基于matlab实现的14节点系统牛顿拉夫逊法潮流计算程序.rar

clear%清除内存变量 format long jc=0; %--------------------------节点原始数据模块--------------------------------% JieDianShuJu=[1 1.0600 0.0000 0.0000 0.0000 0.0000 0.0000; 2 1.0450 0.0000 0.4000 0.4000 0.2170 0.1270; 3 1.0100 0.0000 0.0000 0.2000 0.9420 0.1900; 4 1.0000 0.0000 0.0000 0.0000 0.4780 -0.0390; 5 1.0000 0.0000 0.0000 0.0000 0.0760 0.0160; 6 1.0700 0.0000 0.0000 0.3000 0.1120 0.0750; 7 1.0000 0.0000 0.0000 0.0000 0.0000 0.0000; 8 1.0900 0.0000 0.0000 0.3000 0.0000 0.0000; 9 1.0000 0.0000 0.0000 0.0000 0.2950 0.1660; 10 1.0000 0.0000 0.0000 0.0000 0.0900 0.0580; 11 1.0000 0.0000 0.0000 0.0000 0.0350 0.0180; 12 1.0000 0.0000 0.0000 0.0000 0.0610 0.0160; 13 1.0000 0.0000 0.0000 0.0000 0.1350 0.0580; 14 1.0000 0.0000 0.0000 0.0000 0.1490 0.0500];%输入节点数据 JieDianZongShu=size(JieDianShuJu,1);%节点总数 PVJieDianShu=4;%PV节点数 PQJieDianShu=JieDianZongShu-PVJieDianShu-1;%PQ节点数 PVI=[0 1 1 0 0 1 0 1 0 0 0 0 0 0];%标示某节点是否为PV节点的矩阵 PQI=[0 0 0 1 1 0 1 0 1 1 1 1 1 1];%标示某节点是否为PQ节点的矩阵 PVR=[1 2 3 6 8];%松弛节点与PV节点编号 PV=[2 3 6 8];%PV节点编号 PQ=[4 5 7 9 10 11 12 13 14];%PQ节点编号 SC=[1 2 5 7];%N,L,M中需删除的行与列 U=JieDianShuJu(:,2);%电压有效值 Theta=JieDianShuJu(:,3);%电压相位 PG=JieDianShuJu(:,4);%发电机有功 QG=JieDianShuJu(:,5);%发电机无功 PL=JieDianShuJu(:,6);%负载有功 QL=JieDianShuJu(:,7);%负载无功 DetP=zeros(JieDianZongShu,1);%节点有功功率误差向量 DetQ=zeros(JieDianZongShu,1);%节点无功功率误差向量 sump=0;%临时变量 sumq=0;%临时变量 DetTheta=0;%临时变量 %==========================节点原始数据模块================================% % %--------------------------定义雅可比阵和修正向量---------------------------% XiuZhengXiangLiang=zeros(JieDianZongShu+PQJieDianShu-1,1);%定义修正向量 ThetaXiuZheng=XiuZhengXiangLiang(1:JieDianZongShu-1);%定义幅角修正向量 UXiuZheng=XiuZhengXiangLiang(JieDianZongShu:end);%定义电压修正向量 J=zeros(JieDianZongShu*2-2);%定义雅可比矩阵 H=zeros(JieDianZongShu-1);%初始化H阵 N=zeros(JieDianZongShu-1);%初始化N阵 L=zeros(JieDianZongShu-1);%初始化L阵 M=zeros(JieDianZongShu-1);%初始化M阵 H1=H;N1=N;M1=M;L1=L; %==========================定义雅可比阵和修正向量===========================% % %--------------------------支路原始数据模块--------------------------------% ZhiLuShuJu=[-1 4 7 0.00000 0.20912 0.00000 0.978; -2 4 9 0.00000 0.55678 0.00000 0.969; -3 5 6 0.00000 0.25202 0.00000 0.932; 0 0 9 0.00000 0.00000 0.01900 0.000 1 1 2 0.01938 0.05917 0.02640 0.000; 2 2 3 0.04699 0.19797 0.02190 0.000; 3 2 4 0.05811 0.17632 0.01870 0.000; 4 1 5 0.05403 0.22304 0.02460 0.000; 5 2 5 0.05695 0.17388 0.01700 0.000; 6 3 4 0.06701 0.17103 0.01730 0.000; 7 4 5 0.01335 0.04211 0.00640 0.000; 8 7 8 0.00000 0.17615 0.00000 0.000; 9 7 9 0.00000 0.11001 0.00000 0.000; 10 9 10 0.03181 0.08450 0.00000 0.000; 11 6 11 0.09498 0.19890 0.00000 0.000; 12 6 12 0.12291 0.15581 0.00000 0.000; 13 6 13 0.06615 0.13027 0.00000 0.000; 14 9 14 0.12711 0.27038 0.00000 0.000; 15 10 11 0.08205 0.19207 0.00000 0.000; 16 12 13 0.22092 0.19988 0.00000 0.000; 17 13 14 0.17093 0.34802 0.00000 0.000];%输入支路原始数据 ZhiLuShu=size(ZhiLuShuJu,1);%支路总数 %==========================支路原始数据模块================================% % %--------------------------等值导纳参数模块--------------------------------% DengZhiDaoNaCanShu=zeros(ZhiLuShu,3); for k=1:ZhiLuShu if ZhiLuShuJu(k,1)<0%支路标号为负，即变压器支路 DengZhiDaoNaCanShu(k,1)=1/(ZhiLuShuJu(k,7)*j*ZhiLuShuJu(k,5)); DengZhiDaoNaCanShu(k,2)=DengZhiDaoNaCanShu(k,1)*(ZhiLuShuJu(k,7)-1); DengZhiDaoNaCanShu(k,3)=-DengZhiDaoNaCanShu(k,2)/ZhiLuShuJu(k,7); elseif ZhiLuShuJu(k,1)>0%支路标号为正，即输电线路支路 DengZhiDaoNaCanShu(k,1)=1/(ZhiLuShuJu(k,4)+j*ZhiLuShuJu(k,5)); DengZhiDaoNaCanShu(k,2)=j*ZhiLuShuJu(k,6); DengZhiDaoNaCanShu(k,3)=j*ZhiLuShuJu(k,6); else%并联支路 DengZhiDaoNaCanShu(k,3)=j*ZhiLuShuJu(k,6); end end %==========================等值导纳参数模块================================% % %--------------------------节点自导纳计算模块------------------------------% JieDianZiDaoNa=zeros(JieDianZongShu,1);%定义节点自导纳矩阵 for k=1:JieDianZongShu for l=1:ZhiLuShu if ZhiLuShuJu(l,2)==k%k是支路首节点 JieDianZiDaoNa(k)=JieDianZiDaoNa(k)+DengZhiDaoNaCanShu(l,1)+DengZhiDaoNaCanShu(l,2); elseif ZhiLuShuJu(l,3)==k%k是支路末节点 JieDianZiDaoNa(k)=JieDianZiDaoNa(k)+DengZhiDaoNaCanShu(l,1)+DengZhiDaoNaCanShu(l,3); end end end %==========================节点自导纳计算模块===============================% % %--------------------------节点互导纳计算模块------------------------------% JieDianHuDaoNa=zeros(JieDianZongShu);%初始化节点互导纳矩阵 for k=1:ZhiLuShu if ZhiLuShuJu(k,1)==0; continue; end JieDianHuDaoNa(ZhiLuShuJu(k,2),ZhiLuShuJu(k,3))=-DengZhiDaoNaCanShu(k,1); JieDianHuDaoNa(ZhiLuShuJu(k,3),ZhiLuShuJu(k,2))=-DengZhiDaoNaCanShu(k,1); end %==========================节点互导纳计算模块===============================% % %--------------------------合成节点导纳矩阵--------------------------------% Y=JieDianHuDaoNa; for k=1:JieDianZongShu Y(k,k)=JieDianZiDaoNa(k); end G=real(Y);%生成节点电导矩阵 B=imag(Y);%生成节点电纳矩阵 %==========================合成节点导纳矩阵=================================% % %--------------------------生成节点邻接矩阵---------------------------------% LinJieZhen=zeros(JieDianZongShu); for k=1:ZhiLuShu if ZhiLuShuJu(k,1)==0; continue; end LinJieZhen(ZhiLuShuJu(k,2),ZhiLuShuJu(k,3))=1; LinJieZhen(ZhiLuShuJu(k,3),ZhiLuShuJu(k,2))=1; end for k=1:JieDianZongShu LinJieZhen(k,k)=1; end %==========================生成节点邻接矩阵=================================% % %------------------------------------牛顿迭代循环---------------------------------------% for DieDaiCiShu=1:10 %------------------------------功率误差计算模块-------------------------% for k=2:JieDianZongShu sump=0; sumq=0; for l=1:JieDianZongShu DetTheta=Theta(k)-Theta(l); sump=sump+LinJieZhen(k,l)*U(l)*(G(k,l)*cos(DetTheta)+B(k,l)*sin(DetTheta)); sumq=sumq+LinJieZhen(k,l)*U(l)*(G(k,l)*sin(DetTheta)-B(k,l)*cos(DetTheta)); end sump=sump*U(k); sumq=sumq*U(k); DetP(k)=PG(k)-PL(k)-sump; DetQ(k)=QG(k)-QL(k)-sumq; end YouGongWuCha=DetP; WuGongWuCha=DetQ; YouGongWuCha(1)=[]; WuGongWuCha(PVR,:)=[]; %==============================功率误差计算模块=========================% % %------------------------------迭代循环终止条件—————————————% if max(abs(YouGongWuCha))<=1e-20&max(abs(WuGongWuCha))<=1e-20 jc=1; break; end %=============================迭代循环终止条件==========================% % %-----------------------------雅可比阵计算模块-------------------------% for k=2:JieDianZongShu for l=2:JieDianZongShu DetTheta=Theta(k)-Theta(l); H(k-1,l-1)=-U(k)*U(l)*(G(k,l)*sin(DetTheta)-B(k,l)*cos(DetTheta)); N(k-1,l-1)=-U(k)*U(l)*(G(k,l)*cos(DetTheta)+B(k,l)*sin(DetTheta)); end end L=H; M=-N; for k=2:JieDianZongShu H(k-1,k-1)=U(k)*U(k)*B(k,k)+(QG(k)-QL(k)); N(k-1,k-1)=-U(k)*U(k)*G(k,k)-(PG(k)-PL(k)); M(k-1,k-1)=U(k)*U(k)*G(k,k)-(PG(k)-PL(k)); L(k-1,k-1)=U(k)*U(k)*B(k,k