第三章 非线性规划
§1 非线性规划
-
- 非线性规划的实例与定义
如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。
下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。
例 1 （投资决策问题）某企业有n 个项目可供选择投资，并且至少要对其中一个
项目投资。已知该企业拥有总资金 A 元，投资于第i(i = 1,L, n) 个项目需花资金ai 元， 并预计可收益bi 元。试选择最佳投资方案。
解 设投资决策变量为
⎧1,
xi = ⎨
决定投资第i个项目
， i = 1,L, n ，
⎩0， 决定不投资第i个项目
n n
则投资总额为åai xi ，投资总收益为åbi xi 。因为该公司至少要对一个项目投资，并
i =1 i =1
且总的投资金额不能超过总资金 A ，故有限制条件
0 < åai xi £ A
i =1
另外，由于 xi (i = 1,L, n) 只取值 0 或 1，所以还有
xi (1- xi ) = 0, i = 1,L, n.
|
n |
|
|
最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取 0 或 1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：
åbi xi
max Q = i = 1
åai xi
i =1
n
s.t.
0 < åai xi £ A
i =1
xi (1 - xi ) = 0, i = 1,L, n.
上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式
和 h j ( j = 1,L, q) 称为约束函数。另外， g i ( x) = 0 ( i = 1,L, p) 称为等式约束，h j ( x) £ 0 ( j = 1,L, q) 称为不等式的约束。
对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：
- 确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。
- 提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。
- 给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。
- 寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。
- 线性规划与非线性规划的区别
如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。
-
- 非线性规划的 Matlab 解法
Matlab 中非线性规划的数学模型写成以下形式
min f ( x)
⎧ Ax £ B
⎪ Aeq × x = Beq
,
⎪C( x) £ 0
⎪⎩Ceq( x) = 0
其中 f ( x) 是标量函数， A, B, Aeq, Beq 是相应维数的矩阵和向量，C( x), Ceq( x) 是非
线性向量函数。
Matlab 中的命令是X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
它的返回值是向量 x ，其中 FUN 是用M 文件定义的函数 f ( x) ；X0 是 x 的初始值；
A,B,Aeq,Beq 定义了线性约束 A * X £ B, Aeq * X = Beq ，如果没有线性约束，则
A=[],B=[],Aeq=[],Beq=[]；LB 和 UB 是变量 x 的下界和上界，如果上界和下界没有约束，则 LB=[]，UB=[]，如果 x 无下界，则 LB 的各分量都为-inf，如果 x 无上界，则 UB 的各分量都为 inf；NONLCON 是用M 文件定义的非线性向量函数C( x), Ceq( x) ；OPTIONS
定义了优化参数，可以使用 Matlab 缺省的参数设置。例2 求下列非线性规划
x + 2 x2 = 3x1, x2 , x3 ³ 0
解 （i）编写 M 文件fun1.m 定义目标函数
function f=fun1(x); f=sum(x.^2)+8;
- 编写M文件fun2.m定义非线性约束条件
function [g,h]=fun2(x); g=[-x(1)^2+x(2)-x(3)^2
x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束h=[-x(1)-x(2)^2+2
x(2)+2*x(3)^2-3]; %非线性等式约束
- 编 写 主 程 序 文 件 example2.m 如 下 ： options=optimset('largescale','off'); [x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ... 'fun2', options)
就可以求得当 x1 = 0.5522, x2 = 1.2033, x3 = 0.9478 时，最小值 y = 10.6511 。
-
- 求解非线性规划的基本迭代格式记（NP）的可行域为 K 。
若 x* Î K ，并且
则称 x* 是（NP）的严格局部最优解， f ( x* ) 是(NP)的严格局部最优值。
由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然，有时求出的某个解虽是一部分可行域上的极值点，但却并不一定是整个可行域上的全局最优解。
对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思想是：从一个选定的初始点 x0 Î Rn 出发，按照某一特定的迭代规则产生一个点列
{xk }，使得当{xk } 是有穷点列时，其最后一个点是(NP)的最优解；当{x k } 是无穷点列时，它有极限点，并且其极限点是(NP)的最优解。
设 xk Î Rn 是某迭代方法的第k 轮迭代点， xk +1 Î Rn 是第k + 1 轮迭代点，记
|
（1）
|
就是求解非线性规划模型(NP)的基本迭代格式。
通常，我们把基本迭代格式（1）中的 p k 称为第 k 轮搜索方向， t 为沿 p k 方向的
步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步长。
设 x Î Rn , p ¹ 0 ，若存在d > 0 ，使
f ( x + tp) < f ( x ), "t Î(0,d ) ，
称向量 p 是 f 在点 x 处的下降方向。
设 x Î Rn , p ¹ 0 ，若存在t > 0 ，使
x + tp Î K ，
称向量 p 是点 x 处关于 K 的可行方向。
一个向量 p ，若既是函数 f 在点 x 处的下降方向，又是该点关于区域 K 的可行方向，称之为函数 f 在点 x 处关于 K 的可行下降方向。
现在，我们给出用基本迭代格式（1）求解(NP)的一般步骤如下：
0° 选取初始点 x0 ，令k := 0 。
1° 构造搜索方向，依照一定规划，构造 f 在点 xk 处关于 K 的可行下降方向作为搜索方向 pk 。
|
3° 求出下一个迭代点。按迭代格式（1）求出
x k +1 = x k + t p k 。
若 xk +1 已满足某种终止条件，停止迭代。
4° 以 xk +1 代替 xk ，回到 1°步。
-
- 凸函数、凸规划
设 f ( x) 为定义在 n 维欧氏空间 E( n ) 中某个凸集 R 上的函数，若对任何实数
a (0 < a < 1) 以及 R 中的任意两点 x(1) 和 x(2) ，恒有
f (ax(1) + (1 - a )x( 2) ) £ af ( x(1) ) + (1 - a ) f ( x( 2) )
则称 f ( x) 为定义在 R 上的凸函数。
若对每一个a (0 < a < 1) 和 x(1) ¹ x( 2) Î R 恒有
f (ax(1) + (1 - a )x( 2) ) < af ( x(1) ) + (1 - a ) f ( x( 2) )
则称 f ( x) 为定义在 R 上的严格凸函数。考虑非线性规划
⎧⎪min
⎨ xÎR
f ( x)
⎪⎩R = {x | g j ( x) £ 0, j = 1,2,L, l}
假定其中 f ( x) 为凸函数，g j ( x)( j = 1,2,L,l) 为凸函数，这样的非线性规划称为凸规划。
可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优
解的集合形成一个凸集。当凸规划的目标函数 f ( x) 为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非
线性规划。
§2 无约束问题
-
- 一维搜索方法
当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法， 0.618 法等）；（2）插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线法，二分法等）。
考虑一维极小化问题
min
a £t £b
f (t)
（2）
若 f (t) 是[a, b] 区间上的下单峰函数，我们介绍通过不断地缩短[a, b] 的长度，来搜索得（2）的近似最优解的两个方法。
为了缩短区间[a, b] ，逐步搜索得（2）的最优解t* 的近似值，我们可以采用以下途径：在[a, b] 中任取两个关于[a, b] 是对称的点t1 和t2 （不妨设t2 < t1 ，并把它们叫
做搜索点），计算 f (t1 ) 和 f (t2 ) 并比较它们的大小。对于单峰函数，若 f (t2 ) < f (t1 ) ，
则必有 t* Î[a, t ] ，因而 [a, t ] 是缩短了的单峰区间；若 f (t ) < f (t
) ，则有
1 1 1 2
t* Î[t
, b] ，故[t , b] 是缩短了的单峰区间；若 f (t ) = f (t ) ，则[a, t ] 和[t
, b] 都是
2 2 2 1 1 2
缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行， 在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（2）最优解的近似值。
应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？
-
-
- Fibonacci 法
-
如用 Fn 表示计算n 个函数值能缩短为单位长区间的最大原区间长度，可推出 Fn 满足关系
F0 = F1 = 1
Fn = Fn -2 + Fn -1 ,
n = 2,3,L,
数列{Fn } 称为 Fibonacci 数列， Fn 称为第n 个 Fibonacci 数，相邻两个 Fibonacci 数
Fn
当用斐波那契法以 n 个探索点来缩短某一区间时，区间长度的第一次缩短率为
Fn -1 ，其后各次分别为 Fn-2 , Fn-3 ,L, F1 。由此，若t 和t
< t ) 是单峰区间[a, b]
Fn Fn-1
Fn-2 F2
1 2 2 1
中第 1 个和第 2 个探索点的话，那么应有比例关系
t1 - a = Fn -1 ， t2 - a =
（3）
它们关于[a, b] 确是对称的点。
如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度d > 0 ，这就要求最后区间的长度不超过d ，即
b - a £ d
Fn
（4）
据此，我们应按照预先给定的精度d ，确定使（4）成立的最小整数n 作为搜索次数， 直到进行到第n 个探索点时停止。
用上述不断缩短函数 f (t) 的单峰区间[a, b] 的办法，来求得问题（2）的近似解， 是 Kiefer(1953 年)提出的，叫做 Finbonacci 法，具体步骤如下：
1o 选取初始数据，确定单峰区间[a0 , b0 ] ，给出搜索精度d > 0 ，由（4）确定搜
索次数n 。
=
|
其中 为任意小的数。在t1 和t2 这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小值。并得最终区间[a, t1 ] 或[t2 , b] 。
由上述分析可知，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，它能
以尽量少的函数求值次数，达到预定的某一缩短率。
例 3 试用斐波那契法求函数 f (t) = t 2 - t + 2 的近似极小点，要求缩短后的区间不大于区间[-1,3] 的 0.08 倍。
程序留作习题。
2.1.2 0.618 法
若w > 0 ，满足比例关系
w = lim Fn -1 。
n®¥ Fn
现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（0.618 法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果也相当好，因而易于为人们所接受。
用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短 0.618 倍。计算n 个探索点的函数值可以把原区间[a0 , b0 ] 连续缩短n -1 次，因为每次的缩短率均为 ，故最后的区间长度为
(b0
- a0
)m n -1
这就是说，当已知缩短的相对精度为d 时，可用下式计算探索点个数n ：
m n -1 £ d
当然，也可以不预先计算探索点的数目n ，而在计算过程中逐次加以判断，看是否已满足了提出的精度要求。
0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618
倍和 0.382 倍处。
-
- 二次插值法
对极小化问题（2），当 f (t) 在[a, b] 上连续时，可以考虑用多项式插值来进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。
-
- 无约束极值问题的解法无约束极值问题可表述为
min
f ( x),
x Î E ( n )
（5）
求解问题（5）的迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数， 称为解析法。另一是仅用到函数值，称为直接法。
-
-
- 解析法
- 梯度法（最速下降法） 对基本迭代格式
- 解析法
-
|
（6）
我们总是考虑从点 xk 出发沿哪一个方向 pk ，使目标函数 f 下降得最快。微积分的知识告诉我们，点 xk 的负梯度方向
pk = -Ñf ( xk ) ，
是从点 xk 出发使 f 下降最快的方向。为此，称负梯度方向- Ñf ( xk ) 为 f 在点 xk 处的最速下降方向。
按基本迭代格式（6），每一轮从点 xk 出发沿最速下降方向- Ñf ( xk ) 作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。
这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同
时，用Ñf ( xk ) = 0 或 Ñf ( xk ) £ e 作为停止条件。其具体步骤如下：
1°选取初始数据。选取初始点 x0 ，给定终止误差，令k := 0 。
2°求梯度向量。计算Ñf ( xk ) ， 若 Ñf (xk )
进行 3°。
3° 构造负梯度方向。取
pk = -Ñf ( xk ) .
4° 进行一维搜索。求tk ，使得
£ e ，停止迭代，输出 xk 。否则，
|
pk ) = min f ( xk + tpk )
t ³0
|
例 4 用最速下降法求解无约束非线性规划问题
min f ( x) = x 2 + 25x2
1 2
|
|
编写 M 文件detaf.m，定义函数 f (x) 及其梯度列向量如下
function [f,df]=detaf(x); f=x(1)^2+25*x(2)^2; df=[2*x(1)
50*x(2)];
（ii）编写主程序文件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
x,f0
-
-
-
- Newton 法
-
-
由于Ñ2 f ( xk ) 正定，函数Q 的驻点 xk +1 是Q( x) 的极小点。为求此极小点，令
ÑQ( xk +1 ) = Ñf ( xk ) + Ñ2 f ( xk )( xk +1 - xk ) = 0 ， 即可解得
xk +1 = xk -[Ñ2 f ( xk )]-1 Ñf ( xk ) .
对照基本迭代格式（1），可知从点 xk 出发沿搜索方向。
pk = -[Ñ2 f ( xk )]-1 Ñf ( xk )
|
1°选取初始数据。选取初始点 x0 ，给定终止误差e > 0 ，令k := 0 。
2°求梯度向量。计算Ñf ( xk ) ，若 Ñf ( xk )
进行 3°。
£ e ，停止迭代，输出 xk 。否则，
3°构造 Newton 方向。计算[Ñ2 f ( xk )]-1 ，取
pk = -[Ñ2 f ( xk )]-1 Ñf ( xk ) .
4° 求下一迭代点。令 xk +1 = xk + pk , k := k + 1 ，转 2°。例 5 用 Newton 法求解，
min f ( x) = x 4 + 25x 4 + x 2 x 2
1
选取 x0 = (2,2)T 。
2 1 2
解：（i） Ñf ( x) = [4x3 + 2x x 2 100x 3 + 2x 2 x ]T
1 1 2 2 1 2
⎡12x2 + 2x2 4x x ⎤
Ñ2 f = ⎢ 1 2
1 2 ⎥
4x x
300x2 + 2x2
⎣ 1 2 2 1 ⎦
clc
（ii） 编 写 M 文 件 nwfun.m 如 下 ： function [f,df,d2f]=nwfun(x); f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2;
df=[4*x(1)^3+2*x(1)*x(2)^2;100*x(2)^3+2*x(1)^2*x(2)]; d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2)
4*x(1)*x(2),300*x(2)^2+2*x(1)^2];
（III）编写主程序文件 example5.m 如下：
x=[2;2];
[f0,g1,g2]=nwfun(x); while norm(g1)>0.00001
p=-inv(g2)*g1; x=x+p;
[f0,g1,g2]=nwfun(x); end
x, f0
如果目标函数是非二次函数，一般地说，用 Newton 法通过有限轮迭代并不能保证可求得其最优解。
为了提高计算精度，我们在迭代时可以采用变步长计算上述问题，编写主程序文件
example5_2 如下： clc,clear x=[2;2];
[f0,g1,g2]=nwfun(x); while norm(g1)>0.00001
p=-inv(g2)*g1;p=p/norm(p); t=1.0;f=nwfun(x+t*p); while f>f0
t=t/2;f=nwfun(x+t*p); end
x=x+t*p; [f0,g1,g2]=nwfun(x); end
x,f0
Newton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算-[Ñ2 f ( xk )]-1 的工作量很大。
-
-
-
- 变尺度法
-
-
变尺度法（Variable Metric Algorithm）是近 20 多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变尺度法—DFP 法的基本原理及其计算过程。这一方法首先由 Davidon 在 1959 年提出， 后经 Fletcher 和 Powell 加以改进。
我们已经知道，牛顿法的搜索方向是- [Ñ2 f ( xk )]-1 Ñf ( xk )
，为了不计算二阶
导数矩阵[Ñ2 f ( xk )] 及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵的逆阵[Ñ2 f ( xk )]-1 ，这一类方法也称拟牛顿法（Quasi-Newton Method）。
下面研究如何构造这样的近似矩阵，并将它记为 H ( k ) 。我们要求：每一步都能以现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似矩阵最后应收敛于解点处的 Hesse 阵的逆阵。
当 f ( x) 是二次函数时，其 Hesse 阵为常数阵 A ，任两点 xk 和 xk +1 处的梯度之差为
Ñf ( xk +1 ) - Ñf ( xk ) = A( xk +1 - xk )
或
xk +1 - xk = A-1[Ñf ( xk +1 ) - Ñf ( xk )]
对于非二次函数，仿照二次函数的情形，要求其 Hesse 阵的逆阵的第k + 1 次近似
并令 p k = - H ( k ) Ñ f ( x k ) ，在 p k 方向上进行一维搜索，得 l ，从而可得下一个近似点
|
5°若 xk +1 满足精度要求，则 xk +1 即为所求的近似解，否则，转回 4°，直到求出某点满足精度要求为止。
-
-
- 直接法
-
在无约束非线性规划方法中，遇到问题的目标函数不可导或导函数的解析式难以表示时，人们一般需要使用直接搜索方法。同时，由于这些方法一般都比较直观和易于理解，因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。
这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成，具体步骤如下：
1° 选取初始数据。选取初始点 x0 ， n 个线性无关初始方向，组成初搜索方向组
{p0 , p1 ,L, pn -1}。给定终止误差e > 0 ，令k := 0 。
2°进行基本搜索。令 y 0 := xk ，依次沿{p0 , p1 ,L, pn -1} 中的方向进行一维搜索。对应地得到辅助迭代点 y1 , y 2 ,L, yn ，即
|
k := k + 1 ，转 2°。
6°不调整搜索方向组。令 xk +1 := yn , k := k + 1，转 2°。
-
- Matlab 求无约束极值问题
在 Matlab 工具箱中，用于求解无约束极值问题的函数有 fminunc 和 fminsearch，用法介绍如下。
求函数的极小值
min f (x) ,
x
其中 x 可以为标量或向量。
Matlab 中 fminunc 的基本命令是
[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2, ...)
其中的返回值 X 是所求得的极小点，FVAL 是函数的极小值，其它返回值的含义参见相关的帮助。FUN 是一个 M 文件，当 FUN 只有一个返回值时，它的返回值是函数 f ( x) ； 当 FUN 有两个返回值时，它的第二个返回值是 f ( x) 的梯度向量；当 FUN 有三个返回值时，它的第三个返回值是 f ( x) 的二阶导数阵（Hessian 阵）。X0 是向量 x 的初始值，
OPTIONS 是优化参数，可以使用缺省参数。P1，P2 是可以传递给 FUN 的一些参数。
例 6 求函数 f (x) = 100(x - x 2 )2 + (1 - x )2 的最小值。
2 1 1
解：编写 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)];
编写主函数文件example6.m如下：
options = optimset('GradObj','on'); [x,y]=fminunc('fun2',rand(1,2),options)
即可求得函数的极小值。
在求极值时，也可以利用二阶导数，编写 M 文件 fun3.m 如下：
function [f,df,d2f]=fun3(x); f=100*(x(2)-x(1)^2)^2+(1-x(1))^2;
df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1)
-400*x(1),200];
编写主函数文件example62.m如下：
options = optimset('GradObj','on','Hessian','on'); [x,y]=fminunc('fun3',rand(1,2),options)
即可求得函数的极小值。
求多元函数的极值也可以使用 Matlab 的 fminsearch 命令，其使用格式为：
[X,FVAL,EXITFLAG,OUTPUT]＝FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...)
例 7 求函数 f (x) = sin(x) + 3 取最小值时的 x 值。解 编写 f (x) 的 M 文件 fun4.m 如下：
function f=fun4(x); f=sin(x)+3;
编写主函数文件example7.m如下：
x0=2;
[x,y]=fminsearch(@fun4,x0)
即求得在初值 2 附近的极小点及极小值。
§3 约束极值问题
带有约束条件的极值问题称为约束极值问题，也叫规划问题。
求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作，可采用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题变换为较简单问题的其它方法。
库恩—塔克条件是非线性规划领域中最重要的理论成果之一，是确定某点为最优点的必要条件，但一般说它并不是充分条件（对于凸规划，它既是最优点存在的必要条件， 同时也是充分条件）。
-
- 二次规划
若某非线性规划的目标函数为自变量 x 的二次函数，约束条件又全是线性的，就称这种规划为二次规划。
Matlab 中二次规划的数学模型可表述如下：
⎨ Aeq × x = beq这里 H 是实对称矩阵， f , b 是列向量， A 是相应维数的矩阵。
Matlab 中求解二次规划的命令是
[X,FVAL]= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)
返回值 X 是决策向量 x 的值，返回值 FVAL 是目标函数在 x 处的值。（具体细节可以参看在 Matlab 指令中运行 help quadprog 后的帮助）。
例 8 求解二次规划
⎧ min f (x) = 2x 2 - 4x x
+ 4x 2 - 6x
- 3x
⎪ 1
⎪ x1 + x2 £ 3
1 2 2 1 2
⎨
⎪ 4 x1
+ x2 £ 9
⎩⎪ x1, x2 ³ 0
解 编写如下程序：
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))
求得
|
1.0500
Min f ( x) = -11.0250 。
⎣ ⎦
-
- 罚函数法
利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题， 因而也称这种方法为序列无约束最小化技术，简记为 SUMT (Sequential Unconstrained Minization Technique)。
罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。
考虑问题：
min f ( x)
⎧gi (x) £ 0, i = 1,L, r,
⎪
s.t. ⎨hj (x) ³ 0, j = 1,L, s,
|
|
取一个充分大的数 M > 0 ，构造函数
r s t
P( x, M ) = f ( x) + M åmax(gi ( x),0) - M åmin(hi ( x),0) +M å| ki ( x) |
i=1 i=1 i=1
⎛ ⎛ G(x) ⎞⎞ ⎛ ⎛ H (x) ⎞⎞
（或 P(x, M ) =
f (x) + Msum⎜max⎜
0
⎟⎟ - Msum⎜min⎜
0
⎟⎟ + M || K (x) ||
⎝ ⎝ ⎠⎠ ⎝ ⎝ ⎠⎠
这里G(x) = [g1 (x) L gr (x)]， H (x) = [h1 (x) L hs (x)]，
K (x) = [k1 (x) L kt (t)]，Matlab 中可以直接利用 max 、min 和 sum 函数。）则以
增广目标函数 P( x, M ) 为目标函数的无约束极值问题
min P( x, M )
的最优解 x 也是原问题的最优解。例 9 求下列非线性规划
⎧min f (x) = x 2 + x 2 + 8
⎪ 1 2
|
⎨ 2
2
|
|
⎩ x1 , x2 ³ 0.
解 (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);
或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下：
function g=test(x); M=50000; f=x(1)^2+x(2)^2+8;
g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-x(2),0)+...
M*abs(-x(1)-x(2)^2+2);
我们也可以修改罚函数的定义，编写test.m如下：
function g=test(x); M=50000; f=x(1)^2+x(2)^2+8;
g=f-M*min(min(x),0)-M*min(x(1)^2-x(2),0)+M*(-x(1)-x(2)^2+2)^2;
(ii)在 Matlab 命令窗口输入
[x,y]=fminunc('test',rand(2,1))
即可求得问题的解。
-
- Matlab 求约束极值问题
在 Matlab 优化工具箱中，用于求解约束最优化问题的函数有：fminbnd、fmincon、
quadprog、fseminf、fminimax，上面我们已经介绍了函数 fmincon 和 quadprog。
-
-
- fminbnd 函数
-
求单变量非线性函数在区间上的极小值
min f (x) ,
x
Matlab 的命令为
x Î[ x1 , x2 ]
[X,FVAL] = FMINBND(FUN,x1,x2,OPTIONS)，
它的返回值是极小点 x 和函数的极小值。这里 fun 是用 M 文件定义的函数或 Matlab 中的单变量数学函数。
例 10 求函数
f (x) = (x - 3)2 - 1, x Î[0,5]
的最小值。
解 编写M 文件 fun5.m function f=fun5(x);
f=(x-3)^2-1;
在 Matlab 的命令窗口输入
[x,y]=fminbnd('fun5',0,5)
即可求得极小点和极小值。
-
-
- fseminf 函数求
-
min {F (x) | C(x) £ 0, Ceq(x) = 0, PHI (x, w) £ 0}
x
s.t.
⎧ A * x £ B
|
其中C(x), Ceq(x), PHI (x, w) 都是向量函数； w 是附加的向量变量， w 的每个分量都限定在某个区间内。
上述问题的 Matlab 命令格式为
X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq)
其中 FUN 用于定义目标函数 F (x) ；X0 为 x 的初始值；NTHETA 是半无穷约束
PHI (x, w) 的个数；函数 SEMINFCON 用于定义非线性不等式约束C(x) ，非线性等
|
|
例 11 求函数 f (x) = (x1
- 0.5)2 + (x
- 0.5)2 + (x
- 0.5)2 取最小值时的 x 值,
约束为：
K (x, w ) = sin(w x ) cos(w x ) -
1 (w
- 50)2 - sin(w x
) - x £ 1
1 1 1 1
1 2 1000 1
1 3 3
K 2 (x, w2
) = sin(w2 x2
) cos(w2
x1 ) -
1
1000(w2
- 50)2 - sin(w
x3 ) - x3 £ 1
1 £ w1 £ 100 ，1 £ w2 £ 100
解 （1）编写 M 文件 fun6.m 定义目标函数如下：
function f=fun6(x,s); f=sum((x-0.5).^2);
- 编写M 文件 fun7.m 定义约束条件如下：
function [c,ceq,k1,k2,s]=fun7(x,s); c=[];ceq=[];
if isnan(s(1,1))
s=[0.2,0;0.2 0];
end
%取样值w1=1:s(1,1):100;
w2=1:s(2,1):100;
%半无穷约束
k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1;
k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形plot(w1,k1,'-',w2,k2,'+');
- 调用函数 fseminf
在 Matlab 的命令窗口输入
[x,y]=fseminf(@fun6,rand(3,1),2,@fun7)
即可。
-
-
- fminimax 函数求解min ⎧⎨max F (x)⎫⎬
-
x ⎩ Fi ⎭
s.t.
⎧ A * x £ b
|
⎨C(x) £ 0
|
⎪⎩LB £ x £ UB
其中 F (x) = {F1 (x),L, Fm (x)}。上述问题的 Matlab 命令为
X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)
例12 求函数族{ f1 (x), f 2 (x), f3 (x), f 4 (x), f5 (x)} 取极大极小值时的 x 值。其中：
⎧ f (x) = 2x 2 + x 2 - 48x - 40x
+ 304
⎪ 1 1 2 1 2
⎪ f (x) = -x 2 - 3x 2
⎪ 2 1 2
⎨ f3 (x) = x1 + 3x2 - 18
|
|
⎪⎩ f5 (x) = x1 + x2 - 8
解 （1）编写 M 文件 fun8.m 定义向量函数如下：
function f=fun8(x);
f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304
-x(1)^2-3*x(2)^2
x(1)+3*x(2)-18
-x(1)-x(2)
x(1)+x(2)-8];
（2） 调 用 函 数 fminimax [x,y]=fminimax(@fun8,rand(2,1))
-
-
- 利用梯度求解约束优化问题
-
例 13 已知函数 f (x) = ex1 (4x 2 + 2x 2 + 4x x
-
- 2x
+ 1) ，且满足非线性约束：
1
⎧x1 x2 - x1 - x2 £ -1.5
2 1 2 2
⎨x x ³ -10
⎩ 1 2
求min f (x)
x
分析：当使用梯度求解上述问题时，效率更高并且结果更准确。题目中目标函数的梯度为：
⎡ex1 (4x 2 + 2x 2 + 4x x + 8x + 6x + 1)⎤
⎢ 1 2 1 2 1 2 ⎥
⎢ex1 (4x + 4x + 2) ⎥
⎣ 1 2 ⎦
解 （1）编写 M 文件 fun9.m 定义目标函数及梯度函数： function [f,df]=fun9(x); f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);
df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1))*(4*x(2)
+4*x(1)+2)];
- 编写M 文件 fun10.m 定义约束条件及约束条件的梯度函数：
function [c,ceq,dc,dceq]=fun10(x); c=[x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10];
dc=[x(2)-1,-x(2);x(1)-1,-x(1)];
ceq=[];dceq=[];
- 调用函数 fmincon，编写主函数文件 example13.m 如下：
%采用标准算法options=optimset('largescale','off');
%采用梯度options=optimset(options,'GradObj','on','GradConstr','on'); [x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)
-
- Matlab 优化工具箱的用户图形界面解法
Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。
optimtool 可应用到所有优化问题的求解，计算结果可以输出到 Matlab 工作空间中。
|
|
图 1 优化问题用户图形界面解法示意图
例 14 用 optimtool 重新求解例 1。
利用例 1 已经定义好的函数 fun1 和fun2。在 Matlab 命令窗口运行 optimtool，就打开图形界面，如图 1 所示，填入有关的参数，未填入的参数取值为空或者为默认值，然后用鼠标点一下“start”按钮，就得到求解结果，再使用“file”菜单下的“Export to Workspace…”选项，把计算结果输出到 Matlab 工作空间中去。
§4 飞行管理问题
在约 10，000m 高空的某边长 160km 的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下：
- 不碰撞的标准为任意两架飞机的距离大于 8km;
- 飞机飞行方向角调整的幅度不应超过 30 度；
- 所有飞机飞行速度均为每小时 800km;
- 进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在 60km 以上；
- 最多需考虑 6 架飞机；
- 不必考虑飞机离开此区域后的状况。
请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过 0.01 度），要求飞机飞行方向角调整的幅度尽量小。
设该区域 4 个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。记录数据见表 1。
表 1 飞行记录数据
| 飞机编号 | 横座标 x | 纵座标 y | 方向角（度） |
| 1 | 150 | 140 | 243 |
| 2 | 85 | 85 | 236 |
| 3 | 150 | 155 | 220.5 |
| 4 | 145 | 50 | 159 |
| 5 | 130 | 150 | 230 |
| 新进入 | 0 | 0 | 52 |
注：方向角指飞行方向与 x 轴正向的夹角。为方便以后的讨论，我们引进如下记号： D 为飞行管理区域的边长；
Ù 为飞行管理区域，取直角坐标系使其为[0, D] ´[0, D] ；
a 为飞机飞行速度， a = 800 km/h；
(x0 , y0 ) 为第i 架飞机的初始位置；
i i
(xi (t), yi (t)) 为第i 架飞机在t 时刻的位置；
q 0 为第i 架飞机的原飞行方向角，即飞行方向与 x 轴夹角， 0 £ q 0 < 2p ；
i i
|
i 6 i 6
q = q 0 + Äq 为第i 架飞机调整后的飞行方向角。
i i i
-
- 模型一
根据相对运动的观点在考察两架飞机i 和 j 的飞行时，可以将飞机i 视为不动而飞机 j 以相对速度
vij = v j - vi = (a cosq j - a cosqi , a sinq j - a sinqi )
相对于飞机i 运动，对（19）式进行适当的化约可
（Äqi
£ 30o ，
i = 1,2,L,6
利用如下的程序：
clc,clear
x0=[150 85 150 145 130 0];
y0=[140 85 155 50 150 0];
q=[243 236 220.5 159 230 52];
xy0=[x0; y0];
d0=dist(xy0); %求矩阵各个列向量之间的距离d0(find(d0==0))=inf;
a0=asind(8./d0) %以度为单位的反函数xy1=x0+i*y0
xy2=exp(i*q*pi/180) for m=1:6
for n=1:6
if n~=m
b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n))); end
end
end
b0=b0*180/pi;
dlmwrite('txt1.txt',a0,'delimiter', '\t','newline','PC'); fid=fopen('txt1.txt','a');
fwrite(fid,'~','char'); %往纯文本文件中写 LINGO 数据的分割符
dlmwrite('txt1.txt',b0,'delimiter', '\t','newline','PC','-append','roffset', 1)
|
|
| 1 | 2 | 3 | 4 | 5 | 6 | |
| 1 | 0 | 5.39119 | 32.23095 | 5.091816 | 20.96336 | 2.234507 |
| 2 | 5.39119 | 0 | 4.804024 | 6.61346 | 5.807866 | 3.815925 |
| 3 | 32.23095 | 4.804024 | 0 | 4.364672 | 22.83365 | 2.125539 |
| 4 | 5.091816 | 6.61346 | 4.364672 | 0 | 4.537692 | 2.989819 |
| 5 | 20.96336 | 5.807866 | 22.83365 | 4.537692 | 0 | 2.309841 |
| 6 | 2.234507 | 3.815925 | 2.125539 | 2.989819 | 2.309841 | 0 |
|
|
| 1 | 2 | 3 | 4 | 5 | 6 | |
| 1 | 0 | 109.26 | -128.25 | 24.18 | 173.07 | 14.475 |
| 2 | 109.26 | 0 | -88.871 | -42.244 | -92.305 | 9 |
| 3 | -128.25 | -88.871 | 0 | 12.476 | -58.786 | 0.31081 |
| 4 | 24.18 | -42.244 | 12.476 | 0 | 5.9692 | -3.5256 |
| 5 | 173.07 | -92.305 | -58.786 | 5.9692 | 0 | 1.9144 |
| 6 | 14.475 | 9 | 0.31081 | -3.5256 | 1.9144 | 0 |
上述飞行管理的数学规划模型可如下输入 LINGO 求解：
model:
sets:
plane/1..6/:delta; link(plane,plane):alpha,beta; endsets
data:
alpha=@file('txt1.txt'); !需要在alpha的数据后面加上分隔符"~"; beta=@file('txt1.txt');
enddata min=@sum(plane:@abs(delta)); @for(plane:@bnd(-30,delta,30));
@for(link(i,j)|i#ne#j:@abs(beta(i,j)+0.5*delta(i)+0.5*delta(j))>a lpha(i,j));
end
求得的最优解为Äq3 = 2.83858o ， Äq5 = -21.0138o ， Äq6 = 0.7908o ，其它调整角度为 0。
-
- 模型二
两架飞机i, j 不发生碰撞的条件为
|
（25）
1 £ i £ 5 ， i + 1 £
j £ 6 ， 0 £ t £ min{Ti ,Tj }
其中Ti ,Tj 分别表示第i, j 架飞机飞出正方形区域边界的时刻。这里
|
|
， y (t) = y0 + at sinq
， i = 1,2,L, n ；
|
|
|
|
i i i i 6
|
|
|
li, j
= (xi (t) - x j
(t))2 + ( y (t) - y
(t))2 - 64 = a~(i, j)t 2 + ~
+ c~(i, j)
其中 a~(i, j) = 4a2 sin 2 qi - q j ，
|
b (i, j) = 2a[(xi (0) - x j (0))(cosqi - cosq j ) + ( yi (0) - y j (0))(sinqi - sin q j )]
|
则两架i, j 飞机不碰撞的条件是
|
2 - 4a~(i, j)c~(i, j) < 0
（26）
这样我们建立如下的非线性规划模型
|
2
i
s.t.
i=1
Ä(i, j) < 0 ，1 £ i £ 5 ，
p
i + 1 £ j £ 6
- 用最速下降法（梯度法）求函数：
f (x) = 4x + 6x - 2x 2 - 2x x - 2x 2
1 2 1 1 2 2
的极大点。给定初始点 x0 = (1,1)T 。
- 试用牛顿法求解：
min f (x) = -
1
1 2
取初始点 x(0) = (4,0)T ，并将采用变步长和采用固定步长l = 1.0 时的情形做比较。
- 某工厂向用户提供发动机，按合同规定，其交货数量和日期是：第一季度末交40 台，第二季末交 60 台，第三季末交 80 台。工厂的最大生产能力为每季 100 台，每季的生产费用是 f ( x) = 50x + 0.2x2 （元），此处 x 为该季生产发动机的台数。若工厂生产的多，多余的发动机可移到下季向用户交货，这样，工厂就需支付存贮费，每台发动机每季的存贮费为 4 元。问该厂每季应生产多少台发动机，才能既满足交货合同，又使工厂所花费的费用最少（假定第一季度开始时发动机无存货）。
- 用 Matlab 的非线性规划命令 fmincon 求解飞行管理问题的模型二。
- 用罚函数法求解飞行管理问题的模型二。
图 3 影院剖面示意图
图 3 为影院的剖面示意图，座位的满意程度主要取决于视角 和仰角 b 。视角 是观众眼睛到屏幕上、下边缘视线的夹角， 越大越好；仰角 b 是观众眼睛到屏幕上边缘视线与水平线的夹角， b 太大使人的头部过分上仰，引起不舒服感，
!!!!!!!!!!!!!
点击链接加入群聊
https://qm.qq.com/q/NGl6WD0Bky
一般要求 b 不超过
30o 。
记影院屏幕高h ，上边缘距地面高 H ，地板线倾角q ，第一排和最后一排座位与屏幕水平距离分别为 d 和 D ，观众平均座高为c（指眼睛到地面的距离）。已知参数 h = 1.8 ， H = 5 ， d = 4.5 ， D = 19 ， c = 1.1 （单位：m）。
- 地板线倾角q = 10o ，问最佳座位在什么地方？
- 求地板线倾角（一般不超过20o ），使所有观众的平均满意程度最大。
- 地板线设计成什么形状可以进一步提高观众的满意程度。