threephasese-onlyenquationconstraint-gn/SolveIt.m

function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,PD0,PD,QD0,QD,mVolt,Volt,wVolt,wPD,wQD,Y,Yangle,UAngel,r,c)
% dH_dx
dP_x = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);  %<EFBFBD>γ<EFBFBD><EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
dP_x=-dP_x(Loadi,:);
dQ_x = jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
dQ_x=-dQ_x(Loadi+Busnum*3,:);
dV_x=[eye(Busnum*3,Busnum*3),zeros(Busnum*3,Busnum*3)];
H=[dP_x;dQ_x;dV_x];
cP=jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
cP=cP(setdiff(1:Busnum*3,Loadi),:);
cQ=jacobian_M(Busnum,Volt,Y,Yangle,UAngel,r,c);
cQ=cQ(setdiff(1:Busnum*3,Loadi)+Busnum*3,:);
C=[cP;cQ];
aa=[H'*diag([wPD;wQD;wVolt])*H C'
    C zeros(size(C,1),size(C',2))];
AngleIJ=sparse(r,c,UAngel(r)-UAngel(c)-Yangle,Busnum*3,Busnum*3);
PD=-diag(Volt)*Y.*cos(AngleIJ)*Volt;
PD=PD(Loadi);
QD=-diag(Volt)*Y.*sin(AngleIJ)*Volt;
QD=QD(Loadi);
dP=diag(Volt)*Y.*cos(AngleIJ)*Volt;
dQ=diag(Volt)*Y.*sin(AngleIJ)*Volt;
deltZ=[(PD0-PD);
    (QD0-QD);
    (mVolt-Volt);
    ];
yy=[
    H'*diag([wPD;wQD;wVolt])*deltZ;
    -dP(setdiff(1:Busnum*3,Loadi));
    -dQ(setdiff(1:Busnum*3,Loadi));
];
% yy=[
%     H'*eye(length(H))*deltZ;
% ];
%% ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD> 
t=0;
aa(t+(Balance-1)*3+1,:)=0;
aa(t+(Balance-1)*3+2,:)=0;
aa(t+(Balance-1)*3+3,:)=0;
aa(:,t+(Balance-1)*3+1)=0;
aa(:,t+(Balance-1)*3+2)=0;
aa(:,t+(Balance-1)*3+3)=0;
aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);

% aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);

%%
t=Busnum*3;
aa(t+(Balance-1)*3+1,:)=0;
aa(t+(Balance-1)*3+2,:)=0;
aa(t+(Balance-1)*3+3,:)=0;
aa(:,t+(Balance-1)*3+1)=0;
aa(:,t+(Balance-1)*3+2)=0;
aa(:,t+(Balance-1)*3+3)=0;
aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(C,1)+2*Busnum*3,size(C,1)+2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+1,t+(Balance-1)*3+1,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
% aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),2*Busnum*3,2*Busnum*3);
%%
dxdy=aa\yy;
%% KLU
%spy(aa)
%dxdy = klu(aa,'\',full(yy));
%%
dX=dxdy(1:Busnum*3*2);
dY=dxdy(Busnum*3*2+1:end);
XX=[
    dX;
    dY;    
];
end