大连理工大学2016年秋季优化方法大作业
2016年大连理工大学优化方法上机大作业
2016年理工大学优化方法上机大作业学院：专业：班级：学号：：上机大作业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 + 0.1*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 447.213595--The 2-th iter, the residual is 0.000000x1 =1.00001.00003.BFGS法：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 + 0.1*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*(fa0)')/((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<=1000ak =1; f0 = fun(x0);f1 = fun(x0+ak*dk);slope = dot(gk,dk);while f1 > f0 + 0.1*ak*slopeak = 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+(0.5/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=0.5; %下降标准系数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 4.0: 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 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|1.1e+01|5.1e+00|6.0e+02|-7.000000e+01 0.000000e+00| 0:0:00| chol 1 11|0.912|1.000|9.4e-01|4.6e-02|6.5e+01|-5.606627e+00 -2.967567e+01| 0:0:01| chol 1 12|1.000|1.000|1.3e-07|4.6e-03|8.5e+00|-2.723981e+00 -1.113509e+01| 0:0:01| chol 1 13|1.000|0.961|2.3e-08|6.2e-04|1.8e+00|-4.348354e+00 -6.122853e+00| 0:0:01| chol 1 14|0.881|1.000|2.2e-08|4.6e-05|3.7e-01|-5.255152e+00 -5.622375e+00| 0:0:01| chol 1 15|0.995|0.962|1.6e-09|6.2e-06|1.5e-02|-5.394782e+00 -5.409213e+00| 0:0:01| chol 1 16|0.989|0.989|2.7e-10|5.2e-07|1.7e-04|-5.399940e+00 -5.400100e+00| 0:0:01| chol 1 17|0.989|0.989|5.3e-11|5.8e-09|1.8e-06|-5.399999e+00 -5.400001e+00| 0:0:01| chol 1 18|1.000|0.994|2.8e-13|4.3e-11|2.7e-08|-5.400000e+00 -5.400000e+00| 0:0:01| stop: max(relative gap, infeasibilities) < 1.49e-08-------------------------------------------------------------------number of iterations = 8primal objective value = -5.39999999e+00dual objective value = -5.40000002e+00gap := trace(XZ) = 2.66e-08relative gap = 2.26e-09actual relative gap = 2.21e-09rel. primal infeas (scaled problem) = 2.77e-13rel. dual " " " = 4.31e-11rel. primal infeas (unscaled problem) = 0.00e+00rel. dual " " " = 0.00e+00norm(X), norm(y), norm(Z) = 4.3e+00, 1.3e+00, 1.9e+00norm(A), norm(b), norm(C) = 6.7e+00, 9.1e+00, 5.4e+00Total CPU time (secs) = 0.71CPU time per iteration = 0.09termination code = 0DIMACS: 3.6e-13 0.0e+00 5.8e-11 0.0e+00 2.2e-09 2.3e-09-------------------------------------------------------------------------------------------------------------------------------Status: SolvedOptimal value (cvx_optval): -5.42、>> clear alln=2; c=[-2,-4]'; G=[0.5,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 4.0: 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 0.000 1 0it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime-------------------------------------------------------------------0|0.000|0.000|8.0e-01|6.5e+00|3.1e+02| 1.000000e+01 0.000000e+00| 0:0:00| chol 1 1 1|1.000|0.987|4.3e-07|1.5e-01|1.6e+01| 9.043148e+00 -2.714056e-01| 0:0:00| chol 1 1 2|1.000|1.000|2.6e-07|7.6e-03|1.4e+00| 1.234938e+00 -5.011630e-02| 0:0:00| chol 1 1 3|1.000|1.000|2.4e-07|7.6e-04|3.0e-01| 4.166959e-01 1.181563e-01| 0:0:00| chol 1 1 4|0.892|0.877|6.4e-08|1.6e-04|5.2e-02| 2.773022e-01 2.265122e-01| 0:0:00| chol 1 1 5|1.000|1.000|1.0e-08|7.6e-06|1.5e-02| 2.579468e-01 2.427203e-01| 0:0:00| chol 1 1 6|0.905|0.904|3.1e-09|1.4e-06|2.3e-03| 2.511936e-01 2.488619e-01| 0:0:00| chol 1 1 7|1.000|1.000|6.1e-09|7.7e-08|6.6e-04| 2.503336e-01 2.496718e-01| 0:0:00| chol 1 1 8|0.903|0.903|1.8e-09|1.5e-08|1.0e-04| 2.500507e-01 2.499497e-01| 0:0:00| chol 1 19|1.000|1.000|4.9e-10|3.5e-10|2.9e-05| 2.500143e-01 2.499857e-01| 0:0:00| chol 1 1 10|0.904|0.904|4.7e-11|1.3e-10|4.4e-06| 2.500022e-01 2.499978e-01| 0:0:00| chol 2 2 11|1.000|1.000|2.3e-12|9.4e-12|1.2e-06| 2.500006e-01 2.499994e-01| 0:0:00| chol 2 2 12|1.000|1.000|4.7e-13|1.0e-12|1.8e-07| 2.500001e-01 2.499999e-01| 0:0:00| chol 2 2 13|1.000|1.000|2.0e-12|1.0e-12|4.2e-08| 2.500000e-01 2.500000e-01| 0:0:00| chol 2 2 14|1.000|1.000|2.6e-12|1.0e-12|7.3e-09| 2.500000e-01 2.500000e-01| 0:0:00|stop: max(relative gap, infeasibilities) < 1.49e-08-------------------------------------------------------------------number of iterations = 14primal objective value = 2.50000004e-01dual objective value = 2.49999996e-01gap := trace(XZ) = 7.29e-09relative gap = 4.86e-09actual relative gap = 4.86e-09rel. primal infeas (scaled problem) = 2.63e-12rel. dual " " " = 1.00e-12rel. primal infeas (unscaled problem) = 0.00e+00rel. dual " " " = 0.00e+00norm(X), norm(y), norm(Z) = 3.2e+00, 1.5e+00, 1.9e+00norm(A), norm(b), norm(C) = 3.9e+00, 4.2e+00, 2.6e+00Total CPU time (secs) = 0.36CPU time per iteration = 0.03termination code = 0DIMACS: 3.7e-12 0.0e+00 1.3e-12 0.0e+00 4.9e-09 4.9e-09------------------------------------------------------------------------------------------------------------------------------- Status: SolvedOptimal value (cvx_optval): -3。
(整理)大连理工大学--优化作业----程序.
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。
优化作业流程的三种方法
优化作业流程的三种方法下载温馨提示:该文档是我店铺精心编制而成，希望大家下载以后，能够帮助大家解决实际的问题。
文档下载后可定制随意修改，请根据实际需要进行相应的调整和使用，谢谢!并且，本店铺为大家提供各种各样类型的实用资料，如教育随笔、日记赏析、句子摘抄、古诗大全、经典美文、话题作文、工作总结、词语解析、文案摘录、其他资料等等，如想了解不同资料格式和写法，敬请关注!Download tips: This document is carefully compiled by theeditor. I hope that after you download them,they can help yousolve practical problems. The document can be customized andmodified after downloading,please adjust and use it according toactual needs, thank you!In addition, our shop provides you with various types ofpractical materials,such as educational essays, diaryappreciation,sentence excerpts,ancient poems,classic articles,topic composition,work summary,word parsing,copy excerpts,other materials and so on,want to know different data formats andwriting methods,please pay attention!优化作业流程是提高工作效率和质量的重要手段。
以下是三种常见的优化作业流程的方法：1. 流程再造：步骤：分析现有流程：对当前的作业流程进行全面的分析，包括流程的各个环节、涉及的人员、使用的工具和资源等。
大连理工大学优化方法上机作业
大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除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。
大连市2016～2017学年度第一学期期末考试试卷及答案
大连市2016～2017学年度第一学期期末考试试卷及答案大连市2016～2017学年度第一学期期末考试试卷高二物理（理科）命题人：侯贵民陈晶李辉校对人：侯贵民注意事项：1．请在答题纸上作答，在试卷上作答无效。
2．本试卷分第Ⅰ卷（选择题）和第Ⅱ卷（非选择题）两部分，满分100分，考试时间90分钟。
第Ⅰ卷选择题（共48分）一、选择题（本题共12小题，每小题4分，共48分。
在每小题给出的四个选项中，第1～8题只有一项符合题目要求，第9～12题有多项符合题目要求。
全部选对的得4分，选对但不全的得2分，有选错的得0分。
）1．真空中有甲、乙两个点电荷相距为r ，它们间的静电斥力为F 。
若甲的电荷量变为原来的2倍，乙的电荷量保持不变，它们间的距离变为2r ，则它们之间的静电斥力将变为（）A ．12FB ．14FC ．FD ．2F2．两个固定的等量异种电荷，在它们连线的垂直平分线上有a 、b 两点，如图所示，下列说法正确的是（）A ．a 点场强比b 点场强大B ．a 点场强比b 点场强小C ．a 点电势比b 点电势高D ．a 点电势比b 点电势低3．如图所示，平行板电容器接在电势差恒为U 的电源两端，下极板接地，现将平行板电容器的上极板竖直向上移动一小段距离，则（）A ．电容器的电容减小B ．电容器的电容增大C ．上极板带电荷量将增大D ．下极板带电荷量保持不变ab4．当导线中分别通以如图所示方向的电流，小磁针静止时N 极指向纸面外的是（）5．1930年劳伦斯制成了世界上第一台回旋加速器，其结构如图所示，加速器由两个铜质D 形盒D 1、D 2构成，其间留有空隙，下列说法正确的是（）A ．粒子做圆周运动的周期随运动半径的增大而增大B ．粒子做圆周运动的周期与高频交变电压的周期相同C ．粒子从磁场中获得能量D ．粒子所获得的最大能量与高频交变电压有关6．如图所示四幅图中，匀强磁场磁感应强度大小都为B ，导体中的电流大小都为I ，导体的长度都为L ，A 、B 、C 三个图中的导体都与纸面平行，则四幅图中导体所受安培力大小不等于BIL 的是（）7．如图甲所示，一圆形金属线圈处于匀强磁场中，磁感应强度B 随t 的变化的规律如图乙所示。
大连理工大学c语言大作业
大连理工大学c语言大作业第一篇：大连理工大学c语言大作业程序设计大作业总结报告——<东北大馅饺子馆>的点餐/帐目信息管理系统选题意义；餐厅账目繁多，通过本系统可以实现餐厅管理的自主化。
更加适应这个信息化的社会。
通过对账目、订单的管理与排序。
也能使餐厅管理者更加直观地看出产品之间的优劣与受众。
使其能够更好的调动资源，达到餐厅的快速发展目的。
设计方案；1）任务分析该系统应包括两大界面—用户界面和管理界面。
用户界面包括用户点餐功能并将用户点餐信息存入账单文件中。
管理界面包括记录饺子的单价，库存等信息并将这些信息按一定规律排列供管理者参考，还要记录每天的收入与库存消耗。
2）系统组成框图系统组成如下图所示，点餐/帐目信息管理系统中任务调度模块是信息管理的指挥中心，所有的功能模块均通过该模块集中管理和调用。
数据文件是用于将改变的记录随时保存起来，I/O交互模块是指数据按键和控制按键的响应操作。
数据按键是窗口的输入输出。
系统平台I/O交互任务调度数据文件添加账目查询账目查询原料销量排序当天收入功能模块设计；本管理系统开发的过程中成功地完成很多函数的编写，而且全部通过程序调试。
下面针对与该系统相关的主要功能函数的编写思路和实现方法作总结。
1）数据描述与数据文件（1）数据描述；系统中共定义了三种结构体分别是struct list（账单信息）、struct dump（原材料信息）、struct system（管理信息）其中账单信息成员定义如下； struct list/*账单信息 */ { int num;/*编号*/ intmonth;/*月份*/ int date;/*日期*/ int table;/*桌号*/ int people;/*人数*/ int add;/*收款金额*/ };原材料信息；struct dump/*原材料信息*/ { int num;/*编号*/ char a[40];/*名称*/ int price;/*单价*/ int quantity;/*库存*/ };管理信息；struct system { struct dump data;/*点餐*/ int cash;/*收款*/ int sale;/*销量*/ int p;/*人数*/ };（2）数据文件；共定义四个数据文件；记录编号的文件count、记录原料的文件dumplings、记录账单的文件customer、数据处理文件system。
大连理工大学庞丽萍最优化方法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方法求函数极小点的计算程序。
大连理工大学2016电路分析作业
1-3 按题1-3图所示参考方向及给定值，判别各元件中电流和电压的实际方向。
计算各元件中的功率，说明元件是吸收还是发出功率。
题1-3图1-7 求题1-7所示局部电路中电压U。
题1-7图题1-8图1-8 题1-8图所示电路中，已知I1＝2A，I3＝－3A，U1＝10V，U4＝－5V，试指出各支路电压、电流的实际方向，并计算各元件吸收的功率。
1-10 求题1-10图中U s与I各为多少?题1-10图1-12 题1-12图所示各电路中，求K打开及闭合时的U a，U b，及U ab。
（a）(b) (c)题1-12图1-17 题1-17图中，（1）求图(a)中电压u AB；（2）图(b)中若u AB＝6V，求电流I。
B B题1-17图1-18 求题1-18图所示电路中电压U。
题1-18图题1-19图1-19 求题1-19图所示电路中电压源和电流源发出的功率。
1-22 题1-22图中，求U和I，若1A电流源换成10A电流源，重新求解。
题1-22图题1-25图1-25题1-25图（b），求电压U ab及U cb。
1-26题1-26图(b)中，已知U s =10V ，I 1=2A ，R 1=4.5Ω，R 2=1Ω，求电流I 2。
题1-26图题1-28图1-28 题1-28图中，若U s ＝－15.5V ，U 1＝1V ，求R ＝？1-29 题1-29图所示电路为MOSFET 放大器（Metal Oxide Semiconductor Field Effect Transistor ，将在《模拟电子技术》课中学习）交流特性等效电路，如果g m ＝4mS ，计算u out 。
u out3sin10t题1-29图2-2 求题2-2图所示各电路中a 、b 端等效电阻R ab 。
2-4 求题2-4图所示各电路的等效电阻R ab 。
题2-4图2-5求题2-5图电路的输入电阻in R 。
+ U 2 －题2-5图2-7将题2-7图所示各电路分别对ab 端等效为一个电压源或电流源。
大连理工大学 秋季优化方法大作业
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;
优化作业流程的三种方法
优化作业流程的三种方法一、简化步骤。
1.1 去除冗余操作。
在作业流程里啊，常常有些步骤就像“画蛇添足”一样，没什么实际用处，还拖慢速度。
比如说，有些文件填写，要求填一大堆重复信息，这不是浪费时间嘛。
咱就得把这些多余的部分找出来，大刀阔斧地砍掉。
就像整理房间，把那些没用的杂物都扔掉，房间立马就清爽了。
1.2 合并相似任务。
有些任务啊，就像双胞胎似的，很相似。
那咱们就别分开做了，把它们合并起来。
例如，数据收集和初步整理，这俩事儿关联性很强，完全可以一起做。
这样就像把散落在各处的小珠子串成一条项链，一下子就整齐有序了，效率也能大大提高。
二、引入新技术。
2.1 自动化工具的使用。
现在科技这么发达，有好多自动化的工具就像“得力助手”一样。
像一些办公软件，能自动处理文档格式、进行数据计算。
以前人工做这些得花老长时间了，现在点几下鼠标就搞定。
比如说财务报表制作，有了专门的软件，数据输入进去，各种图表、分析立马就出来了，又快又准，简直是“神来之笔”。
2.2 信息化管理系统。
建立信息化管理系统那可是个好办法。
这就好比给作业流程装上了一个智能大脑。
所有的任务分配、进度跟踪、资源调配都能在这个系统里一目了然。
员工们也能清楚地知道自己该干什么，什么时候干。
就像大家都在一条清晰的轨道上行驶的列车，不会乱套。
三、员工培训与激励。
3.1 针对性培训。
员工要是对作业流程不熟悉，就像盲人摸象一样，只能瞎干。
所以要进行针对性的培训。
根据不同岗位的需求，把流程详细地教给员工。
这就像给战士们配上合适的武器并且教会他们怎么用。
比如生产线上的工人，要让他们清楚每个环节的操作规范，这样才能保证产品质量，提高生产效率。
3.2 激励机制。
人都是需要激励的，有了激励就像给汽车加足了油。
设立一些奖励制度，对于那些在作业流程优化方面有好点子或者执行得特别好的员工，给予奖励。
可以是奖金、荣誉称号或者晋升机会。
这样员工们就会像打了鸡血一样，积极主动地去寻找优化作业流程的方法，整个团队也就充满了活力。
优化方法课程大作业
优化方法课程上机大作业学部：电子信息与电气工程学部专业：生物医学工程班级：电信硕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）代码运行结果如下：。
优化方法
这是一个原料配制问题，是在生产任务确定的条件下， 这是一个原料配制问题，是在生产任务确定的条件下，合 理的组织生产，使所消耗的资源数最少的数学规划问题。 理的组织生产，使所消耗的资源数最少的数学规划问题。 满足一组约束条件的同时， 寻求变量x 满足一组约束条件的同时 ， 寻求变量 1 和 x2 的值使目标函 数取得最小值。 数取得最小值。
max z=5x1+2x2 30x1+20x2≤160 5x1+ x2 ≤15
甲种药品每周产量不应超过4吨的限制 甲种药品每周产量不应超过 吨的限制 计划生产数不可能是负数， 计划生产数不可能是负数，
x1≤4 x1≥0 x2≥0
每吨产品的消耗 甲 维生素（公斤） 维生素（公斤） 设备（ 设备（台） 单位利润(万元 单位利润 万元) 万元 30 5 5 乙 20 1 ２
主要参考书目： 主要参考书目： 理论方面： 理论方面： (1) 解可新、韩健，《最优化方法》，天津大学出版社 解可新、韩健， 最优化方法》 天津大学出版社,2004 (2) 何坚勇， 《最优化方法》, 清华大学出版社， 2007 何坚勇， 最优化方法》 清华大学出版社， 计算方面： 计算方面： (3) 曹卫华，郭正， 《最优化技术方法及 曹卫华，郭正， 最优化技术方法及MATLAB的实现》， 的实现》 的实现 化学工业出版社， 化学工业出版社，2005 (4) 朱德通，《最优化模型与实验》, 同济大学出版社， 2003 朱德通， 最优化模型与实验》 同济大学出版社， 其它参考书： 其它参考书： (5)卢名高、刘庆吉编著，《最优化应用技术》，石油工业出版社 卢名高、 卢名高 刘庆吉编著， 最优化应用技术》 石油工业出版社,2002 (6)唐焕文，秦学志 《实用最优化方法》，大连理工大学出版社，2004 唐焕文， 唐焕文 秦学志,《实用最优化方法》 大连理工大学出版社， (7)钱颂迪，《运筹学》，清华大学出版社，1990 钱颂迪， 运筹学》 清华大学出版社， 钱颂迪 (8)袁亚湘、孙文瑜著，《最优化理论与方法》，科学出版社，2005 袁亚湘、 袁亚湘 孙文瑜著， 最优化理论与方法》 科学出版社，
大连理工大学优化方法上机大作业
学院：专业：班级：学号：姓名：上机大作业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。
最优化方法部分课后习题解答(1-7)
最优化方法部分课后习题解答习题一1.一直优化问题的数学模型为：22121122123142min ()(3)(4)5()02()50..()0()0f x x xg x x x g x x x s t g x x g x x =−+−⎧=−−≥⎪⎪⎪=−−+≥⎨⎪=≥⎪=≥⎪⎩试用图解法求出：（1）无约束最优点，并求出最优值。
（2）约束最优点，并求出其最优值。
（3）如果加一个等式约束，其约束最优解是什么？12()0h x x x =−=解：（1）在无约束条件下，的可行域在整个平面上，不难看出，当=（3，4）()f x 120x x *x 时，取最小值，即，最优点为=（3，4）：且最优值为：=0()f x *x *()f x （2）在约束条件下，的可行域为图中阴影部分所示，此时，求该问题的最优点就是()f x 在约束集合即可行域中找一点，使其落在半径最小的同心圆上，显然，从图示中可12(,)x x 以看出，当时，所在的圆的半径最小。
*155(,)44x =()f x 其中：点为和的交点，令求解得到：1()g x 2()g x 1122125()02()50g x x x g x x x ⎧=−−=⎪⎨⎪=−−+=⎩1215454x x ⎧=⎪⎪⎨⎪=⎪⎩即最优点为：最优值为：=*155(,)44x =*()f x 658（3）.若增加一个等式约束，则由图可知，可行域为空集，即此时最优解不存在。
2.一个矩形无盖油箱的外部总面积限定为S，怎样设计可使油箱的容量最大？试列出这个优化问题的数学模型，并回答这属于几维的优化问题.解：列出这个优化问题的数学模型为：该优化问题属于三维的优化问题。
123122313123max ()220..00f x x x x x x x x x x S x s t x x =++≤⎧⎪>⎪⎨>⎪⎪>⎩32123sx y z v⎛⎞=====⎜⎟⎝⎠习题二3.计算一般二次函数的梯度。
实用最优化方法大连理工课后答案
实用最优化方法大连理工课后答案
1．下列情况引起的误差是系统误差还是偶然误差？
(1)砝码锈蚀（系统误差）
(2)称重时试样吸收了空气中的水分；（系统误差）
(3)滴定管读数时末位数字估计不准嘶；（偶然误差）
(4)滴定剂中台有少量待测组分；睬统误差）
(5)标定用的基准物Na2C03在保存过程中吸收了水分：（系统误整）
(6)称量过程中天平零点由于环境条件的变化稍有变动（偶然误差）
2．什么是误差？什么是偏差？有什么区别和联系？
误差是测量值与真值之差偏差是单次测量值与n次测量平均值之差,误差是用测量位与真实值作比较，衡量准确度的高低，偏差是用测定值与平均位作比较，用于衡量青密度的大小，准确度高则精密度一定高，精密度高准确度不一定高。
大连理工大学庞丽萍最优化方法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方法求函数极小点的计算程序。
关于专科生毕业大作业的规定
附件5关于专科生毕业大作业的规定网络教育学院0 0 6. 2 .大连理工大学网络教育学院关于专科生毕业大作业的规定毕业大作业是专科生培养过程中最后一个教学环节，相当于本科生的毕业论文（设计），是学生学习后期的一次综合训练。
为做好毕业大作业工作，达到综合训练和提高教学质量的目的，特作如下规定。
一、目的要求1．毕业大作业的基本教学目的是培养学生综合运用所学的基础理论、专业知识和基本技能，提高分析与解决实际问题的能力以及科技写作或设计能力。
要求学生在做毕业大作业的过程中要有意识地培养自己的能力。
毕业大作业要按照各专业的教学计划要求进行，按期完成。
二、工作程序2．毕业大作业动员毕业大作业开始前，各学习中心要通过会议或网络进行毕业大作业动员，明确目的要求及有关规定。
同时，聘请好指导教师，并与学生见面。
3．选题原则上，学生可以在自己所从事的工作或所承担的项目中自选或自拟毕业大作业的题目，也可以依据某一门专业性课程的内容拟定一个大型作业作为毕业作业题目。
如自己不能立题，也可在指导教师提出的课题方向或题目中选题。
每位学生选定一个题目。
鼓励学生结合工作实际选题，但其选题必须与所学专业有关，并有一定的份量。
各学习中心要及时填报学生的毕业大作业题目、指导教师及安排情况统计表,并报奥鹏汇总，以便学院教学部备案。
4．毕业作业检查毕业作业进行过程中，各学习中心按要求进行阶段检查。
学院不定期组织抽查。
5．评定成绩学生的毕业大作业完成后，指导教师应予评阅并评出成绩，最后经学院教学部或指定所在的学习中心审查确认。
6．总结及资料保存毕业大作业成绩评定结束后，毕业大作业及成绩单等全部资料上交，由学院教学部安排专人统一保存。
三、写作要求7．毕业大作业的文体（题目的类型）1）．学生可选多种类型的题目，如做“工作研究型”、“技术开发型”、“调研报告型”、“读书报告型”、“大型作业型”、“工程设计型”和“课程设计型”（“设计型”毕业大作业必须带有设计说明书）的毕业大作业，即运用学习过的专业理论知识去解释、解决或回答社会现实生活或本人从事的技术或管理工作或课程学习中的有关问题等，侧重于理论应用，或技术开发与应用，或课程学习的深化。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
//f(x)//
function f=f(x) x1=[1 0]*x; x2=[0 1]*x; f=(1-x1)^2+100*(x2-x1^2)^2; //梯度函数// function g=g(x) x1=[1 0]*x; x2=[0 1]*x; g=[-2*(1-x1)-400*x1*(x2-x1^2); 200*(x2-x1^2)];
function h=h(x) x1=[1 0]*x; x2=[0 1]*x; h=[ 2+1200*x1^2-400*x2-400*x1; -400*x1 200]; //运行过程// >> x=[0;0] x= 0 0 >> di1tinewton(x) after 2 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal "x" is "ans". ans = 1.0000 1.0000 3°BFGS 方法
优化方法
期末上机大作业
姓 名：李岚松 学 部：电信学部电气工程 学 号：31609020
2016 年 11 月 9 日
1°最速下降法 //最速下降法主函数// function lls=di1titidu(x) =di1titidu(x) x0=x; eps=1e-4; 4; k=0; g0=g(x0); while (k>=0) if norm(g0)<eps break; else lanmed=10;c=0.1;i=0; while i>=0&&i<100 100 x=x0+lanmed*s0; if f(x)>(f(x0)+c*lanmed*g0'*s0) lanmed=lanmed/2; i=i+1; else break; end end x=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; end end lls=x; sll=f(x); fprintf('after %d iterations,obtain the optimal solution.\n\ solution. \nThe optimal solution is %f.\n\n n The optimal "x" is "ans".\n',k,sll); "ans". s0= s0=-g0;
//共轭梯度法主函数// function [x,val,k]=frcg(fun,gfun,x0) rho=0.6;sigma=0.4; k=0; epsilon=1e-4; 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; else beta=(g'*g)/(g0'*g0); d=-g+beta*d0; gd=g'*d; if(gd>=0.0) d=-g; end end if(norm(g)<epsilon), break; end m=0; mk=0; while(m<20) if(feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d) mk=m; break; end
//运行过程//
>> x=[0;0] x= 0 0 >> di1tiBFGS(x) after 14 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal "x" is "ans". ans = 1.0000 1.0000 4°共轭梯度法
//运行过程// >> x=[0;0] x= 0 0 >> di1titidu(x) after 4574 iterations,obtain the optimal solution. The optimal solution is 0.000000. The optimal "x" is "ans". ans = 1.0001 1.0002 2°牛顿法 //牛顿法主函数// function lls=di2tinewton(x) x0=x; eps=1e-4; k=0; g0=g(x0); h0=h(x0);s0=-inv(h0)*g0; while (k>=0&&k<1000) if norm(g0)<eps break; else x=x0+s0; x0=x; g0=g(x); h0=h(x); s0=-inv(h0)*g0; k=k+1; end end
lls=x; sll=f(x); fprintf('after %d iterations,obtain the optimal solution.\n\nThe optimal solution is %f.\n\n The optimal "x" is "ans".\n',k,sll); //f(x)//
//梯度函数//
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]' x0 = 0 0 >> [x,val,k]=frcg('fun','gfun',x0) x= 1.0001 1.0002
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=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);
//f(x)//
function f=f(x) x1=[1 0]*x;
x2=[0 1]*x; f=(1-x1)^2+100*(x2-x1^2)^2; //梯度函数// function g=g(x) x1=[1 0]*x; x2=[0 1]*x; g=[-2*(1-x1)-400*x1*(x2-x1^2); 200*(x2-x1^2)];
val = 7.2372e-09 k= 122
//增广拉格朗日法主函数// function [x,mu,lambda,output]=multphr(fun,hf,gf,dfun,dhf,dgf,x0) maxk=500; sigma=2.0; eta=2.0; theta=0.8; k=0; ink=0; epsilon=1e-4; 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); lambda=0.1*ones(m,1); btak=10; btaold=10; while(btak>epsilon (btak>epsilon & k<maxk) [x,ival,ik]=bfgs('mpsi' 'mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,s ,x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,s igma); ink=ink+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),lambda(i)/sigma); btak=btak+temp^2; end btak=sqrt(btak); if btak>epsilon if(k>=2 (k>=2 & btak > theta*btaold) sigma=eta*sigma; end for i=1:l, mu(i)=mu(i)-sigma*he(i); mu(i)=mu(i) end for i=1:m lambda(i)=max(0.0,lambda(i) )=max(0.0,lambda(i)-sigma*gi(i)); end end
//增广拉格朗日函数//
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