牛顿拉夫逊潮流计算Word下载.docx

1、 （5）是阶方阵，其元素为；阶矩阵，其元素为。在这里把节点不平衡功率对节点电压幅值的偏导数都乘以该节点电压，相应地把节点电压的修正量都除以该节点的电压幅值，这样，雅可比矩阵元素的表达式就具有比较整齐的形式。对式（2）和式（3）求偏导数，可以得到雅可比矩阵元素的表达式如下当时 （6） （7）二、 程序流程框图三、 程序清单% B1为支路参数矩阵：1、支路首端号；2、末端号；3、支路阻抗；4、支路对地导纳；5、支路的变比；6、支路首端处于K侧为1，1侧为0% B2为节点参数矩阵：1、该节点发电机功率；2、该节点负荷功率；3、节点电压初始值4、PV节点电压V的给定值；5、节点所接的无功补偿设备的容量

2、 6、节点分类标号：1为平衡节点；2为PQ节点；3为PV节点；clear;n=4;%input（请输入节点数：n=）;nl=4;请输入支路数：nl=isb=1;请输入平衡母线节点号：isb=pr=0.00001;请输入误差精度：pr=B1=1 2 0.1880i -0.6818i 1 1;1 3 0.1302+0.2479i 0.0129i 1 0;1 4 0.1736+0.3306i 0.0172i 1 0;3 4 0.2603+0.4959i 0.0259i 1 0;B2=0 0 1 0 0 2;0 0.5+0.3i 1 0 0 2;0.2 0 1.05 1.05 0 3;0 0.15+0

3、.1i 1.05 1.05 0 1;Y=zeros（n）;e=zeros（1,n）;f=zeros（1,n）;V=zeros（1,n）;sida=zeros（1,n）;S1=zeros（nl）;% % %-for i=1:nl %支路数if B1（i,6）=0 %左节点处于1侧p=B1（i,1）;q=B1（i,2）;else %左节点处于K侧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; %对角

4、元K侧Y（p,p）=Y（p,p）+1./B1（i,3）+B1（i,4）./2; %对角元1侧%求导纳矩阵disp（导纳矩阵 Y=disp（Y）%-G=real（Y）;B=imag（Y）; %分解出导纳阵的实部和虚部n %给定各节点初始电压的实部和虚部e（i）=real（B2（i,3）;f（i）=imag（B2（i,3）;V（i）=B2（i,4）; %PV节点电压给定模值n %给定各节点注入功率S（i）=B2（i,1）-B2（i,2）; %i节点注入功率SG-SLB（i,i）=B（i,i）+B2（i,5）; %i节点无功补偿量%=P=real（S）;Q=imag（S）; %分解出各节点注入的有功

5、和无功功率ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数ICT1、a；不满足收敛要求的节点数IT2while IT2=0 % N0=2*n 雅可比矩阵的阶数；N=N0+1扩展列IT2=0;a=a+1;nif i=isb %非平衡节点C（i）=0;D（i）=0;for j1=1:C（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）P1=C（i）*e（i）+f（i）*D（i）;%节点功率P计算ei（Gij*e

6、j-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）%求i节点有功和无功功率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矩阵 =if j1=isb&j1=i %非平衡节点&非对角元X1=-G（i,j1）*e（i）-B（i,j1）*f（i）; % dP

7、/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; % X3=dQ/de J（p,N）=DQ节点无功功率差J（m,q）=X1;J（m,N）=DP;q=q+1; % X1=dP/de J（m,N）=DP节点有功功率差J（p,q）=X4;J（m,q）=X2; % X4=dQ/df X2=dp/dfelseif j1=i&j1=isb %非平衡节点&对角元

8、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/df%扩展列Q%扩展列Pelse%= 下面是针对PV节点来求取Jacobi矩阵的元素 = % PV节点有功误差DV=V（i）2-V2; % PV节点电压误差 % dP/de % dP/dfX5=0;X6=0;J（p,q）=X5;J（p,N）=DV; % PV节点电压误差J（p,q）=

9、X6; % PV节点有功误差j1=isb %非平衡节点&X5=-2*e（i）;X6=-2*f（i）;%以上为求雅可比矩阵的各个元素及扩展列的功率差或电压差 for k=3:N0 % N0=2*n （从第三行开始，第一、二行是平衡节点）k1=k+1;N1=N; % N=N0+1 即 N=2*n+1扩展列P、Q 或 Ufor k2=k1:N1 % 从k+1列的Jacobi元素到扩展列的P、Q 或 UJ（k,k2）=J（k,k2）./J（k,k）;% 用K行K列对角元素去除K行K列后的非对角元素进行规格化J（k,k）=1; % 对角元规格化K行K列对角元素赋1%= 回代运算=if k=3 % 不是第

10、三行 k 3k4=k-1;for k3=3:k4 % 用k3行从第三行开始到当前行的前一行k4行消去N1 % k3行后各行上三角元素J（k3,k2）=J（k3,k2）-J（k3,k）*J（k,k2）;%消去运算（当前行k列元素消为0）end %用当前行K2列元素减去当前行k列元素乘以第k行K2列元素J（k3,k）=0; %当前行第k列元素已消为0if k=N0 %若已到最后一行break;%= 前代运算=for k3=k1:N0 % 从k+1行到2*n最后一行N1 % 从k+1列到扩展列消去k+1行后各行下三角元素%消去运算else %是第三行k=3%= 第三行k=3的前代运算=N0 %从第四

