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

2016年大连理工大学优化方法上机大作业学院：专业：班级：学号：姓名：上机大作业1：1.最速下降法：functionf=fun(x)f=(1-x(1))^2+100*(x(2)-x(1)^2)^2;endfunctiong=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);endfunctionx_star=steepest(x0,eps)gk=grad(x0);res=norm(gk);k=0;whilereseps&&k=1000dk=-gk;ak=1;f0=fun(x0);f1=fun(x0+ak*dk);slope=dot(gk,dk);whilef1f0+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-thiter,theresidualis%f\n',k,res);endx_star=xk;endclearx0=[0,0]';eps=1e-4;x=steepest(x0,eps)2.牛顿法：functionf=fun(x)f=(1-x(1))^2+100*(x(2)-x(1)^2)^2;endfunctiong=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;endfunctiong=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);endfunctionx_star=newton(x0,eps)gk=grad(x0);bk=[grad2(x0)]^(-1);res=norm(gk);k=0;whilereseps&&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-thiter,theresidualis%f\n',k,res);endx_star=xk;endclearx0=[0,0]';eps=1e-4;x1=newton(x0,eps)--The1-thiter,theresidualis447.213595--The2-thiter,theresidualis0.000000x1=1.00001.00003.BFGS法：functionf=fun(x)f=(1-x(1))^2+100*(x(2)-x(1)^2)^2;endfunctiong=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);endfunctionx_star=bfgs(x0,eps)g0=grad(x0);gk=g0;res=norm(gk);Hk=eye(2);k=0;whilereseps&&k=1000dk=-Hk*gk;ak=1;f0=fun(x0);f1=fun(x0+ak*dk);slope=dot(gk,dk);whilef1f0+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-thiter,theresidualis%f\n',k,res);endx_star=xk;Endclearx0=[0,0]';eps=1e-4;x=bfgs(x0,eps)4.共轭梯度法：functionf=fun(x)f=(1-x(1))^2+100*(x(2)-x(1)^2)^2;endfunctiong=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);endfunctionx_star=CG(x0,eps)gk=grad(x0);res=norm(gk);k=0;dk=-gk;whilereseps&&k=1000ak=1;f0=fun(x0);f1=fun(x0+ak*dk);slope=dot(gk,dk);whilef1f0+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-thiter,theresidualis%f\n',k,res);endx_star=xk;endclearx0=[0,0]';eps=1e-4;x=CG(x0,eps)上机大作业2：functionf=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);endfunctionf=alobj(x)%拉格朗日增广函数%N_equ等式约束个数?%N_inequ不等式约束个数N_equ=1;N_inequ=3;globalr_alpena;%全局变量h_equ=0;h_inequ=0;[h,g]=constrains(x);%等式约束部分?fori=1:N_equh_equ=h_equ+h(i)*r_al(i)+(pena/2)*h(i).^2;end%不等式约束部分fori=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;functionf=compare(x)globalr_alpenaN_equN_inequ;N_equ=1;N_inequ=3;h_inequ=zeros(3,1);[h,g]=constrains(x);%等式部分fori=1:1h_equ=abs(h(i));end%不等式部分fori=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:最优函数值%============================程序开始================================globalr_alpena;%参数（全局变量）pena=10;%惩罚系数r_al=[1,1,1,1];c_scale=2;%乘法系数乘数cta=0.5;%下降标准系数e_al=1e-4;%误差控制范围max_itera=25;out_itera=1;%迭代次数%===========================算法迭代开始=============================whileout_iteramax_iterax_al0=x_al;r_al0=r_al;%判断函数?compareFlag=compare(x_al0);%无约束的拟牛顿法BFGS[X,fmin]=fminunc(@alobj,x_al0);x_al=X;%得到新迭代点%判断停止条件?ifcompare(x_al)e_aldisp('wegettheoptpoint');breakend%c判断函数下降度?ifcompare(x_al)cta*compareFlagpena=1*pena;%可以根据需要修改惩罚系数变量elsepena=min(1000,c_scale*pena);%%乘法系数最大1000disp('pena=2*pena');end%%?更新拉格朗日乘子[h,g]=constrains(x_al);fori=1:1%%等式约束部分r_al(i)=r_al0(i)+pena*h(i);endfori=1:3%%不等式约束部分r_al(i+1)=max(0,(r_al0(i+1)+pena*g(i)));endout_itera=out_itera+1;end%+++++++++++++++++++++++++++迭代结束+++++++++++++++++++++++++++++++++disp('theiterationnumber');k=out_itera;disp('thevalueofconstrains');compare(x_al)disp('theoptpoint');x=x_al;fmin=obj(X);clearx_al=[0,0];[x,fmin,k]=almain(x_al)上机大作业3：1、clearalln=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_beginvariablex(n)minimize(c'*x)subjecttoA*x=bcvx_endCallingSDPT34.0:6variables,3equalityconstraints------------------------------------------------------------num.ofconstraints=3dim.oflinearvar=6*******************************************************************SDPT3:Infeasiblepath-followingalgorithms*******************************************************************versionpredcorrgamexponscale_dataNT10.00010itpstepdsteppinfeasdinfeasgapprim-objdual-objcputime-------------------------------------------------------------------0|0.000|0.000|1.1e+01|5.1e+00|6.0e+02|-7.000000e+010.000000e+00|0:0:00|chol111|0.912|1.000|9.4e-01|4.6e-02|6.5e+01|-5.606627e+00-2.967567e+01|0:0:01|chol11