最优潮流简化梯度法计算程序

1、最优潮流简化梯度法计算程序7节点参数:线路参数IJRX130.100.4140.150.6150.050.2230.050.2250.050.2340.100.4程序:clear all;elf;y=2.1569-8.6274i-1.1765+4.7059i;-0.5882+2.3529i-0.3922+1.5686i0-1.1765+4.7059i;2.3530-9.4118i-1.1765+4.7059i-0.5882+2.3529i-1.1765+4.7059i2.3529-9.4117i-0.5882+2.3529i-0.3922+1.5686i-0.5882+2.3529i0.980

2、4-3.9215i-1.1765+4.7059i-1.1765+4.7059i0;0;02.3530-9.4118i;% 导纳矩阵g=real(y);%各节点初始电压值b=imag(y); u(:,1)=1;1;1;1;1;ct(:,1)=0;0;0;0;0;%初始角度p0=1;0;-0.6;-0.4;-0.6;%各节点的输入有功功率%控制变量的初始值%迭代精度%最优控制变量解对应为q0=0.5;0;-0.3;-0.1;-0.2; % 各节点的输入无功功率 x(:,1)=0;1;0;1;0;1;0;P1、U1、U2、Te=0.0001;zu(:,1)=1;1;1;1;for k=1:1000f

3、or i=1:5a=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;enddp(i,k)=p0(i)-u(i,k)*a;endfor i=3:5z=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;enddq(i,k)=q0(i)-u(i,k)*z;endfor i=1:5for j=1:5if i=j hh(i,j)=u(i,k)*u(j,k)*(g(i,j)*sin(ct(i,k

4、)-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:5if m=i

5、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; endend hh(i,i)=-u(i,k)*c;jj(i,i)=u(i,k)*d;nn (i,i)=u(i,k)*d+2*u(i,kF2*g(i,i);ll(i,i)=u(i,k)*c-2*u(i,k)A2*b(i,i);endendendjkb=hh(3,3)nn(3,3)hh(3,4)nn(3,4) hh(3,5)nn(3,

6、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)

7、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)

8、=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

9、 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)+(

10、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)

11、*(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)

12、-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)-

13、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)=1.2 % 不等式约束条件zu(1,k+1)=1

14、.2;elseif zu(1,k+1)=1.1zu(2,k+1)=1.1;elseif zu(2,k+1)=1.1zu(3,k+1)=1.1;elseif zu(3,k+1)=1.0zu(3,k+1)=1.0;endzu(3,k+1)=1.089;u(1,k+1)=zu(2,k+1);u(2,k+1)=zu(3,k+1);p0(1)=zu(1,k+1);endu=u(:,k+1) %输出结果ct=ct(:,k+1)x=x(:,k+1)zu=zu(:,k)计算结果：15各节点的电压值:u =1.09201.08901.01930.99541.0561 15 各节点的角度值:ct =0.00570-0.0770-0.1290-0