function [z,l,u,w,x,y,PG,QR,Pg,Qg,e,f,T,capK,G,B]=solve(HSB,const,mu,r,nodeNum,balNode,pgNum,pvNum,xNum,transNum,capNum,capI,capBi,capK,...
    pgNode,pvNode,m,x,y,z,l,u,w,Lz,Lw,LZ,UW,dgxT,L1,U1,G,B,transI,transJ,T,transG,transB,transG1,transB1,PG,QR,Pg,Qg)
%程序功能：预计算过程，求解修正方程
%编写时间：2010年11月

%% 求解方程，求出各变量的修正量
dX=-HSB\const;
dx=dX(1:xNum,1);
dy=dX(xNum+1:xNum+m,1);
dl=dgxT*dx+Lz;
du=-dgxT*dx-Lw;
Ll=(l.*z-mu*ones(r,1));
Lu=(u.*w+mu*ones(r,1));
diagZ=sparse(1:r,1:r,z,r,r);
diagW=sparse(1:r,1:r,w,r,r);
dz=-LZ*dgxT*dx-L1*(diagZ*Lz+Ll);
dw=UW*dgxT*dx+U1*(diagW*Lw-Lu);
%% 求最大步长
minusdz=find(dz<0);
dzmin=min(-z(minusdz)./dz(minusdz));
minusdl=find(dl<0);
dlmin=min(-l(minusdl)./dl(minusdl));
minusdw=find(dw>0);
dwmin=min(-w(minusdw)./dw(minusdw));
minusdu=find(du<0);
dumin=min(-u(minusdu)./du(minusdu));

ap=min(dlmin,dumin);
ad=min(dzmin,dwmin);
stepp=0.9995*min(ap,1);
stepd=0.9995*min(ad,1);
%% 更新变量
x=x+stepp*dx;
l=l+stepp*dl;
u=u+stepp*du;
y=y+stepd*dy;
z=z+stepd*dz;
w=w+stepd*dw;
%% 给控制变量、状态变量重新赋值
Tlast=T;                                                                    %保留上一次计算的变比，用于修正节点导纳矩阵
capKlast=capK;
T=x(1:transNum,1);
PG(pgNode,1)=x(transNum+1:transNum+pgNum,1);
QR(pvNode,1)=x(transNum+pgNum+1:transNum+pgNum+pvNum,1);
capK=x(transNum+pgNum+pvNum+1:transNum+pgNum+pvNum+capNum,1);
e=x(transNum+pgNum+pvNum+capNum+1:transNum+pgNum+pvNum+capNum+nodeNum,1);
f=x(transNum+pgNum+pvNum+capNum+nodeNum+1:xNum,1);
Pg(pgNode)=PG(pgNode);
Qg(pvNode)=QR(pvNode);
%% 变比变化后需修正节点导纳矩阵
transG=transG1;
transB=transB1;
G=G-sparse(transI,transI,(T.*T-Tlast.*Tlast).*transG,nodeNum,nodeNum)...
    +sparse(transI,transJ,transG.*(T-Tlast),nodeNum,nodeNum)...
    +sparse(transJ,transI,transG.*(T-Tlast),nodeNum,nodeNum);
B=B-sparse(transI,transI,(T.*T-Tlast.*Tlast).*transB,nodeNum,nodeNum)...
    +sparse(transI,transJ,transB.*(T-Tlast),nodeNum,nodeNum)...
    +sparse(transJ,transI,transB.*(T-Tlast),nodeNum,nodeNum)...
    +sparse(capI,capI,capK-capKlast,nodeNum,nodeNum);
% GTii=sub2ind(size(G),transI,transI);                                        %GTij为G的第i行第j列在G中所对应的位置
% GTij=sub2ind(size(G),transI,transJ);
% GTji=sub2ind(size(G),transJ,transI);
% BTii=sub2ind(size(B),transI,transI);
% BTij=sub2ind(size(B),transI,transJ);
% BTji=sub2ind(size(B),transJ,transI);
% for i=1:length(GTii)
%     G(GTii(i))=G(GTii(i))-(T(i).*T(i)-Tlast(i).*Tlast(i)).*transG(i);
% end
% for i=1:length(GTij)
%     G(GTij(i))=G(GTij(i))+transG(i).*(T(i)-Tlast(i));
% end
% for i=1:length(GTji)
%     G(GTji(i))=G(GTji(i))+transG(i).*(T(i)-Tlast(i));
% end
% % G(GTij)=G(GTij)+transG.*(T-Tlast);
% % G(GTji)=G(GTij);
% for i=1:length(BTii)
%     B(BTii(i))=B(BTii(i))-(T(i).*T(i)-Tlast(i).*Tlast(i)).*transB(i);
% end
% for i=1:length(BTij)
%     B(BTij(i))=B(BTij(i))+transB(i).*(T(i)-Tlast(i));
% end
% for i=1:length(BTji)
%     B(BTji(i))=B(BTji(i))+transB(i).*(T(i)-Tlast(i));
% end
% %  B(BTij)=B(BTij)+transB.*(T-Tlast);
% % B(BTji)=B(BTij);
% deG=G-newG;
% deB=B-newB;
% G(708,708)
% newG(708,708)
% [Grows,Gcols,G1] = find(deG);
% [Brows,Bcols,B1] = find(deB);
% Bcapii=sub2ind(size(B),capI,capI);
% B(Bcapii)=B(Bcapii)+(capK-capKlast);%2011年10月5日修改，去掉*capBi.
% aa=0;
% for i=1:100000
%     aa=aa+1/1.1;
% end
% bb=0;
% bb=bb+sparse(ones(100000,1),ones(100000,1),1/1.1*ones(100000,1),1,1);