现代电力系统——潮流计算作业
0 序章
作业要求（A组）：
0.1 调用matpower中的runpf函数，分析输入文件中各矩阵定义；
0.2 调用某一个算例，输出潮流结果，并分析。
0.3 完成0.1和0.2的基础上，分析matpower中牛顿法和快速解耦法，给出流程图，写出newtonpf和fdpf函数每行程序定义。
0.4 完成0.3的基础上，制造一个病态潮流算例，并跟踪调试，分析病态原因。
1 分析输入文件中各矩阵的定义
1.1 MATPOWER的安装
MATPOWER工具箱的安装步骤如下：
1）下载matpower压缩包。官方下载网址：https://www.doczj.com/doc/4e11064791.html,/matpower/，目前最新版本为6.0b1，稳定版本为5.1，建议下载稳定版本。
2）解压压缩包，得到文件夹matpower5.1，并将文件夹移动到MATLAB所在路径的toolbox文件夹下。我的路径为：C:\Program Files\MATLAB\R2016a\toolbox。
3）添加地址到MATLAB路径。打开MATLAB，点击“文件”→“设置路径”→“添加并包含子文件夹…”，找到matpower5.1所在的位置，点击“确定”，再点“保存”→“关闭”。
4）测试matpower工具是否安装成功。在MATLAB命令行窗口输入“test_matpower”，出现一系列的测试，均显示“ok”，最后显示“All tests successful (3256 passed, 682 skipped of 3938)”，则表示安装成功。
1.2 矩阵的定义
打开文档“caseformat.m”，或者在MATLAB命令行窗口中输入“help caseformat”，可以得到关于输入矩阵的数据定义。当然，也可以参考docs文件夹下的manual文档，其中对matpower工具箱进行了详细说明。
在matpower中，输入矩阵至少包含三种：母线参数矩阵（Bus Data），发电机参数矩阵（Generator Data），支路参数矩阵（Branch Data）。为了进行最优潮流的相关计算，
输入矩阵还包含发电机费用参数矩阵（generator cost data）。以下对三种基本的输入参数矩阵数据格式进行详细说明。
表1.1 母线参数矩阵主要数据格式说明
表1.2 发电机参数矩阵主要数据格式说明
表1.3 支路参数矩阵主要数据格式说明
1.3 case9数据分析
根据以上分析，打开一个算例，比如默认的case9，进行分析。算例case9.m文件包含两个变量和四个矩阵。其中baseMVA=100，表示功率的基准值为100MVA。三个基本的矩阵定义如下。
表1.4 算例case9的母线参数矩阵
表1.6 算例case9的支路参数矩阵
根据参数矩阵，可以推测出case9的电力系统单线图，如图1.1所示。该系统是一个环形网络，包含三个带有发电机的母线，其中母线1是平衡节点，母线2和3均为PV节点，其他的母线都是PQ节点。所有的母线电压初始幅值均设置为1（p.u.），相角为0度，电压基准值为345kV。系统包含3个负荷，分别是母线5上的负荷为90+j30（MVA），母线7上的负荷为100+j35（MVA），母线9上的负荷为125+j50（MVA）。支路1-4，3-6，8-2只有电抗值，电阻和电纳均为0，可以推测该支路为变压器支路的等效。
图1.1 系统单线图
2 计算潮流并分析
2.1 调用runpf计算case9系统的潮流
在MATLAB命令行窗口输入“runpf”或“runpf(‘case9’)”，或者直接运行“runpf.m”，得到case9系统的潮流计算结果：
MATPOWER Version 5.1, 20-Mar-2015 -- AC Power Flow (Newton)
Newton's method power flow converged in 4 iterations.
Converged in 0.03 seconds
================================================================================
| System Summary |
================================================================================
How many? How much? P (MW) Q (MVAr)
--------------------- ------------------- ------------- -----------------
Buses 9 Total Gen Capacity 820.0 -900.0 to 900.0
Generators 3 On-line Capacity 820.0 -900.0 to 900.0
Committed Gens 3 Generation (actual) 320.0 34.9
Loads 3 Load 315.0 115.0
Fixed 3 Fixed 315.0 115.0
Dispatchable 0 Dispatchable -0.0 of -0.0 -0.0
Shunts 0 Shunt (inj) -0.0 0.0
Branches 9 Losses (I^2 * Z) 4.95 51.31
Transformers 0 Branch Charging (inj) - 131.4
Inter-ties 0 Total Inter-tie Flow 0.0 0.0
Areas 1
Minimum Maximum
------------------------- --------------------------------
Voltage Magnitude 0.958 p.u. @ bus 9 1.003 p.u. @ bus 6
Voltage Angle -4.35 deg @ bus 9 9.67 deg @ bus 2
P Losses (I^2*R) - 2.46 MW @ line 8-9
Q Losses (I^2*X) - 16.74 MVAr @ line 8-2
================================================================================ | Bus Data | ================================================================================ Bus Voltage Generation Load
# Mag(pu) Ang(deg) P (MW) Q (MVAr) P (MW) Q (MVAr)
----- ------- -------- -------- -------- -------- --------
1 1.000 0.000* 71.95 24.07 - -
2 1.000 9.669 163.00 14.46 - -
3 1.000 4.771 85.00 -3.65 - -
4 0.987 -2.407 - - - -
5 0.975 -4.017 - - 90.00 30.00
6 1.003 1.926 - - - -
7 0.986 0.622 - - 100.00 35.00
8 0.996 3.799 - - - -
9 0.958 -4.350 - - 125.00 50.00
-------- -------- -------- --------
Total: 319.95 34.88 315.00 115.00
================================================================================ | Branch Data | ================================================================================ Brnch From To From Bus Injection To Bus Injection Loss (I^2 * Z)
# Bus Bus P (MW) Q (MVAr) P (MW) Q (MVAr) P (MW) Q (MVAr) ----- ----- ----- -------- -------- -------- -------- -------- --------
1 1 4 71.95 24.07 -71.95 -20.75 -0.000 3.32
2 4 5 30.7
3 -0.59 -30.55 -13.69 0.17
4 0.94
3 5 6 -59.45 -16.31 60.89 -12.43 1.449 6.31
4 3 6 85.00 -3.6
5 -85.00 7.89 0.000 4.24
5 6 7 24.11 4.54 -24.01 -24.40 0.095 0.81
6 7 8 -75.99 -10.60 76.50 0.26 0.506 4.29
7 8 2 -163.00 2.28 163.00 14.46 0.000 16.74
8 8 9 86.50 -2.53 -84.04 -14.28 2.465 12.40
9 9 4 -40.96 -35.72 41.23 21.34 0.266 2.26
-------- --------
Total: 4.955 51.31
2.2 潮流计算结果分析
Matpower工具箱的潮流计算结果由四部分组成：程序运行信息，系统概述，母线数据，支路数据。
其中，程序运行信号包含潮流计算类型，使用的迭代算法，迭代次数，所用时间。本次潮流计算是交流潮流计算，采用了Newton法，迭代了4次得到了符合精度要求的结果，耗时0.03s。
系统概述描述了系统的基本信息。包括系统元件的数量，元件的功率大小，电压和功率损耗的极值。如case9系统包含9个节点，3个发电机，3个负荷，9条支路。总装机容量820MW，在线容量820MW，实际发电320MW，负荷消耗有功315MW，总网损4.95MW。母线9上电压幅值最小：0.958（p.u.），电压相角也最小：-4.35°，母线6上电压幅值最大：1.003（p.u.），母线2上相角最大：9.67°，支路8-9上消耗了最多的有功功率：2.46MW。
母线数据包含母线电压结果，发电机输出功率，负荷消耗功率，累计功率。
表2.1 潮流计算结果母线数据
支路数据包含起始母线注入功率、终止母线注入功率和支路上的功率损耗。
表2.1 潮流计算结果支路数据
319.95MW，负荷消耗的总有功功率为315MW，网损为4.955MW。另外可以得到系统的潮流分布图如图2.1所示。
71.95
71.95
1.00
∠
0.987-2.407
∠0.975-4.017
∠
0.958-4.350
∠0.9860.622
∠
1.003 1.926
∠
1.0 4.771
∠
1.09.669
∠
0.996 3.799
∠
图2.1 系统潮流分布图
3 迭代算法分析
3.1 牛顿法分析
打开newtonpf.m文档，可以看到matpower的牛顿法的介绍和代码。函数的输入参数包含系统的节点导纳矩阵，节点的注入复功率，初始电压，平衡节点、PV节点和PQ 节点的标号列向量，以及包含终止误差、最大迭代次数和输出选项的结构体。返回节点电压，收敛标志和迭代次数。
通过分析可以得到matpower的牛顿法的程序流程图，如图3.1所示，这和一般的牛顿法潮流计算程序并没有什么区别。
图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节点的电压相角
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%3d %10.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; ]; % 得到雅克比矩阵
%% 计算修正量
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节点的有功功率误差