基于matlab实现简化梯度法进行电力系统最优潮流计算 .rar

clear all; clf; y=[0-66.66667i 0 0+63.49206i 0 0; 0 0-33.33333i 0 0+31.74603i 0; 0+63.49206i 0 1.45390-66.98082i -0.82987+3.11203i -0.62402+3.90015i; 0 0+31.74603i -0.82987+3.11203i 1.58459-35.73786i -0.75471+2.64150i; 0 0 -0.62402+3.90015i -0.75471+2.64150i 1.37874-6.29166i] ;%导纳矩阵 g=real(y); b=imag(y); u(:,1)=[1.05;1.05;1;1;1]; %各节点初始电压值 ct(:,1)=[0;0;0;0;0]; %初始角度 p0=[1;0;-0.4;-0.74;-0.32]; %各节点的输入有功功率 q0=[0;0;-0.2;-0.26;-0.16]; %各节点的输入无功功率 x(:,1)=[0;1;0;1;0;1;0]; %控制变量的初始值 e=0.0001; %迭代精度 zu(:,1)=[1;1;1;1.05]; %最优控制变量解对应为P1、U1、U2、T for k=1:1000 for i=1:5 a=0; for j=1:5 a=u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k)))+a; end dp(i,k)=p0(i)-u(i,k)*a; end for i=3:5 z=0; for j=1:5 z=u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k)))+z; end dq(i,k)=q0(i)-u(i,k)*z; end for i=1:5 for j=1:5 if i~=j hh(i,j)=u(i,k)*u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k))); jj(i,j)=-u(i,k)*u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k))); nn(i,j)=u(i,k)*u(j,k)*(g(i,j)*cos(ct(i,k)-ct(j,k))+b(i,j)*sin(ct(i,k)-ct(j,k))); ll(i,j)=u(i,k)*u(j,k)*(g(i,j)*sin(ct(i,k)-ct(j,k))-b(i,j)*cos(ct(i,k)-ct(j,k))); else c=0; d=0; for m=1:5 if m~=i c=u(m,k)*(g(i,m)*sin(ct(i,k)-ct(m,k))-b(i,m)*cos(ct(i,k)-ct(m,k)))+c; d=u(m,k)*(g(i,m)*cos(ct(i,k)-ct(m,k))+b(i,m)*sin(ct(i,k)-ct(m,k)))+d; end end hh(i,i)=-u(i,k)*c; jj(i,i)=u(i,k)*d; nn(i,i)=u(i,k)*d+2*u(i,k)^2*g(i,i); ll(i,i)=u(i,k)*c-2*u(i,k)^2*b(i,i); end end end jkb=[hh(3,3) nn(3,3) hh(3,4) nn(3,4) hh(3,5) nn(3,5) hh(3,1); jj(3,3) ll(3,3) jj(3,4) ll(3,4) jj(3,5) ll(3,5) jj(3,1); hh(4,3) nn(4,3) hh(4,4) nn(4,4) hh(4,5) nn(4,5) hh(4,1); jj(4,3) ll(4,3) jj(4,4) ll(4,4) jj(4,5) ll(4,5) jj(4,1); hh(5,3) nn(5,3) hh(5,4) nn(5,4) hh(5,5) nn(5,5) hh(5,1); jj(5,3) ll(5,3) jj(5,4) ll(5,4) jj(5,5) ll(5,5) jj(5,1); hh(1,3) nn(1,3) hh(1,4) nn(1,4) hh(1,5) nn(1,5) hh(1,1)]; %牛拉法的雅各比矩阵 dd=[dp(3,k);dq(3,k);dp(4,k);dq(4,k);dp(5,k);dq(5,k);dp(1,k)]; dx(:,k)=inv(jkb)*dd; xzh=[1 0 0 0 0 0 0 ;0 u(3,k) 0 0 0 0 0;0 0 1 0 0 0 0;0 0 0 u(4,k) 0 0 0;0 0 0 0 1 0 0;0 0 0 0 0 u(5,k) 0;0 0 0 0 0 0 1]; x(:,k+1)=x(:,k)+xzh*dx(:,k); %潮流计算解 u(:,k+1)=[u(1,k);u(2,k);x(2,k+1);x(4,k+1);x(6,k+1)]; ct(:,k+1)=[x(7,k+1);0;x(1,k+1);x(3,k+1);x(5,k+1)]; f1=0; for i=1:5 f1=u(j,k+1)*(g(1,j)*cos(ct(1,k+1)-ct(j,k+1))+b(1,j)*sin(ct(1,k+1)-ct(j,k+1)))+f1; end pg1=u(1,k+1)*f1; f2=0; for i=1:5 f2=u(j,k+1)*(g(2,j)*cos(ct(2,k+1)-ct(j,k+1))+b(2,j)*sin(ct(2,k+1)-ct(j,k+1)))+f2; end pg2=u(2,k+1)*f2; fdx=[(702*pg1+100)*u(1,k+1)*u(3,k+1)*(g(1,3)*sin(ct(1,k+1)-ct(3,k+1))-b(1,3)*cos(ct(1,k+1)-ct(3,k+1)))+(778*pg2+100)*u(2,k+1)*u(3,k+1)*(g(2,3)*sin(ct(2,k+1)-ct(3,k+1))-b(2,3)*cos(ct(2,k+1)-ct(3,k+1))); (702*pg1+100)*u(1,k+1)*(g(1,3)*cos(ct(1,k+1)-ct(3,k+1))+b(1,3)*sin(ct(1,k+1)-ct(3,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,3)*cos(ct(2,k+1)-ct(3,k+1))+b(2,3)*sin(ct(2,k+1)-ct(3,k+1))); (702*pg1+100)*u(1,k+1)*u(4,k+1)*(g(1,4)*sin(ct(1,k+1)-ct(4,k+1))-b(1,4)*cos(ct(1,k+1)-ct(4,k+1)))+(778*pg2+100)*u(2,k+1)*u(4,k+1)*(g(2,4)*sin(ct(2,k+1)-ct(4,k+1))-b(2,4)*cos(ct(2,k+1)-ct(4,k+1))); (702*pg1+100)*u(1,k+1)*(g(1,4)*cos(ct(1,k+1)-ct(4,k+1))+b(1,4)*sin(ct(1,k+1)-ct(4,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,4)*cos(ct(2,k+1)-ct(4,k+1))+b(2,4)*sin(ct(2,k+1)-ct(4,k+1))); (702*pg1+100)*u(1,k+1)*u(5,k+1)*(g(1,5)*sin(ct(1,k+1)-ct(5,k+1))-b(1,5)*cos(ct(1,k+1)-ct(5,k+1)))+(778*pg2+100)*u(2,k+1)*u(5,k+1)*(g(2,5)*sin(ct(2,k+1)-ct(5,k+1))-b(2,5)*cos(ct(2,k+1)-ct(5,k+1))); (702*pg1+100)*u(1,k+1)*(g(1,5)*cos(ct(1,k+1)-ct(5,k+1))+b(1,5)*sin(ct(1,k+1)-ct(5,k+1)))+(778*pg2+100)*u(2,k+1)*(g(2,5)*cos(ct(2,k+1)-ct(5,k+1))+b(2,5)*sin(ct(2,k+1)-ct(5,k+1))); (702*pg1+100)*u(1,k+1)*u(1,k+1)*(g(1,1)*sin(ct(1,k+1)-ct(1,k+1))-b(1,1)*cos(ct(1,k+1)-ct(1,k+1)))+(778*pg2+100)*u(2,k+1)*u(1,k+1)*(g(2,1)*sin(ct(2,k+1)-ct(1,k+1))-b(2,1)*cos(ct(2,k+1)-ct(1,k+1)))]; lmd=-inv((jkb)')*fdx; aa=g(2,1)*cos(ct(2,k+1)-ct(1,k+1))+b(2,1)*sin(ct(2,k+1)-ct(1,k+1)); tdf=[702*pg1+100;878*pg2*u(2,k+1)*aa;878*(u(1,k+1)*aa+2*u(2,k+1)*g(2,2));0]; if norm(tdf)<=e %精度判别 pg1=zu(1,k) u1=zu(2,k) u2=zu(3,k) T=zu(4,k) break; end bch=1.2; %步长因子 zu(:,k+1)=zu(:,k)-bch*tdf; if zu(1,k+1)>=1.2 %不等式约束条件 zu(1,k+1)=1.2; elseif zu(1,k+1)<=0.3 zu(1,k+1)=0.3; end zu(1,k+1)=1; if zu(2,k+1)>=1.1 zu(2,k+1)=1.1; elseif zu(2,k+1)<=1.0 end zu(2,k+1)=1.02; if zu(3,k+1)>=1.1 zu(3,k+1)=1.1; elseif zu(3,k+1)<=1.0 zu(3,k+1)=1.0; end zu(3,k+1)=1.02; u(1,k+1)=zu(2,k+1); u(2,k+1)=zu(3,k+1); p0(1)=zu(1,k+1); end u=u(:,k+1) %输出结果 ct=ct(:,k+1) x=x(:,k+1) zu=zu(:,k)