非线性规划模型

非线性规划模型

一、 简介

具有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。非线性规划研究一个n元实函数在一组等式或不等式的约束条件下的极值问题，且目标函数和约束条件至少有一个是未知量的非线性函数。目标函数和约束条件都是线性函数的情形则属于线性规划。

二、 背景

非线性规划是20世纪50年代才开始形成的一门新兴学科。1951年H.W.库恩和A.W.塔克发表的关于最优性条件（后来称为库恩-塔克条件）的论文是非线性规划正式诞生的一个重要标志。在50年代还得出了可分离规划和二次规划的n种解法，它们大都以G.B.丹齐克提出的解线性规划的单纯形法为基础的。50年代末到60年代末出现了许多解非线性规划问题的有效的算法，70年代又得到进一步发展。非线性规划在工程、管理、经济、科研、军事等方面都有广泛的应用，为最优设计提供了有利的工具。20世纪80年代以来，随着计算机技术的快速发展，非线性规划方法取得了长足进步，在信赖域法、稀疏拟牛顿法、并行计算、内点法和有限存储法等领域取得了丰硕的成功。

三、 非线性规划引例

问题1、容器设计问题

问题提出 某公司生产贮藏用容器，订货合同要求该公司制造一种敞口的长方形容器，容积为12立方米，该容器的底为正方形，容器总重量不超过68公斤。已知用作容器四壁的材料为每平方米10元，重3公斤；用作容器底的材料每平方米20元，重2公斤。试问制造该容器所需的最小费用是多少？

模型建立 设该容器的底边长和高分别为x1，x2。 则问题的数学模型为 minf

x()

4x02x1

2

2x0 1

x12x212212xx2x68121

s.t.

x,x012

问题2、营业计划的制定

问题提出 某公司经营两种设备，第一种设备每件售价30元，第二种设备每件售价450元。据统计，售出一件第一种设备所需要的营业时间平均是0.5小时，第二种设备是小时，其中x2是第二种设备的售出数量。已知该公（2+0.25x2）司在这段时间内的总营业时间为800小时。试决定使其营业额最大的营业计划。

模型建立 设该公司计划经营第一种设备x1件，第二种设备x2件。其数学模型为

maxf(x)30x1450x2

8000.5x1(20.25x2)x2

s.t. x,x0

12

四、 非线性规划的基本概念

minfx(

)

1,2,m,1,2,l（1,）.

gi(X)0,i

s.t. h(X)0,jj

Tn

其中X(x1,x2,,xn)E，f,gi,hj是定义在En上的实值函数，简

记：f:E

n

E,gi:EE,hj:EE

1n1n1

其他情况：求目标函数的最大值或约束条件为小于等于零的情况，都可通过取其相反数化为上述一般形式。

定义1 把满足问题（1）中条件的解x(E)称为可行解（或可行点），所有可行点的集合称为可行集（或可行域）。记为D。即

n

DXgi(X)



h0j,X(

)X0E,，问题（1）可简记为minf(X)。

XD

n



*定义2 对于问题（1），设XD，若存在0，使得对一切XD且

XX

*



，都有f(X*)f(X)，则称X*是f(X)在D上的局部极小值点（局

部最优解）。特别地当XX*时，若f(X*)f(X)，则称X*是f(X)在D上的严格局部极小值点（严格局部最优解）。

定义3 对于问题（1），设X*D，对任意的XD，都有f(X*)f(X)，则称X*是f(X)在D上的全局极小值点（全局最优解）。特别地当XX*，若则称Xf(X)f(X)，

*

*

是f(X)在D上的严格全局极小值点（严格全局最优解）。

五、 非线性规划的基本解法

1. 罚函数法

罚函数法基本思想是通过构造罚函数把约束问题转化为一系列无约束最优化问题，进而用无约束最优化方法去求解。这类方法称为序列无约束最小化方法。简称为SUMT法。

其一为SUMT外点法，其二为SUMT内点法。 1) SUMT外点法

对一般的非线性规划：

minfx( )

gi(X)0,i1,2,,m;

s.t. （1）

h(X)0,j1,2,,l.j

m

可设：T(X,M)f(X)M

min(0,gi(X))

i1

2

l

M



j1

hj(X) （2） 

2

将问题（1）转化为无约束问题：

