C++用户的Cplex使用指南（三）——Cutting stock problem求解

1 Cutting Stock Problem

1.1 问题引入及模型建立

以切割下料问题（Cutting Stock Problem）为例介绍Cplex求解优化模型，并引出列生成算法。

假设工厂有标准长度为218cm的钢管，现有客户需要44个长度为81cm的钢管，3个长度为70cm的钢卷，48个长度为68cm的钢卷。请问如何将标准长度为218cm的钢管进行切割，才能保证所使用标准长度钢管的数目最小？

切法1：将1个标准长度的钢管切成1个81cm的钢管；
切法2：将1个标准长度的钢管切成1个70cm的钢管；
切法3：将1个标准长度的钢管切成1个68cm的钢管；
……
切法n:……

上述3种切法虽然浪费材料，但这么切一定能满足要求，所以可以作为求解该问题时的初始解。
还可以有好多种切法，最后会对该问题进行求解。

切割下料问题经典的数学模型如下所示：

假设标准钢卷长度为$W$，$m$个顾客想要$n_i$个长度为$w_i$的钢卷$(i=1,\dots ,m)$，且$(w_i\le W)$。

符号定义：
$K$：可用标准钢卷的集合；
$y_k$：如果标准钢卷$k$被切割，$y_k=1$.否则$y_k=0$；
$x_i^k$：在标准钢卷$k$上将切割出长度为$w_i$的钢卷的数目。

则数学模型（$P_1$）如下：

min ⁡ ∑ k ∈ K y k s.t. ∑ k ∈ K x i k ≥ n i , i = 1 , … , m (需求长度约束) ∑ i = 1 m w i x i k ≤ W y k k ∈ K (长度限制约束) x i k ∈ Z + , y k ∈ { 0 , 1 } \min\quad \sum_{k\in K}y_k\\ \text{s.t.}\quad \sum_{k\in K}x_{i}^{k}\ge n_{i},\quad i=1,\dots,m\quad\text{(需求长度约束)}\\ \sum_{i=1}^{m}w_{i}x_{i}^k\leq W{y_{k}}\quad k\in K\quad\text{(长度限制约束)}\\ x_i^k\in \mathbb{Z}_+,y_k\in\{0,1\} mink∈K∑​yk​s.t.k∈K∑​xik​≥ni​,i=1,…,m(需求长度约束)i=1∑m​wi​xik​≤Wyk​k∈K(长度限制约束)xik​∈Z+​,yk​∈{0,1}

但是这个数学模型从计算角度和理论角度而言效率不高。

主要原因是这个数学模型的线性松弛问题LP很差（即LP问题的解与原问题的解相差很大）。事实上，松弛问题LP的界为$\frac{{\sum }_{i=1}^m{w}_i{n}_i}{W}$:

$$\sum _{k\in K}{y}_k=\sum _{k\in K}\sum _{i=1}^m\frac{w_i{x}_i^k}{W}=\sum _{i=1}^m{w}_i\sum _{k\in K}\frac{x_i^k}{W}=\sum _{i=1}^m\frac{w_i{n}_i}{W}$$

1.2 模型$P_1$的Cplex直接求解

/*cplex header file*/
#include <ilcplex/ilocplex.h>

/*using name space for cplex program*/
ILOSTLBEGIN

typedef IloArray<IloNumVarArray> NumVarMatrix2; // 2-d float variables container

static void readData(const char* filename, IloNum& rollWidth, IloInt& N, IloNumArray& size, IloNumArray& amount){
    ifstream in(filename);
    if(in){
        in >> rollWidth;//input the standard width of the steel roll
        in >> N;//input the costs of rolls
        in >> size;//input the size demand of cilents
        in >> amount;//input the amount demand of cilents
    }
    else{//in case file dose not exist
        cerr << "No such file: " << filename << endl;
        exit(0);
    }
    /*check out the input data*/
    /*
    cout << rollWidth << endl;
    cout << N << endl;
    for (int i = 0; i < size.getSize(); i++) {
        cout << size[i] << "\t";
    }
    cout << endl;
    for (int i = 0; i < amount.getSize(); i++) {
        cout << amount[i] << "\t";
    }
    cout << endl;
    */
}