11、行到2n行（最后一行）N1 %从第四列到2n+1列（即扩展列）%消去运算（当前行3列元素消为0）end %用当前行K2列元素减去当前行3列元素乘以第三行K2列元素 %当前行第3列元素已消为0%=上面是用线性变换方式高斯消去法将Jacobi矩阵化成单位矩阵=2:N0-1L=（k+1）./2;e（L）=e（L）-J（k,N）; %修改节点电压实部f（L）=f（L）-J（k1,N）; %修改节点电压虚部%-修改节点电压-N0DET=abs（J（k,N）;if DET=pr %电压偏差量是否满足要求IT2=IT2+1; %不满足要求的节点数加1ICT2（a）=IT2; %不满足要求的节点数ICT1=I

12、CT1+1; %迭代次数%用高斯消去法解w=-J*V迭代次数：disp（ICT1）;没有达到精度要求的个数：disp（ICT2）;for k=1:V（k）=sqrt（e（k）2+f（k）2）; %计算各节点电压的模值sida（k）=atan（f（k）./e（k）*180./pi; %计算各节点电压的角度E（k）=e（k）+f（k）*j; %将各节点电压用复数表示%= 计算各输出量 =各节点的实际电压标幺值E为（节点号从小到大排列）：disp（E）; %显示各节点的实际电压标幺值E用复数表示-各节点的电压大小V为（节点号从小到大排列）：disp（V）; %显示各节点的电压大小V的模值各节点的电压

13、相角sida为（节点号从小到大排列）：disp（sida）; %显示各节点的电压相角for p=1:C（p）=0;for q=1:C（p）=C（p）+conj（Y（p,q）*conj（E（q）; %计算各节点的注入电流的共轭值S（p）=E（p）*C（p）; %计算各节点的功率 S = 电压 X 注入电流的共轭值各节点的功率S为（节点号从小到大排列）：disp（S）; %显示各节点的注入功率各条支路的首端功率Si为（顺序同您输入B1时一致）：nlif B1（i,6）=0Si（p,q）=E（p）*（conj（E（p）*conj（B1（i,4）./2）+（conj（E（p）*B1（i,5）.-con

14、j（E（q）*conj（1./（B1（i,3）*B1（i,5）;Siz（i）=Si（p,q）;Si（p,q）=E（p）*（conj（E（p）*conj（B1（i,4）./2）+（conj（E（p）./B1（i,5）.disp（Si（p,q）;SSi（p,q）=Si（p,q）;ZF=S（,num2str（p）,num2str（q）,）=,num2str（SSi（p,q）;disp（ZF）;各条支路的末端功率Sj为（顺序同您输入B1时一致）：Sj（q,p）=E（q）*（conj（E（q）*conj（B1（i,4）./2）+（conj（E（q）./B1（i,5）.-conj（E（p）*conj（1.

15、/（B1（i,3）*B1（i,5）;Sjy（i）=Sj（q,p）;Sj（q,p）=E（q）*（conj（E（q）*conj（B1（i,4）./2）+（conj（E（q）*B1（i,5）.disp（Sj（q,p）;SSj（q,p）=Sj（q,p）;,num2str（SSj（q,p）;各条支路的功率损耗DS为（顺序同您输入B1时一致）：DS（i）=Si（p,q）+Sj（q,p）;disp（DS（i）;DDS（i）=DS（i）;DS（,num2str（DDS（i）;figure（1）;subplot（1,2,1）;plot（V）;xlabel（节点号ylabel（电压标幺值grid on;subpl

16、ot（1,2,2）;plot（sida）;电压角度figure（2）;subplot（2,2,1）;bar（P）;节点注入有功subplot（2,2,2）;bar（Q）;节点注入无功subplot（2,2,3）;P1=real（Siz）;Q1=imag（Siz）;bar（P1）;支路号支路首端注入有功subplot（2,2,4）;bar（Q1）;支路首端注入无功运行结果：导纳矩阵 Y= 2.9056 -11.1778i 0 + 5.3191i -1.6606 + 3.1617i -1.2450 + 2.3710i 0 + 5.3191i 0 - 5.6600i 0 0 -1.6606 + 3.

17、1617i 0 2.4904 - 4.7233i -0.8298 + 1.5809i -1.2450 + 2.3710i 0 -0.8298 + 1.5809i 2.0749 - 3.9304i 4 5 5 2 0 1.0000 0.8686 - 0.0940i 1.0500 + 0.0058i 0.9879 - 0.0189i- 1.0000 0.8736 1.0500 0.9881 0 -6.1767 0.3150 -1.0945 0.4585 + 0.8816i -0.5000 - 0.3000i 0.2000 + 0.2171i -0.1500 - 0.1000i 0.5000 + 1

18、.0400iS（1,2）=0.5+1.04i -0.1013 - 0.1549iS（1,3）=-0.10126-0.1549i 0.0598 - 0.0035iS（1,4）=0.059786-0.0034655i 0.0945 + 0.0678iS（3,4）=0.09454+0.067788i -0.5000 - 0.3000iS（2,1）=-0.5-0.3i 0.1055 + 0.1493iS（3,1）=0.10546+0.14934i -0.0592 - 0.0123iS（4,1）=-0.059161-0.012341i -0.0908 - 0.0877iS（4,3）=-0.090839-0.087659i -0.0000 + 0.7400iDS（1,2）=-5.5511e-017+0.73997i 0.0042 - 0.0056iDS（1,3）=0.0042042-0.0055564i 0.0006 - 0.0158iDS（1,4）=0.00062508-0.015806i-