threephasese-huber/SolveIt.m

function XX=SolveIt(ContrlCount,Balance,Busnum,Loadi,deltF,PD0,PD,QD0,QD,mVolt,Volt,Mat_H,wVolt,wPD,wQD,dH_dx)
% dH_dx
dP_x=-dH_dx(Loadi,:);
dQ_x=-dH_dx(Loadi+Busnum*3,:);
dV_x=[sparse(1:Busnum*3,1:Busnum*3,1,Busnum*3,Busnum*3),sparse(Busnum*3,Busnum*3)];
dV_x=dV_x(Loadi,:);
H=[dP_x;dQ_x;dV_x];
c=dH_dx(setdiff(1:Busnum*3,Loadi),:);
c=[c;dH_dx(setdiff(1:Busnum*3,Loadi)+Busnum*3,:)];
wVolt=wVolt(Loadi);
aa=[H'*diag([wPD;wQD;wVolt])*H c'
    c zeros(size(c,1),size(c',2))];
% wVolt(find(wVolt==0))=Inf;
% diag([wPD.^(-1);wQD.^(-1);wVolt.^(-1)])
aa=[
    diag([wPD.^(-1);wQD.^(-1);wVolt.^(-1)]) H zeros(size(H,1),size(c',2))
    H' zeros(size(H',1),size(H,2)) c'
    zeros(size(c,1),size(H',2)) c zeros(size(c,1),size(c',2))
    ];
% aa=H'*eye(length(H))*H;
PD=-Mat_H(Loadi,:);
QD=-Mat_H(Loadi+Busnum*3,:);
deltZ=[(PD0-PD);
    (QD0-QD);
    (mVolt(Loadi)-Volt(Loadi));
    ];
yy=[
   deltZ;
   zeros(Busnum*6,1);
    -Mat_H(setdiff(1:Busnum*3,Loadi));
    -Mat_H(setdiff(1:Busnum*3,Loadi)+Busnum*3);
];
% yy=[
%     H'*eye(length(H))*deltZ;
% ];
%% ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD> 
t=size(H,1);
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(H,1)+size(c,1)+2*Busnum*3,size(H,1)+size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(H,1)+size(c,1)+2*Busnum*3,size(H,1)+size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(H,1)+size(c,1)+2*Busnum*3,size(H,1)+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=size(H,1)+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(H,1)+size(c,1)+2*Busnum*3,size(H,1)+size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+2,t+(Balance-1)*3+2,ones(length(Balance),1),size(H,1)+size(c,1)+2*Busnum*3,size(H,1)+size(c,1)+2*Busnum*3);
aa=aa+sparse(t+(Balance-1)*3+3,t+(Balance-1)*3+3,ones(length(Balance),1),size(H,1)+size(c,1)+2*Busnum*3,size(H,1)+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