clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%B1来表示支路参数,以下是关于B1矩阵的说明
%B1的第一列表示支路的首字母
%B1的第二列表示支路的末字母
%B1的第三列表示支路的电阻参数
%B1的第四列表示支路的电抗参数
%B1的第五列表示支路的电纳参数;注明:表示电纳参数的时候只是用了其数值
%B1的第六列表示支路的变压器变比;注明:变压器阻抗总是折算到变压器的低压侧,支路的首端与变压器的低压侧阻抗相连;
%B1的第七列表示支路是否含有变压器,如果有则为1,否则为0
B1=[1,2,1.099,56.911,0,1.05,1;
2,3,4,32.5,0.0000312,1,0;
2,7,5,40.625,0.000039,1,0;
3,7,3,24.375,0.0000234,1,0;
6,3,1.22,76.963,0,0.956,1;
6,4,0.749,-5.273,0,1.05,1;
6,5,0.791,50.845,0,1,1];
%B2来表示节点参数,以下是关于B2矩阵的说明
%第一列表示该节点电源注入的功率
%第二列表示该节点负荷所需要的功率
%第三列表示该节点电压的初始值
%第四列表示PV节点电压的规定值
%第五列表示该节点的无功补偿容量
%第六列表示该节点类型,1表示平衡节点,2表示PQ节点,3表示PV节点;
B2=[1.5+0.9i,0.0019+0.00736i,1,0,0,2;
0,1+0.936i,1,0,0,2;
0,-0.0499i,1,0,0,2;
0,1+0.8i,1,0,0,2;
0,0.2+0.15i,1,0,0,2;
0,0.00075+0.00079i,1,0,0,2;
0,0,1.05,0,0,1];
%n为网络的节点数
n=7;
%b为网络支路的条数
b=7;
%per为误差的精度
per=0.00001;
%ns为不满足误差精度节点的个数,迭代次数用k表示
%PE,PF,QE,QF是PQ节点在雅格比矩阵中的元素
%PE,PF,VE,VF是PV节点在雅格比矩阵中的元素
PE=zeros(n-1,n-1);PF=zeros(n-1,n-1);QE=zeros(n-1,n-1);QF=zeros(n-1,n-1);
VE=zeros(n-1,n-1);VF=zeros(n-1,n-1);
%定义雅格比矩阵
J=zeros(2*(n-1),2*(n-1));
%雅格比矩阵形成的增广矩阵
JW=zeros(2*(n-1),2*(n-1)+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第一步形成节点导纳矩阵Y
Y=zeros(n,n);
for i=1:n
if B1(i,n)==0 %the road don't have transformer
p=B1(i,1);
q=B1(i,2);
Y(p,q)=-(Y(p,q)+1/(B1(i,3)+sqrt(-1)*B1(i,4)));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+B1(i,5)*sqrt(-1)+1/(B1(i,3)+sqrt(-1)*B1(i,4));
Y(q,q)=Y(q,q)+B1(i,5)*sqrt(-1)+1/(B1(i,3)+sqrt(-1)*B1(i,4));
else %the road have transformer
p=B1(i,1);
q=B1(i,2);
Y(p,q)=-(Y(p,q)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2)));
Y(q,p)=Y(p,q);
Y(p,p)=Y(p,p)+B1(i,5)*sqrt(-1)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2));
Y(q,q)=Y(q,q)+B1(i,5)*sqrt(-1)+1/((B1(i,3)+sqrt(-1)*B1(i,4))*(B1(i,6)^2));
end
end
disp('得出网络的节点导纳矩阵Y');
disp(Y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第二步做一些准备工作,为求出功率的偏移量做准备,这些做准备的偏移量都是以多行一列的形式表达
%1)分解出节点导纳矩阵的实部G和虚部B
G=real(Y);
B=imag(Y);
%求出各节点输入输出功率的初值;注明:流入网络的功率为正
%各节点注入的有功为P,注入的无功为Q;注明此时得出的无功功率是指数值大小,不含有i;
P=zeros(n-1,1);
Q=zeros(n-1,1);
for i=1:n-1
P(i,1)=real(B2(i,1))-real(B2(i,2));
Q(i,1)=imag(B2(i,1))-imag(B2(i,2))+B(i,5);
end
disp(P);
disp(Q);
%2)先定义原始的PV电压值,再分解出电压值的实部e和虚部f,这里要特别清楚复数实部和虚部的定义
for i=1:n
V(i,1)=B2(i,n-3); %PV节点电压的初值
end
disp(V);
E=zeros(n,1);
F=zeros(n,1);
for i=1:n
if B2(i,n-1)==3
e(i,1)=real(B2(i,n-3));
f(i,1)=imag(B2(i,n-3));
else
e(i,1)=real(B2(i,n-4));
f(i,1)=imag(B2(i,n-4));
end
end
disp(e);
disp(f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第三步奏开始判断且求出第一次迭代量
%定义迭代次数k
k=0;
%定义含有不符合误差精度的节点个数ns,先假设ns=1;
ns=1;
%由于下面的程序必须涉及到功率偏移量或者电压的偏移量
%定义有功偏移量XP,无功偏移量XQ,电压的偏移量XV,为了万能程序的形成,此处先不定义明确
XP=zeros(n-1,1);
XQ=zeros(n-1,1);
XV=zeros(n-1,1);
%%%%%%%%%%%%%!!!!说明由于要反复迭代计算,里面的元素不能够赋值为0!!!!!!!!
while ns~=0;
ns=0; %先赋值为0,为后面的判断做铺垫
k=0;
for i=1:n-1
for j=1:n
XP(i,1)=XP(i,1)+e(i,1)*(G(i,j)*e(j,1)-B(i,j)*f(j,1))+f(i,1)*(G(i,j)*f(j,1)+B(i,j)*e(j,1));
XQ(i,1)=XQ(i,1)+f(i,1)*(G(i,j)*e(j,1)-B(i,j)*f(j,1))-e(i,1)*(G(i,j)*f(j,1)+B(i,j)*e(j,1));
XV(i,1)=XV(i,1)+e(i,1)^2+f(i,1)^2;
end
if B(i,n-1)==3
XV(i,1)=V(i,1)^2-XV(i,1);
else
XP(i,1)=P(i,1)- XP(i,1);
XQ(i,1)=Q(i,1)- XQ(i,1);
end
end
disp('功率的偏移量')
disp(XP);
disp(XQ);
disp(XV); %XV的结果都市7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%做几点说明:本来在此处判断三个偏移量是否满足收敛要求,但是为了程序的连续性,先求雅格比矩阵求出电压的偏移量再做判断,并且程序在此处一般不能满足收敛要求%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第四部求雅格比矩阵
for i=1:(n-1)
for j=1:(n-1)
if i~=j
PE(i,j)=-(G(i,j)*e(i,1)+B(i,j)*f(i,1));
QF(i,j)=-PE(i,j);
PF(i,j)=B(i,j)*e(i,1)-G(i,j)*f(i,1);
PF(i,j)=QE(i,j);
VE(i,j)=0;
VF(i,j)=0;
else
for x=1:n
PE(i,j)=PE(i,j)-(G(i,x)*e(x,1)-B(i,x)*f(x,1));
PF(i,j)=PF(i,j)-(G(i,x)*f(x,1)+B(i,x)*e(x,1));
QE(i,j)=QE(i,j)+(G(i,x)*f(x,1)+B(i,x)*e(x,1));
QF(i,j)=QF(i,j)-(G(i,x)*e(x,1)-B(i,x)*f(x,1));
VE(i,i)=-2*e(i,1);
VF(i,i)=-2*f(i,1);
end
PE(i,j)=PE(i,j)-G(i,i)*e(i,1)-B(i,i)*f(i,1);
PF(i,j)=PF(i,j)+B(i,i)*e(i,1)-G(i,i)*f(i,1);
QE(i,j)=QE(i,j)+B(i,i)*e(i,1)-G(i,i)*f(i,1);
QF(i,j)=QF(i,j)+G(i,i)*e(i,1)+B(i,i)*f(i,1);
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:2*(n-1)
for j=1:2*(n-1)
if rem(i,2)~=0 %奇数行
if rem(j,2)~=0 %为奇数行奇数列
a=(i+1)/2;
b=(j+1)/2;
J(i,j)=PE(a,b);
else %为奇数行偶数列
a=(i+1)/2;
b=j/2;
J(i,j)=PF(a,b);
end
else %偶数行
a=i/2;
if rem(j,2)~=0 %偶数行奇数列
b=(j+1)/2;
if B2(a,n-1)==3 %如果是PV节点的
J(i,j)=VE(a,b);
else %如果是PQ节点
J(i,j)=QE(a,b);
end
else %偶数行偶数列
b=j/2;
if B2(a,n-1)==3 %如果是PV节点的
J(i,j)=VF(a,b);
else %如果是PQ节点
J(i,j)=QF(a,b);
end
end
end
end
end
disp('输出雅格比矩阵')
disp(J);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%第五步雅格比矩阵和功率的偏移量形成增广矩阵JW
%说明这里的增广矩阵是定义为JW,在形成时乘以*(-1)
J=(-1)*J;
for i=1:12
for j=1:13
if j<=12
JW(i,j)=J(i,j);
else
if rem(i,2)~=0 %不论是PQ,PV节点,奇数行都是一样的
a=(i+1)/2;
JW(i,j)=XP(a,1);
else
a=i/2;
if B2(a,n-1)==3 %偶数行为PV节点
JW(i,j)=XV(a,1);
else %偶数行为PQ节点
JW(i,j)=XQ(a,1);
end
end
end
end
end
disp('输出增广矩阵JW')
disp(JW);
%%%%%%%%%%%%%%%%%%%%%%%%%