T(Xmin

XE

n

（3） ,M)

其中T(X,M)称为罚函数，M称为罚因子，带M的项称为罚项，这里的罚函数只对不满足约束条件的点实行惩罚：当XD时，满足gi(X)0,hi(X)0，故罚项=0，不受惩罚。XD时，gi(X)0或hi(X)0的约束条件，故罚项>0，要受惩罚。

迭代步骤：

1、任意给定初始点X0，取M11，给定允许误差0，令k1；

k

2、求无约束极值问题minT(X,M)的最优解，设为XX(Mk)，即

n

XE

m

XE

inT(X,M)T(X

n

k

M,k；)

3、若存在i(1im)，使gi(Xk)，则取MkM(Mk1M,10)，

令kk1返回（2），否则，停止迭代。得最优解X*Xk。计算时也可将收敛性判别准则gi(X)改为Mmin(0,gi(X))0。

k

m

2

i1

2) SUNT内点法（障碍函数法）

minf(X)

考虑问题： （1）

s.t.g(X)0,i1,2,,mi

设集合D0Xgi(X)0,i1,2,,m，D0是可行域中所有严格内点的集合。

构造障碍函数

m

m

I(X,r)：I(X,r)f(X)r

m



i1

lngi(X)或I(X,r)f(X)r

i1

1gi(X)

。

m

其中称r

i1

lngi(X)或r

i1

1gi(X)

为障碍项，r为障碍因子。

k

这样问题（1）就转化为求一系列极值问题：minI(X,rk)得X(rk)。

XD

迭代步骤：

1、给定允许误差0，取r10，01； 2、求出约束集合D的一个内点X0D0，令k1；

3、以Xk1D0为初始点，求解minI(X,rk)，其中XD0的最优解，设为

XD

X

k

X(rk)D

；

m

k

m

4、检验是否满足rlngi(X)或rk

i1

1gi(X)

，满足，停止迭代，

i1

*k

XX；否则取rk1rk，令kk1，返回3。

3) 罚函数的缺点

每个近似最优解Xk往往不是容许解，而只能近似满足约束，在实际问题中这种结果可能不能使用；在解一系列无约束问题中，计算量太大，特别是随着Xk的增大，可能导致错误。

2. 近似规划法

基本思想：将问题中的目标函数f(X)和约束条件gi(X)0(i1,,m)；

hj(X)0(j1,,l)

近似为线性函数，并对变量的取值范围加以限制，从而得到

一个近似线性规划问题，再用单纯形法求解，把其符合原始条件的最优解作为问题的解的近似。每得到一个近似解后，都从这点出发，重复以上步骤。

这样，通过求解一系列线性规划问题，产生一个由线性规划最优解组成的序列，经验表明，这样的序列往往收敛于非线性规划问题的解。

迭代步骤：

11

1、给定初始可行点X1x11,x1,,xn，步长限制j(j1,,n)，步长缩小2

系数(0,1)，允许误差，令k1；

2、在点Xk处，将f(X)，gi(X)，hj(X)按泰勒级数展开并取一阶近似，得到近似线性规划问题：

kkT minfX()fX()fX()X(

k

X

i1,m,

kkT

)g(X)(X gi(X)gi(Xi

k

X)

kkT)h(X)(X hj(X)hj(Xj

k

X) j1,l,；

3、在上述近似线性规划问题的基础上增加一组限制步长的线性约束条件。

因为线性近似通常只在展开点附近近似程度较高，故需要对变量的取值范围加以限制，所增加的约束条件是：

xjxkjjk (j1,n, )求解该线性规划问题，得到最优解Xk1；

4、检验Xk1点对原约束是否可行。若Xk1对原约束可行，则转步骤5；否则，缩小步长限制，令jkjk (j13，重解当前的线性规,n,，返回步骤)划问题；

k15、判断精度：若jk (j1为近似最优解；否则，令,n,，则点)X

j

k1

j (j1,n,，)kk1，返回步骤2。

k

六、 用Matlab求解非线性规划问题

1. 二次规划

标准型为：

MinZ

12

XHXcX

T

T

AXb

Xbeq s.t. Aeq

VLBXVUB用Matlab求解，其输入格式如下：

1、x=quadprog(H,c,A,b)； 2、x=quadprog(H,c,A,b,Aeq,beq)；

3、x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)； 4、x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB,X0)； 5、x=quadprog(H,c,A,b,Aeq,beq,VLB,VUB,X0,options)； 6、[x,fval]=quadprog(...)； 7、[x,fval,exitflag]=quadprog(...)； 8、[x,fval,exitflag,output]=quadprog(......)；

