前言
分支定界（branch and bound）和分支定价（branch and price）仅一字之差，这两者之间也有着紧密联系，简单来说分支定价=分支定界+列生成。个人觉得在运筹学领域，分支定价算法也算是比较高级的算法了，要学习分支定价，首先要对分支定界和列生成非常熟悉，分支定界和列生成的基本原理和Cplex实现可参考我的博客分支定界和列生成，到目前为止，我没有在网上找到完整的分支定价代码，因此我决定自己从头开始写一份，前面的几篇博客都是分支定价的铺垫。
分支定价算法
话不多说，首先给出分支定价算法的基本流程图。这个流程图是我在网上找到的清晰且完整的流程图。简单解释一下，首先我们将找到原问题的一个限制主问题，比如在切割问题中，我们初始只选择少数几种可行的切割方案，然后松弛求解，接下来以列的检验数为目标函数建立子问题，并求解子问题，如果有满足条件的进基列，则添加到主模型，到这里为止，都是列生成的内容。这时候我们需要检查解是否为整数解，如果是非整数解，则需要进行分支定界（例如切割问题的解为20.2，这显然是不合适的），分支定界之后并没有结束，因为我们往限制主模型中添加了新的约束，我们需要重新调用列生成算法寻找进基列，等到无法找到进基列后，再分支定界，分支定界后，再寻找进基列，直到解为整数或者解差于当前界或者解不可行。
切割问题
为了与前面所写的列生成相呼应，本文依然选择切割问题作为本文要求解的问题，但是为了保证求解效果的体现，本文选择规模稍微大一些的切割问题。首先再简单回忆一下切割问题。
本文设计的数据如下，现有长度L=100m的纸卷若干，顾客的需求为2500个3m，2000个6米，1800个7米，1000个11m,要求使用最少的纸卷满足顾客需求。
Cplex实现
至此，完整的分支定价代码已经给出，代码中注释详细，在此不再赘述。
#include<ilcplex/ilocplex.h>
#include<queue>
#include<vector>
//=====================================================
//切割问题描述：现有长度L=100m的纸卷若干，顾客的需求为250个3m,
//2000个6米，1800个7米，1000个11m,要求使用最少的纸卷满足顾客需求
//=====================================================
ILOSTLBEGIN
vector<vector<int>>plan;//记录切割方案
//==================================模型类
class c_model{
public:
double obj;//目标值
int add_size;//添加的列数
bool is_feasible;//解是否可行
vector<int>branch_id;//分支变量
vector<double>lb;//下界
vector<double>ub;//上界
//构造函数
c_model(){
obj=-1e8;
add_size=0;
is_feasible=1;
}
//运算符重载
friend bool operator<(const c_model& m1, const c_model& m2)
{
return m1.obj<m2.obj;
}
void get_obj(IloModel model, IloNumVarArray var);
};
//==================================求解RMP
void c_model::get_obj(IloModel model, IloNumVarArray var)
{
for(int i=0;i<add_size;i++)var[branch_id[i]].setBounds(lb[i],ub[i]);//添加约束
for(int j=4+add_size;j<var.getSize();j++)var[j].setBounds(0,0);//当前模型不涉及的变量取值定为0
//求解当前限制主问题
IloEnv env = model.getEnv();
IloCplex cplex(model);
cplex.setOut(env.getNullStream());
if(!cplex.solve())is_feasible=0;
else obj= cplex.getObjValue();
}
//==================================列生成
bool column_generation(c_model cm,IloModel &model, IloCplex &cplex,
IloRangeArray &con, IloNumVarArray &var, IloObjective &obj)
{
//将cm表达成cplex可求解的模型
for(int i=0;i<cm.add_size;i++)var[cm.branch_id[i]].setBounds(cm.lb[i],cm.ub[i]);//添加约束
for(int j=4+cm.add_size;j<var.getSize();j++)var[j].setBounds(0,0);//当前模型不涉及的变量取值定为0
IloEnv env = model.getEnv();
cplex.setOut(env.getNullStream());
if(!cplex.solve())return 0;//如果该模型无解则直接不再进行列生成
//建立子模型
IloModel s_model(env);
IloNumVarArray a(env);
//子模型变量
a.add(IloNumVar(env, 0, IloInfinity, ILOINT));
a.add(IloNumVar(env, 0, IloInfinity, ILOINT));
a.add(IloNumVar(env, 0, IloInfinity, ILOINT));
a.add(IloNumVar(env, 0, IloInfinity, ILOINT));
//子模型目标函数
IloExpr expr(env);
for(int i=0;i<4;i++)expr+=cplex.getDual(con[i])*a[i];
s_model.add(IloMinimize(env, 1-expr));
//子模型约束
s_model.add(3*a[0]+6*a[1]+7*a[2]+11*a[3]<=100);
//求解子模型
IloCplex s_cplex(s_model);
s_cplex.setOut(env.getNullStream());
s_cplex.solve();
//进基
if(s_cplex.getObjValue()<0){
IloNumArray aa(env);
s_cplex.getValues(a, aa);
var.add(IloNumVar(obj(1)+con[0](aa[0])+con[1](aa[1])+con[2](aa[2])+con[3](aa[3])));//添加列
//记录切割方案
vector<int>p;
for(int j=0;j<aa.getSize();j++)p.push_back(aa[j]);
plan.push_back(p);
}
return 1;
}
int main()
{
priority_queue<c_model>c_models;//模型队列
double bound = 0;//当前最优解
//====================求解松弛模型
IloEnv env;
IloNumArray value(env);//当前最优解对应变量值
IloModel model(env);//最初限制主模型
IloNumVarArray var(env);
IloRangeArray con(env);
IloObjective obj = IloMinimize(env);
con.add(IloRange(env, 2500, IloInfinity));
con.add(IloRange(env, 2000, IloInfinity));
con.add(IloRange(env, 1800, IloInfinity));
con.add(IloRange(env, 1000, IloInfinity));
var.add(IloNumVar(obj(1)+con[0](33)+con[1](0)+con[2](0)+con[3](0)));
var.add(IloNumVar(obj(1)+con[0](0)+con[1](16)+con[2](0)+con[3](0)));
var.add(IloNumVar(obj(1)+con[0](0)+con[1](0)+con[2](14)+con[3](0)));
var.add(IloNumVar(obj(1)+con[0](0)+con[1](0)+con[2](0)+con[3](9)));
vector<int>p1;p1.push_back(33);p1.push_back(0);p1.push_back(0);p1.push_back(0);
vector<int>p2;p2.push_back(0);p2.push_back(16);p2.push_back(0);p2.push_back(0);
vector<int>p3;p3.push_back(0);p3.push_back(0);p3.push_back(14);p3.push_back(0);
vector<int>p4;p4.push_back(0);p4.push_back(0);p4.push_back(0);p4.push_back(9);
plan.push_back(p1);plan.push_back(p2);plan.push_back(p3);plan.push_back(p4);
model.add(obj);
model.add(con);
IloCplex cplex(model);
cplex.setOut(env.getNullStream());
cplex.solve();
//====================分支定价（基本原理与分支定界同，有额外的列生成操作）
c_model cm;
cm.obj=cplex.getObjValue();
cm.add_size=0;
c_models.push(cm);
while(!c_models.empty()){
//无进基列则弹出队顶模型
if(!column_generation(c_models.top(), model,cplex,con,var,obj)){
c_models.pop();
continue;
}
//求解模型
cplex.solve();
IloNumArray val(env);
cplex.getValues(var, val);
//寻找非整数变量
int id = -1;
for (int i = 0; i < var.getSize(); i++) {
if (val[i] != (int)val[i]) {
id = i;
break;
}
}
//若当前解为整数解则更新界
if (id == -1) {
if (cplex.getObjValue() > bound) {
bound = cplex.getObjValue();
value = val;
c_models.pop();
}
else {
c_models.pop();
}
}
//若当前解为非整数解则分两支
else{
c_model m1;//左分支
m1=c_models.top();//继承父节点所有约束
m1.branch_id.push_back(id);//分支变量
m1.lb.push_back(-1e8);//定义下界
int ub = (int)val[id];
m1.ub.push_back(ub);//定义上界
m1.get_obj(model,var);//求解当前模型
m1.add_size++;//记录约束数量
c_models.push(m1);//压入队列
c_model m2;//右分支
m2=c_models.top();//继承父节点所有约束
m2.branch_id.push_back(id);//分支变量
int lb = (int)val[id]+1;
m2.lb.push_back(lb);//定义下界
m2.ub.push_back(1e8);//定义上界
m2.get_obj(model,var);//求解当前模型
m2.add_size++;//记录约束数量
c_models.push(m2);//压入队列
c_models.pop();//弹出队顶模型
}
}
//计算卷纸生产量
int pd[4];for(int i=0;i<4;i++)pd[i]=0;
for(int i=0;i<plan.size();i++){
for(int j=0;j<plan[i].size();j++){
pd[j]+=value[i]*plan[i][j];
}
}
env.out()<<bound<<endl;//输出最优解
for(int i=0;i<4;i++)env.out()<<pd[i]<<"\t";cout<<endl;//输出生产量
env.end();
system("pause");
}
分支定价运行结果如下所示，需要的总纸卷数为439，生产长度为3的卷纸2506，长度为6的卷纸2000，长度为7的卷纸1800，长度为11的卷纸1008。
总结
学会分支定价算法后，就可以可以使用算法求解实际的大型整数规划问题了，如车辆路径规划问题、机组人员调度问题。
423
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
