潮流计算的计算机方法

第四章潮流计算的计算机算法第一节概述潮流计算是电力系统最基本、最常用的计算。

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

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

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

电力系统潮流计算分为离线潮流计算和在线潮流计算。

前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。

本章主要讨论离线潮流计算问题，它的基本算法同样适用于在线潮流计算。

潮流计算在数学上是多元非线性方程组的求解问题，求解的方法有很多种。

自从五十年代计算机应用于电力系统以来，当时求解潮流的方法是以节点导纳矩阵为基础的逐次代入法（导纳法），后来为解决导纳法的收敛性较差的问题，出现了以阻抗矩阵为基础的逐次代入法（阻抗法）。

到六十年代，针对阻抗法占用计算机内存大的问题又出现了分块阻抗法及牛顿－拉夫逊（Newton－Raphson）法。

Newton —Raphson法是数学上解非线形方程式的有效方法，有较好的收敛性。

将N－R法用于潮流计算是以导纳矩阵为基础的，由于利用了导纳矩阵的对称性、稀疏性及节点编号顺序优化等技巧，使N－R法在收敛性、占用内存、计算速度方面的优点都超过了阻抗法，成为六十年代末期以后普遍采用的方法。

同时国内外广泛研究了诸如非线形规划法、直流法、交流法等各种不同的潮流计算方法。

七十年代以来，又涌现出了更新的潮流计算方法。

其中有1974年由B、Stott、O、Alsac 提出的快速分解法以及1978年由岩本伸一等提出的保留非线性的高129速潮流计算法。

其中快速分解法（Fast decoupled load flow）从1975年开始已在国内使用，并习惯称之为PQ分解法。

由于PQ分解法在计算速度上大大超过N－R法，不但能应用于离线潮流计算，而且也能应用于在线潮流计算。

由于本人参加我们电气学院的电气小课堂，主讲的是计算机算法计算潮流这章，所以潜心玩了一个星期，下面整理给大家分享下。

本人一个星期以来的汗水，弄清楚了计算机算法计算潮流的基础，如果有什么不懂的可以发信息到邮箱：zenghao616@接下来开始弄潮流的优化问题，吼吼！电力系统的潮流计算的计算机算法：以MATLAB为环境这里理论不做过多介绍，推荐一本专门讲解电力系统分析的计算机算法的书籍---------《电力系统分析的计算机算法》—邱晓燕、刘天琪编著。

这里以这本书上的例题【2-1】说明计算机算法计算的过程，分别是牛顿拉弗逊算法的直角坐标和极坐标算法、P-Q分解算法。

主要是简单的网络的潮流计算，其实简单网络计算和大型网络计算并无本质区别，代码里面只需要修改循环迭代的N即可，这里旨在弄清计算机算法计算潮流的本质。

代码均有详细的注释.其中简单的高斯赛德尔迭代法是以我们的电稳教材为例子讲，其实都差不多，只要把导纳矩阵Y给你，节点的编号和分类给你，就可以进行计算了，不必要找到原始的电气接线图。

理论不多说，直接上代码：简单的高斯赛德尔迭代法：这里我们只是迭代算出各个节点的电压值，支路功率并没有计算。

S_ij=P_ij+Q_ij=V_i(V_i* - V_j*) * y_ij*可以计算出各个线路的功率在显示最终电压幅角的时候注意在MATLAB里面默认的是弧度的形式，需要转化成角度显示。

