2023年优化方法上机大作业

优化措施上机大作业机械工程与材料能源学部能源与动力学院能源与环境工程联络方式:

x0=[0;1］T;%初始值s０＝[-１;1]T;%初始搜索方向c１=0.1;c2=0.5;a＝0;b=ｉｎf;d=１;n＝０;x１=x０＋d*s0；g0=[(x0(2）-x０(1）^2）*ｘ０(1)-２*（1-ｘ0（1));(x０(2）-x０(1）^2)];g１=［(x１(2)-x1(1）＾2）*ｘ1(1)-2＊(１-x1(1));(ｘ1(2)-x1(1)＾2)］;f1=(ｘ1(2)－ｘ１(１)＾2)^２＋(1－x１（１))^2;f０=（x０(2）－x０(1)^2)^２+(1-x0(1）)^２;wｈiｌe((f０-f1<-ｃ１*d*ｇ０'*ｓ0)||(g1'*s0<c2*ｇ0'*ｓ0))if(（f０-f1）<（-ｃ1*d＊ｇ0'*s0))b=ｄ；d=(ｄ+a)/2；x1=x0+d＊s0;g０＝[(ｘ０(2)-x0(1)^２)*x0(１）-２＊(1-x0(1）);（ｘ0（2)－x0（1)^2）];g1＝[（x1（2)-x1(1）＾2)＊x1（1)-2*(１－x1(1));(x1(2)-x1(1）＾2)］;ｆ１=(x1(２)-x1(1）^2）＾２＋(1-x1(1))^２;f0=（x0(2）－ｘ0(1)＾2)^2＋(1-ｘ0(1)）＾２；eｌｓeiｆ（((g1')*s０)<(c2*(g0'）＊s０))a=d;if（２*d<=(d+b）/2)d=2*ｄ;eｌsed=(d＋ｂ)/2;ｅndx1=x0+ｄ*s0;g0＝［(ｘ０(2)－ｘ0(1)^2)*x0(1)-2*（１-x0（1)）;（ｘ0(２)-x０(1)^2)]；ｇ１=[(x1(2)-x1(1）＾2)*x１(１)-2＊（1-ｘ１（1));(ｘ１(２)－x1(1)^2)］;f1=(x1(2)-ｘ1(1)＾２)＾２＋(1-ｘ１（１))^２;f０=(x0（2)－x0（1)^2）^2＋（1-x0(１））^２；ｅnｄeｎｄx1f1=（x1(2)-x1(1)^２)^2＋（1-x1(１))^2x1=-０.0０001.0000ｄ=1．１1０2e-016f1=２functｉonｆ=fｕn（x)%UＮTIＴLＥＤ3Sumｍaryofthisｆｕnｃtiongｏesｈere%Ｄeｔailedeｘｐlanaｔｉoｎgoesheｒef=x（1)^２-2*ｘ(1)＊ｘ（２)+2*x(２）＾２+x(３)^2+x（４）^2－x(2）*x(3)+2*x(1)+3＊ｘ（2）-x(３）;Ｅndfｕnctionｇ＝fun(x）%ＵNTＩTLＥD4Ｓｕmmａrｙｏｆthiｓfunctｉｏngoeｓhere%Ｄetailedexｐlanatｉｏngoeshereg=[2-200;-2４－１0;0-１20;0０02]＊ｘ＋[2；3;-1；０]；ｅnｄx0＝[0;0;0;0];%初始值eps=1.０e-４;％精度ｇ0＝gｆuｎ(ｘ０）;s0=-g０;n=0;symsｄ1;whiｌenorm(g0)>eｐｓifn<３g＝gfuｎ(ｘ0+ｄ1*s０）；d=ｄoubｌｅ（soｌvｅ(s0'*ｇ));x1＝ｘ0＋ｄ*s0;g1＝ｇfuｎ（x1)；ifnoｒｍ(ｇ1)<eｐｓｎ＝n+1;ｘ0=x1；breakelses0=－g1+(norm(ｇ1)^2/norm(g0)＾2)＊s0;x0=x1;ｇ０=g1;endelｓeｉfn=３ｘ0=x1;g0＝gｆun(x0);s0=-g0;ｎ=0;ｅndn=n＋1;endx0ｎfｕn(ｘ0)x0=-４-3-10n＝3anｓ=-８functionf=fun３_1（x)％FUＮ3Summaｒyofthｉｓfｕｎｃtiｏnｇｏeshere%Detaiｌedexｐlａnaｔionｇoeshｅref=ｘ（1)+2*x(2)^2+exp(x(1)^2＋ｘ(2)^2);endfunctioｎｇ=gfｕn３_1(ｘ)%GFUＮ3＿1Sｕmmaｒyｏｆthisfunｃtiｏnｇoeshｅre％Detaｉledeｘｐlａnatｉｏngｏｅsheｒeg＝［1+2*x（1)*eｘp(x(1）＾2+x（２）^2）;4*ｘ（2)+2*x(2)*exｐ（x(１）^２＋x(２)＾2)]；eｎd最速下降法x0＝[0;1］;%初始值eps=１.0ｅ-5；%精度n=0;g0＝gｆuｎ3_1（x0)；sｙｍsd１;whｉlenoｒm(g０)>＝epss0＝-g０;g=ｇｆuｎ3_1（x０+d1*s0);d=douｂle(sｏlve（s0'*g）)；x１=ｘ0＋d*s0;g1=gｆun3＿1(x1);if(norｍ(g1)<ｅps)ｎ=n＋1;x0=x1;bｒｅaｋ；ｅlsｅｘ0＝x1；g0＝gｆun3＿1(x0);ｅnｄn=n+1;eｎdf0=fun3_1(x０）x0nf0=0.7７29ｘ0=－0.41940.０000ｎ=6（2）牛顿法functiong2=hesse(ｘ)%HEＳSEＳｕmmaryｏfthiｓｆuｎcｔiｏngoｅshere％Ｄetailedｅxpｌanationgoeｓherｅg2＝[2*exp(x（1)＾2+x（2）＾2）+4＊x(1)^2*exp(ｘ(1)^2+x（2）＾2)，4＊ｘ(１)*x(2）*exp(x（１)^2+x（2）^2)４*ｘ（1）*x(2）*ｅxp(x(１)^2+x（2)^２)，4+2*exp（x（１)＾２+x(２)^2)＋4*ｘ(2）^2*exp(x(１）^２+x（2)^2)];ｅndx0=[0;１]；％初始值eps=1.0ｅ-5;%精度n＝0;g0=gfuｎ3_１（ｘ0);g20=hessｅ(ｘ０);wｈｉleｎorm(g0)>＝ｅpｓd＝-ｇ２0\g0;x1=ｘ0+ｄ;g1=gfun3_１(ｘ１);if（norｍ（ｇ1）<eps)n=n+1;x0＝x1;break；ｅlsｅx0=x1;ｇ0=ｇfｕn3＿1（x0）;eｎdn=n+１;ｅndf0=ｆun３＿１(x０)x0nf0=0.772９x０=-0.4１9４0.０000n=３５（3)ＢFＧSx0=[0;1];%初始值ｅps=１．0ｅ-５;％精度n=０;g０=ｇfun3_1(x0）;ｓymｓd1;h0＝eye(２);whilenorm(g0)＞=epｓs0＝-h0*g0；g=gｆun３_1(x０+ｄ１*s0)；ｄ=double(soｌvｅ(s０＇*g)）；x1=ｘ0＋ｄ*s0;g1=gfun3_１(x1);iｆ（nｏｒm(g1)<eps)n=n+1;x０=x1;ｂreak;elseh0=h0-(ｈ0*(ｇ１-g0)＊(g１-g０)'*h0）／（(g1-g0）'*h0*(g1-g0))+...((ｘ1-x0)*(ｘ1-x０)＇)/(（ｘ1-x０)'*(g1-g0))+．．.((g1－g0）'*h0*（g1－g0）)*((x1-x0)＊(ｘ1-x0)'）;x０＝x1;g０=gｆuｎ3_1(ｘ0);enｄn＝n＋１;ｅndf０=fun3_1(x0)x０nf0=0.77２9ｘ0=-0.41940．0000n=4function[x,lamk，exiｔflaｇ，ouｔpuｔ]＝ｑpacｔ（H,c,Ae,bｅ,Ａi，bi,x0）ｅｐsilon=1.0ｅ-９;ｅrｒ=1.0e-6;k=0；x=x0;n=ｌenｇｔh(ｘ）;ｋmax=1.0e3;nｅ=lｅｎgｔh(ｂｅ);ni=ｌenｇth(bi);lamｋ=zeroｓ(ne＋ni,1);ｉndｅｘ=onｅｓ(ｎｉ,1);fｏｒ(i＝1:ni)ｉf(Ａｉ(i,:)*ｘ>bi(i)+epｓilon）,index(i）=0;endendwｈile(ｋ<＝ｋmａx)Aｅｅ=[];if(ne>0）,Ａee=Ae;endfor(j=1：ni）if(index(j)>0)，Aee=[Aｅｅ；Ai(j,：)];endendｇk=H*ｘ+c;[ｍ１,n1]=ｓｉｚｅ(Aｅe);[ｄｋ,lａｍk］=qsｕbp(H,gｋ,Aee,zeros(m1,1)）;if(norｍ（dk）<=ｅrr)y=０.0;iｆ(lｅngｔh(laｍk)>ne)[y,jk]=min（lamk(ne＋1:leｎgth(laｍｋ)));endｉf(y>=0）ｅxitflag＝0;elseｅxｉｔfｌaｇ＝1;for(i=1:ni)ｉｆ(indｅｘ(ｉ)&(ｎe＋sum（iｎdex(1：ｉ)))＝=jk)iｎｄex(i)=０;brｅａk;enｄenｄenｄk=k+1;elseexiｔfｌag=1；ａlpha＝1.0；tm=1.０;fｏｒ（ｉ＝1:ni)ｉf((iｎdex(i）＝=０)&(Ai(ｉ，：)＊dk＜0)）ｔｍ1=（bｉ(ｉ）-Ai(i,:)＊x)/(Ai(i，:）*dk）;iｆ(ｔｍ1<tm）tm=tm１;ｔｉ=i;endｅndendaｌpha=mｉn(ａlpha,tm);x=ｘ＋ａｌpｈａ＊dk;iｆ(ｔm<1)，inｄex(ti)=1;enｄｅｎdif(exiｔｆlａg=＝０)，ｂreａk;ｅnｄk=k+1;ｅndｏuｔｐut.fval=０.5＊x＇*H*x+c'*x;ouｔput.itｅr＝k;ｆunction［x,ｌａmbda]=qsubｐ（Ｈ，ｃ,Ae,bｅ)ｇinvＨ=pｉnv（H）；[ｍ，ｎ]=size（Ａｅ);if（m>0)ｒb=Ae*ｇｉnvH＊c＋ｂｅ；lambda＝pinｖ(Ae*ｇinｖH＊Aｅ'）*rb；x＝ginvＨ＊(Ａe'*lａmｂｄa-ｃ);elｓｅx=-ginvH*c;laｍbda=ｚｅros(m,1);ｅndｃａllｑｐａct.ｍ文献functionｃallqpａctＨ＝[２0;02];ｃ=［-2－5]';Ａe＝[］;ｂe=［];Aｉ=[1-2;-1－2;-１2;10;０1];bi＝［－2-６-2００]＇；x0＝[00］＇;［x，lａmbｄa,ｅｘitｆlag,ｏutpｕt]＝qpacｔ(H,ｃ，Ae,be,Aｉ,bi,x0)运行成果:ｃａllqpaｃtx＝１.４0０0１.７000ｌambda＝0．80００exiｔｆlaｇ＝0ｏutpｕt=fvａl：-６.4500iｔｅｒ:7Functｉｏn[x，mｕ，ｌaｍbda,ｏutput]=ｍuｌtｐhｒ(ｆuｎ,hf,ｇf,dfun，ｄhf,dｇｆ,x0)functｉonpsｉ＝mｐｓｉ(ｘ，fｕn,hｆ,gf,dｆun，dｈf,dgf，ｍu，ｌamｂda,sigma）functiｏnｈe=ｈ1(x）he＝-ｘ（1)^2-x(2)^2+２5.０;fｕnctｉｏｎf=f1（ｘ）f=4*x(1)-x（2)^2－１2;ｆunctiongi=ｇ1(x)ｇｉ=10＊x(1)-x(1)^2+１０＊x（2)－ｘ(2）^2－3４;ｆunctiｏnｄhe=dｈ1(x)dhe=[-1*x（１)，-1*ｘ(２)]';ｆｕnctionｄgi=dg1(x)ｄｇi＝[1０-２*ｘ(１）,10－2＊x(2)］';ｆｕnｃｔionｇ=dｆ1(x)g=[4，-2.0*x(2）]';fｕncｔion[x,ｖal,ｋ］=bfｇs(ｆuｎ,gfun,x0,vａrａrgin）maxｋ=５０0;ｓigma＝2.0;eta＝２.0;thｅta=0．8；k=0；iｎk＝０;epsiｌｏｎ=1e-5;x＝x0；ｈe＝feval(hf,x);gｉ=fｅval(gｆ,ｘ);n＝lengｔｈ(x);ｌ=lｅngｔｈ(he)；m=length(gi);mu=0.1＊ｏnes(l,１);lambda=0.１*ones(ｍ，１);btaｋ=10;ｂｔaｏｌd=10；ｗｈｉｌe(btak>eｐｓilon&k＜maxk)［ｘ,ｉval,iｋ]=bｆgs('mpsi＇,＇dｍｐsi＇,x０,fun,hf，ｇf，dfun，ｄhｆ,dgf,mu,lambｄa,siｇma);ｉｎk=inｋ＋ik;he=fｅval(hf,x);gｉ＝fｅval(gf，ｘ);ｂtak=0.0;ｆoｒ(i=1：l）,bｔak＝btak＋ｈｅ(i）^2;ｅndｆori＝1:mteｍｐ=ｍｉｎ(gｉ（i)，lａmｂda(i)/ｓｉgma);bｔak=ｂtak+teｍp^2;endbｔak=sqｒｔ(btak）;iｆｂtak>epsｉlｏnif(k＞=２&ｂtak＞ｔheｔa*ｂtaolｄ）siｇｍａ=ｅta*ｓigma;endｆor(i=1：ｌ),mu（i）=mｕ(i)－sｉgmａ*he(ｉ);eｎｄfor(i＝1:m)lamｂdａ(i)＝mａx（０．0,lambda(i)－ｓｉgｍａ＊gｉ(i));endeｎｄk=k+1;btaolｄ＝btak;ｘ0＝x;endｆ=ｆｅval(ｆun，x）;ｏutput．fｖaｌ＝f;output.itｅr＝k；outpｕt.ｉnｎer_itｅr＝ｉｎｋ；oｕtput.bｔa=btak；ｆ=fevaｌ(fun，x);he=fevａl（ｈｆ,x);gi=fｅval(gf，ｘ）；l=lｅngtｈ(ｈｅ);m＝lenｇth（gi);psi=f;s1＝0．0;ｆｏｒ(i=1:l）ｐsｉ＝psｉ-ｈe(i）*ｍｕ(i)；s1=s1+hｅ(ｉ)＾2;endpｓi＝psi＋0.5*sigｍa*s1；s2=0．０；for（i=1:m）s3=max(0.０，ｌamｂｄa(i）－ｓiｇｍａ*gi(i));s2=ｓ2+s3^2-lamｂda(ｉ)＾2；endpsｉ=psi+s2/（2．0*sigｍａ）;maxk=5０0；rho=0.55;sigｍa1=０.4;epsilｏn1=1e-5；k=0;ｎ＝lengtｈ(x０);Bk=eye(n）；whｉｌｅ(ｋ<maxk）ｇk＝fｅｖaｌ(gｆｕｎ,x0，vaｒargin{:});iｆ(norm(gk)<eｐsilon1),ｂreａk；eｎddk=-Bk\ｇｋ;m=０;mk＝0;whｉｌe(m＜20）newf＝fｅvaｌ(fｕn,x０+ｒho＾ｍ*dｋ,ｖarargin{：}）;ｏldf＝feval(ｆun，ｘ0,varargｉn{:});if(nｅｗf＜oldｆ+sigmａ1*rho^ｍ*gk'*dｋ)ｍk=m;bｒeａk；enｄm＝m+1;enｄx=x0+rho^mｋ*dk;sｋ＝x-x０;ｙk=feval（ｇｆｕn,x,varaｒgin｛:｝）-gk;iｆ（yk＇＊sk>０）Bk＝Bk-(Bｋ*sk*sk'＊Bk)/(sk'*Bｋ*ｓk）＋(yk＊yk')/(yk'*ｓｋ)；ｅndk=k+1；x0=ｘ;endval=feval(ｆｕn，ｘ0,ｖarａrgiｎ{:});functiｏｎdpsi＝dmpｓｉ(x，ｆun，ｈf,gｆ，dｆun,dhf，ｄgf，mu，lambda,sｉgma）dpsi=fｅｖal（dfｕn，ｘ);he＝ｆeval(hｆ,x);ｇi=feｖal(gf,x）；ｄhe=feval(dhf，ｘ);dgi=fｅval(dgｆ，ｘ)；l=lenｇtｈ(ｈｅ);ｍ=lengｔｈ(gi)；foｒ(ｉ=1:l)ｄpｓi=dpsi＋(sｉgma＊hｅ（ｉ)-mu(ｉ))＊dｈe(:,i）；ｅｎdｆor(i=1:ｍ)dｐｓi=ｄｐｓｉ+(siｇma*gi(ｉ)-lａmbda(i）)*ｄｇi(:，i);end>>x0=［1,1]＇;[x,ｍｕ，lambdａ，ouｔｐｕｔ］＝mulｔphｒ('f1','h１','g1','dｆ1＇,＇dh1'，＇ｄg1＇,x０)x=1.０0１３4.８987mｕ=2．０3１2lamｂda=0.7545oｕtｐut＝ｆval:－３1．9９２3iｔｅr:5inneｒ＿ｉter:58bta:4.３１87e-07fuｎｃtioｎ[x，ｍu,ｌam,ｖaｌ,k]=sqpｍ(x0,mｕ0，lａm0)maxｋ=10０；n＝ｌeｎgｔｈ(x0);l=ｌengtｈ(mｕ0)；m=leｎgth（laｍ0);rｈｏ=０．5;eta＝0.1;B0=eyｅ(ｎ);ｘ=x0;mu=mｕ０；ｌａｍ＝laｍ0;Bk=B0;siｇma＝0.8;epｓilｏｎ1=1e-６;eｐsilｏn2=１e－5;[hk,ｇk]＝cons(x);dfk=dｆ1(ｘ);[Ａe,Ａi]=dconｓ(ｘ）;Ak=[Ae;Ai]；k＝0；while（k＜ｍaxk)[dk,ｍu，lam]=qpｓｕbp（dfk，Bk，Aｅ,hk,Ai，gｋ);mp１＝ｎｏrm(hk，1）＋ｎｏrm(max（－gk,0)，１）；if(noｒm（ｄk,1)<epsilon1)&(ｍp1<epｓilｏｎ２)ｂrｅaｋ;enddｅta=0.0５;tau=max(ｎoｒm(mｕ,iｎf),ｎorm(ｌam,inf));ｉf（ｓigｍa＊(tau＋deta)<１)sｉgmａ＝sｉgmａ；elsesigma=1.0／(ｔaｕ+2＊deta);ｅｎdim＝０；whｉle(ｉm＜=２0)ｉf(phi1（x+rｈｏ^ｉm*ｄk,sigｍa）-phi1(x,sigma)<ｅta*ｒhｏ＾ｉm*dphi1(ｘ,sｉgｍa,ｄk）)mk=im；bｒeａk;enｄim＝im+1；if(im==20),ｍk=10;endeｎdalpha＝rho^ｍk；x１=x+alpha＊ｄk;[hk,gk]=cｏns(ｘ１);ｄｆk=df１（x１);[Ae,Aｉ]=dcoｎs(x1)；Ak=［Ae；Ai]；lamu=pｉnv(Ak)＇*dfk;iｆ(l>０&ｍ>0)mu=lａmu(1:l);laｍ＝lamu(l+1:l+m);ｅｎdif（l==0)，mｕ=[]；lam＝lamu；eｎdif（m＝=０)，mu=laｍｕ;lａm=[]；enｄsｋ=alpha*dk;yk＝dlａx(ｘ1,mｕ，ｌam)-dlax(x，mu，lａｍ);ｉｆ(sｋ'*yk>0.2*ｓk'*Bｋ*sｋ）tｈeｔa=１;elseｔheta=0．8*sk＇*Ｂk*ｓk/(sk＇*Bk＊ｓｋ-sk'*yｋ);enｄzk＝theta*yk+（1-theｔa)*Ｂｋ*sk;Bk＝Bk+ｚk＊zk'／(ｓk＇*zk)-(Bk*ｓｋ)*(Bｋ*sk)'／(sk'*Ｂk*sｋ);x=ｘ1;k=k+1;endvａl=f1(x);fuｎｃtioｎp=pｈi1(x，ｓigmａ）f＝f１(x)；［h,g]=cons(x);ｇn=max(－g,０）;l0=lｅｎgtｈ(ｈ）;m０=leｎgtｈ(g);if(ｌ0==0)，p＝ｆ＋1．０/sigｍa*ｎorｍ（ｇn,1)；endｉｆ(m0==０）,ｐ=f+１.０／sigma*norm(ｈ,1）；enｄiｆ（l0＞０&m0>0)p=f+１.０／sigma＊(noｒm(h,1)+norｍ(gn,1）);endfunctｉoｎｄp=ｄphi1(ｘ,ｓigmａ,d)df＝df1(ｘ）；[h，g]=conｓ(x);gn=max（-ｇ,0);l0=length（h)；m０＝leｎgth(g);iｆ(l0==0),dp=ｄf＇＊ｄ-1.0/siｇｍa*ｎｏrｍ(gｎ，１）;endiｆ(m0==０),dp＝dｆ'*d-１．０／sigｍa*nｏｒm(h,1);enｄif（l0>0&ｍ0>0)dｐ＝df＇＊d-1.0／sigma*（ｎoｒm(h,1)+nｏrｍ(ｇｎ,１));enｄfuncｔｉonl=lａ(x，mｕ,laｍ)f＝f１（x);[ｈ，ｇ］=conｓ(ｘ）;l0=lｅmgth(h);m０=lengｔh(g);if(ｌ0==0),ｌ=ｆ-lam*g；endif(ｍ0==0）,l=f-ｍu'*ｈ;endif(l０>0＆m0>0)l=f-mｕ＇*ｈ-lａｍ＇*ｇ;enｄｆunｃｔiondl=dlａx(x，mｕ,lam)df=dｆ1(ｘ);[Ae,Ai］＝dcｏns（ｘ);[ｍ1,m２]=siｚe(Aｉ);[l1,l２]=size(Ae)；if(l1==０）,dｌ=ｄｆ-Ａｉ'＊lam;endｉf（ｍ1==0),ｄl=df-Ae＇*mｕ;endiｆ（l1>0&m1>0）,dl＝df-Ae＇*mｕ-Ａi'*lam;eｎｄfｕnｃtｉonf=f1(x)f=4*x(１）-x(２)＾2-1２；funｃtｉoｎdf＝ｄｆ1（x)df=［4,-２*x(２)］';ｆuｎctｉon[ｈ,ｇ］=ｃons(x）h＝[２5-x(1)^2-x(2）^2]；ｇ＝[x(1);ｘ(２)]；functiｏｎ[ｄｈ,dg]=dｃoｎｓ(x)dh=[-2*ｘ(1）,－2＊ｘ(２)];dg=[１0;0１];Qｐｓubp.ｍ文献ｆunｃtｉon[d,mｕ,laｍ,vaｌ,ｋ]＝qpsubp(dfk,Bk，Ae,hk，Ai,gｋ）ｎ=length(ｄfk);l=lｅｎgtｈ(hｋ);m＝lenｇth(gk)；gamma=0.０５；eｐｓｉlon=1.0e-6;rhｏ=0.5；siｇｍa=0.2;eｐ0=0.0５;mu０=０.05＊zｅros(l,1);lａm０=0.0５*zeｒoｓ(m,１);d0=oｎes（n，1);u0＝［ep0;ｚeros（ｎ+l+m,1）];ｚ0=[ep0;d0；mｕ0；ｌａm0，]；k=0;ｚ=z０;ｅｐ=ｅp0;d=d0;mu＝mu0;lａｍ＝laｍ0;whilｅ(ｋ<=150)dh=ｄah(ep,d,mu,lａm,dfｋ，Ｂk,Ae,hk,Ai,gk);iｆ（norm(dｈ）<ｅｐsilon)ｂｒeａk；ｅndＡ=JacｏbiＨ（ep,ｄ,ｍu,ｌａm，dｆｋ,Ｂk,Ae,hk,Ａi,gｋ）;b=beｔa(ep,d,mu,laｍ,dｆk,Bk,Aｅ,ｈｋ,Ai,gｋ,gａmma)*u０-dh;dz=A\b；iｆ(ｌ＞0＆m>0)ｄｅ=dz(1);dd=ｄz（2：n+1);ｄｕ=dz(ｎ+2:n+l+１)；dl=ｄz(n＋l+2:ｎ+l＋ｍ＋１);eｎdif(ｌ==０)de=ｄz(1);ｄｄ=dz(2:n＋1);dl＝dz(n+2：n+m＋1);ｅndif(ｍ==0)de＝dz(1)；dd＝dｚ(2：n+１);dｕ=ｄｚ(n+２:n+ｌ+1);endｉ=0;ｗhile(i＜=２０)if(ｌ＞０&ｍ>０)dh１=ｄah(ep+rho^i*ｄe,d＋rho^ｉ*dd,mu＋rho＾ｉ*ｄu,lａm＋rｈo^i*dl,dfｋ，Bk,Ae，ｈk,Ai,gk);endif(l=＝0）dh1＝dah（ｅp＋ｒｈo^ｉ*ｄe,d+ｒho^ｉ*ｄd,mu,ｌam+ｒhｏ＾i＊dl,dｆk,Bk，Ae，hｋ,Ai，ｇk);endif（m==0)ｄh１＝dah(eｐ＋rho^ｉ*dｅ,ｄ+ｒhｏ^ｉ*dｄ，mu+ｒｈo^i*dｕ，lam,ｄfk，Ｂk，Ae,ｈk,Ai，gk);eｎdiｆ(nｏrm(dｈ1)<＝（1－sｉgma*(1-gamｍa*ep0)*rho^i)*norm(dh))mｋ=i;break；ｅｎdi=ｉ+１;ｉf(i=＝20）,mk=10;endｅｎｄalpha=rho^ｍk；iｆ（ｌ>0&ｍ＞0)ep=ep＋alpha*de;d=d+alｐha*ｄd；mu=mu+alpha*ｄu;ｌaｍ＝lam＋aｌｐha*dｌ;ｅndif(l=＝０)eｐ=eｐ＋ａｌｐhａ＊de;d＝d+aｌpha*ｄd；ｌam＝ｌam+alｐha*dl;endif(m=＝0)ｅｐ=ep+alpha*de；d=d+ａｌpha*ｄd;mu=mu+alpha＊ｄu;endk=k+1；endｖal=０.5＊ｄ'*Bk＊d+dｆk＇*d;functionp=ｐhi(eｐ,a,b)p=a+b－sｑrt（a^2＋b^2＋2＊ep＾２）;fｕnctionｄh=dah（ep，d,mｕ，lam，dfk,Bk,Ae,ｈk,Aｉ,gｋ)ｎ=lengtｈ(dfk)；l=leｎgth（hｋ）;ｍ＝ｌenｇth(gk）;ｄｈ＝zerｏs(n+ｌ+m+１,1）;dｈ(1)=ep;iｆ(ｌ＞0&m＞0)ｄh(2:n+１)=Ｂk*d-Ａe'*ｍu-Ai'*ｌａm+dfk;dh(n＋２：n+l+1)=hk＋Ae*d;for（i＝1:m）dｈ(n+l+1+ｉ)=ｐhｉ(ep,ｌam(ｉ)，gｋ(ｉ)＋Ａi(i，:)*d);endｅｎｄif（l==0)dh(２:ｎ+1)=Bｋ*ｄ-Ai'*laｍ+dfｋ；for(i=１:ｍ）dh(n+1+ｉ)=phｉ(ep，lam(i)，gｋ(i)＋Ａｉ（i,:）*ｄ);endｅｎｄif(m＝=0）dh（２：n+1)=Ｂk*d-Ae'*mu＋dfk;dh(n+2：n＋ｌ+１)＝ｈｋ+Ae*ｄ;ｅnddh＝dh(:）;ｆuncｔiｏｎｂet=bｅta(ep,d，mu,lａm,dfk,Bｋ,Ae,hk,Ａi,gｋ,ｇamｍａ)dh＝ｄａh（ep，d,mｕ，lａm，dfｋ,Bｋ,Ae，hk,Ai，gk）;ｂet=gaｍmａ*nｏｒm(ｄh)＊min(1,ｎoｒm(ｄh）);