基于MATLAB的潮流计算源程序代码(优.选)

合集下载

基于MATLAB的潮流计算源程序代码(优.选)

%*************************电力系统直角坐标系下的牛顿拉夫逊法潮流计算**********clearclcload E:\data\IEEE014_Node.txtNode=IEEE014_Node;weishu=size(Node);nnum=weishu(1,1); %节点总数load E:\data\IEEE014_Branch.txtbranch=IEEE014_Branch;bwei=size(branch);bnum=bwei(1,1); %支路总数Y=(zeros(nnum));Sj=100;%********************************节点导纳矩阵*******************************for m=1:bnum;s=branch(m,1); %首节点e=branch(m,2); %末节点R=branch(m,3); %支路电阻X=branch(m,4); %支路电抗B=branch(m,5); %支路对地电纳k=branch(m,6);if k==0 %无变压器支路情形Y(s,e)=-1/(R+j*X); %互导纳Y(e,s)=Y(s,e);endif k~=0 %有变压器支路情形Y(s,e)=-(1/((R+j*X)*k));Y(e,s)=Y(s,e);Y(s,s)=-(1-k)/((R+j*X)*k^2);Y(e,e)=-(k-1)/((R+j*X)*k); %对地导纳endY(s,s)=Y(s,s)-j*B/2;Y(e,e)=Y(e,e)-j*B/2; %自导纳的计算情形endfor t=1:nnum;Y(t,t)=-sum(Y(t,:))+Node(t,12)+j*Node(t,13);%求支路自导纳endG=real(Y); %电导B=imag(Y); %电纳%******************节点分类************************************* *pq=0; pv=0; blancenode=0;pqnode=zeros(1,nnum);pvnode=zeros(1,nnum);for m=1:nnum;if Node(m,2)==3blancenode=m; %平衡节点编号else if Node(m,2)==0pq=pq+1;pqnode(1,pq)=m; %PQ 节点编号else if Node(m,2)==2pv=pv+1;pvnode(1,pv)=m; %PV 节点编号endendendend%*****************************设置电压初值********************************** Uoriginal=zeros(1,nnum); %对各节点电压矩阵初始化for n=1:nnumUoriginal(1,n)=Node(n,9); %对各点电压赋初值if Node(n,9)==0;Uoriginal(1,n)=1; %该节点为非PV节点时，将电压值赋为1endendPresion=input('请输入误差精度要求：Presion=');disp('该电力系统节点数:');disp(nnum);xiumax=0.1;counter=0;while xiumax>Presion%****************************计算不平衡量***********************************e=real(Uoriginal); %取初始电压的实部f=imag(Uoriginal); %取初始电压的虚部deta=zeros(2*pq+2*pv,1); %构造储存功率变化量的列矩阵n=1;for m=1:pq;Pi=0;Qi=0;for t=1:nnum;Pi=Pi+e(1,pqnode(1,m))*(G(pqnode(1,m),t)*e (1,t)-B(pqnode(1,m),t)*f(1,t))+f(1,pqnode(1,m ))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t) *e(1,t));%计算该PQ节点的负荷有功Qi=Qi+f(1,pqnode(1,m))*(G(pqnode(1,m),t)*e (1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1,m ))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m),t) *e(1,t));%计算该PQ节点的负荷无功endS1(1,pqnode(1,m))=Pi+j*Qi;P=(Node(pqnode(1,m),7)-Node(pqnode(1,m), 5))/Sj-Pi;%计算该PQ节点的实际有功功率deta(n,1)=P;%在该列向量中储存有功功率n=n+1;Q=(Node(pqnode(1,m),8)-Node(pqnode(1,m),6))/Sj-Qi;%计算该PQ节点的实际无功功率deta(n,1)=Q;%在该列向量中储存无功功率n=n+1;endfor m=1:pv;Pv=0; Qv=0;for t=1:nnum;Pv=Pv+e(1,pvnode(1,m))*(G(pvnode(1,m),t)* e(1,t)-B(pvnode(1,m),t)*f(1,t))+f(1,pvnode(1, m))*(G(pvnode(1,m),t)*f(1,t)+B(pvnode(1,m), t)*e(1,t));%计算该PV节点的负荷有功Ui=e(1,pvnode(1,m))^2+f(1,pvnode(1,m))^2; %计算该节点的负荷电压值Qv=Qv+f(1,pqnode(1,m))*(G(pqnode(1,m),t)* e(1,t)-B(pqnode(1,m),t)*f(1,t))-e(1,pqnode(1, m))*(G(pqnode(1,m),t)*f(1,t)+B(pqnode(1,m), t)*e(1,t));endS1(1,pvnode(1,m))=Pv+j*Qv;P=(Node(pvnode(1,m),7)-Node(pvnode(1,m), 5))/Sj-Pv; %计算该节点的实际有功功率deta(n,1)=P; %储存该有功功率n=n+1;U=Node(pvnode(1,m),3)^2-Ui; %计算电压变化量deta(n,1)=U; %储存该电压变化量n=n+1;enddeta;cerate=zeros(pq+pv,1);for k=1:pqcerate(k,1)=pqnode(1,k);endfor v=1:pvcerate(pq+v,1)=pvnode(1,v);end%******************************雅克比矩阵****************************** Jacob=ones(2*nnum-2);L=0;J=0;H=0;N=0; R=0;S=0;n=1;k=1;form=1:pq; %m表示雅克比矩阵中pq节点的行数for u=1:pq+pv; %u 表示雅克比矩阵中pq节点的列数t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的个数if pqnode(1,m)~=t %非对角元素的情况H=G(pqnode(1,m),t)*f(1,pqnode(1,m))-B(pqn ode(1,m),t)*e(1,pqnode(1,m));N=G(pqnode(1,m),t)*e(1,pqnode(1,m))+B(pq node(1,m),t)*f(1,pqnode(1,m));L=H;J=-N;elseifpqnode(1,m)==t %对角线元素时的情况I=0;for g=1:nnumI=Y(t,g)*Uoriginal(1,g)+I; %计算节点的注入电流endaii=real(I);bii=imag(I);H=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode (1,m))+bii;N=G(t,t)*e(1,pqnode(1,m))+B(t,t)*f(1,pqnode (1,m))+aii;L=-B(t,t)*e(1,pqnode(1,m))+G(t,t)*f(1,pqnode (1,m))-bii; J=-G(t,t)*e(1,pqnode(1,m))-B(t,t)*f(1,pqnode( 1,m))+aii;endendJacob(n,k)=H;k=k+1;Jacob(n,k)=N;k=k-1;n=n+1;Jacob(n,k)=J;k=k+1;Jacob(n,k)=L;n=n-1; k=k+1; %按照雅克比矩阵的排列规则排列pq节点的雅克比元素endk=1; n=2*m+1; %将光标定位于下一个待排列PQ节点元素的第一个位置endn=2*pq+1; k=1; %定位于PV节点的第一个位置处for m=1:pv;for u=1:pq+pv;t=cerate(u,1); %t为中间变量，用来标记雅克比矩阵中指定元素的位置if pvnode(1,m)~=t %非对角线元素情况H=G(pvnode(1,m),t)*f(1,pvnode(1,m))-B(pvno de(1,m),t)*e(1,pvnode(1,m));N=G(pvnode(1,m),t)*e(1,pvnode(1,m))+B(pvn ode(1,m),t)*f(1,pvnode(1,m));R=0; S=0;endifpvnode(1,m)==t %对角线元素情况I=0;for g=1:nnumI=Y(t,g)*Uoriginal(1,g)+I; %计算PV节点的注入电流endaii=real(I);bii=imag(I);H=-B(t,t)*e(1,pvnode(1,m))+G(t,t)*f(1,pvnode (1,m))+bii;N=G(t,t)*e(1,pvnode(1,m))+B(t,t)*f(1,pvnode( 1,m))+aii;R=2*f(1,pvnode(1,m));S=2*e(1,pvnode(1,m));endJacob(n,k)=H; k=k+1;Jacob(n,k)=N; k=k-1;n=n+1;Jacob(n,k)=R; k=k+1;Jacob(n,k)=S;n=n-1;k=k+1; %按照雅克比矩阵的排列规则排列PV节点的雅克比元素endk=1;n=n+2; %定位于下一个待排列PV节点的雅克比元素第一个位置end%*************************电压变化量化的计算与存储************************************ Detau=inv(Jacob)*deta; %构建电压的变化量的列向量f=zeros(1,nnum); %给电压实部赋初值0e=zeros(1,nnum); %给电压虚部赋初值0for p=1:(pq+pv);f(1,cerate(p,1))=j*Detau(2*p-1,1);%将电压变量的奇数行赋值给fe(1,cerate(p,1))=Detau(2*p,1); %将电压变量的偶数行赋值给eendt=e+f;xiumax=abs(Detau(1,1)); %将电压变化量的第一个元素赋值给最大允许误差for n=2:2*nnum-2;if abs(Detau(n,1))>xiumaxxiumax=abs(Detau(n,1)); %找出最大的电压误差endendUoriginal=Uoriginal+t; %迭代修正后的电压值counter=1+counter; %统计迭代次数enddisp('迭代次数counter:');disp(counter);%**************************平衡节点功率及显示**********************************m=blancenode;t=0;for n=1:nnum;t=t+(G(m,n)-j*B(m,n))*(real(Uoriginal(1,n))-j*i mag(Uoriginal(1,n)));endS1(1,m)=Uoriginal(1,m)*t;%**************************直角坐标下各节点电压及显示****************************U=zeros(1,nnum);for n=1:nnumUi(n,1)=Node(n,1);U(1,n)=real(Uoriginal(1,n))+i*imag(Uoriginal(1 ,n)); %将电压值由极坐标转化为直角坐标形式Ui(n,2)=U(1,n);Ui(n,3)=S1(1,n);enddisp('各节点电压直角坐标形式及节点注入功率:');disp(' 节点号节点电压值节点注入功率');disp(Ui);disp('修正电压的最大误差:')disp(xiumax);%**************************功率损耗************************************* **for m=1:bnum %支路功率及损耗startnode=branch(m,1);endnode=branch(m,2); %终止节点y=sum(Y,2);Sij=Uoriginal(1,startnode)*((conj(Uoriginal(1,s tartnode))*conj(y(startnode,1)))+(conj(Uorigi nal(1,startnode))-conj(Uoriginal(1,endnode))) *conj(-Y(startnode,endnode)));Sji=Uoriginal(1,endnode)*((conj(Uoriginal(1,e ndnode))*conj(y(endnode,1)))+(conj(Uoriginal (1,endnode))-conj(Uoriginal(1,startnode)))*co nj(-Y(endnode,startnode)));S(m,1)=startnode; %始节点S(m,2)=endnode; %末节点S(m,3)=Sij; %始节点流入的功率S(m,4)=Sji; %末节点流入的功率S(m,5)=Sij+Sji; %线路损耗enddisp('支路功率及损耗：');disp('节点号(I) 节点号(J) 支路功率（I-J）支路功率（J-I）线路损耗（delta_S)')disp(S); 最新文件---------------- 仅供参考--------------------已改成word文本--------------------- 方便更改。

电力系统潮流计算matlab程序

电力系统潮流计算matlab程序电力系统潮流计算是电力系统运行和规划中的重要环节，它用于计算电力系统中各节点的电压、功率和电流等参数。

随着电力系统规模的不断扩大和复杂性的增加，传统的手工计算方法已经无法满足需求，因此，利用计算机编程进行潮流计算成为了一种必要的选择。

Matlab是一种功能强大的科学计算软件，它提供了丰富的数学函数和工具箱，可以方便地进行电力系统潮流计算。

下面我将介绍一下如何使用Matlab编写电力系统潮流计算程序。

首先，我们需要建立电力系统的节点模型。

节点模型是电力系统中各节点的电压、功率和电流等参数的数学表示。

在Matlab中，我们可以使用矩阵来表示节点模型。

假设电力系统有n个节点，我们可以定义一个n×n的复数矩阵Y来表示节点之间的导纳关系，其中Y(i,j)表示节点i和节点j之间的导纳。

同时，我们还需要定义一个n×1的复数向量V来表示各节点的电压，其中V(i)表示节点i的电压。

接下来，我们需要编写潮流计算的主程序。

主程序的主要功能是根据节点模型和潮流计算算法，计算出各节点的电压、功率和电流等参数。

在Matlab中，我们可以使用循环语句和矩阵运算来实现潮流计算。

具体的计算过程可以参考电力系统潮流计算的算法。

在编写主程序之前，我们还需要定义一些输入参数，如电力系统的节点数、发电机节点和负荷节点等。

这些参数可以通过用户输入或者读取文件的方式获取。

同时，我们还需要定义一些输出参数，如各节点的电压、功率和电流等。

这些参数可以通过矩阵运算和循环语句计算得到，并输出到文件或者显示在屏幕上。

