数学建模-非线性规划模型用MATLAB++LINGO

约束非线性规划非现性规划的基本概念

定义如果目标函数或约束条件中至少有一个是非线性函数时的最优化问题就叫做非线性规划问题．

一般形式:

（1）其中，是定义在En上的实值函数，简记:

其它情况:

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

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

对于问题(1)，设，若存在,使得对一切，且，都有，则称X*是f(X)在D上的局部极小值点（局部最优解）．特别地当时，若则称X*是f(X)在D上的严格局部极小值点（严格局部最优解）．定义3对于问题(1),设,对任意的,都有则称X*是f(X)在D上的全局极小值点（全局最优解）．特别地当时，若，则称X*是f(X)在D上的严格全局极小值点（严格全局最优解）．

返回非线性规划的基本解法SUTM外点法SUTM内点法（障碍罚函数法）1、罚函数法2、近似规划法

返回

罚函数法

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

其中T(X,M)称为罚函数，M称为罚因子，带M的项称为罚项，这里的罚函数只对不满足约束条件的点实行惩罚：当时，满足各，故罚项=0，不受惩罚．当时，必有的约束条件，故罚项>0，要受惩罚．SUTM外点法

罚函数法的缺点是：每个近似最优解Xk往往不是容许解，而只能近似满足约束，在实际问题中这种结果可能不能使用；在解一系列无约束问题中，计算量太大，特别是随着Mk的增大，可能导致错误．1、任意给定初始点X0，取M1>1，给定允许误差，令k=1；2、求无约束极值问题的最优解，设为Xk=X(Mk)，即；3、若存在，使，则取Mk>M()令k=k+1返回（2），否则，停止迭代．得最优解.计算时也可将收敛性判别准则改为.SUTM外点法(罚函数法)的迭代步骤例题.minx2

s.t.x>=1显然本问题的最优解为x*=1用SMT外点法:

T(x,M)=x2+M[min(0,x-1)]2

=求minT(x,M).本题可由

T’(x,M)=2x+2M(x-1)=0,解得:x=M/(1+M),M趋于无穷.可知x从小于1趋于1,罚函数从外部趋于最优解.

SUTM内点法（障碍函数法）

内点法的迭代步骤例题.minx2

s.t.x>=1显然本问题的最优解为x*=1用SMT内点法:罚函数I(x,rk)=x2

–rkln(x-1),求其最值.可见,x从大于1趋于1.

近似规划法的基本思想：将问题(3)中的目标函数和约束条件近似为线性函数，并对变量的取值范围加以限制，从而得到一个近似线性规划问题，再用单纯形法求解之，把其符合原始条件的最优解作为(3)的解的近似．近似规划法每得到一个近似解后，都从这点出发，重复以上步骤．

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

近似规划法的算法步骤如下

返回例题1、用近似规划法求解下列问题。解：第一次迭代：在点x1处将g1(x),g2(x)线性化。步长限制：解（1）（2）（3）（4），得，检验，可得x2不满足原始约束。第二次迭代，减少解（1,2,3,5)，得x2=(4,3)T,f(x2)=-11.所以x2为最小值点.用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]=quaprog(...);7. [x,fval,exitflag]=quaprog(...);8. [x,fval,exitflag,output]=quaprog(...);1、二次规划例1

minf(x1,x2)=-2x1-6x2+x12-2x1x2+2x22

s.t.x1+x2≤2-x1+2x2≤2x1≥0,x2≥01、写成标准形式：2、输入命令：

H=[1-1;-12];c=[-2;-6];A=[11;-12];b=[2;2];

Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3、运算结果为：

x=0.66671.3333z=-8.2222s.t.

1.首先建立M文件fun.m,定义目标函数F（X）:functionf=fun(X);f=F(X);2、一般非线性规划

其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同.用Matlab求解上述问题，基本步骤分三步：3.建立主程序.非线性规划求解的函数是fmincon,命令的基本格式如下：

(1)x=fmincon(‘fun’,X0,A,b)

