📄 zsglcx.m
字号:
%% x=[PG1 PG2 QG1 QG2 A1 A2 A3 A4 A5 V1 V2 V3 V4 V5]
function zsglcx
clc;
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 初始化 p ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%
n=5; %%%%%%%%%% 网络原始数据
Y=[1.37874194213052-6.29166544020251i -0.62402496099844+3.90015600624025i -0.75471698113208+2.64150943396226i 0
0;
-0.62402496099844+3.90015600624025i 1.45390047967064-66.98082109846432i -0.82987551867220+3.11203319502075i
63.49206349206349i 0 ;
-0.75471698113208+2.64150943396226i -0.82987551867220+3.11203319502075i 1.58459249980427-35.73785857758467i 0
31.74603174603175i;
0 63.49206349206349i 0 -66.66666666666667i 0 ;
0 0 31.74603174603175i 0
-33.33333333333334i];
G=real(Y);
B=imag(Y);
Ua=[1 1 1 1 1.05; %%%% 电压 相角 节点编号
0 0 0 0 0;
1 2 3 4 5];
V=Ua(1,:); %%%% 电压
A=Ua(2,:)*pi/180; %%%% 相角
PGG=[5 2.58]; %%%% 发电机有功
QGG=[0 1.45]; %%%% 发电机无功
PG=[0 0 0 PGG];
QG=[0 0 0 QGG];
PL=[1.6 2 3.7 0 0];
QL=[0.8 1 1.3 0 0];
P=PG-PL; %%%% 各节点注入功率
Q=QG-QL;
PT=zeros(1,10); %%%% 节点间传输功率
xxgg=zeros(1,9);
PGGg=zeros(2,1);
QGGg=zeros(2,1);
lyyy=zeros(1,10);
lx=zeros(13,1);
Lxx=zeros(13,1);
%%%%%%%%%%%%%%% 各限制条件
VM=[1.1 1.1 1.1 1.1 1.1;
0.9 0.9 0.9 0.9 0.9];
PGM=[8 8;
1 1];
QGM=[3 5;
-3 -2.1];
PTM=[2 2 0.65 0.65 1 1 6 6 5 5;
1 2 1 3 2 3 2 4 3 5; %%%%% 实际的 i j 编号
2 1 3 1 3 2 4 2 5 3];
for k=1:10
i=PTM(2,k);j=PTM(3,k);
s=sin(A(i)-A(j));c=cos(A(i)-A(j));
PT(k)=V(i)^2*G(i,j)-V(i)*V(j)*(G(i,j)*c+B(i,j)*s);
end
%%%%%%%%%%%%%%% 内点法数学计算参数
for i=1:4
xxgg(2*i-1)=A(i);
xxgg(2*i)=V(i);
end
xxgg(9)=V(5);
x=[PGG QGG xxgg];
x=x';
%%%%% g(x) 和 x
gmax=[PGM(1,:) QGM(1,:) VM(1,:) PTM(1,:)];gmax=gmax';
gmin=[PGM(2,:) QGM(2,:) VM(2,:) -PTM(1,:)];gmin=gmin';
l=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];l=l';
u=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];u=u';
z=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];z=z';
w=[-0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5 -0.5];w=w';
y=zeros(1,10);
for i=1:5
y(2*i-1)=1E-10;
y(2*i)=-1E-10;
end
y=y';
fb=[15 5;
20 8]; %%%发电机增加/削减出力时的报价
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 初始化 d ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 进入迭代 开始计算
g=[PGG QGG V PT];g=g'; %%% 更新g(x)
gap=0; %%% 计算 gap
gap=l'*z-u'*w;
kkk=0;
while gap>0.000001 %%%%%% 迭代公式开始
uu=0.1*gap/38; %disp(gap);
Shx=shxxx(A,V,G,B);
Sgx=sgxx(PTM,A,V,G,B) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 系数矩阵 p ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%%
%%%%%%%%% 对角矩阵 %%%%%%%%%
L=diag(l);
U=diag(u);
Z=diag(z);
W=diag(w);
%%%%%%%%%%% 海森伯矩阵 %%%%%%%%%%%
S2fx=zeros(13,13);
%%%%%%%%%%% 目标函数的海森矩阵
S2hx=s2hxx(A,y,V,G,B,n) ;
S2gx=s2gxx(PTM,A,V,G,B,z,w);
H=-S2fx+S2hx+S2gx;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 系数矩阵 d ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 常数项 p ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
for i=1:5
ppp(i)=0;qqq(i)=0;
for j=1:5
c=cos(A(i)-A(j));s=sin(A(i)-A(j));
ppp(i)=ppp(i)+V(j)*(G(i,j)*c+B(i,j)*s);
qqq(i)=qqq(i)+V(j)*(G(i,j)*s-B(i,j)*c);
end
lyyy(2*i-1)=PG(i)-PL(i)-V(i)*ppp(i);
lyyy(2*i)=QG(i)-QL(i)-V(i)*qqq(i);
end
e=ones(19,1);
Ly=lyyy';
Lz=g-l-gmin;
Lw=g+u-gmax;
Lul=L*Z*e-uu*e;
Luu=U*W*e+uu*e;
fpg=zeros(2,1);
fqr=zeros(2,1);
fx=zeros(9,1);
fpg(1)=-fb(1,2);
fpg(2)=fb(2,1);
Sfx=[fpg
fqr
fx]; %%%disp(Sfx);
lx=Sfx-Shx*y-Sgx*(z+w);
%%%% LX 矩阵形成。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 常数项 d ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ** 修正方程计算 d ** %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I=eye(19);
oo=zeros(19);o1=zeros(19,13);o2=zeros(13,19);o3=zeros(19,10);o4=zeros(10,19);o5=zeros(10);
%NH=[O O O I Shx' O];
NH=[L Z oo oo o1 o3;
oo -I oo oo Sgx' o3;
oo oo U W o1 o3;
oo oo oo I Sgx' o3;
Sgx o2 Sgx o2 H Shx;
o4 o4 o4 o4 Shx' o5] ;
NR=[-Lul
-Lz
-Luu
-Lw
lx
-Ly];
DX=inv(NH)*NR;
Dz=DX(1:19);
Dl=DX((1*19+1):(2*19));
Dw=DX((2*19+1):(3*19));
Du=DX((3*19+1):(4*19));
Dx=DX((4*19+1):89);
Dy=DX(90:99);
ap=1;ad=1;
t=1;
for i=1:19
if Dl(i)<0
t=-l(i)/Dl(i);
if ap>t
ap=t;
end
end
if Du(i)<0
t=-u(i)/Du(i);
if ap>t
ap=t;
end
end
if Dz(i)<0
t=-z(i)/Dz(i);
if ad>t
ad=t;
end
end
if Dw(i)>0
t=-w(i)/Dw(i);
if ad>t
ad=t;
end
end
end
ap=ap*0.9995;
ad=ad*0.9995;
%% ap=0.01;
%% ad=0.01;
x=x+ap*Dx;
l=l+ap*Dl;
u=u+ap*Du;
y=y+ad*Dy;
z=z+ad*Dz;
w=w+ad*Dw;
% disp(Ly);
% disp(Lz);
% disp(Lw);
% disp(Dx);
PGGg=x(1:2);PGG=PGGg';
QGGg=x(3:4);QGG=QGGg';
for i=3:6
pp=2*i-1;qq=2*i;
A(i-2)=x(pp);
V(i-2)=x(qq);
end
A(5)=0;
V(5)=x(13);
for k=1:10
i=PTM(2,k);j=PTM(3,k);
s=sin(A(i)-A(j));c=cos(A(i)-A(j));
PT(k)=V(i)^2*G(i,j)-V(i)*V(j)*(G(i,j)*c+B(i,j)*s);
end
g=[PGG QGG V PT];
g=g';
PG=[0 0 0 PGG];
QG=[0 0 0 QGG];
for i=1:4
xxgg(2*i-1)=A(i);
xxgg(2*i)=V(i);
end
xxgg(9)=V(5);
x=[PGG QGG xxgg];
x=x';
gap=l'*z-u'*w;
kkk=kkk+1;
if kkk>50
disp('计算不收敛');
break;
end
end
disp('发电机有功 发电机无功 角度 幅值 ');
disp(x);
disp('首端 末端 支路功率');
result =[PTM(2,:)',PTM(3,:)',PT'];
disp(result);
disp('阻塞管理费用');
cost=(fb(1,2)*(5-x(1,1))+fb(2,1)*(x(2,1)-2.58))*100;
disp(cost);
function Shx=shxxx(A,V,G,B)
%%%%%%% 等式约束的雅克比矩阵 %%%%%%%
hpg=[0 0 0 0 0 0 1 0 0 0;
0 0 0 0 0 0 0 0 1 0];
hqr=[0 0 0 0 0 0 0 1 0 0;
0 0 0 0 0 0 0 0 0 1];
%%%%常规潮流中的雅克比矩阵
for i=1:5
ha(i)=0;nb(i)=0;
for jj=1:5
c=cos(A(i)-A(jj));
s=sin(A(i)-A(jj));
nb(i)=nb(i)+V(jj)*(G(i,jj)*c+B(i,jj)*s);
ha(i)=ha(i)+V(jj)*(G(i,jj)*s-B(i,jj)*c);
end
for j=1:5
if j~=i
c=cos(A(j)-A(i));s=sin(A(j)-A(i));
pp=2*i-1;qq=2*j-1;
HH(pp,qq)=-V(i)*V(j)*(G(i,j)*s-B(i,j)*c);
qm=qq+1;
HH(pp,qm)=V(i)*V(j)*(G(i,j)*c+B(i,j)*s);
pm=pp+1;
HH(pm,qq)=-V(j)*(G(i,j)*c+B(i,j)*s);
HH(pm,qm)=-V(j)*(G(i,j)*s-B(i,j)*c);
else
pp=2*i-1;
HH(pp,pp)=ha(i)*V(i)+B(i,i)*V(i)^2;
pm=pp+1;
HH(pp,pm)=-nb(i)*V(i)+G(i,i)*V(i)^2;
HH(pm,pp)=-nb(i)-V(i)*G(i,i);
HH(pm,pm)=-ha(i)+V(i)*B(i,i);
end
end
end
hx=HH;
%disp(hx);
Shx=[hpg
hqr
hx];
Shx(13,:)=[];
function Sgx=sgxx(PTM,A,V,G,B)
%%%%%% 不等式约束的雅克比矩阵 %%%%%%%
g1g=[1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
g2q=[0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
g1x=zeros(10,2);
g2x=zeros(10,2);
g3x=zeros(10,5);
g3x0=zeros(5);
g3x1=eye(5);
for i=1:5
pp=2*i-1;qq=2*i;
g3x(pp,:)=g3x0(i,:);
g3x(qq,:)=g3x1(i,:);
end
g4x=zeros(10,10); g4x1=zeros(5,10); g4x2=zeros(5,10);
for k=1:10
i=PTM(2,k);j=PTM(3,k);
s=sin(A(i)-A(j));c=cos(A(i)-A(j));
g4x1(i,k)=V(i)*V(j)*(G(i,j)*s-B(i,j)*c);
g4x1(j,k)=-g4x1(i,k);
g4x2(i,k)=2*V(i)*G(i,j)-V(j)*(G(i,j)*c+B(i,j)*s);
g4x2(j,k)=-V(i)*(G(i,j)*c+B(i,j)*s);
end
for i=1:5
pp=2*i-1;qq=2*i;
g4x(pp,:)=g4x1(i,:);
g4x(qq,:)=g4x2(i,:);
end
%disp(g4x);
Sgx=[g1g;
g2q;
g1x g2x g3x g4x];
% disp(Sgx);
Sgx(13,:)=[];
%%%%%%
%%%%%%
%%%%%
function S2hx=s2hxx(A,y,V,G,B,n)
%%%%%%%% 等式约束海森伯矩阵和拉格朗日乘子的乘积 S2h(x)*y
S2hx=zeros(14,14);
for i=1:5
for j=1:5
if j==i
suma=0;sumb=0;
for k=1:5
if k~=i
s=sin(A(i)-A(k));c=cos(A(i)-A(k));
y1=y(2*i-1);y2=y(2*i);y3=y(2*k-1);y4=y(2*k);
suma=suma+V(k)*(G(i,k)*(c*y1+s*y2+c*y3-s*y4)+B(i,k)*(s*y1-c*y2-s*y3-c*y4));
sumb=sumb+V(k)*(G(i,k)*(s*y1-c*y2+s*y3+c*y4)+B(i,k)*(-c*y1-s*y2+c*y3-s*y4));
end
end
pp=2*i-1;
suma=suma*V(i);
S2hxa(pp,pp)=suma;
pm=pp+1;
S2hxa(pp,pm)=sumb;
S2hxa(pm,pp)=sumb;
S2hxa(pm,pm)=-2*(G(i,i)*y(2*i-1)-B(i,i)*y(2*i));
else
s=sin(A(i)-A(j));c=cos(A(i)-A(j));
y1=y(2*i-1);y2=y(2*i);y3=y(2*j-1);y4=y(2*j);
pp=2*i-1;qq=2*j-1;
S2hxa(pp,qq)=V(i)*V(j)*(G(i,j)*(-c*y1-s*y2-c*y3+s*y4)+B(i,j)*(-s*y1+c*y2+s*y3+c*y4));
qm=qq+1;
S2hxa(pp,qm)=V(i)*(G(i,j)*(s*y1-c*y2+s*y3+c*y4)+B(i,j)*(-c*y1-s*y2+c*y3-s*y4));
pm=pp+1;
S2hxa(pm,qq)=-V(j)*(G(i,j)*(s*y1-c*y2+s*y3+c*y4)+B(i,j)*(-c*y1-s*y2+c*y3-s*y4));
S2hxa(pm,qm)=-(G(i,j)*(c*y1+s*y2+c*y3-s*y4)+B(i,j)*(s*y1-c*y2-s*y3-c*y4));
end
end
end
S2hx0=S2hxa;
for i=5:14
for j=5:14
S2hx(i,j)=S2hx0(i-4,j-4);
end
end
S2hx(13,:)=[];
S2hx(:,13)=[];
%
%%%%%%
%%%%%%
function S2gx=s2gxx(PTM,A,V,G,B,z,w)
%%%%%%%%%%% 不等式约束海森伯矩阵与拉格朗日乘子Z+W的乘积
S2gx=zeros(14,14);
S2gx0=zeros(10,10);
S2gxa=zeros(5,5);S2gxb=zeros(5,5);S2gxc=zeros(5,5);S2gxd=zeros(5,5);
C=z+w;
for kk=1:10
i=PTM(2,kk);j=PTM(3,kk);cc=C(2+2+5+kk);
s=sin(A(i)-A(j));c=cos(A(i)-A(j));
DS2gxa1=V(i)*V(j)*(G(i,j)*c+B(i,j)*s); %%%%%%%%%%% SaSa
DS2gxa2=-DS2gxa1;
S2gxa(i,i)=S2gxa(i,i)+DS2gxa1*cc;
S2gxa(i,j)=S2gxa(i,j)+DS2gxa2*cc;
S2gxa(j,i)=S2gxa(j,i)+DS2gxa2*cc;
S2gxa(j,j)=S2gxa(j,j)+DS2gxa1*cc;
DS2gxb1=V(j)*(G(i,j)*s-B(i,j)*c); %%%%%%%%%%% SaSv
DS2gxb2=V(i)*(G(i,j)*s-B(i,j)*c);
S2gxb(i,i)=S2gxb(i,i)+DS2gxb1*cc;
S2gxb(i,j)=S2gxb(i,j)+DS2gxb2*cc;
S2gxb(j,i)=S2gxb(j,i)-DS2gxb1*cc;
S2gxb(j,j)=S2gxb(j,j)-DS2gxb2*cc;
DS2gxd1=2*G(i,j); %%%%%%%%%%% SvSv
DS2gxd2=-(G(i,j)*c+B(i,j)*s);
S2gxd(i,i)=S2gxd(i,i)+DS2gxd1*cc;
S2gxd(i,j)=S2gxd(i,j)+DS2gxd2*cc;
S2gxd(j,i)=S2gxd(j,i)+DS2gxd2*cc;
end
S2gxc=S2gxb';
for i=1:5
for j=1:5
pp=2*i-1;
qq=2*j-1;
pm=2*i;
qm=2*j;
S2gx0(pp,qq)=S2gxa(i,j);
S2gx0(pp,qm)=S2gxb(i,j);
S2gx0(pm,qq)=S2gxc(i,j);
S2gx0(pm,qm)=S2gxd(i,j);
end
end
for i=5:14
for j=5:14
S2gx(i,j)=S2gx0(i-4,j-4);
end
end
S2gx(13,:)=[];
S2gx(:,13)=[];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -