1整数规划的MATLAB求解方法(一)用MATLAB求解一般混合整数规划问题由于MATLAB优化工具箱中并未提供求解纯整数规划和混合整数规划的函数,因而需要自行根据需要和设定相关的算法来实现。现在有许多用户发布的工具箱可以解决该类问题。这里我们给出开罗大学的Sherif和Tawfik在MATLABCentral上发布的一个用于求解一般混合整数规划的程序,在此命名为intprog,在原程序的基础上做了简单的修改,将其选择分枝变量的算法由自然序改造成分枝变量选择原则中的一种,即:选择与整数值相差最大的非整数变量首先进行分枝。intprog函数的调用格式如下:[x,fval,exitflag]=intprog(c,A,b,Aeq,beq,lb,ub,M,TolXInteger)该函数解决的整数规划问题为:)取整数(MjxnixubxlbbxAbAxtsxcfjieqeqT),,2,1(0..min在上述标准问题中,假设x为n维设计变量,且问题具有不等式约束1m个,等式约束2m个,那么:c、x均为n维列向量,b为1m维列向量,eqb为2m维列向量,A为nm1维矩阵,eqA为nm2维矩阵。在该函数中,输入参数有c,A,b,Aeq,beq,lb,ub,M和TolXInteger。其中c为目标函数所对应设计变量的系数,A为不等式约束条件方程组构成的系数矩阵,b为不等式约束条件方程组右边的值构成的向量。Aeq为等式约束方程组构成的系数矩阵,beq为等式约束条件方程组右边的值构成的向量。lb和ub为设计变量对应的上界和下界。M为具有整数约束条件限制的设计变量的序号,例如问题中设计变量为621,,,xxx,要求32,xx和6x为整数,则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求解整数规划问题:0,121124124..max212212121,且取整数值xxxxxxxtsxxf算法:c=[-1;-1];A=[-42;42;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.50002.5000fval=-4.0000x1=%整数规划的最优解21fval2=-3.0000(二)用MATLAB求解0-1规划问题在MATLAB优化工具箱中,提供了专门用于求解0-1规划问题的函数bintprog,其算法基础即为分枝界定法,在MATLAB中调用bintprog函数求解0-1规划时,需要遵循MATLAB中对0-1规划标准性的要求。0-1规划问题的MATLAB标准型1,0..minxbxAbAxtsxcfeqeqT在上述模型中,目标函数f需要极小化,以及需要满足的约束条件,不等式约束一定要化为形式为“”。假设x为n维设计变量,且问题具有不等式约束1m个,等式约束2m个,那么:c、x均为n维列向量,b为1m维列向量,eqb为2m维列向量,A为nm1维矩阵,eqA为nm2维矩阵。如果不满足标准型的要求,则需要对原问题进行转化,化为标准型之后才能使用相关函数,标准化的方法和线性规划中的类似。0-1规划问题的MATLAB求解函数MATLAB优化工具箱中求解0-1规划问题的命令为bintprogbintprog的调用格式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规划问题1,0..minxtsxcfT2)x=bintprog(c,A,b)该函数调用格式求解如下形式的0-1规划问题1,0..minxbAxtsxcfT3)x=bintprog(c,A,b,Aeq,beq)该函数调用格式求解如下形式的0-1规划问题:1,0..minxbxAbAxtsxcfeqeqT4)x=bintprog(c,A,b,Aeq,beq,x0)该函数调用格式求解如下形式的0-1规划问题1,0..minxbxAbAxtsxcfeqeqT在前一个调用格式的基础上同时设置求解算法的初始解为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规划问题21;214,3,2110211211..max1111,n,,j,n,,ix,n,,jx,n,,ixtsxEfijniijnjijninjijij),(或23321622243113212329152226331220E为了程序的可读性,我们用一维下标来表示设计变量,即用41~xx表示1411~xx,用85~xx表示2421~xx,用129~xx表示3431~xx,用1613~xx表示4441~xx,于是约束条件和目标函数分别为:)16,,2,1(1,0111111111612841511731410621395116151413121110987654321ixxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxi1644622521414313212111xExExExExExExEf算法:c=[20;12;33;26;22;15;29;23;21;13;31;24;22;16;32;23];Aeq=[1111000000000000;0000111100000000;0000000011110000;0000000000001111;1000100010001000;0100010001000100;0010001000100010;0001000100010001;];beq=ones(1,8);[x,fval]=bintprog(c,[],[],Aeq,beq);B=reshape(x,4,4);%由于x是一列元素,为了使结果更加直观,故排成与效率矩阵E相对应的形式B'fval结果:Optimizationterminated.ans=0100001010000001fval=85整数规划的应用——组件配套问题某机械产品需要由三个工厂开工一起生产后组装完成。每件机械需要4个组件1和3个组件2。生产这两种组件需要消耗两种原材料A和B。已知这两种原材料的供应量分别为400kg和600kg。由于三个工厂的生产条件和拥有设备工艺条件不同,每个工厂生产组件的能力和原材料的消耗也不尽相同,且每个工厂开工一次都是配套生产一定数量的组件1和组件2,其具体的数据如表所示。表11-2各工厂生产能力和消耗原材料的数据表每个工厂消耗原材料的数量(公斤)每个工厂各组件的生产能力(件数)A材料B材料组件1组件2工厂1工厂2工厂39647109879695现在的最优化问题是,这三个工厂应当如何安排生产,才能使该产品的配套数达到最大?(Ⅰ)组件配套问题的建模设21xx、和3x是三个工厂分别开工的次数,将其作为该问题的设计变量。由于每个工厂开工一次都是配套生产,故每次开工消耗的原材料一定,且生产的组件数也是一定的。在这个假设的前提之下,我们可以得出该问题的目标函数和约束条件。因为原材料的总量是有限的,根据工厂的开工次数,可得工厂1将消耗A材料19x,工厂2将消耗A材料26x,工厂3将消耗A材料34x,故有约束条件:400469321xxx同理,对于材料B的总量约束条件可以表达为:6009107321xxx然后再来分析三个工厂零件生产的情况,三个工厂生产的组件1的数量分别为2178xx、和39x,故组件1的总数为:3211978xxxQ同理,组件2的总数为:3212596xxxQ下一步是分析目标函数,该问题要求的不是生产的各种组件总数最多,而是要求产品的配套数量最大,即能组成的机械数目最多。问题中已经给出了该种机械中两种组件的配比为4:3,故组件1所能成套的数目1T满足约束条件:4978432111xxxQT同理,组件2所能成套的数目2T满足约束条件:3596332122xxxQT因而,所能组成的成品机械的数目f应该为1T和2T中的较小值,即:),min(21TTf那么,我们求解的目标即是使得f能够尽可能大,可以看出,这种形式并不是有关设计变量的线性函数,我们需要对其进行转化,此时我们可以令一个人工设计变量等于目标函数值,则有:),min(214TTx在该假设下,一定满足关系式:41xT且42xT结合约束关系,对1T的约束可以表示为:49784321114xxxQTx同理,对2T的约束可以表示为:35963321224xxxQTx对1T的上述关系进行整理,可以得到关系:049784321xxxx同理对2T也可以得到不等式关系为:035964321xxxx上面两个式子即为由组件的配比数得到的关于组件数目的约束条件,此时问题的目标就是如何使得4x取到最大值,由于开工的次数一定是一个非负整数,故也是一个约束条件,决定了该问题是一个纯整数规划问题。结合前面给出的原材料约束,可以得到如下的数学模型:4,3,2,1003596049786009107400469maxi432143213213214ixxxxxxxxxxxxxxxs.t.xf且取整数值(Ⅱ)组件配套问题的求解利用§8.1节中给出的MATLAB函数对此问题求解,代码和运行结果如下:算法:%目标函数所对应的设计变量的系数,为转化为极小,故取原目标函数的相反数f=[0;0;0;-1];%不等式约束A=[9640;71090;-8-7-94;-6-9-53];b=[400;600;0;0];%边界约束,由于无上界,故设置ub=[Inf;Inf;Inf;Inf];lb=[0;0;0;0];ub=[Inf;Inf;Inf;Inf];%所有变量均为整数变量,故将