现代电力系统潮流计算作业

3 迭代算法分析

3.1 牛顿法分析

打开newtonpf.m文档，可以看到matpower的牛顿法的介绍和代码。函数的输入参数包含系统的节点导纳矩阵，节点的注入复功率，初始电压，平衡节点、PV节点和PQ节点的标号列向量，以及包含终止误差、最大迭代次数和输出选项的结构体。返回节点电压，收敛标志和迭代次数。

通过分析可以得到matpower的牛顿法的程序流程图，如图3.1所示，这和一般的牛顿法潮流计算程序并没有什么区别。

开始给定精度和最大迭代次数给定电压初始值初始化迭代次数t=0计算有功和无功的误差?P??、?Q??tt检查是否达到迭代次数或精度否迭代次数加1，t=t+1 是返回节点电压，收敛标志和迭代次数结束 计算雅克比矩阵J?t?计算修正量dx修正各节点电压

–9–

图3.1 牛顿法潮流计算程序流程图

以下是newtonpf函数的每行程序的定义。

function [V, converged, i] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt) %NEWTONPF 使用完整的牛顿法求解潮流

% [V, CONVERGED, I] = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT) % 通过分别给定完整系统的导纳矩阵（针对所有节点），节点的注入

% 复功率（针对所有节点），节点电压的初始值，和平衡节点、PV节点和PQ节 % 点标号的列向量，求解节点电压。节点电压矢量包含发电机节点（包括平衡 % 节点）的设定值和平衡节点的参考角度，以及幅度的大小和角度的初始值。 % MPOPT是一个MATPOWER选项结构体，可用于设置终止误差限，最大迭代次数和 % 输出选项（有关详细信息，请参阅MPOPTION）。如果未指定此参数，则使用 % 默认选项。最终返回节点电压相量，收敛标志以及迭代次数。 %

% 参考RUNPF.

%% 缺省参数设置

if nargin < 7 % 如果输入参数少于7项

mpopt = mpoption; % 则设置mpopt的缺省值为mpoption end

%% 求解选项

tol = mpopt.pf.tol; % 终止误差限 max_it = mpopt.pf.nr.max_it; % 最大迭代次数

%% 初始化

converged = 0; % 收敛标志位清零，不收敛 i = 0; V = V0;

% 迭代次数清零 % 初始电压值

Va = angle(V); % 电压相位初始值 Vm = abs(V); % 电压幅值初始值

%% 为了更新电压，建立电压的指针 npv = length(pv); % PV节点数目 npq = length(pq); % PQ节点数目

j1 = 1; j2 = npv; %% PV节点的电压相角 j3 = j2 + 1; j4 = j2 + npq; %% PQ节点的电压相角

–10–

j5 = j4 + 1; j6 = j4 + npq; %% PQ节点的电压幅值

%% 计算修正方程式的常数项

mis = V .* conj(Ybus * V) - Sbus; % 计算误差 F = [ real(mis([pv; pq])); % delta(P) 有功误差 imag(mis(pq)) ]; % delta(Q) 无功误差

%% 判断误差

normF = norm(F, inf); % F的无穷范数（等效于取最大值）

if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\n it max P & Q mismatch (p.u.)'); fprintf('\\n---- ---------------------------'); fprintf('\\n= .3e', i, normF); end% 将进度信息输出到文档中 if normF < tol

% 误差是否小于误差限

converged = 1; % 收敛标志置1

if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\nConverged!\\n');% 将该字符串输出到指定文档中 end end

%% 进行牛顿迭代

while (~converged && i < max_it) % 不收敛并且小于最大迭代次数 %% 更新迭代次数

i = i + 1; % 迭代次数加1

%% 更新雅克比矩阵

[dSbus_dVm, dSbus_dVa] = dSbus_dV(Ybus, V); % 计算雅克比矩阵各元素

j11 = real(dSbus_dVa([pv; pq], [pv; pq])); % 得到雅克比矩阵的子矩阵H j12 = real(dSbus_dVm([pv; pq], pq)); % 得到雅克比矩阵的子矩阵N j21 = imag(dSbus_dVa(pq, [pv; pq])); % 得到雅克比矩阵的子矩阵J j22 = imag(dSbus_dVm(pq, pq)); % 得到雅克比矩阵的子矩阵L

J = [ j11 j12;

j21 j22; ]; % 得到雅克比矩阵

%% 计算修正量

–11–

dx = -(J \\ F); % 得到修正量

%% 更新电压

if npv % 如果是PV节点

Va(pv) = Va(pv) + dx(j1:j2); % 更新PV节点的电压相位 end

if npq % 如果是PQ节点

Va(pq) = Va(pq) + dx(j3:j4); % 更新PQ节点的电压相位 Vm(pq) = Vm(pq) + dx(j5:j6); % 更新PQ节点的电压幅值 end

V = Vm .* exp(1j * Va); % 更新所有的节点电压 Vm = abs(V); %% 所有节点电压幅值-->Vm Va = angle(V); %% 所有节点电压相位-->Va

%% 计算修正方程式的常数项

mis = V .* conj(Ybus * V) - Sbus; % 计算误差 F = [ real(mis(pv)); % PV节点的有功功率误差 real(mis(pq)); % PQ节点的有功功率误差 imag(mis(pq)) ]; % PQ节点的无功功率误差

%% 校验收敛性

normF = norm(F, inf); % F的无穷范数（等效于取最大值）

if mpopt.verbose > 1 % 如果该标志位大于1，则保存进度信息到文档 fprintf('\\n= .3e', i, normF); % 将进度信息存到文档中 end

if normF < tol % 误差是否小于误差限 converged = 1; % 收敛标志置1

if mpopt.verbose % 如果该标志位不等于0，则保存最终信息到文档

fprintf('\\nNewton''s method power flow converged in %d iterations.\\n', i); end % 将该字符串输出到指定文档中 end end

if mpopt.verbose % 如果该标志位不等于0，则保存最终信息到文档

if ~converged % 如果收敛标志为0，即超出了迭代次数，则输出不收敛的信息到文档 fprintf('\\nNewton''s method power flow did not converge in %d iterations.\\n', i); end end

–12–