2016年大连理工大学优化方法上机大作业

合集下载

Removed_大连理工大学工科数学分析上机作业

工科数学分析上机作业说明：以下两道题均是使用Matlab 语言，且在Matlab 7.0中运行通过。

1.（两个重要极限）计算下列函数的函数值并画出图形，观察两个重要极限值。

(1)y=f(x)=; (2)y=f(x)=.sin x x (1+x)1x 解：（1）求解过程如下：>> syms x>> y=limit(sin(x)/x)y =1>> ezplot(sin(x)/x,[-10*pi,10*pi])>> ezplot(sin(x)/x,[-1*pi,1*pi])其图形如下：（2）求解过程如下：>> syms x>> y=(1+x)^(1/x)y =(1+x)^(1/x)>> y=limit((1+x)^(1/x))y =exp(1)>> ezplot((1+x)^(1/x),[-1000,1000]) >> ezplot((1+x)^(1/x),[-10,10]) >> ezplot((1+x)^(1/x),[-1,1])其图像如下：分析如下：（1）当x 取值为[-30,30]时，由该题的第一个图像可以看到，函数值在不断震荡，一会为正数，一会为负数。

而当x 取值为[-3，3]时，函数值始终大于0。

当x 趋近于0时，由该题的第二个图像可以得到函数值为1。

另外，该结论也可以由夹逼法则证明，结果不变，当x 趋近于0时，函数值仍为1。

（2）由该题的三个图像可以知道，该函数在定义域内为单调递减函数。

且由该题的第一和二个图像知道，当x 在[0,10]区间内，函数递减趋势非常迅速。

由该题的第三个图像知道，当x 趋于0 时，函数值为自然对数的底数 e ,即约为2.71828.3.计算f(x)=,12+1√2π∫x 0e ‒t 2/2dt 1≪x ≪3的函数值{f （0.1k ）;k=1,2,…,30}.计算结果取7位有效数字。

教学大纲-大连理工大学教务处

目录《机械设计基础A》 (1)《机械设计基础B》 (8)《**模型设计概论》 (15)阅后删除：请以学部下设学院为单位将全部课程编辑在同一个文档内《机械设计基础A》教学大纲（学分4 学时64）一、课程说明(200字以内，简单说明本课程的地位及教学内容等，阅后删除红色字体）本课程是工科近机械类（包括机械类某些专业）和非机械类专业大类课程之一，是工科学生学习和掌握各种类型的机械中常用机构和通用机械零件的基本知识和基本设计方法的技术基础课。

该课程也是工科学生将来学习专业机械设备课程的理论基础。

本课程在教学内容方面着重基本知识、基本理论和基本设计方法的讲解；在培养实践能力方面着重设计构思和基本设计技能的基本训练。

二、课程目标（对应毕业要求：1-○1、1-○2、1-○3）1. 学习机械工程基础知识和基本理论知识，掌握常用机构的结构、特性等基本知识，了解各种机械的传动原理，具有分析、选用和设计机械设备中基本机构的能力（对应毕业要求：1-○1）；2. 通用机械零件的设计原理、方法和机械设计等的一般规律，具有设计机械传动装置和简单机械的能力（对应毕业要求：1-○1）；3. 掌握基本的机械设计创新方法，培养学生追求创新的态度和意识（对应毕业要求：1-○1）；4. 培养学生树立正确的设计思想，了解机械设计过程中国家有关的经济、环境、法律、安全、健康、伦理等政策和制约因素（对应毕业要求：1-○1）；5. 培养学生的工程实践学习能力，使学生掌握典型零件的实验方法，获得实验技能的基本训练，具有运用标准、规范、手册、图册和查阅有关技术资料的能力（对应毕业要求：1-○1）；6. 了解机械设计的前沿和新发展动向（对应毕业要求：1-○1）。

三、教学内容、基本要求与学时分配序号教学内容教学要求学时教学方式对应课程目标1 一、基本概念1. 研究的对象、内容；2. 机械设计的基本要求和一般设计过程。

1. 了解本课程研究的对象、内容2. 了解机械设计的基本要求、一般设计过程。

(整理)大连理工大学--优化作业----程序.

1.1程序（Java）public class Wolfe_Powell {public static double getFx ( double[] x ) {double x1= x[0]; double x2 = x[1];double Fx= 100 * (x2-x1*x1)* (x2-x1*x1) + (1-x1)* (1-x1) ;return Fx;}public static double[] getDeltFx ( double[] x ) {double x1= x[0]; double x2 = x[1];double [] deltFx = new double[2];deltFx [0] = -400*(x2 - x1* x1) *x1- 2*(1- x1) ;deltFx [1] = 200*(x2- x1 * x1) ;return deltFx ;}public static double getDeltFx_Sk ( double[] deltFx , double[] Sk ) { double a = 0 ;for ( int i = 0 ; i < Sk.length ; i++ ) {a = a + deltFx[ i ] * Sk[ i ] ;}return a ;}public static double getL ( double[] x, double[] s ) {double x1= x[0]; double x2 = x[1];double c1 =0.1 , c2 =0.5 ,a =0 , b=1e8 ,L= 1;double Fx0 , 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 ( int i=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 ;}else if( (Fx0 - Fx1 ) < temp2 ){b = L ;L = (L +a) /2 ;}else if ( deltFx1_Sk < temp ) {a = L ;L = ( L + b ) / 2 >= 2*L ? (2*L):( L + b ) / 2;}}System.out.println(" L= " + L);System.out.println(" 计算次数" + i );return L ;}public static void main(String[] args) {Wolfe_Powell temp =new Wolfe_Powell();double [] X = { -1 ,1 } ; double [] sk = { 1 ,1} ; temp.getL( X ,sk) ;}}1.2实验结果步长L = 0.00390625 x =[-0.9992 , 1.0324] 计算次数82.1程序（Java）public class GongE {public static double getFx ( double[] x ) {double x1= x[0];double x2 = x[1];double Fx= x1*x1 - 2*x1*x2 + 2*x2*x2 +x3*x3 - x2*x3 +2 * x1 +3*x2 -x3 ; return Fx;}public static double[] getDeltFx ( double[] x ) {double x1= x[0];double x2 = x[1];double [] deltFx = new double[x.length];deltFx [0] = 2*x1 - 2*x2+2 ;deltFx [1] = -2*x1 +4*x2 - x3 +3;deltFx [2] = 2*x3 -x2 -1 ;return deltFx ;}public static double[] getX ( double[] x ) {double[] g0,g1;double[] s0= new double[x.length];double[] s1=new double[x.length];double g0_L,g1_L ,L ,temp;double[] x0 =x ;int k =0 ;g0 = getDeltFx ( x0 ) ;for ( int j = 0 ; j < x.length ; j++ ) {s0[ j ] = -g0[ j ] ;}for (int i = 0 ;i<2; i ++,k++){g0 = getDeltFx ( x0 ) ;g0_L = getDeltFx_Sk ( s0 , s0 ) ;L =getL(x0,s0); // 例题一中的方法取得步长Lfor(int j=0;j<x.length ; 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(int j=0;j<x.length ; j++){s0[j] = -g1[j] + temp * s0[j];}} }return x0;}public static void main(String[] args) {GongE temp =new GongE();double [] x = { 1,1 } ;double[] result = temp.getX(x) ;for ( int i = 0 ; i < x.length ; i++ ) {System.out.println ( "result[" + i + "]=" + result[ i ] ) ;} } }2.2实验结果最优点x*=[-4,-3,-1] 最优解f*=-83.1公用程序（Java）public static double getFx ( double[] x ) { //取得Fx 值double x1= x[0]; double x2 = x[1];double Fx = x1 + 2 * x2 * x2 + Math.exp ( x1 * x1 + x2 * x2 ) ;return Fx ;}public static double[] getDeltFx ( double[] x ) { //取得Fx 的梯度值double x1= x[0]; double x2 = x[1];double[] deltFx = new double[ 2 ] ;deltFx[ 0 ] = 1 + 2 * x1 * Math.exp ( x1 * x1 + x2 * x2 ) ;deltFx[ 1 ] = 4 * x2 + 2 * x2 * Math.exp ( x1 * x1 + x2 * x2 ) ;return deltFx ;}3.2.1最速下降法程序（Java）public class FastWay {public static double[] getX ( double[] x ) {double [] deltF0 = new double[2]; double L =0;for ( int i = 0 ; i < 1e1 ; i++ ) {deltF0 = getDeltFx(x);for(int j=0 ;j <deltF0.length ;j++){ //取得负梯度deltF0[j] = - deltF0[j];}L = getL ( x , deltF0 ) ; // 调用习题1的不精确搜索取得步长Lif ( 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 ] ;}return x;}public static void main ( String[] args ) {FastWay temp = new FastWay () ;double[] x0 = { 2 , 2} ; temp.getX(x0) ;}3.2.2最速下降法结果最优点X*=[-0.4194 0] 最优解f*=0.7729 计算次数count=10 3.3.1牛顿法程序（Java）public static double[] getDeltFx ( double[] x ) {double x1 = x[ 0 ] ; double x2 = x[ 1 ] ;double[] one = new double[ 2 ] ;double exp =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 = new double[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 = new double[ 2 ] ;for (int i = 0 ; i < 2 ; i++ ) {deltFx[0] = one[ 0 ]/two[0][0] ;deltFx[1] = one[ 1 ]/two[1][1] ;}return deltFx;}public static void main ( String[] args ) {double[] x = { 1 , 0} ;double[] DeltFx = new double [2] ;for(int i =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.1 BFGS法程序（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<maxk)gk=feval(gfun,x0);if(norm(gk)<epsion),break;end;dk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk)oldf=feval(fun,x0)if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=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;endval=feval(fun,x0);3.4.2 BFGS法结果最优点X*=[-0.4194 0] 最优解f*=0.7729 计算次数count=44.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;endend%算法主程序k=0;while(k<=1e4)%求解子问题Temp=[];if(ne>0),Temp=Ae ; endfor(j=1:ni)if(index(j)>0),Temp=[Temp;Ai(j,:)];endendgk=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 )));endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk)index(i)=0;break;endendendk=k+1;elseexitflag=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)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1),index(ti)=1;endendif(exitflag==0),break;endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;4.1.2 目标函数function f=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);elsex=-ginvH*c;Lagrange =0;end4.1.4 运行函数H=[2 -2;-2 4];c=[-2 -6]';Ae=[ ];be=[ ];Ai=[1 -2;-0.5 -0.5;1 0;0 1];bi=[-2 -1 0 0]';x0=[0 1 ]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0)4.2 有效集法结果内部点初始点x0=[0 0] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=10 边界点初始点x0=[1 1] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=2 检验点初始点x0=[0 1] 最优点X*=[0.8 1.2] 最优解f*=-7.2 迭代次数=75.1 乘子法程序（matlab）5.1.1 chengZi程序---乘子法主程序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; endfor(i=1:m)temp=min(gi(i),Lagrange (i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if btak>1e-6if(count>=2&btak>θ*btaold)sigma=eta*sigma;end%更新乘子向量for(i=1:l),mu(i)=mu(i)-sigma*he(i);endfor(i=1:m)Lagrange (i)=max(0.0,Lagrange (i)-sigma*gi(i));endendcount=count+1;btaold=btak;x0=x;endf=feval(fun,x)output.inner_iter=innerCount;output.iter=count;output.bta=btak;output.fval=f;5.1.2 f1程序---目标函数function f=f1(x)f=4*x(1)-x(2)^2-12;5.1.3 h1程序---等式约束function he=h1(x)he=25-x(1)^2-x(2)^2;5.1.4 g1程序---不等式约束function gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;5.1.5 df1程序---目标函数的梯度文件function g=df1(x)g=[4 ,-2.0*x(2)]';5.1.6 dhe程序---等式约束（向量）函数的Jacobi矩阵（转置）function dhe=dh1(x)dhe=[-2*x(1),-2.0*x(2)]';5.1.7 dgi程序---不等式约束（向量）函数的Jacobi矩阵（转置）function dgi=dg1(x)dgi=[10-2*x(1),10-2*x(2);0,1;1,0]';5.1.8 LagrTiDu程序---增广拉格朗日函数的梯度程序function result=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);精品文档endfor(i=1:m)result=result+(sigma*gi(i)-Lagrange (i))*dgi(:,i);end5.1.9 Lagr程序---增广拉格朗日函数程序function result=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;for(i=1:l)result=result-he(i)*mu(i);s1=s1+he(i)^2;endresult=result+0.5*sigma*s1;s2=0.0;for(i=1:m)s3=max(0.0,Lagrange (i)-sigma*gi(i));s2=s2+s3^2-Lagrange (i)^2;endresult=result+s2/(2.0*sigma);5.2 乘子法结果初始点x0=[0 , 0] 最优点X*=[1.0013，4.8987] 最优解f*= -31.9923 等式乘子向量L hu=1.0156 不等式乘子向量Lg=0.75445精品文档。

优化方法大作业1

优化方法大作业一、Wolfe-Powell法利用MATLAB软件编写，其中初始值x0=-5；其他参数按照已知条件来取。

当b分别等于8、9、10时，均得到如下结果：而当初值x0变化时，则结果变化比较大，如将x0取-6，则计算结果如下：通过比较可以看出，b的取值对计算结果的影响较小，而初始值x0的取值则对结果影响很大。

从中也表明Wolfe-Powell准则的收敛条件比较弱，容易出现当函数还没取极小值而迭代循环已结束的情况。

具体代码见附录1.二、无约束优化1.DFP法精度ε = 0.02，确定λ使用的是一维非精确搜索算法(直接法，Wolfe-Powell准则)。

结果如下图所示：2. BFGS方法同上一种方法，精度ε = 0. 02，确定λ使用的是一维非精确搜索算法(直接法，Wolfe-Powell准则)。

结果如下图所示：比较两种方法，查看每一步的函数值，得出如下结果。

图1：DFP与BFGS法结果对比图通过图1与表1的比较，BFGS比DFP法最初收敛得更快，但是DFP法比BFGS法的最终结果更好。

DFP与BFGS法的代码分别见附录2、3。

三、Rockafellar乘子法文件myfun.m：function y=myfun(x)y=1000-x(1)^2-2*x(2)^2-x(3)^2-x(1)*x(2)-x(1)*x(3);文件mycon.m：function [c,ceq]=mycon(x)c(1)=(-x(1));c(2)=(-x(2));c(3)=(-x(3));ceq(1)=x(1)^2+x(2)^2+x(3)^2-25;ceq(2)=8*x(1)+14*x(2)+7*x(3)-56文件Q3.m：clearclcx0=[2 2 2]';[x,fval,exitflag,output]=fmincon(@myfun,x0,[],[],[],[],[],[],@mycon);附录1：%% Wolfe-Powell 准则方法clear;tic;c1=0.1;c2=0.65;a=0;b=8;syms x;f=x*x-2*x+7; %函数G0=jacobian(f,x);%求梯度x0=-6; %初始取x=-5% 由于只有一个未知数，默认S=1lambda = 1;x=x0;for k=1:1:1000f0 = eval(f); %求出函数值grad1= eval(G0); %求出梯度值figure1(k) = f0; % 用于绘出迭代过程中函数值变化x= x0+lambda;f1 = eval(f); %求出函数值grad2 = eval(G0); %求出梯度值value1 = f0 - f1 + c1* lambda * grad1 ;value2 = grad2 -c2*grad1 ;if value1 >=-1e-6 && value2 >=-1e-6disp('The variable matrix is : ');disp(x);fun=eval(f);fprintf('The minimum value of function is %f \n',fun);break; %如果满足两个条件，则退出循环。

大连理工优化方法大作业MATLAB编程

function [x,dk,k]=fjqx(x,s) flag=0;a=0;b=0;k=0;d=1;while(flag==0)[p,q]=getpq(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;x=x+d*s;flag=1;endk=k+1;if(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend%定义求函数值的函数fun，当输入为x0=（x1，x2）时，输出为f function f=fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;function gf=gfun(x)gf=[-4*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),2*(x(2)-x(1)^2)]; function [p,q]=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s';q=gfun(x+d*s)*s'-0.60*gfun(x)*s';结果：x=[0,1];s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54function f= fun( X )%所求问题目标函数f=X(1)^2-2*X(1)*X(2)+2*X(2)^2+X(3)^2+ X(4)^2-X(2)*X(3)+2*X(1)+3*X(2)-X(3);endfunction g= gfun( X )%所求问题目标函数梯度g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X(3)-X(2)-1,2*X(4)];endfunction [ x,val,k ] = frcg( fun,gfun,x0 )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值，k是迭代次数maxk=5000;%最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while(k<maxk)g=feval(gfun,x0);%计算梯度itern=k-(n+1)*floor(k/(n+1));itern=itern+1;%计算搜索方向if(itern==1)d=-g;elsebeta=(g*g')/(g0*g0');d=-g+beta*d0;gd=g'*d;if(gd>=0.0)d=-g;endendif(norm(g)<eps)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=[0,0,0,0];>> [ x,val,k ] = frcg( 'fun','gfun',x0 )x =-4.0000 -3.0000 -1.0000 0val =-8.0000k =21或者function [x,f,k]=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while(norm(g1)>=0.02)if(k==3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if(k<3)u=((norm(g1))^2)/(norm(g0)^2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x)f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-x(2)*x(3)+2*x(1)+3*x(2)-x(3); function gf=gfun(x)gf=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)];function [p,q]=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss)';q=gfun(x+d*ss)*(ss)'-0.6*gfun(x)*(ss)';function dk=dfun(x)flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endEnd结果：x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132 -1.0090 0 f = -7.9999k = 1function [f,x,k]=third_1(x)k=0;g=gfun(x);while(norm(g)>=0.001)s=-g;dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s)%获取步长flag=0;a=0;d=1;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend结果：x=[0,1];[f,x,k]=third_1(x)f =0.7729x = -0.4196 0.0001k =8（1）程序：function [f,x,k]=third_2(x)k=0;H=inv(ggfun(x));g=gfun(x);while(norm(g)>=0.001)s=(-H*g')';dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2); function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)]; function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)*exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2)]; function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)% 步长获取flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendEnd结果：x=[0,1];[f,x,k]=third_2(x)f =0.7729x = -0.4193 0.0001k =8（2）程序：function [f,x,k]=third_3(x) k=0;X=cell(2);g=cell(2);X{1}=x;H=eye(2);g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});while(norm(g{2})>=0.001)dx=X{2}-X{1};dg=g{2}-g{1};v=dx/(dx*dg')-(H*dg')'/(dg*H*dg'); h1=H*dg'*dg*H/(dg*H*dg');h2=dx'*dx/(dx*dx');h3=dg*H*dg'*v'*v;H=H-h1+h2+h3;k=k+1;X{1}=X{2};g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});norm(g{2});f=fun(x);x=X{2};endfunction f=fun(x)f=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)^2+x(2)^2),4*x(2)+2*x(2)*(x(1)^2+x(2)^2)];function ggf=ggfun(x)ggf=[(4*x(1)^2+2)*exp(x(1)^2+x(2)^2),4*x(1)*x(2)*exp(x(1)^2+x(2)^2);4*x(1)*x(2)* exp(x(1)^2+x(2)^2),4+(4*x(2)^2+2)*exp(x(1)^2+x(2)^2);function [p,q]=con(x,d,s)p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while(flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0) dk=d;flag=1;endif(p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.4195 0.0000 k=6function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[ ];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) function [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilon=1.0e-9; err=1.0e-6;k=0; x=x0; n=length(x); kmax=1.0e3;ne=length(be); ni=length(bi); lamk=zeros(ne+ni,1); index=ones(ni,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon), index(i)=0; endendwhile(k<=kmax)Aee=[];if(ne>0), Aee=Ae; endfor(j=1:ni)if(index(j)>0), Aee=[Aee; Ai(j,:)]; end endgk=H*x+c;[m1,n1] = size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if(norm(dk)<=err)y=0.0;if(length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk))); endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break;endendendk=k+1;elseexitflag=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)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1), index(ti)=1; end endif(exitflag==0), break; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x; output.iter=k;function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H);[m,n]=size(Ae);if(m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0)%功能: 用乘子法解一般约束问题: min f(x), s.t. h(x)=0, g(x).=0%输入: x0是初始点, fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束（向量）函数及其Jacobi矩阵的转置；% gf, dgf分别是不等式约束（向量）函数及其Jacobi矩阵的转置；%输出: x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量; % output是结构变量, 输出近似极小值f, 迭代次数, 内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;ink=0;epsilon=0.00001;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while(btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for i=1:lbtak=btak+he(i)^2;end%更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp^2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endfor i=1:mlambda(i)=max(0,lambda(i)-c*gi(i));endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*c);%不等式约束函数文件g1.mfunction gi=g1(x)gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;%目标函数的梯度文件df1.mfunction g=df1(x)g=[4, -2*x(2)]';%等式约束（向量）函数的Jacobi矩阵（转置）文件dh1.m function dhe=dh1(x)dhe=[-2*x(1), -2*x(2)]'%不等式约束（向量）函数的Jacobi矩阵（转置）文件dg1.m function dgi=dg1(x)dgi=[10-2*x(1), 10-2*x(2)]';function [x,val,k]=bfgs(fun,gfun,x0,varargin)maxk=500;rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while(k<maxk)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon)break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});结果x=[2 2]';[x,mu,lambda,output]=multphr('fun','hf','gf1','df','dh','dg',x0) x =1.00134.8987mu =0.7701lambda =0.9434output =fval: -31.9923iter: 4f=[3,1,1];A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0];x=linprog(f,A,b,zeros(3),[0,0,0]',lb)结果：Optimization terminated.x =0.00000.50000.5000。

