第三章 非线性规划
§1 非线性规划
1.1 非线性规划的实例与定义
如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。
下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。
例1 （投资决策问题）某企业有个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金元，投资于第个项目需花资金元，并预计可收益元。试选择最佳投资方案。
解 设投资决策变量为
，，
则投资总额为，投资总收益为。因为该公司至少要对一个项目投资，并且总的投资金额不能超过总资金，故有限制条件
另外，由于只取值0或1，所以还有
最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取0或1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：
s.t.
上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中目标函数或约束条件中至少有一个非线性函数，这类问题称之为非线性规划问题，简记为（NP）。可概括为一般形式
(NP)
其中称为模型（NP）的决策变量，称为目标函数，和称为约束函数。另外， 称为等式约束， 称为不等式约束。
对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：
（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。
（ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。
（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。
（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。
1.2 线性规划与非线性规划的区别
如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。
1.3 非线性规划的Matlab解法
Matlab中非线性规划的数学模型写成以下形式
,
其中是标量函数，是相应维数的矩阵和向量，是非线性向量函数。
Matlab中的命令是
X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)
它的返回值是向量，其中FUN是用M文件定义的函数；X0是的初始值；A,B,Aeq,Beq定义了线性约束，如果没有等式约束，则A=[],B=[],Aeq=[],Beq=[]；LB和UB是变量的下界和上界，如果上界和下界没有约束，则LB=[]，UB=[]，如果无下界，则LB=-inf，如果无上界，则UB=inf；NONLCON是用M文件定义的非线性向量函数；OPTIONS定义了优化参数，可以使用Matlab缺省的参数设置。
例2 求下列非线性规划问题
（i）编写M文件fun1.m
function f=fun1(x);
f=x(1)^2+x(2)^2+8;
和M文件fun2.m
function [g,h]=fun2(x);
g=-x(1)^2+x(2);
h=-x(1)-x(2)^2+2; %等式约束
（ii）在Matlab的命令窗口依次输入
options=optimset;
[x,y]=fmincon('fun1',rand(2,1),[],[],[],[],zeros(2,1),[], ...
'fun2', options)
就可以求得当时，最小值。
1.4 求解非线性规划的基本迭代格式
记（NP）的可行域为。
若，并且
则称是（NP）的整体最优解，是(NP)的整体最优值。如果有
则称是（NP）的严格整体最优解，是(NP)的严格整体最优值。
若，并且存在的邻域，使
，
则称是（NP）的局部最优解，是(NP)的局部最优值。如果有
则称是（NP）的严格局部最优解，是(NP)的严格局部最优值。
由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然，有时求出的某个解虽是一部分可行域上的极值点，但并不一定是整个可行域上的全局最优解。
对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思想是：从一个选定的初始点出发，按照某一特定的迭代规则产生一个点列，使得当是有穷点列时，其最后一个点是(NP)的最优解；当是无穷点列时，它有极限点，并且其极限点是(NP)的最优解。
设是某迭代方法的第轮迭代点，是第轮迭代点，记
（1）
这里，显然是由点与点确定的方向。式（1）就是求解非线性规划模型(NP)的基本迭代格式。
通常，我们把基本迭代格式（1）中的称为第轮搜索方向，为沿方向的步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步长。
设，若存在，使
，
称向量是在点处的下降方向。
设，若存在，使
，
称向量是点处关于的可行方向。
一个向量，若既是函数在点处的下降方向，又是该点关于区域的可行方向，则称之为函数在点处关于的可行下降方向。
现在，我们给出用基本迭代格式（1）求解(NP)的一般步骤如下：
0° 选取初始点，令。
1° 构造搜索方向，依照一定规则，构造在点处关于的可行下降方向作为搜索方向。
2° 寻求搜索步长。以为起点沿搜索方向寻求适当的步长，使目标函数值有某种意义的下降。
3° 求出下一个迭代点。按迭代格式（1）求出
。
若已满足某种终止条件，停止迭代。
4° 以代替，回到1°步。
1.5 凸函数、凸规划
设为定义在维欧氏空间中某个凸集上的函数，若对任何实数以及中的任意两点和，恒有
则称为定义在上的凸函数。
若对每一个和恒有
则称为定义在上的严格凸函数。
考虑非线性规划
假定其中为凸函数，为凸函数，这样的非线性规划称为凸规划。
可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非线性规划。
§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元。问该厂每季应生产多少台发动机，才能既满足交货合同，又使工厂所花费的费用最少（假定第一季度开始时发动机无存货）。
