fprintf('\\n - = .3e .3e', i, normP, normQ); end % 将进度信息输出到文档中
if normP < tol && normQ < tol % 误差是否都小于误差限 converged = 1;
% 收敛标志置1
if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\nConverged!\\n'); % 将该字符串输出到指定文档中 end end
%% 降低B矩阵维度
Bp = Bp([pv; pq], [pv; pq]); % 通过索引，降低Bp的维度 Bpp = Bpp(pq, pq); % 通过索引，降低Bpp的维度
%% 矩阵B的LU分解
if lu_vec % 选择系统支持的lu语句
[Lp, Up, pp, qp ] = lu(Bp, 'vector'); % 矩阵Bp的LU分解 [Lpp, Upp, ppp, qpp] = lu(Bpp, 'vector');% 矩阵Bpp的LU分解 [junk, iqp ] = sort(qp); % [~, iqp ] = sort(qp); % [~, iqpp] = sort(qpp); else
[Lp, Up, Pp] = lu(Bp); % 矩阵Bp的LU分解 [Lpp, Upp, Ppp] = lu(Bpp);% 矩阵Bpp的LU分解 end
%% 进行P和Q迭代
while (~converged && i < max_it) %不收敛并且小于最大迭代次数 %% 更新迭代次数
i = i + 1; % 迭代次数加1
%%----- 进行P迭代，更新Va -----
if lu_vec % 判断之前进行的lu分解操作，以便进行对应的后续操作 dVa = -( Up \\ (Lp \\ P(pp)) ); % 求电压相角的修正量 dVa = dVa(iqp); % 对电压相角的修正量向量重新排序 else
dVa = -( Up \\ (Lp \\ (Pp * P))); % 求电压相角的修正量 end
% 行向量qp的排序 % junk没有用到，用~代替 % junk没有用到，用~代替
[junk, iqpp] = sort(qpp); % 行向量qpp的排序
–16–
%% 更新电压
Va([pv; pq]) = Va([pv; pq]) + dVa; % 更新电压相角 V = Vm .* exp(1j * Va); % 更新电压值
%% 计算误差
mis = (V .* conj(Ybus * V) - Sbus) ./ Vm; % 计算误差 P = real(mis([pv; pq])); % delta(P) 有功误差 Q = imag(mis(pq)); % delta(Q) 无功误差
%% 校验误差
normP = norm(P, inf); % P的无穷范数（等效于取最大值） normQ = norm(Q, inf); % Q的无穷范数（等效于取最大值）
if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\n P = .3e .3e', i, normP, normQ); end % 将进度信息输出到文档中
if normP < tol && normQ < tol % 误差是否小于误差限 converged = 1; % 收敛标志置1
if mpopt.verbose % 如果该标志位不等于0，则保存最终信息到文档
fprintf('\\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\\n', i, i-1);
end % 将该字符串输出到指定文档中 break; % 跳出循环，结束迭代 end
%%----- 进行Q迭代，更新Vm -----
if lu_vec %判断之前进行的lu分解操作，以便进行对应的后续操作 dVm = -( Upp \\ (Lpp \\ Q(ppp)) ); %求电压幅值的修正量 dVm = dVm(iqpp); % 对电压幅值的修正量向量重新排序 else
dVm = -( Upp \\ (Lpp \\ (Ppp * Q)) ); %求电压幅值的修正量 end
%% 更新电压
Vm(pq) = Vm(pq) + dVm; % 更新电压相角 V = Vm .* exp(1j * Va); % 更新电压值
%% 计算误差
–17–
mis = (V .* conj(Ybus * V) - Sbus) ./ Vm; % 计算误差 P = real(mis([pv; pq])); % delta(P) 有功误差 Q = imag(mis(pq)); % delta(Q) 无功误差
%% 校验误差
normP = norm(P, inf); % P的无穷范数（等效于取最大值） normQ = norm(Q, inf); % Q的无穷范数（等效于取最大值）
if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\n Q = .3e .3e', i, normP, normQ); end % 将进度信息输出到文档中
if normP < tol && normQ < tol % 误差是否小于误差限 converged = 1; % 收敛标志置1
if mpopt.verbose % 如果该标志位不等于0，则保存最终信息到文档
fprintf('\\nFast-decoupled power flow converged in %d P-iterations and %d Q-iterations.\\n', i, i);
end % 将该字符串输出到指定文档中 break; % 跳出循环，结束迭代 end end
if mpopt.verbose % 如果该标志位不等于0，则保存最终信息到文档
if ~converged % 如果收敛标志为0，即超出了迭代次数，则输出不收敛的信息到文档 fprintf('\\nFast-decoupled power flow did not converge in %d iterations.\\n', i); end end
4 病态潮流分析
4.1 重负荷导致的病态潮流
如图4.1所示，将原9节点的系统9#母线上的负荷由原来的125+j50改为125+j300，便构成了一个病态问题。如表4.1所示，逐渐增大9#母线上的负荷，可以清楚地发现节点电压的幅值迅速减小，相位逐渐滞后，牛顿法迭代次数迅速增大，直到问题不收敛，雅克比矩阵的条件数也迅速增大。经过仿真分析，可以发现，雅克比矩阵的条件数最大达到了1?106，所以修正方程是一个病态方程。该问题在给定的运行条件下，潮流方程无解，如果想让系统有解，则必须调整系统的参数，使得系统回到有解区域。
–18–
表4.1 9节点系统9#母线负荷对系统潮流计算的影响
负荷（MVA） 125+j50 125+j100 125+j150 125+j200 125+j250 125+j280 125+j300 节点电压(p.u.) 0.958?-4.350°0.907?-4.408° 0.850?-4.551° 0.780?-4.855° 0.684?-5.579° 0.574?-7.010° - r?0.017x?0.092b?0.158牛顿法迭代次数 4 4 4 5 5 7 迭代200次不收敛 90?j30?MVA?雅克比矩阵最大条件数 101 102 104 107 111 165 ?106 x?0.0586G1x?0.057614r?0.01x?0.085b?0.1765r?0.039x?0.17b?0.358G36r?0.0119x?0.1008b?0.2093P=85MW9125?j300?MVA?r?0.032x?0.161b?0.306r?0.00857x?0.072b?0.1498x?0.0625100?j35?MVA?2G2P=163MW图4.1 重负荷导致的病态问题系统单线图
通过设置断点、单步调试和查看变量等方式，可以清楚地看到9#节点重负荷情况下，系统潮流的计算过程。如表4.2所示。第7次迭代时，雅克比矩阵接近奇异，条件数急剧增大，节点电压也变得非常大，问题无解。
表4.2 9节点系统9#母线电压随迭代次数的变化
迭代次数 0 1 2 3 4 5 6 7 8 9#节点电压（p.u.） 1?0° 0.754?-2.965° 0.607?-5.548° 0.491?-7.970° 0.787?-1.644° 0.628?-4.950° 0.515?-7.425° 4.752?59.738° 2.558?30.695° 雅克比矩阵的条件数 - 99.05 108.28 116.62 472.34 107.17 114.69 8464.3 1596.9 –19–
9 10 11 12 13 14 15 16 17 18 1.416?173.727° 0.626?158.783° 0.125?73.611° 1.151?-120.745° 0.524?174.071° 0.060?48.981° 1.729?19.571° 0.808?17.927° 0.096?102.919° 0.005?-95.546° 9753.4 2367.7 276170 7697.4 598.7 989.3 5227.1 1416.2 13329 4813.0 4.2 小阻抗支路导致的病态潮流
还有一类病态问题是系统有解，但是由于潮流方程自身的一些原因或是算法本身不完善等，导致求不到解。比如系统包含小阻抗支路。如图4.2所示，将原9节点的系统支路4-9的电阻改为0，电抗改为10-9，便构成了一条小阻抗支路，通过潮流计算可以发现这个问题是个病态问题。
G1x?0.057690?j30?MVA?14r?0r?0.017x?0.092b?0.1585r?0.039x?0.17b?0.358x?0.0586G36r?0.0119x?0.1008b?0.2093P=85MWx?10b?0.176-99125?j50?MVA?r?0.032x?0.161b?0.306r?0.00857x?0.072b?0.1498x?0.0625100?j35?MVA?2G2P=163MW图4.2 小阻抗支路导致的病态问题系统单线图
使用牛顿法对该系统进行潮流计算，结果显示迭代了20次，系统不收敛。但是，通过跟踪调试可以发现，系统的电压在第4次迭代时便已经达到了很高的精度，只是因为存在一条小阻抗支路，导致雅克比矩阵的条件数达到了2.26?109（这和1x是同一个量级），导致了牛顿潮流病态。表4.3是节点电压的计算结果，支路功率无法计算。
–20–