§2 无约束问题
2.1 一维搜索方法
当用迭代法求函数的极小点时,常常用到一维搜索,即沿某一已知方向求目标函数的极小点.一维搜索的方法很多,常用的有:(1)试探法("成功—失败",斐波那契法,0.618法等);插值法(抛物线插值法,三次插值法等);(3)微积分中的求根法(切线法,二分法等).
考虑一维极小化问题
(2)
若是区间上的下单峰函数,我们介绍通过不断地缩短的长度,来搜索得(2)的近似最优解的两个方法.
为了缩短区间,逐步搜索得(2)的最优解的近似值,我们可以采用以下途径:在中任取两个关于是对称的点和(不妨设,并把它们叫做搜索点),计算和并比较它们的大小.对于单峰函数,若,则必有,因而是缩短了的单峰区间;若,则有,故是缩短了的单峰区间;若,则和都是缩短了的单峰.因此通过两个搜索点处目标函数值大小的比较,总可以获得缩短了的单峰区间.对于新的单峰区间重复上述做法,显然又可获得更短的单峰区间.如此进行,在单峰区间缩短到充分小时,我们可以取最后的搜索点作为(2)最优解的近似值.
应该按照怎样的规则来选取探索点,使给定的单峰区间的长度能尽快地缩短
Fibonacci法
若数列{}满足关系:
则称为Fibonacci数列,称为第个Fibonacci数,称相邻两个Fibonacci数之比为Fibonacci分数.
当用斐波那契法以个探索点来缩短某一区间时,区间长度的第一次缩短率为,其后各次分别为.由此,若和是单峰区间中第1个和第2个探索点的话,那么应有比例关系
,
从而
, (3)
它们关于确是对称的点.
如果要求经过一系列探索点搜索之后,使最后的探索点和最优解之间的距离不超过精度,这就要求最后区间的长度不超过,即
(4)
据此,我们应按照预先给定的精度,确定使(4)成立的最小整数作为搜索次数,直到进行到第个探索点时停止.
用上述不断缩短函数的单峰区间的办法,来求得问题(2)的近似解,是Kiefer(1953年)提出的,叫做Finbonacci法,具体步骤如下:
1° 选取初始数据,确定单峰区间,给出搜索精度,由(4)确定搜索次数.
2° ,计算最初两个搜索点,按(3)计算和.
3° while
if
else
end
end
4° 当进行至时,
这就无法借比较函数值和的大小确定最终区间,为此,取
其中为任意小的数.在和这两点中,以函数值较小者为近似极小点,相应的函数值为近似极小值.并得最终区间或.
由上述分析可知,斐波那契法使用对称搜索的方法,逐步缩短所考察的区间,它能以尽量少的函数求值次数,达到预定的某一缩短率.
例3 试用斐波那契法求函数的近似极小点,要求缩短后的区间不大于区间的0.08倍.
程序留作习题.
2.1.2 0.618法
若,满足比例关系
称之为黄金分割数,其值为.
黄金分割数和Fibonacci分数之间有着重要的关系,它们是
1° ,为偶数,
,为奇数.
2° .
现用不变的区间缩短率0.618,代替斐波那契法每次不同的缩短率,就得到了黄金分割法(0.618法).这个方法可以看成是斐波那契法的近似,实现起来比较容易,效果也相当好,因而易于为人们所接受.
用0.618法求解,从第2个探索点开始每增加一个探索点作一轮迭代以后,原单峰区间要缩短0.618倍.计算个探索点的函数值可以把原区间连续缩短次,因为每次的缩短率均为,故最后的区间长度为
这就是说,当已知缩短的相对精度为时,可用下式计算探索点个数:
当然,也可以不预先计算探索点的数目,而在计算过程中逐次加以判断,看是否已满足了提出的精度要求.
0.618法是一种等速对称进行试探的方法,每次的探索点均取在区间长度的0.618倍和0.382倍处.
2.2 二次插值法
对极小化问题(2),当在上连续时,可以考虑用多项式插值来进行一维搜索.它的基本思想是:在搜索区间中,不断用低次(通常不超过三次)多项式来近似目标函数,并逐步用插值多项式的极小点来逼近(2)的最优解.
2.3 无约束极值问题的解法
无约束极值问题可表述为
(5)
求解问题(5)的迭代法大体上分为两种:一是用到函数的一阶导数或二阶导数,称为解析法.另一是仅用到函数值,称为直接法.
2.3.1 解析法
2.3.1.1 梯度法(最速下降法)
对基本迭代格式
(6)
我们总是考虑从点出发沿哪一个方向,使目标函数下降得最快.微积分的知识告诉我们,点的负梯度方向
,
是从点出发使下降最快的方向.为此,称负梯度方向为在点处的最速下降方向.
按基本迭代格式(6),每一轮从点出发沿最速下降方向作一维搜索,来建立求解无约束极值问题的方法,称之为最速下降法.
这个方法的特点是,每轮的搜索方向都是目标函数在当前点下降最快的方向.同时,用或作为停止条件.其具体步骤如下:
1°选取初始数据.选取初始点,给定终止误差,令.
2°求梯度向量.计算, 若,停止迭代,输出.否则,进行3°.
3° 构造负梯度方向.取
.
4° 进行一维搜索.求,使得
令转2°.
例4 用最速下降法求解无约束非线性规划问题
其中,要求选取初始点,终止误差.
解:(i)
编写M文件detaf.m如下
function [f,df]=detaf(x);
f=x(1)^2+25*x(2)^2;
df(1)=2*x(1);
df(2)=50*x(2);
(ii)编写M文件zuisu.m
clc
x=[2;2];
[f0,g]=detaf(x);
while norm(g)>0.000001
p=-g''/norm(g);
t=1.0;f=detaf(x+t*p);
while f>f0
t=t/2;f=detaf(x+t*p);
end
x=x+t*p
[f0,g]=detaf(x)
end
2.3.1.2 Newton法
考虑目标函数在点处的二次逼近式
假定Hesse阵
正定.
由于正定,函数的稳定点是的最小点.为求此最小点,令
,
即可解得
.
对照基本迭代格式(1),可知从点出发沿搜索方向.
并取步长即可得的最小点.通常,把方向叫做从点出发的Newton方向.从一初始点开始,每一轮从当前迭代点出发,沿Newton方向并取步长为1的求解方法,称之为Newton法.其具体步骤如下:
1°选取初始数据.选取初始点,给定终止误差,令.
2°求梯度向量.计算,若,停止迭代,输出.否则,进行3°.
3°构造Newton方向.计算,取
.
4° 求下一迭代点.令,转2°.
例5 用Newton法求解,
选取,.
解:(i)
编写M文件nwfun.m如下:
function [f,df,d2f]=nwfun(x);
f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2;
df(1)=4*x(1)^3+2*x(1)*x(2)^2;
df(2)=100*x(2)^3+2*x(1)^2*x(2);
d2f(1,1)=12*x(1)^2+2*x(2)^2;
d2f(1,2)=4*x(1)*x(2);
d2f(2,1)=d2f(1,2);
d2f(2,2)=300*x(2)^2+4*x(1)*x(2);
(ii)编写M文件:
clc
x=[2;2];
[f0,g1,g2]=nwfun(x)
while norm(g1)>0.00001 %dead loop,for i=1:3
p=-inv(g2)*g1'',p=p/norm(p)
t=1.0,f=detaf(x+t*p)
while f>f0
t=t/2,f=detaf(x+t*p),
end
x=x+t*p
[f0,g1,g2]=nwfun(x)
end
如果目标函数是非二次函数,一般地说,用Newton法通过有限轮迭代并不能保证可求得其最优解.
Newton法的优点是收敛速度快;缺点是有时不好用而需采取改进措施,此外,当维数较高时,计算的工作量很大.
2.3.1.3 变尺度法
变尺度法(Variable Metric Algorithm)是近20多年来发展起来的,它不仅是求解无约束极值问题非常有效的算法,而且也已被推广用来求解约束极值问题.由于它既避免了计算二阶导数矩阵及其求逆过程,又比梯度法的收敛速度快,特别是对高维问题具有显著的优越性,因而使变尺度法获得了很高的声誉.下面我们就来简要地介绍一种变尺度法—DFP法的基本原理及其计算过程.这一方法首先由Davidon在1959年提出,后经Fletcher和Powell加以改进.
我们已经知道,牛顿法的搜索方向是,为了不计算二阶导数矩阵及其逆阵,我们设法构造另一个矩阵,用它来逼近二阶导数矩阵的逆阵,这一类方法也称拟牛顿法(Quasi-Newton Method).
下面研究如何构造这样的近似矩阵,并将它记为.我们要求:每一步都能以现有的信息来确定下一个搜索方向;每做一次选代,目标函数值均有所下降;这些近似矩阵最后应收敛于解点处的Hesse阵的逆阵.
当是二次函数时,其Hesse阵为常数阵,任两点和处的梯度之差为
或
对于非二次函数,仿照二次函数的情形,要求其Hesse阵的逆阵的第次近似矩阵满足关系式
(7)
这就是常说的拟Newton条件.
若令
(8)
则式(7)变为
, (9)
现假定已知,用下式求(设和均为对称正定阵);
(10)
其中称为第次校正矩阵.显然,应满足拟Newton条件(9),即要求
或
(11)
由此可以设想, 的一种比较简单的形式是
(12)
其中和为两个待定列向量.
将式(12)中的代入(11),得
这说明,应使
(13)
考虑到应为对称阵,最简单的办法就是取
(14)
由式(13)得
(15)
若和不等于零,则有
(16)
于是,得校正矩阵
(17)
从而得到
(18)
上述矩阵称为尺度矩阵.通常,我们取第一个尺度矩阵为单位阵,以后的尺度矩阵按式(18)逐步形成.可以证明:
(i)当不是极小点且正定时,式(17)右端两项的分母不为零,从而可按式(18)产生下一个尺度矩阵;
(ii)若为对称正定阵,则由式(18)产生的也是对称正定阵;
(iii)由此推出DFP法的搜索方向为下降方向.
现将DFP变尺度法的计算步骤总结如下.
1°给定初始点及梯度允许误差.
2°若,则即为近似极小点,停止迭代,否则,转向下一步.
3°令
(单位矩阵),
在方向进行一维搜索,确定最佳步长:
如此可得下一个近似点
4°一般地,设已得到近似点,算出,若
则即为所求的近似解,停止迭代;否则,计算:
并令,在方向上进行一维搜索,得,从而可得下一个近似点
5°若满足精度要求,则即为所求的近似解,否则,转回4°,直到求出某点满足精度要求为止.
2.3.2 直接法
在无约束非线性规划方法中,遇到问题的目标函数不可导或导函数的解析式难以表示时,人们一般需要使用直接搜索方法.同时,由于这些方法一般都比较直观和易于理解,因而在实际应用中常为人们所采用.下面我们介绍Powell方法.
这个方法主要由所谓基本搜索,加速搜索和调整搜索方向三部分组成,具体步骤如下:
1° 选取初始数据.选取初始点,个线性无关初始方向,组成初搜索方向组.给定终止误差,令.
2°进行基本搜索.令,依次沿中的方向进行一维搜索.对应地得到辅助迭代点,即
3°构造加速方向.令,若,停止迭代,输出.否则进行4°.
4°确定调整方向.按下式
找出.若
成立,进行5°.否则,进行6°.
5°调整搜索方向组.令
.
同时,令
,
,转2°.
6°不调整搜索方向组.令,转2°.
2.4 Matlab求函数的极小值和函数的零点
2.4.1 求单变量有界非线性函数在区间上的极小值
Matlab的命令为
[X,FVAL] = FMINBND(FUN,x1,x2,OPTIONS),
它的返回值是极小点和函数的极小值.这里fun 是用M文件定义的函数或Matlab中的单变量数学函数.
例6 求函数 的最小值.
解 编写M文件fun1.m
function f=fun1(x);
f=(x-3)^2-1;
在Matlab的命令窗口输入
[x,y]=fminbnd(''fun1'',0,5)
即可求得极小点和极小值.
求多变量函数的极小值
其中是一个向量,是一个标量函数.
Matlab中的基本命令是
[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2, ...)
它的返回值是向量的值和函数的极小值.FUN 是一个M文件,当FUN只有一个返回值时,它的返回值是函数;当FUN有两个返回值时,它的第二个返回值是的一阶导数行向量;当FUN有三个返回值时,它的第三个返回值是的二阶导数阵(Hessian阵).X0是向量的初始值,OPTIONS是优化参数,使用确省参数时,OPTIONS为空矩阵.P1,P2是可以传递给FUN的一些参数.
例7 求函数的最小值.
解:编写M文件fun2.m如下:
function [f,g]=fun2(x);
f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
g=[-400*x(1)*(x(2)-x(1)^2)-2(1-x(1)) 200*(x(2)-x(1)^2)];
在Matlab命令窗口输入
fminunc(''fun2'',rand(1,2))
即可求得函数的极小值.
求多元函数的极值也可以使用Matlab的命令
[X,FVAL]= FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...).
§3 约束极值问题
带有约束条件的极值问题称为约束极值问题,也叫约束规划问题.
求解约束极值问题要比求解无约束极值问题困难得多.为了简化其优化工作,可采用以下方法:将约束问题化为无约束问题;将非线性规划问题化为线性规划问题,以及能将复杂问题变换为较简单问题的其它方法.
最优性条件
库恩—塔克条件是非线性规划领域中最重要的理论成果之一,是确定某点为最优点的必要条件,但一般说它并不是充分条件(对于凸规划,它既是最优点存在的必要条件,同时也是充分条件).
二次规划
若某非线性规划的目标函数为自变量的二次函数,约束条件又全是线性的,就称这种规划为二次规划.
Matlab中二次规划的数学模型可表述如下:
这里是实对称矩阵,是列向量,是相应维数的矩阵.
Matlab中求解二次规划的命令是
[X,FVAL]= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
X的返回值是向量,FVAL的返回值是目标函数在X处的值.(具体细节可以参看在Matlab指令中运行help quadprog后的帮助).
例8 求解二次规划
解 编写如下程序:
h=[4,-4;-4,8];
f=[-6;-3];
a=[1,1;4,1];
b=[3;9];
[x,value]=quadprog(h,f,a,b,[],[],zeros(2,1))
求得
.
罚函数法
利用罚函数法,可将非线性规划问题的求解,转化为求解一系列无约束极值问题,因而也称这种方法为序列无约束最小化技术,简记为 SUMT (Sequential Unconstrained Minimization Technique).
罚函数法求解非线性规划问题的思想是,利用问题中的约束函数作出适当的罚函数,由此构造出带参数的增广目标函数,把问题转化为无约束非线性规划问题.主要有两种形式,一种叫外罚函数法,另一种叫内罚函数法,下面介绍外罚函数法.
考虑如下问题:
s.t.
取一个充分大的数 ,构造函数
(或
这里 ,,,为适当的行向量,Matlab中可以直接利用 和 函数.)则以增广目标函数为目标函数的无约束极值问题
的最优解也是原问题的最优解.
例9 求下列非线性规划
解 (i)编写 M 文件 test.m
function g=test(x);
M=50000;
f=x(1)^2+x(2)^2+8;
g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)^2-x(2),0)...
+M*abs(-x(1)-x(2)^2+2);
(ii)在Matlab命令窗口输入
[x,y]=fminunc(''test'',rand(2,1))
即可求得问题的解.
§4 飞行管理问题
在约10,000m高空的某边长160km的正方形区域内,经常有若干架飞机作水平飞行.区域内每架飞机的位置和速度向量均由计算机记录其数据,以便进行飞行管理.当一架欲进入该区域的飞机到达区域边缘时,记录其数据后,要立即计算并判断是否会与区域内的飞机发生碰撞.如果会碰撞,则应计算如何调整各架(包括新进入的)飞机飞行的方向角,以避免碰撞.现假定条件如下:
1)不碰撞的标准为任意两架飞机的距离大于8km;
2)飞机飞行方向角调整的幅度不应超过30度;
3)所有飞机飞行速度均为每小时800km;
4)进入该区域的飞机在到达区域边缘时,与区域内飞机的距离应在60km以上;
5)最多需考虑6架飞机;
6)不必考虑飞机离开此区域后的状况.
请你对这个避免碰撞的飞行管理问题建立数学模型,列出计算步骤,对以下数据进行计算(方向角误差不超过0.01度),要求飞机飞行方向角调整的幅度尽量小.
设该区域4个顶点的座标为(0,0),(160,0),(160,160),(0,160).记录数据为:
飞机编号 横座标 纵座标 方向角(度)
1 150 140 243
2 85 85 236
3 150 155 220.5
4 145 50 159
5 130 150 230
新进入 0 0 52
注:方向角指飞行方向与轴正向的夹角.
试根据实际应用背景对你的模型进行评价与推广.
提示:
,,
其中为飞机的总架数,为时刻第架飞机的坐标,分别表示第架飞机飞出正方形区域边界的时刻.这里
,,;
,,;
其中为飞机的速度,分别为第架飞机的初始方向角和调整后的方向角.
令
其中,
则两架飞机不碰撞的条件是.
习 题 三
1. 用Matlab的非线性规划命令fmincon求解飞行管理问题.
2. 用罚函数法求解飞行管理问题.
3. 某工厂向用户提供发动机,按合同规定,其交货数量和日期是:第一季度末交40台,第二季末交60台,第三季末交80台.工厂的最大生产能力为每季100台,每季的生产费用是(元),此处为该季生产发动机的台数.若工厂生产的多,多余的发动机可移到下季向用户交货,这样,工厂就需支付存贮费,每台发动机每季的存贮费为4元.问该厂每季应生产多少台发动机,才能既满足交货合同,又使工厂所花费的费用最少(假定第一季度开始时发动机无存货). |
发表于 2006-3-27 02:31
提升卡
变色卡
显身卡
楼主