(2)x=fmincon(‘fun’,X0,A,b,Aeq,beq)

(3)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB)

(4)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’)(5)x=fmincon(‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options)

(6)[x,fval]=

fmincon(...)

(7)[x,fval,exitflag]=

fmincon(...)(8)[x,fval,exitflag,output]=

fmincon(...)输出极值点M文件迭代的初值参数说明变量上下限注意：[1]fmincon函数提供了大型优化算法和中型优化算法。默认时，若在fun函数中提供了梯度（options参数的GradObj设置为’on’），并且只有上下界存在或只有等式约束，fmincon函数将选择大型算法。当既有等式约束又有梯度约束时，使用中型算法。[2]fmincon函数的中型算法使用的是序列二次规划法。在每一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日Hessian矩阵。[3]fmincon函数可能会给出局部最优解，这与初值X0的选取有关。1、写成标准形式：

s.t.

2x1+3x26

s.tx1+4x25x1,x20例22、先建立M-文件fun3.m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^23、再建立主程序youh2.m：

x0=[1;1];A=[23;14];b=[6;5];

Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)4、运算结果为：

x=0.76471.0588

fval=-2.02941．先建立M文件fun4.m,定义目标函数:

functionf=fun4(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);例32．再建立M文件mycon.m定义非线性约束：

function[g,ceq]=mycon(x)g=x(1)+x(2);

ceq=[1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10];3．主程序youh3.m为:x0=[-1;1];A=[];b=[];Aeq=[11];beq=[0];vlb=[];vub=[];[x,fval]=fmincon('fun4',x0,A,b,Aeq,beq,vlb,vub,'mycon')3.运算结果为：

x=-1.22501.2250

fval=1.8951

例4

1．先建立M-文件fun.m定义目标函数:functionf=fun(x);f=-2*x(1)-x(2);2．再建立M文件mycon2.m定义非线性约束：

function[g,ceq]=mycon2(x)g=[x(1)^2+x(2)^2-25;x(1)^2-x(2)^2-7];

3.主程序fxx.m为:x0=[3;2.5];VLB=[00];VUB=[510];[x,fval,exitflag,output]=fmincon('fun',x0,[],[],[],[],VLB,VUB,'mycon2')4.运算结果为:x=4.00003.0000fval=-11.0000exitflag=1output=iterations:4

funcCount:17

stepsize:1algorithm:[1x44char]

firstorderopt:[]

cgiterations:[]应用实例：供应与选址

某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系a，b表示，距离单位：千米）及水泥日用量d(吨)由下表给出。目前有两个临时料场位于A(5,1)，B(2,7)，日储量各有20吨。假设从料场到工地之间均有直线道路相连。（1）试制定每天的供应计划，即从A，B两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。（2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数有多大？（一）、建立模型

记工地的位置为(ai，bi)，水泥日用量为di，i=1,…,6;料场位置为(xj，yj)，日储量为ej，j=1,2；从料场j向工地i的运送量为Xij。当用临时料场时决策变量为：Xij，当不用临时料场时决策变量为：Xij，xj，yj。（二）使用临时料场的情形

使用两个临时料场A(5,1)，B(2,7).求从料场j向工地i的运送量为Xij，在各工地用量必须满足和各料场运送量不超过日储量的条件下，使总的吨千米数最小，这是线性规划问题.线性规划模型为：设X11=X1,X21=X2,,X31=X3,X41=X4,X51=X5,,X61=X6X12=X7,X22=X8,,X32=X9,X42=X10,X52=X11,,X62=X12

编写程序gying1.mgying1.mcleara=[1.258.750.55.7537.25];b=[1.250.754.7556.57.75];d=[3547611];x=[52];y=[17];e=[2020];fori=1:6forj=1:2

aa(i,j)=sqrt((x(j)-a(i))^2+(y(j)-b(i))^2);endendCC=[aa(:,1);aa(:,2)]';A=[111111000000;000000111111];B=[20;20];Aeq=[100000100000%从第一＼二料场运到工地一的料

