1 整数规划的MATLAB求解方法
(一) 用MATLAB求解一般混合整数规划问题
由于MATLAB优化工具箱中并未提供求解纯整数规划和混合整数规划的函数,因而需要自行根据需要和设定相关的算法来实现。现在有许多用户发布的工具箱可以解决该类问题。这里我们给出开罗大学的Sherif和Tawfik在MATLAB Central上发布的一个用于求解一般混合整数规划的程序,在此命名为intprog,在原程序的基础上做了简单的修改,将其选择分枝变量的算法由自然序改造成分枝变量选择原则中的一种,即:选择与整数值相差最大的非整数变量首先进行分枝。intprog函数的调用格式如下:
[x,fval,exitflag]=intprog(c,A,b,Aeq,beq,lb,ub,M,TolXInteger) 该函数解决的整数规划问题为:
mins.t.
fcTxAxbAeqxbeqlbxubxi0
(i1,2,,n)
xj取整数(jM)
在上述标准问题中,假设x为n维设计变量,且问题具有不等式约束m1个,等式约束m2个,那么:c、x均为n维列向量,b为m1维列向量,beq为m2维列向量,A为m1n维矩阵,Aeq为m2n维矩阵。
在该函数中,输入参数有c,A,b,Aeq,beq,lb,ub,M和TolXInteger。其中c为目标函数所对应设计变量的系数,A为不等式约束条件方程组构成的系数矩阵,b为不等式约束条件方程组右边的值构成的向量。Aeq为等式约束方程组构成的系数矩阵,beq为等式约束条件方程组右边的值构成的向量。lb和ub为设计变量对应的上界和下界。M为具有整数约束条件限制的设计变量的序号,例如问题中设计变量为x1,x2,,x6,要求x2,x3和x6为整数,则M=[2;3;6];若要求全为整数,则M=1:6,或者M=[1;2;3;4;5;6]。TolXInteger为判定整数的误差限,即若某数x和最邻近整数相差小于该误差限,则认为x即为该整数。
在该函数中,输出参数有x, fval和exitflag。其中x为整数规划问题的最优解向量,fval为整数规划问题的目标函数在最优解向量x处的函数值,exitflag为函数计算终止时的状态指示变量。 例1 求解整数规划问题:
fx1x2max
s.t. 4x12x21
4x12x211 2x1
2
x1,x20,且取整数值
算法:
c=[-1;-1]; A=[-4 2;4 2;0 -2]; b=[-1;11;-1]; lb=[0;0]; M=[1;2]; Tol=1e-8;
%均要求为整数变量 %判断是否整数的误差限
%求解原问题松弛线性规划
[x,fval]=linprog(c,A,b,[],[],lb,[])
[x1,fval1]=intprog(c,A,b,[],[],lb,[],M,Tol) %求最优解整数解 结果:
x =
%松弛线性规划问题的最优解
1.5000 2.5000 fval = -4.0000 x1 = 2 1 fval2 = -3.0000
%整数规划的最优解
(二) 用MATLAB求解0-1规划问题
在MATLAB优化工具箱中,提供了专门用于求解0-1规划问题的函数bintprog,其算法基础即为分枝界定法,在MATLAB中调用bintprog函数求解0-1规划时,需要遵循MATLAB中对0-1规划标准性的要求。 0-1规划问题的MATLAB标准型
mins.t.
fcTxAxbAeqxbeqx0,1
在上述模型中,目标函数f 需要极小化,以及需要满足的约束条件,不等式约束一定要化为形式为“”。
假设x为n维设计变量,且问题具有不等式约束m1个,等式约束m2个,那么:
c、x均为n维列向量,b为m1维列向量,beq为m2维列向量,A为m1n维矩阵,
Aeq为m2n维矩阵。
如果不满足标准型的要求,则需要对原问题进行转化,化为标准型之后才能使用相关函数,标准化的方法和线性规划中的类似。 0-1规划问题的MATLAB求解函数
MATLAB优化工具箱中求解0-1规划问题的命令为bintprog bintprog的调用格式
x = bintprog(f) x = bintprog(f,A,b) x = bintprog(f,A,b,Aeq,beq) x = bintprog(f,A,b,Aeq,beq,x0) x = bintprog(f,A,b,Aeq,Beq,x0,options) [x,fval] = bintprog(...) [x,fval,exitflag] = bintprog(...) [x,fval,exitflag,output] = bintprog(...)
命令详解
1)x = bintprog(f)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
fcTxx0,1
2)x = bintprog(c,A,b)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
3)x = bintprog (c,A,b,Aeq,beq)
fcTxAxbx0,1
该函数调用格式求解如下形式的0-1规划问题:
mins.t.
fcTxAxbAeqxbeqx0,1
4)x = bintprog (c,A,b,Aeq,beq,x0)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
fcTxAxbAeqxbeqx0,1
在前一个调用格式的基础上同时设置求解算法的初始解为x0,如果初始解x0不在0-1规划问题的可行域中,算法将采用默认的初始解 5)x = bintprog (c,A,b,Aeq,beq,x0,options)
用options指定的优化参数进行最小化。可以使用optimset来设置这些参数
上面的函数调用格式仅设置了最优解这一输出参数,如果需要更多的输出参数,则可以参照下面的调用格式:
[x,fval] = bintprog(...)
在优化计算结束之时返回整数规划问题在解x处的目标函数值fval
[x,fval,exitflag] = bintprog(...)
在优化计算结束之时返回exitflag值,描述函数计算的退出条件。
[x,fval,exitflag,output] = bintprog(...) 在优化计算结束之时返回结构变量output 例2:求解0-1规划问题
nn
fEijxij
max
i1j1
20n22s.t. xij1 i1,2,,nEj1
21n
xij1 j1,2,,n22
i1
x0或1( 1,2,3,4) i1,2,,n;j1,2,,n ij
123326
152923 133124
163223
为了程序的可读性,我们用一维下标来表示设计变量,即用x1~x4表示x11~x14,用x5~x8表示x21~x24,用x9~x12表示x31~x34,用x13~x16表示x41~x44,于是约束条件和目标函数分别为:
x1x2x3x41
xxxx1
6785
x9x10x11x121
x13x14x15x161
x1x5x9x131 xxxx1
61014
2
x3x7x11x151
x4x8x12x161xi0,1 (i1,2,,16)
fE11x1E12x2E13x3E14x4E21x5E22x6E44x16
算法:
c=[20;12;33;26;22;15;29;23;21;13;31;24;22;16;32;23]; Aeq=[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 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 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0; 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0; 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1; ]; beq=ones(1,8);
[x,fval]=bintprog(c,[],[],Aeq,beq);
B=reshape(x,4,4); %由于x是一列元素,为了使结果更加直观,故排成与效率矩阵E相
对应的形式
B' fval 结果:
Optimization terminated. ans =
0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 fval = 85
整数规划的应用——组件配套问题
某机械产品需要由三个工厂开工一起生产后组装完成。每件机械需要4个组件1和3个组件2。生产这两种组件需要消耗两种原材料A和B。已知这两种原材料的供应量分别为400kg和600kg。
由于三个工厂的生产条件和拥有设备工艺条件不同,每个工厂生产组件的能力和原材料的消耗也不尽相同,且每个工厂开工一次都是配套生产一定数量的组件1和组件2,其具体的数据如表所示。
表11-2 各工厂生产能力和消耗原材料的数据表
现在的最优化问题是,这三个工厂应当如何安排生产,才能使该产品的配套数达到最大?
(Ⅰ)组件配套问题的建模
设x1、x2和x3是三个工厂分别开工的次数,将其作为该问题的设计变量。由于每个工厂开工一次都是配套生产,故每次开工消耗的原材料一定,且生产的组件数也是一定的。在这个假设的前提之下,我们可以得出该问题的目标函数和约束条件。
因为原材料的总量是有限的,根据工厂的开工次数,可得工厂1将消耗A材料9x1,工厂2将消耗A材料6x2,工厂3将消耗A材料4x3,故有约束条件:
9x16x24x3400
同理,对于材料B的总量约束条件可以表达为:7x110x29x3600 然后再来分析三个工厂零件生产的情况,三个工厂生产的组件1的数量分别为8x1、7x2和9x3,故组件1的总数为:Q18x17x29x3
同理,组件2的总数为:Q26x19x25x3
下一步是分析目标函数,该问题要求的不是生产的各种组件总数最多,而是要求产品的配套数量最大,即能组成的机械数目最多。问题中已经给出了该种机械中两种组件的配比为4:3,故组件1所能成套的数目T1满足约束条件:
T1
Q18x17x29x3
44
同理,组件2所能成套的数目T2满足约束条件:T2
Q26x19x25x3
33
因而,所能组成的成品机械的数目f应该为T1和T2中的较小值,即:
fmin(T1,T2)
那么,我们求解的目标即是使得f能够尽可能大,可以看出,这种形式并不是有关设计变量的线性函数,我们需要对其进行转化,此时我们可以令一个人工设计变量等于目标函数值,则有:x4min(T1,T2)
在该假设下,一定满足关系式:T1x4且T2x4 结合约束关系,对T1的约束可以表示为:x4T1
Q18x17x29x3
44
同理,对T2的约束可以表示为:x4T2
Q26x19x25x3
33
对T1的上述关系进行整理,可以得到关系:8x17x29x34x40 同理对T2也可以得到不等式关系为:6x19x25x33x40
上面两个式子即为由组件的配比数得到的关于组件数目的约束条件,此时问题的目标就是如何使得x4取到最大值,由于开工的次数一定是一个非负整数,故也是一个约束条件,决定了该问题是一个纯整数规划问题。结合前面给出的原材料约束,可以得到如下的数学模型:
max fx4
s.t. 9x6x4x400
123
7x110x29x3600
8x17x29x34x40 6x19x25x33x40 xi0且取整数值i1,2,3,4
(Ⅱ)组件配套问题的求解
利用§8.1节中给出的MATLAB函数对此问题求解,代码和运行结果如下:
算法:
%目标函数所对应的设计变量的系数,为转化为极小,故取原目标函数的相反数 f=[0;0;0;-1]; %不等式约束 A=[ 9 6 4 0;
7 10 9 0; -8 -7 -9 4; -6 -9 -5 3]; b=[400;600;0;0];
%边界约束,由于无上界,故设置ub=[Inf;Inf;Inf;Inf]; lb=[0;0;0;0]; ub=[Inf;Inf;Inf;Inf];
%所有变量均为整数变量,故将所有序号组成向量M M=[1;2;3;4];
%判定为整数的误差限 Tol=1e-8;
%求最优解x和目标函数值fval,并返回状态指示 [x,fval,exitflag]=intprog(f,A,b,[],[],lb,ub,M,Tol) 结果:
x = %最优解向量x 18 15 36 141
fval = %在最优解向量x处,原目标函数值的相反数 -141.000
exitflag= %算法终止指示变量,说明问题收敛到了最优解x 1
由运行结果可知,工厂1、2和3需要分别开工18、15和36次,这样所能生产出来的成套的机械为141件。
2 多目标规划的MATLAB求解方法
(一) 多目标规划的MATLAB求解
由于多目标规划中的求解涉及到的方法非常多,故在MATLAB中可以利用
不同的函数进行求解,例如在评价函数法中我们所得最后的评价函数为一线性函数,且约束条件也为线性函数,则我们可以利用MATLAB优化工具箱中提供的linprog函数进行求解,如果我们得到的评价函数为非线性函数,则可以利用MATLAB优化工具箱中提供的fmincon函数进行求解,如果我们采用最大最小法进行求解,则可以利用MATLAB优化工具箱中提供的fminimax函数进行求解。下面我们就结合理论求解的几种方法,讲解一下典型多目标规划问题的MATLAB求解方法。 例1 利用理想点法求解
min f1(x)2x13x2min f(x)5x3x
212
s.t. 3x12x212 xx8
12
x1,x20
我们首先根据评价函数法中的理想点法的理论基础,按照理想点法的求解思路分别对两个单目标规划问题P1,P2进行求解:
min f1(x)2x13x2min f2(x)5x13x2
s.t. 3x2x12s.t. 3x2x121212
P1,P2
xx8 xx81212 x1,x20 x1,x20
求解P1的MATLAB的代码和相应的运行结果为:
算法:
c=[2;-3]; A=[3 2;1 1]; b=[12;8] lb=[0;0]
[x,fval]=linprog(c,A,b,[],[],lb,[]) 结果:
x = 0.0000 6.0000
-18.0000
*
于是可知。当x10,6时,单目标线性规划P 1的最优函数值为f118。
T
求解P2的MATLAB的代码和相应的运行结果为:
算法:
c=[-5;-3]; A=[3 2;1 1]; b=[12;8] lb=[0;0]
[x,fval]=linprog(c,A,b,[],[],lb,[]) 结果:
Optimization terminated. x = 4.0000 0.0000 fval = -20.0000
于是可知,当x24,0时,单目标线性规划P2的最优函数值为f2*20。
T
由上述两个单目标线性规划的求解结果可知x2x2,因而
f
*1
,f2*18,20是一个不可能达到的理想点,因而我们构造如下评价函数:
minhf(x)
xR
f1(x)182f2(x)202
2x13x21825x13x202
构造描述该评价函数的M-函数文件objfun.m如下: function f=objfun(x)
f=sqrt((2*x(1)-3*x(2)+18)^2+(5*x(1)+3*x(2)-20)^2); 然后用非线性规划的方式求解上述问题:
算法:
A=[3 2;1 1]; b=[12;8];
x0=[0;0];
[x,fval,exitflag]=fmincon(@objfun,x0,A,b,[],[],lb,[]) 结果:
Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1
x = 0.0235 5.9647
fval = 1.9941
exitflag = 5
由运行结果可知在该评价函数标准之下,问题的最优解为:
x*0.0235,5.9647此时,各目标函数的取值为:f1*17.8471,f2*18.0118。
T
它与理想点f1*,f2*18,20在评价函数标准下的最小距离为1.9941。 例2 利用评价函数中的线性加权和法求解如下多目标规划问题:
22
min f1(x)x12x2x3222min f2(x)x12x23x3
s.t. x1x2x33 x,x,x0
123
其中权系数为10.8,20.2。 建立线性加权和法的评价函数为:
2222
minhf(x)1x12x2x32x122x23x3
将相应的权系数代入上式即整理出目标函数f(x)为:
22f(x)x121.2x21.4x3
于是建立目标函数的M-函数文件objfun.m:
function f=objfun(x)
f=x(1)^2+1.2*x(2)^2+1.4*x(3)^2;
由于目标函数非线性函数且具有线性等式约束和边界约束,因而我们调用MATLAB中求解非线性规划的fmincon函数对此问题进行求解,同时注意如果只考虑第一个目标函数,由这种特殊形式,即在设计变量的和为一定值,需要求其平方和的最小值时,最优解必然是当这几个设计变量的值相等时取得,于是我们可以将这个解设为问题的初始点,开始迭代。
算法:
Aeq=[1 1 1]; beq=[3]; lb=[0;0;0]; x0=[1;1;1];
[x,fval]=fmincon(@objfun,x0,[],[],Aeq,beq,lb,[]) 结果:
No active inequalities. x = 1.1776 0.9812 0.8412 fval = 3.5327
*x11.1776*,f(x*)3.5327 结果说明,问题的最优解为:x*x20.9812*x30.8412
我们在求解多目标规划问题时,可以采用评价函数法中的最大最小法,而MATLAB也为这种方法提供了专门的求解函数fminimax,在讲解这方面的例题之前,我们首先介绍一下MATLAB优化工具箱中所提供的最大最小法的求解函数fminimax。
最大最小法问题的MATLAB标准形式为:
min max fi(x)
i
x
s.t. c(x)0 c(x)0eq
Axb Axb
eqeq
lbxub
函数fminimax的调用方式和其他的最优化函数类似,其中所涉及的输入参数和输出参数的含义与非线性规划的求解函数fmincon类似,使用方法也基本相同,细节问题读者可以参考MATLAB的帮助文件。 例3 求解最大最小问题:
min max fi(x)i
x
2
s.t. f1(x)2x12x248x140x23042 f2(x)x123x2
f3(x)x13x218
f4(x)x1x2 f5(x)x1x28
首先建立描述目标函数的M-函数文件objfun.m,注意到一共有五个目标函数,所求目标为这五个函数最大值中的最小值,代码如下:
function f = objfun(x)
f(1)= 2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304; f(2)= -x(1)^2 - 3*x(2)^2; f(3)= x(1) + 3*x(2) -18; f(4)= -x(1)- x(2); f(5)= x(1) + x(2) - 8;
然后设置求解的初始点为x0=[0;0],调用fminimax求解该问题。
算法:
x0 = [0; 0];
[x,fval,maxfval] = fminimax(@objfun,x0) 结果:
Local minimum possible. Constraints satisfied.
fminimax stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints were satisfied to within the default value of the constraint tolerance.
Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1 5 x = 4.0000 4.0000 fval =
0.0000 -64.0006 -1.9999 -8.0000 -0.0000 maxfval = 2.6897e-008
上述结果说明当x14,x24时,目标函数fi(x) i1,2,,5的最大值达到最小,这一组的函数值为0.0000,-64.0006,-1.9999,-8.0000,-0.0000,于是最大值为0。
多目标规划的应用——投资问题(全国大学生数学建模竞赛试题)
假设市场上有n种资产,比如股票、债券等可以供投资者选择,某公司有数额为M的一笔相当大的资金可用作一个时间的投资。通过财务人员对Si种资产进行评估,大概可以估算出在这一时期内购买资产Si的平均收益率为ri,并预测出购买Si的损失率为qi。考虑到投资越分散,总的风险越小,公司决定当用这笔资金购买若干种资产时,总体风险可用所投资的Si中的最大一个风险来度量。
购买Si要付交易费,费率为pi,并且当购买额不超过给定值ui时,交易费按购买ui计算(不买当然无须付费)。另外,假定同期银行存款利率是r0,且既无交易费又无风险(r0=5%)。已知n4时的相关数据如下表所示:
表1 投资各种资产的参数值
试给该公司设计一种投资组合方案,即用给定的资金M,有选择地购买若干种资产或存银行生息,使净收益尽可能大,而总体风险尽可能小。 (Ⅰ)投资问题的建模
为了建立数学模型,首先对模型进行一些必要的假设及符号说明: (1) 模型的假设
① 在一个时期内所给出的ri,qi,pi保持不变。
②在一个时间内所购买的各种资产(如股票、证券等)不进行买卖交易,即在买入后不再卖出。
③ 每种投资是否收益是相互独立的。
④ 在投资过程中,无论盈利与否必须先付交易费。 (2)符号说明
M (元):公司现有投资总金额;
Si,i0,1,,n:欲购买的第i种资产种类(其中i0表示存入银行): xi,i0,1,,n:公司购买Si的金额; ri,i0,1,,n:公司购买Si的平均收益率; qi,i0,1,,n:公司购买Si的平均损失率;
pi,i0,1,,n:公司购买Si超过ui时所付交易费率。
下面来建立模型。设购买Si的金额为xi,所付的交易费ci(xi),则
0 xi0
pu 0xu(i1,2,,n)iiii
ci(xi)
pixi xiuic0(x0)0
由于投资额M相当大,所以总可以假定对每个Si的投资xiui。这时上面的大括号公式可简化为:ci(xi)pixi (i1,2,,n)
对Si投资的净收益为:Ri(xi)rixici(xi)ripixi 对Si的风险为:Qi(xi)qixi
对Si投资所需资金为投资金额xi与所需的手续费ci(xi)之和,即:
fi(xi)xici(xi)1pixi
当购买Si的金额为xi (i0,1,2,,n),投资组合x(x0,x1,,xn)的净收益总额为:
R(x)Ri(xi)
i0n
整体风险为:Q(x)maxQi(xi)maxqi(xi)
1in
1in
资金约束为:fi(xi)M
i0
n
根据题目要求,以净收益总额R(x)最大,而整体风险Q(x)最小为目标建立模型如下:
n
min rpxmaxqxiiiiii0
n
s.t. 1pixiM
i0
xi0,i1,2,,n
很显然,这是一个多目标规划问题。 (Ⅱ)投资问题的求解
在此我们采用主要目标法对该问题进行求解,即根据问题的实际情况,确
定一个目标为主要目标,而把其余目标作为次要目标,并且根据经验,选取一定的界限值。这样就可以把次要目标作为约束来处理,于是就将原来的多目标问题转化为一个在新的约束下的单目标最优化问题。
在上述例子中如果以收益为主要目标,则可以固定风险水平,给定风险一个界限a,讲问题转化称为求最大风险不超过a时的最大收益,即下面的线性规划模型:
n
max ripixi
i0
s.t. qixiMa i1,2,,n
n
1pxM
ii
i0
xi0,i1,2,,n
(1)
若投资者希望总盈利至少达到水平K以上,则可以在风险最小的情况下寻找相应的投资组合,从而将原模型转化成为下列的线性规划模型进行求解:
min maxqixiix
n
s.t. ripixiKi0
n
1pixiM
i0
xi0,i1,2,,n
(2)
根据上面的分析,我们利用主要目标法建立了该问题的多目标规划模型,进而转化成为了线性规划模型,下面我们利用MATLAB对此问题进行求解并进行投资分析。将n4时的数据代入模型(1)中,我们可以得到MATLAB的标准型为:
min f0.05x00.27x10.19x20.185x30.185x4
s.t. x1.01x1.02x1.045x1.065x1
01234
0.025x1a
0.015x2a
0.055xa
3
0.026x4a
xi0,i1,2,3,4
由于a是任意给定的风险度,因而实际上没有一个选择的准则,不同的投资者可能有自身不同的判断,因而使得a的取值不同。我们在求解的过程中不妨用试探的方法,从a0开始,以步长a0.001进行搜索,通过实验来分析和总结风险度a和收益Q之间的数量关系。利用MATLAB编程求解。
算法:
clear i=1; a(i)=0; x=cell(1,101); while(1.1-a(i))>1 i=i+1; a(i)=a(i-1)+0.001;
f=[-0.05 -0.27 -0.19 -0.185 -0.185]; Aeq=[1 1.01 1.02 1.045 1.065]; beq=[1]; A=[0 0.025 0 0 0; 0 0 0.015 0 0; 0 0 0 0.055 0; 0 0 0 0 0.026]; b=[a(i);a(i);a(i);a(i)]; lb=[0;0;0;0;0];
[y,fval(i)]=linprog(f,A,b,Aeq,beq,lb,[]); x{i}=y; Q(i)=-fval(i); end plot(a,Q,'.'); xlabel('a'); ylabel('Q'); 结果:
Q
a
风险度与收益关系图
模型(1)的结果图分析:
(1) 由上图可知,收益随风险增大而增大,就是说,风险越大,收益也越大。 (2) 由数据表可以看出,投资越分散,投资者承担的风险越小。这与实际情况相符,就是说,冒险的投资者会出现集中投资的情况,保守的投资者则尽量分散投资。
(3) 曲线上的任一点都表示该风险水平的最大可能收益和该收益要求的最小风险。投资者应根据对不同风险的承受能力,选择该风险水平下的最优投资组合。
(4) 由局部放大上图可以看出,在a(7)0.0060附近有一个转折点。在这一点左边,风险增加很少时,利润增长很快;在这一点右边,风险增加很大时,利润增长很缓慢。所以对于风险和收益没有特殊偏好的投资者来说,应该选择曲线的该转折点作为最优投资组合,大约是a(7)0.0060,Q*Q(7)0.2019,所对应投资方案为:
x[0,0.2400,0.4000,0.1091,0.2212]T
1 整数规划的MATLAB求解方法
(一) 用MATLAB求解一般混合整数规划问题
由于MATLAB优化工具箱中并未提供求解纯整数规划和混合整数规划的函数,因而需要自行根据需要和设定相关的算法来实现。现在有许多用户发布的工具箱可以解决该类问题。这里我们给出开罗大学的Sherif和Tawfik在MATLAB Central上发布的一个用于求解一般混合整数规划的程序,在此命名为intprog,在原程序的基础上做了简单的修改,将其选择分枝变量的算法由自然序改造成分枝变量选择原则中的一种,即:选择与整数值相差最大的非整数变量首先进行分枝。intprog函数的调用格式如下:
[x,fval,exitflag]=intprog(c,A,b,Aeq,beq,lb,ub,M,TolXInteger) 该函数解决的整数规划问题为:
mins.t.
fcTxAxbAeqxbeqlbxubxi0
(i1,2,,n)
xj取整数(jM)
在上述标准问题中,假设x为n维设计变量,且问题具有不等式约束m1个,等式约束m2个,那么:c、x均为n维列向量,b为m1维列向量,beq为m2维列向量,A为m1n维矩阵,Aeq为m2n维矩阵。
在该函数中,输入参数有c,A,b,Aeq,beq,lb,ub,M和TolXInteger。其中c为目标函数所对应设计变量的系数,A为不等式约束条件方程组构成的系数矩阵,b为不等式约束条件方程组右边的值构成的向量。Aeq为等式约束方程组构成的系数矩阵,beq为等式约束条件方程组右边的值构成的向量。lb和ub为设计变量对应的上界和下界。M为具有整数约束条件限制的设计变量的序号,例如问题中设计变量为x1,x2,,x6,要求x2,x3和x6为整数,则M=[2;3;6];若要求全为整数,则M=1:6,或者M=[1;2;3;4;5;6]。TolXInteger为判定整数的误差限,即若某数x和最邻近整数相差小于该误差限,则认为x即为该整数。
在该函数中,输出参数有x, fval和exitflag。其中x为整数规划问题的最优解向量,fval为整数规划问题的目标函数在最优解向量x处的函数值,exitflag为函数计算终止时的状态指示变量。 例1 求解整数规划问题:
fx1x2max
s.t. 4x12x21
4x12x211 2x1
2
x1,x20,且取整数值
算法:
c=[-1;-1]; A=[-4 2;4 2;0 -2]; b=[-1;11;-1]; lb=[0;0]; M=[1;2]; Tol=1e-8;
%均要求为整数变量 %判断是否整数的误差限
%求解原问题松弛线性规划
[x,fval]=linprog(c,A,b,[],[],lb,[])
[x1,fval1]=intprog(c,A,b,[],[],lb,[],M,Tol) %求最优解整数解 结果:
x =
%松弛线性规划问题的最优解
1.5000 2.5000 fval = -4.0000 x1 = 2 1 fval2 = -3.0000
%整数规划的最优解
(二) 用MATLAB求解0-1规划问题
在MATLAB优化工具箱中,提供了专门用于求解0-1规划问题的函数bintprog,其算法基础即为分枝界定法,在MATLAB中调用bintprog函数求解0-1规划时,需要遵循MATLAB中对0-1规划标准性的要求。 0-1规划问题的MATLAB标准型
mins.t.
fcTxAxbAeqxbeqx0,1
在上述模型中,目标函数f 需要极小化,以及需要满足的约束条件,不等式约束一定要化为形式为“”。
假设x为n维设计变量,且问题具有不等式约束m1个,等式约束m2个,那么:
c、x均为n维列向量,b为m1维列向量,beq为m2维列向量,A为m1n维矩阵,
Aeq为m2n维矩阵。
如果不满足标准型的要求,则需要对原问题进行转化,化为标准型之后才能使用相关函数,标准化的方法和线性规划中的类似。 0-1规划问题的MATLAB求解函数
MATLAB优化工具箱中求解0-1规划问题的命令为bintprog bintprog的调用格式
x = bintprog(f) x = bintprog(f,A,b) x = bintprog(f,A,b,Aeq,beq) x = bintprog(f,A,b,Aeq,beq,x0) x = bintprog(f,A,b,Aeq,Beq,x0,options) [x,fval] = bintprog(...) [x,fval,exitflag] = bintprog(...) [x,fval,exitflag,output] = bintprog(...)
命令详解
1)x = bintprog(f)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
fcTxx0,1
2)x = bintprog(c,A,b)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
3)x = bintprog (c,A,b,Aeq,beq)
fcTxAxbx0,1
该函数调用格式求解如下形式的0-1规划问题:
mins.t.
fcTxAxbAeqxbeqx0,1
4)x = bintprog (c,A,b,Aeq,beq,x0)
该函数调用格式求解如下形式的0-1规划问题
mins.t.
fcTxAxbAeqxbeqx0,1
在前一个调用格式的基础上同时设置求解算法的初始解为x0,如果初始解x0不在0-1规划问题的可行域中,算法将采用默认的初始解 5)x = bintprog (c,A,b,Aeq,beq,x0,options)
用options指定的优化参数进行最小化。可以使用optimset来设置这些参数
上面的函数调用格式仅设置了最优解这一输出参数,如果需要更多的输出参数,则可以参照下面的调用格式:
[x,fval] = bintprog(...)
在优化计算结束之时返回整数规划问题在解x处的目标函数值fval
[x,fval,exitflag] = bintprog(...)
在优化计算结束之时返回exitflag值,描述函数计算的退出条件。
[x,fval,exitflag,output] = bintprog(...) 在优化计算结束之时返回结构变量output 例2:求解0-1规划问题
nn
fEijxij
max
i1j1
20n22s.t. xij1 i1,2,,nEj1
21n
xij1 j1,2,,n22
i1
x0或1( 1,2,3,4) i1,2,,n;j1,2,,n ij
123326
152923 133124
163223
为了程序的可读性,我们用一维下标来表示设计变量,即用x1~x4表示x11~x14,用x5~x8表示x21~x24,用x9~x12表示x31~x34,用x13~x16表示x41~x44,于是约束条件和目标函数分别为:
x1x2x3x41
xxxx1
6785
x9x10x11x121
x13x14x15x161
x1x5x9x131 xxxx1
61014
2
x3x7x11x151
x4x8x12x161xi0,1 (i1,2,,16)
fE11x1E12x2E13x3E14x4E21x5E22x6E44x16
算法:
c=[20;12;33;26;22;15;29;23;21;13;31;24;22;16;32;23]; Aeq=[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 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 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0;
0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0; 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0; 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1; ]; beq=ones(1,8);
[x,fval]=bintprog(c,[],[],Aeq,beq);
B=reshape(x,4,4); %由于x是一列元素,为了使结果更加直观,故排成与效率矩阵E相
对应的形式
B' fval 结果:
Optimization terminated. ans =
0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 fval = 85
整数规划的应用——组件配套问题
某机械产品需要由三个工厂开工一起生产后组装完成。每件机械需要4个组件1和3个组件2。生产这两种组件需要消耗两种原材料A和B。已知这两种原材料的供应量分别为400kg和600kg。
由于三个工厂的生产条件和拥有设备工艺条件不同,每个工厂生产组件的能力和原材料的消耗也不尽相同,且每个工厂开工一次都是配套生产一定数量的组件1和组件2,其具体的数据如表所示。
表11-2 各工厂生产能力和消耗原材料的数据表
现在的最优化问题是,这三个工厂应当如何安排生产,才能使该产品的配套数达到最大?
(Ⅰ)组件配套问题的建模
设x1、x2和x3是三个工厂分别开工的次数,将其作为该问题的设计变量。由于每个工厂开工一次都是配套生产,故每次开工消耗的原材料一定,且生产的组件数也是一定的。在这个假设的前提之下,我们可以得出该问题的目标函数和约束条件。
因为原材料的总量是有限的,根据工厂的开工次数,可得工厂1将消耗A材料9x1,工厂2将消耗A材料6x2,工厂3将消耗A材料4x3,故有约束条件:
9x16x24x3400
同理,对于材料B的总量约束条件可以表达为:7x110x29x3600 然后再来分析三个工厂零件生产的情况,三个工厂生产的组件1的数量分别为8x1、7x2和9x3,故组件1的总数为:Q18x17x29x3
同理,组件2的总数为:Q26x19x25x3
下一步是分析目标函数,该问题要求的不是生产的各种组件总数最多,而是要求产品的配套数量最大,即能组成的机械数目最多。问题中已经给出了该种机械中两种组件的配比为4:3,故组件1所能成套的数目T1满足约束条件:
T1
Q18x17x29x3
44
同理,组件2所能成套的数目T2满足约束条件:T2
Q26x19x25x3
33
因而,所能组成的成品机械的数目f应该为T1和T2中的较小值,即:
fmin(T1,T2)
那么,我们求解的目标即是使得f能够尽可能大,可以看出,这种形式并不是有关设计变量的线性函数,我们需要对其进行转化,此时我们可以令一个人工设计变量等于目标函数值,则有:x4min(T1,T2)
在该假设下,一定满足关系式:T1x4且T2x4 结合约束关系,对T1的约束可以表示为:x4T1
Q18x17x29x3
44
同理,对T2的约束可以表示为:x4T2
Q26x19x25x3
33
对T1的上述关系进行整理,可以得到关系:8x17x29x34x40 同理对T2也可以得到不等式关系为:6x19x25x33x40
上面两个式子即为由组件的配比数得到的关于组件数目的约束条件,此时问题的目标就是如何使得x4取到最大值,由于开工的次数一定是一个非负整数,故也是一个约束条件,决定了该问题是一个纯整数规划问题。结合前面给出的原材料约束,可以得到如下的数学模型:
max fx4
s.t. 9x6x4x400
123
7x110x29x3600
8x17x29x34x40 6x19x25x33x40 xi0且取整数值i1,2,3,4
(Ⅱ)组件配套问题的求解
利用§8.1节中给出的MATLAB函数对此问题求解,代码和运行结果如下:
算法:
%目标函数所对应的设计变量的系数,为转化为极小,故取原目标函数的相反数 f=[0;0;0;-1]; %不等式约束 A=[ 9 6 4 0;
7 10 9 0; -8 -7 -9 4; -6 -9 -5 3]; b=[400;600;0;0];
%边界约束,由于无上界,故设置ub=[Inf;Inf;Inf;Inf]; lb=[0;0;0;0]; ub=[Inf;Inf;Inf;Inf];
%所有变量均为整数变量,故将所有序号组成向量M M=[1;2;3;4];
%判定为整数的误差限 Tol=1e-8;
%求最优解x和目标函数值fval,并返回状态指示 [x,fval,exitflag]=intprog(f,A,b,[],[],lb,ub,M,Tol) 结果:
x = %最优解向量x 18 15 36 141
fval = %在最优解向量x处,原目标函数值的相反数 -141.000
exitflag= %算法终止指示变量,说明问题收敛到了最优解x 1
由运行结果可知,工厂1、2和3需要分别开工18、15和36次,这样所能生产出来的成套的机械为141件。
2 多目标规划的MATLAB求解方法
(一) 多目标规划的MATLAB求解
由于多目标规划中的求解涉及到的方法非常多,故在MATLAB中可以利用
不同的函数进行求解,例如在评价函数法中我们所得最后的评价函数为一线性函数,且约束条件也为线性函数,则我们可以利用MATLAB优化工具箱中提供的linprog函数进行求解,如果我们得到的评价函数为非线性函数,则可以利用MATLAB优化工具箱中提供的fmincon函数进行求解,如果我们采用最大最小法进行求解,则可以利用MATLAB优化工具箱中提供的fminimax函数进行求解。下面我们就结合理论求解的几种方法,讲解一下典型多目标规划问题的MATLAB求解方法。 例1 利用理想点法求解
min f1(x)2x13x2min f(x)5x3x
212
s.t. 3x12x212 xx8
12
x1,x20
我们首先根据评价函数法中的理想点法的理论基础,按照理想点法的求解思路分别对两个单目标规划问题P1,P2进行求解:
min f1(x)2x13x2min f2(x)5x13x2
s.t. 3x2x12s.t. 3x2x121212
P1,P2
xx8 xx81212 x1,x20 x1,x20
求解P1的MATLAB的代码和相应的运行结果为:
算法:
c=[2;-3]; A=[3 2;1 1]; b=[12;8] lb=[0;0]
[x,fval]=linprog(c,A,b,[],[],lb,[]) 结果:
x = 0.0000 6.0000
-18.0000
*
于是可知。当x10,6时,单目标线性规划P 1的最优函数值为f118。
T
求解P2的MATLAB的代码和相应的运行结果为:
算法:
c=[-5;-3]; A=[3 2;1 1]; b=[12;8] lb=[0;0]
[x,fval]=linprog(c,A,b,[],[],lb,[]) 结果:
Optimization terminated. x = 4.0000 0.0000 fval = -20.0000
于是可知,当x24,0时,单目标线性规划P2的最优函数值为f2*20。
T
由上述两个单目标线性规划的求解结果可知x2x2,因而
f
*1
,f2*18,20是一个不可能达到的理想点,因而我们构造如下评价函数:
minhf(x)
xR
f1(x)182f2(x)202
2x13x21825x13x202
构造描述该评价函数的M-函数文件objfun.m如下: function f=objfun(x)
f=sqrt((2*x(1)-3*x(2)+18)^2+(5*x(1)+3*x(2)-20)^2); 然后用非线性规划的方式求解上述问题:
算法:
A=[3 2;1 1]; b=[12;8];
x0=[0;0];
[x,fval,exitflag]=fmincon(@objfun,x0,A,b,[],[],lb,[]) 结果:
Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1
x = 0.0235 5.9647
fval = 1.9941
exitflag = 5
由运行结果可知在该评价函数标准之下,问题的最优解为:
x*0.0235,5.9647此时,各目标函数的取值为:f1*17.8471,f2*18.0118。
T
它与理想点f1*,f2*18,20在评价函数标准下的最小距离为1.9941。 例2 利用评价函数中的线性加权和法求解如下多目标规划问题:
22
min f1(x)x12x2x3222min f2(x)x12x23x3
s.t. x1x2x33 x,x,x0
123
其中权系数为10.8,20.2。 建立线性加权和法的评价函数为:
2222
minhf(x)1x12x2x32x122x23x3
将相应的权系数代入上式即整理出目标函数f(x)为:
22f(x)x121.2x21.4x3
于是建立目标函数的M-函数文件objfun.m:
function f=objfun(x)
f=x(1)^2+1.2*x(2)^2+1.4*x(3)^2;
由于目标函数非线性函数且具有线性等式约束和边界约束,因而我们调用MATLAB中求解非线性规划的fmincon函数对此问题进行求解,同时注意如果只考虑第一个目标函数,由这种特殊形式,即在设计变量的和为一定值,需要求其平方和的最小值时,最优解必然是当这几个设计变量的值相等时取得,于是我们可以将这个解设为问题的初始点,开始迭代。
算法:
Aeq=[1 1 1]; beq=[3]; lb=[0;0;0]; x0=[1;1;1];
[x,fval]=fmincon(@objfun,x0,[],[],Aeq,beq,lb,[]) 结果:
No active inequalities. x = 1.1776 0.9812 0.8412 fval = 3.5327
*x11.1776*,f(x*)3.5327 结果说明,问题的最优解为:x*x20.9812*x30.8412
我们在求解多目标规划问题时,可以采用评价函数法中的最大最小法,而MATLAB也为这种方法提供了专门的求解函数fminimax,在讲解这方面的例题之前,我们首先介绍一下MATLAB优化工具箱中所提供的最大最小法的求解函数fminimax。
最大最小法问题的MATLAB标准形式为:
min max fi(x)
i
x
s.t. c(x)0 c(x)0eq
Axb Axb
eqeq
lbxub
函数fminimax的调用方式和其他的最优化函数类似,其中所涉及的输入参数和输出参数的含义与非线性规划的求解函数fmincon类似,使用方法也基本相同,细节问题读者可以参考MATLAB的帮助文件。 例3 求解最大最小问题:
min max fi(x)i
x
2
s.t. f1(x)2x12x248x140x23042 f2(x)x123x2
f3(x)x13x218
f4(x)x1x2 f5(x)x1x28
首先建立描述目标函数的M-函数文件objfun.m,注意到一共有五个目标函数,所求目标为这五个函数最大值中的最小值,代码如下:
function f = objfun(x)
f(1)= 2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304; f(2)= -x(1)^2 - 3*x(2)^2; f(3)= x(1) + 3*x(2) -18; f(4)= -x(1)- x(2); f(5)= x(1) + x(2) - 8;
然后设置求解的初始点为x0=[0;0],调用fminimax求解该问题。
算法:
x0 = [0; 0];
[x,fval,maxfval] = fminimax(@objfun,x0) 结果:
Local minimum possible. Constraints satisfied.
fminimax stopped because the predicted change in the objective function
is less than the default value of the function tolerance and constraints were satisfied to within the default value of the constraint tolerance.
Active inequalities (to within options.TolCon = 1e-006): lower upper ineqlin ineqnonlin 1 5 x = 4.0000 4.0000 fval =
0.0000 -64.0006 -1.9999 -8.0000 -0.0000 maxfval = 2.6897e-008
上述结果说明当x14,x24时,目标函数fi(x) i1,2,,5的最大值达到最小,这一组的函数值为0.0000,-64.0006,-1.9999,-8.0000,-0.0000,于是最大值为0。
多目标规划的应用——投资问题(全国大学生数学建模竞赛试题)
假设市场上有n种资产,比如股票、债券等可以供投资者选择,某公司有数额为M的一笔相当大的资金可用作一个时间的投资。通过财务人员对Si种资产进行评估,大概可以估算出在这一时期内购买资产Si的平均收益率为ri,并预测出购买Si的损失率为qi。考虑到投资越分散,总的风险越小,公司决定当用这笔资金购买若干种资产时,总体风险可用所投资的Si中的最大一个风险来度量。
购买Si要付交易费,费率为pi,并且当购买额不超过给定值ui时,交易费按购买ui计算(不买当然无须付费)。另外,假定同期银行存款利率是r0,且既无交易费又无风险(r0=5%)。已知n4时的相关数据如下表所示:
表1 投资各种资产的参数值
试给该公司设计一种投资组合方案,即用给定的资金M,有选择地购买若干种资产或存银行生息,使净收益尽可能大,而总体风险尽可能小。 (Ⅰ)投资问题的建模
为了建立数学模型,首先对模型进行一些必要的假设及符号说明: (1) 模型的假设
① 在一个时期内所给出的ri,qi,pi保持不变。
②在一个时间内所购买的各种资产(如股票、证券等)不进行买卖交易,即在买入后不再卖出。
③ 每种投资是否收益是相互独立的。
④ 在投资过程中,无论盈利与否必须先付交易费。 (2)符号说明
M (元):公司现有投资总金额;
Si,i0,1,,n:欲购买的第i种资产种类(其中i0表示存入银行): xi,i0,1,,n:公司购买Si的金额; ri,i0,1,,n:公司购买Si的平均收益率; qi,i0,1,,n:公司购买Si的平均损失率;
pi,i0,1,,n:公司购买Si超过ui时所付交易费率。
下面来建立模型。设购买Si的金额为xi,所付的交易费ci(xi),则
0 xi0
pu 0xu(i1,2,,n)iiii
ci(xi)
pixi xiuic0(x0)0
由于投资额M相当大,所以总可以假定对每个Si的投资xiui。这时上面的大括号公式可简化为:ci(xi)pixi (i1,2,,n)
对Si投资的净收益为:Ri(xi)rixici(xi)ripixi 对Si的风险为:Qi(xi)qixi
对Si投资所需资金为投资金额xi与所需的手续费ci(xi)之和,即:
fi(xi)xici(xi)1pixi
当购买Si的金额为xi (i0,1,2,,n),投资组合x(x0,x1,,xn)的净收益总额为:
R(x)Ri(xi)
i0n
整体风险为:Q(x)maxQi(xi)maxqi(xi)
1in
1in
资金约束为:fi(xi)M
i0
n
根据题目要求,以净收益总额R(x)最大,而整体风险Q(x)最小为目标建立模型如下:
n
min rpxmaxqxiiiiii0
n
s.t. 1pixiM
i0
xi0,i1,2,,n
很显然,这是一个多目标规划问题。 (Ⅱ)投资问题的求解
在此我们采用主要目标法对该问题进行求解,即根据问题的实际情况,确
定一个目标为主要目标,而把其余目标作为次要目标,并且根据经验,选取一定的界限值。这样就可以把次要目标作为约束来处理,于是就将原来的多目标问题转化为一个在新的约束下的单目标最优化问题。
在上述例子中如果以收益为主要目标,则可以固定风险水平,给定风险一个界限a,讲问题转化称为求最大风险不超过a时的最大收益,即下面的线性规划模型:
n
max ripixi
i0
s.t. qixiMa i1,2,,n
n
1pxM
ii
i0
xi0,i1,2,,n
(1)
若投资者希望总盈利至少达到水平K以上,则可以在风险最小的情况下寻找相应的投资组合,从而将原模型转化成为下列的线性规划模型进行求解:
min maxqixiix
n
s.t. ripixiKi0
n
1pixiM
i0
xi0,i1,2,,n
(2)
根据上面的分析,我们利用主要目标法建立了该问题的多目标规划模型,进而转化成为了线性规划模型,下面我们利用MATLAB对此问题进行求解并进行投资分析。将n4时的数据代入模型(1)中,我们可以得到MATLAB的标准型为:
min f0.05x00.27x10.19x20.185x30.185x4
s.t. x1.01x1.02x1.045x1.065x1
01234
0.025x1a
0.015x2a
0.055xa
3
0.026x4a
xi0,i1,2,3,4
由于a是任意给定的风险度,因而实际上没有一个选择的准则,不同的投资者可能有自身不同的判断,因而使得a的取值不同。我们在求解的过程中不妨用试探的方法,从a0开始,以步长a0.001进行搜索,通过实验来分析和总结风险度a和收益Q之间的数量关系。利用MATLAB编程求解。
算法:
clear i=1; a(i)=0; x=cell(1,101); while(1.1-a(i))>1 i=i+1; a(i)=a(i-1)+0.001;
f=[-0.05 -0.27 -0.19 -0.185 -0.185]; Aeq=[1 1.01 1.02 1.045 1.065]; beq=[1]; A=[0 0.025 0 0 0; 0 0 0.015 0 0; 0 0 0 0.055 0; 0 0 0 0 0.026]; b=[a(i);a(i);a(i);a(i)]; lb=[0;0;0;0;0];
[y,fval(i)]=linprog(f,A,b,Aeq,beq,lb,[]); x{i}=y; Q(i)=-fval(i); end plot(a,Q,'.'); xlabel('a'); ylabel('Q'); 结果:
Q
a
风险度与收益关系图
模型(1)的结果图分析:
(1) 由上图可知,收益随风险增大而增大,就是说,风险越大,收益也越大。 (2) 由数据表可以看出,投资越分散,投资者承担的风险越小。这与实际情况相符,就是说,冒险的投资者会出现集中投资的情况,保守的投资者则尽量分散投资。
(3) 曲线上的任一点都表示该风险水平的最大可能收益和该收益要求的最小风险。投资者应根据对不同风险的承受能力,选择该风险水平下的最优投资组合。
(4) 由局部放大上图可以看出,在a(7)0.0060附近有一个转折点。在这一点左边,风险增加很少时,利润增长很快;在这一点右边,风险增加很大时,利润增长很缓慢。所以对于风险和收益没有特殊偏好的投资者来说,应该选择曲线的该转折点作为最优投资组合,大约是a(7)0.0060,Q*Q(7)0.2019,所对应投资方案为:
x[0,0.2400,0.4000,0.1091,0.2212]T