【数学建模】数学建模学习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 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)) 

即可求得问题的解。