010000010000001000001000000100000100

000010000010000001000001];beq=[d(1);d(2);d(3);d(4);d(5);d(6)];VLB=[000000000000];VUB=[];x0=[123010010101];[xx,fval]=linprog(CC,A,B,Aeq,beq,VLB,VUB,x0)计算结果为：x=[3.00005.00000.00007.00000.00001.00000.00000.00004.00000.00006.000010.0000]’fval=136.2275（三）改建两个新料场的情形

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

用lingo解此题MODEL:设X11=X1,X21=X2,,X31=X3,X41=X4,X51=X5,,X61=X6X12=X7,X22=X8,,X32=X9,X42=X10,X52=X11,,X62=X12

x1=X13,y1=X14,x2=X15,y2=X16

（1）先编写M文件liaoch.m定义目标函数。(2)取初值为线性规划的计算结果及临时料场的坐标:x0=[35070100406105127]';编写主程序gying2.m.functionf=liaoch(x)

a=[1.258.750.55.7537.25];

b=[1.250.754.7556.57.75];

d=[3547611];

e=[2020];

f1=0;

fori=1:6

s(i)=sqrt((x(13)-a(i))^2+(x(14)-b(i))^2);

f1=s(i)*x(i)+f1;

end

f2=0;

fori=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;

clear%x0=[35070100406105127]';%x0=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435.75117.1867]';%x0=[3.00005.00000.30947.00000.01080.6798003.690605.989210.32025.53694.91945.82917.2852]';x0=[35471000005115.63484.86877.24797.7499]';A=[11111100000000000000001111110000];B=[20;20];原程序gying2.mAeq=[100000100000000001000001000000000010000010000000000100000100000000001000001000000000010000010000];beq=[3547611]';vlb=[zeros(12,1);-inf;-inf;-inf;-inf];vub=[];[x,fval,exitflag]=fmincon('liaoch',x0,A,B,Aeq,beq,vlb,vub)(3)计算结果为：x=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435.75117.1867]’fval=105.4626exitflag=1(4)若修改主程序gying2.m,取初值为上面的计算结果:x0=[3.00005.00000.07077.000000.9293003.929306.000010.07076.38754.39435.75117.1867]’得结果为:x=[3.00005.00000.30947.00000.01080.6798003.690605.989210.32025.53694.91945.82917.2852]’fval=103.4760exitflag=1总的吨千米数比上面结果略优.(5)若再取刚得出的结果为初值,却计算不出最优解.(6)若取初值为:x0=[35471000005115.63484.86877.24797.7499]',

则计算结果为:x=[3.00005.00004.00007.00001.0000000005.000011.00005.69594.92857.25007.7500]’fval=89.8835exitflag=1总的吨千米数89.8835比上面结果更好.

gying2.m通过此例可看出fmincon函数在选取初值上的重要性.用lingo解此题的方法:Lingo的使用－－数学函数LINGO提供了大量的标准数学函数：@abs(x)返回x的绝对值@sin(x)返回x的正弦值，x采用弧度制@cos(x)返回x的余弦值@tan(x)返回x的正切值@exp(x)返回常数e的x次方@log(x)返回x的自然对数@lgm(x)返回x的gamma函数的自然对数@sign(x)如果x<0返回-1；否则，返回1@floor(x)返回x的整数部分。当x>=0时，返回不超过x的最大整数；当x<0时，返回不低于x的最大整数。@smax(x1,x2,…,xn)返回x1，x2，…，xn中的最大值@smin(x1,x2,…,xn)返回x1，x2，…，xn中的最小值例4.3

给定一个直角三角形，求包含该三角形的最小正方形。

LINGO代码如下：model:Titletrangle;sets:object/1..3/:f;endsetsdata:a,b=3,4;!两个直角边长，修改很方便;Enddataf(1)=a*@sin(x);f(2)=b*@cos(x);f(3)=a*@cos(x)+b*@sin(x);min=@smax(f(1),f(2),f(3));@bnd(0,x,1.57);end例题:

