本文介绍了Benders分解算法，主要用于解决混合整数规划问题。算法通过将问题拆分为主问题和子问题，逐步优化求解。文中给出了算法原理、流程，并通过Cplex库展示了如何构建主模型、子模型及其对偶模型，以及迭代求解过程。最终，通过示例展示了算法的具体应用。
算法介绍
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,t
solutionMP=new double[varNumMp];
//α1,α2,α3,β1,β2,β3,β4,β5
solutionDS=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技术的学习过程。
以上学习笔记，如有侵犯，请立即联系并删除！
扫一扫
1466
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
