大连理工大学优化方法上机大作业
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。
大连理工优化方法大作业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。
《机械优化设计》大作业
高等流体力学班级：机设15学硕班学号： ********** *名：***授课老师：毕新胜日期： 2016年7月 1日一、研究报告内容：1、λ＝0.618的证明、一维搜索程序作业；2、单位矩阵程序作业；3、连杆机构问题＋自行选择小型机械设计问题或其他工程优化问题；（1）分析优化对象，根据设计问题的要求，选择设计变量，确立约束条件，建立目标函数，建立优化设计的数学模型并编制问题程序；（2）选择适当的优化方法，简述方法原理，进行优化计算；（3）进行结果分析，并加以说明。
4、写出课程实践心得体会，附列程序文本。
5、为响应学校2014年度教学工作会议的改革要求，探索新的课程考核评价方法，特探索性设立一开放式考核项目，占总成绩的5%。
试用您自己认为合适的方式（书面）表达您在本门课程学习方面的努力、进步与收获。
（考评将重点关注您的独创性、简洁性与可验证性）。
二、研究报告要求1、报告命名规则：学号－姓名－《机械优化设计》课程实践报告.doc2、报告提交邮址：*****************.cn（收到回复，可视为提交成功）。
追求：问题的工程性，格式的完美性，报告的完整性。
不追求：问题的复杂性，方法的惟一性。
评判准则：独一是好，先交为好；切勿拷贝。
目录：λ＝0.618的证明、一维搜索程序作业① 关于618.0=λ的证明……………………………………………………4 ② 一维搜索的作业采用matlab 进行编程…………………………………………… 5 采用C 语言进行编程……………………………………………… 7 单位矩阵程序作业① 采用matlab 的编程………………………………………………… 9 ② 采用c 语言进行编程………………………………………………… 9 机械优化工程实例① 连杆机构...........................................................................11 ② 自选机构...........................................................................16 课程实践心得.............................................................................. 20 附列程序文本.............................................................................. 21 进步，努力，建议 (25)一、λ＝0.618的证明、一维搜索程序作业①关于618.0=λ的证明黄金分割法要求插入点1α，2α的位置相对于区间],[b a 两端具有对称性，即)(1a b b --=λα)(2a b a -+=λα其中λ为待定常数。
大连理工矩阵上机作业
第一题Lagrange插值函数function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endx0=[1:10];y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777];lagrange(x0,y0,17)ans=-1.9516e+12x=[1:0.1:10];x=x';plot(x0,y0,'r');hold onplot(x,y,'k');legend('原函数','拟合函数')拟合图像如下拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断累积，导致龙格现象的发生。
第二题function fun =nihe(n)m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063 *10^6,74.669*10^6,74.487*10^6,74.065*10^6,76.777*10^6];w=[1,2,3,4,5,6,7,8,9,10];d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p))-m(p))^2;d2=d2+(subs(y2,w(p))-m(p))^2;d3=d3+(subs(y3,w(p))-m(p))^2;endd1=sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=[f1 f2 f3;d2 d3 d4];return;结果三次函数的均方误差最小，拟合的最好。
大连理工大学优化方法上机作业
大连理工大学优化方法上机作业本页仅作为文档页封面，使用时可以删除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。
最优化方法大作业模板
命题人：审核人：大作业学期：至学年度第学期课程：最优化方法课程代号：签到序号：使用班级：姓名：学号：题号一二三四五六七八九十总分得分一、（目标1）请从以下6种算法中任选一种，说明算法的来源、定义、基本思想和优缺点，并给出算法步骤（包含算法流程图）和例子（包含程序与运算结果）。
①禁忌搜索算法；②模拟退火算法；③遗传算法；④神经网络算法；⑤粒子群算法；⑥蚁群算法。
二、（目标1）某工厂生产甲、乙两种产品，已知生产这两种产品需要消耗三种材料A 、B 和C ，其中生产过程中材料的单位产品消耗量和总量，以及单位产品的利润如下表所示。
该如何配置安排生产计划，使得工厂所获得的利润最大？材料甲乙资源总量材料A （Kg ）3265材料B （Kg ）2140材料C （Kg ）0375单位利润（元/件）15002500-（1）要保证工厂利润的最大化，写出相应的生产计划数学模型；（2）根据对偶理论，直接写出该线性规划的对偶问题；（3）采用单纯形表法对该该线性规划问题进行求解，写出详细的计算过程；（4）采用Matlab 软件对该线性规划问题进行求解，写出完整的源程序，并给出程序运行结果；（5）讨论当材料B 的资源总量发生变化时，该线性规划问题的最优解会如何变化？课程目标目标1……题号一、二、三、四、五……分值20、25、20、20、15……得分得分三、（目标1）求解下列无约束非线性规划问题（1）采用黄金分割法求解：min 4()24f x x x =++。
初始区间为[-1.0]，精度为ε=10-4。
（要求：采用黄金分割法进行Matlab 编程求解，写出源程序，并给出运行结果，列出迭代过程的数据表格）（2）采用阻尼牛顿法求解：222121212min (,)4f x x x x x x =+-。
分别取两个初始点：x A =(1,1)T ，x B =(3,4)T 。
（要求：采用阻尼牛顿法进行Matlab 编程求解，并给出运行结果，列出迭代过程的数据表格）四、（目标1）求解下列约束非线性规划问题：22112212121212min ()23532..00f x x x x x x x x x x x s t x x =-+--+≤⎧⎪-≤⎪⎨≥⎪⎪≥⎩（1）采用罚函数法进行求解，需写出具体计算过程；（2）采用二次规划方法进行求解，需写出具体计算过程，并进行MATLAB 编程，写出源程序和运算结果；五、（目标1）（1）某商店在未来的4个月里，准备利用它的一个仓库来专门经营某种商品，仓库的最大容量为1000单位，而且该商店每月只能出卖仓库现有的货。
大连理工大学庞丽萍最优化方法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;
大连理工大学--优化作业----程序共9页文档
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。
优化方法课程大作业
优化方法课程上机大作业学部：电子信息与电气工程学部专业：生物医学工程班级：电信硕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．下列情况引起的误差是系统误差还是偶然误差？
(1)砝码锈蚀（系统误差）
(2)称重时试样吸收了空气中的水分；（系统误差）
(3)滴定管读数时末位数字估计不准嘶；（偶然误差）
(4)滴定剂中台有少量待测组分；睬统误差）
(5)标定用的基准物Na2C03在保存过程中吸收了水分：（系统误整）
(6)称量过程中天平零点由于环境条件的变化稍有变动（偶然误差）
2．什么是误差？什么是偏差？有什么区别和联系？
误差是测量值与真值之差偏差是单次测量值与n次测量平均值之差,误差是用测量位与真实值作比较，衡量准确度的高低，偏差是用测定值与平均位作比较，用于衡量青密度的大小，准确度高则精密度一定高，精密度高准确度不一定高。
实用最优化方法编程大作业
实用最优化方法编程大作业班级：姓名：指导老师：学号：大连理工大学2015年11月27日版本号：MATLAB7.11.0(R2010b)【文件名WP.m】function x=WP(f,x0,var,s0,eps)clcsyms x1x2;c1=0.1;c2=0.5;b=inf;lambda=1;if nargin==4eps=1.0e-6;endgradf=jacobian(f,var);g0=subs(gradf,var,x0);f0=subs(f,var,x0);gs=g0*s0';a=0;j=0;while j<1000new_x=x0+lambda*s0;new_f=subs(f,var,new_x);left=f0-new_f;new_g=subs(gradf,var,new_x);new_gs=new_g*s0';right=-1*c1*lambda*gs;if left<right%不满足第一个条件j=j+1;b=lambda;lambda=0.5*(lambda+a);else%满足第一个条件left2=new_gs;right2=c2*g0;if left2<right2%不满足第二个条件j=j+1;a=lambda;lambda=min(2*lambda,0.5*(lambda+b));elsex=lambda;break;endendend在Command Window输入：syms x1x2WP(100*(x2-x1^2)^2+(1-x1)^2,[-1,1],[x1x2],[11])程序运行后可得出结果：ans=0.0039【文件名minGETD.m】function[x,minf]=minGETD(f,x0,var,eps)clcsyms x1x2x3format long;n=length(x0);if nargin==3eps=1.0e-3;endx0=x0';syms lambda;gradf=jacobian(f,var);g=subs(gradf,var,x0);s=-g';k=0;while1tol=norm(double(g));if tol<=epsx=x0;break;endx1=x0+lambda*s;f1=subs(f,var,x1);dy1=diff(f1);lambda0=solve(dy1);x1=x0+lambda0*s;g1=subs(gradf,var,x1)tol=norm(double(g1))if tol<=epsx=x1;break;endif k+1==nx0=x1;continue;elsemiu=dot(g1,g1)/dot(g,g)s=-g1'+miu*s;k=k+1;x0=x1;endendx在Command Window输入：syms x1x2x3x0=[000];var=[x1x2x3];f=x1^2-2*x1*x2+2*x2^2+x3^2-x2*x3+2*x1+3*x2-x3;minGETD(f,x0,var,eps)程序运行后可得出结果：x=[-236894/59711,-178563/59711,-59465/59711]可认为最终解为[-4,-3,-1]。
大连理工大学概率上机作业
大连理工大学概率上机作业————————————————————————————————作者: ————————————————————————————————日期：ﻩ第一次上机作业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个随机数的均值和方差。
教学内容基本要求与学时分配-大连理工大学教务处
《机械设计基础A》教学大纲（学分4 学时64）一、课程说明本课程是工科近机械类（包括机械类某些专业）和非机械类专业大类课程之一，是工科学生学习和掌握各种类型的机械中常用机构和通用机械零件的基本知识和基本设计方法的技术基础课。
该课程也是工科学生将来学习专业机械设备课程的理论基础。
本课程在教学内容方面着重基本知识、基本理论和基本设计方法的讲解；在培养实践能力方面着重设计构思和基本设计技能的基本训练。
二、课程目标（需对应于本专业2013级培养方案中的毕业生能力进行细化分解）1．学习机械工程基础知识和基本理论知识，掌握常用机构的结构、特性等基本知识，了解各种机械的传动原理，具有分析、选用和设计机械设备中基本机构的能力；2．通用机械零件的设计原理、方法和机械设计等的一般规律，具有设计机械传动装置和简单机械的能力；3．掌握基本的机械设计创新方法，培养学生追求创新的态度和意识；4．培养学生树立正确的设计思想，了解机械设计过程中国家有关的经济、环境、法律、安全、健康、伦理等政策和制约因素；5．培养学生的工程实践学习能力，使学生掌握典型零件的实验方法，获得实验技能的基本训练，具有运用标准、规范、手册、图册和查阅有关技术资料的能力；6．了解机械设计的前沿和新发展动向。
三、教学内容、基本要求与学时分配序号教学内容教学要求学时教学方式对应能力1 一、基本概念1.研究的对象、内容；2.机械设计的基本要求和一般设计过程。
1.了解本课程研究的对象、内容2.了解机械设计的基本要求、一般设计过程。
2 讲授2、42 二、平面机构的自由度和速度分析1.机构运动简图2.平面机构自由度3.速度瞬心1.了解平面机构运动简图的绘制。
2.掌握平面机构自由度的计算以及机构具有确定运动的条件。
3.掌握速度瞬心的求解方法以及在平面机构运动分析上的应用。
3讲授、上机1、53 三、平面连杆机构及其设计1.平面四杆机构2.平面四杆机构的设计方法1.了解平面四杆机构的特点、应用、基本型式及其演化。
优化方法MATLAB编程——大连理工大学
The optimal solution is 0.000000.
The optimal "x" is "ans".
ans =
1.0000 1.0000 1.0000 1.0000 可以看出，用Newton法经过6次迭代就能求出最优解。
3. BFGS法 源程序如下： function zy_x=di2tiBFGS(x) %该函数用来解大作业第二题。 x0=x; yimuxulong=1e-5; k=0; g0=g(x0); H0=eye(4);s0=-H0*g0;
3
s0=-g0;
大连理工大学优化方法上机大作业
end end x=x0+lanmed*s0; x0=x; g0=g(x); s0=-g0; k=k+1; end end zy_x=x; zyj=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,zyj);
4
大连理工大学优化方法上机大作业
>> x=[-1.2 1 -1.2 1]'; >> di2titidu(x) after 5945 iterations,obtain the optimal solution.
The optimal solution is 0.000000.
The optimal "x" is "ans".
大连理工大学优化方法上机大作业
大连理工大学优化方法上机大作业
学院：专业：班级：学号：姓名：上机大作业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、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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 + *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(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);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 is x1 =法：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)')/((f a0)'*(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 + *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+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); %%乘法系数最大1000 disp('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。