【数学建模】数学建模学习3---非线性规划（例题+matlab代码实现）

1 非线性规划

1.1 非线性规划的实例与定义

如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。

下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。

例 1 （投资决策问题）某企业有 n 个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金 A 元，投资于第i(i = 1,L,n) 个项目需花资金ai 元，并预计可收益bi 元。试选择最佳投资方案。

解 设投资决策变量为

最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取 0 或 1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：

上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式

对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：

（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。

（ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。

（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。

（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。

1.2 线性规划与非线性规划的区别

如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。

1.3 非线性规划的 Matlab 解法

Matlab 中非线性规划的数学模型写成以下形式

其中 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 求下列非线性规划

解 （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)+2*x(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 凸函数、凸规划

设 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)为凸函数，g j(x)( j =1,2,L,l) 为凸函数，这样的非线性规划称为凸规划。

可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数 f (x)为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非线性规划。

2 无约束问题

2.1 一维搜索方法

当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法，0.618 法等）；（2）插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线法，二分法等）。

考虑一维极小化问题

若 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 满足关系

数列{Fn }称为 Fibonacci 数列， Fn 称为第n 个 Fibonacci 数，相邻两个 Fibonacci 数之比 称为 Fibonacci 分数。

如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度δ > 0 ，这就要求最后区间的长度不超过δ ，即

据此，我们应按照预先给定的精度δ ，确定使（4）成立的最小整数n 作为搜索次数，直到进行到第n 个探索点时停止。

用上述不断缩短函数 f (t) 的单峰区间[a,b] 的办法，来求得问题（2）的近似解，是 Kiefer(1953 年)提出的，叫做 Finbonacci 法，具体步骤如下：

这就无法借比较函数值 f (t1 )和 f (t2 ) 的大小确定最终区间，为此，取

其中ε 为任意小的数。在t1 和t2 这两点中，以函数值较小者为近似极小点，相应的函数值为近似极小值。并得最终区间[a,t1 ]或[t2 ,b] 。

由上述分析可知，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，它能以尽量少的函数求值次数，达到预定的某一缩短率。

例 3 试用斐波那契法求函数 f (t) = t 2 − t + 2 的近似极小点，要求缩短后的区间不大于区间[−1,3]的 0.08 倍。

程序留作习题。

2.1.2 0.618 法

若ω > 0 ，满足比例关系

称之为黄金分割数，其值为 0.6180339887L

黄金分割数ω 和 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 梯度法（最速下降法）

对基本迭代格式

我们总是考虑从点 xk 出发沿哪一个方向 pk ，使目标函数 f 下降得最快。微积分的知识告诉我们，点 xk 的负梯度方向pk = −∇f (xk ) ，是从点 xk 出发使 f 下降最快的方向。为此，称负梯度方向 − ∇f (xk ) 为 f 在点 xk 处的最速下降方向。

按基本迭代格式（6），每一轮从点 xk 出发沿最速下降方向 − ∇f (xk ) 作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。

这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，用∇f (xk ) = 0或 ∇f (xk ) ≤ ε 作为停止条件。其具体步骤如下：

1°选取初始数据。选取初始点 x 0 ，给定终止误差，令k := 0 。

2°求梯度向量。计算∇f (xk ) ， 若 ∇f (xk ) ≤ ε ，停止迭代，输出 xk 。否则，进行 3°。

3° 构造负梯度方向。取pk = −∇f (xk ).

4° 进行一维搜索。求tk ，使得

例 4 用最速下降法求解无约束非线性规划问题

解 （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=[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

2.3.1.2 Newton 法

考虑目标函数 f 在点 xk 处的二次逼近式

假定 Hesse 阵

由于∇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 )

并取步长 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° 求下一迭代点。令 xk +1 = xk + pk ,k := k +1，转 2°。

例 5 用 Newton 法求解，

（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 如下：

clc
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 的工作量很大。

2.3.2 直接法

在无约束非线性规划方法中，遇到问题的目标函数不可导或导函数的解析式难以表示时，人们一般需要使用直接搜索方法。同时，由于这些方法一般都比较直观和易于理解，因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。

这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成，具体步骤

如下：

1° 选取初始数据。选取初始点 x0 ， n 个线性无关初始方向，组成初搜索方向组{p0 , p1,L, pn−1}。给定终止误差ε > 0，令k := 0 。

2°进行基本搜索。令 y 0 := xk ，依次沿{p0 , p1 ,L, pn−1}中的方向进行一维搜索。对应地得到辅助迭代点 y1 , y 2 ,L, y n ，即

3°构造加速方向。令 pn = y n − y 0 ，若 pn ≤ ε ，停止迭代，输出 xk +1 = y n 。否则进行 4°。

4°确定调整方向。按下式

找出m 。若

成立，进行 5°。否则，进行 6°。

5°调整搜索方向组。令

同时，令

k := k +1，转 2°。

6°不调整搜索方向组。令 xk +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=[-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 约束极值问题

带有约束条件的极值问题称为约束极值问题，也叫规划问题。

求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作，可采用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题变换为较简单问题的其它方法。

库恩—塔克条件是非线性规划领域中最重要的理论成果之一，是确定某点为最优点的必要条件，但一般说它并不是充分条件（对于凸规划，它既是最优点存在的必要条件，同时也是充分条件）。

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 求解二次规划

解 编写如下程序：

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)。

罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。

例 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)+...-47- 
 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)) 

即可求得问题的解。