java调用cplex实现经典Benders分解算法求解混合整数规划问题

算法介绍

Benders分解算法是由Jacques F. Benders在1962年首先提出，目的是用于解决混合整数规划问题(MIP问题)，即连续变量与整数变量同时出现的极值问题。随着分解算法的不断演变，不严格要求问题形式的广义Benders分解算法被A.M. Geoffrion提出，实现了对具有Benders分解基本形式的非线性问题求解，也就是说在广义Benders分解算法对子问题的求解方法也不必一定是线性的。

这里仅对经典的Benders分解算法进行简要说明。

经典的Benders分解算法将决策变量分为简单变量和复杂变量两类。其中复杂变量一般是指整数变量，而简单变量一般是指连续变量。
在Benders考虑的一类特殊问题中，先将复杂变量的值固定，从而将问题规约为一个一般的线性规划问题（子问题），然后利用割平面的方式向主问题添加极点约束和极射线约束，以达到主问题缩小可行域的目的。因而，该分解算法的基本思想是将大问题分解为小问题（主问题和子问题），并且通过子问题构造主问题的可行约束，以使主问题的逐渐收敛。实质是一种行生成方法。

算法原理

考虑以下的规划模型：

模型当中，x是连续变量，y是整数变量。与x相比，y则更为复杂。因而当规模较大时，变量y导致模型难以求解。那么将变量y的值进行固定，然后对变量x进行求解，相对而言就要简单许多。

如果给定 y 的值，假定为 g(y)，那么剩余部分的模型就可以写成如下主问题：

该模型当中g(y)是一个关于y的变量，因而主问题是一个关于变量y的整数规划问题。其子问题如下：

子问题模型当中y是一个固定值，因而子问题是一个关于变量x的线性规划问题。

将原问题划分为主问题和子问题以后，两个模型之间又存在什么关联呢？答案就是主问题的g(y)是子问题的目标值，也就是说g(y)表示当y值固定时子问题的最优解。那么通过联合主问题和子问题，就能够得到原问题的解，必定存在以下特征（原问题是最小化）：
（1）如果子问题无界，那么主问题也必定无界，此时原问题也无界，即原问题没有最优解；
（2）如果子问题有界，那么主问题也必定有界，此时原问题也有界，即原问题存在最优解。
那么问题又来了，子问题的范围x与固定值y的有关（即x是关于y的变量值），不同的y值对应了不同的x，如何求解子问题？
可以通过求解子问题的对偶问题来计算g(y)，子问题的对偶问题可以写成：

通过将子问题转换为对偶问题，将其中的y值放置在目标函数当中，那么这样就使得无论y取什么值，都不会影响到问题的可行范围。
如果对偶问题的可行域为空，那么原问题无界或可行域也为空。如果对偶问题不为空，可以枚举对偶问题所有的极点（extreme points）和极射线(extreme rays)。关于使用极点和极射线构造问题的可行域，可以参考博文【DW分解介绍】。
对于一个组合问题来说，极点集合Ω_p和极射线集合Ω_r至少存在一个甚至是无法枚举，但是当问题存在最优解时至少存在一个极点使得问题得到最优解；否则，至少存在一条极射线使得问题无界。
使用极点和极射线将子问题进行重写，可以得到一下新子模型：

进而主问题满足以下模型：

这样就可以通过对偶问题获得可行割集和最优割集，并使用其对主问题进行改善。详细的理论介绍可以参考视频【Benders 分解-李纪柳】。

算法流程

示例模型

原模型：

主模型：

//建立主模型private void buildModelMP() throws IloException {modelMP=new IloCplex();modelMP.setOut(null);varsMP=new IloNumVar[varNumMp];for(int i=0;i<varNumMp;i++){if(i==varNumMp-1) {varsMP[i] = modelMP.numVar(0, Double.MAX_VALUE, IloNumVarType.Float, "t");}else{varsMP[i]=modelMP.numVar(0,1, IloNumVarType.Int,"y["+(i+1)+"]");}}//设置目标函数IloNumExpr obj = modelMP.numExpr();for(int i=0;i<varNumMp;i++){if(i==varNumMp-1) {obj = modelMP.sum(obj, varsMP[i]);}else{obj=modelMP.sum(obj,modelMP.prod(7,varsMP[i]));}}modelMP.addMinimize(obj);}