clear;clc;%电稳书Page 102 例题3-5%计算网络的潮流分布 --- 高斯-赛德尔算法%其中节点1是平衡节点%节点2、3是PV节点，其余是PQ节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；y=zeros(5,5);y(1,2)=1/(0.0194+0.0592*1i);y(1,5)=1/(0.054+0.223*1i);y(2,3)=1/(0.04699+0.198*1i);y(2,4)=1/(0.0581+0.1763*1i);%由于电路网络的互易性，导纳矩阵为对称的矩阵for i=1:1:5for j=1:1:5y(j,i)=y(i,j);endend%节点导纳矩阵的形成Y=zeros(5,5);%求互导纳for i=1:1:5for j=1:1:5if i~=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:1:5%这句话是说将y矩阵的第i行的所有元素相加，得到自导纳的值 Y(i,i)=sum(y(i,:));end%上面求得的自导纳不包含该节点的对地导纳数值，需要加上Y(2,2)=Y(2,2)+0.067*1i;Y(3,3)=Y(3,3)+0.022*1i;Y(4,4)=Y(4,4)+0.0187*1i;Y(5,5)=Y(5,5)+0.0246*1i;%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);Qc2=0;Qc3=0;%原始节点功率%这里电源功率为正，负荷功率为负S(1)=0;S(2)=-0.217-0.121*1i+Qc2*1i;S(3)=-0.749-0.19*1i+Qc3*1i;S(4)=-0.658+0.039*1i;S(5)=-0.076-0.016*1i;%节点功率的P QP = real(S);Q = imag(S);%下面是两个PV节点的无功初始值Q(2) = 0;Q(3) = 0;U=ones(5,1); %1列5行的‘1’矩阵%节点电压初始值U(1)=1.06;U(2)=1.045;U(3)=1.01;U_reg=U;Sum_YU0=0;%中间变量Sum_YU1=0;%中间变量for cont=1:1:6 %这里的cont是迭代次数for i=2:1:5for j=1:1:iif i~=jSum_YU0 = Sum_YU0 + Y(i,j)*U_reg(j);endendfor j=i+1:1:5Sum_YU1 = Sum_YU1 + Y(i,j)*U(j);endU(i)=( (P(i)-Q(i)*1i ) / conj(U(i)) - Sum_YU0 - Sum_YU1 ) / Y(i,i); U_reg(i)=U(i);%PV节点计算%下面是把求出的U2、U3只保留其相位，幅值不变if i==2angle_U2 = angle(U(2));U(2)=1.045*cos(angle_U2)+1.045*sin(angle_U2)*1i;Q(2)=imag( U(2)*( conj(Sum_YU0) + conj(Sum_YU1) + conj(Y(2,2)*U(2)) ) );endif i==3angle_U3 = angle(U(3));U(3)=1.01*cos(angle_U3)+1.01*sin(angle_U3)*1i;Q(3)=imag( U(3)*( conj(Sum_YU0) + conj(Sum_YU1) + conj(Y(3,3)*U(3)) ) );end% 下面做越界检查%if Q(4)>Q_Max% Q(4) = Q_Max;%end%if Q(4)<Q_Min% Q(4) = Q_Min;%end%下面可以做PV节点收敛判断Sum_YU0 = 0;Sum_YU1 = 0;endend%节点注入无功，流入为正，流出为负Qc2=Q(2)+0.121-1.045^2 * 0.067;Qc3=Q(3)+0.19-1.01^2 * 0.022;%电压幅值和相角angle_U=angle(U)*180/pi;U=abs(U);S_Line=zeros(5,5);%计算平衡节点功率S_BalanceNode=0;for j=1:1:5S_BalanceNode = S_BalanceNode + U(1) * conj(Y(1,j)*U(j));end%下面由上面算出的电压值求线路的功率%这里计算出来的线路功率的有功、无功%for i=1:1:5% for j=i:1:5% if i~=j% S_Line(i,j)=U(i)*( conj(U(i))-conj(U(j)) ) * conj(y(i,j));% end% if i==2% %S_Line(2,j)=S_Line(2,j)+U(2)*conj(0.067*1i);% end% if i==3% %S_Line(3,j)=S_Line(3,j)+U(3)*conj(0.022*1i);% end% end%end计算网络的潮流分布 ---- Newton算法(直角坐标)clear;clc;%电稳书Page 102 例题3-5%计算网络的潮流分布 ---- Newton算法(直角坐标)%其中节点1是平衡节点%节点2、3是PV节点，其余是PQ节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；y=zeros(5,5);y(1,2)=1/(0.0194+0.0592*1i);y(1,5)=1/(0.054+0.223*1i);y(2,3)=1/(0.04699+0.198*1i);y(2,4)=1/(0.0581+0.1763*1i);%由于电路网络的互易性，导纳矩阵为对称的矩阵for i=1:1:5for j=1:1:5y(j,i)=y(i,j);endend%节点导纳矩阵的形成Y=zeros(5,5);%求互导纳for i=1:1:5for j=1:1:5if i~=jY(i,j)=-y(i,j);endendend%求自导纳for i=1:1:5%这句话是说将y矩阵的第i行的所有元素相加，得到自导纳的值Y(i,i)=sum(y(i,:));end%上面求得的自导纳不包含该节点的对地导纳数值，需要加上Y(2,2)=Y(2,2)+0.067*1i;Y(3,3)=Y(3,3)+0.022*1i;Y(4,4)=Y(4,4)+0.0187*1i;Y(5,5)=Y(5,5)+0.0246*1i;%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);%节点2、3需补偿的无功Qc2=0;Qc3=0;%原始节点功率%这里电源功率为正，负荷功率为负S(1)=0;S(2)=-0.217-0.121*1i+Qc2*1i;S(3)=-0.749-0.19*1i+Qc3*1i;S(4)=-0.658+0.039*1i;S(5)=-0.076-0.016*1i;%节点功率的P QP = real(S);Q = imag(S);%下面是两个PV节点的无功初始值Q(2) = 0;Q(3) = 0;%给点电压初始值e=[1.06,1.045,1.01,1,1];f=[0,0,0,0,0];U=e+f*1i;delta_U=zeros(1,5);delta_P=zeros(1,5);delta_Q=zeros(1,5);delta_PQV=ones(8,1);Sum_GB1=0;Sum_GB2=0;cont=0;while max(delta_PQV > 1e-6),cont=cont+1;%for cont=1:1:3%下面开始计算delta_P/delta_Q/delta_Ufor i=2:1:5for j=1:1:5Sum_GB1=Sum_GB1 + ( G(i,j)*e(j) - B(i,j)*f(j) );Sum_GB2=Sum_GB2 + ( G(i,j)*f(j) + B(i,j)*e(j) );enddelta_P(i)=P(i)-e(i)*Sum_GB1-f(i)*Sum_GB2;if i~=2 && i~=3 %不为节点2,3则计算无功delta_Q(i)=Q(i)-f(i)*Sum_GB1+e(i)*Sum_GB2;endif i==2 || i==3 %这里计算delta_U的值，始终为零delta_U(i)=U(i)^2-( e(i)^2 + f(i)^2 );endSum_GB1=0;Sum_GB2=0;end%___________________________________%%下面计算雅克比矩阵J=zeros(8,8);for ii=2:1:5i=ii-1;for j=1:1:5Sum_GB1=Sum_GB1 + ( G(ii,j)*e(j) - B(ii,j)*f(j) );Sum_GB2=Sum_GB2 + ( G(ii,j)*f(j) + B(ii,j)*e(j) );endfor jj=2:1:5j=jj-1;if ii~=2 && ii~=3 %PQ节点if ii==jjJ(2*i-1,2*i-1)=-Sum_GB1-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);J(2*i-1,2*i)=-Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i-1)=Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i)=-Sum_GB1+G(ii,ii)*e(ii)+B(ii,ii)*f(ii);elseJ(2*i-1,2*j-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii)); J(2*i-1,2*j)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);J(2*i,2*j-1)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);J(2*i,2*j)=(G(ii,jj)*e(ii)+B(ii,jj)*f(ii));endelse%PV节点if ii==jjJ(2*i-1,2*i-1)=-Sum_GB1-G(ii,ii)*e(ii)-B(ii,ii)*f(ii);J(2*i-1,2*i)=-Sum_GB2+B(ii,ii)*e(ii)-G(ii,ii)*f(ii);J(2*i,2*i-1)=-2*e(ii);J(2*i,2*i)=-2*f(ii);elseJ(2*i-1,2*j-1)=-(G(ii,jj)*e(ii)+B(ii,jj)*f(ii)); J(2*i-1,2*j)=B(ii,jj)*e(ii)-G(ii,jj)*f(ii);J(2*i,2*j-1)=0;J(2*i,2*j)=0;endendendSum_GB1=0;Sum_GB2=0;end%在求解修正方程之前建议把delta_P和delta_Q,delta_U全部放在一个矩阵delta_PQV=[delta_P(2);delta_U(2);delta_P(3);delta_U(3);delta_P(4) ;delta_Q(4);delta_P(5);delta_Q(5)];%下面求解修正方程;注意矩阵运算时候的左除和右除的区别delta_ef=-J\delta_PQV;%下面修正各个节点的电压for i=2:1:5e(i)=e(i)+delta_ef(2*(i-1)-1);f(i)=f(i)+delta_ef(2*(i-1));end%到这里第一轮迭代完成end%电压幅值和相角U=e+f*1i;angle_U=angle(U)*180/pi;%节点注入无功，流入为正，流出为负Sum_YU=0;for i=2:1:3for j=1:1:5Sum_YU = Sum_YU + Y(i,j)*U(j);endQ(i)=imag( U(i)*conj( Sum_YU ) );Sum_YU=0;endQc2=Q(2)+0.121-1.045^2 * 0.067;Qc3=Q(3)+0.19-1.01^2 * 0.022;U=abs(U);disp(['Iteration times : ' num2str(cont)]);%显示最终的迭代次数牛顿算法求解潮流 (极坐标):clear;clc;%牛顿算法求解潮流 (极坐标)%计算网络的潮流分布%其中节点5是平衡节点%节点1、2、3是PQ节点，节点4是PV节点% 如果节点有对地导纳支路%需将对地导纳支路算到自导纳里面%------------------------------------------------%%输入原始数据，每条支路的导纳数值，包括自导和互导纳；Y=[0.8381-3.7899*1i,-0.4044+1.6203*1i,0,0,-0.4337+2.2586*1i;...-0.4044+1.6203*1i,0.7769-3.3970*1i,-0.3726+1.8557*1i,0,0;...0,-0.3726+1.8557*1i,1.1428-7.0210*1i,-0.5224+4.1792*1i,-0.2739+1. 2670*1i;...0,0,-0.5224+4.1792*1i,0.5499-4.3591*1i,0;...-0.4337+2.2586*1i,0,-0.2739+1.2670*1i,0,0.7077-3.4437*1i];%导纳矩阵的实部和虚部G = real(Y);B = imag(Y);%给点电压初始值U = [1,1,1,1,1.05];angle_U=[0,0,0,0,0];%for i=1:1:5% U(i)=U_abs(i)*cos(angle_U(i))+U_abs(i)*sin(angle_U(i))*1i;%end%原始节点功率%这里电源功率为正，负荷功率为负%下面给点PQ PV节点功率值S=[-0.22-0.14*1i,-0.18-0.09*1i,-0.27-0.13*1i,0.35,0];%节点功率的P QP = real(S);Q = imag(S);%下面是PV节点的无功初始值Q(4) = 0;delta_P=zeros(1,5);delta_Q=zeros(1,5);%delta_angleU=zeros(1,4);%delta_absU=zeros(1,4);delta_PQ=ones(8,1);Sum_GB1=0;Sum_GB2=0;cont=0;%最外层循环，cont代表迭代的次数，这里可以用约束条件来代替%for cont=1:1:4while max(delta_PQ)>1e-6,%下面计算delta_P/delta_Q/delta_Ucont=cont+1;for i=1:1:4for j=1:1:5Sum_GB1=Sum_GB1 + U(j)*( G(i,j)*cos(angle_U(i)-angle_U(j)) + B(i,j)*sin(angle_U(i)-angle_U(j)) );Sum_GB2=Sum_GB2 + U(j)*( G(i,j)*sin(angle_U(i)-angle_U(j)) - B(i,j)*cos(angle_U(i)-angle_U(j)) );enddelta_P(i)=P(i)-U(i)*Sum_GB1;if i~=4 %不为节点四则计算无功delta_Q(i)=Q(i)-U(i)*Sum_GB2;endSum_GB1=0;Sum_GB2=0;end%_______________________________________________________%%下面计算雅克比矩阵J=zeros(7,7);for ii=1:1:4for jj=1:1:4if ii ~= 4 %PQ节点if ii==jjJ(2*ii-1,2*ii-1)=U(ii)^2*B(ii,ii)+Q(ii);J(2*ii-1,2*ii)=-U(ii)^2*G(ii,ii)-P(ii);J(2*ii,2*ii-1)=U(ii)^2*G(ii,ii)-P(ii);J(2*ii,2*ii)=U(ii)^2*B(ii,ii)-Q(ii);elseJ(2*ii-1,2*jj-1)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U( jj)) - B(ii,jj)*cos(angle_U(ii)-angle_U(jj)) );J(2*ii-1,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)) + B(ii,jj)*sin(angle_U(ii)-angle_U(jj)) );J(2*ii,2*jj-1)=U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)) + B(ii,jj)*sin(angle_U(ii)-angle_U(jj)) );J(2*ii,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U(jj)) - B(ii,jj)*cos(angle_U(ii)-angle_U(jj)) );endelse%PV节点if ii==jjJ(2*ii-1,2*ii-1)=U(ii)^2*B(ii,ii)+Q(ii);J(2*ii-1,2*ii)=-U(ii)^2*G(ii,ii)-P(ii);elseJ(2*ii-1,2*jj-1)=-U(ii)*U(jj)*( G(ii,jj)*sin(angle_U(ii)-angle_U( jj)) - B(ii,jj)*cos(angle_U(ii)-angle_U(jj)) );J(2*ii-1,2*jj)=-U(ii)*U(jj)*( G(ii,jj)*cos(angle_U(ii)-angle_U(jj)) + B(ii,jj)*sin(angle_U(ii)-angle_U(jj)) );endendendend%在求解修正方程之前建议把delta_ef和delta_ef全部放在一个矩阵delta_PQ=[delta_P(1);delta_Q(1);delta_P(2);delta_Q(2);delta_P(3); delta_Q(3);delta_P(4)];%下面求解修正方程;注意矩阵运算时候的左除和右除的区别J=J(1:7,1:7);delta_ef=-J\delta_PQ;%下面修正各个节点的电压for i=1:1:4if i~=4U(i)=U(i)+delta_ef(2*i)*U(i);endangle_U(i)=angle_U(i)+delta_ef(2*i-1);end%到这里第一轮迭代完成end%下面显示出满足条件后的迭代的次数disp(['Iteration times : ' num2str(cont)]);%下面计算平衡节点5的功率PQfor j=1:1:5Sum_GB1=Sum_GB1 + U(j)*( G(5,j)*cos(angle_U(5)-angle_U(j)) + B(5,j)*sin(angle_U(5)-angle_U(j)) );Sum_GB2=Sum_GB2 + U(j)*( G(5,j)*sin(angle_U(5)-angle_U(j)) - B(5,j)*cos(angle_U(5)-angle_U(j)) );endP(5)=U(5)*Sum_GB1;Q(5)=U(5)*Sum_GB2;%下面将相角用角度表示for i=1:1:5angle_U(i)=angle_U(i)*180/pi;End计算计算法P-Q算法计算潮流：这个算法是由牛顿算法的极坐标形式简化而来。