大连理工大学结构优化复习总结

⼤连理⼯⼤学结构优化复习总结结构优化设计－基于结构分析技术，在给定的设计空间实现满⾜使⽤要求且具有最佳性能或最低成本的⼯程结构设计的技术优化设计的三要素:设计变量;约束条件;⽬标函数凸域:基于n维空间的区域s⾥，如果取任意两点x1和x2，连接这两点的线段也属于s，该区域称凸域（=αx1+(1-α)x2 ）凸函数:如果函数f(x)定义在n维空间的凸域s上，⽽且对s中的任意两点x1和x2和任意常数α,0.0<=α<=1.0,有f[αx1+(1- α)x2]<=αf(x1)+(1- α)f(x2),则f(x)称为s上的凸函数严格凸函数:上式⼩于严格成⽴凸规划:如果可⾏域是凸域，⽬标函数是凸函数，这样构成的数学规划问题为凸规划问题。

准则设计法:依靠⼯程经验;效率⾼;缺乏严格数学基础最优准则法基于库塔克（K-T）条件:需构造迭代求解算法;通⽤性不强数学规划⽅法:有严格的数学基础，有较好的通⽤性，计算效率要考虑。

结构优化问题的求解布骤I. 建⽴优化模型。

给定初始设计⽅案。

II. 结构分析(有限元）III.优化（收敛性）检验。

满⾜则结束程序，否则继续IVIV. 灵敏度分析V. 求解优化问题，修改结构模型，返回II。

优化求解的两⼤类⽅法：准则法；数学规划法准则设计⽅法：⽤优化准则代替原来的优化问题同步失效准则设计的评价：{优点：简单、⽅便，特别是独⽴约束个数n=m时；⼯程实⽤；适合于构件设计。

缺点：只能处理简单构件设计；缩⼩了设计空间，不能保证最优解；若n < m ，可能⽆解；当n > m时，确定哪些破坏模式应同时发⽣⽐较困难。

改进：为了弥补等式约束代替不等式约束的缺陷，引⼊松弛因⼦ψiσi (X ) =ψiσip , 0 ≤ψi ≤1, i =1,2,......n启发：⽤准则代替原来的优化问题，准则法的基本思想；如果将桁架的每根杆看作⼀种可能的破坏模式，桁架看作⼀个元件。

可以得到满应⼒准则满应⼒⽅法的缺点：完全⽆视重量会漏掉最轻设计；中间点⼀般是不可⾏设计，对⼯程实际不利。

大连理工大学优化方法上机作业

大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除This document is for reference only-rar21year.March优化方法上机大作业学院：电子信息与电气工程学部姓名：学号：指导老师：上机大作业（一）%目标函数function f=fun(x)f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;end%目标函数梯度function gf=gfun(x)gf=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; End%目标函数Hess矩阵function He=Hess(x)He=[1200*x(1)^2-400*x(2)+2,-400*x(1);-400*x(1), 200;];end%线搜索步长function mk=armijo(xk,dk)beta=0.5; sigma=0.2;m=0; maxm=20;while (m<=maxm)if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) mk=m; break;endm=m+1;endalpha=beta^mknewxk=xk+alpha*dkfk=fun(xk)newfk=fun(newxk)%最速下降法function [k,x,val]=grad(fun,gfun,x0,epsilon)%功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值maxk=5000; %最大迭代次数beta=0.5; sigma=0.4;k=0;while(k<maxk)gk=feval(gfun,x0); %计算梯度dk=-gk; %计算搜索方向if(norm(gk)<epsilon), break;end%检验终止准则m=0;mk=0;while(m<20) %用Armijo搜索步长if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+beta^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);>> x0=[0;0];>> [k,x,val]=grad('fun','gfun',x0,1e-4)迭代次数：k =1033x =0.99990.9998val =1.2390e-008%牛顿法x0=[0;0];ep=1e-4;maxk=10;k=0;while(k<maxk)gk=gfun(x0);if(norm(gk)<ep)x=x0miny=fun(x)k0=kbreak;elseH=inv(Hess(x0));x0=x0-H*gk;k=k+1;endendx =1.00001.0000miny =4.9304e-030迭代次数k0 =2%BFGS方法function [k,x,val]=bfgs(fun,gfun,x0,varargin) %功能：梯度法求解无约束优化问题：minf(x)%输入：fun，gfun分别是目标函数及其梯度，x0是初始点，% epsilon为容许误差%输出：k是迭代次数，x，val分别是近似最优点和最优值N=1000;epsilon=1e-4;beta=0.55;sigma=0.4;n=length(x0);Bk=eye(n);k=0;while(k<N)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon), break;enddk=-Bk\gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+beta^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<=oldf+sigma*beta^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+beta^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});>> x0=[0;0];>> [k,x,val]=bfgs('fun','gfun',x0)k =20x =1.00001.0000val =2.2005e-011%共轭梯度法function [k,x,val]=frcg(fun,gfun,x0,epsilon,N)if nargin<5,N=1000;endif nargin<4, epsilon=1e-4;endbeta=0.6;sigma=0.4;n=length(x0);k=0;while(k<N)gk=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)dk=-gk;elsebetak=(gk'*gk)/(g0'*g0);dk=-gk+betak*d0; gd=gk'*dk;if(gd>=0),dk=-gk;endendif(norm(gk)<epsilon),break;endm=0;mk=0;while(m<20)if(feval(fun,x0+beta^m*dk)<=feval(fun,x0)+sigma*beta^m*gk'*dk) mk=m;break;endm=m+1;endx=x0+beta^m*dk;g0=gk; d0=dk;x0=x;k=k+1;endval=feval(fun,x);>> x0=[0;0];[k,x,val]=frcg('fun','gfun',x0,1e-4,1000)k =122x =1.00011.0002val =7.2372e-009上机大作业（二）%目标函数function f_x=fun(x)f_x=4*x(1)-x(2)^2-12;%等式约束条件function he=hf(x)he=25-x(1)^2-x(2)^2;end%不等式约束条件function gi_x=gi(x,i)switch icase 1gi_x=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;case 2gi_x=x(1);case 3gi_x=x(2);otherwiseend%求目标函数的梯度function L_grad=grad(x,lambda,cigma)d_f=[4;2*x(2)];d_g(:,1)=[-2*x(1);-2*x(2)];d_g(:,2)=[10-2*x(1);10-2*x(2)];d_g(:,3)=[1;0];d_g(:,4)=[0;1];L_grad=d_f+(lambda(1)+cigma*hf(x))*d_g(:,1);for i=1:3if lambda(i+1)+cigma*gi(x,i)<0L_grad=L_grad+(lambda(i+1)+cigma*gi(x,i))*d_g(:,i+1);continueendend%增广拉格朗日函数function LA=lag(x,lambda,cee)LA=fun(x)+lambda(1)*hf(x)+0.5*cee*hf(x)^2;for i=1:3LA=LA+1/(2*cee)*(min(0,lambda(i+1)+cee*gi(x,i))^2-lambda(i+1)^2); endfunction xk=BFGS(x0,eps,lambda,cigma)gk=grad(x0,lambda,cigma);res_B=norm(gk);k_B=0;a_=1e-4;rho=0.5;c=1e-4;length_x=length(x0);I=eye(length_x);Hk=I;while res_B>eps&&k_B<=10000dk=-Hk*gk;m=0;while m<=5000if lag(x0+a_*rho^m*dk,lambda,cigma)-lag(x0,lambda,cigma)<=c*a_*rho^m*gk'*dkmk=m;break;endm=m+1;endak=a_*rho^mk;xk=x0+ak*dk;delta=xk-x0;y=grad(xk,lambda,cigma)-gk;Hk=(I-(delta*y')/(delta'*y))*Hk*(I-(y*delta')/(delta'*y))+(delta*delta')/(delta'*y);k_B=k_B+1;x0=xk;gk=y+gk;res_B=norm(gk);end%增广拉格朗日法function val_min=ALM(x0,eps)lambda=zeros(4,1);cigma=5;alpha=10;k=1;res=[abs(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);while res>eps&&k<1000xk=BFGS(x0,eps,lambda,cigma);lambda(1)=lambda(1)+cigma*hf(xk);for i=1:3lambda(i+1)=lambda(i+1)+min(0,lambda(i+1)+gi(x0,1)); endk=k+1;cigma=alpha*cigma;x0=xk;res=[norm(hf(x0)),0,0,0];for i=1:3res(1,i+1)=norm(min(gi(x0,i),-lambda(i+1)/cigma)); endres=max(res);endval_min=fun(xk);fprintf('k=%d\n',k);fprintf('fmin=%.4f\n',val_min);fprintf('x=[%.4f;%.4f]\n',xk(1),xk(2));>> x0=[0;0];>> val_min=ALM(x0,1e-4)k=10fmin=-31.4003x=[1.0984;4.8779]val_min =-31.4003上机大作业（三）A=[1 1;-1 0;0 -1];n=2;b=[1;0;0];G=[0.5 0;0 2];c=[2 4];cvx_solver sdpt3cvx_beginvariable x(n)minimize (x'*G*x-c*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -2.40.40000.6000A=[2 1 1;1 2 3;2 2 1;-1 0 0;0 -1 0;0 0 -1]; n=3;b=[2;5;6;0;0;0];C=[-3 -1 -3];cvx_solver sdpt3cvx_beginvariable x(n)minimize (C*x)subject toA*x<=bcvx_enddisp(x)Status: SolvedOptimal value (cvx_optval): -5.40.20000.00001.600011。

大连理工优化方法大作业MATLAB编程

fun ctio n [x,dk,k]=fjqx(x,s) flag=0;a=0;b=0;k=0;d=1;while (flag==0)[p,q]=getpq(x,d,s);if (P<0)b=d;d=(d+a)/2;endif(p>=0) &&( q>=0)dk=d;x=x+d*s;flag=1;endk=k+1;if (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endend%定义求函数值的函数 fun ，当输入为 x0= (x1 ， x2 )时，输出为 f function f=fun(x)f=(x(2)-x(1)A2)A2+(1-x(1)F2;function gf=gfun(x)gf=[-4*x(1)*(x (2) -x(1)A2)+2*(x(1)-1),2*(x(2)-x(1)A2)];function [p,q]=getpq(x,d,s)p=fun(x)-fun(x+d*s)+0.20*d*gfun(x)*s';q=gfun(x+d*s)*s'-0.60*gfun(x)*s';结果：x=[0,1];s=[-1,1];[x,dk,k]=fjqx(x,s)x =-0.0000 1.0000dk =1.1102e-016k =54取初始= (0.0. 0,0)r^'l用兵柜梯皮法求解下面无约東优化问题：min f (x) = x孑—2x^X2 十2x孑 + x孑H-爲—X2天3 十 2xj + 3|X2 —*3,其中步长g的选取可利用习題1戎精确一维披索.注:通过比习题验证共範梯度法求辉门无二次西数极小点至多需要“次迭代.fun ctio n f= fun( X )%所求问题目标函数f=X(1)A2-2*X(1)*X (2)+2*X(2)A2+X(3)A2+ X(4) A2-X( 2)*X(3)+2*X(1)+3*X(2)-X(3);end function g= gfun( X )%所求问题目标函数梯度g=[2*X(1)-2*X(2)+2,-2*X(1)+4*X(2)-X(3)+3,2*X (3) -X (2)-1,2*X(4)];end function [ x,val,k ] = frcg( fun,gfun,xO )%功能：用FR共轭梯度法求无约束问题最小值%输入：x0是初始点，fun和gfun分别是目标函数和梯度%输出：x、val分别是最优点和最优值，k是迭代次数maxk=5000; %最大迭代次数rho=0.5;sigma=0.4;k=0;eps=10e-6;n=length(x0);while (k<maxk)g=feval(gfun,x0); % 计算梯度 itern=k-(n+1)*floor(k/(n+1));itern=itern+1;%计算搜索方向if (itern==1)d=-g;elsebeta=(g*g')/(g0*g0');d=-g+beta*d0;gd=g'*d;if (gd>=0.0)d=-g;endendif (norm(g)<eps)break ;endm=0;mk=0;while (m<20)if(feval(fu n,xO+rhoAm*d)<feval(fu n,xO)+sigma*rhoAm*g'*d) mk=m; break ;endm=m+1;endx0=x0+rho A mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x0);end结果：>> x0=[0,0,0,0];>> [ x,val,k ] = frcg( 'fun','gfun',x0 ) x =-4.0000 -3.0000 -1.0000 0val =-8.0000k =或者function [x,f,k]=second(x)k=0;dk=dfun(x);g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);while (norm(g1)>=0.02)if (k==3)k=0;g0=gfun(x);s=-g0;x=x+dk*s;g1=gfun(x);else if (k<3)u=(( norm(g1))A2)/( norm(gO)A2); s=-g1+u*s;k=k+1;g0=g1;dk=dfun(x);x=x+dk*s;g1=gfun(x);endendf=fun(x);endfunction f=fun(x)f=x(1F2-2*x(1)*x (2)+2*x (2)A2+x(3)A2+x(4)A2-x (2) *x (3)+2*x(1)+3*x(2)-x(3); function gf=gfun(x)gf=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)];function [p,q]=con(x,d)ss=-gfun(x);p=fun(x)-fun(x+d*ss)+0.2*d*gfun(x)*(ss)';q=gfun(x+d*ss)*(ss)'-0.6*gfun(x)*(ss)';function dk=dfun(x)flag=0;a=0;d=1;while (flag==0)[p,q]=con(x,d);if (p<0)b=d;d=(d+a)/2;endif (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2};endEnd结果： x=[0,0,0,0];>> [x,f,k]=second(x)x =-4.0147 -3.0132-1.0090 0 f = -7.9999k = 1取初始点3 = （0」）二考虑下面无约東优化问题：min f（x）二冷 + 2x2 + exp（xf + 天孑），其中歩长Qk的选取可別用习题1或精确一维搜索•搜索方向为一HNW ♦取垃=b•取皿=R2f防）］"9耳丈啟为BFG5公式亠通过此习题体会上述三种算法的收敛速度.fun ctio n [f,x,k]=third_1(x) k=0;g=gfu n(x);while (norm(g)>=0.001) s=-g;dk=dfu n( x,s);x=x+dk*s;k=k+1;g=gfu n(x);f=fun( x);endfun ctio n f=fun(x)f=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);fun ctio n gf=gfu n(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)]; function[j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)'; j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s) % 获取步长 flag=0;a=0;d=1;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;d=min{2*d,(d+b)/2}; end结果：x=[0,1];[f,x,k]=third_1(x)f =0.7729x = -0.4196 0.0001k =8(1 ) 程序：function [f,x,k]=third_2(x)k=0;H=inv(ggfun(x));g=gfun(x);while (norm(g)>=0.001)s=(-H*g')';dk=dfun(x,s);x=x+dk*s;k=k+1;g=gfun(x);f=fun(x);endfunction f=fun(x)f=x(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x) gf=[1+2*x(1)*exp(x(1F2+x(2)A2),4*x(2)+2*x(2)*(x(1F2+x(2)A2)]; function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x (2) A2),4*x(1)*x (2) *exp(x(1)A2+x(2)A2);4*x(1)*x(2)*exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2)];function [j_1,j_2]=con(x,d,s)j_1=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';j_2=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)'; function dk=dfun(x,s) % 步长获取flag=0;a=0;d=1;b=10000;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;endif(p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;endendEnd结果：x=[0,1];[f,x,k]=third_2(x)f =0.7729x = -0.4193 0.0001k =8(2) 程序：function [f,x,k]=third_3(x) k=0;X=cell(2);g=cell(2);X{1}=x;H=eye(2);g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});while (norm(g{2})>=0.001)dg=g{2}-g{1};v=dx/(dx*dg')-(H*dg')'/(dg*H*dg');h1=H*dg'*dg*H/(dg*H*dg');h2=dx'*dx/(dx*dx');h3=dg*H*dg'*v'*v;H=H-h1+h2+h3;k=k+1;X{1}=X{2};g{1}=gfun(X{1});s=(-H*g{1}')';dk=dfun(X{1},s);X{2}=X{1}+dk*s;g{2}=gfun(X{2});norm(g{2});f=fun(x);x=X{2};endfunction f=fun(x)f=x(1)+2*x(2)A2+exp(x(1F2+x(2)A2);function gf=gfun(x)gf=[1+2*x(1)*exp(x(1)A2+x(2)A2),4*x(2)+2*x(2)*(x(1)A2+x(2)A2)];function ggf=ggfun(x)ggf=[(4*x(1)A2+2)*exp(x(1)A2+x(2)A2),4*x(1)*x(2)*exp(x(1)A2+x(2)A2);4*x(1)*x(2)* exp(x(1)A2+x(2)A2),4+(4*x(2)A2+2)*exp(x(1)A2+x(2)A2);function [p,q]=con(x,d,s)p=fun(x)-fun(x+d*s)+0.1*d*gfun(x)*(s)';q=gfun(x+d*s)*(s)'-0.5*gfun(x)*(s)';function dk=dfun(x,s)flag=0;a=0;d=1;b=10000;while (flag==0)[p,q]=con(x,d,s);if (p<0)b=d;d=(d+a)/2;if (p>=0)&&(q>=0)dk=d;flag=1;endif (p>=0)&&(q<0)a=d;if 2*d>=(d+b)/2d=(d+b)/2;else d=2*d;endendend结果：x=[0,1];[f,x,k]=third_3(x)f =0.7729x = -0.41950.0000 k=6*U 用有效集法求解下面勺勺二次规划问题:(XI 一 I)2 + (x 2 一 2.5)2 X1 - 2X2 + 2 > 0-Xi — 2>(2 + 6 > 0-Xi + 2X2 + 2 > 0xi,x 2 > 0function callqpactH=[2 0; 0 2];c=[-2 -5]';Ae=[ ]; be=[];Ai=[1 -2; -1 -2; -1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';[x,lambda,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,xO)fun ctio n [x,lamk,exitflag,output]=qpact(H,c,Ae,be,Ai,bi,x0) epsilo n=1.0e-9; err=1.0e-6;k=0; x=x0; n=len gth(x); kmax=1.0e3;n e=le ngth(be); ni=le ngth(bi); lamk=zeros( ne+n i,1); in dex=ones(n i,1);for (i=1:ni)if(Ai(i,:)*x>bi(i)+epsil on), i ndex(i)=0; end while (k<=kmax)mmSi.Aee=[];if (ne>0), Aee=Ae; endfor (j=1:ni)if (index(j)>0), Aee=[Aee; Ai(j,:)]; end endgk=H*x+c;[m1,n1] = size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1)); if (norm(dk)<=err)y=0.0;if (length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk))); endif (y>=0)exitflag=0;elseexitflag=1;for (i=1:ni)if (index(i)&(ne+sum(index(1:i)))==jk) index(i)=0; break ;endendk=k+1;elseexitflag=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)tm=tm1; ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if (tm<1), index(ti)=1; endendif (exitflag==0), break ; endk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;function [x,lambda]=qsubp(H,c,Ae,be) ginvH=pinv(H); [m,n]=size(Ae);if (m>0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae')*rb; x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;end结果>>callqpactx =1.40001.7000lambda =0.8000exitflag =output =fval: -6.4500iter: 7function [x,mu,lambda,output]=multphr(fu n, hf,gf,dfu n, dhf,dgf,xO)%功能：用乘子法解一般约束问题：min f(x), s.t. h(x)=0, g(x).=0%输入：x0是初始点，fun, dfun分别是目标函数及其梯度；% hf, dhf分别是等式约束(向量)函数及其 Jacobi矩阵的转置；% gf, dgf分别是不等式约束(向量)函数及其 Jacobi矩阵的转置；%输出：x是近似最优点，mu, lambda分别是相应于等式约束和不等式约束的乘子向量% output是结构变量，输出近似极小值f,迭代次数，内迭代次数等maxk=500;c=2.0;eta=2.0;theta=0.8;k=0;i nk=0;epsilo n=0.00001;x=xO;he=feval(hf,x);gi=feval(gf,x);n=len gth(x);l=le ngth(he);m=le ngth(gi);mu=zeros(l,1);lambda=zeros(m,1);btak=10;btaold=10;while (btak>epsilon&&k<maxk)%调用BFGS算法程序求解无约束子问题[x,ival,ik]=bfgs( 'mpsi' ,'dmpsi' ,x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for i=1:lbtak=btak+he(y2;end% 更新乘子向量for i=1:mtemp=min(gi(i),lambda(i)/c);btak=btak+temp A2;endbtak=sqrt(btak);if btak>epsilonif k>=2&&btak>theta*btaoldc=eta*c;endfor i=1:lmu(i)=mu(i)-c*he(i);endlambda(i)=max(0,lambda(i)-c*gi(i));endk=k+1;btaold=btak;x0=x;endendf=feval(fun,x);output.fval=f;output.iter=k;%增广拉格朗日函数function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,c) f=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for i=1:lpsi=psi-he(i)*mu(i);s仁 s1+he(y2;psi=psi+0.5*c*s1;s2=0;for i=1:ms3=max(0,lambda(i)-c*gi(i));s2=s2+s3A2-lambda(i)A2;endpsi=psi+s2/(2*c);% 不等式约束函数文件 g1.mfunction gi=g1(x)gi=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;% 目标函数的梯度文件df1.mfunction g=df1(x)g=[4, -2*x(2)]';% 等式约束(向量)函数的Jacobi 矩阵(转置)文件 dh1.m function dhe=dh1(x)dhe=[-2*x(1), -2*x(2)]'% 不等式约束(向量)函数的Jacobi 矩阵(转置)文件 dg1.m function dgi=dg1(x)dgi=[10-2*x(1), 10-2*x(2)]';function [x,val,k]=bfgs(fun,gfun,x0,varargin) maxk=500; rho=0.55;sigma=0.4;epsilon=0.00001;k=0;n=length(x0);Bk=eye(n);while (k<maxk)gk=feval(gfun,x0,varargin{:});if (norm(gk)<epsilon)break ;enddk=-Bk\gk;m=0;mk=0;while (m<20)n ewf=feval(fu n, x0+rho A m*dk,vararg in {:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rhoAm*gk'*dk) mk=m;break ;endm=m+1;endx=x0+rhoAmk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if (yk'*sk>0)Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});结果x=[2 2]';[x,mu,lambda,output]=multphr( 'fun' ,'hf' ,'gf1' ,'df' ,'dh' ,'dg' ,x0) x =1.00134.8987mu =0.7701lambda =0.9434output =fval: -31.9923iter: 4利用序列二次规划方法求解习题5中的约束优化问题:min 4xi 一好一 12s.t. 25 - x? —x孑=Q10x一召 + 10旳-xj - 34 > 0 X1,X2 > 0tf=[3,1,1]；A=[2,1,1;1,-1,-1];b=[2;-1];lb=[0,0,0]； x=li nprog(f,A,b,zeros(3),[0,0,0]',lb)结果：Optimization terminated.0.00000.50000.5000。

大连理工高等物理化学大作业

共价有机骨架（COFs）材料的合成与应用引言：在过去的十年里，由于纳米多孔材料在气体存储、气体分离、超疏水相互作用、催化、能量转换、能量存储和光电性能等方面的出色表现和广泛应用而引起了众多化学家的兴趣。

化学家们已经发现并合成了多种多样的孔材料。

但是，在网状化学的概念提出以前，对于共价有机骨架多孔材料的合成是十分困难的。

网状化学运用拓扑学的方法，设计并构建模块以指导孔材料的合成1。

运用网状化学成功合成的第一类孔材料是金属有机骨架（MOFs）材料，在这种材料中，金属离子作为中心离子，有机化合物作为配体，以配位键形成有机无机杂化多孔材料。

MOFs材料因其中心离子和有机配合物的多样性以及材料的可修饰性，展现出广阔的发展空间和应用潜力。

2005年，Yaghi小组首先报道了利用拓扑学理论设计合成出由共价键连接的孔材料，即共价有机骨架（COFs）材料2。

这篇标志性文章的发表对COFs材料的化学合成产生了深远的影响，并且该材料在功能探索方面也具有巨大的潜力1。

因为COFs材料是由轻元素通过共价键连接的，所以具有较低的分子质量和较好的热力学稳定性。

依据COFs材料的维度，可以将其分为2D和3D的COFs 材料。

2D的COFs材料是层状结构，类似于石墨，其在光电性能上表现出很好的潜力，由于其中碳原子采用sp2杂化，形成离域电子体系，是很好的导电材料。

而3D的COFs材料是碳原子通过sp3杂化成键，形成三维立体结构，具有极高的比表面积，而其质量又很低，估计在储氢上会有应用。

本文将从COFs材料的合成和应用两方面来介绍这种新型材料。

1.COFs材料的合成对于可逆反应，在化学反应进程中，不仅包含反应物化学键的断裂和生成物化学键的形成，还包括生成物的化学键断裂和再形成。

COFs材料的合成过程便是可逆过程，这种键的形成、断裂、再形成的过程，给予材料在合成过程中自我校验的能力，修复不规整的部分并形成较为完美的晶体多孔聚合物。

所以实现热力学控制的可逆反应就成为一个最重要的问题。

大连理工大学庞丽萍最优化方法MATLAB程序

班级：优化1班授课老师：庞丽萍姓名：学号：第二章12.（1）用修正单纯形法求解下列LP问题:>>clear>>A=[121100;123010;215001];[m,n]=size(A);b=[10;15;20];r=[-1-2-31];c=[-1-2-31];bs=[3:3];nbs=[1:4];a1=A(:,3);T=A(:,bs);a2=inv(T)*a1;b=inv(T)*b;A=[eye(m),a2];B=eye(m);xb=B\b;cb=c(bs);cn=c(nbs);con=1;M=zeros(1);while conM=M+1;t=cb/B;r=c-t*A;if all(r>=0)x(bs)=xb;x(nbs)=0;fx=cb*xb;disp(['当前解是最优解，minz=',num2str(fx)])disp('对应的最优解为,x=')disp(x)breakendrnbs=r(nbs);kk=find(rnbs==min(rnbs));k=kk(1);Anbs=A(:,nbs);yik=B\Anbs(:,k);xb=B\b;%yi0if all(yik<=0)disp('此LP问题无有限的最优解，计算结束',x)disp(xb)breakelsei=find(yik>0);w=abs(xb(i,1)./yik(i,1));l=find(w==min(w));rr=min(l);yrrk=yik(rr,1);Abs=A(:,bs);D=Anbs(:,k);Anbs(:,k)=Abs(:,rr);Abs(:,rr)=D;F=bs(rr);bs(rr)=nbs(k);nbs(k)=F;AA=[Anbs,Abs];EE=eye(m);EE(:,rr)=-yik./yrrk;Errk=EE;Errk(rr,rr)=1/yrrk;BB=Errk/B;B=inv(BB);cb=c(:,bs);xb=Errk*xb;x(bs)=xb;x(nbs)=0;fx=cb*xb;endif M>=1000disp('此问题无有限最优解')breakendend%结果当前解是最优解，minz=-15对应的最优解为,x=2.5000 2.5000 2.50000第三章30题DFP算法求函数极小点的计算程序function[x,val,k]=dfp(fun,gfun,x0)%功能:用DFP算法求解无约束问题:minf(x)%输入:x0是初始点,fun,gfun分别是目标函数及其梯度%输出:x,val分别是近似最优点和最优值,k是迭代次数.maxk=1e5;%给出最大迭代次数rho=0.55;sigma=0.4;epsilon=1e-5;k=0;n=length(x0);Hk=inv(feval('Hess',x0));%Hk=eye(n);while(k<maxk)gk=feval(gfun,x0);%计算梯度if(norm(gk)<epsilon),break;end%检验终止准则dk=-Hk*gk;%解方程组,计算搜索方向m=0;mk=0;while(m<20)%用Armijo搜索求步长if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk’*dk)mk=m;break;endm=m+1;end%DFP校正x=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x)-gk;if(sk'*yk>0)Hk=Hk-(Hk*yk*yk'*Hk)/(yk'*Hk*yk)+(sk*sk')/(sk'*yk);endk=k+1;x0=x;endval=feval(fun,x0);%习题26的程序调用方式及结果：function y=fun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=(x(1)-1)^2+5*(x2-x(1)^2)^2endfunction y=gfun(x)%UNTITLED Summary of this function goes here%Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[20]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=1.000001.00000val=k=6%习题27的程序调用方式及结果：function y=fun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=x1+2*x（2）^2+exp(x(1)^2+x(2)^2)endfunction y=gfun(x)%UNTITLED Summary of this function goes here %Detailed explanation goes herey=[diff(y,x1)diff(y,x2)]endx0=[10]’;[x,val,k]=dfp(fun,gfun,x0)%结果x=-0.419360val=0.77291k=536题编写Hooke-Jeeves方法求函数极小点的计算程序。

大连理工大学秋季优化方法大作业

m=m+1; end x0=x0+rho^mk*d; val=feval(fun,x0); g0=g; d0=d; k=k+1; end x=x0; val=feval(fun,x);
//f(x)//
function f=fun(x) x1=[1 0]*x; x2=[0 1]*x; f=(1-x1)^2+100*(x2-x1^2)^2; //梯度函数// function g=gfun(x) x1=[1 0]*x; x2=[0 1]*x; g=[-2*(1-x1)-400*x1*(x2-x1^2); 200*(x2-x1^2)]; //运行过程// >> x0=[0 0]'
k=k+1; btaold=btak; x0=x; end f=feval(fun,x); output.fval=f; output.iter=k; output.inner_iter=ink; output.bta=btak;
//增广拉格朗日函数//
function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x); l=length(he); m=length(gi); psi=f; s1=0.0; for(i=1:l) psi=psi-he(i)*mu(i); s1=s1+he(i)^2; end psi=psi+0.5*sigma*s1; s2=0.0; for(i=1:m) s3=max(0.0, lambda(i) - sigma*gi(i)); s2=s2+s3^2-lambda(i)^2; end psi=psi+s2/(2.0*sigma); //增广拉格朗日函数// function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma) dpsi=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) dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i); end for(i=1:m) dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i); end //f(x)// function f=f1(x) f=4*x(1)-x(2)^2-12; //等式约束// function he=h1(x) he=25-x(1)^2-x(2)^2;

最优化方法大作业

发动机空燃比控制器引言：我主要从事自动化相关研究。

这里介绍我曾经接触过的发动机空燃比控制器设计中的优化问题。

发动机空燃比控制器设计中的最优化问题AFR =afm m && (1)空燃比由方程（1）定义，在发动机运行过程中如果控制AFR 稳定在14.7可以获得最好的动力性能和排放性能。

如果假设进入气缸的空气流量am &可以由相关单元检测得到，则可以通过控制进入气缸的燃油流量f m &来实现空燃比的精确控制。

由于实际发动机的燃油喷嘴并不是直接对气缸喷燃油，而是通过进气歧管喷燃油，这么做会在进气歧管壁上液化形成油膜，因此不仅是喷嘴喷出的未液化部分燃油会进入气缸，油膜蒸发部分燃油也会进入气缸，如方程（2）。

这样如何更好的喷射燃油成为了一个问题。

1110101122211ττττ⎡⎤⎡⎤-⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥⎢⎥=+⎢⎥⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥-⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤⎢⎥⎣⎦⎢⎥⎣⎦f ff v X x x u x x X x y =x && (2)其中12、,==ff fv x m x m &&=f y m &,=fi u m &这里面，表示油膜蒸发量ff m &、fvm &表示为液化部分燃油、fim &表示喷嘴喷射的燃油，在τf 、τv 、X 都已知的情况下，由现代控制理论知识，根据系统的增广状态空间模型方程（3）00000011011011114.70ττττ⎡⎤⎡⎤-⎡⎤⎢⎥⎢⎥⎡⎤⎡⎤⎢⎥=-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦⎡⎤⎡⎤=⎢⎥⎣⎦⎣⎦ff v v a X X u +q q m y q x x x &&& (3)其中()014.7⎰taq =y -m&。

由极点配置方法，只要设计控制器方程（4），就可以使得y 无差的跟踪阶跃输入，那么y 也能较好的跟踪AFR *am /&。

优化方法课程大作业

优化方法课程上机大作业学部：电子信息与电气工程学部专业：生物医学工程班级：电信硕1303学号：21309210姓名：史益新大连理工大学Dalian University of Technology解：（1）MATLAB代码如下：clc;clear all;close all;[x,y]=meshgrid(-2:0.1:2,-1:0.1:3);z=(y-x.^2).^2+(1-x).^2;mesh(x,y,z)hold on;xk=[0 1]';epsilon=1e-5;plot(xk(1),xk(2),'ro');text(xk(1),xk(2),'start point');hold on;[ x,val,k ]=Newton('fun','gfun','Hess',xk,epsilon)plot(x(1),x(2),'ro');text(x(1),x(2),'end point');function [ x,val,k ] = Newton( fun,gfun,Hess,xk,epsilon ) k=0;while(1)gk=feval(gfun,xk);hk=feval(Hess,xk);sk=-inv(hk)*gk;if(norm(gk)<epsilon)break;endxk=xk+sk;k=k+1;endx=xk;val=feval(fun,x);endfunction f = fun(x)f=(x(2)-x(1)^2)^2+(1-x(1))^2;endfunction g=gfun(x)g=[-4*x(1)*(x(2)-x(1)^2)+2*(1-x(1)),2*(x(2)-x(1)^2)]'; endfunction He = Hess( x )n=length(x);He=zeros(n,n);He=[12*x(1)^2-4*x(2)+2,-4*x(1);-4*x(1), 2 ];end（2）代码运行结果如下：解：（1）MATLAB代码如下：clc;clear all;close all;x0=[0 0 0 0]';epsilon=1e-5;[ x,val,k ]=Frcg('fun','gfun',x0,epsilon)function [ x,val,k ] = Frcg( fun,gfun,x0,epsilon )rho=0.6;sigma=0.5;k=0;n=length(x0);while(1)g=feval(gfun,x0);itern=k-(n+1)*floor(k/(n+1));itern=itern+1;if(itern==1)d=-g;elsebeta=(g'*g)/(g0'*g0);d=-g+beta*d0;gd=g'*d;if(gd>=0.0)d=-g;endendif(norm(g)<epsilon)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^mk*d;val=feval(fun,x0);g0=g;d0=d;k=k+1;endx=x0;val=feval(fun,x);function f = fun( x )f=x(1)^2-2*x(1)*x(2)+2*x(2)^2+x(3)^2+x(4)^2-x(2)*x(3)+2*x(1)+3*x(2)-x(3); endfunction g = gfun( x )g=[2*x(1)-2*x(2)+2,-2*x(1)+4*x(2)-x(3)+3,2*x(3)-x(2)-1,2*x(4)]';end（2）代码运行结果如下：解：（1）MATLAB代码如下：clc;clear all;close all;[x,y]=meshgrid(-2:0.1:2,-1:0.1:3);z=5*(y-x.^2).^2+(x-1).^2;mesh(x,y,z)hold on;x0=[2 0]';epsilon=1e-5;plot(x0(1),x0(2),'ro');text(x0(1),x0(2),'start point');hold on;[ x1,val1,k1 ]=grad('fun','gfun',x0,epsilon)plot(x1(1),x1(2),'ro');text(x1(1),x1(2),'end point1');hold on;[ x2,val2,k2 ]=znNewton('fun','gfun','Hess',x0,epsilon) plot(x2(1),x2(2),'bo');text(x2(1),x2(2),'end point2');hold on;[ x3,val3,k3 ]=BFGS('fun','gfun',x0,epsilon)plot(x3(1),x3(2),'go');text(x3(1),x3(2),'end point3');hold on;function [ x,val,k ] = grad( fun,gfun,x0,epsilon )rho=0.5;sigma=0.4;k=0;while(1)g=feval(gfun,x0);d=-g;if(norm(d)<epsilon)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m;break;endm=m+1;endx0=x0+rho^mk*d;k=k+1;endx=x0;val=feval(fun,x0);endfunction [ x,val,k ] = znNewton( fun,gfun,Hess,x0,epsilon )rho=0.5;sigma=0.4;k=0;while(1)gk=feval(gfun,x0);hk=feval(Hess,x0);dk=-inv(hk)*gk;if(norm(gk)<epsilon)break;endm=0;mk=0;while(m<20)if(feval(fun,x0+rho^m*dk)<feval(fun,x0)+sigma*rho^m*gk'*dk) mk=m;break;endm=m+1;endx0=x0+rho^mk*dk;k=k+1;endx=x0;val=feval(fun,x0);endfunction [ x,val,k ] = BFGS( fun,gfun,x0,epsilon )rho=0.5;sigma=0.4;k=0;n=length(x0);bk=eye(n);while(1)gk=feval(gfun,x0);if(norm(gk)<epsilon)break;enddk=-inv(bk)*gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk);oldf=feval(fun,x0);if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=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);endk=k+1;x0=x;endval=feval(fun,x0);endfunction f = fun(x)f=5*(x(2)-x(1)^2)^2+(x(1)-1)^2;endfunction g=gfun(x)g=[-20*x(1)*(x(2)-x(1)^2)+2*(x(1)-1),10*(x(2)-x(1)^2)]'; endfunction He = Hess( x )n=length(x);He=zeros(n,n);He=[60*x(1)^2-20*x(2)+2,-20*x(1);-20*x(1), 10 ];end（2）代码运行结果如下：解：（1）MATLAB程序如下：clc;close all;clear all;x0=[1,0]';epsilon=1e-5;[ x,mu,lambda,output ] = multphr( 'f1','h1','g1','df1','dh1','dg1',x0,epsilon )function [ x,mu,lambda,output ] = multphr( fun,hf,gf,dfun,dhf,dgf,x0,epsilon )%MULTPHR Summary of this function goes here% Detailed explanation goes heresigma=2.0;eta=2.0;theta=0.8;k=0;ink=0;x=x0;he=feval(hf,x);gi=feval(gf,x);n=length(x);l=length(he);m=length(gi);%initial of multi-vectormu=0.1*ones(l,1);lambda=0.1*ones(m,1);btak=10;btaold=10;while(btak>epsilon)[x,ival,ik]=BFGS('mpsi','dmpsi',x0,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma);ink=ink+ik;he=feval(hf,x);gi=feval(gf,x);btak=0;for(i=1:l)btak=btak+he(i)^2;endfor(i=1:m)temp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);if(btak>epsilon)if(k>=2 & btak>theta*btaold)sigma=eta*sigma;endfor(i=1:l)mu(i)=mu(i)-sigma*he(i);endfor(i=1:m)lambda(i)=max(0,lambda(i)-sigma*gi(i));endendk=k+1;btaold=btak;x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=ink;output.bta=btak;endfunction [ x,val,k ] = BFGS( fun,gfun,x0,varargin )rho=0.5;epsilon=1e-5;sigma=0.4;k=0;n=length(x0);bk=eye(n);while(1)gk=feval(gfun,x0,varargin{:});if(norm(gk)<epsilon)break;enddk=-inv(bk)*gk;m=0;mk=0;while(m<20)newf=feval(fun,x0+rho^m*dk,varargin{:});oldf=feval(fun,x0,varargin{:});if(newf<oldf+sigma*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(gfun,x,varargin{:})-gk;if(yk'*sk>0)bk=bk-(bk*sk*sk'*bk)/(sk'*bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(fun,x0,varargin{:});endfunction psi = mpsi( x,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma ) %MPSI Summary of this function goes here% Detailed explanation goes heref=feval(fun,x);he=feval(hf,x);gi=feval(gf,x);l=length(he);m=length(gi);psi=f;s1=0;for(i=1:l)psi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*sigma*s1;s2=0;for(i=1:m)s3=max(0,lambda(i)-sigma*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2*sigma);endfunction dpsi = dmpsi( x,epsilon,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma ) %DMPSI Summary of this function goes here% Detailed explanation goes heredpsi=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)dpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);endfor(i=1:m)dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);endendfunction f = f1( x )f=4*x(1)-x(2)^2-12;endfunction gi = g1( x )gi=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;%unequation constrainendfunction he = h1( x )he=25-x(1)^2-x(2)^2;%equation constrainendfunction g = df1( x )g=[4,-2*x(2)]';endfunction dgi = dg1( x )dgi=[-2*x(1)+10,-2*x(2)+10]';endfunction dhe = dh1( x )dhe=[-2*x(1),-2*x(2)]';end（2）代码运行结果如下：解：（1）MATLAB代码如下：clc;close all;clear all;H=[2 0;0 2];c=[-2 -5]';Ae=[];be=[];Ai=[1 -2;-1 -2;-1 2;1 0;0 1];bi=[-2 -6 -2 0 0]';x0=[0 0]';epsilon=1e-9;[ x,lamk,exitflag,output ] = qpact( H,c,Ae,be,Ai,bi,x0,epsilon )function [ x,lamk,exitflag,output ] = qpact( H,c,Ae,be,Ai,bi,x0,epsilon ) %QPACT Summary of this function goes here% Detailed explanation goes hereerr=1e-6;k=0;x=x0;n=length(x);kmax=1e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for(i=1:ni)if(Ai(i,:)*x>bi(i)+epsilon)index(i)=0;endendwhile(k<=kmax)Aee=[];if(ne>0)Aee=Ae;endfor(j=1:ni)if(index(j)>0)Aee=[Aee;Ai(j,:)];endendgk=H*x+c;[m1,n1]=size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1));if(norm(dk)<=err)y=0;if(length(lamk)>ne)[y,jk]=min(lamk(ne+1:length(lamk)));endif(y>=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk)index(i)=0;break;endendendk=k+1;elseexitflag=1;alpha=1;tm=1;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk<0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1<tm)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1)index(ti)=1;endendif(exitflag==0)break;end %updata the setk=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k;endfunction [ x,lambda ] = qsubp( H,c,Ae,be )%QSUBP Summary of this function goes here % Detailed explanation goes hereginvH=pinv(H);[m,n]=size(Ae);if(m>0)rb=Ae*ginvH*c+be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;endend（2）代码运行结果如下：解：（1）MATLAB代码如下：clc;close all;clear all;x0=[0.5 0.2]';mu0=[ ]';lam0=[0 0 0 0]';epsilon=1e-6;[ x,mu,lam,val,k ] = sqpm( x0,mu0,lam0,epsilon)function [ x,mu,lam,val,k ] = sqpm( x0,mu0,lam0,epsilon ) %SQPM Summary of this function goes here% Detailed explanation goes heren=length(x0);l=length(mu0);m=length(lam0);rho=0.5;eta=0.1;B0=eye(n);x=x0;mu=mu0;lam=lam0;Bk=B0;sigma=0.8;[hk,gk]=cons(x);dfk=df1(x);[Ae,Ai]=dcons(x);Ak=[Ae;Ai];k=0;while(1)[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk,epsilon);mp1=norm(hk,1)+norm(max(-gk,0),1);if (norm(dk,1)<epsilon) & (mp1<1e-5)break;enddeta=0.05;tau=max(norm(mu,inf),norm(lam,inf));if(sigma*(tau+deta)<1)sigma=sigma;elsesigma=1.0/(tau+2*deta);endim=0;while(im<=20)if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)<eta*rho^im*dphi1(x,sigma,dk)) mk=im;break;endim=im+1;if(im==20)mk=10;endendalpha=rho^mk;x1=x+alpha*dk;[hk,gk]=cons(x1);dfk=df1(x1);[Ae,Ai]=dcons(x1);Ak=[Ae;Ai];lamu=pinv(Ak)'*dfk;if(l>0 & m>0)mu=lamu(1:l);lam=lamu(l+1:l+m);endif(l==0)mu=[];lam=lamu;endif(m==0)mu=lamu;lam=[];endsk=alpha*dk;yk=dlax(x1,mu,lam)-dlax(x,mu,lam);if(sk'*yk>0.2*sk'*Bk*sk)theta=1;elsetheta=0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);endzk=theta*yk+(1-theta)*Bk*sk;Bk=Bk+zk*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);x=x1;k=k+1;endval=f1(x);endfunction [ d,mu,lam,val,k ] = qpsubp( dfk,Bk,Ae,hk,Ai,gk,epsilon ) %QPSUBP Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);gamma=0.05;rho=0.5;sigma=0.2;ep0=0.05;mu0=0.05*zeros(l,1);lam0=0.05*zeros(m,1);d0=ones(n,1);u0=[ep0;zeros(n+l+m,1)];z0=[ep0;d0;mu0;lam0,];k=0;z=z0;ep=ep0;d=d0;mu=mu0;lam=lam0;while(1)dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);if(norm(dh)<epsilon)break;endA=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);b=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)*u0-dh;dz=pinv(A)*b;if(l>0 & m>0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);dl=dz(n+l+2:n+l+m+1);endif(l==0)de=dz(1);dd=dz(2:n+1);dl=dz(n+2:n+m+1);endif(m==0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);endi=0;while(i<=20)if(l>0 & m>0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(l==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(m==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);endif(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))mk=i;break;endi=i+1;if(i==20)mk=10;endendalpha=rho^mk;if(l>0 & m>0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;lam=lam+alpha*dl;endif(l==0)ep=ep+alpha*de;d=d+alpha*dd;lam=lam+alpha*dl;endif(m==0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;endk=k+1;endendfunction dh = dah( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk )%DAH Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);dh=zeros(n+l+m+1,1);dh(1)=ep;if(l>0 & m>0)dh(2:n+1)=Bk*d-Ae'*mu-Ai'*lam+dfk;dh(n+2:n+l+1)=hk+Ae*d;for(i=1:m)dh(n+l+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(l==0)dh(2:n+1)=Bk*d-Ai'*lam+dfk;for(i=1:m)dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(m==0)dh(2:n+1)=Bk*d-Ae'*mu+dfk;dh(n+2:n+l+1)=hk+Ae*d;enddh=dh(:);endfunction bet = beta( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma )%BETA Summary of this function goes here% Detailed explanation goes heredh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);bet=gamma*norm(dh)*min(1,norm(dh));endfunction [ h,g ] = cons( x )%CONS Summary of this function goes here% Detailed explanation goes hereh=[ ];g=[-x(1)^2+6*x(1)-4*x(2)+11,x(1)*x(2)-3*x(2)-exp(x(1)-1)+1,x(1),x(2)]'; endfunction [ dh,dg ] = dcons( x )%DCONS Summary of this function goes here% Detailed explanation goes heredh=[ ];dg=[-2*x(1)+6,-4;x(2)-exp(x(1)-1),x(1)-1;1,0;0,1];endfunction [ dd1,dd2,v1 ] = ddv( ep,d,lam,Ai,gk)%DDV Summary of this function goes here% Detailed explanation goes herem=length(gk);dd1=zeros(m,m);dd2=zeros(m,m);v1=zeros(m,1);for(i=1:m)fm=sqrt(lam(i)^2+(gk(i)+Ai(i,:)*d)^2+2*ep^2);dd1(i,i)=1-lam(i)/fm;dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm;v1(i)=-2*ep/fm;endendfunction df = df1( x )%DF1 Summary of this function goes here% Detailed explanation goes heredf=[2*x(1)-16,2*x(2)-10]';endfunction dl = dlax( x,mu,lam )%DLAX Summary of this function goes here% Detailed explanation goes heredf=df1(x);[Ae,Ai]=dcons(x);[m1,m2]=size(Ai);[l1,l2]=size(Ae);if(l1==0)dl=df-Ai'*lam;endif(m1==0)dl=df-Ae'*mu;endif(l1>0 & m1>0)dl=df-Ae'*mu-Ai'*lam;endendfunction f = f1( x )%F1 Summary of this function goes here% Detailed explanation goes heref=x(1)^2+x(2)^2-16*x(1)-10*x(2);endfunction p = phi( ep,a,b )%PHI Summary of this function goes here% Detailed explanation goes herep=a+b-sqrt(a^2+b^2+2*ep^2);endfunction p = phi1( x,sigma )%PHI1 Summary of this function goes here% Detailed explanation goes heref=f1(x);[h,g]=cons(x);gn=max(-g,0);l0=length(h);m0=length(g);if(l0==0)p=f+1.0/sigma*norm(gn,1);endif(m0==0)p=f+1.0/sigma*norm(h,1);endif(l0>0 & m0>0)p=f+1.0/sigma*(norm(h,1)+norm(gn,1)); endendfunction dp = dphi1( x,sigma,d )%DPHI1 Summary of this function goes here% Detailed explanation goes heredf=df1(x);[h,g]=cons(x);gn=max(-g,0);l0=length(h);m0=length(g);if(l0==0)dp=df'*d-1.0/sigma*norm(gn,1);endif(m0==0)dp=df'*d-1.0/sigma*norm(h,1);endif(l0>0 & m0>0)dp=df'*d-1.0/sigma*(norm(h,1)+norm(gn,1)); endendfunction A = JacobiH( ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk )%JACOBIH Summary of this function goes here% Detailed explanation goes heren=length(dfk);l=length(hk);m=length(gk);A=zeros(n+l+m+1,n+l+m+1);[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk);if(l>0 & m>0)A=[1, zeros(1,n), zeros(1,l), zeros(1,m);zeros(n,1), Bk, -Ae', -Ai';zeros(l,1), Ae, zeros(l,l), zeros(l,m);v1, dd2*Ai, zeros(m,l), dd1]; endif(l==0)A=[1, zeros(1,n), zeros(1,m);zeros(n,1), Bk, -Ai';v1, dd2*Ai, dd1];endif(m==0)A=[1, zeros(1,n), zeros(1,l);zeros(n,1), Bk, -Ae';zeros(l,1), Ae, zeros(l,l)];endend（2）代码运行结果如下：解：（1）MATLAB代码如下：clc;close all;clear all;f=[1 1 1 1 1 1 1 1];A=[1 0 0 0.5;0 1 0.2 0.3;0 0.1 1 0.2];Aeq=[A,-A];beq=[-1 0.2 1]';lb=zeros(8,1);x=linprog(f,[],[],Aeq,beq,lb);x=[x(1)-x(5) x(2)-x(6) x(3)-x(7) x(4)-x(8)]' （2）代码运行结果如下：。

大连理工大学概率上机作业

大连理工大学概率上机作业————————————————————————————————作者: ————————————————————————————————日期：ﻩ第一次上机作业1.利用Matlａb自带命令产生100０个均匀随机变量服从Ｕ(０，1）。

>＞unｉfrnｄ(0，1，2０,5０)aｎs＝Ｃoｌｕｍns 1 tｈrouｇh 1０0.81４７0．6５５7０．4387 ０.75130.３517 0.1６22０．106７0.85３00.7８0３0．5４700．9058 0.03570.3８１6 0.25510.8308 0.7９４3 0．９619 0.6221 ０．3897 0.２963０.1270 0.849１０.7６55 0.５06０0．585３0.３1１2 0.０046 0.３5100.241７0.74４70.９134 0.93４00.7952０.6991 0.5４97 ０．52８5 ０.７749 0.51３2 0．４039 0.18900.６3２4 0.６７87 0.1869 0.8909 0.91７2 0.1656 0．８173 0.4０１８０.0965 0.6868０.09750.757７0.48９８0.95９3 0.28５８0．6020 ０．86870．0760０.13２0 0.１8350.2785 ０.743１0.４4５60.５472 0．7５7２0.2６3０0．084４0.23990.９42１0．368５０.５469 0.3９220.64630．138６0.7５37０.65４1 0.３998 ０.1233 0．95６1 ０.６２560．9575 0．6555 ０.7０９4 0．1493 0.3８04 ０.6８９2 0．25990.１83９0.5７５2 0.78０20.９６49 0.1712 0.７５470．257５０.5６7８0．74８2 0．80010.2４00０．05９８0.081１0.1５760.7060 0.27６0 0.８４07 ０.0７59 0.４505 0.４３14 0.4１730.2348 0.929４０.970６0．0３1８０.679７０.２543 0.０5400.083８0．9106 0.０497 0．35３20.７７57０.９572 0．2769 ０.655１0.８１43 0.5308 ０．22９0０.18１8０．９0２7 0.8212 ０．4８６８0.4854 0．0462 0.１６26 ０.2435 0.7792 0.91３3 0.２638 0.94480.01540.４35９0.8003 ０.０971 0.１19０0．92９３0．9340 ０．1524 ０．1455 0．４909 0.０４30 0.4468０.1419 0.8２3５0.49８4 ０．35０0 0.1２9９0.８2５８0.１36１0.4８93 0.1690 0.３０63０.4２18 0.6948０.9597 0.１966０.５６8８０.5383 0.8693 ０.3３７7 0.6491 0.508５0.9157 0.3171０.３4０4 0.2５１1 0．46９4 0.996１0.579７0.9０0１0．7３17 0.51０80.7922 0．９50２0.5853 0.61600.０1190.07８20．54９９０.３692 0.6４77 ０．81760.9５950.0344 ０.2２38 ０.4733 0.３3７1 0.4427０.1450 ０．1１120.4509 0.794８Colｕmns １1 throｕｇh 200.6443 0.31１10．08５5 ０.0377 0.03050.0596 ０．173４0.９51６0.０3２6 0．25180.３7８6 ０.9２3４0.262５0．８8５2 ０.74４1０.68200．39０9 0.920３0.５6120.29040.８116 0.4302 0.8010 ０.9１3３0．500０0．04２4 0.831４0.0527０.８819 0.6171０．５328 0．18480．0292 ０.79620.４7990.07１４0.8034 ０.７379 0.6６9２0.265３0.3５07 0.9049 ０.92８9 0.0987 0.９047０．５21６0.06050.26９10.１９0４0.824４0.9390 0．9７97 0．７303 0.26１９０.609９0.０９6７0．39９3０.４228０．3689 0．982７0.８759 ０.4389 0.4886 ０.３354 0.6177 0.８18１0．5269 0.54790.4607 ０.730２0.55020.11１1 ０.5785 0．67９7 ０.8５94 0．81750．4１680.94２７0.9816 0.３4390.６２２５0.25８1 0.23７30.1366 0.8055 0.７２24 0.６5690.4177 0.1５6４0.58４1０．５870 0.4087 0.458８０．7212 0.5７6７0.1499０．6２80 ０.９83１０．8555 0.1０780．20７70.5949 ０.963１0.106８0.18２９０．６5９6 0.2920 0.3０１5 0.6448 0.９06３0．3012 ０．2622 0．5４680.6５38 ０.２３99０.５186 0.4３1７0.７0１1 0.３763 0.8797０.4709 0.6０2８0．52110.494２0．88６5 0.９７3００.0155 0.６663 0.190９0．817８0.230５0．7112 0.231６0.７79１０.02８70.6490 ０.9841 0．5３91 0．4283 ０.260７0．844３０.2217 ０.488９０.７１50 0．4899 ０.8003 0.1６72 0.6９８１0.４８20 0.59４4０．19４8 0.1174 ０.６241 0．90370．１6790．4５38 0.106２0.666５0．1206 0.02２５0.2２5９0.2９67 0.67９1 0.89０9 0．97８7 0.４３2４0.3724 0.1781 0.５8９５0．425３０.1707 ０．３188 ０.３９55 0.3342 ０.7127 ０.82５3 ０.1981 ０.１280 0．2262 0.31２7０.2277 0．4242 0.3674 0.６987 0.5００5 0.０835 ０.48970.99９1 ０.３846 0.１6150.43５7 0.5０７9 ０.988００.197８0.471１0.1332 ０.3395０.171１0．5830 0.１７88Coluｍnｓ２1thrｏugh 300.4２２９０.７7８8 0.２54８0．17５9 0．６476 0.5８22 ０．4046 0.3477 0．８2１７0.5１４4０．0942 0．42３５0.22４0 ０．7２18 ０．67900.5４070.44８4 ０．1500 ０.42９９0.884３0.598５0．09080.6６７80．473５0.63５8 0．86990．3658 0．5861 ０.88780．58８00.47090.2６65 0.8４４4 0.１527 0.945２0.26480.７63５0.2621 ０．3912 0.1548０.69５9 0.1５3７０.３44５0．34１1０.２089 0.３１８1 ０.６27９0.044５0．7691 0.１999０.６999０．28１0 ０.780５０.６07４0.70930.119２0．77２0 ０.7５49 0．39６8 ０．40７0０.638５0.440１０.6753 0.19170.２３6２0.9398 0.932９0．2428 ０．8085 ０.74870.033６0.5２7１0.０0６7 0.73840．11940.64560.9727 0.442４０．7551 0.8２5６０.０６88 0.４５7４０.6022 0.2４２80.６073 0.4795 0.192００.687８０.3７740.79000.31９6 0.875４０．386８0.9１74 ０.４５０1 0．639３0.13８9０．359２0.2160 0.3185０.53090.5181 ０.9１600.２6９1 ０.4５87０.54４7 0.69６３０.7３63 0.7904 0.53410.6５４4 0.9436 ０.00１2 0．76５5 ０.6619 0.64730.09３8 0.3947 0.９４9３0.090０0.４0７6 0.6377 0.４6２4 0.1887 0．７70３０．5４39 0.52５4 0．6834 0.３2７60．１11７0.8200 0.9５770.4２4３0.287５0.3502 0.７2１0 0.５3030．7040 0．6713 0.13６３0.７18４0.24０７0.46090.0911 0.66２0 0.5225 ０.8611 0.4423 0.4３8６0.67８70．9６8６0.6７61 0.7７0２0.５762 0．41６20.９937 0.４849 0.０196 ０．8３35 0.49520.53１3 0.２89１0.3２25 0．６83４０.８４１９0.２18７０．39３50.33０9 ０.７689 ０．１8970.3251 0.671８0.7847 ０.５466 ０.８32９０．１０58 ０.67１４0．4243 ０.１67３0.495００.1０5６0．6９51０.４7１4 0.4257 ０.２56４０．109７0．７４１3 ０．２7０3 0.８620 0.14７60．61１0 ０．0６８０0.0358０．64４4 ０.６13５0．063６0．520１0．1971 0.9899 0．0５５0Cｏlumｎs 31 thｒｏugh 400．8５０７0．7３8６0.55２３0．１２３９０．737８0.55９0 ０.1781 0.894９0．6３11 ０.6９25０．５60６0.5８６00.62９９０.4904 0．0６3４０.8541 ０.３5９6 0.0715０.0８9９0.55６70.９２96 0.246７0.０320０.8530 0.860４0.3４79 0.０567 0.２425 ０.080９０.3９65０.696７0.６664 ０．6１4７0.873９0.93４４0.4460 0.5２19 ０.05３8 0.777２0.0６1６０．５８２80.０8３50.３624 0.2703 0.9844 0.0542 ０.３358 0．441７0．９０51 0．78020.８1５4 0．62６00.０４９50．2085 0.８589 0.1771０．１７57０.0１330．53380.3３760.87９0 0.６６09 0．４8９6 0.５650 0．7856 0.6６28 0．２08９０.８９7２0．109２0.６0790.988９0.7298 0.19２50．6403 0．51340.３３0８０．90520．19６7 0.8２5８0.7４13０.０0050.89０80.12310.４17０0.17760.8985 ０.67５4 0.09３４0.3３81 0.１048０．8６5４０.98230.２05５0.2060 0.39860.1182 0.4６８5 0．3０74 0．２９40 0.1２79０.６126 0.76９０0．１４6５0.9479０.１３3９0．9884 0.9１2１0.4５61 0.74６3 0.54９50.99０0０.５81４0.1891 0.0８21 ０．030９０.54000.1０４００.１017 0.0103 0.4８520.5277 0.９２８3 ０.0４２7 0．１０57０．9３９1 0.７０69 ０．74５５0.9954 0.０484 0.89０50.４795 0.5801 0．6352 0.142０0.３０1３0.9９95 ０．7363 0.3321 0.６6７９0.7990０.8013 0.０170 0.2819 0．１665 0.295５0.287８0.5619 0.２973 0．6035 0.7343０.2278 0.1209 ０.５386 0.62１00.3３２9 ０.41４5 0.18４20.0620０.526１０.05１30.49８1 0．86２7 0.6952 ０．5７3７0.4671 0．4648 ０.5972 0.２982 ０.7２97０.07２90.9009０.4843 0.４９91 0.0521 0.648２0.7640 ０.2９99 0.0464 0.7073０．088５０.57４70.８44９0.53５８0．9312 0.02５2 ０.81820.134１0.505４０.78１4 0.7984０.８452 0.20940．44５2 0.7287 0.8４22 0.１０0２0.２12６0．76１4０.288０0.9４30Coｌuｍns 41 ｔhｒough 5０0．６8３7 ０.789４0．1123 ０.６733 0．09８６0.9879 0．5９７5 0．７5９3０.809２0.7５190.1３２1 0.367７0.７8４４0.4２9６0.142００.1704 0.33５3 ０.7４０6 0.7486 0.２2870.7２27 ０．2060 0.2916 0.4517 0．１６83 0.2578 ０.２9９2 0．74370.120２０．０6420．110４0.0８67 0.6０350.6０99 0.1９620.3968 ０.4526 ０.１0590.５２50 0．7６730.11７５０.77１９0.9644 0．0594 0.３17５0.0740 0.42２6 ０.68１６０.３２５8 0.6712０.6407 0.2057 0.４3250．3１５8 0.３1640．６841 0.3５9６0.46330.５４64 0．７15２0.3２88 ０.３8８３0.6948 0.7727 0.2176 0.4024 0.5583 0.21２２0.39８9 0.64210.65380.55１8 0.7５810.６９64 0.2５１00.９8２8 0.74250.0985０.4１51 0.41900.74９1 0.2290 0.43２6 0.１２5３0．8929 0.4022 0.42４3 ０．8２360.18０7 0．39080.５８320.６４１9 0.６5550．１３02 0.703２0.６２０7 0．4294 0.17５0 0.25５4 0.81６1０.7400０.48４５０.１0980．0９24 0.５557 0.1544 0．1２49 0.1636 0.02０5 0.317４0.23４8 0.15180．933８0.０078０．1844 0．3813 0.0244 0．66６００.9237 0.８１45０.7350 0.78１9０．18７5 0．４2３１0.２１200．1611 0.２902 ０.８944 0．６５3７0.７89１0.97060.10060．2662 0.655６0.0７7３０.7581０．31７5 ０.5166 0.9３２６0.８5２3０.８669 ０.29410．７978 0．7229 ０.91３８0.８７11 0.６5370.70２7０．1６３5 0.５056０.08620.2３7４0.48760.5312０.7０6７0.３5080.9569 0.1536 0.９２11 0.63570．３664 0.５309 0.769０0.10880.55７8 ０.68550.9３57 0．95３50.７94７0.９509０．３6９2 0.0915 0.39６0 0.6318０.３13４０.2９４1 0.４５79 ０.540９0.577４0.４44０0.68５0 0．405３0.2729０.126５0.1662 ０.530６０.240５0.6797０.4400 ０.06000.5９７9 ０.１04８0.0３72 0.13４3 0．６225 0.832４０.76390.03660．257６0.86672.参考课本综合例题2.5.４和2．5.5中的方法，模拟产生1０00个随机变量，使其服从参数为2的指数分布,进而计算这1000个随机数的均值和方差。

大连理工大学数据结构(一)上机作业答案——张老师

#define OVERFLOW -2
typedef int ElemType;
typedef int Status;
#define LIST_INIT_SIZE 100
#define LISTTINCREMENT 10
typedef struct{
ElemType *elem;
int length;
for(i=0;i<n;i++)
{
printf("input elem:");
scanf("%d",&L.elem[i]);
}
}
//输出顺序表中的元素
void ListOutput_Sq(SqList L){
int i,n;
n=L.length;
for(i=0;i<n;i++)
printf("%2d",L.elem[i]);
else
return -2;
}//优先级比较
void transform(char suffix[], char exp[])
{
char *p, ch, c;
p = exp;
ch = *p;
SqStack S;
InitStack (S);
Push (S, '#');
while(!StackEmpty(S))
}
}
}
void DeleteList(LinkList &L,ElemType mink,ElemType maxk){
LNode *p=L,*q;
while(p->next&&p->next->data<=mink)

大连理工大学优化方法上机大作业

学院：专业：班级：学号：姓名：上机大作业1：1.最速下降法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = steepest(x0,eps) gk = grad(x0);res = norm(gk);k = 0;while res > eps && k<=1000dk = -gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;gk = grad(xk);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x=steepest(x0,eps)2.牛顿法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad2(x)g = zeros(2,2);g(1,1)=2+400*(3*x(1)^2-x(2));g(1,2)=-400*x(1);g(2,1)=-400*x(1);g(2,2)=200;endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = newton(x0,eps)gk = grad(x0);bk = [grad2(x0)]^(-1);res = norm(gk);k = 0;while res > eps && k<=1000dk=-bk*gk;xk=x0+dk;k = k+1;x0 = xk;gk = grad(xk);bk = [grad2(xk)]^(-1);res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x1=newton(x0,eps)--The 1-th iter, the residual is--The 2-th iter, the residual isx1 =法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star = bfgs(x0,eps) g0 = grad(x0);gk=g0;res = norm(gk);Hk=eye(2);k = 0;while res > eps && k<=1000dk = -Hk*gk;ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slopeak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;fa0=xk-x0;x0 = xk;go=gk;gk = grad(xk);y0=gk-g0;Hk=((eye(2)-fa0*(y0)')/((fa0)'*(y0)))*((eye(2)-(y0)*(fa0)')/((fa0)'*(y0)))+(fa0*(fa 0)')/((fa0)'*(y0));res = norm(gk);fprintf('--The %d-th iter, the residual is %f\n',k,res);endx_star = xk;End>> clear>> x0=[0,0]';>> eps=1e-4;>> x=bfgs(x0,eps)4.共轭梯度法：function f = fun(x)f = (1-x(1))^2 + 100*(x(2)-x(1)^2)^2; endfunction g = grad(x)g = zeros(2,1);g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); g(2) = 200*(x(2)-x(1)^2);endfunction x_star =CG(x0,eps) gk = grad(x0);res = norm(gk);k = 0;dk = -gk;while res > eps && k<=1000 ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + *ak*slope ak = ak/4;xk = x0 + ak*dk;f1 = fun(xk);endk = k+1;x0 = xk;g0=gk;gk = grad(xk);res = norm(gk);p=(gk/g0)^2;dk1=dk;dk=-gk+p*dk1;fprintf('--The %d-th iter, the residual is %f\n',k,res); endx_star = xk;end>> clear>> x0=[0,0]';>> eps=1e-4;>> x=CG(x0,eps)上机大作业2：function f= obj(x)f=4*x(1)-x(2)^2-12;endfunction [h,g] =constrains(x) h=x(1)^2+x(2)^2-25;g=zeros(3,1);g(1)=-10*x(1)+x(1)^2-10*x(2)+x(2)^2+34;g(2)=-x(1);g(3)=-x(2);endfunction f=alobj(x) %拉格朗日增广函数%N_equ等式约束个数?%N_inequ不等式约束个数N_equ=1;N_inequ=3;global r_al pena;%全局变量h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式约束部分?for i=1:N_equh_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i).^2;end%不等式约束部分for i=1:N_inequh_inequ=h_inequ+pena)*(max(0,(r_al(i)+pena*g(i))).^2-r_al(i).^2); end%拉格朗日增广函数值f=obj(x)+h_equ+h_inequ;function f=compare(x)global r_al pena N_equ N_inequ;N_equ=1;N_inequ=3;h_inequ=zeros(3,1);[h,g]=constrains(x);%等式部分for i=1:1h_equ=abs(h(i));end%不等式部分for i=1:3h_inequ=abs(max(g(i),-r_al(i+1)/pena));endh1 = max(h_inequ);f= max(abs(h_equ),h1); %sqrt(h_equ+h_inequ);function [ x,fmin,k] =almain(x_al)%本程序为拉格朗日乘子算法示例算法%函数输入：% x_al:初始迭代点% r_al：初始拉格朗日乘子N-equ:等式约束个数N_inequ:不等式约束个数?%函数输出% X：最优函数点FVAL:最优函数值%============================程序开始================================ global r_al pena ; %参数（全局变量）pena=10; %惩罚系数r_al=[1,1,1,1];c_scale=2; %乘法系数乘数cta=; %下降标准系数e_al=1e-4; %误差控制范围max_itera=25;out_itera=1; %迭代次数%===========================算法迭代开始============================= while out_itera<max_iterax_al0=x_al;r_al0=r_al;%判断函数?compareFlag=compare(x_al0);%无约束的拟牛顿法BFGS[X,fmin]=fminunc(@alobj,x_al0);x_al=X; %得到新迭代点%判断停止条件?if compare(x_al)<e_aldisp('we get the opt point');breakend%c判断函数下降度?if compare(x_al)<cta*compareFlagpena=1*pena; %可以根据需要修改惩罚系数变量elsepena=min(1000,c_scale*pena); %%乘法系数最大1000disp('pena=2*pena');end%%?更新拉格朗日乘子[h,g]=constrains(x_al);for i=1:1%%等式约束部分r_al(i)= r_al0(i)+pena*h(i);endfor i=1:3%%不等式约束部分r_al(i+1)=max(0,(r_al0(i+1)+pena*g(i)));endout_itera=out_itera+1;end%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++ disp('the iteration number');k=out_itera;disp('the value of constrains'); compare(x_al)disp('the opt point');x=x_al;fmin=obj(X);>> clear>> x_al=[0,0];>> [x,fmin,k]=almain(x_al)上机大作业3： 1、>> clear alln=3; c=[-3,-1,-3]'; A=[2,1,1;1,2,3;2,2,1;-1,0,0;0,-1,0;0,0,-1];b=[2,5,6,0,0,0]';cvx_beginvariable x(n)minimize( c'*x)subject toA*x<=bcvx_endCalling SDPT3 : 6 variables, 3 equality constraints------------------------------------------------------------num. of constraints = 3dim. of linear var = 6*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0|||+01|+00|+02|+01 +00| 0:0:00| chol 1 11|||||+01|+00 +01| 0:0:01| chol 1 12|||||+00|+00 +01| 0:0:01| chol 1 13|||||+00|+00 +00| 0:0:01| chol 1 14||||||+00 +00| 0:0:01| chol 1 15||||||+00 +00| 0:0:01| chol 1 16||||||+00 +00| 0:0:01| chol 1 17||||||+00 +00| 0:0:01| chol 1 18||||||+00 +00| 0:0:01|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 8primal objective value = +00dual objective value = +00gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval):2、>> clear alln=2; c=[-2,-4]'; G=[,0;0,1]; A=[1,1;-1,0;0,-1]; b=[1,0,0]'; cvx_beginvariable x(n)minimize( x'*G*x+c'*x)subject toA*x<=bcvx_endCalling SDPT3 : 7 variables, 3 equality constraintsFor improved efficiency, SDPT3 is solving the dual problem.------------------------------------------------------------num. of constraints = 3dim. of socp var = 4, num. of socp blk = 1dim. of linear var = 3*******************************************************************SDPT3: Infeasible path-following algorithms*******************************************************************version predcorr gam expon scale_dataNT 1 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime -------------------------------------------------------------------0||||+00|+02| +01 +00| 0:0:00| chol 1 11|||||+01| +00 | 0:0:00| chol 1 12|||||+00| +00 | 0:0:00| chol 1 13|||||| | 0:0:00| chol 1 14|||||| | 0:0:00| chol 1 15|||||| | 0:0:00| chol 1 16|||||| | 0:0:00| chol 1 17|||||| | 0:0:00| chol 1 18|||||| | 0:0:00| chol 1 19|||||| | 0:0:00| chol 1 110|||||| | 0:0:00| chol 2 211|||||| | 0:0:00| chol 2 212|||||| | 0:0:00| chol 2 213|||||| | 0:0:00| chol 2 214|||||| | 0:0:00|stop: max(relative gap, infeasibilities) <------------------------------------------------------------------- number of iterations = 14primal objective value =dual objective value =gap := trace(XZ) =relative gap =actual relative gap =rel. primal infeas (scaled problem) =rel. dual " " " =rel. primal infeas (unscaled problem) = +00rel. dual " " " = +00norm(X), norm(y), norm(Z) = +00, +00, +00norm(A), norm(b), norm(C) = +00, +00, +00Total CPU time (secs) =CPU time per iteration =termination code = 0DIMACS: +00 +00-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -3。

浮式生产储油船(FPSO)设计建造研究船舶与海洋工程论文

大连理工大学博士学位论文浮式生产储油船（FPSO）设计建造研究姓名：马延德申请学位级别：博士专业：船舶与海洋结构物设计制造指导教师：王言英20061201 大连理工大学博士学位论文摘要ＦＰＳＯ是ＦｌｏａｔｉｎｇＰｒｏｄｕｃｔｉｏｎＳｔｏｒａｇｅａｎｄＯｆｆｉｏａｄｉｎｇ的英文缩写，即浮式生产储油卸油装置，习惯上我们称为浮式生产储油船。

它是集生产、储油、外输、生活、动力于一体的多功能采油设施，是海洋石油开发中非常重要、也是最有应用前景的装备之一。

国外ＦＰＳＯ的设计建造始于二十世纪七十年代，经过多年的发展，国夕｝公司对于ＦＰＳＯ关键技术的研究日趋成熟。

国内对ＦＰＳＯ设计建造的研究起步相对较晚，虽然相关单位也对ＦＰＳＯ设计建造的部分技术进行了多年研究并取得了一定成果，但是对ＦＰＳＯ的总体设计和ＦＰＳＯ建造过程中的特别之处尚缺乏系统分析，对于ＦＰＳＯ设计中关键技术之一“系泊系统的设计”尚缺乏理论研究。

针对这些不足之处，本论文对影响浮式生产储油设施ＦＰＳＯ设计建造的因素进行了综合分析，并对系泊系统的设计进行了重点研究。

基于对ＦＰＳＯ相关规范的研究，结合我国自行设计建造的１５万吨级ＰＦＳＯ的实际经验，采用了层式分析法和模糊评判法，对该船型的方案论证、总布置、可靠性评估等方面做了详细分析，总结得出ＦＰＳ０的总体设计框架与原则，设计单位可以根据该原则对承接的ＦＰＳ０进行初步设计。

通过将浮式生产储油设旋ＦＰＳＯ和普通油船进行对比分析，首次给出一系列反映两者区别的直观表格，并在此基础上归纳了设计建造ＦＰＳＯ所必须考虑的影响因素，可供船厂建造ＦＰｓＯ过程中结合已有的油船建造经验进行参考分析。

应用流体力学理论和数学工具开发了ＦＰＳＯ环境载荷的计算方法，对系泊系统的设计和模型实验进行了研究。

如何确定外部环境载荷，是本部分研究的重点。

采用线性化Ｗｅｉｂｕｌｌ概率密度函数分析得到设计波参数，并根据三维源汇分布方法建立浮体运动与波浪荷载计算方法，完成了一浮式生产储油船（ＦＰＳＯ）在波浪中的运动响应和船体表面水动力压力分布以及总体荷载的概率特性的计算，并以此为基础对ＦＰＳＯ的系泊系统的设计和模型实验开展研究。

2016年大连理工大学优化方法上机大作业

Removed_大连理工大学工科数学分析上机作业

教学大纲-大连理工大学教务处

(整理)大连理工大学--优化作业----程序.

优化方法大作业1

大连理工优化方法大作业MATLAB编程

大连理工大学结构优化复习总结

大连理工大学优化方法上机作业

大连理工优化方法大作业MATLAB编程

大连理工高等物理化学大作业

大连理工大学庞丽萍最优化方法MATLAB程序

大连理工大学 秋季优化方法大作业

最优化方法大作业

优化方法课程大作业

大连理工大学概率上机作业

大连理工大学数据结构(一)上机作业答案——张老师

大连理工大学优化方法上机大作业

浮式生产储油船(FPSO)设计建造研究 船舶与海洋工程论文

大连理工大学庞丽萍最优化方法MATLAB程序

大连理工大学秋季优化方法大作业

浮式生产储油船(FPSO)设计建造研究船舶与海洋工程论文