一奶制品加工厂用牛奶生产A1,A2两种奶制品，1桶牛奶可以在甲车间用12小时加工成3公斤A1，或者在乙车间用8小时加工成4公斤A2。根据市场需求，生产的A1,A2全部能售出，且每公斤A1获利24元，每公斤A2获利16元。现在加工厂每天能得到50桶牛奶的供应，每天正式工人总的劳动时间480小时，并且甲车间每天至多能加工100公斤A1，乙车间的加工能力没有限制。试为该厂制订一个生产计划，使每天获利最大，并进一步讨论以下3个附加问题：1）若用35元可以买到1桶牛奶，应否作这项投资？若投资，每天最多购买多少桶牛奶？2）若可以聘用临时工人以增加劳动时间，付给临时工人的工资最多是每小时几元？3）由于市场需求变化，每公斤A1的获利增加到30元，应否改变生产计划？模型代码如下：max=72*x1+64*x2;x1+x2<=50;12*x1+8*x2<=480;3*x1<=100;求解这个模型,结果如下:Globaloptimalsolutionfoundatiteration:0Objectivevalue:3360.000

VariableValueReducedCostX120.000000.000000X230.000000.000000RowSlackorSurplusDualPrice13360.0001.00000020.00000048.0000030.0000002.000000440.000000.000000Rangesinwhichthebasisisunchanged:

ObjectiveCoefficientRanges

CurrentAllowableAllowable

VariableCoefficientIncreaseDecrease

X172.0000024.000008.000000

X264.000008.00000016.00000

RighthandSideRanges

RowCurrentAllowableAllowable

RHSIncreaseDecrease

250.0000010.000006.666667

3480.000053.3333380.00000

4100.0000INFINITY40.00000做灵敏性分析，结果如下:用lingo解此题

MODEL:

TitleLocationProblem;

sets:

demand/1..6/:a,b,d;

supply/1..2/:x,y,e;

link(demand,supply):c;

endsets

data:

!locationsforthedemand(需求点的位置);

a=1.25,8.75,0.5,5.75,3,7.25;

b=1.25,0.75,4.75,5,6.5,7.75;

!quantitiesofthedemandandsupply（供需量）;

d=3,5,4,7,6,11;e=20,20;

enddata

init:

!initiallocationsforthesupply（初始点）;

x,y=5,1,2,7;

endinit

!Objectivefunction（目标）;

[OBJ]min=@sum(link(i,j):c(i,j)*((x(j)-a(i))^2+(y(j)-b(i))^2)^(1/2));

!demandconstraints（需求约束）;

@for(demand(i):[DEMAND_CON]@sum(supply(j):c(i,j))=d(i););

!supplyconstraints（供应约束）;

@for(supply(i):[SUPPLY_CON]@sum(demand(j):c(j,i))<=e(i););

@for(supply:@bnd(0.5,X,8.75);@bnd(0.75,Y,7.75););

END

运行的结果为:X(1)5.695967X(2)7.249997Y(1)4.928560Y(2)7.750000Feasiblesolutionfound.Objectivevalue:89.88349Totalsolveriterations:81

1.3应用举例(二)---整数规化解先看有多少种裁料方案，再进行组合和选择。方案：

例1合理利用线材问题现要做一百套钢管,每套要长为2.9m、2.1m和1.5m的钢管各一根。已知原料长7.4m，问应如何下料，使用的原料最省。

设用方案Ⅰ,Ⅱ,…,Ⅷ分别裁原料钢管x1,x2,…,x8根,则：Minz=x1+x2+x3+x4+x5+x6+x7+x82x1+x2+x3+x4≥

1002x2+x3+3x5+2x6+x7≥100x1+x3+3x4+2x6+3x7+4x8≥100x1,x2,x3,x4,x5,x6,x7,x8≥0

上题中,如果每套1.8m的钢管要70根,要求使用的切割模式不超过3种.问应如何下料，使用的原料最省。

