电力系统状态估计/Matlab编程 MATLAB编程计算电力系统状态估计
电力系统状态估计这玩意儿就像给电网做体检，得从一堆带噪声的量测数据里算出系统真实状态。咱今天就用Matlab撸个最基础的加权最小二乘法（WLS）状态估计，手把手教你怎么把教科书算法变成跑得通的代码。
先上硬菜——算法核心就三句话：构建量测方程、拼雅可比矩阵、迭代求解修正量。但真写起代码来，魔鬼全在细节里。比如下面这个量测数据的结构体，用表格形式存节点参数和量测值最省事：
% 节点数据 [节点编号 类型 电压幅值 电压相角] bus = [ 1 1 1.06 0; 2 3 1.04 0; 3 2 1.01 0; ]; % 量测数据 [起始节点 结束节点 类型 值 标准差] measure = [ 1 2 1 0.5 0.02; % P12 2 3 2 0.3 0.01; % Q23 1 0 3 1.02 0.005; % V1 ];这里类型1是线路有功，2是无功，3是电压量测。注意量测的方差处理——咱直接取标准差的平方构造权重矩阵：
W = diag(1./measure(:,5).^2); % 权重矩阵状态变量初始化也有讲究，直接拿电压量测值当初始值能减少迭代次数。比如用第三个量测的1.02作为节点1电压初值：
V = bus(:,3); % 电压幅值初始 theta = zeros(size(bus,1),1); % 相角初始构建雅可比矩阵H是重头戏。这里有个暴力但管用的办法：遍历每个量测单独计算偏导。比如处理线路有功量测时：
for m = 1:size(measure,1) i = measure(m,1); % 起始节点 j = measure(m,2); % 结束节点 type = measure(m,3); if type == 1 % 线路有功 G = 1.0; B = -0.5; % 假设导纳 H(m,i) = -B*V(i)*V(j)*cos(theta(i)-theta(j)) + G*V(i)*V(j)*sin(theta(i)-theta(j)); H(m,j) = B*V(i)*V(j)*cos(theta(i)-theta(j)) - G*V(i)*V(j)*sin(theta(i)-theta(j)); end % 其他量测类型处理... end这种逐个元素填充的方式虽然看着笨，但胜在可读性强，改起来方便。实际工程中当然要用稀疏矩阵加速，但教学代码就别整那么复杂了。
电力系统状态估计/Matlab编程 MATLAB编程计算电力系统状态估计
迭代过程才是见证奇迹的时刻。注意每次更新状态后要重新计算雅可比矩阵和残差：
max_iter = 10; tol = 1e-5; for iter = 1:max_iter [H, h] = build_jacobian(V, theta); % 构造H矩阵和量测函数 r = z - h; % 残差计算 dx = (H'*W*H) \ (H'*W*r); % 关键求解步骤 theta = theta + dx(1:end/2); % 更新相角 V = V + dx(end/2+1:end); % 更新电压 if norm(dx) < tol break; end end这里有个小陷阱：状态变量排列顺序必须是相角在前、电压在后，否则雅可比矩阵维度会乱套。要是遇到矩阵奇异的情况，加个正则化项（比如H'WH + λI）能救命。
跑完迭代后做个残差分析，用个简单粗暴的方法检测异常数据：
bad_index = find(abs(r) > 3*sqrt(diag(inv(H'*W*H)))); disp(['可疑量测编号：', num2str(bad_index')]);最后把估计结果和真值对比下，通常IEEE14节点系统迭代5次内就能收敛到1e-4的精度。不过实际系统中拓扑变化、量测丢包这些破事才是真正的挑战，那得用到鲁棒估计或者人工智能这些高级货了。
完整代码传到了Github（假装有个链接），包含了可视化电压分布的功能。下次可以试试加入PMU量测，搞个混合估计玩玩——那才是现代电网的状态估计该有的样子。