3.3复杂电力网潮流计算的计算机解法3.3.1 导纳矩阵的形成1．自导纳节点i的自导纳，亦称输入导纳，在数值上等于在节点i施加单位电压，其他节点全部接地时，经节点i注入网络的电流。

主对角线元素，更具体地说，就等于与节点连接的所有支路导纳的和。

2．互导纳节点i、j间的互导纳，在数值上等于在节点i施加单位电压，其他节点全部接地时，经节点j注入网络的电流。

非对角线元素。

更具体地说，是连接节点j和节点i支路的导纳之和再加上负号而得。

3．导钠矩阵的特点：（1）因为，导纳矩阵Y是对称矩阵；（2）导纳矩阵是稀疏矩阵，每一非对角元素是节点i和j间支路导纳的负值，当i和j间没有直接相连的支路时，即为零，根据一般电力系统的特点，每一节点平均与3-5个相邻节点有直接联系，所以导纳矩阵是一高度稀疏的矩阵；（3）导纳矩阵能从系统网络接线图直观地求出。

4．节点导纳矩阵的修改（1）从原有网络引出一支路，同时增加一节点，设i为原有网络结点，j为新增节点，新增支路ij的导纳为y ij。

如图3-17（a）所示。

因新增一节点，新的节点导纳阵需增加一阶。

且新增对角元Y jj=y ij，新增非对角元Y ij=Y ji=－y ij，同时对原阵中的对角元Y ii进行修改，增加ΔY ii＝y ij。

（2）在原有网络节点i、j间增加一支路。

如图3-17（b）所示。

设在节点i增加一条支路，由于没有增加节点数，节点导纳矩阵Y阶次不变，节点的自导纳Y ii、Y jj和互导纳Y ij分别变化量为（3-57）图 3-17 网络接线的变化图（a）网络引出一支路，（b）节点间增加一支路，（c）节点间切除一支路，（d）节点间导纳改变（3）在原有网络节点i、j间切除一支路。

如图3-17（c）所示。

设在节点i切除一条支路，由于没有增加节点数，节点导纳矩阵Y阶次不变，节点的自导纳Y ii、Y jj和互导纳Y ij分别发生变化，其变化量为（3-58）（4）原有网络节点i、j间的导纳改变为。

电力系统潮流计算机算法电力系统潮流计算是电力系统分析中最基本的一项计算，其目的是确定电力系统中各母线电压的幅值和相角、各元件中的功率以及整个系统的功率损耗等。

随着计算机技术的发展，电力系统潮流计算算法也在不断更新和完善。

以下是电力系统潮流计算的一些常用算法：1. 牛顿-拉夫逊法（Newton-Raphson Method）：这是一种求解非线性方程组的方法，应用于电力系统潮流计算中。

该方法在多数情况下没有发散的危险，且收敛性较强，可以大大节约计算时间，因此得到了广泛的应用。

2. 快速迪科法（Fast Decoupled Method）：这是一种高效的电力系统潮流计算方法，将电力系统分为几个子系统进行计算，从而提高了计算速度。

3. 最小二乘法（Least Squares Method）：这是一种用于求解线性方程组的方法，通过最小化误差平方和来获得最优解。

在电力系统潮流计算中，可用于优化电压幅值和相角。

4. 遗传算法（Genetic Algorithm）：这是一种全局优化搜索算法，应用于电力系统潮流计算中，可以解决一些复杂和非线性问题。

5. 粒子群优化算法（Particle Swarm Optimization）：这是一种启发式优化算法，通过模拟鸟群觅食行为来寻找最优解。

在电力系统潮流计算中，可用于优化网络参数和运行条件。

6. 模拟退火算法（Simulated Annealing）：这是一种全局优化搜索算法，应用于电力系统潮流计算中，可以在较大范围内寻找最优解。

7. 人工神经网络（Artificial Neural Network）：这是一种模拟人脑神经网络的计算模型，可用于电力系统潮流计算。

通过训练神经网络，可以实现对电力系统中复杂非线性关系的建模和预测。

以上所述算法在电力系统潮流计算中起着重要作用，为电力系统运行、设计和优化提供了有力支持。

同时，随着计算机技术的不断发展，未来还将出现更多高效、精确的电力系统潮流计算算法。

高斯——赛德尔法潮流计算潮流计算高斯——赛德尔迭代法(Gauss一Seidel method)是求解电力系统潮流的方法。

潮流计算高斯——赛德尔迭代法又分导纳矩阵迭代法和阻抗矩阵迭代法两种。

前者是以节点导纳矩阵为基础建立的赛德尔迭代格式;后者是以节点阻扰矩阵为基础建立的赛德尔迭代格式。

高斯——赛德尔迭代法这是数学上求解线性或非线性方程组的一种常用的迭代方法。

本实验通过对电力网数学模型形成的计算机程序的编制与调试，获得形成电力网数学模型：高斯---赛德尔法的计算机程序，使数学模型能够由计算机自行形成，即根据已知的电力网的接线图及各支路参数由计算程序运行形成该电力网的节点导纳矩阵和各节点电压、功率。

通过实验教学加深学生对高斯---赛德尔法概念的理解，学会运用数学知识建立电力系统的数学模型，掌握数学模型的形成过程及其特点，熟悉各种常用应用软件，熟悉硬件设备的使用方法，加强编制调试计算机程序的能力，提高工程计算的能力，学习如何将理论知识和实际工程问题结合起来。