例 minf(x1,x2)2x16x2x12x1x22x2

2

2

x1x22

x12x22

s.t.

x,x012

1

（1）写成标准形式：minz(x1,x2)

1

1x1

2x2

2x1

6x2

T

1

1

s.t. 

0

0

1x12



2x22x1x2

（2）输入命令：

H=[1 -1;-1 2]; c=[-2;-6]; A=[1 1;-1 2]; b=[2;2]; Aeq=[]; beq=[]; VLB=[0;0]; VUB=[];

[x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)

（3）运算结果为：

x=0.6667 1.3333 z=-8.2222 2. 一般非线性规划 标准型为：

minFX(

)

AXb

AeqXbeq

G(X)0

s.t. 

Ceq(X)0VLBXVUB

其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同。

用Matlab求解，基本步骤分三步：

（1）首先建立M文件fun.m，定义目标函数F(X)： function

ffun(X);

fF(X);

（2）若约束条件中有非线性约束：G(X)或Ceq(X)0，则建立M文件

nonlcon.m定义函数G(X)与Ceq(X)：

function [G,Cep]nonlcon(X)

G

Ceq

（3）建立主程序。非线性规划求解的函数是fmincon，命令的基本格式如下：

a) x=fmincon('fun',X0,A,b) b) x=fmincon('fun',X0,A,b,Aeq,beq)

c) x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB) d) x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB,'nonlcon') e) x=fmincon('fun',X0,A,b,Aeq,beq,VLB,VUB，'nonlcon',options) f) [x,fval]=fmincon()

g) [x,fval,exitflag]=fmincon()

h) [x,fval,exitflag,putput]=fmincon()

其中，x为输出极值点，fun为M文件，X0为迭代的初值，VLB、VUB为变量上下限，options为参数说明。

注意：

1) fmincon函数提供了大型优化算法和中型优化算法。默认时，若在fun函数中提供了梯度（options参数的GradObj设置为'on'），并且只有上下界存在或只有等式约束，fmincon函数将选择大型化。当既有等式约束又有梯度约束时，使用中型算法。

2) fmincon函数的中型算法使用的是序列二次规划法。在每一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日Hessian矩阵。 3) fmincon函数可能会给出局部最优解，这与初值X0的选取有关。 例1 minfx12x2

12

x

21

12

x2

2

2x13x26

x14x25

s.t.

x,x012

（1）写成标准形式：

minfx12x2

12

x

21

12

x2

2

2x13x260x14x250

s.t. 0x

1



0x2（2）建立M文件fun3.m：

function f=fun3(x);

f=-x(1)-2*x(2)+(1/2)*x(1)^x(1)^2+(1/2)*x(2)^2

（3）建立主程序youh2.m:

x0=[1;1]; A=[2 3;1 4]; b=[6;5]; Aeq=[]; beq=[]; VLB=[0;0]; VUB=[];

[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)

（4）运算结果为：

x=0.7647 1.0588 fval=-2.0294

1

例2 f(x)e(4x12x24x1x22x21)

x22

x1x20

1.5x1x2x1x20

s.t. 

xx10012

（1）建立M文件fun4.m，定义目标函数：

function f=fun4(x);

f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

（2）建立M文件mycon.m定义非线性约束：

function[g,cep]=mycon(x)

g=[x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];

（3）主程序youh3.m为：

x0=[-1;1]; A=[]; b=[]; Aeq=[1 1]; beq=[0]; VLB=[]; VUB=[];

[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,VLB,VUB,'mycon')

（4）运算结果为：

x=-1.2250 1.2250 fval=1.8951 七、 非线性规划应用实例

某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系a，b表示，

距离单位：千米）及水泥日用量d（吨）由下表给出。目前有两个临时料场位于A（5,1），B（2,7），日储量各有20吨。假设从料场到工地之间均有直线道路相连。

（1）试制定每天的供应计划，即从A，B两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。

