最优化大作业

最优化方法大作业

---------用优化算法求解函数最值问题

摘要

最优化(optimization) 是应用数学的重要研究领域.它是研究在给定约束之下如何寻求某些因素(的量)，以使某一(或某些)指标达到最优的一些学科的总称。最优化问题一般包括最小化问题和最大化问题，而最大化问题可以通过简单的转化使之成最最小化问题。最小化问题分为两类，即约束最小化和无约束最小化问题。在此报告中，前两个问题属于无约束最小化问题的求解，报告中分别使用了“牛顿法”和“共轭梯度法”。后两个问题属于有约束最小化问题的求解，报告中分别用“外点法”和“内点法”求解。虽然命名不一样，其实质都是构造“惩罚函数”或者“障碍函数”，通过拉格朗日乘子法将有约束问题转化为无约束问题进行求解。再此报告中，“外点法”和“内点法”分别用了直接求导和调用“牛顿法”来求解无约束优化问题。

在此实验中，用“共轭梯度法”对“牛顿法”所解函数进行求解时出现错误，报告中另取一函数用“共轭梯度法”求解得到正确的结果。此实验中所有的函数其理论值都是显见的，分析计算结果可知程序正确，所求结果误差处于可接受范围内。

报告中对所用到的四种方法在其使用以前都有理论说明，对“外点法”中惩罚函数和“内点法”中障碍函数的选择也有相应的说明，另外，对此次试验中的收获也在报告的三部分给出。

本报告中所用程序代码一律用MATLAB编写。

【关键字】函数最优化牛顿法共轭梯度法内点法外点法 MATLAB

一，问题描述

1，

分别用共轭梯度发法和牛顿法来求解一下优化问题

()()()()()4

41432243221102510min x x x x x x x x x f -+-+-++=

2， 分别用外点法和内点发求解一下优化问题

??

?≥-++0

1.min 212

231x x t s x x

二、问题求解

1.1 用牛顿法求解

()()()()()4

414

322

432

21102510min x x x x x x x x x f -+-+-++=

1.1.1问题分析：

取步长为1而沿着牛顿方向迭代的方法称为牛顿法，在牛顿法中，初始点的取值随意，在以后的每次迭代中，()[]

()k k k k x f x f x x ??-=-+1

21，直到终止条件成立时停止。

1.1.2 问题求解

注：本程序编程语言为MATLAB ，终止条件为()162

110-≤?x f ，初始取值

为[10 10 10 10]

M 文件（求解函数）如下：

function s=newton1(f,c,eps)

%c 是初值，eps 为允许误差值 if nargin==2 eps=1.0e-16; elseif nargin<1 error('') % return end

syms x1 x2 x3 x4 x=[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)>=eps

x=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;

end

disp('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 x4

f=(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 为

1.0e-005 * [-0.1111 0.0111 0.0095 0.0095], 显然，此题的理论解为[0 0 0 0],分析上述结果，与理论解的误差处于可接受范围之内。求解完成。

1.2 用共轭梯度法求解函数

()()()()()4

414

322

432

21102510min 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 右乘上式

得

()k

T

k k k k k k Ap p Ap x f Ap p λ+-?=++11，由

1=+k T

k Ap p 得

()2,...,1,01-=?=+n k Ap p Ap x f k

T p k

T

k k λ，若不满足条件，进行下一次迭代。

1.1.2 问题求解

注：程序所用语言为MATLAB ，精度为1610-=eps

syms x1 x2 x3 x4 t0 t1

f=(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)>=eps

p0=-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)>=eps

if k+1~=n

gradk2=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;

else

k=0;

p0=-subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]); end

clear 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 t1

gradk=subs(a,[x1 x2 x3 x4],[x(1) x(2) x(3) x(4)]);

end

