最优化方法大作业---------用优化算法求解函数最值问题摘要最优化(optimization) 是应用数学的重要研究领域.它是研究在给定约束之下如何寻求某些因素(的量),以使某一(或某些)指标达到最优的一些学科的总称。
最优化问题一般包括最小化问题和最大化问题,而最大化问题可以通过简单的转化使之成最最小化问题。
最小化问题分为两类,即约束最小化和无约束最小化问题。
在此报告中,前两个问题属于无约束最小化问题的求解,报告中分别使用了“牛顿法”和“共轭梯度法”。
后两个问题属于有约束最小化问题的求解,报告中分别用“外点法”和“内点法”求解。
虽然命名不一样,其实质都是构造“惩罚函数”或者“障碍函数”,通过拉格朗日乘子法将有约束问题转化为无约束问题进行求解。
再此报告中,“外点法”和“内点法”分别用了直接求导和调用“牛顿法”来求解无约束优化问题。
在此实验中,用“共轭梯度法”对“牛顿法”所解函数进行求解时出现错误,报告中另取一函数用“共轭梯度法”求解得到正确的结果。
此实验中所有的函数其理论值都是显见的,分析计算结果可知程序正确,所求结果误差处于可接受范围内。
报告中对所用到的四种方法在其使用以前都有理论说明,对“外点法”中惩罚函数和“内点法”中障碍函数的选择也有相应的说明,另外,对此次试验中的收获也在报告的三部分给出。
本报告中所用程序代码一律用MATLAB编写。
【关键字】函数最优化牛顿法共轭梯度法内点法外点法 MATLAB一,问题描述1,分别用共轭梯度发法和牛顿法来求解一下优化问题()()()()()441432243221102510min x x x x x x x x x f -+-+-++=2, 分别用外点法和内点发求解一下优化问题⎩⎨⎧≥-++01.min 212231x x t s x x二、问题求解用牛顿法求解()()()()()441432243221102510min x x x x x x x x x f -+-+-++=1.1.1问题分析:取步长为1而沿着牛顿方向迭代的方法称为牛顿法,在牛顿法中,初始点的取值随意,在以后的每次迭代中,()[]()k k k k x f x f x x ∇∇-=-+121,直到终止条件成立时停止。
1.1.2 问题求解注:本程序编程语言为MATLAB ,终止条件为()162110-≤∇x f ,初始取值为[10 10 10 10]M 文件(求解函数)如下:function s=newton1(f,c,eps) %c 是初值,eps 为允许误差值 if nargin==2 eps=;elseif nargin<1 error('') % return endsyms x1 x2 x3 x4x=[x1 x2 x3 x4].';grad = jacobian(f).';hesse = jacobian(grad);a=grad;b=hesse;i=1;gradk=subs(a,[x1 x2 x3 x4],[c(1) c(2) c(3) c(4)]); hessek=subs(b,[x1 x2 x3 x4],[c(1) c(2) c(3) c(4)]);pk=-1*(hessek\gradk);x=tihuan(c);while norm(gradk)>=epsx=x+pk;gradk=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]); hessek=subs(b,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]); pk=-hessek\gradk;i=i+1;enddisp('the times of iteration is:')disp(i)disp('The grad is:')disp(gradk)disp('and the result is:')x=x.';disp(x)return“tihuan”子函数:function x=tihuan(x)x(1)=x(1);x(2)=x(2);x(3)=x(3);x(4)=x(4);end调用方式如下:syms x1 x2 x3 x4f=(x1+10*x2)^2+5*(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4;c=[10 10 10 10]';%初始值newton1(f,c,eps);1.1.3 计算结果如下:由上述结果可知,当迭代次数达到47次时满足终止条件,此时x 为* [ ],显然,此题的理论解为[0 0 0 0],分析上述结果,与理论解的误差处于可接受范围之内。
求解完成。
用共轭梯度法求解函数()()()()()441432243221102510min x x x x x x x x x f -+-+-++=用共轭梯度法求解上述函数的程序代码如下:1.2.1问题分析: 取()00x f p -∇=,当搜索到1+k x 时,共轭方向()2,...,1,0,11-=+-∇=++n k p x f p k k k k λ,此时,1+k p 与k p A 共轭,用k Ap 右乘上式得()kTk k k k k k Ap p Ap x f Ap p λ+-∇=++11,由1=+k Tk Ap p 得()2,...,1,01-=∇=+n k Ap p Ap x f kT p kTk k λ,若不满足条件,进行下一次迭代。
1.1.2 问题求解注:程序所用语言为MATLAB ,精度为1610-=epssyms x1 x2 x3 x4 t0 t1f=(x1+10*x2)^2+5*(x3-x4)^2+(x2-2*x3)^4+10*(x1-x4)^4;c=[10;10;10;10];grad1 = diff(f,x1);grad2=diff(f,x2);grad3 = diff(f,x3);grad4=diff(f,x4);grad=[grad1;grad2;grad3;grad4];a=grad;i=1;n=40;gradk=subs(a,[x1 x2 x3 x4],[c(1) c(2) c(3) c(4)]);x=tihuan(c);p0=0;while norm(gradk)>=epsp0=-gradk;y=x;x=x+t0*p0;k=0;gradk=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]);w=solve(gradk(1)+gradk(2)+gradk(3)+gradk(4));t0=real(w);t0=eval(t0);t0=t0(1);x=y+t0*p0;gradk=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]);while norm(gradk)>=epsif k+1~=ngradk2=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]); gradk1=subs(a,[x1 x2 x3 x4],[y(1) y(2) y(3) y(4)]); lamda=norm(gradk2).^2/norm(gradk1).^2;p0=-gradk2+lamda*p0;k=k+1;elsek=0;p0=-subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]); endclear y; y=x;x=x+t1*p0;gradk=subs(f,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]);m=solve(gradk);t1=real(m); t1=eval(t1(1));x=x+t1*p0;x=eval(x);clear m;clear t1;syms t1gradk=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]);enddisp(x.') return; enddisp(x.')此程序为一初步程序,假设初值为[10;10;10;10],则第一次运算得t0=,lamda=,迭代后的x=NaN 。
现用共轭梯度法求解另一函数 ()222125min x x x f +=对上述程序稍加改动来求解本题的代码如下: 注:程序所用语言为MATLAB ,精度为1610-=epsfunction s=gongegrad2(f,c,eps) %c 是初值,eps 为允许误差值 if nargin==2 %eps=;elseif nargin<1 error('') return end ticsyms x1 x2 t0 t1 grad1 = diff(f,x1); grad2=diff(f,x2); grad=[grad1;grad2]; a=grad; i=1;n=40;gradk=subs(a,[x1 x2],[c(1) c(2)]); x=tihuan(c); p0=0;while norm(gradk)>=eps p0=-gradk; y=x;x=x+t0*p0; k=0;gradk=subs(f,[x1 x2],[x(1) x(2)]); w=solve(gradk); t0=real(w); t0=eval(t0); t0=t0(1); x=y+t0*p0;gradk=subs(a,[x1 x2],[x(1) x(2)]);while norm(gradk)>=epsif k+1~=ngradk2=subs(a,[x1 x2],[x(1) x(2)]);gradk1=subs(a,[x1 x2],[y(1) y(2)]);lamda=norm(gradk2)^2/norm(gradk1)^2;p0=-gradk2+lamda*p0;k=k+1;elsek=0;p0=-subs(a,[x1 x2],[x(1) x(2)]);endclear y; y=x;x=x+t1*p0;gradk=subs(f,[x1 x2],[x(1) x(2)]);m=solve(gradk);t1=real(m); t1=eval(t1(1));x=y+t1*p0;clear m;clear t1;syms t1gradk=subs(a,[x1 x2],[x(1) x(2)]);enddisp(sprintf('the last point we want is [%f %f]',x(1),x(2))); disp(sprintf('the times used to recursion is %f',k));disp(sprintf('the function value is %f',x(1)^2+25*x(2)^2));tocreturn;enddisp(sprintf('the last point we want is [%f %f]',x(1),x(2))); disp(sprintf('the times used to recursion is %f',k));disp(sprintf('the function value is %f',x(1)^2+25*x(2)^2)); toc“tihuan”子函数为:function x=tihuan(x)[v,g]=size(x);for i=1:vx(i)=x(i);end程序调用方式为:clear allclcsyms x1 x2 t0 t1f=x1^2+25*x2^2;c=[2;2];%初值gongegrad2(f,c,eps)程序结果如下:由上述结果知,用共轭梯度法对上述函数求解需要迭代三次得到最优解0,此时x 为[0 0];符合理论分析的结果,求解完成。