facat
/
ncp_sigmoid
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
ncp_sigmoid/NCP函数法/linesearch.m

188 lines
8.6 KiB
Mathematica
Raw Permalink Normal View History

备份 Signed-off-by: dmy@lab <dmy@lab.lab>
2015-01-29 21:32:27 +08:00
function [mk,meritF,stepp,stepd]=linesearch(r,nodeNum,balNode,pgNum,pvNum,pgvNum,xNum,transNum,pgNode,pvNode,m,x,y,z,l,u,w,G,B,transI,transJ,transG,...
transB,transG1,transB1,Tn,capKn,capNum,capI,Tn1,capKn1,PG,QR,Pg,Qg,dx,dl,du,dy,dz,dw,opfGoal,Gama,Arlfa,meritF,k2,mu,T,capK,lineNum,lineI,...
lineJ,Pd,Qd,gxx,gsx,mut,mub,stepp,stepd)
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ܣ<EFBFBD><EFBFBD>⻬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>̣<EFBFBD><EFBFBD>󲽳<EFBFBD>
%<EFBFBD><EFBFBD>дʱ<EFBFBD>䣺2010<EFBFBD><EFBFBD>11<EFBFBD><EFBFBD>
mk=0;
kk=k2;
isOK=0;
while isOK==0&&mk<10
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
xk=x+Arlfa.^mk*dx;
lk=l+Arlfa.^mk*dl;
uk=u+Arlfa.^mk*du;
yk=y+Arlfa.^mk*dy;
zk=z+Arlfa.^mk*dz;
wk=w+Arlfa.^mk*dw;
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʱ<EFBFBD><EFBFBD><EFBFBD>Ʊ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>״̬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
Tlastk=T;
capKlastk=capK;
T=xk(1:transNum,1);
capK=xk(transNum+pgvNum+1:transNum+pgvNum+capNum);
PG(pgNode,1)=xk(transNum+1:transNum+pgNum,1);
QR(pvNode,1)=xk(transNum+pgNum+1:transNum+pgNum+pvNum,1);
e=xk(transNum+pgNum+pvNum+capNum+1:transNum+pgNum+pvNum+capNum+nodeNum,1);
f=xk(transNum+pgNum+pvNum+capNum+nodeNum+1:xNum,1);
%<EFBFBD><EFBFBD><EFBFBD>ȱ仯<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ㵼<EFBFBD>ɾ<EFBFBD><EFBFBD><EFBFBD>
G=G-sparse(transI,transI,(T.*T-Tlastk.*Tlastk).*transG1,nodeNum,nodeNum)...
+sparse(transI,transJ,transG1.*(T-Tlastk),nodeNum,nodeNum)...
+sparse(transJ,transI,transG1.*(T-Tlastk),nodeNum,nodeNum);
B=B-sparse(transI,transI,(T.*T-Tlastk.*Tlastk).*transB1,nodeNum,nodeNum)...
+sparse(transI,transJ,transB1.*(T-Tlastk),nodeNum,nodeNum)...
+sparse(transJ,transI,transB1.*(T-Tlastk),nodeNum,nodeNum)...
+sparse(capI,capI,capK-capKlastk,nodeNum,nodeNum);
%% <EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD>
Ezu=exp(-(T-Tn)./mut);
Ezu_1=exp(-(Tn1-T)./mut);
dT_Zt=1;
dNCPT=-mut.*log(Ezu+Ezu_1); %<EFBFBD><EFBFBD><EFBFBD>ȶ<EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ۺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
dNCPT_dt=(Ezu-Ezu_1)./(Ezu+Ezu_1); %<EFBFBD><EFBFBD><EFBFBD>ۺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>T<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD>
transGmut=transG.*dT_Zt;
transBmut=transB.*dT_Zt;
Ezbu=exp((capKn-capK)./mub);
Ezbu_1=exp((capK-capKn1)./mub);
dNCPcap=-mub.*log(Ezbu+Ezbu_1); %<EFBFBD><EFBFBD><EFBFBD>ݵ翹<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ۺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
dNCP_dcap=(Ezbu-Ezbu_1)./(Ezbu+Ezbu_1); %<EFBFBD><EFBFBD><EFBFBD>ۺ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>capK<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD>
ef1=e(transI).*e(transJ)+f(transI).*f(transJ);
ef2=e(transI).*f(transJ)-e(transJ).*f(transI);
ef3=e(transI).*e(transI)+f(transI).*f(transI);
ef1G=ef1.*transGmut;
ef1B=ef1.*transBmut;
ef2G=ef2.*transGmut;
ef2B=ef2.*transBmut;
dhPTi=-ef1G+ef2B+2*ef3.*transGmut.*T;
dhPTj=-ef1G-ef2B;
dhPT=sparse(1:transNum,transI,dhPTi,transNum,nodeNum)+sparse(1:transNum,transJ,dhPTj,transNum,nodeNum);%<EFBFBD>й<EFBFBD>ƽ<EFBFBD>ⷽ<EFBFBD><EFBFBD>P<EFBFBD>Ա<EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD>
dhQTi=ef1B+ef2G-2*ef3.*transBmut.*T;
dhQTj=ef1B-ef2G;
dhQT=sparse(1:transNum,transI,dhQTi,transNum,nodeNum)+sparse(1:transNum,transJ,dhQTj,transNum,nodeNum);%<EFBFBD>޹<EFBFBD>ƽ<EFBFBD>ⷽ<EFBFBD><EFBFBD>Q<EFBFBD>Ա<EFBFBD>ѹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD>
dhNCP_T=sparse(1:transNum,1:transNum,dNCPT_dt,transNum,transNum+capNum);
dhP=sparse(1:pgNum,pgNode,ones(pgNum,1),pgNum,m);
dhQ1=sparse(1:pvNum,pvNode,ones(pvNum,1),pvNum,m-nodeNum);
dhQ=[sparse(pvNum,nodeNum) dhQ1];
dhcapK1=sparse(1:capNum,capI,(e(capI).^2+f(capI).^2),capNum,nodeNum); %2011<EFBFBD><EFBFBD>10<EFBFBD><EFBFBD>5<EFBFBD><EFBFBD><EFBFBD>޸ģ<EFBFBD>ȥ<EFBFBD><EFBFBD>.*capBi_dBci
dhNCP_capK=sparse(1:capNum,transNum+1:transNum+capNum,dNCP_dcap,capNum,transNum+capNum);
dhcapK=[sparse(capNum,nodeNum) dhcapK1 dhNCP_capK];
%<EFBFBD>γɹ<EFBFBD><EFBFBD>ʷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD>
H1=G*e-B*f;
N1=G*f+B*e;
diagE=sparse(1:nodeNum,1:nodeNum,e,nodeNum,nodeNum);
diagF=sparse(1:nodeNum,1:nodeNum,f,nodeNum,nodeNum);
deG=diagE*G;
dfB=diagF*B;
deB=diagE*B;
dfG=diagF*G;
dH1=sparse(1:nodeNum,1:nodeNum,H1,nodeNum,nodeNum);
dN1=sparse(1:nodeNum,1:nodeNum,N1,nodeNum,nodeNum);
Hij=-dH1-deG-dfB;
Nij=-dN1+deB-dfG;
Jij=dN1-dfG+deB;
Lij=-dH1+dfB+deG;
dhJacb=[Hij Nij;Jij Lij]'; %<EFBFBD>ϲ<EFBFBD><EFBFBD>γɹ<EFBFBD><EFBFBD>ʷ<EFBFBD><EFBFBD>̵<EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD>
dhT=[dhPT dhQT dhNCP_T]; %<EFBFBD>ϲ<EFBFBD><EFBFBD>γɵ<EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD>h(x)
dhx=[dhT;dhP;dhQ;dhcapK;dhJacb sparse(2*nodeNum,transNum+capNum)];
%% <EFBFBD>γɲ<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD>ſɱȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>g(x)
gtransNum=0; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>ά<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>NCP<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
gcapNum=0;
dgT1=sparse(transNum,r-lineNum);
dglT=(e(transI).^2+f(transI).^2-e(transI).*e(transJ)-f(transI).*f(transJ)).*transGmut+(e(transI).*f(transJ)-e(transJ).*f(transI)).*transBmut;
dgT2=sparse(1:transNum,lineNum-transNum+1:lineNum,dglT,transNum,lineNum);
dgT=[dgT1 dgT2]; %<EFBFBD><EFBFBD><EFBFBD>в<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>T<EFBFBD>󵼵<EFBFBD>ģ<EFBFBD><EFBFBD>
dgP=sparse(1:pgNum,gtransNum+1:gtransNum+pgNum,ones(pgNum,1),pgNum,r); %<EFBFBD><EFBFBD><EFBFBD>в<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>P<EFBFBD>󵼵<EFBFBD>ģ<EFBFBD><EFBFBD>
dgQ=sparse(1:pvNum,gtransNum+pgNum+1:gtransNum+pgvNum,ones(pvNum,1),pvNum,r);%<EFBFBD><EFBFBD><EFBFBD>в<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Q<EFBFBD>󵼵<EFBFBD>ģ<EFBFBD><EFBFBD>
dgcapK=sparse(capNum,r); %<EFBFBD><EFBFBD><EFBFBD>в<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>capK<EFBFBD>󵼵<EFBFBD>ģ<EFBFBD><EFBFBD>
dgve=2*diagE;
dgvf=2*diagF;
sizeG=size(G);
GIJ=G(sub2ind(sizeG,lineI,lineJ));
BIJ=B(sub2ind(sizeG,lineI,lineJ));
dgleI=(2*e(lineI)-e(lineJ)).*GIJ+f(lineJ).*BIJ;
dgleJ=-e(lineI).*GIJ-f(lineI).*BIJ;
dglfI=(2*f(lineI)-f(lineJ)).*GIJ-e(lineJ).*BIJ;
dglfJ=-f(lineI).*GIJ+e(lineI).*BIJ;
dgle=sparse(lineI,(1:lineNum),dgleI,nodeNum,lineNum)+sparse(lineJ,(1:lineNum),dgleJ,nodeNum,lineNum);
dglf=sparse(lineI,(1:lineNum),dglfI,nodeNum,lineNum)+sparse(lineJ,(1:lineNum),dglfJ,nodeNum,lineNum);
dg0=sparse(2*nodeNum,gtransNum+pgvNum+gcapNum);
dgvl=[dgve dgle;dgvf dglf];
dgef=[dg0 dgvl]; %<EFBFBD><EFBFBD><EFBFBD>в<EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD>Ե<EFBFBD>ѹʵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>鲿<EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD>
dgx=[dgT;dgP;dgQ;dgcapK;dgef]; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD>ſ˱Ⱦ<EFBFBD><EFBFBD><EFBFBD>
%% <EFBFBD>γɶԽ<EFBFBD><EFBFBD><EFBFBD>
Z=sparse(1:r,1:r,zk,r,r);
W=sparse(1:r,1:r,wk,r,r);
L1=sparse(1:r,1:r,1./lk,r,r);
U1=sparse(1:r,1:r,1./uk,r,r);
LZ=L1*Z;
UW=U1*W;
Ll=(lk.*zk-mu*ones(r,1));
Lu=(uk.*wk+mu*ones(r,1));
%% <EFBFBD>γ<EFBFBD>Ŀ<EFBFBD>꺯<EFBFBD><EFBFBD>һ<EFBFBD><EFBFBD>ƫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>f(.)
if(opfGoal==2)
dPG=ones(pgNum,1);
dfx=sparse(transNum+1:transNum+pgNum,ones(pgNum,1),dPG,xNum,1);
else if(opfGoal==1)
dPG=2*a(pgNode).*PG(pgNode)+b(pgNode);
dfx=sparse(transNum+1:transNum+pgNum,ones(pgNum,1),dPG,xNum,1);
else if(opfGoal==3)
dfT=-1*transG.*dT_Zt.*((f(transI)-f(transJ)).^2+(e(transI)-e(transJ)).^2);
dfPQ=zeros(pgvNum,1);
dfcapK=zeros(capNum,1);
EE=e*ones(1,nodeNum);
FF=f*ones(1,nodeNum);
dfe=2*sum(G.*(EE-EE'));
dff=2*sum(G.*(FF-FF'));
dfx=[dfT;dfPQ;dfcapK;dfe';dff'];
else if(opfGoal==4)
dfT=1*transB.*dT_Zt.*((f(transI)-f(transJ)).^2+(e(transI)-e(transJ)).^2);
dfPQ=zeros(pgvNum,1);
dfcapK=zeros(capNum,1);
EE=e*ones(1,nodeNum);
FF=f*ones(1,nodeNum);
dfe=-2*sum(B.*(EE-EE'));
dff=-2*sum(B.*(FF-FF'));
dfx=[dfT;dfPQ;dfcapK;dfe';dff'];
end
end
end
end
%% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
gxv=[PG(pgNode);QR(pvNode);e.^2+f.^2];
gxl=(e(lineI).*e(lineI)+f(lineI).*f(lineI)-e(lineI).*e(lineJ)-f(lineI).*f(lineJ)).*GIJ+(e(lineI).*f(lineJ)-e(lineJ).*f(lineI)).*BIJ;%<EFBFBD><EFBFBD>·<EFBFBD><EFBFBD><EFBFBD>ʵ<EFBFBD>ǰֵ
gx=[gxv;gxl]; %Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ǰֵ
Lz=gx-lk-gxx;
Lw=gx+uk-gsx;
muUL=mu*(U1-L1);
BB1=UW*Lw-LZ*Lz-muUL*ones(r,1);
BB=dhx*yk-dfx+dgx*BB1;
BB(transNum+pgvNum+capNum+nodeNum+balNode,1)=0; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>ѹ<EFBFBD>鲿<EFBFBD><EFBFBD><EFBFBD>䣬ʼ<EFBFBD><EFBFBD>Ϊ0<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD>ڵ㴦<EFBFBD><EFBFBD>0
BB(transNum+pgvNum+capNum+balNode,1)=0;
%% <EFBFBD><EFBFBD>ʽԼ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
Pg(pgNode)=PG(pgNode);
Qg(pvNode)=QR(pvNode);
P=Pg-Pd;
Q=Qg-Qd;
dP=P-diagE*H1-diagF*N1;
dQ=Q-diagF*H1+diagE*N1;
hx=[dP;dQ;dNCPT;dNCPcap]; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>̳<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>h(x)
const=[BB;hx]; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>̳<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
L3=-L1*(Z*Lz+Ll);
L4=U1*(W*Lw-Lu);
L5=Lz;
L6=Lw;
f_NCP=[const;L3;L4;L5;L6];
f_NCPsum=(f_NCP'*f_NCP)/2;
maxfk=max(meritF(kk-min(10,kk)+1:end));
f3=maxfk-Gama*(Arlfa.^mk)*meritF(kk);
if f_NCPsum<=f3
isOK=1; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
else
mk=mk+1; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
end
end
stepp=stepp*Arlfa.^mk;
stepd=stepd*Arlfa.^mk;