高斯---赛德尔法潮流计算框图开始输入数据，定义数组给定PQ节点电压初值给定PV节点电压实部（或虚部）置迭代计数b=0计算PQ节点电压实部和虚部先计算PV节点无功功率再用其计算PV节点电压实部和虚部计算平衡节点的有功和无功NY［1］系统节点的分类根据给定的控制变量和状态变量的不同分类如下①P 、Q 节点（负荷节点），给定Pi 、Qi 求Vi 、Si ，所求数量最多；②负荷节点，变电站节点（联络节点、浮游节点），给定P Gi 、QGi 的发电机节点，给定Q Gi 的无功电源节点；③PV 节点（调节节点、电压控制节点），给定P i 、Q i 求Q n 、S n ，所求数量少，可以无有功储备的发电机节点和可调节的无功电源节点；④平衡节点（松弛节点、参考节点（基准相角）、S 节点、VS 节点、缓冲节点），给定V i ，δi =0，求P n 、Q n (V s 、δs 、P s 、Q s )。

电力系统潮流计算算法设计及实现潮流计算是电力系统分析中的一种最基本的计算，它的任务是对给定的运行条件确定系统的运行状态，如各母线上的电压（幅值及相角）、网络中的功率分布以及功率损耗等。

建模是用数学的方法建立的数学模型，但它严格依赖于物理系统。

根据电力系统的实际运行条件，按给定的变量不同，一般将节点分为PQ节点，PV节点，平衡节点三种类型。

当这三个节点与潮流计算的约束条件结合起来时，便是潮流计算的数学模型。

PQ节点：有功功率P和无功功率Q是已知的，节点电压（V，δ）是待求量。

通常变电所都是这一类型的节点。

PV节点：有功功率P和电压复制V是已知的，节点的无功功率Q和电压相位δ是待求量。

一般选择有一定无功储备的发电厂和具有可调无功电源设备的变电所作为PV节点。

平衡节点：在潮流分布算出之前，网络中的功率损失是未知的，所以，网络中至少有一个节点的有功功率P不能给定，这个节点承担了系统的有功功率平衡，所以称为平衡节点。

一般选择主调频发电厂为平衡节点。

潮流计算的约束条件是：1、所有的节点电压必须满足：这一约束主要是对PQ节点而言。

2、2、所有电源节点的有功功率和无功功率必须满足：对平衡节点的P和Q以及PV节点的Q按以上条件进行检验。

3、某些节点之间电压的相位差应满足：稳定运行的一个重要条件。

功率方程的非线性雅可比矩阵的特点：●各元素是各节点电压的函数●不是对称矩阵●因为Y =0，所以H =N =J =L =0，另R =S =0，故稀疏两种常见的求解非线性方程的方法：1）高斯-赛德尔迭代法；2）牛顿-拉夫逊迭代法。

高斯-赛德尔迭代法潮流计算1、方程表示：①用高斯-赛德尔计算电力系统潮流首先要将功率方程改写成能收敛的迭代形式；②Q：设系统有n个节点，其中m个PQ节点，n-(m+1)个是PV节点，一个平衡节点,平衡节点不参加迭代；③功率方程改写成:2、求解的步骤：1)上述迭代公式假设n个节点全部为PQ节点。

2)始终等号右边采用第k次迭代结果，当j<i时，采用经(k+1)次迭代后的值，当j>i时，采用第k次迭代结果。

潮流计算的基本算法及使用方法Company number：【0089WT-8898YT-W8CCB-BUUT-202108】潮流计算的基本算法及使用方法一、 潮流计算的基本算法1.牛顿－拉夫逊法1．1 概述牛顿－拉夫逊法是目前求解非线性方程最好的一种方法。

这种方法的特点就是把对非线性方程的求解过程变成反复对相应的线性方程求解的过程，通常称为逐次线性化过程，就是牛顿－拉夫逊法的核心。

牛顿-拉夫逊法的基本原理是在解的某一邻域内的某一初始点出发，沿着该点的一阶偏导数——雅可比矩阵，朝减小方程的残差的方向前进一步，在新的点上再计算残差和雅可矩阵继续前进，重复这一过程直到残差达到收敛标准，即得到了非线性方程组的解。

因为越靠近解，偏导数的方向越准，收敛速度也越快，所以牛顿法具有二阶收敛特性。

而所谓“某一邻域”是指雅可比方向均指向解的范围，否则可能走向非线性函数的其它极值点，一般来说潮流由平电压即各母线电压(相角为0，幅值为1)启动即在此邻域内。

1．2 一般概念对于非线性代数方程组即 ()0,,,21=n i x x x f ()n i ,2,1= (1－1)在待求量x 的某一个初始计算值()0x 附件，将上式展开泰勒级数并略去二阶及以上的高阶项，得到如下的线性化的方程组()()()()()0000=∆'+x x f x f (1－2)上式称之为牛顿法的修正方程式。

由此可以求得第一次迭代的修正量()()()[]()()0100x f x f x -'-=∆ (1－3)将()0x ∆和()0x 相加，得到变量的第一次改进值()1x 。

接着再从()1x 出发，重复上述计算过程。

因此从一定的初值()0x 出发，应用牛顿法求解的迭代格式为()()()()()k k k x f x x f -=∆' (1－4)()()()k k k x x x ∆+=+1 (1－5)上两式中：()x f '是函数()x f 对于变量x 的一阶偏导数矩阵，即雅可比矩阵J ；k 为迭代次数。

潮流计算的主要方法
最近几年，随着计算机仿真技术和复杂系统全面发展，潮流计算也受到越来越多的重视。

潮流计算是研究不同电力网络的物理特性和操作规律的一项重要工作。

针对潮流计算的主要方法，总结如下：
一、基于动力学的方法
1. 碰撞模型：根据动力学方法，计算电力系统的运行稳定性。

基于动力学的碰撞模型能够快速而精确地预测两个潮流的变化情况。

2. 时变快速收敛：在碰撞模型的基础上，为快速求解电力系统潮流，提出了时变快速收敛算法。

可以更快地获得潮流解。

二、基于牛顿迭代法的方法
1.牛顿迭代潮流计算方法：根据牛顿迭代法，采用迭代算法，求解电力系统潮流运行状态。

2. 功率流计算方法：计算机基于牛顿迭代法，快速求解节点电能的功率流公式。

可以有效的缩短潮流计算的时间，提高计算效率。

三、基于模糊聚类算法的方法
1. 基于模糊聚类的潮流计算方法：采用模糊聚类算法，对潮流计算进行多维度分析，可以得出最优的潮流结果。

2. 基于模糊划分的多目标模糊控制：根据模糊聚类理论，对潮流算法进行最佳控制，以满足电力网不同优化目标。

四、基于期望最大化的方法
1、基于粒子群优化的潮流计算方法：采用粒子群优化算法，将电力网潮流计算定义为多目标最优化问题，以期望最大化来求解潮流值，提高计算效率。

2、基于遗传算法的潮流计算方法：遗传算法利用进化过程来搜索全局最优解，使用遗传变异原则来改变候选解，以期望最大化来求解潮流计算问题。

一、潮流计算的计算机方法对于复杂网络的潮流计算，一般必须借助电子计算机进行。

其计算步骤是：建立电力网络的数学模型，确定计算方法、制定框图和编制程序。

本章重点介绍前两部分，并着重阐述在电力系统潮流实际计算中常用的、基本的方法。

1，电力网络的数学模型电力网络的数学模型指的是将网络有关参数相变量及其相互关系归纳起来所组成的．可以反映网络性能的数学方程式组。

也就是对电力系统的运行状态、变量和网络参数之间相互关系的—种数学描述。

电力网络的数学模型有节点电压方程和回路电流方程等，前者在电力系统潮流计算中广泛采用。

节点电压方程又分为以节点导纳矩阵表示的节点电压方程和以节点阻抗矩阵表示的节点电压方程。

（1）节点导纳矩阵在电路理论课中。

已讲过了用节点导纳矩阵表示的节点电压方程：对于n个节点的网络其展开为：上式中，I是节点注入电流的列向量。

在电力系统计算中，节点注入电流可理解为节点电源电流与负荷电流之和，并规定电源向网络节点的注人电流为正。

那么，只有负荷节点的注入电流为负，而仅起联络作用的联络节点的注入电流为零。

U是节点电压的列向量。

网络中有接地支路时，通常以大地作参考点，节点电压就是各节点的对地电压。

并规定地节点的编号为0。

y是一个n×n阶节点导纳矩阵，其阶数n就等于网络中除参考节点外的节点数。

