第三章 非线性规划
§ 1 非线性规划
1.1 非线性规划的实例与定义
如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且， 也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。
下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。
例 1 （投资决策问题）某企业有n 个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金 A 元，投资于第i(i = 1, … , n) 个项目需花资金ai 元，并预计可收益bi 元。试选择最佳投资方案。
解 设投资决策变量为
xi 目，i = 1, … , n ，
则投资总额为 ai xi ，投资总收益为 bi xi 。因为该公司至少要对一个项目投资，并且总的投资金额不能超过总资金 A ，故有限制条件
另外，由于xi (i = 1, … , n) 只取值 0 或 1，所以还有
xi (1− xi ) = 0, i = 1, … , n.
最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取 0 或 1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：
s.t. ai xi ≤ A
xi (1− xi ) = 0, i = 1, … , n.
上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问
题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式min f (x)
s.t. hj (x) ≤ 0, j = 1, … , q (NP)
gi (x) = 0, i = 1, … , p
-32-
其中x = [x1 … xn ]T 称为模型（NP）的决策变量，f 称为目标函数，gi (i = 1, … , p)和 hj (j = 1, … , q) 称为约束函数。另外， gi (x) = 0 (i = 1, … , p) 称为等式约束， hj (x) ≤ 0 (j = 1, … , q) 称为不等式的约束。
对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：
（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。
（ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。
（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。
（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。
1.2 线性规划与非线性规划的区别
如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。
1.3 非线性规划的 Matlab 解法
Matlab 中非线性规划的数学模型写成以下形式
min f (x)
其中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 求下列非线性规划
min f (x) = x12 + x22 + x + 8
s.t. x − x2 + x ≥ 0
x1 + x + x ≤ 20
− x1 − x22 + 2 = 0
-33-
x2 + 2x = 3
x1 , x2 , x3 ≥ 0
解 （i）编写 M 文件 fun1.m 定义目标函数function f=fun1(x);
f=sum(x.^2)+8;
（ii）编写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)+2x(3)^2-3]; %非线性等式约束
（iii）编写主程序文件 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。
1.4 求解非线性规划的基本迭代格式
记（NP）的可行域为K 。
若x * ∈ K ，并且
f (x * ) ≤ f (x), ∀x ∈ K
则称x是（NP）的整体最优解，f (x * ) 是(NP)的整体最优值。如果有f (x * ) < f (x), ∀x ∈ K, x ≠ x *
则称x* 是（NP）的严格整体最优解，f (x * ) 是(NP)的严格整体最优值。
若x * ∈ K ，并且存在 x* 的邻域Nδ (x * ) ，使f (x* ) ≤ f (x), ∀x ∈ Nδ (x* ) ∩ K ，
则称x* 是（NP）的局部最优解，f (x * ) 是(NP)的局部最优值。如果有f (x* ) < f (x), ∀x ∈ Nδ (x* ) ∩ K
则称x* 是（NP）的严格局部最优解，f (x * ) 是(NP)的严格局部最优值。
由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然， 有时求出的某个解虽是一部分可行域上的极值点，但却并不一定是整个可行域上的全局最优解。
对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思想是：从一个选定的初始点 x0 ∈ R n 出发，按照某一特定的迭代规则产生一个点列{xk } ，使得当{xk } 是有穷点列时，其最后一个点是(NP)的最优解；当 {xk } 是无穷点列时，它有极限点，并且其极限点是(NP)的最优解。
设x k ∈ Rn 是某迭代方法的第k 轮迭代点，x k +1 ∈ Rn 是第k +1 轮迭代点，记
x k +1 = x k + tk pk （1）
这里tk ∈ R1 , pk ∈ R n , pk = 1 ，并且pk 的方向是从点x k 向着点xk +1 的方向。式（1）就是求解非线性规划模型(NP)的基本迭代格式。
-34-
通常，我们把基本迭代格式（1）中的pk 称为第k 轮搜索方向，tk 为沿pk 方向的步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步长。
设x ∈ R n , p ≠ 0 ，若存在δ > 0 ，使
f (x + tp) < f (x ), ∀t ∈(0, δ) ，称向量p 是f 在点x 处的下降方向。
设x ∈ R n , 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 在点 x k 处关于K 的可行下降方向作为搜索方向pk 。
2° 寻求搜索步长。以xk 为起点沿搜索方向pk 寻求适当的步长tk ，使目标函数值有某种意义的下降。
3° 求出下一个迭代点。按迭代格式（1）求出
x k +1 = x k + tk pk 。
若x k +1 已满足某种终止条件，停止迭代。
4° 以xk +1 代替x k ，回到 1°步。
1.5 凸函数、凸规划
设 f(x) 为定义在 n 维欧氏空间 E(n ) 中某个凸集 R 上的函数，若对任何实数α(0 < α < 1) 以及R 中的任意两点x( 1) 和x(2) ，恒有
f ( αx( 1) + (1 − α)x( 2) ) ≤ αf (x( 1) ) + (1 − α)f (x( 2) )则称f(x) 为定义在R 上的凸函数。
若对每一个 α(0 < α < 1) 和 x( 1) ≠ x( 2) ∈ R 恒有
f ( αx( 1) + (1 − α)x( 2) ) < αf (x( 1) ) + (1 − α)f (x( 2) )则称f(x) 为定义在R 上的严格凸函数。
考虑非线性规划
假定其中f (x) 为凸函数，gj (x)(j = 1,2, … , l) 为凸函数，这样的非线性规划称为凸规划。
可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数f(x) 为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见， 凸规划是一类比较简单而又具有重要理论意义的非
-35-
线性规划。
§2 无约束问题
2.1 一维搜索方法
当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法， 0.618 法等）；（2）插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线法，二分法等）。
考虑一维极小化问题
min 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, t1 ] ，因 而 [a, t1 ] 是缩短了 的单峰区间 ；若 f (t1 ) < f (t2 ) ，则有t* ∈ [t2, b] ，故 [t2, b] 是缩短了的单峰区间；若f (t2 ) = f (t1 ) ，则 [a, t1 ] 和 [t2, b] 都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较， 总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行，在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（2）最优解的近似值。
应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？
2.1.1 Fibonacci 法
如用Fn 表示计算n 个函数值能缩短为单位长区间的最大原区间长度，可推出Fn 满足关系
F0 = F1 = 1
Fn = Fn −2 +Fn −1 , n = 2,3, … ,
数列{Fn } 称为 Fibonacci 数列，Fn 称为第n 个 Fibonacci 数，相邻两个 Fibonacci 数之比 称为 Fibonacci 分数。
当用斐波那契法以 n 个探索点来缩短某一区间时，区间长度的第一次缩短率为 , 其后各次分别为 。由此， 若t1 和t2 (t2 < t1 ) 是单峰区间 [a, b]中第 1 个和第 2 个探索点的话，那么应有比例关系
从而
它们关于 [a, b] 确是对称的点。
-36-
如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度δ> 0 ，这就要求最后区间的长度不超过δ ,即
（4）
据此，我们应按照预先给定的精度 δ ,确定使（4）成立的最小整数n 作为搜索次数，直到进行到第n 个探索点时停止。
用上述不断缩短函数f(t) 的单峰区间 [a, b] 的办法，来求得问题（2）的近似解，是 Kiefer(1953 年)提出的，叫做 Finbonacci 法，具体步骤如下：
1o 选取初始数据，确定单峰区间 [a0, b0 ] ，给出搜索精度 δ > 0 ，由（4）确定搜索次数n 。
2o k = 1, a = a0, b = b0 ，计算最初两个搜索点，按（3）计算t1 和 t2 。
3o while k < n −1
f 1 = f (t1 ), f2 = f (t2 )
if f1 < f2
else
end
k = k +1
end
4o 当进行至k = n −1时，
这就无法借比较函数值f (t1 ) 和f(t2 ) 的大小确定最终区间，为此，取
其中ε 为任意小的数。在t1 和t2 这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小值。并得最终区间 [a, t1 ] 或 [t2, b] 。
由上述分析可知，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，它能以尽量少的函数求值次数，达到预定的某一缩短率。
例 3 试用斐波那契法求函数f(t) = t2 − t + 2 的近似极小点，要求缩短后的区间不大于区间[−1,3] 的 0.08 倍。
程序留作习题。
2.1.2 0.618 法
若 ω > 0 ，满足比例关系
-37-
称之为黄金分割数，其值为 黄金分割数 ω 和 Fibonacci 分数之间有着重要的关系
现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（0.618 法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果也相当好，因而易于为人们所接受。
用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短 0.618 倍。计算n 个探索点的函数值可以把原区间 [a0, b0 ] 连续缩短n −1次，因为每次的缩短率均为 μ ,故最后的区间长度为
(b0 − a0 ) μ n −1
这就是说，当已知缩短的相对精度为 δ 时，可用下式计算探索点个数n ：
μ n −1 ≤ δ
当然，也可以不预先计算探索点的数目n ，而在计算过程中逐次加以判断，看是否已满足了提出的精度要求。
0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618倍和 0.382 倍处。
2.2 二次插值法
对极小化问题（2），当 f(t) 在 [a, b] 上连续时，可以考虑用多项式插值来进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。
2.3 无约束极值问题的解法
无约束极值问题可表述为
min f (x), x ∈ E(n ) （5）
求解问题（5）的迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数，称为解析法。另一是仅用到函数值，称为直接法。
2.3.1 解析法
2.3.1.1 梯度法（最速下降法）
对基本迭代格式
x k +1 = x k + tk pk （6）
我们总是考虑从点xk 出发沿哪一个方向pk ，使目标函数f 下降得最快。微积分的知识告诉我们，点xk 的负梯度方向
pk = −▽f (xk ) ，
是从点xk 出发使f 下降最快的方向。为此，称负梯度方向 − ▽f(xk ) 为f 在点xk 处的最速下降方向。
-38-
按基本迭代格式（6），每一轮从点xk 出发沿最速下降方向 − ▽f(xk ) 作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。
这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，用 ▽f (xk ) = 0 或作为停止条件。其具体步骤如下：
1°选取初始数据。选取初始点x0 ，给定终止误差，令 k := 0 。
2°求梯度向量。计算 ▽f (xk ) ， 若 ▽f (xk ) ≤ ε ,停止迭代，输出 x k 。否则，进行 3°。
3° 构造负梯度方向。取
pk = −▽f (xk ) .
4° 进行一维搜索。求tk ，使得
令xk +1 = xk + tk pk , k := k +1, 转 2°。
例 4 用最速下降法求解无约束非线性规划问题
min f (x) = x + 25x
其中x = (x1, x2 )T ，要求选取初始点 x0 = (2,2)T 。解 （i） ▽f (x) = (2x1,50x2 )T
编写 M 文件 detaf.m，定义函数f (x) 及其梯度列向量如下function [f,df]=detaf(x);
f=x(1)2+25*x(2)2;
df=[2x(1)
50x(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+tp);
while f>f0
t=t/2;
f=detaf(x+tp);
end
x=x+t*p;
[f0,g]=detaf(x);
end
x,f0
2.3.1.2 Newton 法
考虑目标函数f 在点xk 处的二次逼近式
假定 Hesse 阵
-39-
正定。
由于 ▽2 f (xk ) 正定，函数Q 的驻点xk +1 是Q(x) 的极小点。为求此极小点，令
▽Q(xk +1) = ▽f (xk ) + ▽2 f (xk )(xk +1 − xk ) = 0 ，即可解得
x k +1 = xk − [▽2 f (xk )]−1 ▽f (xk ) .
对照基本迭代格式（1），可知从点 xk 出发沿搜索方向。
pk = −[▽2 f (xk )]−1 ▽f (xk )
并取步长 tk = 1 即可得 Q(x) 的最小点 xk +1 。通常， 把方向 pk 叫做从点 xk 出发的Newton 方向。从一初始点开始，每一轮从当前迭代点出发，沿 Newton 方向并取步长为 1 的求解方法，称之为 Newton 法。其具体步骤如下：
1°选取初始数据。选取初始点x0 ，给定终止误差 ε > 0 ，令 k := 0 。
2°求梯度向量。计算 ▽f (xk ) ，若▽f (xk ) ≤ ε ,停止迭代，输出 xk 。否则，进行 3°。
3°构造 Newton 方向。计算 [▽2 f (xk )]−1 ，取
pk = −[▽2 f (xk )]−1 ▽f (xk ) .
4° 求下一迭代点。令x k +1 = x k +pk , k := k +1，转 2°。
例 5 用 Newton 法求解，
min f (x) = x + 25x + xx选取x0 = (2,2)T 。
解：（i） ▽f (x) = [4x + 2x1x 100x + 2xx2 ]T
（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=[4x(1)3+2*x(1)*x(2)2;100x(2)3+2*x(1)2x(2)];
d2f=[2x(1)2+2*x(2)2,4*x(1)x(2)
4x(1)x(2),300x(2)2+2*x(1)2];
（III）编写主程序文件 example5.m 如下：
clc
x=[2;2];
[f0,g1,g2]=nwfun(x);
while norm(g1)>0.00001
p=-inv(g2)*g1;
x=x+p;
-40-
[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§;
t=1.0;f=nwfun(x+tp);
while f>f0
t=t/2;f=nwfun(x+tp);
end
x=x+tp;
[f0,g1,g2]=nwfun(x);
end
x,f0
Newton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算 − [▽2 f (xk )]−1 的工作量很大。
2.3.1.3 变尺度法
变尺度法（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）。
下面研究如何构造这样的近似矩阵，并将它记为 。我们要求：每一步都能以现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似矩阵最后应收敛于解点处的 Hesse 阵的逆阵。
当f (x) 是二次函数时，其 Hesse 阵为常数阵 A ，任两点xk 和xk +1 处的梯度之差为
▽f (xk +1) − ▽f (xk ) = A(xk +1 − xk )或
对于非二次函数，仿照二次函数的情形，要求其 Hesse 阵的逆阵的第k +1 次近似矩阵H (k +1) 满足关系式
-41-
这就是常说的拟 Newton 条件。
若令
则式（7）变为
现假定H (k ) 已知，用下式求H (k +1) （设 H (k ) 和H (k +1) 均为对称正定阵）；
其中 称为第k 次校正矩阵。显然，H (k +1) 应满足拟 Newton 条件（9），即要求
或
由此可以设想， ΔH (k ) 的一种比较简单的形式是
ΔH (k ) = Δxk (Q(k ) )T - H (k ) ΔG(k ) (W (k ) )T (12)其中Q(k ) 和W(k ) 为两个待定列向量。
将式（12）中的 ΔH (k ) 代入（11），得
这说明，应使
(Q(k ) )T ΔG(k ) = (W (k ) )T ΔG(k ) = 1 （13）
考虑到 ΔH (k ) 应为对称阵，最简单的办法就是取
由式（13）得
若(Δxk )T ΔG(k ) 和 不等于零，则有
于是，得校正矩阵
从而得到
上述矩阵称为尺度矩阵。通常，我们取第一个尺度矩阵H(0) 为单位阵，以后的尺度矩
-42-
阵按式（18）逐步形成。可以证明：
（i）当 xk 不是极小点且H (k ) 正定时，式（17）右端两项的分母不为零，从而可按式（18）产生下一个尺度矩阵
（ii）若H (k ) 为对称正定阵，则由式（18）产生的 H (k +1) 也是对称正定阵；
（iii）由此推出 DFP 法的搜索方向为下降方向。
现将 DFP 变尺度法的计算步骤总结如下。
1°给定初始点x0 及梯度允许误差 ε > 0 。
2°若 ，则 x0 即为近似极小点，停止迭代，否则，转向下一步。 3°令
H ( 0) = I （单位矩阵），
p0 = −H ( 0) ▽f (x0 )
在p0 方向进行一维搜索，确定最佳步长 λ0 ：
如此可得下一个近似点
x1 = x0 + λ0 p0
4 °一般地，设已得到近似点xk ，算出 ▽f (xk ) ，若
▽f (x0 ) ≤ ε
则x k 即为所求的近似解，停止迭代；否则，计算H (k ) ：
并令pk = −H(k ) ▽f (xk ) ，在 pk 方向上进行一维搜索，得 λk ，从而可得下一个近似点x k +1 = xk + λk pk
5°若xk +1 满足精度要求，则xk +1 即为所求的近似解，否则，转回 4°,直到求出某点满足精度要求为止。
2.3.2 直接法
在无约束非线性规划方法中，遇到问题的目标函数不可导或导函数的解析式难以表示时，人们一般需要使用直接搜索方法。同时， 由于这些方法一般都比较直观和易于理解，因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。
这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成， 具体步骤如下：
1° 选取初始数据。选取初始点x0 ，n 个线性无关初始方向，组成初搜索方向组{p0 , p1 , … , pn −1 } 。给定终止误差 ε > 0 ，令 k := 0 。
2°进行基本搜索。令 y0 := xk ，依次沿 {p0 , p1 , … , pn −1 } 中的方向进行一维搜索。对应地得到辅助迭代点y1 , y2 , … , y n ，即
yj = yj −1 + tj −1pj−1
f j = 1, … , n
-43-
3°构造加速方向。令p n = y n − y 0 ，若 p n ≤ ε ,停止迭代，输出 x k +1 = y n 。否则进行 4°。
4°确定调整方向。按下式
f (y m−1) − f (y m ) = max{f (yj −1) − f (yj ) | 1 ≤ j ≤ n}找出m 。若
f (y 0 ) − 2f (y n) + f (2y n − y 0 ) < 2[f (y m−1) − f (y m )]成立，进行 5°。否则，进行 6°。
5°调整搜索方向组。令
x k +1 = y n + tn p n : f (y n + tn p n) = min f (y n + tpn ) .
t ≥0
同时，令
{p0 , p1 , … , pn +1 }k +1 := {p0 , … , pm−1 , pm+1 , … , pn −1 , pn} ， k := k +1，转 2°。
6°不调整搜索方向组。令x k +1 := y n , k := k +1，转 2°。
2.4 Matlab 求无约束极值问题
在 Matlab 工具箱中，用于求解无约束极值问题的函数有 fminunc 和 fminsearch，用法介绍如下。
求函数的极小值
其中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(x2 − x12 )2 + (1 − x1 )2 的最小值。
解：编写 M 文件 fun2.m如下：
function [f,g]=fun2(x);
f=100*(x(2)-x(1)2)2+(1-x(1))^2;
g=[-400x(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);
-44-
编写主函数文件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 约束极值问题
带有约束条件的极值问题称为约束极值问题，也叫规划问题。
求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作， 可采用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题变换为较简单问题的其它方法。
库恩—塔克条件是非线性规划领域中最重要的理论成果之一，是确定某点为最优点的必要条件，但一般说它并不是充分条件（对于凸规划，它既是最优点存在的必要条件，同时也是充分条件）。
3.1 二次规划
若某非线性规划的目标函数为自变量x 的二次函数，约束条件又全是线性的，就称这种规划为二次规划。
Matlab 中二次规划的数学模型可表述如下：
这里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 求解二次规划
解 编写如下程序：
-45-
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))
求得
3.2 罚函数法
利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题，因而也称这种方法为序列无约束最小化技术，简记为 SUMT (Sequential Unconstrained Minization Technique)。
罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。
考虑问题：
min f (x)
取一个充分大的数 M > 0 ，构造函数
这里 G(x) = [g1 (x) … gr (x)] ，H(x) = [h1 (x) … hs (x)]，
K(x) = [k1 (x) … kt (t)]，Matlab 中可以直接利用 max 、min 和 sum 函数。）则以增广目标函数P(x, M) 为目标函数的无约束极值问题
min P(x, M)的最优解x 也是原问题的最优解。
例 9 求下列非线性规划
解 (i)编写 M 文件 test.m
function g=test(x);
M=50000;
f=x(1)2+x(2)2+8;
g=f-Mmin(x(1),0)-Mmin(x(2),0)-M*min(x(1)^2-x(2),0)+…
-46-
Mabs(-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-Msum(min([x’;zeros(1,2)]))-Mmin(x(1)^2-x(2),0)+…
Mabs(-x(1)-x(2)^2+2);
我们也可以修改罚函数的定义，编写test.m如下： function g=test(x);
M=50000;
f=x(1)2+x(2)2+8;
g=f-Mmin(min(x),0)-Mmin(x(1)2-x(2),0)+M*(-x(1)-x(2)2+2)^2;
(ii)在 Matlab 命令窗口输入
[x,y]=fminunc(‘test’,rand(2,1))即可求得问题的解。
3.3 Matlab 求约束极值问题
在 Matlab 优化工具箱中，用于求解约束最优化问题的函数有：fminbnd 、fmincon、 quadprog 、fseminf、fminimax，上面我们已经介绍了函数 fmincon 和 quadprog。
3.3.1 fminbnd 函数
求单变量非线性函数在区间上的极小值
Matlab 的命令为
[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)即可求得极小点和极小值。
3.3.2 fseminf 函数
求
min {F(x) | C(x) ≤ 0, Ceq(x) = 0, PHI(x, w) ≤ 0}
x
其中 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) ，非线性等
-47-
式约束 Ceq(x) 和半无穷约束PHI(x, w) 的每一个分量函数，函数 SEMINFCON 有两个输入参量 X 和 S ，S 是推荐的取样步长，也许不被使用。
例 11 求函数 f (x) = (x1 - 0.5)2 + (x2 - 0.5)2 + (x3 - 0.5)2 取最小值时的 x 值,约束为：
1 ≤ w1 ≤ 100 ，1 ≤ w2 ≤ 100
解 （1）编写 M 文件 fun6.m 定义目标函数如下：
function f=fun6(x,s);
f=sum((x-0.5).^2);
（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(w1x(1)).cos(w1x(2))-1/1000(w1-50).^2-sin(w1x(3))-x(3)-1;
k2=sin(w2x(2)).cos(w2x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1;
%画出半无穷约束的图形
plot(w1,k1,‘-’,w2,k2,‘+’);
（3）调用函数 fseminf
在 Matlab 的命令窗口输入
[x,y]=fseminf(@fun6,rand(3,1),2,@fun7)即可。
3.3.3 fminimax 函数求解 x F
其中F(x) = {F1 (x),… , Fm (x)}。
上述问题的 Matlab 命令为
X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)
-48-
例 12 求函数族{f1 (x), f2 (x), f3 (x), f4 (x), f5 (x)} 取极大极小值时的x 值。其中：
解 （1）编写 M 文件 fun8.m 定义向量函数如下：
function f=fun8(x);
（2）调用函数 fminimax
[x,y]=fminimax(@fun8,rand(2,1))
3.3.4 利用梯度求解约束优化问题
例 13 已知函数f (x) = ex1 (4x + 2x + 4x1x2 + 2x2 + 1) ，且满足非线性约束：
求min f (x)
x
分析：当使用梯度求解上述问题时，效率更高并且结果更准确。
题目中目标函数的梯度为：
解 （1）编写 M 文件 fun9.m 定义目标函数及梯度函数：
function [f,df]=fun9(x);
f=exp(x(1))(4x(1)2+2*x(2)2+4x(1)x(2)+2x(2)+1);
df=[exp(x(1))(4x(1)2+2*x(2)2+4x(1)x(2)+8x(1)+6x(2)+1);exp(x(1))(4x(2) +4x(1)+2)];
（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=[];
（3）调用函数 fmincon，编写主函数文件 example13.m 如下：
%采用标准算法
options=optimset(‘largescale’,‘off’);
%采用梯度
options=optimset(options,‘GradObj’,‘on’,‘GradConstr’,‘on’);
[x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options)
-49-
3.4 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 的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞， 则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下：
1）不碰撞的标准为任意两架飞机的距离大于 8km;
2）飞机飞行方向角调整的幅度不应超过 30 度；
3）所有飞机飞行速度均为每小时 800km;
4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在 60km 以上；
5）最多需考虑 6 架飞机；
6）不必考虑飞机离开此区域后的状况。
请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过 0.01 度），要求飞机飞行方向角调整的幅度尽量小。
设该区域 4 个顶点的座标为(0,0) ，(160,0) ，(160,160) ，(0,160)。记录数据见表 1。
-50-
表 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；
(x, y) 为第i 架飞机的初始位置；
(xi (t), yi (t)) 为第i 架飞机在t 时刻的位置；
θi0 为第i 架飞机的原飞行方向角，即飞行方向与x 轴夹角，0 ≤ θi0 < 2π ;
Δ θi 为第i 架飞机的方向角调整，
θi = θi0 + Δθi 为第i 架飞机调整后的飞行方向角。
4.1 模型一
根据相对运动的观点在考察两架飞机i 和 j 的飞行时，可以将飞机i 视为不动而飞机j 以相对速度
vij = vj - vi = (a cos θj - a cos θi, a sinθj - a sinθi ) （19）相对于飞机i 运动，对（19）式进行适当的化约可得
不妨设θj ≥ θi ，此时相对飞行方向角为 ，见图 2。
图 2 相对飞行方向角由于两架飞机的初始距离为
-51-
则只要当相对飞行方向角 βij 满足
α < βij < 2π − α （23）
时，两架飞机不可能碰撞（见图 2）。
记 β为调整前第 j 架飞机相对于第i 架飞机的相对速度（矢量）与这两架飞机连
线（从j 指向i 的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负）。则由式（23）知，两架飞机不碰撞的条件为
其中 ＝相对速度vmn 的幅角－从n 指向m 的连线矢量的幅角
（注意 βn 表达式中的i 表示虚数单位）这里我们利用复数的幅角，可以很方便地计算角度 βn （ m, n = 1,2, … ,6 ）。
本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小；所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数，可以得到如下的数学规划模型：
Δ θi ≤ 30o ， i = 1,2, … ,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+iy0
xy2=exp(iq*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
-52-
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)
求得 α 的值如表 2 所示。
表 2 α 的值
1 2 3 4 5 61 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
求得 β 的值如表 3 所示。
表 3 β 的值
1 2 3 4 5 61 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.5delta(i)+0.5delta(j))>a
lpha(i,j));
end
求得的最优解为 Δθ3 = 2.83858o ， Δ θ5 = -21.0138o ， Δ θ6 = 0.7908o ，其它调整角度为 0。
4.2 模型二
-53-
两架飞机i, j 不发生碰撞的条件为
(xi (t) − xj (t))2 + (yi (t) − yj (t))2 > 64 （25）
1 ≤ i ≤ 5 ，i +1 ≤ j ≤ 6 ，0 ≤ t ≤ min{Ti, Tj }
其中 Ti, Tj 分别表示第i, j 架飞机飞出正方形区域边界的时刻。这里
xi (t) = x + at cos θi ，yi (t) = y + at sinθi ，i = 1,2, … , n ；
θi = θi0 + Δθi ， | Δ θi i = 1,2, … , n ；
下面我们把约束条件（25）加强为对所有的时间t 都成立，记
li ,j = (xi (t) − xj (t))2 + (yi (t) − yj (t))2 − 64 = (i, j)t 2 + (i, j)t + (i, j)
其中 i, j) = 4a2 sin
(i, j) = 2a[(xi (0) − xj (0))(cosθi − cos θj ) + (yi (0) − yj (0))(sinθi − sinθj )](i, j) = (xi (0) − xj (0))2 + (yi (0) − yj (0))2 − 64
则两架i, j 飞机不碰撞的条件是
Δ (i, j) = (i, j)2 − 4(i, j)(i, j) < 0 （26）
这样我们建立如下的非线性规划模型
s.t. Δ (i, j) < 0 ，1 ≤ i ≤ 5 ， i + 1 ≤ j ≤ 6
i = 1,2, … ,6
习 题 三
- 用最速下降法（梯度法）求函数：
f (x) = 4x1 + 6x2 − 2x − 2x1x2 − 2x的极大点。给定初始点x 0 = (1,1)T 。 - 试用牛顿法求解：
取初始点x (0) = (4,0)T ，并将采用变步长和采用固定步长 λ = 1.0 时的情形做比较。
3. 某工厂向用户提供发动机，按合同规定，其交货数量和日期是：第一季度末交40 台，第二季末交 60 台，第三季末交 80 台。工厂的最大生产能力为每季 100 台，每季的生产费用是f (x) = 50x + 0.2x2 （元），此处x 为该季生产发动机的台数。若工厂生产的多，多余的发动机可移到下季向用户交货，这样，工厂就需支付存贮费，每台发动机每季的存贮费为 4 元。问该厂每季应生产多少台发动机， 才能既满足交货合同，又使工厂所花费的费用最少（假定第一季度开始时发动机无存货）。
4. 用 Matlab 的非线性规划命令 fmincon 求解飞行管理问题的模型二。
5. 用罚函数法求解飞行管理问题的模型二。
-54-
- 求下列问题的解
max f (x) = 2x1 + 3x + 3x2 + x + x3
s.t. x1 + 2x + x2 + 2x + x3 ≤ 10
x1 + x + x2 + x − x3 ≤ 50 2x1 + x + 2x2 + x3 ≤ 40 x + x3 = 2
x1 + 2x2 ≥ 1
x1 ≥ 0, x2, x3不约束7．
图 3 影院剖面示意图
图 3 为影院的剖面示意图，座位的满意程度主要取决于视角 α 和仰角 β 。视角 α 是观众眼睛到屏幕上、下边缘视线的夹角， α 越大越好；仰角 β 是观众眼睛到屏幕上边缘视线与水平线的夹角， β 太大使人的头部过分上仰，引起不舒服感，一般要求 β 不超过30o 。
记影院屏幕高h ，上边缘距地面高H ，地板线倾角θ , 第一排和最后一排座位与屏幕水平距离分别为d 和D ，观众平均座高为c（指眼睛到地面的距离）。已知参数h = 1.8 ， H = 5 ，d = 4.5 ，D = 19 ，c = 1.1 （单位：m）。
（1）地板线倾角 θ = 10o ，问最佳座位在什么地方？
（2）求地板线倾角（一般不超过20o ），使所有观众的平均满意程度最大。
（3）地板线设计成什么形状可以进一步提高观众的满意程度。
)
)