function[,S_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);
loss_marix=diag([S_loss_active/p2,S_loss_active/p3,S_loss_active/p5]);
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'
%如通過支路功率求和來求取網絡的總損耗，容易出現正負號的問題，所以選擇用網絡所有節點注入功率之和來求網損