（2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，

记工地的位置为（ai，bi），水泥日用量为di，i=1，···，6；料场位置为（xj，yj），日储量为ej，j=1，2；从料场j向工地i的运送量为Xij。

2

6

目标函数为：minf



j1i1

Xij2

Xijdi,i1,2,,6j

16

约束条件为：

Xijej,j1,2i1

当用临时料场时决策变量为：Xij， 当不用临时料场时决策变量为：Xij,xj,yj 2、使用临时料场的情形

使用两个临时料场A（5,1），B（2,7）。求从料场j向工地i的运送量为Xij，在各工地用量必须满足和各料场运送量不超过日储量的条件下，使总的吨千米数

最小，这是线性规划问题。线性规划模型为：

2

6

ij

minf

aa(i,j)X

j1i1

s.t.

2

i1,2,Xijd,i

j16j1,2Xije,ji1

,6

其中，aa(i,j)i1,2,,6.j1,2，为常数。

设X11X1，X21X2，X31X3，X41X4，X51X5，X61X6，X12X7，

X22X8，X32X9，X42X10，X52X11，X62X12。

编写matlab程序：

a=[1.25 8.75 0.5 5.75 3 7.25]; b=[1.25 0.75 4.75 5 6.5 7.25]; d=[3 5 4 7 6 11]; x=[5 2]; y=[1 7]; e=[20 20]; for i=1:6 for j=1:2

aa(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2); end end

CC=[aa(:1);aa(:2)]';

A=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1]; B=[20;20];

Aeq=[1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 ]

beq=[d(1);d(2);d(3);d(4);d(5);d(6)]; VLB=[0 0 0 0 0 0 0 0 0 0 0 0]; VUB=[];

x0=[1 2 3 0 1 0 0 1 0 1 0 1];

[xx,fval]=linprog(CC ,A,B,Aeq,beq,VLB,VUB,x0) 计算结果为：

x=[3.0000 5.0000 0.0000.7.0000 0.0000 1.0000 0.0000 0.0000 4.0000 0.0000 6.0000 10.0000]'

fval=136.2275

3、改建两个新料场的情形

改建两个新料场，要同时确定料场的位置（xj，yj）和运送量Xij，在同样条件下使总吨千米数最小。这是非线性规划问题。非线性规划模型为：

26

ij

minf

aa(i,j)X

j1i1

s.t.

2

i1,2,Xijd,i

j16j1,2Xije,ji1

,6

设X11X1，X21X2，X31X3，X41X4，X51X5，X61X6，X12X7，

x1X13，X22X8，X32X9，X42X10，X52X11，X62X12，y1X14

，x2X15，

y2X16。

先编写M文件liaoch.m定义目标函数：

function f=liaoch(x)

a=[1.25 8.75 0.5 5.75 3 7.27]; b=[1.25 0.75 4.75 5 6.5 7.75]; d=[3 5 4 7 6 11]; e=[20 20]; for i=1:6

s(i)=sqrt((x(13)-a(i))^2+(x(14)-b(i))^2); f1=s(i)*x(i)+f1; end f2=0; for i=7:12

s(i)=sqrt((x(15)-a(i-6))^2+(x(16)-b(i-6))^2); f2=s(i)*x(i)+f2; end

f=f1+f2;

取初值为线性规划的计算结果及临时料场的坐标： x0=[3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]', 编写主程序gying2.m

%x0=[2 2 2 2 2 2 2 2 2 2 2 2]';

x0=[3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7]'; A=[1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 ]; B=[20;20];

Aeq=[1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0];

beq=[3 5 4 7 6 11]';

vlb=[zeros(12,1);-inf;-inf;-inf;-inf]; vub=[];

[x,fval,exitflag]=fmincon('liaoch',x0,A,B,Aeq,beq,vlb,vub) 计算结果为：

x=[3.0000 5.0000 0.0707 7.0000 0 0.9293 0 0 3.9293 0 6.0000 10.0707 6.3875 4.3943 5.7511 7.1867]'

fval=105.4626 exitflag=1

即两个新料场的坐标分别为（6.3875,4.3943），（5.7511,7.1867），由料场A、

4、结果分析

由于迭代初值只能凭经验取得，所以有可能所取的初值并不合适，由此造成误差甚至错误。