物理意义：节点i单位电压，其余节点接地，此时各节点向网络注入的电流就是节点i 的自导纳和其余节点的与节点i之间的互导纳。

特点：对称矩阵，稀疏矩阵，对角占优(2) 节点阻抗矩阵对导纳阵求逆，得：其中称为节点阻抗矩阵，是节点导纳矩阵的逆阵。

物理意义：节点i注入单位电流，其余节点不注入电流，此时各节点的电压就是节点i 的自阻抗和其余节点的与节点i之间的互阻抗。

特点：满阵，对称，对角占优2，功率方程、变量和节点分类（1）功率方程已知的是节点的注入功率，因此，需要重新列写方程： **==B B B B B U S I U Y其展开式为： i i i nj j ij U jQ P U Y ~1-=∑= 所以：∑=**=+nj jij i i i U Y U jQ P 1 展开写成极坐标方程的形式：)cos sin ()sin cos (11ij ij ij ij n j j i i ij ij ij ij n j j i i B G U U Q B G U U P δδδδ-=+=∑∑==所以节点的功率方程为：)cos sin ()sin cos (11ij ij ij ij n j j i di Gi i ij ij ij ij nj j i di Gi i B G U U Q Q Q B G U U P P P δδδδ---=∆+--=∆∑∑==（2） 变量分类负荷消耗的有功、无功功率取决于用户，因而是无法控制的，故称为不可控变量或扰动变量。

潮流计算方法 —计算机算法
————
需求分析：
使用MATLAB软件计算美国电气服务公司 IEEE30节点系统（图见课本P136页图6.16） 各节点的电压幅值和相角的潮流分布，并比较 各种算法的优劣。
总体思路：
1.直接使用“牛顿——拉夫逊”法，通过图像 显示每次迭代后的误差并记录程序运行时间。 2.直接使用“PQ分解法”，通过图像显示每次 迭代后的误差并记录程序运行时间。 3.首先使用直流潮流计算，计算各节点电压的 相角，并将计算结果作为 “牛顿——拉夫逊” 法的初值，通过图像显示每次迭代后的误差并 记录程序运行时间。
迭代步骤：
1.根据各节点电压幅值及相角初值，计算电压 相量。 2.计算各PQ节点的复功率。 3.计算各节点流出的电流初值。 4.根据S=VI*计算出各节点复功率，并与已知复 程序 功率相减得到误差。 5.当误差小于设定值时，迭代结束。
雅克比矩阵的形成
当 j≠i 时
Pi H ij Vi V j (Gij sin ij Bij cos ij ) j Pi N ij V j Vi V j (Gij cos ij Bij sin ij ) Vj
形成导纳矩阵



2.对角线元素： for n=1:nbus for k=1:nbr if n1(k)==n Y(n,n) = Y(n,n)+Y1(k)/(a(k)^2) + Bc(k); elseif n2(k)==n Y(n,n) = Y(n,n)+Y1(k) +Bc(k); else, end end end
Qi Lii Vi Vi V j (Gij sin ij Bij cos ij ) 2Vi 2 Bii Vi 2 Bii Qi ji i j i

高等电力系统分析（潮流计算的计算机算法）PQ分解法潮流计算（IEEE14）目录一、MATLAB源程序二、对支路参数（B1）、节点参数（B2）的说明三、带入数据，运行结果非对角元 对角元MATLAB 源程序 clearclose alln=in put(' 请输入节点数:n=');n 仁in put(' 请输入支路数：n 1=') )isb=in put(' 请输入平衡节点号 :isb=')pr=i nput(' 请输入误差精度:pr=');B1=i nput(' 请输入支路参数：B1 =');B2=i nput(' 请输入节点参数:B2=');n2=in put(' 请输入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)% --------------------------------------------% -------------- 下面是求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:n所有节点数 节点数 求 DP,DQ 对V 矩阵求逆 △ P/V △ P/V/B3 △角=-△ P/V/V/B3△ V=- △ Q/V/Bif i~=isbC(i)=O;D(i)=O;for j1=1:n C(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;endendendif 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 0 J1 4 0 0.55618 0 0 J1 10 0.05811 0.17632 0 0.0341 11 0.06701 0.17103 0 0.01282 10 0.05695 0.17388 0 0.0346 2 12 0 0.25202 0 0 J2 14 0.05403 0.22304 0 0.04923 4 0 0.11001 0 0 J3 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 07 8 0.22092 0.19988 0 0 ;[1 1 -0.4782 1 -0.0763 10 04 1 -0.2955 1 -0.096 1 -0.0357 1 -0.0618 1 -0.1359 1 -0.14910 2 0.18311 2 -0.94212 2 -0.11213 214 00.039 1 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.047 1.70;0 0.174 1.9 0;0 0 1.06 0;]7 12 0.12291 0.25581 0 08 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 alln=input('请输入节点数：n=');n1=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 节点% --------------------------------------------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 角=-A% -------------- 下面是求△ 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/Bmax2=max(abs(DQ));for i=1:m2if max2<prMq=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(0 (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 J1 4 0 0.55618 0 0 J1 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 0 J2 14 0.05403 0.22304 0 0.049210.5130 -38.2963i -6.8410 +21.5786i 0.0000 + 4.7819i0.0000 + 1.7980i0.0000 +0.0000i-6.8410 +21.5786i 9.5680 -34.8916i0.0000 + 0.0000i 0.0000 + 0.0000i0.0000 +0.0000i0.0000 + 4.7819i 0.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.3654i 5.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 4 0 0.11001 00 ; 3 13 0 0.17615 00 J 4 5 0.03181 0.08450 ; 4 9 0.12711 0.27038 0 0 ; 5 60.08205 0.19207 00 ; 6 12 0.09498 0.1989 0 0 ; 7 8 0.22092 0.19988 0 0 ; 7 12 0.12291 0.25581 00 ; 8 9 0.17093 0.34802 0 0 ; 8 12 0.06615 0.13027 0 0 ;10 11 0.04699 0.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.9420 1.01 0;12 2 -0.112 0.047 1.7 0;13 20 0.1741.9 0;14 0 0 0 1.06 0;]请输入PQ 节点个数:n2=9导纳矩阵 Y=1Columns 1 through 5；0.0000 +O.OOOOiO.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.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 + 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i-1.4240 + 3.0291i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000iColu mns 6 through 100.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.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.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.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.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 +0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i -4.9991-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.0000i 0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 3.9679i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-1.9550 + 4.0941i-1.5260 + 3.1760i-3.0989 + 6.1028i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i6.5799 -17.3407i0.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.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 - 5.6770i0.0000 + 0.0000i0.0000 + 0.0000i-1.0259 + 4.2350i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i-4.9991 +15.2631i0.0000 + 0.0000i0.0000 + 0.0000i0.0000 + 0.0000i6.0250 -19.3961i+15.2631iColu mns 11 through 14迭代次数为:42修正后各节点电压标么值为Colu mns 1 through 10（节点号从小到大排列）1.2128 1.21481.5809 1.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-0.0475 --1.5950 -各节点的功率为：Columns 1 through 5 1.4718 - 1.8800i 0.9041 - 1.1423i0.0000 - O.OOOOi -0.1314 + 0.1399i0.0646iColu mns 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.2396i3.2005 + 3.0414i1.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各条支路的首端功率为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.6677i19.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+02i21.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.3933i39.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。

(0) xn xn ) 0
在初值处按泰勒级数展开，并略去高阶项，得到：
(0) f1 (x1 x1 , x (0) 2 x 2 , (0) (0) , x (0) n x n ) f1 (x 1 , x 2 ,
, x (0) n )
f1 f x1 1 x 2 x1 0 x 2 0
2 2 20 2 1 12 2 3
23
U y (U U ) y (U U )y I 3 3 30 3 2 23 3 1 13
电力系统分析
根据KCL定律 U y (U U ) y (U U )y I 1 1 10 1 2 12 1 3 13 U y (U U ) y (U U )y I
2 2 20 2 1 12 2 3
23
U y (U U ) y (U U )y I 3 3 30 3 2 23 3 1 13
I 1 Y11 U 1 Y12 U 2 Y13 U 3 I 2 Y21 U 1 Y22 U 2 Y23 U 3 I 3 Y31 U 1 Y32 U 2 Y33 U 3
已 知
相量形式
f1 f1 f1 ... x1 0 x2 0 xn f 2 f 2 ... f 2 x1 0 x2 0 xn . . . f n f n ...f n x1 x2 xn 0 0
电力系统分析
牛顿—拉夫逊法计算非线性方程
一、N-R原理
1. 非线性方程的求解：
f(x)=0
设：x(0)为的初始近似解，x(0)为与真实解的偏差 则：x= x(0) +x(0) f(x(0) +x(0))=0

