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