子模型：

//建立子模型private void buildModelSP(double[] solutionMP) throws IloException {modelSP=new IloCplex();modelSP.setOut(null);varsSP=new IloNumVar[varNumSp];for(int i=0;i<varNumSp;i++){varsSP[i] = modelSP.numVar(0, Double.POSITIVE_INFINITY, IloNumVarType.Float, "x["+(i+1)+"]");}//设置目标函数IloNumExpr obj = modelSP.numExpr();for(int i=0;i<varNumSp;i++){obj = modelSP.sum(obj, varsSP[i]);}modelSP.addMinimize(obj);//添加约束IloNumExpr expr = modelSP.numExpr();expr=modelSP.sum(varsSP[0],modelSP.sum(varsSP[3],varsSP[4]));modelSP.addEq(expr,8);expr = modelSP.numExpr();expr=modelSP.sum(varsSP[1],varsSP[4]);modelSP.addEq(expr,3);expr = modelSP.numExpr();expr=modelSP.sum(varsSP[2],varsSP[3]);modelSP.addEq(expr,5);modelSP.addLe(varsSP[0],8*solutionMP[0]);modelSP.addLe(varsSP[1],8*solutionMP[1]);modelSP.addLe(varsSP[2],8*solutionMP[2]);modelSP.addLe(varsSP[3],8*solutionMP[3]);modelSP.addLe(varsSP[4],8*solutionMP[4]);}

子模型的对偶模型：

//建立子模型的对偶模型private void buildModelDS(double[] solutionMP) throws IloException {modelDS=new IloCplex();modelDS.setOut(null);varsDS=new IloNumVar[varNumDs];for(int i=0;i<varNumDs;i++){if(i<3) {varsDS[i] = modelDS.numVar(Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, IloNumVarType.Float, "α["+(i+1)+"]");}else{varsDS[i]=modelDS.numVar(Double.NEGATIVE_INFINITY,0, IloNumVarType.Float,"β["+(i-3+1)+"]");}//System.out.println(varsDS[i].getLB()+"\t"+varsDS[i].getUB());}//设置目标函数IloNumExpr obj = modelDS.numExpr();obj = modelDS.sum(obj, modelDS.prod(8,varsDS[0]));obj = modelDS.sum(obj, modelDS.prod(3,varsDS[1]));obj = modelDS.sum(obj, modelDS.prod(5,varsDS[2]));if(solutionMP[0]!=0)obj = modelDS.sum(obj, modelDS.prod(8*solutionMP[0],varsDS[3]));if(solutionMP[1]!=0)obj = modelDS.sum(obj, modelDS.prod(3*solutionMP[1],varsDS[4]));if(solutionMP[2]!=0)obj = modelDS.sum(obj, modelDS.prod(5*solutionMP[2],varsDS[5]));if(solutionMP[3]!=0)obj = modelDS.sum(obj, modelDS.prod(5*solutionMP[3],varsDS[6]));if(solutionMP[4]!=0)obj = modelDS.sum(obj, modelDS.prod(3*solutionMP[4],varsDS[7]));modelDS.addMaximize(obj);//添加约束IloNumExpr expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],varsDS[3]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[1],varsDS[4]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[2],varsDS[5]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],modelDS.sum(varsDS[2],varsDS[6]));modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],modelDS.sum(varsDS[1],varsDS[7]));modelDS.addLe(expr,1);}

实现流程：
（1）初始化主模型；
（2）主模型求解；
（3）如果主模型的变量t和子模型的目标值相同，则终止循环；
（4）根据主模型的解方案构造子问题的对偶模型并求解；
（5）若对偶模型无界，调用getRay()函数生成极射线；否则，将模型最优解设置为极点；
（6）根据极射线或者极点生成主模型约束；
（7）将约束添加到主模型并更新主模型，转至（2）。

