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

1、 现代电力系统潮流计算作业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压缩包。官方下载网址：，目前最新版本为 6.0b1，稳定版本为 5.1，建议下载稳定

2、版本。2）解压压缩包，得到文件夹 matpower5.1，并将文件夹移动到 MATLAB所在路径的toolbox文件夹下。我的路径为： C:Program FilesMATLABR2016atoolbox。3）添加地址到 MATLAB路径。打开 MATLAB，点击“文件”“设置路径”“添加并包含子文件夹”，找到 matpower5.1所在的位置，点击“确定”，再点“保存”“关闭”。4）测试 matpower工具是否安装成功。在 MATLAB命令行窗口输入“test_matpower”，出现一系列的测试，均显示“ ok”，最后显示“ All tests successful (3256 pass

3、ed, 682skipped of 3938)”，则表示安装成功。1.2矩阵的定义打开文档“caseformat.m”，或者在 MATLAB命令行窗口中输入“help caseformat”，可以得到关于输入矩阵的数据定义。当然，也可以参考docs文件夹下的 manual文档，其中对 matpower工具箱进行了详细说明。在 matpower中，输入矩阵至少包含三种：母线参数矩阵（ Bus Data），发电机参数矩阵（Generator Data），支路参数矩阵（Branch Data）。为了进行最优潮流的相关计算，1 输入矩阵还包含发电机费用参数矩阵（ generator cost data

4、）。以下对三种基本的输入参数矩阵数据格式进行详细说明。表 1.1母线参数矩阵主要数据格式说明含义列名称单位1bus number母线编号，从 1开始，正整数节点类型， 1为 PQ节点， 2为 PV节点， 3为平衡(参考 )节点， 4为孤立节点2bus type34PdQd注入负荷的有功功率注入负荷的无功功率与母线并联的电导与母线并联的电纳MWMVarMW5Gs6BsMVar7area numberVm电网断面号，一般设置为母线电压的幅值118p.u.度9Va母线电压的相角10111213baseKVzone基准电压kV损耗分区号，一般设置为VmaxVmin工作时，母线的最高电压幅值工作时，母线

5、的最低电压幅值p.u.p.u.表 1.2发电机参数矩阵主要数据格式说明列名称含义单位1bus发电机所在母线的编号234567PgQg接入发电机的有功功率接入发电机的无功功率发电机的最大输出无功功率发电机的最小输出无功功率发电机的工作电压幅值MWMVarMVarMVarp.u.QmaxQminVgmBase发电机的功率基准值，默认为baseMVAMVA发电机是否投入运行，大于 0表示投入运行，小于或等于 0表示退出运行8status9PmaxPmin发电机的最大输出有功功率发电机的最小输出有功功率MWMW102 表 1.3支路参数矩阵主要数据格式说明含义列名称单位1fbus起始母线编号终止母线编

6、号支路电阻23tbusrp.u.4x支路电抗p.u.5b支路充电电纳p.u.6rateArateBrateCratioanglestatusangmin支路长期运行允许的功率MVAMVAMVA7支路短期运行允许的功率8支路紧急运行允许的功率9支路上变压器的变比，如果支路元件不是变压器，则该值为支路上变压器的转角，如果支路元件不是变压器，则该值为支路的初始工作状态， 1表示投入运行， 0表示退出运行支路最小相角差0010111213度度度angmax支路最大相角差1.3 case9数据分析根据以上分析，打开一个算例，比如默认的 case9，进行分析。算例 case9.m文件包含两个变量和四个矩阵

7、。其中 baseMVA=100，表示功率的基准值为 100MVA。三个基本的矩阵定义如下。表 1.4 算例 case9的母线参数矩阵12type33Pd04Qd05Gs06Bs07area18Vm19Va010baseKV3451112Vmax1.11.11.11.11.11.11.11.11.113Vmin0.90.90.90.90.90.90.90.90.9bus_izone1234567891111111112000011034520000110345100001103451900300001103451001103451100035000110345100110345112550001

8、103453 表 1.5 算例 case9的发电机参数矩阵123456789101112bus Pg Qg Qmax Qmin Vg mBase statusPmax Pmin Pc1 Pc21230000300300300-300-300-30011110010010011125030027010101000000016385131415161718192021Qc1min Qc1max Qc2min Qc2max ramp_agc ramp_10 ramp_30 ramp_q apf000000000000000000000000000表 1.6 算例 case9的支路参数矩阵1fbus1

9、2tbus43r4x5b06rateA2502501503001502502502502507rateB2502501503001502502502502508rateC2502501503001502502502502509ratio010111213anglestatusangmin angmax00.05760.0920.17000000000111111111-360-360-360-360-360-360-360-360-360360360360360360360360360360450.0170.03900.1580.35800560360.05860.10080.0720.062