最后，我们需要进行程序的测试和调试。

可以通过输入一些测试数据，运行程序并检查输出结果是否正确。

如果发现程序有错误或者结果不准确，可以通过调试工具和打印调试信息的方式进行调试。

总之，利用Matlab编写电力系统潮流计算程序可以提高计算效率和准确性，为电力系统的运行和规划提供有力的支持。

当然，编写一个完整的潮流计算程序需要考虑很多细节和特殊情况，这需要有一定的电力系统和编程知识。

潮流计算matlab程序

clear;%各节点参数：节点编号，类型，电压幅值，电压相位，注入有功，注入无功%类型：1＝PQ节点，2＝PV节点，3＝平衡节点%本程序中将最后一个节点设为平衡节点R_1=[1 1 1.0 0 0.2 0.2j;2 1 1.0 0 -0.45 -0.15j;3 1 1.0 0 -0.45 -0.05j;4 1 1.0 0 -0.6 -0.1j;5 3 1.0 0 0 0];%支路号首端节点末端节点支路导纳R_2=[1 5 2 1.25-3.75j;2 23 10.00-30.00j;3 34 1.25-3.75j;4 1 4 2.50-7.50j;5 1 5 5.00-15.00j;6 1 2 1.667-5.00j];n=5;L=6;%需要改变的到此为止i=0;j=0;a=0;precision=1;k=0;Y=zeros(n,n);u=zeros(1,n);delt=zeros(1,n);P=zeros(1,n);Q=zeros(1,n);G=[];B=[];PP=[];uu=[];U=[];dp=[];dq=[];for a=1:Li=R_2(a,2);j=R_2(a,3);Y(i,j)=-R_2(a,4);Y(j,i)=Y(i,j);endfor a=1:nfor b=1:nif a~=bY(a,a)=Y(a,a)+Y(a,b);endendendfor i=1:nfor j=1:nif i==jY(i,j)=-Y(i,j);endendendY %形成导纳矩阵for i=1:nfor j=1:nG(i,j)=real(Y(i,j));B(i,j)=imag(Y(i,j));endendfor a=1:nu(a)=R_1(a,3);P(a)=R_1(a,5);Q(a)=R_1(a,6);delt(a)=R_1(a,4);endwhile precision>0.0001 %判断是否满足精度要求for a=1:n-1for b=1:npt(b)=u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));qt(b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b))); endpt,qtpi(a)=sum(pt);qi(a)=sum(qt); %计算PQ节点的注入功率dp(a)=P(a)-pi(a);dq(a)=Q(a)-qi(a); %计算PQ节点的功率不平衡量endfor a=1:n-1for b=1:n-1if a==bH(a,a)=-qi(a)-u(a)^2*B(a,a); N(a,a)=pi(a)+u(a)^2*G(a,a);J(a,a)=pi(a)-u(a)^2*G(a,a); L(a,a)=qi(a)-u(a)^2*B(a,a);JJ(2*a-1,2*a-1)=H(a,a); JJ(2*a-1,2*a)=N(a,a);JJ(2*a,2*a-1)=J(a,a); JJ(2*a,2*a)=L(a,a);elseH(a,b)=u(a)*u(b)*(G(a,b)*sin(delt(a)-delt(b))-B(a,b)*cos(delt(a)-delt(b)));J(a,b)=-u(a)*u(b)*(G(a,b)*cos(delt(a)-delt(b))+B(a,b)*sin(delt(a)-delt(b)));N(a,b)=-J(a,b);L(a,b)=H(a,b);JJ(2*a-1,2*b-1)=H(a,b);JJ(2*a-1,2*b)=N(a,b);JJ(2*a,2*b-1)=J(a,b); JJ(2*a,2*b)=L(a,b);endendend %计算jocbi各项,并放入统一矩阵JJ中，对JJ下标统一编号JJfor a=1:n-1PP(2*a-1)=dp(a);PP(2*a)=dq(a);end %按统一矩阵形成功率不平衡uu=inv(JJ)*PP';precision=max(abs(uu)); %判断是否收敛for b=1:n-1delt(b)=delt(b)+uu(2*b-1);u(b)=u(b)+uu(2*b)*u(b); %将结果分解为电压幅值和角度end %求解修正方程，得电压幅值变化量（标幺值）和角度变化量k=k+1;endfor a=1:nU(a)=u(a)*(cos(delt(a))+j*sin(delt(a)));endfor b=1:nI(b)=Y(n,b)*U(b);%求平衡节点的注入电流endS5=U(n)*sum(conj(I))%求平衡节点的注入功率for a=1:nfor b=1:nS(a,b)=U(a)*(conj(U(a))-conj(U(b)))*conj(-Y(a,b));endend %求节点i,j节点之间的功率，方向为由i指向j,S %显示支路功率。

基于MATLAB进行潮流计算

基于MATLAB进行潮流计算本文介绍了基于MATLAB软件的潮流计算方法。

电力系统潮流计算方法分为手算潮流和计算机潮流计算两类。

手算潮流主要适用于规模较小的辐射型电力潮流计算，而计算机潮流计算有两种途径：编程实现网络方程的迭代求解和借助电力系统分析仿真软件搭建系统模型完成潮流计算。

MATLAB具有强大的矩阵运算功能和电力系统仿真平台，可以为实现潮流计算提供更便捷的手段。

本文采用极坐标形式牛顿─拉夫逊法进行潮流计算，为其他形式的潮流计算提供借鉴。

Abstract: The power flow n method can be divided into two categories: hand n of tidal current and computer power flow XXX simplified equivalent circuits。

making it XXX: programming XXX ns。

or using power system XXX system model for power flow n。

MATLAB are has strong matrix ns and its power system XXX-Raphson method of power flow n in polar coordinates with MATLAB are。

and can serve as a reference for other forms of power flow n.1.电力系统中的牛顿法潮流计算是一种常用的电力系统分析方法。

该方法基于节点电压的相等条件和潮流方程的等式条件，通过迭代求解电压和相位的不平衡量，最终得到各节点的电压、相位和功率等参数。

2.牛顿法潮流计算的步骤包括输入系统原始数据、形成节点导纳矩阵、给定各节点电压初值、计算功率偏差向量、判断收敛条件、计算雅克比矩阵、解修正方程、计算节点电压和相位的修正值、迭代计算直至满足收敛条件、计算各节点功率等参数并输出计算结果。

matlab电力系统潮流计算程序

matlab电力系统潮流计算程序电力系统潮流计算是电力系统分析的关键步骤之一，用于确定电力系统各节点的电压和相角分布。

以下是一个简单的MATLAB电力系统潮流计算的基本步骤和代码示例：1.定义电力系统参数：-定义系统节点数量、支路数据、发电机数据、负荷数据等电力系统参数。

```matlab%电力系统参数busdata=[1,1.05,0,0,0,0,0,0;2,1.02,0,0,0,0,0,0;%...其他节点数据];linedata=[1,2,0.02,0.06,0.03;%...其他支路数据];gendata=[1,2,100,0,999,1.05,0.95;%...其他发电机数据];loaddata=[1,50,20;%...其他负荷数据];```2.构建潮流计算矩阵：-利用节点支路导纳、节点负荷和发电机功率等信息构建潮流计算的阻抗矩阵。

```matlabYbus=buildYbus(busdata,linedata);```3.迭代求解潮流方程：-利用迭代算法（如牛顿-拉夫森法）求解潮流方程，更新节点电压和相角。

```matlab[V,delta]=powerflow(Ybus,gendata,loaddata,busdata);```4.结果分析和可视化：-分析计算结果，可视化电压和相角分布。

```matlabplotVoltageProfile(busdata,V,delta);```这只是一个简单的潮流计算示例。

具体的程序实现可能涉及更复杂的算法和工程细节，取决于电力系统的复杂性和精确性要求。

您可能需要根据实际情况和数据格式进行调整和改进。

在实际工程中，也可以考虑使用专业的电力系统仿真软件。

潮流计算MATLAB 粗略程序

%========================================================================== %========================================================================== %========================================================================== %潮流计算MATLAB 粗略程序 C.zhou 2009.3.27%========================================================================== %========================================================================== %========================================================================== %creat a new_datat=0;s=0;r=0;w=0;number=input('How many node are there=');% Convert Pq to a new arrayfor ii=1:numberif data(ii,4)==1t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;s=s+1; %record the number of the PQ node end;end;%Convert pv to a new arrayfor ii=1:numberif data(ii,4)==2t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;r=r+1; %record the number of the PV node end;end;%Convert set_v to a new arrayfor ii=1:numberif data(ii,4)==3t=t+1;for jj=1:14new_data1(t,jj)=data(ii,jj);end;a(1,t)=ii;w=w+1;end;end;%creat a new_data2[x,y]=size(data2)for ii=1:xfor jj=1:2for mm=1:numberif data2(ii,jj)==a(1,mm)new_data2(ii,jj)=mm;end;end;end;end;for ii=1:xfor jj=3:14new_data2(ii,jj)=data2(ii,jj);end;end;%creat a YY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,number);for ii=1:x% for jj=1:14iii=new_data2(ii,1);jjj=new_data2(ii,2);if new_data2(ii,5)==2sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6 ))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;Y(iii,jjj)=-sub./new_data2(ii,14);YY(iii,jjj)=sub./new_data2(ii,14);Y(jjj,iii)=-sub/new_data2(ii,14);YY(jjj,iii)=sub./new_data2(ii,14);yy(iii,jjj)=(1.-new_data2(ii,14))./(new_data2(ii,14).*new_data2(ii,14)).*sub;yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14)).*sub;elseY(iii,jjj)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2 (ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data 2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2 (ii,6))+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6))*i;yy(iii,jjj)=new_data2(ii,8)./2.*i;yy(jjj,iii)=new_data2(ii,8)./2.*i;end;%end;end;for iii=1:numberY(iii,iii)=0;end;%for ii=1:x% for jj=1:14for iii=1:numberfor jj=1:number% if iii~=jjY(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);% end;end;end;%creat B, Gfor ii=1:numberfor jj=1:numberG(ii,jj)= real(Y(ii,jj));B(ii,jj)= imag(Y(ii,jj));end;end;%creat Initial_P Initial_Q Initial_Vfor ii=1:(s+r)set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7))./100;end;for ii=1:s;set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8))./100;end;for ii=1:rset_V(ii,1)=new_data1(ii+s,12).*new_data1(ii+s,12);%try to modify for sike of correcting end;Initial_p_q_v=[set_P;set_Q;set_V];disp(Initial_p_q_v);%creat Initial_e,Initial_ffor ii=1:number-1e(ii,1)=1;f(ii,1)=0.0;%change f to test used to be 1.0end;e(number,1)=new_data1(number,12);f(number,1)=0;% e(64,1)=0.88;%test 118ieee% f(64,1)=0.39395826829394;% f(14,1)=0;% e(10,1)=1.045;%e(11,1)=1.01;%e(12,1)=1.07;%e(13,1)=1.09;%//////////////////////////////////////////////////////////////////////////%/////////////////////////////////////////////////////////////////////////%//////////////////////////////////////////////////////////////////////////%//////////////////////////////////////////////////////////////////////////% Start NEWTOWN CALULATIONfor try_time=1:25%Creat every node consume P Q and Un=s;m=r;for ii=1:(n+m)sum1=0;for jj=1:(n+m+1)sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));end;p(ii,1)=sum1;end;for ii=1:nsum2=0;for jj=1:(n+m+1)sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1))-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1));end;q(ii,1)=sum2;end;disp('q=');disp(q);u=zeros((n+m),1);for ii=(n+1):(n+m)u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;for ii=n+1:(n+m)extra_u((ii-n),1)=u(ii,1);end;disp('extra_u=');disp(extra_u);sum=[p;q;extra_u];disp(sum)disp(s);disp(p);%creat Jacobiandisp(n);disp(m);for ii=1:(n+m)for jj=1:(n+m)if (ii~=jj)PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);%TEST+1end;end;end;copy=3.14159;disp('================copy================')for ii=1:nfor jj=1:m+nif (ii~=jj)QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1end;end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros((n+m),(n+m));UF=zeros((n+m),(n+m));for ii=n+1:n+mfor jj=1:(n+m)if (ii~=jj)UE(ii,jj)=0;UF(ii,jj)=0;elsess=0;qq=0;for num=1:(n+m+1)ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);end;UF(ii,jj)=-2.*f(ii,1);UE(ii,jj)=-2.*e(ii,1);end;end;end;for ii=(n+1):(n+m)for jj=1:(n+m)extra_UE((ii-n),jj)=UE(ii,jj);extra_UF((ii-n),jj)=UF(ii,jj);end;end;%disp('extra_UE');%disp(extra_UE);%disp('extra_Uf');%disp(extra_UF);Jacobian=[PF,PE;QF,QE;extra_UF,extra_UE];%disp('Jacobian=');%disp(Jacobian);%creat substract resultsubstract_result=Initial_p_q_v-sum;%disp('substract_result');%disp(substract_result);%calculate delta_f_edelta_f_e=-inv(Jacobian)*substract_result; %disp(delta_f_e);for ii=1:number-1;f(ii,1)=f(ii,1)+delta_f_e(ii,1);e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1); end;if max(substract_result)<1e-4break;end ;end;%disp('substract_result');%disp(substract_result);%disp('e=');%disp(e);%disp('f=');%disp(f);for ii=1:numberuuu(ii,1)= e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1));end;for ii=1:numberfor jj=1:numberif ii==a(1,jj)Old_Uresult(ii,1)=U_RESULT(jj,1)end;end;end;for ii=1:numberOld_Uresult(ii,2)=ii;end;%disp('U_result');%disp(U_RESULT);disp('====================================='); disp('The last result is :')disp('===========U===================BUS-NO.'); disp('U=')disp(Old_Uresult);%calculate the anglePI=3.141592for ii=1:numberAngle(ii,1)=atan(f(ii,1)./e(ii,1))./PI*180;end;for ii=1:numberfor jj=1:numberif ii==a(1,jj)Old_Angle(ii,1)=Angle(jj,1);Old_Angle(ii,2)=ii;end;end;end;disp('=======Angle===================BUS-NO.'); disp('Angle=');disp(Old_Angle);disp('=====Try-times=======================') disp('Try-times=')disp(try_time);%disp('p====================');%disp(p);% for jj=1:number% if a(1,jj)==118% man=jj% end;%end;%disp('man=========');%disp(man)sum4=0;for jj=1:numberY_conj(number,jj)=conj(Y(number,jj));sum4=sum4+Y_conj(number,jj).*(e(jj,1)-f(jj,1)*i);end;%sum4=sum4*e(number,1);disp('===============Balance P Q=========BUS-NO');%disp(sum4);Blance_Q(1,1)=imag(sum4)*100;Blance_Q(1,2)=a(1,number);Blance_P(1,1)=real(sum4)*100;Blance_P(1,2)=a(1,number);disp('Q Of the Balance node= ');disp(Blance_Q);disp('P Of the Balance node= ')disp(Blance_P);disp('=================================BUS-NO');%calculate the Q of the P-V nodeQ_PV_node=zeros(number,2);Y_conj=conj(Y);for ii=(s+1):(s+r)for jj=1:numberQ_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)*i).*(Y_conj(ii,jj).*(e(jj,1)-f(jj,1)*i));end;end;for ii=(s+1):(s+r);Q_PV_node(ii,1)=Q_PV_node(ii,1).*100+new_data1(ii,8)*i;end;disp('This program is from /breadwinner') ;for ii=1:numberold_number=a(1,ii);Q_PV_node_old(old_number,1)=Q_PV_node(ii,1);end;for ii=1:numberQ_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));end;for ii=1:numberQ_PV_node_old(ii,2)=ii;end;disp('Q gen=');disp(Q_PV_node_old);。

