clc
clear
tic
%% *************************读取数据*************************
mpc=case5;
ng=size(mpc.gen,1); %发电机数2
nn=size(mpc.bus,1); %节点数5
nb=size(mpc.branch,1); %支路数5
%节点数据提取
Pd=mpc.bus(:,3);
Qd=mpc.bus(:,4);
V=mpc.bus(:,5);
delta=mpc.bus(:,6);
%发电机数据提取
PG=zeros(nn,1);
QG=zeros(nn,1);
PG(nn-ng+1:nn,1)=mpc.gen(:,3);
QG(nn-ng+1:nn,1)=mpc.gen(:,4);
a2=mpc.gen(:,9); %燃料耗费曲线二次系数
al=mpc.gen(:,10); %燃料耗费曲线一次系数
a0=mpc.gen(:,11); %燃料耗费曲线常数项
%% ************************参数初始化*************************
equ_cons=2*nn; %等式约束个数10
nonequ_cons=2*ng+nn+nb;%不等式约束个数14
xn=2*(ng+nn);%系统总变量个数14
%系统总变量
x=zeros(xn,1);
x(1:ng)=PG(nn-ng+1:nn,1);
x(ng+1:2*ng)=QG(nn-ng+1:nn,1);
x((2*ng+2):2:2*(ng+nn))=V;
x((2*ng+1):2:2*(ng+nn)-1)=delta;
%松弛变量
l=ones(nonequ_cons,1);
u=ones(nonequ_cons,1);
%拉格朗日乘子
y=zeros(equ_cons,1);
y(1:2:2*nn-1)=1e-10;
y(2:2:2*nn)=-1e-10;
z=ones(nonequ_cons,1);
w=-0.5*ones(nonequ_cons,1);
%gmax
gmax=zeros(nonequ_cons,1);
gmax(1:ng)=mpc.gen(:,5);
gmax(ng+1:2*ng)=mpc.gen(:,6);
gmax(2*ng+1:2*ng+nn)=mpc.bus(:,7);
gmax(2*ng+nn+1:2*ng+nn+nb)=mpc.branch(:,8);
%gmin
gmin=zeros(nonequ_cons,1);
gmin(1:ng)=mpc.gen(:,7);
gmin(ng+1:2*ng)=mpc.gen(:,8);
gmin(2*ng+1:2*ng+nn)=mpc.bus(:,8);
gmin(2*ng+nn+1:2*ng+nn+nb)=-mpc.branch(:,8);
%% ********************形成节点导纳矩阵Y***********************
Y=Ybus_build(mpc.branch);
G=real(Y);
B=imag(Y);
% ****关联矩阵****
INM=zeros(nn,nb);
for i=1:5
INM(mpc.branch(i,2),i)=1;
INM(mpc.branch(i,3),i)=-1;
end
%%
k=0;%迭代次数
Gap_save=zeros(1,50);
dx_save=zeros(50,14);
Gap=l'*z-u'*w;
while k<50&&Gap>=1e-6
L=diag(l);
U=diag(u);
Z=diag(z);
W=diag(w);
Gap=l'*z-u'*w; %计算互补间隙Gap
Gap_save(1,k+1)=Gap;
miu=0.1*Gap/(2*nonequ_cons);
theta=zeros(nn,nn);
for i=1:nn
for j=1:nn
theta(i,j)=delta(i)-delta(j);
end
end
%% *********************形成系数矩阵**********************
%%----- ①等式约束的雅克比矩阵 -------%%
J=zeros(2*nn,2*nn);
for n=1:nn
for m=1:nn
if m~=n
J(2*n-1,2*m-1)=-V(m)*V(n)*(G(m,n)*sin(theta(m,n))-B(m,n)*cos(theta(m,n))); %
J(2*n,2*m-1)=-V(m)*(G(m,n)*cos(theta(m,n))+B(m,n)*sin(theta(m,n)));
J(2*n-1,2*m)=V(m)*V(n)*(G(m,n)*cos(theta(m,n))+B(m,n)*sin(theta(m,n)));%
J(2*n,2*m)=-V(m)*(G(m,n)*sin(theta(m,n))-B(m,n)*cos(theta(m,n)));
J(2*n-1,2*n-1)=J(2*n-1,2*n-1)+V(m)*V(n)*(G(n,m)*sin(theta(n,m))-B(n,m)*cos(theta(n,m)));
J(2*n-1,2*n)=J(2*n-1,2*n)-V(m)*V(n)*(G(n,m)*cos(theta(n,m))+B(n,m)*sin(theta(n,m)));
J(2*n,2*n-1)=J(2*n,2*n-1)-V(m)*(G(n,m)*cos(theta(n,m))+B(n,m)*sin(theta(n,m)));
J(2*n,2*n)=J(2*n,2*n)-V(m)*(G(n,m)*sin(theta(n,m))-B(n,m)*cos(theta(n,m)));
end
end
J(2*n,2*n-1)=J(2*n,2*n-1)-2*V(n)*G(n,n);
J(2*n,2*n)=J(2*n,2*n)+2*V(n)*B(n,n);
end
%形成hPG，hQR
hPG=zeros(ng,2*nn);
hQR=zeros(ng,2*nn);
for qq=1:nn
for ee=1:ng
if (ee==1&&qq==4)||(ee==2&&qq==5)
hPG(ee,2*qq-1)=1;
hQR(ee,2*qq)=1;
end
end
end
hx1=vertcat(hPG,hQR,J);
%% %%----- ②不等式约束的雅克比矩阵 -------%%
g1PG=eye(ng);
g2QR=eye(ng);
g3x=zeros(2*nn,nn);
for m=1:nn
g3x(2*m,m)=1;
end
g4x=zeros(2*nn,nn);
for j=1:nb
m=mpc.branch(j,2);
n=mpc.branch(j,3);
g4x(2*m-1,j)=-V(m)*V(n)*(G(m,n)*sin(theta(m,n))-B(m,n)*cos(theta(m,n)));
g4x(2*m,j)=V(n)*(G(m,n)*cos(theta(m,n))+B(m,n)*sin(theta(m,n)))-2*V(m)*G(m,n);
g4x(2*n-1,j)=V(m)*V(n)*(G(m,n)*sin(theta(m,n))-B(m,n)*cos(theta(m,n)));
g4x(2*n,j)=V(m)*(G(m,n)*cos(theta(m,n))+B(m,n)*sin(theta(m,n)));
end
gx1=horzcat(blkdiag(g1PG,g2QR,g3x),vertcat(zeros(4,5),g4x));
%% %%----- ③对角矩阵 -------%%
L_Z=L\Z;
U_W=U\W;
%% %%----- ④海森伯矩阵 -------%%
%a.目标函数的海森伯矩阵--HH1
A2=diag(a2);
aafx=zeros(xn,xn);
aafx(1:ng,1:ng)=2*A2;
HH1=aafx;
%b.等式约束海森伯矩阵hx2与拉格朗日y的乘积--HH2
A=zeros(2*nn,2*nn);
for i=1:nn
for n=1:nn
if i==n %%%对角矩阵块
for j=1:nn
if j~=i
A(2*i-1,2*n-1)=A(2*i-1,2*n-1)+V(i)*V(j)*(G(i,j)*(cos(angle(i)-angle(j))*y(2*i-1)+sin(angle(i)-angle(j))*y(2*i)+cos(angle(i)-angle(j))*y(2*j-1)-sin(angle(i)-angle(j))*y(2*j)) ...
+B(i,j)*(sin(angle(i)-angle(j))*y(2*i-1)-cos(angle(i)-angle(j))*y(2*i)-sin(angle(i)-angle(j))*y(2*j-1)-cos(angle(i)-angle(j))*y(2*j)));
end
end
for j=1:nn
if j~=i
A(2*i-1,2*n)=A(2*i-1,2*n)+V(j)*(G(i,j)*(sin(angle(i)-angle(j))*y(2*i-1)-cos(angle(i)-angle(j))*y(2*i)+sin(angle(i)-angle(j))*y(2*j-1)+cos(angle(i)-angle(j))*y(2*j)) ...
+B(i,j)*(-cos(angle(i)-angle(j))*y(2*i-1)-sin(angle(i)-angle(j))*y(2*i)+cos(angle(i)-angle(j))*y(2*j-1)-sin(angle(i)-angle(j))*y(2*j)));
end
end
A(2*i,2*n)=-2*(G(i,i)*y(2*i-1)-B(i,i)*y(2*i));
A(2*i,2*n-1)= A(2*i-1,2*n);
else %非对角块
A(2*i-1,2*n-1)=V(i)*V(n)*(G(i,n)*(-cos(angle(i)-angle(n))*y(2*i-1)-sin(angle(i)-angle(n))*y(2*i)-cos(angle(i)-angle(n))*y(2*n-1)+sin(angle(i)-angle(n))*y(2*n)) ...
+B(i,n)*(-sin(angle(i)-angle(n))*y(2*i-1)+cos(angle(i)-angle(n))*y(2*i)+sin(angle(i)-angle(n))*y(2*n-1)+cos(angle(i)-angle(n))*y(2*n)));
A(2*i-1,2*n)=V(i)*(G(i,n)*(sin(angle(i)-angle(n))*y(2*i-1)-cos(angle(i)-angle(n))*y(2*i)+sin(angle(i)-angle(n))*y(2*n-1)+cos(angle(i)-angle(n))*y(2*n)) ...
+B(i,n)*(-cos(angle(i)-angle(n))*y(2*i-1)-sin(angle(i)-angle(n))*y(2*i)+cos(angle(i)-angle(n))*y(2*n-1)-sin(angle(i)-angle(n))*y(2*n)));
A(2*i,2*n-1)=V(n)*(G(i,n)*(-sin(angle(i)-angle(n))*y(2*i-1)+cos(angle(i)-angle(n))*y(2*i)-sin(angle(i)-angle(n))*y(2*n-1)-cos(angle(i)-angle(n))*y(2*n)) ...
+B(i,n)*(cos(angle(i)-angle(n))*y(2*i-1)+sin(angle(i)-angle(n))*y(2*i)-cos(angle(i)-angle(n))*y(2*n-1)+sin(angle(i)-angle(n))*y(2*n)));
A(2*i,2*n)=-(G(i,n)*(cos(angle(i)-angle(n))*y(2*i-1)+sin(angle(i)-angle(n))*y(2*i)+cos(angle(i)-angle(n))*y(2*n-1)+-sin(angle(i)-angle(n))*y(2*n)) ...
+B(i,n)*(sin(angle(i)-angle(n))*y(2*i-1)-cos(angle(i)-angle(n))*y(2*i)-sin(angle(i)-angle(n))*y(2*n-1)-cos(angle(i)-angle(n))*y(2*n)));
end
end
end
HH2=blkdiag(zeros(4,4),A);
%c.不等式约束海森伯矩阵gx2与拉格朗日c=z+w的乘积--HH3
HH3=zeros(xn,xn);
AA=zeros(2*nn,2*nn);
c=z+w;
a1=mpc.branch(:,2);
b1=mpc.branch(:,3);
for i=1:nb
for j=1:nb
AA(2*i-1,2*i-1)=AA(2*i-1,2*i-1)-V(b1(j))*V(a1(j))*(G(a1(j),b1(j))*cos(angle(a1(j))-angle(b1(j)))+B(a1(j),b1(j))*sin(angle(a1(j))-angle(b1(j))))*c(j+9)*abs(INM(i,j));
end
for j=1:nb
if INM(i,j)==1
AA(2*i-1,2*i)=AA(2*i-1,2*i)-V(b1(j))*(G(a1(j),b1(j))*sin(angle(a1(j))-angle(b1(j)))-B(a1(j),b1(j))*cos(angle(a1(j))-angle(b1(j))))*c(j+9)*INM(i,j);
else
AA(2*i-1,2*i)=AA(2*i-1,2*i)-V(a1(j))*(G(a1(j),b1(j))*sin(angle(a1(j))-angle(b1(j)))-B(a1(j),b1(j))*cos(angle(a1(j))-angle(b1(j))))*c(j+9)*INM(i,j);
end
end
AA(2*i,2*i-1)=AA(2*i-1,2*i);
for j=1:nb
AA(2*i,2*i)=AA(2*i,2*i)-2*G(a1(j),b1(j))*c(j+9)*(INM(i,j)>0);
end
end
for j=1:nb %非对角矩阵快
AA(2*a1(j)-1,2*b1(j)-1)=V(b1(j))*V(a1(j))*(G(a1(j),b1(j))*cos(angle(a1(j))-angle(b1(j)))+B(a1(j),b1(j))*sin(angle(a1(j))-angle(b1(j))))*c(j+9);
AA(2*a1(j)-1,2*b1(j))=V(a1
评论1