10、50.1610.0850670.01190.008500.2090.14900780820890.0320.010.3060.1760940根据参数矩阵，可以推测出 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

11、，3-6，8-2只有电抗值，电阻和电纳均为 0，可以推测该支路为变压器支路的等效。4 90 j30 MVA0.039x 0.17b 0.358x 0.0576G1x 0.05860.0119rx 0.092b 0.1580.017rG3514P=85MW63r0.01rx 0.085b 0.176x 0.1008b 0.209r0.032rx0.00850.07279x 0.161b 0.306b 0.149125 j50 MVA100 j35 MVA8x 0.06252P=163MWG2图 1.1系统单线图2计算潮流并分析2.1调用 runpf计算 case9系统的潮流在 MATLAB命令行

12、窗口输入“runpf”或“runpf(case9”，或者直接运行“runpf.m”，得到 case9系统的潮流计算结果：MATPOWER Version 5.1, 20-Mar-2015 - AC Power Flow (Newton)Newtons method power flow converged in 4 iterations.Converged in 0.03 seconds=|System Summary|=How many?-BusesHow much?P (MW)Q (MVAr)- - -9Total Gen CapacityOn-line CapacityGeneratio

13、n (actual)Load820.0820.0320.0-900.0 to 900.0GeneratorsCommitted GensLoads3-900.0 to 900.034.933315.0115.0115.0-0.0Fixed3Fixed315.0DispatchableShunts0DispatchableShunt (inj)-0.0 of -0.009-0.04.95-0.0BranchesLosses (I2 * Z)Branch Charging (inj)51.31131.4Transformers05 Inter-tiesAreas0Total Inter-tie F

14、lowMinimum0.00.01Maximum- -0.958 p.u. bus 9 1.003 p.u. bus 6-4.35 deg bus 9Voltage MagnitudeVoltage Angle9.67 deg2.46 MW bus 2P Losses (I2*R)Q Losses (I2*X)- line 8-9-16.74 MVAr line 8-2=|Bus Data|=Bus#VoltageGenerationLoadMag(pu) Ang(deg)P (MW)Q (MVAr)P (MW)Q (MVAr)- - - - - - -1 1.0002 1.0003 1.00

15、04 0.9875 0.9756 1.0037 0.9868 0.9969 0.9580.000*9.6694.771-2.407-4.0171.9260.6223.799-4.35071.9524.07-163.0014.46-85.00-3.65-30.00-90.00-100.00-35.00-125.0050.00- - - -Total:319.9534.88315.00115.00=|Branch Data|=Brnch#FromBusToFrom Bus InjectionP (MW) Q (MVAr)To Bus InjectionLoss (I2 * Z)BusP (MW)Q

16、 (MVAr)P (MW)Q (MVAr)- - - - - - - - -12345678914536788945667829471.9530.7324.07-0.59-16.31-3.654.54-71.95-30.5560.89-20.75-13.69-12.437.89-0.0000.1741.4490.0000.0950.5060.0002.4650.2663.320.946.314.240.814.2916.7412.402.26-59.4585.00-85.00-24.0176.5024.11-24.400.26-75.99-163.0086.50-10.602.28163.00

17、-84.0441.2314.46-14.2821.34-2.53-35.72-40.966 - -Total:4.95551.312.2潮流计算结果分析Matpower工具箱的潮流计算结果由四部分组成：程序运行信息，系统概述，母线数据，支路数据。其中，程序运行信号包含潮流计算类型，使用的迭代算法，迭代次数，所用时间。本次潮流计算是交流潮流计算，采用了 Newton法，迭代了 4次得到了符合精度要求的结果，耗时 0.03s。系统概述描述了系统的基本信息。包括系统元件的数量，元件的功率大小，电压和功率损耗的极值。如 case9系统包含 9个节点， 3个发电机， 3个负荷， 9条支路。总装机容量 8

18、20MW，在线容量 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潮流计算结果母线数据母线电压幅值pu发电机有功发电机无功负荷消耗有负荷消耗无功电压相角 deg编号1MWMVar功 MWMVar1.0001.0001.0000.9870.9751.0030.9860.996

19、0.9580.000*9.6694.771-2.407-4.0171.9260.6223.799-4.350累计：71.9524.07-2163.0014.46-385.00-3.65-4-5-90.00-30.00-6-7-100.00-35.00-8-9-125.00315.0050.00115.00319.9534.88支路数据包含起始母线注入功率、终止母线注入功率和支路上的功率损耗。表 2.1潮流计算结果支路数据支路编号1起始母线1终止母线4起始母线注起始母线注终止母线注终止母线注支路有功支路无功入有功 MW入无功 MVar 入有功 MW入无功 MVar损耗 MW损耗 MVar71.9

20、524.07-71.95-20.75-0.0003.327 23456789453678895667829430.73-59.4585.00-0.59-16.31-3.654.54-30.5560.89-85.00-24.0176.50163.00-84.0441.23-13.69-12.437.890.1741.4490.0000.0950.5060.0002.4650.2664.9550.946.314.2424.11-24.400.260.81-75.99-163.0086.50-10.602.284.2914.4616.7412.402.26-2.53-35.72-14.2821.34

21、-40.96累计损耗：51.31通过以上数据可以发现系统的有功功率守恒：发电机发出的总有功功率为319.95MW，负荷消耗的总有功功率为 315MW，网损为 4.955MW。另外可以得到系统的潮流分布图如图 2.1所示。0.975 -4.01759.451.003 1.9260.987 -2.4071.0 030.73 j0.5990 j30 MVA1.0 4.77171.95 j24.07j16.31G185.0 j7.8930.55+j13.6960.89 j12.43G371.95 j20.7571.95 j24.07585.0 j3.651441.23 j21.3463 85.0 j3

22、.6524.11 j4.5424.01 j24.400.986 0.62240.96 j35.72984.04 j14.28775.99 j10.600.958 -4.35086.50j2.53125 j50 MVA100 j35 MVA76.50 j0.268163.0 j2.280.996 3.799163.0 j14.4621.0 9.669G2 163.0 j14.46图 2.1系统潮流分布图3迭代算法分析3.1牛顿法分析打开 newtonpf.m文档，可以看到 matpower的牛顿法的介绍和代码。函数的输入参数包含系统的节点导纳矩阵，节点的注入复功率，初始电压，平衡节点、 PV节点

23、和 PQ节点的标号列向量，以及包含终止误差、最大迭代次数和输出选项的结构体。返回节点电压，收敛标志和迭代次数。8 通过分析可以得到 matpower的牛顿法的程序流程图，如图 3.1所示，这和一般的牛顿法潮流计算程序并没有什么区别。开始给定精度和最大迭代次数给定电压初始值初始化迭代次数 t=0计算有功和无功的误差ttP、 Q返回节点电压，收敛标志和迭代次数检查是否达到迭代次数或精度是否迭代次数加 1，t=t+1结束计算雅克比矩阵tJ计算修正量 dx修正各节点电压图 3.1牛顿法潮流计算程序流程图以下是 newtonpf函数的每行程序的定义。function V, converged, i =

24、newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt)%NEWTONPF使用完整的牛顿法求解潮流%V, CONVERGED, I = NEWTONPF(YBUS, SBUS, V0, REF, PV, PQ, MPOPT)通过分别给定完整系统的导纳矩阵（针对所有节点），节点的注入复功率（针对所有节点），节点电压的初始值，和平衡节点、PV节点和 PQ节点标号的列向量，求解节点电压。节点电压矢量包含发电机节点（包括平衡9 %节点）的设定值和平衡节点的参考角度，以及幅度的大小和角度的初始值。MPOPT是一个 MATPOWER选项结构体，可用于设置终止误差限，最大迭代