while(count++<ITER){System.out.println("****************************");System.out.println("主模型：");System.out.println(modelMP);//主模型求解boolean typeMP=solveModelMP();//终止条件if(solutionMP[varNumMp-1]==bestObjDS)break;//主模型无解if(!typeMP){break;}//主模型有解else{//构造子模型的对偶模型buildModelDS(solutionMP);System.out.println("子模型的对偶模型：");System.out.println(modelDS);//设置子模型求解参数modelDS.setParam(IloCplex.Param.Preprocessing.Reduce, 0);//modelDS.setParam(IloCplex.Param.RootAlgorithm, IloCplex.Algorithm.Primal);//对偶模型状态判断boolean typeDS=solveModelDS();//根据状态选择可行割和最优割if ( modelDS.getStatus().equals(IloCplex.Status.Unbounded) ) {//获取可行割向量feasibleCut();}else{//获取最优割向量if(typeDS)opyimalCut(solutionDS);}}}

第一次迭代

第二次迭代

第三次迭代

第四次迭代

循环终止得到问题的最优解：22.其中y=[0,0,0,1,1],x=[0,0,0,5,3]

完整代码

以主模型和子模型对偶模型为基础的代码如下：

import ilog.concert.*;
import ilog.cplex.IloCplex;import java.util.HashMap;public class BendersMLPDemo {//变量数int varNumMp,varNumDs,varNumSp;//主模型对象IloCplex modelMP;//主模型变量IloNumVar[] varsMP;//主模型方案double[] solutionMP;//主模型目标值double bestObjMP;//子模型对偶对象IloCplex modelDS;//子模型对偶变量IloNumVar[] varsDS;//子模型对偶方案double[] solutionDS;//子模型对偶目标值double bestObjDS;//主模型对象IloCplex modelSP;//主模型变量IloNumVar[] varsSP;//记录模型变量HashMap<IloNumVar,IloNumExpr>subObj;//迭代次数int ITER;//构造函数public BendersMLPDemo(){varNumMp=6;varNumDs=8;varNumSp=5;//y1,y2,y3,y4,y5,tsolutionMP=new double[varNumMp];//α1,α2,α3,β1,β2,β3,β4,β5solutionDS=new double[varNumDs];ITER=100;bestObjDS=-1e15;}//建立主模型private void buildModelMP() throws IloException {modelMP=new IloCplex();modelMP.setOut(null);varsMP=new IloNumVar[varNumMp];for(int i=0;i<varNumMp;i++){if(i==varNumMp-1) {varsMP[i] = modelMP.numVar(0, Double.MAX_VALUE, IloNumVarType.Float, "t");}else{varsMP[i]=modelMP.numVar(0,1, IloNumVarType.Int,"y["+(i+1)+"]");}}//设置目标函数IloNumExpr obj = modelMP.numExpr();for(int i=0;i<varNumMp;i++){if(i==varNumMp-1) {obj = modelMP.sum(obj, varsMP[i]);}else{obj=modelMP.sum(obj,modelMP.prod(7,varsMP[i]));}}modelMP.addMinimize(obj);}private boolean solveModelMP()throws IloException{boolean rstType=false;if (modelMP.solve()){bestObjMP=modelMP.getObjValue();System.out.println("主模型目标值："+bestObjMP);System.out.println("主模型变量值：");for(int i=0;i< varNumMp;i++){solutionMP[i]=modelMP.getValue(varsMP[i]);System.out.print(solutionMP[i]+"\t");}System.out.println();rstType=true;}else{System.out.println("模型不可解");}return rstType;}//建立子模型private void buildModelSP(double[] solutionMP) throws IloException {modelSP=new IloCplex();modelSP.setOut(null);varsSP=new IloNumVar[varNumSp];for(int i=0;i<varNumSp;i++){varsSP[i] = modelSP.numVar(0, Double.POSITIVE_INFINITY, IloNumVarType.Float, "x["+(i+1)+"]");}//设置目标函数IloNumExpr obj = modelSP.numExpr();for(int i=0;i<varNumSp;i++){obj = modelSP.sum(obj, varsSP[i]);}modelSP.addMinimize(obj);//添加约束IloNumExpr expr = modelSP.numExpr();expr=modelSP.sum(varsSP[0],modelSP.sum(varsSP[3],varsSP[4]));modelSP.addEq(expr,8);expr = modelSP.numExpr();expr=modelSP.sum(varsSP[1],varsSP[4]);modelSP.addEq(expr,3);expr = modelSP.numExpr();expr=modelSP.sum(varsSP[2],varsSP[3]);modelSP.addEq(expr,5);modelSP.addLe(varsSP[0],8*solutionMP[0]);modelSP.addLe(varsSP[1],8*solutionMP[1]);modelSP.addLe(varsSP[2],8*solutionMP[2]);modelSP.addLe(varsSP[3],8*solutionMP[3]);modelSP.addLe(varsSP[4],8*solutionMP[4]);}private void solveModelSP()throws IloException{if (modelSP.solve()){System.out.println("子模型目标值："+modelSP.getObjValue());System.out.println("子模型变量值：");for(int i=0;i< varNumSp;i++){System.out.print(modelSP.getValue(varsSP[i])+"\t");}System.out.println();}else{System.out.println("模型不可解");}}//建立子模型的对偶模型private void buildModelDS(double[] solutionMP) throws IloException {modelDS=new IloCplex();modelDS.setOut(null);varsDS=new IloNumVar[varNumDs];for(int i=0;i<varNumDs;i++){if(i<3) {varsDS[i] = modelDS.numVar(Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, IloNumVarType.Float, "α["+(i+1)+"]");}else{varsDS[i]=modelDS.numVar(Double.NEGATIVE_INFINITY,0, IloNumVarType.Float,"β["+(i-3+1)+"]");}//System.out.println(varsDS[i].getLB()+"\t"+varsDS[i].getUB());}//设置目标函数IloNumExpr obj = modelDS.numExpr();obj = modelDS.sum(obj, modelDS.prod(8,varsDS[0]));obj = modelDS.sum(obj, modelDS.prod(3,varsDS[1]));obj = modelDS.sum(obj, modelDS.prod(5,varsDS[2]));if(solutionMP[0]!=0)obj = modelDS.sum(obj, modelDS.prod(8*solutionMP[0],varsDS[3]));if(solutionMP[1]!=0)obj = modelDS.sum(obj, modelDS.prod(3*solutionMP[1],varsDS[4]));if(solutionMP[2]!=0)obj = modelDS.sum(obj, modelDS.prod(5*solutionMP[2],varsDS[5]));if(solutionMP[3]!=0)obj = modelDS.sum(obj, modelDS.prod(5*solutionMP[3],varsDS[6]));if(solutionMP[4]!=0)obj = modelDS.sum(obj, modelDS.prod(3*solutionMP[4],varsDS[7]));modelDS.addMaximize(obj);//添加约束IloNumExpr expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],varsDS[3]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[1],varsDS[4]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[2],varsDS[5]);modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],modelDS.sum(varsDS[2],varsDS[6]));modelDS.addLe(expr,1);expr = modelDS.numExpr();expr=modelDS.sum(varsDS[0],modelDS.sum(varsDS[1],varsDS[7]));modelDS.addLe(expr,1);}private boolean solveModelDS()throws IloException{boolean rstType=false;if (modelDS.solve()){bestObjDS=modelDS.getObjValue();System.out.println("子对偶模型目标值："+bestObjDS);System.out.println("子对偶模型变量值：");for(int i=0;i< varNumDs;i++){solutionDS[i]=modelDS.getValue(varsDS[i]);System.out.print(solutionDS[i]+"\t");}System.out.println();rstType=true;}else{bestObjDS=-1;System.out.println("模型不可解");}return rstType;}private void recordSubObj()throws IloException{subObj=new HashMap<>(varNumDs);IloNumExpr expr = modelMP.numExpr();expr=modelMP.sum(expr,8);subObj.put(varsDS[0],expr);expr = modelMP.numExpr();expr=modelMP.sum(expr,3);subObj.put(varsDS[1],expr);expr = modelMP.numExpr();expr=modelMP.sum(expr,5);subObj.put(varsDS[2],expr);expr=modelMP.prod(8,varsMP[0]);subObj.put(varsDS[3],expr);expr=modelMP.prod(3,varsMP[1]);subObj.put(varsDS[4],expr);expr=modelMP.prod(5,varsMP[2]);subObj.put(varsDS[5],expr);expr=modelMP.prod(5,varsMP[3]);subObj.put(varsDS[6],expr);expr=modelMP.prod(3,varsMP[4]);subObj.put(varsDS[7],expr);}//获取可行割private void feasibleCut()throws IloException{//记录模型变量recordSubObj();//获取极射线IloLinearNumExpr rayExpr = modelDS.getRay();IloLinearNumExprIterator iter = rayExpr.linearIterator();IloNumExpr expr= modelDS.numExpr();while ( iter.hasNext() ){IloNumVar var = iter.nextNumVar();double value=iter.getValue();expr=modelDS.sum(expr,modelDS.prod(value,subObj.get(var)));}//添加可行割进入主模型modelMP.addLe(expr,0);System.out.println("可行割：");System.out.println(expr);}//获取最优割private void opyimalCut(double[] solutionDS)throws IloException{IloNumExpr expr = modelMP.numExpr();expr=modelMP.sum(expr,8*solutionDS[0]);expr=modelMP.sum(expr,3*solutionDS[1]);expr=modelMP.sum(expr,5*solutionDS[2]);expr=modelMP.sum(expr,modelMP.prod(8*solutionDS[3],varsMP[0]));expr=modelMP.sum(expr,modelMP.prod(3*solutionDS[4],varsMP[1]));expr=modelMP.sum(expr,modelMP.prod(5*solutionDS[5],varsMP[2]));expr=modelMP.sum(expr,modelMP.prod(5*solutionDS[6],varsMP[3]));expr=modelMP.sum(expr,modelMP.prod(3*solutionDS[7],varsMP[4]));modelMP.addLe(expr,varsMP[5]);System.out.println("最优割：");System.out.println(expr);}//benders分解方法private void bendersMethod()throws IloException{//建立模型buildModelMP();//记录迭代次数int count=1;//循环迭代while(count++<ITER){System.out.println("****************************");System.out.println("主模型：");System.out.println(modelMP);//主模型求解boolean typeMP=solveModelMP();//终止条件if(solutionMP[varNumMp-1]==bestObjDS)break;//主模型无解if(!typeMP){break;}//主模型有解else{//构造子模型的对偶模型buildModelDS(solutionMP);System.out.println("子模型的对偶模型：");System.out.println(modelDS);//设置子模型求解参数modelDS.setParam(IloCplex.Param.Preprocessing.Reduce, 0);//modelDS.setParam(IloCplex.Param.RootAlgorithm, IloCplex.Algorithm.Primal);//对偶模型状态判断boolean typeDS=solveModelDS();//根据状态选择可行割和最优割if ( modelDS.getStatus().equals(IloCplex.Status.Unbounded) ) {//获取可行割向量feasibleCut();}else{//获取最优割向量if(typeDS)opyimalCut(solutionDS);}}}//求解子问题buildModelSP(solutionMP);solveModelSP();}public static void main(String[] args) throws IloException {BendersMLPDemo lp=new BendersMLPDemo();lp.bendersMethod();}
}

以主模型和子模型为基础的代码如下：
如何通过子问题获取对偶问题的极射线（函数DualFarka使用方式）暂不会，以后再考虑补充。

========================================
今天到此为止，后续记录其他cplex技术的学习过程。
以上学习笔记，如有侵犯，请立即联系并删除！