整理大连理工大学优化作业程序.docx
- 文档编号:13818843
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:17
- 大小:289.18KB
整理大连理工大学优化作业程序.docx
《整理大连理工大学优化作业程序.docx》由会员分享,可在线阅读,更多相关《整理大连理工大学优化作业程序.docx(17页珍藏版)》请在冰点文库上搜索。
整理大连理工大学优化作业程序
1.1程序(Java)
publicclassWolfe_Powell{
publicstaticdoublegetFx(double[]x){
doublex1=x[0];doublex2=x[1];
doubleFx=100*(x2-x1*x1)*(x2-x1*x1)+(1-x1)*(1-x1);
returnFx;
}
publicstaticdouble[]getDeltFx(double[]x){
doublex1=x[0];doublex2=x[1];
double[]deltFx=newdouble[2];
deltFx[0]=-400*(x2-x1*x1)*x1-2*(1-x1);
deltFx[1]=200*(x2-x1*x1);
returndeltFx;
}
publicstaticdoublegetDeltFx_Sk(double[]deltFx,double[]Sk){
doublea=0;
for(inti=0;i a=a+deltFx[i]*Sk[i]; } returna; } publicstaticdoublegetL(double[]x,double[]s){ doublex1=x[0];doublex2=x[1]; doublec1=0.1,c2=0.5,a=0,b=1e8,L=1; doubleFx0,Fx1,deltFx1_Sk,deltFx0_Sk,temp,temp2; double[]deltFx0,deltFx1; Fx0=getFx(x); deltFx0=getDeltFx(x); deltFx0_Sk=getDeltFx_Sk(deltFx0,s); temp=c2*getDeltFx_Sk(deltFx0,s); for(inti=0;i<1e8;i++){ temp2=-c1*L*deltFx0_Sk; x[0]=x1+L*s[0]; x[1]=x2+L*s[1]; Fx1=getFx(x); deltFx1=getDeltFx(x); deltFx1_Sk=getDeltFx_Sk(deltFx1,s); if((Fx0-Fx1)>=temp2&&deltFx1_Sk>=temp){ break; } elseif((Fx0-Fx1) b=L; L=(L+a)/2; } elseif(deltFx1_Sk a=L; L=(L+b)/2>=2*L? (2*L): (L+b)/2; } } System.out.println("L="+L); System.out.println("计算次数"+i); returnL; } publicstaticvoidmain(String[]args){ Wolfe_Powelltemp=newWolfe_Powell(); double[]X={-1,1};double[]sk={1,1};temp.getL(X,sk); } } 1.2实验结果 步长L=0.00390625x=[-0.9992,1.0324]计算次数8 2.1程序(Java) publicclassGongE{ publicstaticdoublegetFx(double[]x){ doublex1=x[0]; doublex2=x[1]; doubleFx=x1*x1-2*x1*x2+2*x2*x2+x3*x3-x2*x3+2*x1+3*x2-x3; returnFx; } publicstaticdouble[]getDeltFx(double[]x){ doublex1=x[0]; doublex2=x[1]; double[]deltFx=newdouble[x.length]; deltFx[0]=2*x1-2*x2+2; deltFx[1]=-2*x1+4*x2-x3+3; deltFx[2]=2*x3-x2-1; returndeltFx; } publicstaticdouble[]getX(double[]x){ double[]g0,g1; double[]s0=newdouble[x.length]; double[]s1=newdouble[x.length]; doubleg0_L,g1_L,L,temp; double[]x0=x; intk=0; g0=getDeltFx(x0); for(intj=0;j s0[j]=-g0[j]; } for(inti=0;i<2;i++,k++){ g0=getDeltFx(x0); g0_L=getDeltFx_Sk(s0,s0); L=getL(x0,s0);//例题一中的方法取得步长L for(intj=0;j x0[j]=x0[j]+s0[j]*L; } g1=getDeltFx(x0); g1_L=getDeltFx_Sk(g1,g1); if(Math.sqrt(g1_L)<=1e-2){ break; } else{ temp=g1_L/g0_L; for(intj=0;j s0[j]=-g1[j]+temp*s0[j]; } }} returnx0; } publicstaticvoidmain(String[]args){ GongEtemp=newGongE(); double[]x={1,1}; double[]result=temp.getX(x); for(inti=0;i System.out.println("result["+i+"]="+result[i]); }}} 2.2实验结果 最优点x*=[-4,-3,-1]最优解f*=-8 3.1公用程序(Java) publicstaticdoublegetFx(double[]x){//取得Fx值 doublex1=x[0];doublex2=x[1]; doubleFx=x1+2*x2*x2+Math.exp(x1*x1+x2*x2); returnFx; } publicstaticdouble[]getDeltFx(double[]x){//取得Fx的梯度值 doublex1=x[0];doublex2=x[1]; double[]deltFx=newdouble[2]; deltFx[0]=1+2*x1*Math.exp(x1*x1+x2*x2); deltFx[1]=4*x2+2*x2*Math.exp(x1*x1+x2*x2); returndeltFx; } 3.2.1最速下降法程序(Java) publicclassFastWay{ publicstaticdouble[]getX(double[]x){ double[]deltF0=newdouble[2];doubleL=0; for(inti=0;i<1e1;i++){ deltF0=getDeltFx(x); for(intj=0;j deltF0[j]=-deltF0[j]; } L=getL(x,deltF0);//调用习题1的不精确搜索取得步长L if(Math.sqrt(getDeltFx_Sk(deltF0,deltF0))<=1e-3){ System.out.println("最终计算次数"+i); System.out.println("x1="+x[0]+"x2="+x[1]); break; } x[0]=x[0]+L*deltF0[0];x[1]=x[1]+L*deltF0[1]; } returnx; } publicstaticvoidmain(String[]args){ FastWaytemp=newFastWay(); double[]x0={2,2};temp.getX(x0); } 3.2.2最速下降法结果 最优点X*=[-0.41940]最优解f*=0.7729计算次数count=10 3.3.1牛顿法程序(Java) publicstaticdouble[]getDeltFx(double[]x){ doublex1=x[0];doublex2=x[1]; double[]one=newdouble[2]; doubleexp=Math.exp(Math.pow(x1,2)+Math.pow(x2,2)); one[0]=1+2*x1*exp;one[1]=4*x2+2*x2*exp; double[][]two=newdouble[2][2]; two[0][0]=2*exp*(1+2*Math.pow(x1,2)); two[1][1]=2*exp*(1+2*Math.pow(x2,2))+4; double[]deltFx=newdouble[2]; for(inti=0;i<2;i++){ deltFx[0]=one[0]/two[0][0]; deltFx[1]=one[1]/two[1][1]; } returndeltFx; } publicstaticvoidmain(String[]args){ double[]x={1,0}; double[]DeltFx=newdouble[2]; for(inti=0;i<1e3;i++){ DeltFx=getDeltFx(x); x[0]=x[0]-DeltFx[0]; x[1]=x[1]-DeltFx[1]; if(Math.sqrt(getDeltFx_Sk(DeltFx,DeltFx))<=1e-4){ System.out.println("计算次数为"+i); break; } } System.out.println("x1="+x[0]+"x2="+x[1]+"\n"); System.out.println("Fx="+getFx(x)); } 3.3.2牛顿法结果 最优点X*=[-0.4194,0]最优解f*=0.7729计算次数count=5 3.4.1BFGS法程序(matlab) function[x,val,k]=bfgs(fun,gfun,x0) maxk=1000;sigma=0.4;rho=0.55;epsion=1e-5;k=0;n=length(x0); Bk=eye(n);%Bk=feval('Hess',x0); while(k gk=feval(gfun,x0); if(norm(gk) dk=-Bk\gk; m=0;mk=0; while(m<20) newf=feval(fun,x0+rho^m*dk) oldf=feval(fun,x0) if(newf mk=m;break; end m=m+1; end x=x0+rho^mk*dk; sk=x-x0; yk=feval(gfun,x)-gk; if(yk'*sk>0) Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk); end; k=k+1;x0=x; end val=feval(fun,x0); 3.4.2BFGS法结果 最优点X*=[-0.41940]最优解f*=0.7729计算次数count=4 4.1有效集法(matlab) 4.1.1主程序 function[x,Lagrange,exitflag,output]=TwoProg(H,c,Ae,be,Ai,bi,x0) n=length(x0);x=x0;ni=length(bi);ne=length(be);Lagrange=zeros(ne+ni,1);index=ones(ni,1); for(i=1: ni) if(Ai(i,: )*x>bi(i)+1e-9),index(i)=0;end end %算法主程序 k=0; while(k<=1e4) %求解子问题 Temp=[]; if(ne>0),Temp=Ae;end for(j=1: ni) if(index(j)>0),Temp=[Temp;Ai(j,: )];end end gk=H*x+c; [m1,n1]=size(Temp); [dk,Lagrange]=SubPro(H,gk,Temp,zeros(m1,1)); if(norm(dk)<=1.0e-6) y=0.0; if(length(Lagrange)>ne) [y,jk]=min(Lagrange(ne+1: length(Lagrange))); end if(y>=0) exitflag=0; else exitflag=1; for(i=1: ni) if(index(i)&(ne+sum(index(1: i)))==jk) index(i)=0;break; end end end k=k+1; else exitflag=1; %求步长 alpha=1.0;tm=1.0; for(i=1: ni) if((index(i)==0)&(Ai(i,: )*dk<0)) tm1=(bi(i)-Ai(i,: )*x)/(Ai(i,: )*dk); if(tm1 tm=tm1;ti=i; end end end alpha=min(alpha,tm); x=x+alpha*dk; if(tm<1),index(ti)=1;end end if(exitflag==0),break;end k=k+1; end output.fval=0.5*x'*H*x+c'*x; output.iter=k; 4.1.2目标函数 functionf=fun(x) x1=x (1);x2=x (2);f=eval('x1+2*x2^2+exp(x1^2+x2^2)'); 4.1.3子问题函数 function[x,Lagrange]=SubPro(H,c,Ae,be) [m,n]=size(Ae); ginvH=pinv(H); if(m>0) rb=Ae*ginvH*c+be; Lagrange=pinv(Ae*ginvH*Ae')*rb; x=ginvH*(Ae'*Lagrange-c); else x=-ginvH*c; Lagrange=0; end 4.1.4运行函数 H=[2-2;-24]; c=[-2-6]'; Ae=[];be=[]; Ai=[1-2;-0.5-0.5;10;01]; bi=[-2-100]'; x0=[01]'; [x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) 4.2有效集法结果 内部点初始点x0=[00]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=10 边界点初始点x0=[11]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=2 检验点初始点x0=[01]最优点X*=[0.81.2]最优解f*=-7.2迭代次数=7 5.1乘子法程序(matlab) 5.1.1chengZi程序---乘子法主程序 function[x,mu,Lagrange,output]=chengZi(fun,hf,gf,dfun,dhf,dgf,x0) sigma=2.0;count=0;innerCount=0; eta=2.0;θ=0.8;%PHR算法中的实参数θ x=x0;he=feval(hf,x);gi=feval(gf,x); n=length(x);l=length(he);m=length(gi); %选取乘子向量的初始值 mu=0.1*ones(l,1);Lagrange=0.1*ones(m,1); btak=10;btaold=10;%用来检验终止条件的两个值 while(btak>1e-6&count<1e3) %调用BFGS算法程序求解无约束子问题 [x,ival,ik]=bfgs('Lagr','LagrTiDu',x0,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma); innerCount=innerCount+ik; he=feval(hf,x);gi=feval(gf,x); btak=0.0; for(i=1: l),btak=btak+he(i)^2;end for(i=1: m) temp=min(gi(i),Lagrange(i)/sigma); btak=btak+temp^2; end btak=sqrt(btak); ifbtak>1e-6 if(count>=2&btak>θ*btaold) sigma=eta*sigma; end %更新乘子向量 for(i=1: l),mu(i)=mu(i)-sigma*he(i);end for(i=1: m) Lagrange(i)=max(0.0,Lagrange(i)-sigma*gi(i)); end end count=count+1; btaold=btak; x0=x; end f=feval(fun,x) output.inner_iter=innerCount; output.iter=count; output.bta=btak; output.fval=f; 5.1.2f1程序---目标函数 functionf=f1(x) f=4*x (1)-x (2)^2-12; 5.1.3h1程序---等式约束 functionhe=h1(x) he=25-x (1)^2-x (2)^2; 5.1.4g1程序---不等式约束 functiongi=g1(x) gi=10*x (1)-x (1)^2+10*x (2)-x (2)^2-34; 5.1.5df1程序---目标函数的梯度文件 functiong=df1(x) g=[4,-2.0*x (2)]'; 5.1.6dhe程序---等式约束(向量)函数的Jacobi矩阵(转置) functiondhe=dh1(x) dhe=[-2*x (1),-2.0*x (2)]'; 5.1.7dgi程序---不等式约束(向量)函数的Jacobi矩阵(转置) functiondgi=dg1(x) dgi=[10-2*x (1),10-2*x (2);0,1;1,0]'; 5.1.8LagrTiDu程序---增广拉格朗日函数的梯度程序 functionresult=LagrTiDu(x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma) result=feval(dfun,x); he=feval(hf,x);gi=feval(gf,x); dhe=feval(dhf,x);dgi=feval(dgf,x); l=length(he);m=length(gi); for(i=1: l) result=result+(sigma*he(i)-mu(i))*dhe(: i); end for(i=1: m) result=result+(sigma*gi(i)-Lagrange(i))*dgi(: i); 1)规划实施可能对相关区域、流域、海域生态系统产生的整体影响。 end (4)预防或者减轻不良环境影响的对策和措施的合理性和有效性;5.1.9Lagr程序---增广拉格朗日函数程序 functionresult=Lagr(x,fun,hf,gf,dfun,dhf,dgf,mu,Lagrange,sigma) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x); (一)环境影响经济损益分析概述l=length(he);m=length(gi); result=f;s1=0.0; 仍以
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 整理 大连理工大学 优化 作业 程序