int main(){
    IloEnv env;//declare the cplex environment

    /*try to establish the model and slove it*/
    try{
        IloInt i, j;

        IloNum rollWidth;//the width of rolls
        IloInt N;//the maximum number of rolls
        IloNumArray amount(env);//the wanted amount of cilents
        IloNumArray size(env);//the wanted size of cilents

        readData("ncutstock.dat", rollWidth, N, size, amount);//input the data of the problem instance
        //system("pause");//check out the data input process

        /*cutting stock problem*/
        IloModel cutOpt(env);//declare the model for cutting stock problem
        /*declare and create variables*/
        IloBoolVarArray y(env, N);//y_i
        NumVarMatrix2 x(env, N);//x_{ij}
        for (i = 0; i < N; i++) {
            y[i] = IloBoolVar(env);
        }
        for (i = 0; i < N; i++) {
            x[i] = IloNumVarArray(env, size.getSize());
            for (j = 0; j < size.getSize(); j++) {
                x[i][j] = IloNumVar(env, 0.0, IloInfinity, ILOINT);
            }
        }

        //construct objective function
        IloExpr obj(env);
        for (i = 0; i < N; i++) {
            obj += y[i];
        }
        //add the objective function to the solving model
        //IloObjective obj = IloMinimize(env, sum);
        cutOpt.add(IloMinimize(env, obj));

        //declare and construct the constraints
        IloRangeArray cons(env);

        IloExprArray con1(env);
        con1 = IloExprArray(env, size.getSize());
        for (j = 0; j < size.getSize(); j++) {
            con1[j] = IloExpr(env);
        }
        for (j = 0; j < size.getSize(); j++) {
            for (i = 0; i < N; i++) {
                con1[j] += x[i][j];
            }
            cons.add(con1[j] >= amount[j]);
        }

        IloExprArray con2(env);
        con2 = IloExprArray(env, N);
        for (i = 0; i < N; i++) {
            con2[i] = IloExpr(env);
        }
        IloExprArray con3(env);
        con3 = IloExprArray(env, N);
        for (i = 0; i < N; i++) {
            con3[i] = IloExpr(env);
        }
        for (i = 0; i < N; i++) {
            for (j = 0; j < size.getSize(); j++) {
                con2[i] += size[j] * x[i][j];
            }
        }
        for (i = 0; i < N; i++) {
            con3[i] = rollWidth * y[i];
        }
        for (i = 0; i < N; i++) {
            cons.add(con3[i] - con2[i] >= 0);
        }

        /*add the constraints to the solving model*/
        cutOpt.add(cons);




        /*sovle the model and refer to the results*/
        IloCplex cutSolver(cutOpt);
        cutSolver.solve();//solve the master proble
        cutSolver.exportModel("cuttingstock.lp");//export the model as lp file
    
        if (cutSolver.getStatus() == IloAlgorithm::Infeasible) {
            env.out() << "No Solution" << endl;
        }
        else {
            cout << endl;
            cout << "The total costs of used rolls: " << cutSolver.getObjValue() << endl; //output the total costs of used rolls
            cout << endl;
        }
    }
    catch(IloException& ex){
        cerr << "Error: " << ex << endl;
    }
    catch(...){
        cerr << "Error" << endl;
    }

    env.end();
    return 0;
}

假设钢管长度为120米，总共有240根，现在有10个顾客，顾客的要求随机生成如下：
demands = [10, 11, 11, 12, 12, 12, 10, 11, 12, 10]
lengths = [92, 59, 97, 32, 38, 55, 80, 75, 108, 57]
求解结果如下：

---

Reference

1. 运筹学（最优化理论）学习笔记 | 列生成法
2. 【学界/编码】从下料问题看整数规划中的列生成方法（附Gurobi求解器源代码）