disp(x.') return; end

disp(x.')

此程序为一初步程序，假设初值为[10;10;10;10]，则第一次运算得t0=0.0011,lamda=0.9291,迭代后的x=NaN 。现用共轭梯度法求解另一函数 ()2

22

125min x x x f +=

对上述程序稍加改动来求解本题的代码如下： 注：程序所用语言为MATLAB ，精度为1610-=eps

function s=gongegrad2(f,c,eps) %c 是初值，eps 为允许误差值 if nargin==2

%eps=1.0e-16; elseif nargin<1 error('') return end tic

syms 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)>=eps

if k+1~=n

gradk2=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;

else

k=0;

p0=-subs(a,[x1 x2],[x(1) x(2)]);

end

clear 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 t1

gradk=subs(a,[x1 x2],[x(1) x(2)]);

end

disp(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

return;

end

disp(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:v

x(i)=x(i);

end

程序调用方式为：

clear all

clc

syms x1 x2 t0 t1

f=x1^2+25*x2^2;

c=[2;2];%初值

gongegrad2(f,c,eps)

程序结果如下：

由上述结果知，用共轭梯度法对上述函数求解需要迭代三次得到最优解0，此时x 为[0 0];符合理论分析的结果，求解完成。

2.1 用外点法法求解函数

?

??≥-++01.min 212

231x x t s x x

2.1.1 问题分析

外点法为序列无约束最优化方法，其基本思想为将条件函数吸收到目标函数中进行求解。其在数学上的直观理解是拉格朗日乘子法：

[]()()()[]?

??

???+=∑=m

i i x g M x f M X T 12,0min min ;min ，

[]M X T ;min 为总代价，()x f 为价格，()()[]∑=m

i i x g M 1

2

,0min 为罚款。即在经济学上总代价为价格和罚款的和。

此时

()()[]()

()()()

m i x g x g x g x g i i i i ~1000,0min 2

2

=?≥???=，当，

当

称[]()()∑=++=m

i i x g M x f M X T 1

;为增广目标函数，通常取

()()

()()???≤≥=+0

00

2

x g x g x g x g i i i i

当当

2.1.2 问题求解

两种方法求解程序如下：

，程序中 2.1.2.1注：程序所用语言为MATLAB，终止条件为()16

-x

g

?

10-

i

无约束优化部分通过求导实现。

M文件如下：

tic

clc

%c 是初值，eps为允许误差值

if nargin==1

eps=1.0e-16;

elseif nargin<1

error('')

return

end

syms m

[x1,x2]=solve('3*x1^2+2*m*(x1+x2-1)=0','2*x2+2*m*(x1+x2-1)=0');

t=1;k=1;

x1=limit(x1,m,t);

x2=limit(x2,m,t);

bound=max(eval(x1(1)+x2(1))-1,eval(x1(2)+x2(2))-1);

while -bound>eps

t=10*t;k=k+1;

[x1,x2]=solve('3*x1^2+2*m*(x1+x2-1)=0','2*x2+2*m*(x1+x2-1)=0');

x1=limit(x1,m,t);

x2=limit(x2,m,t);

bound=max(eval(x1(1)+x2(1))-1,eval(x1(2)+x2(2))-1);

end

x1=eval(x1);x2=eval(x2);

f1=x1(1)^3+x2(1)^2;

f2=x1(2)^3+x2(2)^2;

if f1

disp(sprintf('the final x is [%f %f]',x1(1),x2(1)));

disp(sprintf('the final function value is %f',f1));

else

disp(sprintf('the final x is [%f %f]',x1(2),x2(2)));

disp(sprintf('the final function value is %f',f2));

end

disp(sprintf('the times used to recursion is %f',k))

toc

调用方式如下

syms x1 x2 m

T=x1^3+x2^2+m*(x1+x2-1)^2;

waidian(T,eps);

function s=waidian(T,eps)

实验结果如下

由上述结果可知当迭代17次时能够达到终止条件，此时,

???==451416.0548584

.02

1x x

函数得到最优解0.368870.求解完成。

2.1.2.2，程序中无约束优化部分通过调用“牛顿法”完成

代码如下；

tic

clear all;close all;clc syms x1 x2

m=1;x0=[2;2];eps=1.0e-6; T=x1^3+x2^2+m*(x1+x2-1)^2; c=10;k=1; while 1

s=newton1(T,x0,eps);%调用newton 法 x1=s(1);x2=s(2);

if m*(x1+x2-1)^2>eps

m=c*m;k=k+1;syms x1 x2

T=x1^3+x2^2+m*(x1+x2-1)^2; else

disp(sprintf('the final result is [%f %f]',x1,x2)); disp(sprintf('the function value is %f',x1^3+x2^2)); disp(sprintf('the times used to recursion is %f',k)); break; end end toc

实验结果如下：

由上述结果可知当迭代7次时能够达到终止条件，此时,

???==451416.0548584

.02

1x x

函数得到最优解0.368869.求解完成。 2.2 用内点法法求解函数

???≥-++0

1.m i n 212

231x x t s x x

2.2.1 问题分析

同外点法，内点法也为序列无约束最优化方法，其基本思想为将条件函数吸收到目标函数中进行求解。 ()()()()()()()x f x B x x B x x wB x f w x I →∞++=不为边界时，为为边界时;

;

称()x B 为障碍函数、围墙函数或者惩罚函数。通常取 ()()()()

()

∑

∑

∑∞

=∞

=∞

==1

2

1

1

i 11

,ln -i i i i i x g x g x g x B 或者

若在上一次迭代中x 处于可解区域外部，则“围墙”不够高，在下一次迭代时加高“围墙函数”。 2.2.2 问题求解 求解程序如下：

注：程序所用语言为MATLAB ，本程序出于对运行时间和算法简单的考

虑，选用()()

∑

∞

==1

1i i x g x B

终止条件为()610?x g i ，此时可以认为当x 接近求解边界时()x B （围墙）为无穷大，以确保函数的自变量迭代发生在可解区域内部 2.2.2.1，程序中无约束优化部分通过求导实现。 tic

clear all; close all;clc; eps=10^6; syms w x1 x2

I=x1^3+x2^2+w*(1/(x1+x2-1)); x10=2;x20=2;

[x1,x2]=solve('3*x1^2-w/(x1+x2-1)^2=0','2*x2-w/(x1+x2-1)^2=0'); c=10;k=1;

x1=limit(x1,w,c); x2=limit(x2,w,c);

x1=eval(x1);x2=eval(x2);

x1=real(x1(1));x2=real(x2(1)); bound=(1/(x1+x2-1)); while bound<=eps c=c/10;k=k+1;

[x1,x2]=solve('3*x1^2-w/(x1+x2-1)^2=0','2*x2-w/(x1+x2-1)^2=0'); x1=limit(x1,w,c); x2=limit(x2,w,c);

x1=eval(x1);x2=eval(x2);

x1=real(x1(1));x2=real(x2(1)); bound=(1/(x1+x2-1)); end

disp(sprintf('the final x is [%f %f]',x1,x2));

disp(sprintf('the final function value is %f',x1(1)^3+x2(1)^2)); disp(sprintf('the times used to recursion is %f',k)) toc

实验结果如下：

由上述结果可知当迭代15次时能够达到终止条件，此时,

???==451416

.0548584

.021x x

函数得到最优解0.368870.求解完成，结果与“外点法”相同。 2.2.2.2，程序中无约束优化部分通过调用“牛顿法”完成（eps x x c k k ≤-=-收敛准则为,）

代码如下；

tic

clear all;close all;clc syms x1 x2

w=10;x0=[2;2];eps=1.0e-6;

I=x1^3+x2^2+w*(1/(x1+x2-1)^2); c=0.3;k=1; while 1

s=newton1(I,x0,eps);%调用newton 法 x1=s(1);x2=s(2);

if sqrt((x1-x0(1))^2+(x2-x0(2))^2)>eps w=c*w;k=k+1;x0(1)=x1;x0(2)=x2; syms x1 x2

I=x1^3+x2^2+w*(1/(x1+x2-1)^2); else

disp(sprintf('the final result is [%f %f]',x1,x2)); disp(sprintf('the function value is %f',x1^3+x2^2)); disp(sprintf('the times used to recursion is %f',k)); break; end end toc

实验结果如下：

由上述结果可知当迭代36次时能够达到终止条件，此时,

???==451418.0548585

.02

1x x

函数得到最优解0.368872.求解完成，结果与“外点法”十分近似。

三、此次实验的收获

此次试验是用优化算法求解一些函数最值。通过此次实验，我的收获主要有以下几点；

1，

函数求最值在以前的学习中基本都是从函数本身出发，研究函数的本身性质，进而得出理论最值，对于一些性质不明确的函数涉猎不多。在这门课程中，我们更多的是从一个初始点出发，通过一些函数解集的迭代过程来求解函数的最值。虽然有的算法也牵扯到函数的性质（函数是否可导，是否有Hesse 矩阵等），但相对于以前的求解方式，对函数本身的要求已经薄弱很多，这就意味着可以对更多的函数运用算法求解最值，扩大了可求解函数的范围。

2，

这门课程的学习多是理论为主，老师在课堂上讲了很多实效有用的算法，对算法的优缺点也多有说明，但没在实战中运用的话理解不会很深。这次实验给了我一个运用所学理论知识的完美机会，在此次实验中能够通过自己的努力将课堂上讲的理论算法

逐步运用到实战中，无疑对理论的理解是大有裨益的。

3，这次实验中另外一个比较重要的收获是学习了MATLAB的运用。以前没有学习过MATLAB，但在这次试验中，能够自己学

习MATLAB并且很拙劣的完成实验内容，这种收获对以后的科目

学习也是很重要的。