下载文档原格式
含风电场电力系统的潮流计算方法综述【摘要】含风电场电力系统的潮流计算对分析风电场并网后对电力系统稳定运行的影响具有十分重要的意义。
本文在介绍在风电场潮流计算模型的基础上，分析和总结了目前含风电场电力系统的确定性潮流计算方法和不确定性潮流计算方法的研究现状和特点。
【关键词】风电场；电力系统；确定性潮流；不确定性潮流0.引言随着能源结构的调整，风力发电的比重日益加大，我国政府一直积极支持风力发电的发展，制定了一系列鼓励风力发电的政策，支持风力发电的快速发展，使得风力发电的成本已大幅下降，成为了可再生能源中发展速度最快和最有前途的发电方式之一[1]。
但风能所具有的随机性和不可控性决定了风电机组的出力具有波动性和间歇性的特点。
与传统发电方式相比，风电场容量可信度较低，随着风电场规模的不断扩大和风力发电装置容量的显著增加，风电并网后对原有电力系统的影响也加大了，因此对其的研究也更加迫切。
当风电机组装机容量在电网总容量的比例较大时，风力发电将改变输电系统中的网损及其原有的潮流分布，输电网运行的安全性将受到较大的挑战，其运行的经济性也可能受到一定的影响[2-5]。
因此，为了研究风电机组接入电网以后对整个电力系统的影响，就必须对风电并网前后的系统潮流分布进行计算。
目前，对风电场潮流计算的研究已经具有一定的基础，风电场潮流计算主要包括含普通异步电机的风电场潮流计算和含双馈异步电机的风电场潮流计算[6]。
从上世纪80年代起，随着并网风电场的出现，人们就开始关注含风电场电力系统的潮流计算问题。
在电力系统潮流计算中，传统节点主要分为PV节点、PQ节点和平衡节点。
一般异步电机本身没有励磁调节装置，不能有效地调整节点电压，因此不能与常规的同步电机一样看作电压幅值恒定的PV节点。
异步电机向系统注入有功功率时也要从系统吸收一定的无功功率，吸收无功大小与发电机发出的有功功率、滑差率和机端电压等有着紧密的联系，因此不能简单的处理为恒功率的PQ节点[7]。
含风电场的电力系统潮流计算一、本文概述随着全球能源结构的转型和可再生能源的大力发展，风电作为一种清洁、可再生的能源形式，其在电力系统中的比重日益增加。
风电场的大规模接入对电力系统的运行和控制带来了新的挑战，尤其是风电场出力的随机性和波动性对电力系统的潮流分布、电压稳定性以及保护控制等方面产生了显著影响。
因此，对含风电场的电力系统进行准确的潮流计算，对于电力系统的规划、设计、运行和控制具有重要的理论价值和现实意义。
本文旨在研究含风电场的电力系统潮流计算方法，分析风电场接入对电力系统潮流分布的影响，提出相应的潮流计算模型和算法。
文章首先介绍了风电场的基本特性及其在电力系统中的接入方式，然后详细阐述了含风电场的电力系统潮流计算的基本原理和方法，包括风电场出力模型的建立、潮流计算的基本方程和求解算法等。
在此基础上，文章进一步探讨了风电场接入对电力系统潮流分布的影响，包括风电场出力波动对电压稳定性、线路潮流和节点功率分布的影响等。
文章提出了针对含风电场的电力系统潮流计算的一些改进措施和优化策略，为提高电力系统的运行效率和稳定性提供参考。
通过本文的研究，可以为含风电场的电力系统潮流计算提供理论支持和实践指导，有助于更好地理解和解决风电场接入带来的电力系统运行问题，推动可再生能源在电力系统中的广泛应用和持续发展。
二、风电场特性及建模风电场作为可再生能源的重要组成部分，具有随机性、间歇性和不可预测性等特点。
这些特性使得风电场在电力系统中的建模和潮流计算变得复杂。
风电场的出力受到风速、风向、湍流等多种因素的影响，因此，准确描述风电场的特性并建立合适的模型是电力系统潮流计算的关键。
在风电场建模中，通常将风电场看作一个由多个风电机组组成的集合。
每个风电机组的出力取决于其装机容量、风速以及控制策略等因素。
为了简化计算，通常将风电场视为一个等效的电源，其出力等于所有风电机组出力的总和。
等效电源的出力特性可以通过统计方法得到，如威布尔分布、贝塔分布等。
含风力发电机组的配电网潮流计算一、概述随着全球能源结构的转型和可再生能源的大力发展，风电作为一种清洁、可再生的能源形式，其在电力系统中的比重日益增加。
风电场的大规模接入为电力系统带来了新的活力，但同时也带来了诸多挑战。
尤其在配电网层面，风力发电机组的接入使得配电网从一个无源网络转变为有源网络，其潮流特性、电压分布以及网损情况都发生了显著变化。
含风力发电机组的配电网潮流计算，是电力系统分析与控制领域的重要课题。
通过潮流计算，可以准确描述风力发电机组接入后配电网的运行状态，分析其对系统电压稳定性、潮流分布以及网损的影响。
这不仅有助于电力系统的规划与设计，更对于电力系统的安全稳定运行和优化调度具有重要意义。
在含风力发电机组的配电网潮流计算中，风电场的特性建模是关键环节。
由于风速的随机性、间歇性和不可预测性，风电场的出力具有极大的不确定性。
在建模过程中需要充分考虑这些因素，建立准确的风电场出力模型。
配电网的结构特点、负荷分布以及控制策略等也是影响潮流计算的重要因素。
针对含风力发电机组的配电网潮流计算已有多种方法，如前推回代法、牛顿拉夫逊法等。
这些方法各有优缺点，需要根据实际情况进行选择和优化。
随着智能电网和分布式发电技术的不断发展，配电网潮流计算也面临着新的挑战和机遇。
本文旨在深入研究含风力发电机组的配电网潮流计算方法，分析风力发电机组接入对配电网潮流分布的影响，提出相应的优化策略和建议。
通过本文的研究，可以为含风力发电机组的配电网潮流计算提供理论支持和实践指导，有助于推动可再生能源在电力系统中的广泛应用和持续发展。
1. 风力发电机组在配电网中的应用背景随着全球能源结构的转型和可再生能源的大力发展，风力发电作为一种清洁、可再生的能源形式，其在配电网中的应用愈发广泛。
风力发电机组，作为风力发电的核心设备，在配电网中发挥着举足轻重的作用。
环境问题日益严重，化石燃料燃烧导致的碳排放量不断增加，加剧了全球气候变暖的速度。
班级：姓名：学号：一、作业要求编写程序计算图1所示算例系统的潮流及三相短路电流..潮流计算：方法不限;计算系统的节点电压和相角..短路电流：4号母线发生金属性三相短路时z f=0;分别按照精确算法和近似算法计算短路电流、系统中各节点电压以及网络中各支路的电流分布;并对两种情况下的计算结果进行比较..二、电路图及参数图1 3机9节点系统表1 9节点系统支路参数表2 9节点系统发电机参数表3 9节点系统负荷参数三、计算步骤(1) 进行系统正常运行状态的潮流计算;求得(0)i U (2) 形成不含发电机和负荷的节点导纳矩阵Y N ；(3) 将发电机表示为电流源i I /i diE jx ''=和导纳i y 1/di jx '=的并联组合；节点负荷用恒阻抗的接地支路表示；形成包括所有发电机支路和负荷支路的节点导纳矩阵Y;即在Y N 中的发电机节点和负荷节点的自导纳上分别增加发电机导纳i y 和负荷导纳,LD i y *,,22LD i LDi LDiLD i i i S P jQ y V V -==； (4) 利用1Z Y-=;计算节点阻抗矩阵;从而得到阻抗矩阵中的第f 列；(5) 利用公式6-7或6-10计算短路电流；(6)利用公式6-8或6-11计算系统中各节点电压；(7)利用公式6-9计算变压器支路的电流；对输电线路利用П型等值电路计算支路电流..四、计算结果节点导纳矩阵Yn：Columns 1 through 50 -17.3611i 0 0 0 +17.3611i 00 0 -16.0000i 0 0 00 0 0 -17.0648i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 2.5528 -17.3382i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 03.2242 -15.8409i 0 0 -1.2820 + 5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 2.7722 -23.3032i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i电压幅值：1.0400 1.0250 1.0250 1.0258 0.9956 1.0127 1.0258 1.0159 1.0324电压相角：0 0.1620 0.0814 -0.0387 -0.0696 -0.0644 0.0649 0.0127 0.0343节点有功：0.7164 1.6300 0.8500 0.0000 -1.2500 -0.9000 -0.0000 -1.0000 -0.0000节点无功：0.2705 0.0665 -0.1086 0.0000 -0.5000 -0.3000 -0.0000 -0.3500 -0.0000修正后的节点导纳矩阵Y：Columns 1 through 50 -20.6944i 0 0 0 +17.3611i 00 0 -19.3333i 0 0 00 0 0 -20.3982i 0 00 +17.3611i 0 0 3.3074 -39.3089i -1.3652 +11.6041i0 0 0 -1.3652 +11.6041i 3.8716 -17.6627i0 0 0 -1.9422 +10.5107i 00 0 +16.0000i 0 0 -1.1876 + 5.9751i0 0 0 0 00 0 0 +17.0648i 0 0 Columns 6 through 90 0 0 00 0 +16.0000i 0 00 0 0 0 +17.0648i-1.9422 +10.5107i 0 0 00 -1.1876 + 5.9751i 0 04.1321 -16.0184i 0 0 -1.2820 +5.5882i0 2.8047 -35.4456i -1.6171 +13.6980i 00 -1.6171 +13.6980i 3.7323 -23.6669i -1.1551 + 9.7843i-1.2820 + 5.5882i 0 -1.1551 + 9.7843i 2.4371 -32.1539i节点阻抗矩阵Z的第4列：0.0463 + 0.1252i0.0329 + 0.0693i0.0316 + 0.0707i0.0552 + 0.1493i0.0589 + 0.1204i0.0562 + 0.1226i0.0397 + 0.0838i0.0416 + 0.0814i0.0378 + 0.0845i精确计算结果：短路电流：模值：6.4459相角：-71.9365节点电压模值：0.1831 0.5687 0.5427 0.0000 0.1466 0.1506 0.4537 0.4463 0.4495支路电流：i j Iij1 4 0.5779-3.1264i2 7 1.3702-1.4433i3 9 0.64294-1.4808i4 5 -0.77968+1.5248i4 6 -0.6411+1.477i5 7 -0.89528+1.6436i6 9 -0.73353+1.5487i7 8 0.50734+0.10234i8 9 0.062766+0.056451i近似计算结果：短路电流：模值：6.2838相角：-69.7198节点电压模值：0.1611 0.5214 0.5157 0.0000 0.1827 0.1675 0.4227 0.4348 0.4217五、程序流程图六、程序及输入文件input_data.xls 文件：powerflow_cal.m 文件：l=9;%支路数n=9;%节点数m=6;%PQ节点数Yn=zerosn;%初始化节点导纳矩阵Y DATA1=xlsread'input_data.xls';1; %计算节点导纳矩阵Yfor k=1:li=DATA1k;1;j=DATA1k;2;R=DATA1k;3;X=DATA1k;4;B2=DATA1k;5;Yni;i=Yni;i+1i*B2+1/R+1i*X; Ynj;j=Ynj;j+1i*B2+1/R+1i*X; Yni;j=Yni;j-1/R+1i*X;Ynj;i=Ynj;i-1/R+1i*X;enddisp'节点导纳矩阵Yn：';dispYn;G=realYn;B=imagYn;DATA2=xlsread'input_data.xls';2;P=zeros1;n;Q=zeros1;n;U=ones1;n;P2:n=DATA22:n;3;Q4:n=DATA24:n;4;U1:3=DATA21:3;5;%设置节点电压初值e1=DATA21;5;e2:n=1.0;f1:n=0.0;%设置迭代次数t=0;tmax=10;while t<=tmax%计算fxa1:n=0.0;c1:n=0.0;for i=2:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj;ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=2:ndeltaPi=Pi-ei*ai-fi*ci;endfor j=4:ndeltaQj=Qj-fj*aj+ej*cj;endfor k=2:3deltaU2k=Uk*Uk-ek*ek-fk*fk;endfx=deltaP2:n deltaQ4:n deltaU22:3';%计算雅克比矩阵Jfor i=2:nfor j=2:nif i~=jHi;j=-Gi;j*ei+Bi;j*fi;Ni;j=Bi;j*ei-Gi;j*fi;elseHi;j=-ai-Gi;i*ei+Bi;i*fi;Ni;j=-ci+Bi;i*ei-Gi;i*fi;endendendfor i=4:nfor j=2:nif i~=jMi;j=Bi;j*ei-Gi;j*fi;Li;j=Gi;j*ei+Bi;j*fi;elseMi;j=ci+Bi;i*ei-Gi;i*fi;Li;j=-ai+Gi;i*ei+Bi;i*fi;endendendfor i=2:3for j=2:nif i~=jRi;j=0;Si;j=0;elseRi;j=-2*ei;Si;j=-2*fi;endendendJ=H2:n;2:n N2:n;2:n;M4:n;2:n L4:n;2:n;R2:3;2:n S2:3;2:n;if maxabsfx<0.0001%输出结果break;else%求解修正方程获得dxdx=-J^-1*fx;dx=dx';e2:n=e2:n+dx1:n-1;f2:n=f2:n+dxn:2*n-1;t=t+1;endendif t>tmaxstr='潮流计算不收敛';dispstr;elsea1:n=0.0;c1:n=0.0;for i=1:nfor j=1:nai=ai+Gi;j*ej-Bi;j*fj; ci=ci+Gi;j*fj+Bi;j*ej;endendfor i=1:nUi=ei+1i*fi;ampi=absUi;argi=angleUi;Pi=ei*ai+fi*ci;Qi=fi*ai-ei*ci;enddisp'电压幅值：';dispamp;disp'电压相角：';disparg;disp'节点有功：';dispP;disp'节点无功：';dispQ;end%计算短路电流f=4;zf=0.0;%修正节点导纳矩阵Xd=DATA21:3;6;E=DATA21:3;7;for i=1:3Iii=Ei/1i*Xdi;endY=Yn;for i=1:3Yi;i=Yi;i+1/1i*Xdi;endfor j=4:nYj;j=Yj;j+-Pj+1i*Qj/Uj*Uj;enddisp'修正后的节点导纳矩阵Y：';Z=Y^-1;disp'节点阻抗矩阵Z的第4列：';dispZ:;4;%精确计算disp'精确计算结果：';U0=U;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';disparg;for i=1:nUi=U0i-Zi;f*If;amp=absU;enddisp'节点电压模值：';dispamp;disp'支路电流： ';str='i ''j '' Iij';dispstr;for k=1:li=DATA1k;1;j=DATA1k;2;r=DATA1k;3;x=DATA1k;4;z=r+1i*x;I=Ui-Uj/z;str=num2stri ' ' num2strj ' ' num2strI; dispstr;end%近似计算disp'近似计算结果：';U01:n=1.0;If=U0f/Zf;f+zf;amp=absIf;arg=atandimagIf/realIf;disp'短路电流：';disp'模值：';dispamp;disp'相角：';for i=1:nUi=U0i-Zi;f*If; amp=absU;enddisp'节点电压模值：'; dispamp;。
含集群风电的潮流计算的两种方法潮流计算是电力系统分析中的关键问题之一，主要用于估计输电网中各节点的电压、功率和相角等状态量，以及计算电网输电能力。
对于含集群风电的潮流计算，由于集群风电的特殊性，需采用一些特殊的方法来完成潮流计算。
下面将介绍两种常用的方法：1.静态等值法：静态等值法是一种将复杂的集群风电系统简化为等效节点和等效负荷的方法。
其基本思想是将集群风电系统中的大量风机等效为一个或多个节点，并根据风电系统与主网的连接关系确定等效节点的等值功率和等值阻抗。
等效节点的功率可以通过对风电机组的功率曲线进行相应的处理得到。
等效阻抗则可以通过对风电机组的等效阻抗进行计算得到。
在等影响因素确定之后，通过将这些等效节点与传统的负荷节点一起放入插值微分方程中进行潮流计算。
优点：（1）较好的处理了集群风电系统的复杂性；（2）能够较快地得到潮流计算结果。
缺点：（1）简化过程中可能会引入一定的误差；（2）对于大规模的集群风电系统，等效节点的选取比较困难。
2.迭代法：迭代法是一种通过迭代求解的方法来完成含集群风电的潮流计算。
其基本思想是通过迭代解算节点的电压和相角等状态量，直到满足潮流计算的收敛条件为止。
迭代法可分为高斯-赛德尔迭代法和牛顿-拉夫逊迭代法。
高斯-赛德尔迭代法通过按照节点顺序来更新节点的状态量，并将新的状态量作为旧状态量的近似值进行下一次迭代，直到达到收敛条件。
牛顿-拉夫逊迭代法是一种基于牛顿法的改进方法，通过构建雅可比矩阵对节点电压和相角进行一次迭代，并根据迭代结果修正节点的状态量，直到达到收敛条件。
优点：（1）准确度较高，能够更好地反映集群风电系统的实际情况；（2）对于大规模集群风电系统有较好的可行性。
缺点：（1）计算量较大，耗时较长；（2）可能会出现收敛问题，需要进行合理的参数选择和初始值设定。
总体而言，静态等值法适用于对集群风电系统的初步分析和初步评估，而迭代法适用于对集群风电系统进行深入研究，以及对系统进行详细的潮流计算。
含风电场的电力系统概率潮流计算摘要：由于对环保的关注，主要收获可再生能源（RES）的分布式能源（分布式能源）得到空前上升的关注。
这种类型的能源的天生不确定性增加电力系统中的不确定性，因此，就必须对系统性能进行概率分析。
此外，除了他们的不确定性，不确定参数具有相当水平的相关性。
两点估计法（2PEM）被公认为是适当的解决小规模甚至中等规模问题概率方法。
本文通过两点估计法计算概率潮流问题。
为了证明该方法的效果，用Mathpower14节点系统验证该方法。
然后，将得到的结果与蒙特卡罗模拟（MCS）的结果相比较。
关键字：概率潮流；两点法；风力发电引言最优概率潮流是电力市场中的重要工具，通过最优潮流模拟市场竞价过程，可获得交易量和节点电价等重要指标。
传统最优潮流研究大都基于确定性模型，即市场报价、负荷分布和元件参数等条件固定不变。
从宏观上看，一定时期内发电商报价和用户消费电能具有一定的确定性，但从微观角度来看，每个时段内发电商报价和用户消费的电能又会在各种因素影响下产生变化，这将引起交易量和节点电价的不断波动。
因此采用确定性模型进行最优潮流计算得到的结果，不能全面反映不确定因素对市场交易的影响。
计及发电报价、负荷分布中存在的不确定因素，采用概率最优潮流对市场交易进行模拟，能揭示出随机性和概率性后面隐藏的规律，为市场运营提供更多信息，降低交易风险，更好地引导市场交易的开展。
2．不确定模型负荷作为最显眼的不确定变量对电力系统运行起着至关重要的作用。
它的波动与时间，天气条件和电价等有关。
对于负载一种常见的做法到通过正态分布特定平均值和STD值，从历史数据获得的模型。
在这项研究中，负荷通过正态分布函数来模拟，平均值等于基本负载并且STD等于其平均值5％。
为了模拟风力发电的不确定性，一些节点被认为具有风电场和不确定的输出功率。
风速随着时间和地点的变化而变化和它的PDF遵循weibull分布。
因此，风速用weibull分布函数建模。
潮流计算，电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。
潮流计算是根据给定的电网结构、参数和发电机、负荷等元件的运行条件，确定电力系统各部分稳态运行状态参数的计算。
通常给定的运行条件有系统中各电源和负荷点的功率、枢纽点电压、平衡点的电压和相位角。
待求的运行状态参量包括电网各母线节点的电压幅值和相角，以及各支路的功率分布、网络的功率损耗等。
潮流计算程序，相关的原始数据数据数据输入格式如下：%B1是支路参数矩阵，第一列和第二列是起始节点编号和终点节点编号%第三列、第四列、第五列、第六列、第七列、第八列分别为：支路电阻、电抗、电导、电纳、变压器变比、是否有变压器（1为有、0为无）。
%B2为节点参数矩阵，其中第一列到第六列为节点编号；为节点类型；注入有功、注入无功、电压幅值、电压相位。
%“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。
n=input('请输入节点数:n=');n1=input('请输入支路数:n1=');isb=input('请输入平衡节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入支路参数:B1=');B2=input('请输入节点参数:B2=');Y=zeros(n);N=1;%建立节点导纳矩阵for i=1:n1if B1(i,8)==0p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4))-B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4))+0.5*B1(i,6)+B1(i,5);elsep=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/((B1(i,3)+B1(i,4))*B1(i,7)) -B1(i,5);Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4))+B1(i,5);Y(q,q)=Y(q,q)+1/(B1(i,7)^2*(B1(i,3)+B1(i,4)))+B1(i,5);endendYG=real(Y);B=imag(Y);PriS=zeros(2*n-2,1);ImbS=zeros(2*n-2,1);%创建PriS，用于存储初始功率参数h=0;j=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i, j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B 2(j,5));endendendfor i=1:nif i~=isb&B2(i,2)==2h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5))endendendPriSU3=zeros(n-h-1,2);%U3存储PV节点的初始电压t=0;for i=1:nif B2(i,2)==2t=t+1;U3(t,1)=B2(i,5);U3(t,2)=B2(i,6);endendU3%ImbS于存储有功功率、无功功率和电压幅值的不平衡量h=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=B2(i,4)-PriS(2*h,1);endendt=0;for i=1:nif i~=isb&B2(i,2)==2h=h+1;t=t+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=U3(t,1)^2+U3(t,2)^2-B2(i,5)^2-B2(i,6)^2;endendImbSI=zeros(n-1,1);%I，存储节点电流参数h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(PriS(2*h-1,1)-PriS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendIJacbi=zeros(2*n-2);%Jacbi(雅可比矩阵)h=0;k=0;for i=1:nif B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);elseJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbi%求解修正方程，获取节点电压的不平衡量ImbU=zeros(2*n-2,1);ImbU=inv(Jacbi)*ImbS;ImbU%修正节点电压j=0;for i=1:nif B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendfor i=1:nif B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendB2while abs(max(ImbU))>prPriS=zeros(2*n-2,1);h=0;j=0;for i=1:nif i~=isb&B2(i,2)==1h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i, j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B 2(j,5));endendendfor i=1:nif i~=isb&B2(i,2)==2h=h+1;for j=1:nPriS(2*h-1,1)=PriS(2*h-1,1)+B2(i,5)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))+B2(i,6)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5));PriS(2*h,1)=PriS(2*h,1)+B2(i,6)*(G(i,j)*B2(j,5)-B(i,j)*B2(j,6))-B2(i,5)*(G(i,j)*B2(j,6)+B(i,j)*B2(j,5))endendendPriS%创建ImbSh=0;for i=1:n %对PQ节点的处理if i~=isb&B2(i,2)==1h=h+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=B2(i,4)-PriS(2*h,1);endendt=0;for i=1:n %对PV节点的处理if i~=isb&B2(i,2)==2h=h+1;t=t+1;ImbS(2*h-1,1)=B2(i,3)-PriS(2*h-1,1);ImbS(2*h,1)=U3(t,1)^2+U3(t,2)^2-B2(i,5)^2-B2(i,6)^2;endendImbS%创建II=zeros(n-1,1);h=0;for i=1:nif i~=isbh=h+1;I(h,1)=(PriS(2*h-1,1)-PriS(2*h,1)*sqrt(-1))/conj(B2(i,5)+B2(i,6)*sqrt(-1));endendI%创建JacbiJacbi=zeros(2*n-2);h=0;k=0;for i=1:nif B2(i,2)==1h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k)+2*real(I(h,1));Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1)-2*imag(I(h,1));elseJacbi(2*h-1,2*k-1)=-B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)=G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=-Jacbi(2*h-1,2*k);Jacbi(2*h,2*k)=Jacbi(2*h-1,2*k-1);endif k==(n-1)k=0;endendendendendk=0;for i=1:nif B2(i,2)==2h=h+1;for j=1:nif j~=isbk=k+1;if i==jJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6)+imag(I(h,1));Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6)+real(I(h,1));Jacbi(2*h,2*k-1)=2*B2(i,6);Jacbi(2*h,2*k)=2*B2(i,5);elseJacbi(2*h-1,2*k-1)= -B(i,j)*B2(i,5)+G(i,j)*B2(i,6);Jacbi(2*h-1,2*k)= G(i,j)*B2(i,5)+B(i,j)*B2(i,6);Jacbi(2*h,2*k-1)=0;Jacbi(2*h,2*k)=0;endif k==(n-1)k=0;endendendendendJacbiImbU=zeros(2*n-2,1);ImbU=inv(Jacbi)*ImbS;ImbU%修正节点电压j=0;for i=1:nif B2(i,2)==1j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendfor i=1:nif B2(i,2)==2j=j+1;B2(i,5)=B2(i,5)+ImbU(2*j,1);B2(i,6)=B2(i,6)+ImbU(2*j-1,1);endendB2N=N+1; %迭代次数加1endN。
%本程序的功能是用牛拉法进行含风电场的电力系统潮流计算function s1=pf(a,B1,B2)n=a(1);%节点数nl=a(2);%支路数isb=a(3);%平衡节点号pr=a(4);%误差精度for i=1:nfor j=1:nG(i,j)=0;B(i,j)=0;endend%求导纳矩阵%B1为支路参数矩阵，其每一行格式为首节点号，末节点号，支路电导，支路电纳，首节点对地电纳，末节点对地电纳for i=1:nlp=B1(i,1);q=B1(i,2)G(p,q)=G(p,q)-B1(i,3);B(p,q)=B(p,q)- B1(i,4);G(q,p)=G(p,q);B(q,p)=B(p,q);G(p,p)=G(p,p)+ B1(i,3);B(p,p)=B(p,p)+B1(i,4);G(q,q)=G(q,q)+ B1(i,3);B(q,q)=B(q,q)+B1(i,4);endfor i=1:nB(i,i)=B(i,i)+B2(i,5);end%求导纳矩阵%B2为节点参数矩阵，每一行格式为节点注入有功，注入无功，电压实部，电压虚部，对地电纳，节点类型%节点类型：1为平衡节点，2为PQ节点，3为PV节点，4为风电场节点%Bf为风电场参数，格式为：有功功率，定子电抗，转子漏抗，转子电阻，励磁电抗for i=1:nif B2(i,6)==4p(i)=Bf(1);a1=2*p(i)^2*(Bf(2)+Bf(3))^2;a2=-Bf(4)^3*Bf(5);a3=Bf(4)^2*Bf(5);a4=Bf(4)^2;a5=4*p(i)^2*(Bf(2)+Bf(3))^2*Bf(4)^2;elseP(i)=B2(i,1);Q(i)=B2(i,2);ende(i)= B2(i,3);f(i)=B2(i,4);V(i)=sqrt(e(i)^2+f(i)^2);endICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:nif i~=isbC(i)=0;D(i)=0;for j1=1:nC(i)= C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);D(i)= D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);endP1=C(i)*e(i)+f(i)*D(i);Q1=f(i)*C(i)-D(i)*e(i);V2=e(i)^2+f(i)^2;if B2(i,6)==2DP=P(i)-P1;DQ=Q(i)-Q1;for j1=1:nif j1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*e(i)-G(i,j1)*f(i);X3=X2;X4=-X1;p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1==i&j1~=isbX1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;endendelseif B2(i,6)==3DP=P(i)-P1;DV=V(i)^2-V2;for j1=1:nif j1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*e(i)-G(i,j1)*f(i);X5=0;X6=0;p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;elseif j1==i&j1~=isbX1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);X5=-2*e(i);X6=-2*f(i);p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6;J(m,q)=X2;endendelseDP=P(i)-P1;for j1=1:nif j1~=isb&j1~=iX1=-G(i,j1)*e(i)-B(i,j1)*f(i);X2=B(i,j1)*e(i)-G(i,j1)*f(i);X3=X2;X4=-X1;p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;elseif j1==i&j1~=isbX1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);x3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);x4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);x5=2*a1*(a2+a3*a4*V2^2/(sqrt(V2^4*a4-a5)))/(a2*V2^2+a3*(sqrt(V2^4*a4-a5)))^2;X3=x3-e(i)*x5;X4=x4-f(i)*x5;p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;m=p+1;J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X4;J(m,q)=X2;endendendendend%求雅克比矩阵for k=3:N0k1=k+1;N1=N;for k2=k1:N1J(k,k2)=J(k,k2)./J(k,k);endJ(k,k)=1;if k~=3;k4=k-1;for k3=3:k4for k2=k1:N1J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endif k==N0,break;endfor k3=k1:N0for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endelsefor k3=k1:N0for k2=k1:N1J(k3,k2)= J(k3,k2)-J(k3,k)*J(k,k2);endJ(k3,k)=0;endendendfor k=3:2:N0-1L=(k+1)./2;e(L)=e(L)-J(k,N);k1=k+1;f(L)=f(L)-J(k1,N);endfor k=3:N0DET=abs(J(k,N));if DET>=prIT2=IT2+1;endendICT2(a)=IT2;ICT1=ICT1+1;end%用高斯消去法解w=-J*Vfid1=fopen('out1.txt','wt')fprintf(fid1,‘各节点的实际电压标幺值E和电压大小V为：\n’);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);fprintf(fid1,'E(%d)=%8.5f+j%8.5f, V(%d)=%8.5f\n',k,e(k),f(k),k,V(k)); endfor p=1:nif p==isbfprintf(fid1,'平衡节点的功率S为:\n');C(p)=0;D(p)=0;for q=1:nC(p)=C(p)+G(p,q)*e(q)-B(p,q)*f(q);D(p)=D(p)-G(p,q)*f(q)-B(p,q)*e(q);endSp(p)=e(p)*C(p)-f(p)*D(p);Sq(p)=(e(p)*D(p)+f(p)*C(p));fprintf(fid1,'S=%8.5f+j%8.5f\n',Sp(p),Sq(p));endendfprintf(fid1,'各条支路的首端功率为:\n');for i=1:nlp=B1(i,1);q=B1(i,2);sip(i)=-e(p)*B1(i,5)+(e(q)-e(p))*G(p,q)-(f(q)-f(p))*B(p,q);siq(i)=(f(p)*B1(i,5)+(e(p)-e(q))*B(p,q)-(f(q)-f(p))*G(p,q));Sip(i)=e(p)*sip(i)-f(p)*siq(i);Siq(i)=e(p)*siq(i)+f(p)*sip(i);fprintf(fid1,'S(%d,%d)=%8.5f+j%8.5f\n',p,q,Sip(i),Siq(i));endfprintf(fid1,'¸各条支路的末端功率为:\n');for i=1:nlp=B1(i,1);q=B1(i,2);sjp(i)=-e(q)*B1(i,6)+(e(p)-e(q))*G(p,q)-(f(p)-f(q))*B(p,q);sjq(i)=(f(q)*B1(i,6)+(e(q)-e(p))*B(p,q)-(f(p)-f(q))*G(p,q));Sjp(i)=e(q)*sjp(i)-f(q)*sjq(i);Sjq(i)=e(q)*sjq(i)+f(q)*sjp(i);fprintf(fid1,'S(%d,%d)=%8.5f+j%8.5f\n',q,p,Sjp(i),Sjq(i));endfprintf(fid1,'各条支路的功率损耗为:\n');for i=1:nlDSp(i)=Sip(i)+Sjp(i);DSq(i)=Siq(i)+Sjq(i);fprintf(fid1,'DS(%d)=%8.5f+j%8.5f\n',i,DSp(i),DSq(i)); endfclose(fid1);。
