最优化大作业.docx
《最优化大作业.docx》由会员分享，可在线阅读，更多相关《最优化大作业.docx（20页珍藏版）》请在冰点文库上搜索。
最优化大作业
最优化方法大作业
---------用优化算法求解函数最值问题
摘要
最优化（optimization）是应用数学的重要研究领域.它是研究在给定约束之下如何寻求某些因素（的量），以使某一（或某些）指标达到最优的一些学科的总称。
最优化问题一般包括最小化问题和最大化问题，而最大化问题可以通过简单的转化使之成最最小化问题。
最小化问题分为两类，即约束最小化和无约束最小化问题。
在此报告中，前两个问题属于无约束最小化问题的求解，报告中分别使用了“牛顿法”和“共轭梯度法”。
后两个问题属于有约束最小化问题的求解，报告中分别用“外点法”和“内点法”求解。
虽然命名不一样，其实质都是构造“惩罚函数”或者“障碍函数”，通过拉格朗日乘子法将有约束问题转化为无约束问题进行求解。
再此报告中，“外点法”和“内点法”分别用了直接求导和调用“牛顿法”来求解无约束优化问题。
在此实验中，用“共轭梯度法”对“牛顿法”所解函数进行求解时出现错误，报告中另取一函数用“共轭梯度法”求解得到正确的结果。
此实验中所有的函数其理论值都是显见的，分析计算结果可知程序正确，所求结果误差处于可接受范围内。
报告中对所用到的四种方法在其使用以前都有理论说明，对“外点法”中惩罚函数和“内点法”中障碍函数的选择也有相应的说明，另外，对此次试验中的收获也在报告的三部分给出。
本报告中所用程序代码一律用MATLAB编写。
【关键字】函数最优化牛顿法共轭梯度法内点法外点法MATLAB
一，问题描述
1，分别用共轭梯度发法和牛顿法来求解一下优化问题
2，分别用外点法和内点发求解一下优化问题
二、问题求解
用牛顿法求解
1.1.1问题分析：
取步长为1而沿着牛顿方向迭代的方法称为牛顿法，在牛顿法中，初始点的取值随意，在以后的每次迭代中，
，直到终止条件成立时停止。
1.1.2问题求解
注：
本程序编程语言为MATLAB，终止条件为
，初始取值为[10101010]
M文件（求解函数）如下：
functions=newton1（f,c,eps）
%c是初值，eps为允许误差值
ifnargin==2
eps=;
elseifnargin<1
error（''）
%return
end
symsx1x2x3x4
x=[x1x2x3x4].';
grad=jacobian（f）.';
hesse=jacobian（grad）;
a=grad;
b=hesse;
i=1;
gradk=subs（a,[x1x2x3x4],[c
（1）c
（2）c（3）c（4）]）;
hessek=subs（b,[x1x2x3x4],[c
（1）c
（2）c（3）c（4）]）;
pk=-1*（hessek\gradk）;
x=tihuan（c）;
whilenorm（gradk）>=eps
x=x+pk;
gradk=subs（a,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
hessek=subs（b,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
pk=-hessek\gradk;
i=i+1;
end
disp（'thetimesofiterationis:
'）
disp（i）
disp（'Thegradis:
'）
disp（gradk）
disp（'andtheresultis:
'）
x=x.';
disp（x）
return
“tihuan”子函数:
functionx=tihuan（x）
x
（1）=x
（1）;
x
（2）=x
（2）;
x（3）=x（3）;
x（4）=x（4）;
end
调用方式如下：
symsx1x2x3x4
f=（x1+10*x2）^2+5*（x3-x4）^2+（x2-2*x3）^4+10*（x1-x4）^4;
c=[10101010]';%初始值
newton1（f,c,eps）;
1.1.3计算结果如下：
由上述结果可知，当迭代次数达到47次时满足终止条件，此时x为
*[],
显然，此题的理论解为[0000],分析上述结果，与理论解的误差处于可接受范围之内。
求解完成。
用共轭梯度法求解函数
用共轭梯度法求解上述函数的程序代码如下：
1.2.1问题分析：
取
，当搜索到
时，共轭方向
，此时，
与
A共轭，用
右乘上式得
，由
得
，若不满足条件，进行下一次迭代。
1.1.2问题求解
注：
程序所用语言为MATLAB，精度为
symsx1x2x3x4t0t1
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,[x1x2x3x4],[c
（1）c
（2）c（3）c（4）]）;
x=tihuan（c）;
p0=0;
whilenorm（gradk）>=eps
p0=-gradk;
y=x;
x=x+t0*p0;
k=0;
gradk=subs（a,[x1x2x3x4],[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,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
whilenorm（gradk）>=eps
ifk+1~=n
gradk2=subs（a,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
gradk1=subs（a,[x1x2x3x4],[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,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
end
cleary;y=x;
x=x+t1*p0;
gradk=subs（f,[x1x2x3x4],[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）;clearm;cleart1;symst1
gradk=subs（a,[x1x2x3x4],[x
（1）x
（2）x（3）x（4）]）;
end
disp（x.'）
return;
end
disp（x.'）
此程序为一初步程序，假设初值为[10;10;10;10]，则第一次运算得t0=,lamda=,迭代后的x=NaN。
现用共轭梯度法求解另一函数
对上述程序稍加改动来求解本题的代码如下：
注：
程序所用语言为MATLAB，精度为
functions=gongegrad2（f,c,eps）
%c是初值，eps为允许误差值
ifnargin==2
%eps=;
elseifnargin<1
error（''）
return
end
tic
symsx1x2t0t1
grad1=diff（f,x1）;
grad2=diff（f,x2）;
grad=[grad1;grad2];
a=grad;
i=1;n=40;
gradk=subs（a,[x1x2],[c
（1）c
（2）]）;
x=tihuan（c）;
p0=0;
whilenorm（gradk）>=eps
p0=-gradk;
y=x;
x=x+t0*p0;
k=0;
gradk=subs（f,[x1x2],[x
（1）x
（2）]）;
w=solve（gradk）;
t0=real（w）;
t0=eval（t0）;
t0=t0
（1）;
x=y+t0*p0;
gradk=subs（a,[x1x2],[x
（1）x
（2）]）;
whilenorm（gradk）>=eps
ifk+1~=n
gradk2=subs（a,[x1x2],[x
（1）x
（2）]）;
gradk1=subs（a,[x1x2],[y
（1）y
（2）]）;
lamda=norm（gradk2）^2/norm（gradk1）^2;
p0=-gradk2+lamda*p0;
k=k+1;
else
k=0;
p0=-subs（a,[x1x2],[x
（1）x
（2）]）;
end
cleary;y=x;
x=x+t1*p0;
gradk=subs（f,[x1x2],[x
（1）x
（2）]）;
m=solve（gradk）;
t1=real（m）;t1=eval（t1
（1））;
x=y+t1*p0;
clearm;cleart1;symst1
gradk=subs（a,[x1x2],[x
（1）x
（2）]）;
end
disp（sprintf（'thelastpointwewantis[%f%f]',x
（1）,x
（2）））;
disp（sprintf（'thetimesusedtorecursionis%f',k））;
disp（sprintf（'thefunctionvalueis%f',x
（1）^2+25*x
（2）^2））;
toc
return;
end
disp（sprintf（'thelastpointwewantis[%f%f]',x
（1）,x
（2）））;
disp（sprintf（'thetimesusedtorecursionis%f',k））;
disp（sprintf（'thefunctionvalueis%f',x
（1）^2+25*x
（2）^2））;
toc
“tihuan”子函数为：
functionx=tihuan（x）
[v,g]=size（x）;
fori=1:
v
x（i）=x（i）;
end
程序调用方式为：
clearall
clc
symsx1x2t0t1
f=x1^2+25*x2^2;
c=[2;2];%初值
gongegrad2（f,c,eps）
程序结果如下：
由上述结果知，用共轭梯度法对上述函数求解需要迭代三次得到最优解0，此时x为[00];符合理论分析的结果，求解完成。
用外点法法求解函数
2.1.1问题分析
外点法为序列无约束最优化方法，其基本思想为将条件函数吸收到目标函数中进行求解。
其在数学上的直观理解是拉格朗日乘子法：
，
为总代价，
为价格，
为罚款。
即在经济学上总代价为价格和罚款的和。
此时
称
为增广目标函数，通常取
2.1.2问题求解
两种方法求解程序如下：
2.1.2.1注：
程序所用语言为MATLAB，终止条件为
，程序中无约束优化部分通过求导实现。
M文件如下：
tic
clc
%c是初值，eps为允许误差值
ifnargin==1
eps=;
elseifnargin<1
error（''）
return
end
symsm
[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;
iff1disp（sprintf（'thefinalxis[%f%f]',x1
（1）,x2
（1）））;
disp（sprintf（'thefinalfunctionvalueis%f',f1））;
else
disp（sprintf（'thefinalxis[%f%f]',x1
（2）,x2
（2）））;
disp（sprintf（'thefinalfunctionvalueis%f',f2））;
end
disp（sprintf（'thetimesusedtorecursionis%f',k））
toc
调用方式如下
symsx1x2m
T=x1^3+x2^2+m*（x1+x2-1）^2;
waidian（T,eps）;
functions=waidian（T,eps）
实验结果如下
由上述结果可知当迭代17次时能够达到终止条件，此时,
函数得到最优解.求解完成。
2.1.2.2，程序中无约束优化部分通过调用“牛顿法”完成
代码如下；
tic
clearall;closeall;clc
symsx1x2
m=1;x0=[2;2];eps=;
T=x1^3+x2^2+m*（x1+x2-1）^2;
c=10;k=1;
while1
s=newton1（T,x0,eps）;%调用newton法
x1=s
（1）;x2=s
（2）;
ifm*（x1+x2-1）^2>eps
m=c*m;k=k+1;symsx1x2
T=x1^3+x2^2+m*（x1+x2-1）^2;
else
disp（sprintf（'thefinalresultis[%f%f]',x1,x2））;
disp（sprintf（'thefunctionvalueis%f',x1^3+x2^2））;
disp（sprintf（'thetimesusedtorecursionis%f',k））;
break;
end
end
toc
实验结果如下：
由上述结果可知当迭代7次时能够达到终止条件，此时,
函数得到最优解.求解完成。
用内点法法求解函数
2.2.1问题分析
同外点法，内点法也为序列无约束最优化方法，其基本思想为将条件函数吸收到目标函数中进行求解。
称
为障碍函数、围墙函数或者惩罚函数。
通常取
若在上一次迭代中x处于可解区域外部，则“围墙”不够高，在下一次迭代时加高“围墙函数”。
2.2.2问题求解
求解程序如下：
注：
程序所用语言为MATLAB，本程序出于对运行时间和算法简单的考虑，选用
终止条件为
，此时可以认为当x接近求解边界时
（围墙）为无穷大，以确保函数的自变量迭代发生在可解区域内部
2.2.2.1，程序中无约束优化部分通过求导实现。
tic
clearall;
closeall;clc;
eps=10^6;
symswx1x2
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））;
whilebound<=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（'thefinalxis[%f%f]',x1,x2））;
disp（sprintf（'thefinalfunctionvalueis%f',x1
（1）^3+x2
（1）^2））;
disp（sprintf（'thetimesusedtorecursionis%f',k））
toc
实验结果如下：
由上述结果可知当迭代15次时能够达到终止条件，此时,
函数得到最优解.求解完成，结果与“外点法”相同。
2.2.2.2，程序中无约束优化部分通过调用“牛顿法”完成（
）
代码如下；
tic
clearall;closeall;clc
symsx1x2
w=10;x0=[2;2];eps=;
I=x1^3+x2^2+w*（1/（x1+x2-1）^2）;
c=;k=1;
while1
s=newton1（I,x0,eps）;%调用newton法
x1=s
（1）;x2=s
（2）;
ifsqrt（（x1-x0
（1））^2+（x2-x0
（2））^2）>eps
w=c*w;k=k+1;x0
（1）=x1;x0
（2）=x2;
symsx1x2
I=x1^3+x2^2+w*（1/（x1+x2-1）^2）;
else
disp（sprintf（'thefinalresultis[%f%f]',x1,x2））;
disp（sprintf（'thefunctionvalueis%f',x1^3+x2^2））;
disp（sprintf（'thetimesusedtorecursionis%f',k））;
break;
end
end
toc
实验结果如下：
由上述结果可知当迭代36次时能够达到终止条件，此时,
函数得到最优解.求解完成，结果与“外点法”十分近似。
三、此次实验的收获
此次试验是用优化算法求解一些函数最值。
通过此次实验，我的收获主要有以下几点；
1，函数求最值在以前的学习中基本都是从函数本身出发，研究函数的本身性质，进而得出理论最值，对于一些性质不明确的函数涉猎不多。
在这门课程中，我们更多的是从一个初始点出发，通过一些函数解集的迭代过程来求解函数的最值。
虽然有的算法也牵扯到函数的性质（函数是否可导，是否有Hesse矩阵等），但相对于以前的求解方式，对函数本身的要求已经薄弱很多，这就意味着可以对更多的函数运用算法求解最值，扩大了可求解函数的范围。
2，这门课程的学习多是理论为主，老师在课堂上讲了很多实效有用的算法，对算法的优缺点也多有说明，但没在实战中运用的话理解不会很深。
这次实验给了我一个运用所学理论知识的完美机会，在此次实验中能够通过自己的努力将课堂上讲的理论算法逐步运用到实战中，无疑对理论的理解是大有裨益的。
3，这次实验中另外一个比较重要的收获是学习了MATLAB的运用。
以前没有学习过MATLAB，但在这次试验中，能够自己学习MATLAB并且很拙劣的完成实验内容，这种收获对以后的科目学习也是很重要的。
升级会员 