牛顿—拉夫逊法潮流计算MATLAB程序

牛顿—拉夫逊法潮流计算程序By Yuluo%牛顿--拉夫逊法进行潮流计算n=input('请输入节点数：n=');n1=input('请输入支路数：n1=');isb=input('请输入平衡母线节点号：isb=');pr=input('请输入误差精度：pr=');B1=input('请输入由支路参数形成的矩阵：B1=');B2=input('请输入各节点参数形成的矩阵：B2=');X=input('请输入由节点参数形成的矩阵：X=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=seros(1,n);O=zeros(1,n);S1=zeros(n1);for i=1:nif X(i,2)~=0;p=X(i,1);Y(p,p)=1./X(i,2);endendfor i=1:n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5);Y(p,q)=Y(p,q);Y(p,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;end %求导纳矩阵G=real(Y);B=imag(Y);for i=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1:nS(i)=B2(i,1)-B2(i,2);B(i,i)=B(i,i)+B2(i,5);endP=rea(S);Q=imag(S);ICT1=0;IT2=1;NO=2*n;N=NO+1;a=0;while IT2~=0IT2=0;a=a+1;for i=1:n;if 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); %求'P,Q'V2=e(i)^2+f(i)^2;if B2(i,6)~=3DP=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; endendelseDP=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;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~=3k4=k-1;for k3=3:k4for k2=k1:N1J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2); endJ(k3,k)=0;endendfor 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+1endendICT2(a)=IT2ICT1=ICT1+1;for k=1:ndy(k)=sqrt(e(k)^2+f(k)^2);endfor i=1:nDy(k)=sqrt(e(k)^2+f(k)^2);endfor i=1:nDy(ICT1,i)=dy(i);endend %用高斯消去法解“w=-J*V”disp('迭代次数');disp(ICT1);disp('没有达到精度要求的个数');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);O(k)=atan(f(k)./e(k))*180./pi;endE=e+f*j;disp('各节点的实际电压标么值E为（节点号从小到大的排列）：');disp(E);disp('各节点的电压大小V为（节点号从小到大的排列）：');disp(V);disp('各节点的电压相角O为（节点号从小到大的排列）：');disp(O);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);enddisp('各节点的功率S为（节点号从小到大排列）：‘）；disp(S);disp('各条支路的首端功率Si为（顺序同您输入B1时一样）：‘）；for i=1:n1if B1（i,6）==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*con j(1./(B1(i,3)*B1(i,5))));disp(Si(p.q));enddisp('各条支路的末端功率Sj为（顺序同您的输入B1时一样）：‘）；for i=1:n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(xonj(E(q)./B1(i,5))-conj(E(p)))*xo nj(1./(B1(i,3)*B1(i,5))));disp(Sj(q,p));enddisp('各条支路的功率损耗DS为（顺序同您输入B1时一样）：’）；for i=1:n1if B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));endfor i=1:ICT1Cs(i)=i;enddisp('以下是每次迭代后各节点的电压值（如图所示）‘）；plot(Cs,Dy),xlabel('迭代次数’），ylabel('电压'),title('电压迭代次数曲线');。

Matlab实现潮流计算程序

程序代码如下:111111.%读入数据clcclearfilename='123.txt';a=textread(filename)n=a(1,1);pinghengjd=a(1,2);phjddianya=a(1,3);jingdu=a(1,4);b=zeros(1,9);j1=0;[m1,n1]=size(a);for i1=1:m1if a(i1,1)==0j1=j1+1;b(j1)=i1;endendb;%矩阵分块a1=a(b(1)+1:b(2)-b(1)+1,1:n1);a2=a(b(2)+1:b(3)-1,1:n1);a3=a(b(3)+1:b(4)-1,1:n1);a4=a(b(4)+1:b(5)-1,1:n1);a5=a(b(5)+1:b(6)-1,1:n1);%设置初值vcz=1;dcz=0;kmax=20;k1=0;%求节点导纳矩阵a11=zeros(4,6);for i0=1:3for j0=1:6a11(i0,j0)=a1(i0,j0);a11(4,j0)=a2(1,j0);endenda11;linei=a11(1:4,2);linej=a11(1:4,3);liner=a11(1:4,4);linex=a11(1:4,5);lineb=a11(1:4,6);branchi=0;branchj=0;branchb=0;G=zeros(4,4);B=zeros(4,4);for k=1:4i2=linei(k,1);j2=linej(k,1);r=liner(k,1);x=linex(k,1);b=0;GIJ=r/(r*r+x*x);BIJ=-x/(r*r+x*x);if k>=4 & lineb(k)~=0k0=lineb(k);G(i2,j2)=-GIJ/k0;G(j2,i2)=G(i2,j2);B(i2,j2)=-BIJ/k0;B(j2,i2)=B(i2,j2);G(i2,i2)=G(i2,i2)+GIJ/k0/k0; B(i2,i2)=B(i2,i2)+BIJ/k0/k0;elseG(j2,i2)=-GIJ;G(i2,j2)=G(j2,i2);B(j2,i2)=-BIJ;B(i2,j2)=B(j2,i2);G(i2,i2)=G(i2,i2)+GIJ;b=lineb(k);B(i2,i2)=B(i2,i2)+BIJ+b;endG(j2,j2)=G(j2,j2)+GIJ;B(j2,j2)=B(j2,j2)+BIJ+b;endG;B;B=B.*i;Yf=G+BY=abs(Yf);alf=angle(Yf);%赋Jacobian矩阵参数P=zeros(n,1);Q=zeros(n,1);Pd=zeros(1,n);Qd=zeros(1,n);dP=zeros(1,n);dQ=zeros(1,n);PG=a4(:,3);PD=a4(:,5);QG=a4(:,4);QD=a4(:,6);i8=a4(:,2);for j8=1:length(i8)P(i8(j8))=PG(i8(j8))-PD(i8(j8));Q(i8(j8))=QG(i8(j8))-QD(i8(j8));enddelt=zeros(n,1);V=ones(n,1);V(3)=1.10;V(4)=1.05;ddelt=zeros(n,1);dV=zeros(n,1);A=zeros(2*n,2*n);B=zeros(2*n,1);Jacobian=Jaco(V,delt,n,Y,alf)%求取矩阵功率for j5=1:kmaxdisp(['第' int2str(j5) '次计算结果'])if k>=kmaxbreakendfor i10=1:4Pd(i10)=0;Qd(i10)=0;for j10=1:nPd(i10)=Pd(i10)+V(i10)*Y(i10,j10)*V(j10)*cos(d elt(i10)-delt(j10)-alf(i10,j10));Qd(i10)=Qd(i10)+V(i10)*Y(i10,j10)*V(j10)*sin(d elt(i10)-delt(j10)-alf(i10,j10));endendfor i4=1:3dP(i4)=P(i4)-Pd(i4);endfor j4=1:2dQ(j4)=Q(j4)-Qd(j4);endA=Jaco(V,delt,n,Y,alf)for i14=1:nB(i14*2-1)=-dP(i14);B(i14*2)=-dQ(i14);endif max(abs(B))>jingduX=A\B;for i16=1:nddelt(i16)=X(2*i16-1);dV(i16)=X(2*i16)*V(i16);endV=V+dVdelt=delt+ddeltelsebreakenddisp('----------------')end%流氓算法% for ii=1:2% V(ii)=V(ii)+dV(ii);% end% V222222.function A=Jaco(V,delt,n,Y,alf)%计算Jacobian矩阵for i7=1:nHd1(i7)=0;Jd1(i7)=0;for j7=1:nHd1(i7)=Hd1(i7)+V(i7)*Y(i7,j7)*V(j7)*sin(delt(i7)-delt(j7)-alf(i7,j7));Jd1(i7)=Jd1(i7)+V(i7)*Y(i7,j7)*V(j7)*cos(delt(i7)-delt(j7)-alf(i7,j7));endendfor i6=1:nfor j6=1:nif i6~=j6H(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*sin(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-N(i6,j6);L(i6,j6)=H(i6,j6);elseH(i6,i6)=Hd1(i6)-V(i6)*Y(i6,i6)*V(i6)*sin(delt(i6)-delt(j6)-alf(i6,j6));J(i6,j6)=-Jd1(i6)+V(i6)*Y(i6,j6)*V(j6)*cos(delt(i6)-delt(j6)-alf(i6,j6));N(i6,j6)=-Jd1(i6)-V(i6)*Y(i6,i6)*V(i6)*cos(alf(i6,i6));L(i6,i6)=-Hd1(i6)+V(i6)*Y(i6,i6)*V(i6)*sin(alf(i6,i6));endendend%修正Jacobian矩阵for j9=3for i9=1:nN(i9,j9)=0;L(i9,j9)=0;J(j9,i9)=0;L(j9,i9)=0;endendL(j9,j9)=1;for j9=4for i9=1:nH(i9,j9)=0;N(i9,j9)=0;J(i9,j9)=0;L(i9,j9)=0;H(j9,i9)=0;N(j9,i9)=0;J(j9,i9)=0;L(j9,i9)=0;endendH(j9,j9)=1;L(j9,j9)=1;%Jaco=[H N;J L];%Jaco=zeros(2*n,2*n);for i11=1:nfor j11=1:nJaco(2*i11-1,2*j11-1)=H(i11,j11); Jaco(2*i11-1,2*j11)=N(i11,j11); Jaco(2*i11,2*j11-1)=J(i11,j11);Jaco(2*i11,2*j11)=L(i11,j11);endendA=Jaco;33333.数据:4 4 1.05 0.000011 12 0.1 0.40 0.015282 1 4 0.12 0.50 0.019203 24 0.08 0.40 0.014131 1 3 0 0.3 0.909090911 1 0 0 0.30 0.182 2 0 0 0.55 0.133 3 0.5 0 0 01 3 1.10 0 0。

(完整版)潮流计算的MATLAB源程序

％简单潮流计算的牛顿拉夫逊程序,相关的原始数据数据数据输入格式如下：％B1是支路参数矩阵，第一列和第二列是节点编号。

节点编号由小到大编写%对于含有变压器的支路，第一列为低压侧节点编号，第二列为高压侧节点%编号,将变压器的串联阻抗置于低压侧处理。

％第三列为支路的串列阻抗参数。

%第四列为支路的对地导纳参数。

%第五列为含变压器支路的变压器的变比％第六列为变压器是否含有变压器的参数,其中“1"为含有变压器，%“0”为不含有变压器。

％B2为节点参数矩阵,其中第一列为节点注入发电功率参数；第二列为节点负荷功率参数;第三列为节点电压参数；第六列为节点类型参数,其中“1”为平衡节点，“2”为PQ节点，“3”为PV节点参数。

％X为节点号和对地参数矩阵.其中第一列为节点编号，第二列为节点对地％参数。

n=input(’请输入节点数:n=');n1=input（’请输入支路数:n1=');isb=input（’请输入平衡节点号:isb='）；pr=input('请输入误差精度：pr=');B1=input（’请输入支路参数:B1=’);B2=input(’请输入节点参数：B2=’)；X=input('节点号和对地参数:X=');Y=zeros（n);Times=1; %置迭代次数为初始值％创建节点导纳矩阵for i=1:n1if B1(i,6)==0 %不含变压器的支路p=B1（i,1)；q=B1（i,2）;Y（p，q）=Y（p，q)-1/B1(i,3)；Y(q,p)=Y(p，q);Y(p，p)=Y（p，p)+1/B1(i，3)+0。

5＊B1（i,4）;Y(q，q）=Y(q，q)+1/B1(i，3)+0。

5＊B1(i，4)；else %含有变压器的支路p=B1(i,1)；q=B1(i，2);Y(p，q）=Y（p,q）—1/（B1(i,3）＊B1(i，5));Y（q，p）=Y（p,q)；Y（p，p)=Y（p,p）+1/B1（i，3)；Y（q,q)=Y(q，q）+1/(B1(i，5）^2*B1（i，3））;endendYOrgS=zeros(2＊n—2，1);DetaS=zeros(2＊n—2,1); %将OrgS、DetaS初始化％创建OrgS，用于存储初始功率参数h=0；j=0;for i=1：n ％对PQ节点的处理if i～=isb＆B2（i,6）==2h=h+1;for j=1：nOrgS(2*h-1，1）=OrgS(2*h—1，1）+real（B2（i，3））*（real（Y （i，j）)＊real(B2(j，3）)-imag（Y(i，j））*imag(B2(j,3）)）+imag(B2(i,3))＊(real（Y(i,j））*imag(B2(j，3）)+imag（Y(i，j)）*real(B2(j,3））)；OrgS（2＊h,1）=OrgS(2＊h,1）+imag（B2（i，3））＊（real(Y（i，j）)＊real(B2（j，3））—imag(Y(i，j））＊imag(B2（j,3)）)-real（B2（i,3））＊（real(Y(i，j)）＊imag(B2(j,3）)+imag（Y(i，j))＊real（B2（j，3））);endendendfor i=1:n ％对PV节点的处理，注意这时不可再将h初始化为0 if i～=isb&B2（i，6)==3h=h+1;for j=1:nOrgS（2*h—1,1）=OrgS(2*h—1,1)+real(B2（i，3）)＊(real(Y（i，j）)＊real（B2（j，3）)-imag(Y（i,j))＊imag(B2(j，3）））+imag(B2(i,3）)*(real(Y（i，j)）*imag（B2(j，3））+imag（Y(i，j))＊real(B2（j，3））);OrgS（2＊h,1)=OrgS(2*h，1）+imag（B2（i，3)）*(real（Y(i，j））*real(B2（j,3）)-imag（Y(i，j)）＊imag(B2（j,3)))—real（B2(i,3））*(real（Y（i,j)）*imag(B2(j,3)）+imag（Y（i，j))＊real(B2(j，3)))；endendendOrgS％创建PVU 用于存储PV节点的初始电压PVU=zeros(n-h—1，1）；t=0;for i=1:nif B2（i,6)==3t=t+1；PVU（t,1)=B2(i,3）;endendPVU%创建DetaS，用于存储有功功率、无功功率和电压幅值的不平衡量h=0；for i=1：n ％对PQ节点的处理if i～=isb&B2(i,6)==2h=h+1;DetaS(2*h—1，1)=real（B2（i,2））-OrgS(2*h—1，1）;DetaS(2＊h,1）=imag(B2(i，2)）—OrgS（2*h，1)；endendt=0；for i=1:n ％对PV节点的处理，注意这时不可再将h初始化为0 if i～=isb&B2（i,6）==3h=h+1;t=t+1;DetaS(2*h-1，1）=real(B2(i,2））-OrgS(2*h-1，1);DetaS（2＊h,1)=real（PVU（t，1））^2+imag（PVU（t,1）)^2-real(B2(i,3）)^2—imag(B2（i,3）)^2；endendDetaS%创建I,用于存储节点电流参数i=zeros(n—1，1）；h=0；for i=1：nif i～=isbh=h+1;I（h，1)=(OrgS(2＊h—1，1)-OrgS(2＊h，1)＊sqrt（-1))/conj（B2(i,3））;endendI％创建Jacbi（雅可比矩阵)Jacbi=zeros(2*n—2);h=0;k=0；for i=1:n ％对PQ节点的处理if B2（i,6）==2h=h+1;for j=1:nif j～=isbk=k+1;if i==j ％对角元素的处理Jacbi(2＊h-1,2＊k—1)=—imag(Y（i，j））*real(B2(i,3))+real（Y（i,j）)*imag（B2(i，3））+imag(I（h，1));Jacbi(2＊h-1,2＊k）=real（Y(i,j）)＊real(B2（i，3)）+imag(Y(i,j)）＊imag（B2(i，3））+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）)；else %非对角元素的处理Jacbi（2＊h—1,2*k-1）=—imag(Y(i,j））＊real(B2(i,3）)+real(Y（i,j）)＊imag（B2(i，3)）;Jacbi(2*h—1,2＊k)=real(Y（i，j))*real(B2（i,3))+imag(Y(i，j)）＊imag(B2（i，3））；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：n %对PV节点的处理if B2(i，6）==3h=h+1;for j=1：nif j～=isbk=k+1;if i==j ％对角元素的处理Jacbi（2＊h—1,2＊k-1）=—imag(Y(i，j））*real(B2(i，3）)+real（Y（i,j））*imag(B2（i，3)）+imag(I(h，1））；Jacbi(2＊h—1,2*k）=real（Y(i，j）)＊real(B2（i,3)）+imag(Y（i,j））*imag(B2(i，3)）+real(I(h,1））;Jacbi（2*h,2*k-1)=2＊imag(B2（i,3）)；Jacbi（2＊h，2*k)=2*real（B2（i,3）)；else %非对角元素的处理Jacbi（2*h-1，2*k—1)=-imag（Y(i,j））＊real(B2（i,3））+real(Y(i，j))＊imag(B2（i，3)）;Jacbi（2*h—1,2*k)=real(Y（i,j))＊real(B2（i，3))+imag(Y(i,j））*imag(B2（i,3））;Jacbi（2＊h，2*k—1)=0;Jacbi（2＊h,2＊k）=0;endif k==(n-1) %将用于内循环的指针置于初始值，以确保雅可比矩阵换行k=0；endendendendendJacbi％求解修正方程，获取节点电压的不平衡量DetaU=zeros(2＊n—2,1)；DetaU=inv（Jacbi)＊DetaS；DetaU％修正节点电压j=0；for i=1:n ％对PQ节点处理if B2（i,6)==2j=j+1;B2(i，3）=B2(i，3）+DetaU(2＊j,1)+DetaU(2*j—1,1）*sqrt（—1);endendfor i=1:n %对PV节点的处理if B2(i,6)==3j=j+1；B2(i，3)=B2(i,3)+DetaU（2＊j,1）+DetaU（2＊j-1，1)*sqrt(-1）;endendB2％开始循环＊**＊*******＊*******＊*＊*＊**＊＊**＊**＊＊＊＊****＊＊＊＊**＊＊***＊******＊＊＊*****＊*＊while abs（max(DetaU))>prOrgS=zeros（2＊n-2,1); ％初始功率参数在迭代过程中是不累加的，所以在这里必须将其初始化为零矩阵h=0;j=0;for i=1:nif i~=isb＆B2（i,6）==2h=h+1;for j=1：nOrgS(2*h—1,1)=OrgS(2＊h-1，1)+real（B2(i，3））＊（real(Y(i，j)）*real(B2(j，3）)—imag（Y（i,j))*imag（B2(j,3)）)+imag(B2(i,3))*（real（Y（i，j)）*imag（B2(j，3）)+imag(Y（i,j))＊real（B2(j,3）））；OrgS（2＊h,1）=OrgS（2＊h,1）+imag（B2(i,3))*（real（Y（i，j））＊real（B2(j，3））—imag(Y(i，j)）*imag(B2(j,3））)-real(B2（i，3）)*（real(Y（i,j））*imag(B2（j,3)）+imag(Y（i,j）)＊real（B2(j，3）)）；endendendfor i=1:nif i~=isb&B2(i,6）==3h=h+1；for j=1：nOrgS(2＊h-1,1)=OrgS(2*h-1,1)+real（B2（i,3））＊（real(Y(i，j））*real(B2(j，3))-imag（Y（i,j））＊imag（B2(j,3)))+imag(B2(i，3））*（real（Y（i，j））*imag(B2(j，3))+imag(Y（i，j)）*real(B2(j，3）））;OrgS(2＊h,1）=OrgS（2*h,1)+imag(B2（i，3）)＊（real（Y（i,j))*real （B2（j，3））-imag(Y（i，j）)＊imag(B2(j，3）））—real（B2（i,3）)＊(real(Y(i,j）)＊imag(B2(j，3)）+imag（Y(i，j）)*real（B2（j，3)））；endendendOrgS%创建DetaSh=0；for i=1:nif i~=isb&B2（i，6)==2h=h+1；DetaS(2＊h-1,1）=real（B2（i，2））—OrgS(2＊h—1，1）；DetaS（2＊h，1)=imag(B2（i,2））—OrgS(2＊h,1）;endendt=0;for i=1:nif i～=isb&B2（i,6)==3h=h+1;t=t+1；DetaS(2*h—1,1）=real(B2（i，2）)—OrgS（2＊h—1,1)；DetaS（2*h,1)=real（PVU（t，1）)^2+imag（PVU(t，1)）^2-real（B2（i,3）)^2-imag(B2(i，3）)^2；endendDetaS%创建Ii=zeros(n-1,1)；h=0;for i=1:nif i~=isbh=h+1;I(h,1）=（OrgS（2*h—1，1）-OrgS（2＊h,1)*sqrt（—1)）/conj（B2(i，3)）；endendI%创建JacbiJacbi=zeros(2*n—2）;h=0;k=0;for i=1:nif B2(i，6)==2h=h+1；for j=1：nif j～=isbk=k+1；if i==jJacbi（2＊h—1，2＊k—1）=-imag（Y(i，j）)＊real（B2(i，3）)+real(Y(i，j))＊imag(B2（i,3）)+imag(I(h，1)）;Jacbi(2*h-1,2*k)=real(Y（i,j））＊real（B2（i,3））+imag（Y(i,j））*imag(B2（i，3))+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)=-imag(Y（i，j))＊real(B2（i，3）)+real(Y(i,j）)＊imag（B2（i，3)）；Jacbi（2＊h—1，2＊k)=real（Y(i，j)）＊real （B2(i,3））+imag（Y（i，j）)＊imag(B2(i,3))；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，6）==3h=h+1；for j=1:nif j～=isbk=k+1;if i==jJacbi(2*h—1,2*k-1）=—imag(Y（i，j))*real(B2(i，3）)+real（Y(i,j))＊imag（B2（i，3））+imag(I(h,1）);Jacbi（2*h-1，2*k）=real（Y（i，j)）*real(B2（i,3)）+imag(Y(i，j))＊imag(B2（i,3))+real（I（h,1))；Jacbi(2*h,2＊k-1)=2＊imag(B2（i，3）);Jacbi（2＊h，2*k)=2＊real(B2(i,3））；elseJacbi(2*h-1,2＊k—1)=-imag（Y（i，j））＊real （B2(i，3））+real（Y（i，j）)＊imag（B2(i，3）)；Jacbi（2＊h-1,2*k)=real(Y(i,j））＊real（B2(i,3)）+imag（Y(i,j）)*imag（B2（i，3））;Jacbi(2＊h，2*k-1）=0；Jacbi（2＊h，2＊k）=0;endif k==（n—1）k=0；endendendendendJacbiDetaU=zeros（2*n-2，1）；DetaU=inv(Jacbi）＊DetaS；DetaU％修正节点电压j=0;for i=1:nif B2(i,6）==2j=j+1;B2（i,3)=B2(i，3)+DetaU(2*j，1)+DetaU(2＊j—1,1)＊sqrt(-1）;endendfor i=1：nif B2（i，6)==3j=j+1；B2(i,3)=B2(i,3)+DetaU（2*j,1）+DetaU（2＊j-1,1）*sqrt(—1)；endendB2Times=Times+1; %迭代次数加1endTimes一个原始数据的例子节点数 5支路数 5平衡节点编号 5精度pr 0。

潮流计算的计算机算法

高等电力系统分析（潮流计算的计算机算法）PQ分解法潮流计算（IEEE14）目录一、MATLAB源程序二、对支路参数（B1）、节点参数（B2）的说明三、带入数据，运行结果、MATLAB源程序clearclose alln=in put(' 请输入节点数:n=');n仁in put(' 请输入支路数：n 1=') Jisb=in put(' 请输入平衡节点号:isb=')pr=i nput(' 请输入误差精度:pr=');B1=i nput(' 请输入支路参数：B1 =');B2=i nput(' 请输入节点参数:B2=');n2=in put(' 请输入PQ节点个数：:n 2=');Y=zeros (n);for i=1: n1P=B1(i,1)；q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4)*1j); %Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; % Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; end 非对角元对角元disp(' 导纳矩阵Y=');disp(Y)% -------------------------------% -------------- 下面是求P,Q,V,0矩阵------------------------- V=zeros(1, n);O=zeros(1, n);P=zeros(1, n);Q=zeros(1, n);G=real(Y);B=imag(Y);for i=1: nP(i)=B2(i,3);Q(i)=B2(i,4);V(i)=B2(i,5);O(i)=B2(i,6);endB3=B(1: n-1,1: n-1); %B4=B(1: n2,1: n2); % % --------------------------- 不含平衡节点，由节点导纳虚部构成所有PQ节点% -------------- 下面是求△ P, △ Q矩阵--------------------- DX=0;ICT=1;Mp=1;Mq=1;while ICT~=0m1=1;m2=1;for i=1: nif i~=isbC(i)=O;D(i)=O;for j1=1:nC(i)=C(i)+V(i)*V(j1)*(G(i,j1)*cos(O(i)-O(j1))+B(i,j1)*si n( O(i)-O(j1)))；D(i)=D(i)+V(i)*V(j1)*(G(i,j1)*si n(O(i)-O(j1))-B(i,j1)*cos(O(i)-O( j1)))；endDP(m1)=P(i)-C(i); m1=m1+1;if B2(i,2)==1 DQ(m2)=Q(i)-D(i); m2=m2+1;endendendm1=m1-1; %m2=m2-1; %PQ DPQ=[DP';DQ']; %V仁V(:,1:m1);V2=diag(V1);V3=i nv(V2); %H=V3*DP'; %K=-i nv(B3)*H; %-deltO=V3*K; % max1=max(abs(DP)); for i=1:m1if max1<prMp=0;elseO(i)=O(i)+deltO(i)'; Mq=1;endendV4=V(:,1:m2);V5=diag(V4);V6=i nv(V5);L=V6*DQ';N=-i nv(B4)*L; deltV=N; %max2=max(abs(DQ));for i=1:m2if max2<prMq=0;elseif B2(i,2)==1;V(i)=V(i)+deltV(i)：Mp=1;end所有节点数节点数求DP,DQ对V矩阵求逆△ P/V△ P/V/B3△角=-△ P/V/V/B3 △ V=- △ Q/V/Bendendif Mp==0&&Mq==0ICT=0;elseICT=1;endDX=DX+1;end% ---------------------------------------% --------------- 迭代结束，开始输出结果disp(' ---------------------------------------- ');disp(' 迭代次数为：');disp(DX);for i=1: nE(i)=V(i)*cos(O(i))+1j*V(i)*si n( O(i)); o(i)= 180*a ngle(E(i))/pi;enddisp(' ---------------------------------------- ');disp(' 修正后各节点电压标么值为(节点号从小到大排列) :'); disp(V);disp(' ---------------------------------------- ');disp(' 修正后各节点电压相角为(节点号从小到大排列) :'); disp(o);% ----------- 计算各个节点的功率disp(' ---------------------------------------- ');disp(' 各节点的功率为：')；for p=1: nC(p)=0;for q=1: nC(p)=C(p)+conj(Y(p,q)*conj(E(q)));endS(p)=E(p)*C(p);enddisp(S);% ----------- 计算各支路的功率-------------------------for i=1: n1 p=B1(i,1);q=B1(i,2);Si(p,q)=E(p)*(conj(E(p))*conj(Y(p,p)-Y(p,q))+(conj(E(p))-conj(E(q )))*conj(Y(p,q)));disp(' ----------------------------------------- ');disp(' 各条支路的首端功率为:');disp(Si(p,q));Si(q,p)=E(q)*(conj(E(q))*conj(Y(q,q)-Y(p,q))+(conj(E(q))-conj(E(p )))*conj(Y(p,q))); disp(' ----------------------------------------- ');disp(' 各条支路的末端功率为:');disp(Si(q,p));DS(i)=Si(p,q)+Si(q,p);disp(' ----------------------------------------- ');disp(' 各条支路的功率损耗为:');disp(DS(i));end% --------------- 计算平衡节点功率----------------------Sp=0;for i=1: nSp=Sp+V( n)*conj(Y( n, i))*conj(V(i));enddisp(' ----------------------------------------- ');disp(' 平衡节点功率为：')；disp(Sp);、对支路参数（B1）、节点参数（B2）的说明1•节点数：142•支路数：203•支路矩阵B1的各支路参数：起点编号，终点编号，电阻，电抗，电导，电纳[1 2 0.01335 0.04211 0 0 ；1 3 0 0.20912 0 0J1 4 0 0.55618 0 0J1 10 0.05811 0.17632 0 0.0341 11 0.06701 0.17103 0 0.01282 10 0.05695 0.17388 0 0.03462 12 0 0.25202 0 0J2 14 0.05403 0.22304 0 0.04923 4 0 0.11001 0 0J3 13 0 0.17615 0 0 J4 5 0.03181 0.0845 0 0 ；4 9 0.12711 0.27038 0 0 ；5 6 0.08205 0.19207 0 0 ；6 12 0.09498 0.1989 0 0 ；7 8 0.22092 0.19988 0 0 ;7 12 0.12291 0.25581 0 0 ;8 9 0.17093 0.34802 0 0 ;8 12 0.06615 0.13027 0 0 ;10 11 0.04699 0.19797 0 0.043810 14 0.01938 0.05917 0 0.05284.节点参数矩阵B2的各节点参数：（对应的每一列为）节点编号，类型，注入有功，注入无功，电压幅值，电压相位其中节点类型：仁PQ 节点，2=PV 节点，0二平衡节点三、带入数据，运行结果>> clearclose all n=input （'请输入节点数：n='）;[1 1 -0.478 21 -0.076 310 0 41 -0.295 51 -0.09 61 -0.035 71 -0.061 81 -0.135 9 1 -0.149 102 0.18311 2 -0.942 12 2 -0.11213 214 00.0391 0; -0.016 1 0; 1 0; -0.166 1 0; -0.058 1 0; -0.018 1 0; -0.016 1 0; -0.058 1 0; -0.05 1 0;0 1.045 0; 0 1.01 0;0 0.1741.9 0; 0 0 1.06n1=input（'请输入支路数：n仁'）;isb=i nput（'请输入平衡节点号：isb='）;pr=input（'请输入误差精度:pr='）;B1= input（'请输入支路参数:B仁'）；B2=input（'请输入节点参数:B2='）;n2=input（'请输入PQ节点个数:n2='）;Y=zeros （n）;for i=1: n1p=B1(i,1);q=B1(i,2);Y(p,q)=Y(p,q)-1/(B1(i,3)+B1(i,4)*1j); %非对角元Y(q,p)=Y(p,q);Y(p,p)=Y(p,p)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j; % 对角元Y(q,q)=Y(q,q)+1/(B1(i,3)+B1(i,4)*1j)+B1(i,6)*1j;enddisp('导纳矩阵Y=');disp(Y)% --------------------------------% ------------- 下面是求PQ,V,0矩阵---------------V=zeros(1, n);O=zeros(1, n);P=zeros(1, n);Q=zeros(1, n); G=real(Y);B=imag(Y); for i=1: nP(i)=B2(i,3);Q(i)=B2(i,4);V(i)=B2(i,5);O(i)=B2(i,6);endB3=B（1: n-1,1: n-1）; %不含平衡节点，由节点导纳虚部构成B4=B（1:n2,1:n2）; %所有PQ 节点% -----------------------------------------% -------------- 下面是求△ P,A Q 矩阵 -------------DX=0;ICT=1;Mp=1;Mq=1;while ICT~=0m1=1;m2=1;for i=1: nif i~=isbC(i)=0;D(i)=0;for j1=1: nC(i)=C(i)+V(i)*V(j1)*(G(i,j1)*cos(O(i)-O(j1))+B(i,j1)*si n(O(i)-O(j1))); D(i)=D(i)+V(i)*V(j1)*(G(i,j1)*si n(O(i)-O(j1))-B(i,j1)*cos(O(i)-O(j1))); endDP(m1)=P(i)-C(i); m1=m1+1;if B2(i,2)==1DQ(m2)=Q(i)-D(i); m2=m2+1;max1=max(abs(DP));for i=1:m1if max1<prMp=0;elseO(i)=O(i)+deltO(i)';Mq=1;endendV4=V(:,1:m2);V5=diag(V4); V6=i nv(V5);L=V6*DQ';N=-i nv(B4)*L;deltV=N;%△ V=-A Q/V/B max2=max(abs(DQ));for i=1:m2if max2<pr endendend m1=m1-1; m2=m2-1;DPQ=[DP';DQ'];V1=V(:,1:m1);V2=diag(V1);V3=i nv(V2);H=V3*DP';K=-i nv(B3)*H;deltO=V3*K; % 所有节点数 %PQ 节点数 %求 DPQQ%对V 矩阵求逆 %A P/V %-A P/V/B3%A 角=-AMq=0;elseif B2(i,2)==1;V(i)=V(i)+deltV(i)： Mp=1;endendendif Mp==0&&M q==0ICT=0;elseICT=1;endDX=DX+1;end% ---------------------------------------% --------------- 迭代结束，开始输出结果disp(' ---------------------------------------- ');disp('迭代次数为：');disp(DX);for i=1: nE(i)=V(i)*cos(O(i))+1j*V(i)*si n(O (i));o(i)= 180*a ngle(E(i))/pi;enddisp(' ---------------------------------------- ');disp('修正后各节点电压标么值为(节点号从小到大排列) :'); disp(V);disp(' ---------------------------------------- ');disp('修正后各节点电压相角为(节点号从小到大排列) :'); disp(o);% ----------- 计算各个节点的功率------------------disp(' ---------------------------------------- ');disp('各节点的功率为:');for p=1: nC(p)=0;for q=1: nC(p)=C(p)+conj(Y(p,q)*conj(E(q)));endS(p)=E(p)*C(p); end disp(S);% ----------- 计算各支路的功率------------------for i=1: n1 p=B1(i,1);q=B1(i,2);Si(p,q)=E(p)*(conj(E(p))*conj(Y(p,p)-Y(p,q))+(conj(E(p))-conj(E(q)))*conj(Y(p,q))); disp(' ----------------------------------------- ');disp('各条支路的首端功率为:');disp(Si(p,q));Si(q,p)=E(q)*(conj(E(q))*conj(Y(q,q)-Y(p,q))+(conj(E(q))-conj(E(p)))*conj(Y(p,q)));disp(' ----------------------------------------- ');disp('各条支路的末端功率为:');disp(Si(q,p));DS(i)=Si(p,q)+Si(q,p);disp(' ----------------------------------------- ');disp('各条支路的功率损耗为:');disp(DS(i));end% --------------- 计算平衡节点功率------------Sp=0;for i=1: nSp=Sp+V( n)*conj(Y( n,i))*conj(V(i));enddisp(' ----------------------------------------- ');disp('平衡节点功率为:');disp(Sp);请输入节点数：n=14请输入支路数：n 1=20请输入平衡节点号:isb=14请输入误差精度：pr=0.00001请输入支路参数:B1=[1 2 0.01335 0.04211 0 0 ;1 3 0 0.20912 0 0 51 4 0 0.55618 0 051 10 0.05811 0.176320.0341 11 0.06701 0.17103 0 0.01282 10 0.05695 0.17388 0 0.03462 12 0 0.25202 0 052 14 0.05403 0.22304 0 0.049210.5130 -38.2963i -6.8410 +21.5786i0.0000 + 4.7819i 0.0000 + 1.7980i0.0000 +0.0000i-6.8410 +21.5786i9.5680 -34.8916i0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 +0.0000i0.0000 + 4.7819i0.0000 + 0.0000i0.0000 -19.5490i 0.0000 + 9.0901i0.0000 +0.0000i0.0000 + 1.7980i0.0000 + 0.0000i0.0000 + 9.0901i5.3261 -24.2825i-3.9020+10.3654i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-3.9020 +10.3654i5.7829-14.7683i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i 0.0000 + 0.0000i-1.8809 +4.4029i0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 +O.OOOOi3 40 0.11001 00 ; 3 13 0 0.17615 0 0 54 5 0.03181 0.0845 0 0 ; 4 9 0.12711 0.27038 0 0 ; 5 6 0.08205 0.19207 0 0 ; 6 12 0.094980.1989 0 0 ; 7 8 0.22092 0.19988 0 0 ; 7 12 0.12291 0.25581 0 0 ; 8 9 0.17093 0.34802 0 0 ; 8 12 0.06615 0.130270 0 ; 10 11 0.046990.19797 00.0438;10 14 0.01938 0.05917 0 0.0528 ;请输入节点参数：B2=[1 1 -0.478 0.039 2 1 -0.076 -0.016 10;3 1 0 0 1 0;4 1 -0.295 -0.166 1 0; 5 1 -0.09 -0.058 1 0; 6 1 -0.035 -0.018 1 0;7 1 -0.061 -0.016 1 0;8 1 -0.135 -0.058 1 0;9 1 -0.149 -0.05 10; 10 2 0.183 01.0450;11 2 -0.942 0 1.01 0;12 2 -0.112 0.0471.7 0;13 2 0 0.174 1.90;14 00 01.06 0;]请输入PQ 节点个数:n2=9导纳矩阵 Y=Columns 1 through 5；0.0000 + O.OOOOi O.OOOOi0.0000 + 0.0000i 0.0000i-1.6860 + 5.1158i 0.0000i-1.9860 + 5.0688i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 5.115 8i0.0000 + 0.0000i 5.1939i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i-1.8809 + 4.4029i 0.0000i3.8359 - 8.4970i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i -30.1895i0.0000 + 0.0000i 4.7819i-1.9550 + 4.0941i 0.0000i0.0000 + 0.0000i 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.7011 + 5.1939i0.0000 + 0.0000i0.0000 + 3.9679i0.0000 + 0.0000i-1.0259 + 4.2350i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i4.0150 -5.4279i-2.4890 + 2.2520i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.5260 + 3.1760i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-2.4890 + 2.2520i6.7249 -10.6697i-1.1370 + 2.3150i0.0000 + 0.0000i0.0000 + 0.0000i-3.0989 + 6.1028i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.4240 + 3.0291i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.4240 + 3.0291i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.1370 + 2.3150i2.5610 - 5.3440i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000iColumns 6 through 100.0000 +0.0000 +0.0000 +0.0000 + 0.0000 +0.0000 +0.0000 + -1.6860 + -1.7011 +0.0000 +0.0000 +0.0000 +0.0000 +0.0000 + 0.0000 +0.0000 +9.5213 -1.1350 +0.0000 +0.0000 ++15.2631iColu mns 11 through 14迭代次数为:42修正后各节点电压标么值为 Colu mns 1 through 10（节点号从小到大排列）1.21281.21481.58091.04501.5627 1.5364 1.5602 1.6264 1.6792 1.6654Colu mns 11 through 141.0100 1.7000 1.90001.0600修正后各节点电压相角为（节点号从小到大排列）Colu mns 1 through 10-13.1377 -11.9998-14.4540 -15.0006 -15.3011 -15.7590 -16.4954 -16.3362-16.0220 -5.2081Colu mns 11 through 14-12.1568 -16.1918 -14.4540-1.9860 + 5.0688i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -1.1350 + 4.7819i 3.1210 - 9.7941i 0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 3.9679i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -1.9550 + 4.0941i -1.5260 + 3.1760i -3.0989 + 6.1028i 0.0000 + 0.0000i0.0000 + 0.0000i 0.0000 + 0.0000i 6.5799 -17.3407i 0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 5.6770i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 - 5.6770i 0.0000 + 0.0000i0.0000 + 0.0000i-1.0259 + 4.2350i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i-4.9991 +15.2631i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 6.0250 -19.3961i各节点的功率为：Columns 1 through 51.4718 - 1.8800i 0.9041 - 1.1423i 0.0000 - 0.0000i -0.1314 + 0.1399i 0.0646iColumns 6 through 10-0.0300 + 0.0035i 0.0039 - 0.1228i -0.0656 - 0.0724i -0.0352 - 0.1437i1.3905iColu mns 11 through 14-0.4870 - 1.2396i 3.2005 + 3.0414i 1.7586 + 3.1847i -2.7057 + 0.3992i各条支路的首端功率为26.1711 +87.9128i各条支路的末端功率为23.5667 +83.4812i各条支路的功率损耗为：4.9738e+01 + 1.7139e+02i各条支路的首端功率为15.2545 +65.3872i各条支路的末端功率为0.2082 +56.8008i各条支路的功率损耗为：1.5463e+01 + 1.2219e+02i -0.0475 - -1.5950 -各条支路的首端功率为15.3538 +59.6752i各条支路的末端功率为12.6813 +60.6685i各条支路的功率损耗为：2.8035e+01 + 1.2034e+02i各条支路的首端功率为18.4735 +62.4536i各条支路的末端功率为11.6194 +39.6840i各条支路的功率损耗为：3.0093e+01 + 1.0214e+02i各条支路的首端功率为18.0012 +62.4930i各条支路的末端功率为5.5097 +16.2405i各条支路的功率损耗为23.5109 +78.7335i各条支路的首端功率为17.0452 +57.7870i各条支路的末端功率为11.7622 +39.7706i各条支路的功率损耗为28.8075 +97.5576i各条支路的首端功率为13.5220 +59.6677i各条支路的末端功率为19.6150 +58.2875i各条支路的功率损耗为：3.3137e+01 + 1.1796e+02i各条支路的首端功率为16.5471 +56.5546i各条支路的末端功率为6.9281 +27.4025i各条支路的功率损耗为23.4752 +83.9571i各条支路的首端功率为-0.2082 +69.5644i各条支路的末端功率为12.7806 +79.1441i各条支路的功率损耗为：1.2572e+01 + 1.4871e+02i各条支路的首端功率为0.0000 +64.5962i各条支路的末端功率为-0.0000 +37.3498i各条支路的功率损耗为：-1.7764e-15 + 1.0195e+02i各条支路的首端功率为21.7957 +82.2160i各条支路的末端功率为23.5612 +60.7482i各条支路的功率损耗为：4.5357e+01 + 1.4296e+02i各条支路的首端功率为15.8994 +64.7377i各条支路的末端功率为9.9896 +20.6498i各条支路的功率损耗为25.8891 +85.3875i各条支路的首端功率为18.7607 +47.1609i各条支路的末端功率为15.0088 +33.6107i各条支路的功率损耗为33.7695 +80.7716i各条支路的首端功率为15.4666 +33.8366i各条支路的末端功率为24.5068 +61.3933i各条支路的功率损耗为39.9734 +95.2299i各条支路的首端功率为18.2997 +21.5842i各条支路的末端功率为25.5956 +35.9107i各条支路的功率损耗为43.8954 +57.4949i各条支路的首端功率为15.7256 +24.3488i各条支路的末端功率为23.3240 +59.2040i各条支路的功率损耗为39.0496 +83.5528i各条支路的首端功率为21.6792 +35.6719i各条支路的末端功率为9.3603 +19.4667i各条支路的功率损耗为31.0395 +55.1385i各条支路的首端功率为27.4696 +46.8495i各条支路的末端功率为27.7461 +67.4149i各条支路的功率损耗为：5.5216e+01 + 1.1426e+02i各条支路的首端功率为10.9761 +38.1226i各条支路的末端功率为4.9835 +14.8560i各条支路的功率损耗为15.9596 +52.9785i各条支路的首端功率为17.4469 +49.3022i各条支路的末端功率为10.7497 +39.1332i各条支路的功率损耗为28.1966 +88.4354i平衡节点功率为：-0.0889 - 0.5670i。

matlab loadflow代码

一、介绍loadflow代码在Matlab中的应用在电力系统中，loadflow（潮流计算）是一种非常重要的分析方法，它用来计算电力系统中各个节点的电压、相角和功率等参数。

在Matlab中，我们可以编写loadflow代码来实现对电力系统进行潮流计算，以便进行系统的稳定性分析、负荷分配等工作。

二、loadflow代码的基本原理loadflow代码的基本原理是通过迭代计算来求解电力系统中各节点的电压和相角。

一般来说，loadflow代码的实现需要根据节点的功率平衡方程、节点电压幅值方程、支路功率方程等，采用牛顿-拉夫逊法或高斯-赛德尔法等迭代方法来求解。

三、编写loadflow代码的步骤1. 收集电力系统的基本参数：首先需要收集电力系统的基本参数，包括节点的功率负荷、支路的阻抗和导纳等。

这些参数将作为loadflow 代码的输入。

2. 建立loadflow模型：根据电力系统的实际情况，构建电力系统的潮流计算模型，包括节点的功率平衡方程、节点电压幅值方程、支路功率方程等。

3. 编写loadflow代码：根据建立的loadflow模型，使用Matlab语言编写loadflow代码，实现节点电压和相角的迭代计算。

在编写过程中，需要注意代码的结构，保证其逻辑清晰，易于阅读和维护。

4. 验证loadflow代码：完成loadflow代码的编写后，需要进行验证，确保代码的准确性和稳定性。

通过对实际电力系统的数据进行潮流计算，对比计算结果与实际情况，来验证loadflow代码的准确性。

四、loadflow代码的优化1. 代码效率优化：在编写loadflow代码时，需要考虑代码的运行效率。

可以通过优化算法、并行计算等手段来提高代码的运行速度，以适应大规模电力系统的潮流计算需求。

2. 界面优化：为方便用户使用，可以在Matlab中设计图形用户界面（GUI），以便用户输入电力系统参数、进行潮流计算，并可视化输出计算结果。

(完整word版)基于MATLAB进行潮流计算

基于MATLAB进行潮流计算学生：王仕龙2011148213指导老师：李咸善摘要：电力系统潮流计算方法有两类，即手算潮流和计算机潮流计算。

手算潮流主要借助于形成简化的等值电路来实现，这种方法尤其适用于规模不大的辐射型电力潮流计算。

计算机潮流计算的实现有两种途径：其一是编程实现网络方程的迭代求解；其二是借助与电力系统分析仿真软件，搭建系统模型来完成潮流计算。

MATLAB具有强大的矩阵运算功能，同时其具有电力系统仿真平台也为直观地实现潮流计算提供了更便捷的手段[1]。

本文是基于MATLAB软件，采用极坐标形式牛顿─拉夫逊法进行潮流计算，为其他形式的潮流计算有借鉴的作用。

关键词：电力系统；计算机潮流计算；MATLAB ；牛顿─拉夫逊法Abstract:The power flow calculation method has two kinds,which are the hand calculation of tidal current and computer power flow calculation.Hand calculation tidal current is mainly realized by means of the formation of simplified equivalent circuit.This method is especially suitable for small scale radiation power flow calculation.There are two ways to realize the computer power flow calculation.The first one is through the programming iteration for solving network equation,the second one is with the help of analysis of power system simulation software to build the system model complete the power flow calculation.The software of MATLAB has strong matrix function,.At the same time,It’s power system simulation platform provides a more convenient means to realize power flow calculation intuitively[1].This paper is based on the software of MATLAB to calculate the power flow calculation by adopting the form of Newton-Raphson method of power flow calculation of polar coordinates.And it can be the role of reference of other forms of power flow calculation.Key words: power system computer; power flow calculation;MATLAB;Newton-Raphson1.计算原理电力系统潮流是指系统中所有运行参数的总体，包括各个母线电压的大小和相位，各个发电机和负荷的功率及电流，以及各个变压器和线路等元件所通过的功率,电流和其中的损耗。

基于matlab的潮流计算

课程设计任务书第学期学院：信息商务学院专业：电气工程及其自动化学生姓名：学号课程设计题目：基于MATLAB的潮流计算起迄日期：月日～月日课程设计地点：电气工程系工程实验室指导教师：系主任：课程设计任务书课程设计任务书1 引言潮流计算是电力系统最基本最常用的计算。

根据系统给定的运行条件，网络接线及元件参数，通过潮流计算可以确定各母线的电压（幅值和相角），各支路流过的功率，整个系统的功率损耗。

潮流计算是实现电力系统安全经济发供电的必要手段和重要工作环节。

因此，潮流计算在电力系统的规划计算，生产运行，调度管理及科学计算中都有着广泛的应用。

潮流计算在数学上是多元非线性方程组的求解问题，牛顿—拉夫逊Newton-Raphson法是数学上解非线性方程组的有效方法，有较好的收敛性。

运用电子计算机计算一般要完成以下几个步骤：建立数学模型，确定解算方法，制订计算流程，编制计算程序。

1.1 设计任务与要求利用MATLAB/SIMULINK仿真软件，搭建多电压等级电力网络潮流计算模型。

掌握各种潮流计算方法，掌握MATLAB/SPS工具箱使用方法及其在潮流计算中的应用。

MATLAB作为先进的仿真软件，在电力系统潮流计算方面，尤其是矩阵运算上具有突出的优秀性。

本设计要求利用MATLAB编写程序，输入数据并完成潮流分布计算。

1.2实用价值与理论意义在电力系统的正常运行中，随着用电负荷的变化和系统运行方式的改变，网络中的损耗也将发生变化。

要严格保证所有的用户在任何时刻都有额定的电压是不可能的，因此系统运行中个节点出现电压的偏移是不可避免的。

为了保证电力系统的稳定运行，要进行潮流调节。

随着电力系统及在线应用的发展，计算机网络已经形成，为电力系统的潮流计算提供了物质基础。

电力系统潮流计算是电力系统分析计算中最基本的内容，也是电力系统运行及设计中必不可少的工具。

根据系统给定的运行条件、网络接线及元件参数，通过潮流计算可以确定各母线电压的幅值及相角、各元件中流过的功率、整个系统的功率损耗等。

(完整版)基于MATLAB牛顿拉夫逊法进行潮流计算

>> %本程序的功能是用牛顿拉夫逊法进行潮流计算n=input('请输入节点数:n=');nl=input('请输入支路数:nl=');isb=input('请输入平衡母线节点号:isb=');pr=input('请输入误差精度:pr=');B1=input('请输入由各支路参数形成的矩阵:B1='); B2=input('请输入各节点参数形成的矩阵:B2=');Y=zeros(n); e=zeros(1,n);f=zeros(1,n);V=zeros(1,n); O=zeros(1,n);S1=zeros(nl);for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endY(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));Y(q,p)=Y(p,q);Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;end%求导纳矩阵disp('导纳矩阵Y=');disp(Y);G=real(Y);B=imag(Y);for i=1:ne(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4);endfor i=1:nS(i)=B2(i,1)-B2(i,2);B(i,i)=B(i,i)+B2(i,5);endP=real(S);Q=imag(S);ICT1=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;%117页malihong打if B2(i,6)~=3DP=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;endendelseDP=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);%118页孟打印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;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;endendend %119页zhengtong打%for 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*V”disp('迭代次数');disp(ICT1);disp('没有达到精度要求的个数');disp(ICT2);for k=1:nV(k)=sqrt(e(k)^2+f(k)^2);shita(k)=atan(f(k)./e(k))*180/pi;E(k)=e(k)+f(k)*j;enddisp('各节点的实际电压标么值E为(节点号从小到大排列):');disp(E);disp('各节点的电压大小V为(节点号从小到大排列):');disp(V);disp('各节点的电压相角时shita为(节点号从小到大排列):');disp(shita);for p=1:nC(p)=0;for q=1:nC(p)=C(p)+conj(Y(p,q))*conj(E(q));endS(p)=E(p)*C(p);end %120页打disp('各节点的功率S为(节点号从小到大排列):');disp(S);disp('各条支路的首端功率Si为(顺序同您输入B1时一样):');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSi(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))-conj(E(q)))*conj(1./(B1(i,3)*B1(i, 5))));disp(Si(p,q));enddisp ('各条支路的末端功率Sj为(顺序同您输入B1时一样):');for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endSj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))-conj(E(p)))*conj(1./(B1(i,3)*B1(i, 5))));disp(Sj(q,p));enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一样):' );for i=1:nlif B1(i,6)==0p=B1(i,1);q=B1(i,2);else p=B1(i,2);q=B1(i,1);endDS(i)=Si(p,q)+Sj(q,p);disp(DS(i));end。

基于Matlab编程的NR法潮流程序

Line=[1,7,1, 0.05293, 0.050815,0;2,7,4, 0.10397, 0.09981,0;3,7,5, 0.04348, 0.04174,0;4,6,2, 0.071835, 0.06896,0;5,6,3, 0.0397, 0.03796,0;6,6,4, 0.04537,0.04356,0;];%线路数据 Line[1-线路编号 2-线路首端节点号 3-线路末端节点号 4-支路电阻 5-支路电抗 6-支路电纳(此处取B/2)] Nbus=7;%节点数Nline=6;%线路个数Slack=7;%平衡节点号Npq=5;%PQ节点数%以下计算节点导纳矩阵Y=zeros(Nbus);%节点导纳矩阵初始化for k=1:Nline %根据每条线路来处理节点导纳t1=Line(k,2);t2=Line(k,3);b=Line(k,6);%分别提取出线路的首端节点号t1、末端节点号t2和线路导纳bY1=1/(Line(k,4)+j*Line(k,5));%计算线路的支路电导Y1Y(t1,t1)=Y(t1,t1)+Y1+j*b;%Y(t1,t2)=Y(t1,t2)-Y1;%Y(t2,t1)=Y(t2,t1)-Y1;%Y(t2,t2)=Y(t2,t2)+Y1+j*b;%endG=real(Y);B=imag(Y);%%以下给节点电压赋初值u=[1,1,1,1,1,1.1,1.1;];%us=u;delt=[0,0,0,0,0,0,0;];%for m=1:Nbuse(m)=u(m)*cos(delt(m));f(m)=u(m)*sin(delt(m));v(m)=e(m)+i*f(m);endp=[-0.075,-0.065,-0.092,-0.085,-0.094,0.114];%q=[-0.035,-0.038,-0.045,-0.038,-0.06,0.07364];%e=real(u);f=imag(u);k=0;a=1;%while a>0.00001 %if k>100 ,break;elsefor m=1:Nbusu(m)=sqrt(e(m)^2+f(m)^2);delt(m)=atan(f(m)/e(m));v(m)=e(m)+i*f(m);endfor m=1:Npqfor n=1:Nbuspt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt( n)));qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt( n)));endpp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);endfor m=Npq+1:Nbus-1for n=1:Nbuspt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt( n)));endpp(m)=p(m)-sum(pt);uu(m)=us(m)^2-u(m)^2;end%以下计算m=n情况下的Jacobi矩阵中的子矩阵元素for m=1:Nbus-1for n=1:Nbush(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))-B(m,n)*sin(delt(m)-delt(n )));n0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))+B(m,n)*cos(delt(m)-delt( n)));endfor n=1:Nbusl(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))+B(m,n)*cos(delt(m)-delt(n )));j(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))-B(m,n)*sin(delt(m)-delt(n )));endH(m,m)=-sum(h)-u(m)^2*(G(m,m)*cos(delt(m)-delt(m))-B(m,m)*sin(delt(m) -delt(m)));N(m,m)=-sum(n0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m )-delt(m)));L(m,m)=sum(l)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));J(m,m)=-sum(j)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m) -delt(m)));endfor m=1:Nbus-1Jacobi(2*m-1,2*m-1)=H(m,m);%P对e求导Jacobi(2*m-1,2*m)=N(m,m);%P对f求导Jacobi(2*m,2*m-1)=L(m,m);%Q对e求导Jacobi(2*m,2*m)=J(m,m);%Q对f求导end%以下计算m！=n情况下的Jacobi矩阵中的子矩阵元素for m=1:Nbus-1for n=1:Nbus-1if m==nelseH(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-del t(n)));J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt (n)));N(m,n)=-J(m,n);L(m,n)=H(m,n);Jacobi(2*m-1,2*n-1)=N(m,n);%P对eJacobi(2*m-1,2*n)=L(m,n);%P对fJacobi(2*m,2*n-1)=H(m,n);%Q对eJacobi(2*m,2*n)=J(m,n);%Q对fendendendfor m=Npq+1:Nbus-1for n=1:Nbus-1if n==mJacobi(2*m,2*n)=-2*imag(v(m));Jacobi(2*m,2*n-1)=-2*real(v(m));elseJacobi(2*m,2*n)=0;Jacobi(2*m,2*n-1)=0;endendendfor m=1:NpqUnbanlance(2*m-1)=pp(m);Unbanlance(2*m)=qq(m);endfor m=Npq+1:Nbus-1Unbanlance(2*m-1)=pp(m);Unbanlance(2*m)=uu(m);endCorrection=-Jacobi\(Unbanlance');a=max(abs(Correction));for m=1:Nbus-1e(m)=e(m)+Correction(2*m-1);f(m)=f(m)+Correction(2*m);endk=k+1;endendk,u,delt,Jacobi,afor m=1:NbusU(m)=u(m)*(cos(delt(m))+i*sin(delt(m)));I(m)=Y(Nbus,m)*U(m);endfor m=1:Nbusfor n=1:NbusS(m,n)=U(m)*(conj(U(m))-conj(U(n)))*conj(-Y(m,n)); endendS。

matlab潮流计算

matlab潮流计算%线路数据line=[1-节点编号2-线路首端节点号3-线路末端节点号4-支路电阻5-支路电抗6-支路电纳(注意：此处取的是b/2)]line=[1120.040.250.252130.10.3503230.080.300.25];%变压器数据transform=[1-支路编号2-支路首节点编号3-支路末节点编号4-支路电阻(p.u.)5-支路电抗(p.u.)6-变压器变比(p.u.)]transform=[12400.0151.0523500.031.05];%数据预处理nbus=5;%节点数nline=size(line,1);%线路个数ntrans=size(transform,1);%变压器个数slack=5;%均衡节点号npq=4;%pq节点的个数ss=0;%计算节点导纳矩阵y=zeros(nbus);ifnline>=1%推论与否存有线路fork=1:nline%以下处理线路t1=line(k,2);t2=line(k,3);b2=line(k,6);%分别抽出线路的首端节点编号t1、末端节点编号t2和对地电纳b2yl=1/(line(k,4)+j*line(k,5));%计算线路的支路电导yly(t1,t1)=y(t1,t1)+yl+j*b2;%修正第k条线路首端节点的自导纳y(t1,t2)=y(t1,t2)-yl;%修正第k条线路首端节点与末端节点之间的互磁滞y(t2,t1)=y(t2,t1)-yl;%修正第k条线路末端节点与首端节点之间的互磁滞y(t2,t2)=y(t2,t2)+yl+j*b2;%修正第k条线路末端节点的自磁滞endendifntrans>=1%判断是否存在变压器fork=1:ntrans%以下处理变压器t1=transform(k,2);t2=transform(k,3);t3=transform(k,6);%分别抽出变压器的首端节点编号t1、末端节点编号t2和变比t3yt=1/(transform(k,4)+j*transform(k,5));yt1=yt/t3;yt2=yt*(1-t3)/(t3*t3);yt3=yt*(t3-1)/t3;y(t1,t1)=y(t1,t1)+yt1+yt2;y(t1,t2)=y(t1,t2)-yt1;y(t2,t1)=y(t2,t1)-yt1;y(t2,t2)=y(t2,t2)+yt1+yt3;endendg=real(y);b=imag(y);%区分节点磁滞矩阵的实部和虚部gb%赋初值delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;u(1)=1;u(2)=1;u(3)=1;u(4)=1;p(1)=-0.30;q(1)=-0.18;p(2)=-0.55;q(2)=-0.13;p(3)=0;q(3)=0;p(4)=0.5;q(4)=1.10;p(5)=0.8;q(5)=0.50;k=0;precision=1;npq=4;%npq分别就是网络中的pq节点数%[unbalance]=-[jacobi][correction]whileprecision>0.00001%设定误差上限，判断是否继续迭代u(5)=1.06;delt(5)=0;%设定平衡节点电压相角与幅值k;u;delt;form=1:npqforn=1:nbuspt(n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));%由节点电压求出的pq节点转化成军功功率qt(n)=u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));%由节点电压求得的pq节点注入无功功率endunbalance(2*m-1)=p(m)-sum(pt);%排序pq节点军功功率不平来衡量unbalance(2*m)=q(m)-sum(qt);%排序pq节点无功功率不平来衡量end%[unbalance]是节点不平衡量矩阵form=1:npqforn=1:nbush0(n)=u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));n0(n)=-u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));j0(n)=-u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));l0(n)=-u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));endh(m,m)=sum(h0)-u(m)^2*(g(m,m)*sin(delt(m)-delt(m))-b(m,m)*cos(delt(m)-delt(m)));n(m,m)=sum(n0)+u(m)^2*(g(m,m)*cos(delt(m)-delt(m))+b(m,m)*sin(delt(m)-delt(m)))-2*u(m)^2*g(m,m);j(m,m)=sum(j0)+u(m)^2*(g(m,m)*cos(delt(m)-delt(m))+b(m,m)*sin(delt(m)-delt(m)));l(m,m)=sum(l0)+u(m)^2*(g(m,m)*sin(delt(m)-delt(m))-b(m,m)*cos(delt(m)-delt(m)))+2*u(m)^2*b(m,m);jacobi(2*m-1,2*m-1)=h(m,m);jacobi(2*m-1,2*m)=n(m,m);jacobi(2*m,2*m-1)=j(m,m);jacobi(2*m,2*m)=l(m,m);end%计算m=n情况下的jacobi矩阵中的子矩阵元素form=1:npqforn=1:npqifm==nelseh(m,n)=-u(m)*u(n)*(g(m,n)*sin(delt(m)-delt(n))-b(m,n)*cos(delt(m)-delt(n)));j(m,n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));n(m,n)=-j(m,n);l(m,n)=h(m,n);jacobi(2*m-1,2*n-1)=h(m,n);jacobi(2*m-1,2*n)=n(m,n);jacobi(2*m,2*n-1)=j(m,n);jacobi(2*m,2*n)=l(m,n);endendend%排序m≠n情况下的jacobi矩阵中的子矩阵元素correction=-jacobi\\(unbalance');%计算电压相角和幅值的修正量precision=max(abs(correction));%取误差最大值form=1:npqdelt(m)=delt(m)+correction(2*m-1);%修正pq节点电压相角u(m)=u(m)+correction(2*m);%修正pq节点电压幅值endk=k+1;%运算轮数+1endk,u,delt,jacobi,precisionform=1:nbusu(m)=u(m)*(cos(delt(m))+j*sin(delt(m)));%采用直角坐标系表示电压i(m)=y(nbus,m)*u(m);%计算注入平衡节点的电流endsslack=u(nbus)*sum(conj(i))%排序转化成均衡节点的功率form=1:nbusforn=1:nbuss(m,n)=u(m)*(conj(u(m))-conj(u(n)))*conj(-y(m,n));%排序线路功率ploss(m,n)=u(m)*u(n)*(g(m,n)*cos(delt(m)-delt(n))+b(m,n)*sin(delt(m)-delt(n)));endends%表明线路功率矩阵，s(m,n)则表示假设功率由节点m流向节点n%若数值为+则说明实际功率流向与假设方向相同，若数值为-则说明实际功率流向与假设方向相反。

matlab110kv潮流计算仿真代码

matlab110kv潮流计算仿真代码以下是一个基于 MATLAB 的 110KV 潮流计算仿真代码，仅供参考:```matlab% 导入必要的文件net = read_net("110kv_power_grid.txt"); % 读取电网模型文件load("grid_voltage.mat"); % 加载电压等级load("grid_current.mat"); % 加载电流等级% 设置潮流计算参数dt = 0.01; % 时间步长t = 0:dt:10; % 时间序列N = 100; % 点数omega = 6.283; % 谐波频率f = omega/(2*pi); % 谐波频率N_oil = 30; % 油浸变压器个数t_oil = 20; % 油浸变压器时间步长t_oil = round(t_oil/dt); % 油浸变压器时间步长omega_oil = omega/(2*pi)*t_oil; % 油浸变压器谐波频率i_oil = zeros(1,N_oil); % 油浸变压器电流for i = 1:N_oilbeta_oil = 0.1; % 油浸变压器 beta 系数omega_oil_prime = omega_oil*beta_oil; % 油浸变压器谐波频率t_oil_prime = t_oil/omega_oil_prime; % 油浸变压器时间步长i_oil(i) =exp(-t_oil_prime*omega_oil)*sin(2*pi*f*t_oil(i) +omega_oil_prime*t(i) + theta_oil(i));end% 计算潮流net_out = zeros(size(i_oil));for i = 1:Nnet_out(i) =sum(i_oil(i)*exp(-t_oil_prime*omega_oil)*sin(2*pi*f*t_oil(i) + omega_oil_prime*t(i) + theta_oil(i)));end% 可视化结果fplot("i_oil",net_out);hold on;fplot("theta_oil",net_out);hold off;title("110KV 潮流计算");```在此代码中，我们使用了 MATLAB 的`read_net`函数读取电网模型文件，`load`函数加载电压等级和电流等级，`grid_voltage`和`grid_current`函数分别加载了电压等级和电流等级，设置了潮流计算参数，并使用`fplot`函数绘制了潮流结果。

MATLAB 极坐标求解潮流计算

电力系统分析大作业——潮流计算班级：XXX姓名：XXX学号：XXX程序清单：%直角坐标法求解潮流计算clear;close;clc;Branch=[1 2 0.1+0.4j 0.01528j 1 03 1 0.3j inf 1.1 11 4 0.12+0.50j 0.01920j 1 02 4 0.08+0.40j 0.01413j 1 0]%Branch矩阵：1、支路首端号；2、支路末端号；3、支路阻抗；4、支路对地导纳；5、支路的变化；6、支路首端处于K侧为1,1侧为0Y=zeros(4); %节点导纳矩阵for i=1:4if Branch(i,6)==0 %不含变压器的支路j=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-1/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3)+Branch(i,4);Y(k,k)=Y(k,k)+1/Branch(i,3)+Branch(i,4);elsej=Branch(i,1);k=Branch(i,2);Y(j,k)=Y(j,k)-Branch(i,5)/Branch(i,3);Y(k,j)=Y(j,k);Y(j,j)=Y(j,j)+1/Branch(i,3);Y(k,k)=Y(k,k)+Branch(i,5)^2/Branch(i,3);endenddisp('节点导纳矩阵');YG=real(Y);B=imag(Y);V=[1;1;1.1;1.05];%给定V的初始计算值e=real(V);f=imag(V);disp('节点电压的实部：')edisp('节点电压的虚部：')fdisp('节点注入有功功率：')Ps=[-0.3;-0.55;0.5;0]disp('节点注入无功功率：')Qs=[-0.18;-0.13;0;0]%由各节点电压向量（状态变量）可得各节点注入功率：P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];n=0;%辅助循环计数变量while (any(del_W>1e-5))&(n<5)n=n+1;disp('迭代次数：');disp(n)%--------------------------------------------------------%修正方程式：% [-0.30-P(1) ] [del_e(1)]% [-0.18-Q(1) ] [del_f(1)]% [-0.55-P(2) ] [del_e(2)]% [-0.13-Q(2) ] + J * [del_f(2)] = 0 % [0.5-P(3) ] [del_e(3)]% [1.1^2-e(3)^2-f(3)^2] [del_f(3)]%雅克比矩阵为：J=[P1_e1 P1_f1 P1_e2 P1_f2 P1_e3 P1_e3;% Q1_e1 Q1_f1 Q1_e2 Q1_f2 Q1_e3 Q1_f3% P2_e1 P2_f1 P2_e2 P2_f2 P2_e3 P2_f3% Q2_e1 Q2_f1 Q2_e2 Q2_f2 Q2_e3 Q2_f3% P3_e1 P3_f1 P3_e2 P3_f2 P3_e3 P3_f3% V3_e1 V3_f1 V3_e2 V3_f2 V3_e3 V3_f3]%--------------------------------------------------------%-------------------求雅可比矩阵的参数-----------------------%GeBf=G*e-B*f;%辅助计算函数GfBe=G*f+B*e;%辅助计算函数for i=1:2for j=1:3if i==jP_e(i,j)=-GeBf(i)-G(i,j)*e(i)-B(i,j)*f(i);P_f(i,j)=-GfBe(i)-G(i,j)*f(i)+B(i,j)*e(i);Q_e(i,j)=GfBe(i)+B(i,j)*e(i)-G(i,j)*f(i);Q_f(i,j)=-GeBf(i)+G(i,j)*e(i)+B(i,j)*f(i);elseP_e(i,j)=-G(i,j)*e(i)-B(i,j)*f(i);Q_f(i,j)=-P_e(i,j);P_f(i,j)=B(i,j)*e(i)-G(i,j)*f(i);Q_e(i,j)=P_f(i,j);endendendfor j=1:2P_e(3,j)=-G(3,j)*e(3)-B(3,j)*f(3);P_f(3,j)=B(3,j)*e(3)-G(3,j)*f(3);V_e(3,j)=0;V_f(3,j)=0;endP_e(3,3)=-GeBf(3)-G(3,3)*e(3)-B(3,3)*f(3);P_f(3,3)=-GfBe(3)-G(3,3)*f(3)+B(3,3)*e(3);V_e(3,3)=-2*e(3);V_f(3,3)=-2*f(3);%----------------------------------------------------%%--------------------求出雅克比矩阵-------------------% J=zeros(6,6);for i=1:3for j=1:3J(2*i-1,2*j-1)=P_e(i,j);J(2*i-1,2*j)=P_f(i,j);endendfor i=1:2for j=1:3J(2*i,2*j-1)=Q_e(i,j);J(2*i,2*j)=Q_f(i,j);endendfor j=1:3J(6,2*j-1)=V_e(3,j);J(6,2*j)=V_f(3,j);enddisp('雅可比矩阵为：')J%-----------------------------------------------------%%------------解修正方程并得出修正后的电压向量-----------% del_V=-J\del_W;clear i;for k=1:3V(k)=V(k)+del_V(2*k-1)+i*del_V(2*k);enddisp('节点电压为：')disp(V)%-----------------------------------------------------%%------------------计算节点不平衡量--------------------%e=real(V);f=imag(V);P=e.*(G*e-B*f)+f.*(G*f+B*e);Q=f.*(G*e-B*f)-e.*(G*f+B*e);del_W=[-0.30-P(1);-0.18-Q(1);-0.55-P(2);-0.13-Q(2);0.5-P(3);1.1^2-e(3)^2-f(3)^2];disp('节点不平衡量为：')disp(del_W)end%------------------------最终结论-------------------------%disp('最终各节点电压幅值为：')disp(abs(V))disp('最终各节点电压相角(度)为：')disp(180*angle(V)/pi)disp('最终各节点注入功率：')S=P+i*Qdisp('最终各节点注入有功功率为：')Pdisp('最终各节点注入无功功率为：')Q运行结果：Branch =1.00002.0000 0.1000 + 0.4000i 0 + 0.0153i 1.0000 03.0000 1.0000 0 + 0.3000i Inf 1.1000 1.00001.0000 4.0000 0.1200 + 0.5000i 0 + 0.0192i1.0000 02.0000 4.0000 0.0800 + 0.4000i 0 + 0.0141i 1.0000 0节点导纳矩阵Y =1.0421 - 8.2429i -0.5882 +2.3529i 0 +3.6667i -0.4539 + 1.8911i-0.5882 + 2.3529i 1.0690 - 4.7274i 0 -0.4808 + 2.4038i0 + 3.6667i 0 0 - 3.3333i 0-0.4539 + 1.8911i -0.4808 + 2.4038i 0 0.9346 - 4.2616i节点电压的实部：e =1.00001.00001.10001.0500节点电压的虚部：f =节点注入有功功率：Ps =-0.3000-0.55000.5000节点注入无功功率：Qs =-0.1800-0.1300迭代次数：1雅可比矩阵为：J =-1.0194 -8.3719 0.5882 2.3529 0 3.6667 -8.1138 1.0648 2.3529 -0.5882 3.6667 00.5882 2.3529 -1.0450 -4.8770 0 02.3529 -0.5882 -4.5778 1.0930 0 00 4.0333 0 0 0 -3.66670 0 0 0 -2.2000 0节点电压为：0.9935 - 0.0088i0.9763 - 0.1078i1.1000 + 0.1267i1.0500节点不平衡量为：-0.0013-0.0028-0.0135-0.05470.0030-0.0160迭代次数：2雅可比矩阵为：J =-0.8091 -8.3613 0.6052 2.3325 0.0324 3.6429 -7.9992 1.4071 2.3325 -0.6052 3.6429 -0.03240.8280 2.2338 -1.0190 -4.6364 0 02.2338 -0.8280 -4.3641 2.0878 0 0-0.4644 4.0333 0 0 -0.0324 -3.64290 0 0 0 -2.2000 -0.2533节点电压为：0.9847 - 0.0086i0.9590 - 0.1084i1.0924 + 0.1289i1.0500节点不平衡量为：-0.0000-0.0000-0.0003-0.00110.0001-0.0001迭代次数：3雅可比矩阵为：J =-0.7940 -8.2936 0.5995 2.3120 0.0315 3.6107 -7.9228 1.4000 2.3120 -0.5995 3.6107 -0.03150.8191 2.1927 -0.9865 -4.6144 0 02.1927 -0.8191 -4.2210 2.0885 0 0-0.4728 4.0056 0 0 -0.0315 -3.61070 0 0 0 -2.1849 -0.2579。