matlab潮流计算程序
MATpower潮流计算软件的使用
MATpower软件的使用一、MATpower软件的使用方法在MATLAB软件中的命令窗口输入runpf（‘程序名’）就可以通过MATpower已经编好的程序计算潮流，而函数runpf的参数是相应需计算潮流的数据文件。
数据文件主要用来定义和返回一下四个变量。
1、baseMVA baseMVA是一个标量，用来设置基准容量。
对于计算中采用有名值，可以根据实际情况设置，如设置100MV A；对于计算中采用标幺值，一般设置为12、bus bus变量是一个矩阵，用来设置电网中各节点参数，该矩阵内的参数如下：%% bus data%bus_i type Pd Qd Gs Bs area Vm Va baseKV zone Vmax Vmin其中，第1列参数即bus_i用来设置母线编号，范围为1~29997；第2列参数type用来设置母线类型，1为PQ节点，2为PV节点，3为平衡节点；第3列参数Pd用来设置母线注入负荷的有功功率；第4列参数Qd用来设置母线注入负荷的无功功率；第5列参数Gs用来设置与母线并联的电导；第6列参数Bs用来设置与母线并联的电纳；第7列参数area用来设置电网断面号，可设置范围为1~100，一般设置为1；第8列参数Vm用来设置母线电压的幅值初值；第9列参数Va用来设置母线电压的相角初值；第10列参数baseKV用来设置该母线的基准电压；第11列参数zone用来设置省耗分区号，可设置范围为1~999，一般设置为1；第12列参数Vmax用来设置工作时母线电压最高幅值；第13列参数Vmin用来设置工作时母线电压最低幅值。
3、gen gen变量是一个矩阵，用来设置接入电网的发电机参数，该矩阵的参数如下：%% generator data%bus Pg Qg Qmax Qmin Vg mBase status Pmax Pmin 其中，第1列参数bus用来设置接入发电机的母线编号；第2列参数Pg用来设置接入发电机的有功功率，注意功率输入的是有名值；第3列参数Qg用来设置接入发电机的无功功率；第4列参数Qmax用来设置接入发电机的无功功率的最大允许值；第5列参数Qmin用来设置接入发电机的无功功率的最小允许值；第6列Vg用来设置接入发电机的工作电压，注意输入的是标幺值；第7列mBase用来设置接入发电机的功率基准；第8列status用来设置发电机的工作状态，1表示投入运行，2表示投出运行；第9列Pmax用来设置接入发电机的无功功率的最大允许值；第10列参数Pmin用来设置接入发电机的无功功率的最小允许值。
基于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程序
0 0.343+0.21256i 1 0 0 2]
;%input('请输入各节点参数形成的矩阵： B2=');
Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
J(p,q)=X4;J(m,q)=X2; % X4=dQ/df X2=dp/df
elseif j1==i&j1~=isb %非平衡节点&对角元
X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/de
IT2=0;a=a+1;
for i=1:n
if i~=isb %非平衡节点
C(i)=0;D(i)=0;
for j1=1:n
C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)
J(m,q)=X2;
end
end
else
%=============== 下面是针对PV节点来求取Jacobi矩阵的元素 ===========
8 9 0.03512+0.08306i 0.13455i 1 0]
B2=[0 0 1.1 1.1 0 1;
0 0 1 0 0 2;
Байду номын сангаас%================= 求取Jacobi矩阵 ===================
for j1=1:n
if j1~=isb&j1~=i %非平衡节点&非对角元
X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/df
MATLAB电力系统PQ潮流计算程序设计
MATLAB电力系统PQ潮流计算程序设计1 绪论1．1潮流计算1.1.1 潮流计算概述电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。
在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性、可靠性和经济性。
此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。
所以潮流计算是研究电力系统的一种很重要也很基础的计算。
电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的随时监视及实时控制。
利用计算机进行电力系统潮流计算从50年代中期就已经开始。
在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。
对潮流计算的要求可以归纳为下面几点：（1）计算方法的可靠性或收敛性；（2）对计算机内存量的要求；（3）计算速度；（4）计算的方便性和灵活性。
电力系统潮流计算问题在数学上是一组多元非线性方程式求解问题，其解法都离不开迭代。
因此，对潮流计算方法，首先要求它能可靠地收敛，并给出正确答案。
由于电力系统结构及参数的一些特点，并且随着电力系统不断扩大，潮流计算方程式的阶数也越来越高，对这样的方程式并不是任何数学方法都能保证给出正确答案的。
这种情况成为促使电力系统计算人员不断寻求新的更可靠方法的重要因素。
在用数字计算机解电力系统潮流问题的开始阶段，普遍采取以节点导纳矩阵为基础的逐次代入法。
这个方法的原理比较简单，要求的数字计算机内存量比较低，适应50年代电子计算机制造水平和当时电力系统理论水平。
但它的收敛性较差，当系统规模变大时，迭代次数急剧上升，在计算中往往出现迭代不收敛的情况。
这就迫使电力系统计算人员转向以阻抗矩阵为基础的逐次代入法。
基于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辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。
此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。
现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法.高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。
lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。
(读者可以到MATPOWER主页下载MATPOWER4.1，然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算）一、高斯—赛德尔法潮流计算使用的程序：高斯—赛德法的具体使用方法读者可参考后面的实例,这里仅介绍各程序的编写格式：lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。
程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。
根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。
对于PV节点，如发电机节点,要提供一个无功功率限定值。
当给定电压过高或过低时，无功功率可能超出功率限定值。
在几次迭代之后(高斯—塞德尔迭代为10次），需要检查一次发电机节点的无功出力，如果接近限定值,电压幅值进行上下5%的调整,使得无功保持在限定值内。
lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。
并将这些参数放在名为linedata的文件中。
这个程序将阻抗转换为导纳,并得到节点导纳矩阵.busout：该程序以表格形式输出结果,节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。
潮流计算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);。
22节点潮流计算matlab仿真程序
【22节点潮流计算matlab仿真程序详解】引言主题指定的22节点潮流计算matlab仿真程序是电力系统领域中的重要内容。
在本文中，我将以深入、广泛的方式来探讨这一主题，帮助您全面了解和掌握相关知识。
一、22节点潮流计算的概念1. 22节点潮流计算的基本介绍22节点潮流计算是电力系统分析中的重要部分，它是对电力系统节点之间功率、电压、相角等参数进行计算和分析的过程。
通过22节点潮流计算，可以有效评估电力系统的稳定性和可靠性。
2. 22节点潮流计算的基本原理在进行22节点潮流计算时，需要考虑节点之间的功率平衡方程、电压平衡方程等基本原理。
这些原理是潮流计算的理论基础，也是实际仿真程序设计的重要依据。
二、22节点潮流计算matlab仿真程序的设计与实现1. matlab在电力系统仿真中的应用matlab作为一种功能强大的科学计算软件，在电力系统领域有着广泛的应用。
通过matlab，可以方便地进行潮流计算、稳定性分析、拓扑优化等工作。
2. 22节点潮流计算matlab仿真程序的设计要点在设计22节点潮流计算matlab仿真程序时，需要考虑程序的模块化、可扩展性、运行效率等要点。
通过合理的程序设计，可以提高仿真程序的稳定性和可靠性。
三、22节点潮流计算matlab仿真程序的应用实例1. 22节点潮流计算matlab仿真程序的基本框架通过一份完整的22节点潮流计算matlab仿真程序，来实际演示程序的结构和实现过程。
这样可以让读者更直观地理解程序的设计和应用。
2. 22节点潮流计算matlab仿真程序的实际应用案例以一个具体的电力系统实例为例，演示22节点潮流计算matlab仿真程序的实际应用过程。
这样可以让读者更清晰地了解程序在实际工程中的价值和作用。
结束语通过本文的深度探讨和广度展示，相信您已经对22节点潮流计算matlab仿真程序有了更全面的了解。
也希望您能够在实际工作中灵活运用所学知识，不断提升自己在电力系统领域的技术水平。
电力系统潮流计算的MATLAB辅助程序设计,潮流计算程序(精编文档).doc
【最新整理，下载后即可编辑】电力系统潮流计算的MATLAB辅助程序设计潮流计算，通常指负荷潮流，是电力系统分析和设计的主要组成部分，对系统规划、安全运行、经济调度和电力公司的功率交换非常重要。
此外，潮流计算还是其它电力系统分析的基础，比如暂态稳定，突发事件处理等。
现代电力系统潮流计算的方法主要:高斯法、牛顿法、快速解耦法和MATLAB的M语言编写的MATPOWER4.1，这里主要介绍高斯法、牛顿法和快速解耦法。
高斯法的程序是lfgauss，其与lfybus、busout和lineflow程序联合使用求解潮流功率。
lfybus、busout和lineflow程序也可与牛顿法的lfnewton程序和快速解耦法的decouple程序联合使用。
(读者可以到MATPOWER主页下载MATPOWER4.1,然后将其解压到MATLAB目录下，即可使用该软件进行潮流计算)一、高斯-赛德尔法潮流计算使用的程序：高斯-赛德法的具体使用方法读者可参考后面的实例，这里仅介绍各程序的编写格式:lfgauss：该程序是用高斯法对实际电力系统进行潮流计算，需要用到busdata和linedata两个文件。
程序设计为输入负荷和发电机的有功MW和无功Mvar，以及节点电压标幺值和相角的角度值。
根据所选复功率为基准值将负荷和发电机的功率转换为标幺值。
对于PV节点，如发电机节点，要提供一个无功功率限定值。
当给定电压过高或过低时，无功功率可能超出功率限定值。
在几次迭代之后(高斯-塞德尔迭代为10次)，需要检查一次发电机节点的无功出力，如果接近限定值，电压幅值进行上下5%的调整，使得无功保持在限定值内。
lfybus：这个程序需要输入线路参数、变压器参数以及变压器分接头参数。
并将这些参数放在名为linedata的文件中。
这个程序将阻抗转换为导纳，并得到节点导纳矩阵。
busout：该程序以表格形式输出结果，节点输出包括电压幅值和相角，发电机和负荷的有功和无功功率，以及并联电容器或电抗器的有功和无功功率。
matlab潮流计算程序
基于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初始化为0if 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初始化为0if 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.000001B1(支路参数矩阵)【1 2 0.04+0.25i 0.5i 1 0;1 3 0.1+0.35i 0 1 0;2 3 0.08+0.30i 0.5i 1 0;4 2 0.015i 0 1.05 1;5 3 0.03i 0 1.05 1】B2(节点参数矩阵)【0 -1.6-0.8i 1 0 0 2;0 -2-1i 1 0 0 2;0 -3.7-1.3i 1 0 0 2;0 5+0i 1.05 1.05 0 3;0 0 1.05 1.05 0 1】X(节点号和对地参数)【1 0;2 0;3 0;4 0;5 0】亲爱的朋友，上文已完，为感谢你的阅读，特加送另一篇范文，如果下文你不需要，可以下载后编辑删除，谢谢！矿井水灾事故专项应急演练方案1 应急演练目的、意义和目标1.1应急演练目的①评估我矿水灾事故的应急准备状态，发现并修改我矿水灾事故专项应急预案和执行程序中存在的缺陷和不足；②评估我矿在发生水灾事故时的应急能力，识别处理水灾事故的资源需求，澄清相关单位和人员的应急职责，改善水灾事故应急救援中的组织协调问题；③检验应急响应人员对水灾事故应急预案及执行程序的了解程度和实际操作技能；同时，通过调整演练难度，进一步培训和提高应急响应人员的业务素质和能力；④提高全员安全意识。
完整版本基于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 时同样 ):' ); fori=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实现潮流计算程序
程序代码如下: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电力系统潮流计算
电力系统潮流计算是电力系统分析和运行中的重要环节。
潮流计算主要用于确定电力系统中各个节点的电压、功率和潮流方向，以便进行功率平衡、电压稳定和线路负荷等方面的分析和评估。
MATLAB作为一种强大的数学建模和仿真工具，被广泛应用于电力系统潮流计算的研究和实际工程中。
本文将介绍MATLAB在电力系统潮流计算中的应用，包括算法原理、建模步骤和实例分析等内容。
一、潮流计算的基本原理潮流计算是指在给定电力网拓扑结构、负荷信息和发电机功率的情况下，通过迭代计算求解节点电压的复数值，以确定各节点的电压幅值和相角，进而计算各支路和各节点上的有功和无功功率。
潮流计算的基本原理是基于电力系统的潮流方程和节点功率平衡等基本理论，通过建立节点电压的复数方程组，利用迭代计算方法求解该方程组，从而得到节点的电压和功率信息。
二、MATLAB在潮流计算中的应用MATLAB作为一种功能强大的数学建模和仿真工具，具有丰富的数学计算函数和图形显示功能，适合于电力系统潮流计算的建模和仿真。
在MATLAB环境下，可以利用其矩阵运算、方程求解和数据可视化等功能，实现电力系统潮流计算的数学模型和算法的实现。
下面将介绍MATLAB在电力系统潮流计算中的具体应用步骤。
1. 建立电力系统潮流计算的数学模型在MATLAB环境下，首先需要建立电力系统潮流计算的数学模型，包括节点电压方程、支路潮流方程、节点功率平衡方程等。
利用MATLAB的矩阵运算和符号计算工具，可以将电力系统的节点和支路参数、负荷信息、发电机功率等数据表示为矩阵形式，建立电力系统潮流计算的数学模型。
2. 编写潮流计算的求解算法在建立电力系统潮流计算的数学模型后，需要编写潮流计算的求解算法。
在MATLAB环境下，可以利用其丰富的数学计算函数和优化工具，实现潮流计算的迭代求解算法，包括高斯-赛德尔迭代法、牛顿-拉夫逊迭代法等。
通过编写求解算法，可以实现电力系统潮流计算的数值求解过程。
3. 进行潮流计算的仿真实验在完成潮流计算的求解算法后，可以利用MATLAB进行潮流计算的仿真实验。
电力网潮流计算的MATLAB总程序(课程组自编)
电力网潮流计算的MATLAB总程序cleardisp('此程序适用可带有双绕组变压器的任意节点数目的通用网络的潮流计算')n=input('请输入节点数目n; n=');m=input('请输入PQ节点数目m; m=');disp('请将节点按如下规则排序,n个节点中前m个节点为PQ节点,第m+1至n-1个节点为PV节点,第n个节点为平衡节点.')t1=menu('请选择系统的给定参数形式:有名值选1,标幺值选2;','1','2')switch t1case 1disp('给定参数为有名值')t2=menu('选择所求系统是否含有变压器;是选1,否选2.','1','2')switch t2case 1disp('系统含有变压器')SB=input('请输入基准功率(MVA);SB=');disp('请输入系统共有几段电压等级Nv,单位:千伏.')Nv=input('Nv=');disp('为了输入快捷,在以下程序中要求输入"数组"或"矩阵"时请按如下方法输入')disp('输入的数据必须方在方括号内,一行中有多个元素时用空格或逗号隔开,行与行之间用分号或回车隔开.')disp('请为段排序后按序输入各段额定电压VN,单位:千伏;数组形式,如[220,110,10]') VN=input('VN=');[a,c1]=size(VN);while a~=1|c1~=Nvdisp('请为段排序后正确输入各段额定电压VN,单位:千伏;数组形式,如[220,110,10]') VN=input('VN=');[a,c1]=size(VN);enddisp('请按序对应输入各段电压的基准电压VB,数组形式,单位:千伏;如[230,115,10.5]') VB=input('VB=');[a,c2]=size(VB);while a~=1| c2~=Nvdisp('请正确输入基准电压VB,数组形式,单位:千伏;如[230,115,10.5]')VB=input('VB=');[a,c2]=size(VB);endb=1;for a=1:Nvdisp('请输入运行在第'),disp(a),disp('段中的节点号.数组形式,如:[2,3,5]') Nvn=input('请输入:'); %运行在第I,II,III,...段的节点编号[c1,c2]=size(Nvn);for b=1:c2NVN(a,b)=Nvn(b);endendNVNYT=zeros(n);for a=1:ndisp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a) disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);[x1,x2]=find(a==NVN);[x3,x4]=find(b==NVN);if x1==x3disp('请输入节点'),disp(a),disp('与),disp(b),disp('之间的线路阻抗Zij单位:欧姆,电纳B,单位:S')disp('注:可以输入算式,如:(a+j*b)*c/d;复数虚部必须输入j*a,a 代表数值;电纳只输入数值即可')Zij(a,b)=input('Zij=');BB(a,b)=input('B=');while BB(a,b)^2<0disp('输入错误,此值不必输入虚数单位')BB(a,b)=input('B=');endZij(a,b)=Zij(a,b)*SB/VB(x1)^2; % 线路阻抗标幺值BB(a,b)=BB(a,b)*VB(x1)^2/SB; % 对地导纳标幺值Y(a,b)=-1/Zij(a,b);Y(b,a)=Y(a,b); %互导纳y0(a,b)=j*BB(a,b)/2; %对地导纳矩阵元素(虚数)elsedisp('节点'),disp(a),disp('与'),disp(b),disp('不在同一段电压等级,通过变压器相连;请输入它们之间的变压器参数')ST=input('额定容量,单位:兆伏安;S=');Vs100=input('短路电压百分数Vs%=');disp('请输入变压器两侧额定电压VTN;单位:千伏;数组形式,如:[10.5 121]')VTN=input('VTN=');while [1 2]~=size(VTN)disp('输入错误,请输入变压器两侧额定电压,单位:千伏;数组形式,如:[10.5 121]')VTN=input('VTN=');endif VN(x1)<VN(x3)if VTN(1)<VTN(2)ZT(a,b)=j*Vs100*VTN(2)^2*SB/(100*ST*VB(x3)^2); %变压器等效阻抗标幺值kT(b,a)=VTN(2)/VTN(1)/(VB(x3)/VB(x1)); %变压器变比标幺值elseZT(a,b)=j*Vs100*VTN(1)^2*SB/(100*ST*VB(x3)^2) ;%变压器等效阻抗标幺值kT(b,a)=VTN(1)/VTN(2)/(VB(x3)/VB(x1));%变压器变比标幺值 endkT(a,b)=kT(b,a)y0(b,a)=(1-kT(a,b))/ZT(a,b);y0(a,b)=kT(a,b)*(kT(a,b)-1)/ZT(a,b);elseif VTN(1)<VTN(2)ZT(a,b)=j*Vs100*VTN(2)^2*SB/(100*ST*VB(x1)^2);kT(a,b)=VTN(2)/VTN(1)/(VB(x1)/VB(x3));elseZT(a,b)=Vs100*VTN(1)^2*SB/(100*ST*VB(x1)^2);kT(a,b)=VTN(1)/VTN(2)/(VB(x1)/VB(x3));endy0(a,b)=(1-kT(a,b))/ZT(a,b);y0(b,a)=kT(a,b)*(kT(a,b)-1)/ZT(a,b);endY(a,b)=-kT(a,b)/ZT(a,b);Y(b,a)=Y(a,b);YT(a,b)=Y(a,b);YT(b,a)=YT(a,b);endendendY(a,a)=sum(y0(a,:))-sum(Y(a,:));endcase 2disp('系统不含变压器')SB=0;disp('为了输入快捷,在以下程序中要求输入"数组"或"矩阵"时请按如下方法输入')disp('输入的数据必须方在方括号内,一行中有多个元素时用空格或逗号隔开,行与行之间用分号或回车隔开.') ;for a=1:ndisp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a) disp('的节点号,数组形式如"[3 5 6]"或"[3,5,6]",若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);disp('请输入节点'),disp(a),disp('与'),disp(b),disp('之间的线路阻抗Zij单位:欧姆,电纳B,单位:S')disp('注:可以输入算式,如:(r+j*x)*l;复数虚部必须输入j*a,a 代表数值;电纳只输入数值即可')Zij(a,b)=input('Zij=');BB(a,b)=input('B=');while BB(a,b)^2<0disp('输入错误,此值不必输入虚数单位')BB(a,b)=input('B=');endY(a,b)=-1/Zij(a,b);Y(b,a)=Y(a,b);y0(a,b)=j*BB(a,b)/2;endY(a,a)=sum(y0(a,:))-sum(Y(a,:))endendendcase 2disp('给定参数为标幺值')SB=0;t3=menu('请选择所求系统是否含有变压器,是选1,否选2', '1', '2');switch t3case 1 %对含变支路的输入和处理disp('系统含有变压器')nT=input('请输入含有变压器支路数目nT;nT=');YT=zeros(n); %产生一个n*n的零矩阵for a=1:nTdisp('请输入第'),disp(a),disp('条含有变压器支路的节点编号,其形式如"[x y]"或"[x,y]"(元素间用空格或逗号隔开)')disp('x为低压侧节点,y为高压侧节点')Tij=input('Tij=');x1=Tij(1);x2=Tij(2);disp('请输入节点'),disp(x1),disp('与'),disp(x2),disp('之间的阻抗Z(虚部j*a,a 代表数值)及变压器变比k')Zij(x1,x2)=input('Z=');kT=input('k=');t4=menu('请选择变压器阻抗Z是折算到高压侧的值还是低压侧,高压侧选1,低压侧选2','1','2');switch t4case 1Y(x1,x2)=-kT/Zij(x1,x2); %%此为含变支路互导纳Y(x2,x1)=Y(x1,x2); %对称性y0(x1,x2)=kT*(kT-1)/Zij(x1,x2);y0(x2,x1)=(1-kT)/Zij(x1,x2);case 2Y(x1,x2)=-1/(kT*Zij(x1,x2)); %%此为含变支路互导纳Y(x2,x1)=Y(x1,x2); %对称性y0(x1,x2)=(kT-1)/kT*Zij(x1,x2);y0(x2,x1)=(1-kT)/(kT^2*Zij(x1,x2));endendcase 2 %不含变disp('系统不含有变压器')YT=zeros(n); %产生一个n*n的零矩阵endfor a=1:n %导纳矩阵计算主程序disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",(变压器支路不必再输);若无则输入0')F=input('请输入:');while (min(F)<=a|max(F)>n)&F~=0disp('请正确输入')disp('请输入与节点'),disp(a),disp('有支路直接相连且编号大于'),disp(a)disp('的节点号,其形式如(冒号之内)"[3 5 6]"或"[3,5,6]",(变压器支路不必再输);若无则输入0')F=input('请输入:');endif F~=0[f1,f]=size(F); %求F中有几个元素for c=1:fb=F(c);disp('请输入节点'),disp(a),disp('与'),disp(b),,disp('之间的线路阻抗Zij,线路对地支路导纳y0(虚部j*a,a 代表数值)')Zij(a,b)=input('Zij=');y0(a,b)=input('y0=');while y0(a,b)^2>0disp('输入错误,请重新输入;虚数须加虚数单位"j*"')y0(a,b)=input('y0=');endy0(b,a)=y0(a,b);Y(a,b)=-1/Zij(a,b); Y(b,a)=Y(a,b); %互导纳endendY(a,a)=sum(y0(a,:))-sum(Y(a,:));endendYG=real(Y);B=imag(Y);disp('以数组形式(形如"[50 1.5 75]")按节点顺序输入各节点注入的有功功率,无功功率和节点电压') disp('请按节点顺序输入PQ,PV节点的有功功率,可以输入算式形如[3+2,5+8]')Ps=input('Pis=');[x1,x2]=size(Ps);while x1~=1 | x2~=n-1disp('输入错误!请按节点顺序输入PQ,PV节点的有功功率,也可以输入算式形如[3+2,5+8]')Ps=input('Pis=');[x1,x2]=size(Ps);enddisp('请按节点顺序输入PQ节点的无功功率,也可以输入算式形如[3+2,5+8]')Qs=input('Qis=');[x1,x2]=size(Qs);while x1~=1 | x2~=mdisp('输入错误!请按节点顺序输入PQ节点的无功功率,也可以输入算式形如[3+2,5+8]')Qs=input('Qis=');[x1,x2]=size(Qs);enddisp('请按节点顺序输入PV节点及平衡节点电压幅值(从第m+1节点开始)')V=input('V=');[x1,x2]=size(V);while x1~=1 | x2~=n-mdisp('输入错误!请按节点顺序输入PV节点及平衡节点电压幅值')V=input('V=');[x1,x2]=size(V);enddisp('请输入平衡节点电压相角')r=input('r=');[x1,x2]=size(r);while x1~=1 | x2~=1disp('输入错误!请输入平衡节点电压相角')r=input('r=');[x1,x2]=size(r);endt10=menu('是否考虑PV节点无功越界?,是选1,否选2','1','2');switch t10case 1Qv=input('请按节点顺序输入PV节点可提供的最大无功功率,Qmax=');case 2Qv=inf;endswitch t1case 1switch t2case 1Ps=Ps/SB;Qs=Qs/SB;Vsd=ones(1,m);V=[Vsd,V];rsd=zeros(1,n-1);r=[rsd,r];for a=m+1:n-1[x1,x2]=find(a==NVN);V(a)=V(a)/VB(x1);endQv=Qv/SB;case 2[c1,c2]=size(V);Vsd=V(c2)*ones(1,m);V=[Vsd,V];rsd=zeros(1,n-1);r=[rsd,r];endcase 2Vsd=ones(1,m);V=[Vsd,V]rsd=zeros(1,n-1);r=[rsd,r]endk=0;jdgd=input('输入要求精度;jd=');jd=inf;while jd>jdgdfor a=1:n-1for b=1:nPn(b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));Qn(b)=V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endPi(a)=sum(Pn(:)); %PiQH(a)=sum(Qn(:));if a<=mQi(a)=sum(Qn(:)); %求QiendendDPi=Ps-Pi;DQi=Qs-Qi; %求得DPifor a=1:n-1for b=1:n-1if a~=bH(a,b)=-V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));if b<=mN(a,b)=-V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));endif a<=mK(a,b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));if b<=mL(a,b)=-V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endendelseH(a,a)=V(a)^2*B(a,a)+QH(a);if b<=mN(a,a)=-V(a)^2*G(a,a)-Pi(a);endif a<=mK(a,a)=V(a)^2*G(a,a)-Pi(a);if b<=mL(a,a)=V(a)^2*B(a,a)-QH(a);endendendendendfor a=1:n-1for b=1:n-1if a<=m & b<=mJJ(2*a-1,2*b-1)=H(a,b);JJ(2*a-1,2*b)=N(a,b);JJ(2*a,2*b-1)=K(a,b); JJ(2*a,2*b)=L(a,b);endif a<=m & b>mJJ(2*a-1,m+b)=H(a,b);JJ(2*a,m+b)=K(a,b);endif a>m & b<=mJJ(m+a,2*b-1)=H(a,b);JJ(m+a,2*b)=N(a,b);endif a>m & b>mJJ(m+a,m+b)=H(a,b);endendendJJfor a=1:n-1if a<=mDPQ(2*a-1)=DPi(a);DPQ(2*a)=DQi(a);elseDPQ(a+m)=DPi(a);endendDPQDrV=-inv(JJ)*(DPQ');for a=1:n-1if a<=mDr(a)=DrV(2*a-1);DVV(a)=DrV(2*a);elseDr(a)=DrV(a+m);endendDrj=Dr/pi*180for a=1:mVd2(a,a)=V(a);endDV=Vd2*(DVV')for a=1:n-1r(a)=r(a)+Dr(a);if a<=mV(a)=V(a)+DV(a);endendjd=max(abs(DPQ))switch t10case 1for a=m+1:n-1if Qv(a-m)<QH(a)disp(a),disp('节点无功越限')pvwgyj=1;m=m+1;for b=1:nif Y(a,b)~=0Yc(a,b)=Y(a,b);Y(a,b)=Y(m,b);Y(m,b)=Yc(a,b);Y(b,a)=Y(a,b);Y(b,m)=Y(m,b);endendQ(a)=Q(m);Qs(m)=Qv(a);Ps(m)=P(a);P(a)=P(m)Vc(a)=V(a);V(a)=V(m);V(m)=Vc(a);rc(a)=r(a);r(a)=r(m);r(m)=rc(a) elsepvwgyj=0endendendk=k+1;k=k+1;enddisp('迭代次数')k=k-1format short edisp('节点功率不平衡量')DPiDQiformat shortdisp('各节点电压相角; 单位:度;')rd=r/pi*180for a=1:nfor b=1:nPn(b)=V(a)*V(b)*(G(a,b)*cos(r(a)-r(b))+B(a,b)*sin(r(a)-r(b)));Qn(b)=V(a)*V(b)*(G(a,b)*sin(r(a)-r(b))-B(a,b)*cos(r(a)-r(b)));endPN=sum(Pn(:)); %P(n)QN=sum(Qn(:)); %Q(n)endPi(n)=PN;QH(n)=QN;Sn=Pi+QH*j;for a=1:nU(a)=V(a)*(cos(r(a))+j*sin(r(a)));endfor a=1:nfor b=1:nS(a,b)=V(a)^2*conj(y0(a,b))+U(a)*(conj(U(a))-conj(U(b)))*conj(-Y(a,b));endendif SB==0disp('各节点电压幅值')Vdisp('各节点的注入功率')Sndisp('各支路的功率分布')Selsefor a=1:n[x1,x2]=find(a==NVN);V(a)=V(a)*VB(x1);enddisp('各节点电压幅值,单位:千伏')Vdisp('各节点的注入功率,单位:兆伏安')Sn=Sn/SBdisp('各支路的功率分布,单位:兆伏安')S=S/SBend。
潮流计算(matlab)实例计算
潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。
2.在给定的电力网络上画出等值电路图。
3.运用计算机进行潮流计算。
4.编写设计说明书。
一、设计原理1．牛顿-拉夫逊原理牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。
为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
牛顿—拉夫逊迭代法的一般步骤：（1）形成各节点导纳矩阵Y。
（2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。
（3）计算各个节点的功率不平衡量。
（4）根据收敛条件判断是否满足，若不满足则向下进行。
（5）计算雅可比矩阵中的各元素。
（6）修正方程式个节点电压（7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。
（8）计算平衡节点输出功率和各线路功率2．网络节点的优化１）静态地按最少出线支路数编号这种方法由称为静态优化法。