25、次数和输出选项（有关详细信息，请参阅MPOPTION）。如果未指定此参数，则使用默认选项。最终返回节点电压相量，收敛标志以及迭代次数。参考 RUNPF.%缺省参数设置if nargin 1 %如果该标志位大于 1，则保存进度信息到文档10 fprintf(n itmax P & Q mismatch (p.u.);-);%10.3e, i, normF);fprintf(n-fprintf(n%3dend%将进度信息输出到文档中if normF 1 %如果该标志位大于 1，则保存进度信息到文档fprintf(nConverged!n);%将该字符串输出到指定文档中endend%进行牛顿迭代wh

26、ile (converged & i Vm%所有节点电压相位 -Va%计算修正方程式的常数项mis = V .* conj(Ybus * V) - Sbus; %计算误差F = real(mis(pv); % PV节点的有功功率误差real(mis(pq);imag(mis(pq)% PQ节点的有功功率误差; % PQ节点的无功功率误差%校验收敛性normF = norm(F, inf); % F的无穷范数（等效于取最大值）if mpopt.verbose 1 %如果该标志位大于 1，则保存进度信息到文档fprintf(n%3d%10.3e, i, normF); %将进度信息存到文档中end

27、if normF tol %误差是否小于误差限converged = 1; %收敛标志置 1if mpopt.verbose %如果该标志位不等于 0，则保存最终信息到文档fprintf(nNewtons method power flow converged in %d iterations.n, i);end %将该字符串输出到指定文档中endendif mpopt.verbose%如果该标志位不等于 0，则保存最终信息到文档if converged %如果收敛标志为 0，即超出了迭代次数，则输出不收敛的信息到文档fprintf(nNewtons method power flow did

28、 not converge in %d iterations.n, i);endend3.2快速解耦法分析在 matpower中，快速解耦法分为 XB型和 BX型。调用方式如下：runpf(case9,mpoption(pf.alg,FDXB); % XB型runpf(case9,mpoption(pf.alg,FDBX); % BX型快速解耦法、牛顿法和高斯 -塞德尔法的比较如表 3.1所示，对于不同节点的算例，牛顿法和快速解耦法都具有很强的收敛性，迭代次数并没有明显的增加，而高斯尔法收敛性不佳，迭代次数增加明显，甚至超过了预设的最大迭代次数。-塞德12 表 3.1各种潮流计算方法的比较算例

29、节点数牛顿法FDXB8-7FDBX7-7高斯 -塞德尔法944335212535518NCNC248-89-8577-79-81183008-77-69-99-8接下来对快速解耦法进行分析，得到流程图如图3.2所示。开始给定精度和最大迭代次数给定电压初始值初始化迭代次数 t=0计算有功和无功的误差计算有功和无功的误差ttttQP、 QP、返回节点电压，收敛标志和迭代次数检查是否达到迭代次数或精度检查是否达到迭代次数或精度是是否否迭代次数加 1，t=t+1进行无功功率迭代计算电压幅值修正量 Vm结束进行有功功率迭代计算电压相位修正量 Va修正各节点电压修正各节点电压图 3.2快速解耦法潮流计算程

30、序流程图以下是 fdpf函数的每行程序的定义。function V, converged, i = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt)%FDPF使用快速解耦法求解电力系统潮流.13 %通过分别给定完整系统的导纳矩阵（针对所有节点），节点的注入复功率（针对所有节点），节点电压的初始值，和平衡节点、PV节点和 PQ节点标号的列向量，求解节点电压。节点电压矢量包含发电机节点（包括平衡节点）的设定值和平衡节点的参考角度，以及幅度的大小和角度的初始值。MPOPT是一个 MATPOWER选项结构体，可用于设置终止误差限，最大迭代次数和输出选项

31、（有关详细信息，请参阅MPOPTION）。如果未指定此参数，则使用默认选项。最终返回节点电压相量，收敛标志以及迭代次数。参考 RUNPF.%缺省参数设置if nargin 7 %如果输入参数少于 7项mpopt = mpoption; %则设置 mpopt的缺省值为 mpoptionend%求解选项tol= mpopt.pf.tol; %终止误差限max_it = mpopt.pf.fd.max_it; %最大迭代次数if have_fcn(matlab) & have_fcn(matlab, vnum) 1 %如果该标志位大于 1，则保存进度信息到文档fprintf(niterationfp

32、rintf(ntypemax mismatch (p.u.)P);#Q);fprintf(n- - -);fprintf(n- %3dend %将进度信息输出到文档中if normP tol & normQ 1 %如果该标志位大于 1，则保存进度信息到文档fprintf(nConverged!n); %10.3e%10.3e, i, normP, normQ);将该字符串输出到指定文档中endend%降低 B矩阵维度Bp = Bp(pv; pq, pv; pq);%通过索引，降低 Bp的维度Bpp = Bpp(pq, pq); %通过索引，降低 Bpp的维度%矩阵 B的 LU分解if lu_v

33、ec %选择系统支持的 lu语句Lp, Up, pp, qp = lu(Bp, vector); %矩阵 Bp的 LU分解Lpp, Upp, ppp, qpp = lu(Bpp, vector);%矩阵 Bpp的 LU分解junk, iqp = sort(qp);junk, iqpp = sort(qpp);% , iqp = sort(qp);% , iqpp = sort(qpp);%行向量 qp的排序%行向量 qpp的排序% junk没有用到，用 代替% junk没有用到，用 代替elseendLp, Up, Pp = lu(Bp); %矩阵 Bp的 LU分解Lpp, Upp, Ppp = lu(Bpp);%矩阵 Bpp的 LU分解%进行 P和 Q迭代while (converged & i 1 %如果该标志位大于 1，则保存进度信息到文档fprintf(nP %3d%10.3e%10.3e, i, normP, normQ);end %将进度信息输出到文档中if normP tol & normQ 1 %如果该标志位大于 1，则保存进度信息到文档Q %3d %10.3e, i, normP, normQ);end %将进度信息输出到文档中% Q的无穷范数（等效于取最大值）fprintf(n%10.3eif normP tol & normQ to