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

3.2 快速解耦法分析

在matpower中，快速解耦法分为XB型和BX型。调用方式如下： runpf('case9',mpoption('pf.alg','FDXB')); % XB型 runpf('case9',mpoption('pf.alg','FDBX')); % BX型

快速解耦法、牛顿法和高斯-塞德尔法的比较如表3.1所示，对于不同节点的算例，牛顿法和快速解耦法都具有很强的收敛性，迭代次数并没有明显的增加，而高斯-塞德尔法收敛性不佳，迭代次数增加明显，甚至超过了预设的最大迭代次数。

表3.1 各种潮流计算方法的比较

算例节点数 9 24 57 118 300 牛顿法 4 4 3 3 5 FDXB 8-7 8-8 7-7 8-7 9-9 FDBX 7-7 9-8 9-8 7-6 9-8 高斯-塞德尔法 212 535 518 NC NC 接下来对快速解耦法进行分析，得到流程图如图3.2所示。

–13–

开始给定精度和最大迭代次数给定电压初始值初始化迭代次数t=0计算有功和无功的误差?P??、?Q??tt计算有功和无功的误差?P??、?Q??tt检查是否达到迭代次数或精度否迭代次数加1，t=t+1进行有功功率迭代计算电压相位修正量Va 是返回节点电压，收敛标志和迭代次数是检查是否达到迭代次数或精度否结束进行无功功率迭代计算电压幅值修正量Vm修正各节点电压修正各节点电压

图3.2 快速解耦法潮流计算程序流程图

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

function [V, converged, i] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt) yPF 使用快速解耦法求解电力系统潮流.

% 通过分别给定完整系统的导纳矩阵（针对所有节点），节点的注入

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

% 参考RUNPF.

–14–

%% 缺省参数设置

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

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

%% 求解选项

tol = mpopt.pf.tol; % 终止误差限

max_it = mpopt.pf.fd.max_it; % 最大迭代次数

if have_fcn('matlab') && have_fcn('matlab', 'vnum') < 7.3 % 测试matlab的版本信息和功能

lu_vec = 0; %% 不支持lu(..., 'vector') 语法 else

lu_vec = 1; %% 支持lu(..., 'vector') 语法 end

%% 初始化

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

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

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

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

%% 校验初始误差

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('\\niteration max mismatch (p.u.) '); fprintf('\\ntype # P Q '); fprintf('\\n---- ---- ----------- -----------');

–15–

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–