潮流计算是电力系统非常重要的分析计算，用以研究系统规划和运行中提出的各种问题。

对规划中的电力系统，通过潮流计算可以检验所提出的电力系统规划方案能否满足各种运行方式的要求；对运行中的电力系统，通过潮流计算可以预知各种负荷变化和网络结构的改变会不会危及系统的安全，系统中所有母线的电压是否在允许的范围以内，系统中各种元件(线路、变压器等)是否会出现过负荷，以及可能出现过负荷时应事先采取哪些预防措施等。

潮流计算是电力系统分析最基本的计算。

除它自身的重要作用之外，在《电力系统分析综合程序》(PSASP)中，潮流计算还是网损计算、静态安全分析、暂态稳定计算、小干扰静态稳定计算、短路计算、静态和动态等值计算的基础。

传统的潮流计算程序缺乏图形用户界面，结果显示不直接难与其他分析功能集成。

网络原始数据输入工作大量且易于出错。

本文采用MATLAB语言运行WINDOWS操作系统的潮流计算软件。

而采用MATLAB界面直观，运行稳定，计算准确。

关键词：电力系统潮流计算；牛顿—拉夫逊法潮流计算；MATLAB一、概述1.1设计目的与要求.................................................1.1.1 设计目的......................................................1.1.2 设计要求.....................................................1.2 设计题目......................................................1.3 设计内容.....................................................二、电力系统潮流计算概述.....................2.1 电力系统简介..........................................2.2 潮流计算简介..........................................2.3 潮流计算的意义及其发展..................... ..............三、潮流计算设计题目..........................3.1 潮流计算题目........................................3.2 对课题的分析及求解思路........................四、潮流计算算法及手工计算...........................4.1 变压器的∏型等值电路..............................4.2 节点电压方程..............................4.3节点导纳矩阵.............................4.4 导纳矩阵在潮流计算中的应用.......................4.5 潮流计算的手工计算..........................五、Matlab概述....................................5.1 Matlab简介............................................5.2 Matlab的应用............................................5.3 矩阵的运算...........................................5.3.1 与常数的运算.............................................5.3.2 基本数学运算..................................5.3.3 逻辑关系运算....................................5.4 Matlab中的一些命令.................................六、潮流计算流程图及源程序................................6.1 潮流计算流程图..............................6.2 潮流计算源程序图...............................6.3 运行计算结果.......................................总结参考文献1.1 设计目的与要求1.1.1设计目的1.掌握电力系统潮流计算的基本原理；2.掌握并能熟练运用一门计算机语言（MATLAB语言或C语言或C++语言）；3.采用计算机语言对潮流计算进行计算机编程。

