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'
%如通过支路功率求和来求取网络的总损耗,容易出现正负号的问题,所以选择用网络所有节点注入功率之和来求网损