069用matlab编写的经典电力系统经济调度程序.rar

function[S_loss_active,loss_active,loss_marix]=loss(p2,p3,p5) A=[1,0,0,0,1.06,0; 2,1,p2-0.20,0.20,1,0; 3,1,p3-0.45,-0.15,1,0; 4,1,-0.40,-0.05,1,0; 5,1,p5-0.60,-0.10,1,0] tic; [dimention,dd]=size(A); y=zeros(dimention,dimention); y(1,3)=1/(0.08+j*0.24);y(3,4)=1/(0.01+j*0.03);y(4,5)=1/(0.08+j*0.24);y(2,5)=1/(0.04+j*0.12); y(2,4)=1/(0.06+j*0.18);y(2,3)=1/(0.06+j*0.18);y(1,2)=1/(0.02+j*0.06); y1=transpose(y); Y1=y+y1; [c,d]=size(Y1); %Y=[]; for m=1:c for n=1:d if(m~=n) Y(m,n)=-Y1(m,n); else sum=0; for k=1:d sum=Y1(m,k)+sum; Y(m,n)=sum; end end end end Y G=real(Y); B_temp=imag(Y); eps=0.00001; Akk_Ak=0.1; K=0; while (Akk_Ak>eps) K=K+1; Ak=A; [dP_temp]=make_in_P(A,G,B_temp);%求取多阶的有功功率不平衡量 B1=B_temp(2:5,2:5);% 电纳矩阵减阶 dP=dP_temp(2:5,:)%有功不平衡量减去平衡节点的不平衡量 B=inv(B1); U=A(2:5,5); for i=1:4 dP_U(i,1)=dP(i,1)*(1/U(i,1)); end U_a=-B*dP_U; for i=1:4 dHD(i,1)=U_a(i,1)/U(i,1);%相角的修正量 end dHD %显示弧度的修正量 for i=2:5 A(i,6)=A(i,6)+dHD(i-1,1);%修正后的相角值更新在矩阵A中 end [dQ_temp]=make_in_Q(A,G,B_temp); dQ=dQ_temp(2:5,1)%无功功率的不平衡量 for i=1:4 dQ_U(i,1)=dQ(i,1)*(1/U(i,1)); end dU=-B*dQ_U %电压幅值的修正量 for i=2:5 A(i,5)=A(i,5)+dU(i-1,1);%修正后的电压幅值更新在矩阵A中 end fprintf('第%d次迭代结果\n',K) Akk=A for m=2:5 for n=5:6 Akk_Ak_temp1(m,n)=abs(Akk(m,n)-Ak(m,n)); end end Akk_Ak_temp=max(Akk_Ak_temp1); Akk_Ak=max(Akk_Ak_temp'); end t1=toc; fprintf('迭代耗时%f\n秒',t1) %以下是计算平衡节点功率 [balance_P]=make_in_P(Akk,G,B_temp); [balance_Q]=make_in_Q(Akk,G,B_temp); i=sqrt(-1); S_balance=-balance_P(1,1)-balance_Q(1,1)*i %以下是计算线路功率 S_mn=[]; [k3,k4]=size(Y); for m=1:k3 for n=1:k4 S_mn(m,n)=-Akk(m,5)*(cos(Akk(m,6))+sin(Akk(m,6))*i)*((Akk(m,5)*(cos(-Akk(m,6))+sin(-Akk(m,6))*i)-Akk(n,5)*(cos(-Akk(n,6))+sin(-Akk(n,6))*i)))*conj(Y(m,n)); %注意公式互导纳加负号 % S_mn(m,n)=G(m,n)*(Akk(m,5)*cos(Akk(m,6))-Akk(n,5)*cos(Akk(n,6)))-B_temp(m,n)*(Akk(m,5)*sin(Akk(m,6))-Akk(n,5)*sin(Akk(n,6)))-(G(m,n)*(Akk(m,5)*sin(Akk(m,6))-Akk(n,5)*sin(Akk(n,6)))+B_temp(m,n)*(Akk(m,5)*cos(Akk(m,6))-Akk(n,5)*cos(Akk(n,6))))*i; end end S_mn %以下是通过各节点注入功率求和得到网络的总损耗S_loss， [k5,k6]=size(Akk); S_loss_temp=0; for m=1:k5 S_loss_temp=S_loss_temp+Akk(m,3)+Akk(m,4)*i; end S_loss=S_loss_temp+S_balance S_loss_active=real(S_loss); P=p2+p3+p5; loss_marix=diag([p2/P,p3/P,p5/P]); loss_active=[S_loss_active;S_loss_active;S_loss_active]; t2=toc fprintf('PQ分解法运算共耗时%f\n秒',t2) %S_mn_1=S_mn+S_mn' %如通过支路功率求和来求取网络的总损耗，容易出现正负号的问题，所以选择用网络所有节点注入功率之和来求网损