∆x ( k +1) = − J −1[ y ( x ( 0 ) ) − y s + y∆x ( k ) ]
x
( k +1)
=x
(0)
+ ∆x
( k +1)
式中： 表示迭代次数； (0)估计而 式中：k表示迭代次数；J为按x＝x(0)估计而 得
2012-3-7
与牛拉法比较： 与牛拉法比较：
雅克比矩阵是定值， 雅克比矩阵是定值，而牛拉法雅克比矩阵需重新计算 ) 的修正量， 修正量 ∆x (k 是相对初始估计值 ∆x ( 0 ) 的修正量，而牛拉法 ∆x (k )是相对上一次迭代所得到的迭代点的 修正量 保留达到收敛所需的迭代次数比牛拉法要多， 保留达到收敛所需的迭代次数比牛拉法要多，但由于每次 迭代所需的计算量要节省很多， 迭代所需的计算量要节省很多，总的计算速度是提高很多 的 对初始值的要求更高 对于具有大R/X比值元件或有串联支路的系统，保留非线 比值元件或有串联支路的系统， 对于具有大 比值元件或有串联支路的系统 性法有更好的收敛可靠性
线路特别重载以致两节点间相角差特别大
2012-3-7
保留非线性潮流算法
为何提出？ 为何提出？ 牛拉法迭代时，采用的是逐次线性化， 牛拉法迭代时，采用的是逐次线性化，略去了泰 勒级数的高阶项， 勒级数的高阶项，出于对精确数学模型可能会提高算 法的收敛性能及计算速度的考虑而提出此算法。 法的收敛性能及计算速度的考虑而提出此算法。 算法的发展 保留非线性的快速潮流算法（极坐标形式） 保留非线性的快速潮流算法（极坐标形式） 保留非线性的快速潮流算法（直角坐标形式） 保留非线性的快速潮流算法（直角坐标形式） 采用直角坐标形式的包括二阶项的快速潮流算法
2012-3-7

电力系统三种潮流计算方法的比较 一、高斯-赛德尔迭代法：以导纳矩阵为基础,并应用高斯--塞德尔迭代的算法是在电力系统中最早得到应用的潮流计算方法,目前高斯一塞德尔法已很少使用; 将所求方程改写为 不能直接得出方程的根,给一个猜测值 得 又可取x1为猜测值,进一步得：反复猜测 则方程的根 优点：1. 原理简单,程序设计十分容易;2. 导纳矩阵是一个对称且高度稀疏的矩阵,因此占用内存非常节省;3. 就每次迭代所需的计算量而言,是各种潮流算法中最小的,并且和网络所包含的节点数成正比关系;缺点：1. 收敛速度很慢;2. 对病态条件系统,计算往往会发生收敛困难：如节点间相位角差很大的重负荷系统、包含有负电抗支路如某些三绕组变压器或线路串联电容等的系统、具有较长的辐射形线路的系统、长线路与短线路接在同一节点上,而且长短线路的长度比值又很大的系统;3. 平衡节点所在位置的不同选择,也会影响到收敛性能;二、牛顿-拉夫逊法：求解 设 ,则按牛顿二项式展开：当△x 不大,则取线性化仅取一次项则可得修正量对 得：作变量修正： ,求解修正方程 牛顿法是数学中求解非线性方程式的典型方法,有较好的收敛性;自从20世纪60年代中期采用了最佳顺序消去法以后,牛顿法在收敛性、内存要求、计算速度方面都超过了其他方法,成为直到目前仍被广泛采用的方法;优点：1. 收敛速度快,若选择到一个较好的初值,算法将具有平方收敛特性,一般迭代4—5次便可以收敛到一个非常精确的解;而且其迭代次数与所计算网络的规模基本无关;2. 具有良好的收敛可靠性,对于前面提到的对以节点导纳矩阵为基础的高斯一塞德尔法呈病态的系统,牛顿法均能可靠地收敛;3. 牛顿法所需的内存量及每次迭代所需时间均较前述的高斯一塞德尔法为多,并与程序设计技巧有密切关系;缺点：牛顿法的可靠收敛取决于有一个良好的启动初值;如果初值选择不当,算法有可能根本不收敛或收敛到一个无法运行的解点上;()0f x =10()x x ϕ=迭代 0x 21()x x ϕ=1()k k x x ϕ+=()x x ϕ=()0f x =0x x x =+∆1k k k x x x +=+∆解决方法：对于正常运行的系统,各节点电压一般均在额定值附近,偏移不会太大,并且各节点间的相位角差也不大,所以对各节点可以采用统一的电压初值也称为“平直电压”,“平直电压”法假定：︒==0100i i U θ 或 );,...,2,1(0100s i n i f e i i ≠===这样一般能得到满意的结果;但若系统因无功紧张或其它原因导致电压质量很差或有重载线路而节点间角差很大时,仍用上述初始电压就有可能出现问题;可以先用高斯一塞德尔法迭代1-2次；以此迭代结果作为牛顿法的初值,也可以先用直流法潮流求解一次以求得一个较好的角度初值,然后转入牛顿法迭代;三、P-Q 分解法：电力系统中常用的PQ 分解法派生于以极坐标表示的牛顿—拉夫逊法,其基本思想是把节点功率表示为电压向量的极坐标形式,以有功功率误差作为修正电压向量角度的依据,以无功功率误差作为修正电压幅值的依据,把有功和无功分开进行迭代其主要特点是以一个n-1阶和一个m 阶不变的、对称的系数矩阵B ,B '''代替原来的n+m-1阶变化的、不对称的系数矩阵M,以此提高计算速度,降低对计算机贮存容量的要求;P-Q 分解法在计算速度方面有显着的提高,迅速得到了推广;原理：修正方程为：⎥⎥⎦⎤⎢⎢⎣⎡∆∆⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡∆∆V V δL K N H Q P 雅克比矩阵元素的表达如下：a) 当i ≠j 时b) 当i ＝j 时对修正方程的第一个简化是：上式可分别写成以下两式在一般情况下,线路两端电压的相角差是不大的不超过100~200,因此可以认为δδij ij ij G sin ,1cos ≈B ij因此可得：B V V H ij j i ij = i,j=1,2,…,n-1B V V L ij j i ij = i,j=1,2,…,m经一系列化简得P —Q 分解法的修正方程式： ⎭⎬⎫∆''=∆∆'=∆V B Q B P δ 原P —Q 分解法的修正方程的简化形式为： ⎪⎭⎪⎬⎫∆''=∆∆'=∆V B V Q V B V PδPQ分解法的修正方程式的特点：'、替代原有的系数矩阵J,提高了计算速度, 1.以一个n-1阶和一个m-1阶系数矩阵BB''降低了对贮存容量的要求;'、替代原有的系数矩阵J,显着的提高了计算2.以迭代过程中保持不变的系数矩阵BB''速度;'、替代原有的系数矩阵J,使求逆等运算量和所需的储存容量3.以对称的系数矩阵BB''都大为减少;P-Q分解法两个主要特点：1.降阶在潮流计算的修正方程中利用了有功功率主要与节点电压相位有关,无功功率主要与节点电压幅值有关的特点,实现P-Q分解,使系数矩阵由原来的2N×2N阶降为N×N 阶,N为系统的节点数不包括缓冲节点;2.因子表固定化利用了线路两端电压相位差不大的假定,使修正方程系数矩阵元素变为常数,并且就是节点导纳的虚部;由于以上两个特点,使快速分解法每一次迭代的计算量比牛顿法大大减少;P-Q分解法只具有一次收敛性,因此要求的迭代次数比牛顿法多,但总体上快速分解法的计算速度仍比牛顿法快;快速分解法只适用于高压网的潮流计算,对中、低压网,因线路电阻与电抗的比值大,线路两端电压相位差不大的假定已不成立,用快速分解法计算,会出现不收敛问题;。

附表1：计算机计算潮流程序：%本程序的功能是用牛顿——拉夫逊法进行潮流计算% B1矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地电纳% 5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值% 4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量% 6、节点分类标号clear;n=13;%input('请输入节点数：n=');nl=13;%input('请输入支路数：nl=');isb=1;%input('请输入平衡母线节点号：isb=');pr=0.00001;%input('请输入误差精度：pr=');B1=[];%input('请输入由支路参数形成的矩阵：B1=');B2=[];%input('请输入各节点参数形成的矩阵：B2=');Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl); %－－－－－－－修改部分－－－－－－－－－－－－ym=0;SB=100;UB=220;%ym=input('您输入的参数是标么值？（若不是则输入一个不为零的数值）'); if ym~=0%SB=input('请输入功率基准值:SB=');%UB=input('请输入电压基准值:UB=');YB=SB./UB./UB;BB1=B1;BB2=B2;for i=1:nlB1(i,3)=B1(i,3)*YB;B1(i,4)=B1(i,4)./YB;enddisp('B1矩阵B1=');disp(B1)for i=1:nB2(i,1)=B2(i,1)./SB;B2(i,2)=B2(i,2)./SB;B2(i,3)=B2(i,3)./UB;B2(i,4)=B2(i,4)./UB;B2(i,5)=B2(i,5)./SB;enddisp('B2矩阵B2=');disp(B2)end% % %---------------------------------------------------for i=1:nl %支路数if B1(i,6)==0 %左节点处于低压侧p=B1(i,1);q=B1(i,2);elsep=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; %对角元K侧Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元1侧end%求导纳矩阵disp('导纳矩阵Y=');disp(Y)%----------------------------------------------------------G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部for i=1:n %给定各节点初始电压的实部和虚部e(i)=real(B2(i,3));f(i)=imag(B2(i,3));V(i)=B2(i,4); %PV节点电压给定模值endfor i=1:n %给定各节点注入功率S(i)=B2(i,1)-B2(i,2); %i节点注入功率SG-SLB(i,i)=B(i,i)+B2(i,5); %i节点无功补偿量end%=========================================================== ========P=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~=isb %非平衡节点C(i)=0;D(i)=0;for j1=1:nC(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ(Gij*ej-Bij*fj)D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ(Gij*fj+Bij*ej)endP1=C(i)*e(i)+f(i)*D(i);%节点功率P计算eiΣ(Gij*ej-Bij*fj)+fiΣ(Gij*fj+Bij*ej)Q1=C(i)*f(i)-e(i)*D(i);%节点功率Q计算fiΣ(Gij*ej-Bij*fj)-eiΣ(Gij*fj+Bij*ej)%求P',Q'V2=e(i)^2+f(i)^2; %电压模平方%========= 以下针对非PV节点来求取功率差及Jacobi矩阵元素=========if B2(i,6)~=3 %非PV节点DP=P(i)-P1; %节点有功功率差DQ=Q(i)-Q1; %节点无功功率差%=============== 以上为除平衡节点外其它节点的功率计算=================%================= 求取Jacobi矩阵===================for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/de=-dQ/dfX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/df=dQ/deX3=X2; % X2=dp/df X3=dQ/deX4=-X1; % X1=dP/de X4=dQ/dfp=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~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX3=D(i)+B(i,i)*e(i)-G(i,i)*f(i); % dQ/deX4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);% dQ/dfp=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Qm=p+1;J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△PJ(m,q)=X2;endendelse%=============== 下面是针对PV节点来求取Jacobi矩阵的元素===========DP=P(i)-P1; % PV节点有功误差DV=V(i)^2-V2; % PV节点电压误差for j1=1:nif j1~=isb&j1~=i %非平衡节点&非对角元X1=-G(i,j1)*e(i)-B(i,j1)*f(i); % dP/deX2=B(i,j1)*e(i)-G(i,j1)*f(i); % dP/dfX5=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~=isb %非平衡节点&对角元X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);% dP/deX2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);% dP/dfX5=-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:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即N=2*n+1扩展列△P、△Qfor k2=k1:N1 % 扩展列△P、△QJ(k,k2)=J(k,k2)./J(k,k); % 非对角元规格化endJ(k,k)=1; % 对角元规格化if k~=3 % 不是第三行%======================================================== ====k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行前的k4行消去for k2=k1:N1 % k3行后各行下三角元素J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算endJ(k3,k)=0;endif k==N0break;end%==========================================for 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%====上面是用线性变换方式将Jacobi矩阵化成单位矩阵=====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); %修改节点电压虚部end%------修改节点电压-----------for k=3:N0DET=abs(J(k,N));if DET>=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1endendICT2(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);sida(k)=atan(f(k)./e(k))*180./pi;E(k)=e(k)+f(k)*j;end%=============== 计算各输出量=========================== disp('各节点的实际电压标幺值E为(节点号从小到大排列)：');disp(E);EE=E*UB;disp(EE);disp('-----------------------------------------------------');disp('各节点的电压大小V为(节点号从小到大排列)：');disp(V);VV=V*UB;disp(VV);disp('-----------------------------------------------------');disp('各节点的电压相角sida为(节点号从小到大排列)：');disp(sida);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('~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~');SS=S*SB;disp(SS);disp('-----------------------------------------------------');disp('各条支路的首端功率Si为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Si(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))));Siz(i)=Si(p,q);elseSi(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))));Siz(i)=Si(p,q);enddisp(Si(p,q));SSi(p,q)=Si(p,q)*SB;ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];disp(ZF);%disp(SSi(p,q));disp('-----------------------------------------------------');enddisp('各条支路的末端功率Sj为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);if B1(i,6)==0Sj(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))));Sjy(i)=Sj(q,p);elseSj(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))));Sjy(i)=Sj(q,p);enddisp(Sj(q,p));SSj(q,p)=Sj(q,p)*SB;ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];disp(ZF);%disp(SSj(q,p));disp('-----------------------------------------------------');enddisp('各条支路的功率损耗DS为(顺序同您输入B1时一致)：');for i=1:nlp=B1(i,1);q=B1(i,2);DS(i)=Si(p,q)+Sj(q,p);disp(DS(i));DDS(i)=DS(i)*SB;ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];disp(ZF);%disp(DDS(i));disp('-----------------------------------------------------');endfigure(1);subplot(2,2,1);plot(V);xlabel('节点号');ylabel('电压标幺值');grid on;subplot(2,2,2);plot(sida);xlabel('节点号');ylabel('电压角度');grid on;subplot(2,2,3);bar(real(S));ylabel('节点注入有功');grid on;subplot(2,2,4);bar(Siz);ylabel('支路首端无功');grid on;1.冬季最大运行方式潮流计算结果:计算机运行的B1,B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.88+0.545*i 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.77+0.4772*i 1 0 0 20 0 1 0 0 20 0.88+0.545*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.2+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：节点号电压值相角值支路标号首端功率末端功率支路功率损耗1 242.0000 0 1-2 148.9391+3.327637i -143-27.45476i 5.93913-24.1271i2 231.0000 -3.0341 1-11 89.1956+37.193i -88.19803-61.4687i 0.997519-24.2757i3 227.4154 -4.4383 11-3 88.19803+61.4687i -88-54.5i 0.19803+6.9687i4 226.3304 -4.9004 1-12 77.8364+36.4186i -77.2404-55.7062i 0.596039-19.2876i5 229.8318 -3.9547 12-4 77.2404+55.7062i -77-47.72i 0.24036+7.9862i6 217.7226 -8.2720 1-5 63.5438+12.2204i -61.8773-32.0321i 1.66656-19.8116i7 234.9326 -3.1047 5-6 77.2597+56.3501i -77-47.72i 0.25974+8.6301i8 226.0991 -6.2429 5-7 15.3825-24.3181i 15.5354+0.274108i 0.152961-24.044i9 231.0000 0.5851 7-8 88.20085+61.55009i -88-54.5i 0.20085+7.0501i10 231.0000 3.4950 1-7 40.806-6.05737i -40.0647-22.0555i 0.741264-28.1128i11 236.1931 -1.3348 7-13 -63.6716-39.7687i 64.0965+17.5832i 0.424934-22.1855i12 237.9346 -0.8892 13-9 -64.0965-17.5832i 64.2+21.0187i 0.10348+3.4355i13 238.1868 -2.2028 1-10 79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果图形：2.冬季最小运行方式潮流计算结果：计算机运行的B1B2矩阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.01+0.033*i 0.204*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.43+0.886*i 1.05 1.05 0 30 0.616+0.3817*i 1 0 0 20 0.539+0.3817*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20 0 1 0 0 20 0.539+0.334*i 1 0 0 20.642+0.3817*i 0 1.05 1.05 0 31+0.75*i 0.1+0.06197*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]电压调整措施：变电所1、4变压器变比：+2.5% 水电厂变压器变比：-2.5%5 234.2241 -2.0051 (12--4) 54.0234+42.2706i -53.9-38.17i 0.12342+4.1006i6 226.2159 -4.8577 (1--5) 34.1669+4.22856 -33.6427-28.276i 0.524165-24.0474i7 235.1128 -0.4737 (5--6) 54.0179+37.3169i -53.9-33.4i 0.11789+3.9169i8 223.9941 -2.4603 (5--7) -20.3752-9.04094i 20.5371-15.4558i 0.161947-24.4968i9 231.0000 3.4429 (1--7) 11.0013+1.45186i -10.8258-31.4513i 0.175442-29.9995i10 231.0000 3.9321 (7--8) 53.9768+36.0957i -53.9-33.4i 0.076797+2.6957i11 238.4187 -0.9561 (7--13) -63.6881+10.8115i 64.0874-32.7763i 0.39932-21.9648i12 239.1742 -0.6185 (13--9) -64.0874+32.7763i 64.2-29.0386i 0.11258+3.7377i13 234.9468 0.6942 (1--10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机运行结果的图形：3.夏季最大运行方式计算机计算结果：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.776+0.481*i 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.543+0.3367*i 1 0 0 20 0 1 0 0 20 0.776+0.481*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]计算机运行结果如下表：10 231.0000 3.4950 (7,8) 77.7585+53.6634i -77.6-48.1i0.1585+5.5634i11 237.0938 -1.1858 (7,13) -113.4984+9.476406i 114.6926-28.38885i 1.19422-18.9124i12 239.1742 -0.6185 (13,9) -114.6926+28.38885 115-18.18369i0.307384+10.2052i13 233.3912 3.2303 (1,10) -79.8613+4.23546i 80+0.641048i 0.13875+4.8765i 计算机计算结果如图：4.夏季最小运行方式：计算机运行B1B2阵如下：B1=[ 1 2 0.0318+0.0454*i 0.282*i 1 01 11 0.0114+0.0374*i 0.2332*i 1 011 3 0.001975+0.0695*i 0 1.025:1.1 11 12 0.0087+0.029*i 0.1788*i 1 012 4 0.0031+0.103*i 0 1:1.05 11 5 0.043+0.142*i 0.22*i 1 05 6 0.0031+0.103*i 0 1:1.05 15 7 0.043+0.142*i 0.22*i 1 01 7 0.051+0.168*i 0.26*i 1 07 8 0.00198+0.0695*i 0 1:1.1 17 13 0.02+0.066*i 0.102*i 1 013 9 0.0025+0.083*i 0 0.9956 11 10 0.00239+0.084*i 0 1.048 1]B2=[0 0 1.1 0 0 10 1.26+0.7814*i 1.05 1.05 0 30 0.543+0.336*i 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.4753+0.2947*i 1 0 0 20 0 1 0 0 20 0.543+0.336*i 1 0 0 21.35+0.6538*i 0.2+0.124*i 1.05 1.05 0 31+0.75*i 0.25+0.1549*i 1.05 1.05 0 30 0 1 0 0 20 0 1 0 0 20 0 1 0 0 2]10 231.0000 3.9321 (7,8) 54.3735+36.2509i -54.3-33.67i 0.073528+2.5809i11 239.0099 -0.8518 (7,13) -113.4428+24.39707i 114.6754-43.79301i 1.23255-19.3959i12 239.8726 -0.5689 (13,9) -114.6754+43.79301i 115-33.01613i 0.324605+10.7769i13 235.9565 4.4510 (1,10) -89.8244+5.1673i 90+1.0049i 0.17561+6.1722i 计算机计算结果如图：5.夏季故障运行状态:调压及无功补偿措施如下:变电所3的变压器变比为-2.5%,无功补偿容量为20Mvar。

在节点 k 单独注入电流，所 有其它节点的注入电流都等 于 0 时，在节点 k 产生的电 压同注入电流之比
从节点 k 向整个网络看进去 的对地总阻抗
11
一、节点电压方程
2、节点阻抗矩阵
Z 矩阵元素的物理意义互阻抗
if k i
在节点 k 单独注入电流，所
Z ik

U i Ik
有其它节点的注入电流都等 于 0 时，在节点 i 产生的电 I j 0, jk 压同注入电流之比
求PGs、QGs。
23
3－2 潮流方程及其迭代解法
二、高斯－赛德尔迭代法（既可解线性，
设有方程组
也可解非线性方程）
a11x1 a12 x2 a13x3 b1 a21x1 a22 x2 a23x3 b2 a31x1 a32 x2 a33x3 b3
24
3－2 潮流方程及其迭代解法
电力网
的投切或变比的调整等）
Y Y (0) Y
Yij Yij(0) Yij
14
三、节点导纳矩阵
Y 矩阵的修改
Y11 Y12 Y1i Y1 j Y1n
Y21
Y22

Y2i
Y2 j

Y2n


电力网
Y
(0)

Yi1
(2) PV节点： PLi、 PGi ，从而Pi给定； ULi 、UGi给定。 即相应的Pi、Ui给定，待求QGi、δi。如有一定无功储备 电源变电所母线（很少，甚至没有）。
(3) 平衡节点(Vδ节点）： 一般只有一个。设s节点
为平衡节点，则： Us 、 δs 给定， Us ＝1.0， δs ＝0。待