解:设xi—按第i种模式切割的原料钢管根数(i=1,2,3)r1i,r2i,r3i,r4i—第i种切割模式下,每根原料生产2.9m、2.1m和1.5m,1.8m的钢管的数量.目标函数(总根数)

(套数约束)按第i种切割模式下,每根钢管的长度限制:这是一个整数非线性规划模型，用LINGO要运行很长时间也难以得到最优解。因三种切割模式的排列顺序无关紧要，所以不妨增加以下约束：又钢管的总根数有明显的上界和下界。首先，原料钢管的总根数不可能少于其次，考虑第一种切割模式下只生产2.9m钢管，一根原料钢管切割成2根2.9m钢管，100套钢管需要50根原料；如此可计算出钢管的上界：

所以可以增加以下约束：Lingo源程序model:title钢管下料；min=x1+x2+x3;x1*r11+x2*r12+x3*r13>=100;x1*r21+x2*r22+x3*r23>=100;x1*r31+x2*r32+x3*r33>=100;x1*r41+x2*r42+x3*r43>=70;x1+x2+x3>=105;x1+x2+x3<=152;2.9*r11+2.1*r21+1.5*r31+1.8*r41>=5.9;2.9*r12+2.1*r22+1.5*r32+1.8*r42>=5.9;2.9*r13+2.1*r23+1.5*r33+1.8*r43>=5.9;2.9*r11+2.1*r21+1.5*r31+1.8*r41<=7.4;2.9*r12+2.1*r22+1.5*r32+1.8*r42<=7.4;2.9*r13+2.1*r23+1.5*r33+1.8*r43<=7.4;x1>=x2;x2>=x3;@gin(x1);@gin(x2);@gin(x3);@gin(r11);@gin(r21);@gin(r31);@gin(r41);@gin(r12);@gin(r22);@gin(r32);@gin(r42);@gin(r13);@gin(r23);@gin(r33);@gin(r43);end运行结果如下：Feasiblesolutionfoundatiteration:26565

ModelTitle:钢管下料；

VariableValue

X172.00000

R111.000000

X229.00000

R120.000000

X314.00000

R132.000000

R211.000000

R221.000000

R230.000000

R310.000000

R323.000000

R331.000000

R411.000000

R420.000000

R430.000000Lingo源程序二model:Title钢管下料-最小化钢管根数的LINGO模型;SETS:NEEDS/1..4/:LENGTH,NUM; !定义基本集合NEEDS及其属性LENGTH,NUM;CUTS/1..3/:X; !定义基本集合CUTS及其属性X;PATTERNS(NEEDS,CUTS):R;!定义派生集合PATTERNS（这是一个稠密集合）及其属性R;ENDSETSDATA: LENGTH=2.92.11.51.8; NUM=100100

10070; CAPACITY=7.4;ENDDATAmin=@SUM(CUTS(I):X(I)); !目标函数;Lingo源程序二@FOR(NEEDS(I):@SUM(CUTS(J):X(J)*R(I,J))>NUM(I));!满足需求约束;@FOR(CUTS(J):@SUM(NEEDS(I):LENGTH(I)*R(I,J))<CAPACITY);!合理切割模式约束;@FOR(CUTS(J):@SUM(NEEDS(I):LENGTH(I)*R(I,J))>CAPACITY-@MIN(NEEDS(I):LENGTH(I))); !合理切割模式约束;@SUM(CUTS(I):X(I))>105;@SUM(CUTS(I):X(I))<115;!人为增加约束;@FOR(CUTS(I)|I#LT#@SIZE(CUTS):X(I)>X(I+1)); !人为增加约束;@FOR(CUTS(J):@GIN(X(J)));@FOR(PATTERNS(I,J):@GIN(R(I,J)));end运行结果如下：Localoptimalsolutionfoundatiteration:64408

Objectivevalue:113.0000

ModelTitle:钢管下料-最小化钢管根数的LINGO模型

VariableValueReducedCost

CAPACITY