C++用户的Cplex使用指南（四）——Cutting stock problem与Column Generation

2022-06-13 23:04:44

1 Cutting stock problem

1.1 模型改进

符号定义：
x j x_j xj​：第 j j j种切法使用的次数（即按这种切法切割的钢卷的数量）；
a i j a_{ij} aij​：在第 j j j种切法中切割出长度为 w i w_i wi​的钢卷的数目。

例如，标准钢卷长度为 W = 100 W=100 W=100,需要的钢卷长度 w i = 25 , 35 , 45 w_i=25,35,45 wi​=25,35,45，每种长度钢卷的需求量 n i = 100 , 200 , 300 , ( i = 1 , … , 3 ) n_i=100,200,300,(i=1,\dots,3) ni​=100,200,300,(i=1,…,3)。
则可以有以下几种切法：
切法1:将1个标准长度的钢卷切成4个 w 1 = 25 w_1=25 w1​=25的，则 a 11 = 4 , a 21 = 0 , a 31 = 0 a_{11}=4,a_{21}=0,a_{31}=0 a11​=4,a21​=0,a31​=0；
切法2:将1个标准长度的钢卷切成1个 w 1 = 25 w_1=25 w1​=25的和2个 w 2 = 35 w_2=35 w2​=35的，则 a 12 = 1 , a 22 = 2 , a 32 = 0 a_{12}=1,a_{22}=2,a_{32}=0 a12​=1,a22​=2,a32​=0；
切法3:将1个标准长度的钢卷切成2个 w 3 = 45 w_3=45 w3​=45的，则 a 13 = 0 , a 23 = 0 , a 33 = 2 a_{13}=0,a_{23}=0,a_{33}=2 a13​=0,a23​=0,a33​=2。
……

则数学模型（ P 2 P_2 P2​）如下：

min ⁡ ∑ j = 1 n x j s.t. ∑ j = 1 n a i j x j ≥ n i , i = 1 , … , m (需求数量约束) x j ∈ Z + , j = 1 , … , n \min\quad \sum^n_{j=1}x_j\\ \text{s.t.}\quad \sum_{j=1}^na_{ij}x_j\ge n_i,i=1,\dots,m\quad\text{(需求数量约束)}\\ x_j\in\mathbb{Z}_+,j=1,\dots,n minj=1∑n​xj​s.t.j=1∑n​aij​xj​≥ni​,i=1,…,m(需求数量约束)xj​∈Z+​,j=1,…,n

其中 n n n为所有切法的数目，同时满足 ∑ i = 1 m w i a i j ≤ W \sum^m_{i=1}w_ia_{ij}\le W ∑i=1m​wi​aij​≤W。

A = [ a 11 a 12 a 13 ⋯ a 21 a 22 a 23 ⋯ ⋮ ⋮ ⋮ ⋮ a m 1 a m 2 a m 3 ⋯ ] A= \begin{bmatrix} a_{11}&a_{12}&a_{13}&\cdots\\ a_{21}&a_{22}&a_{23}&\cdots\\ \vdots&\vdots&\vdots&\vdots\\ a_{m1}&a_{m2}&a_{m3}&\cdots\\ \end{bmatrix} A=⎣⎢⎢⎢⎡​a11​a21​⋮am1​​a12​a22​⋮am2​​a13​a23​⋮am3​​⋯⋯⋮⋯​⎦⎥⎥⎥⎤​

A中的每一列代表一种切法。

1.2 切法数量的估计

C m k ˉ = m ! k ˉ ! ( m − k ˉ ) ! C^{\bar{k}}_m=\frac{m!}{\bar{k}!(m-\bar{k})!} Cmkˉ​=kˉ!(m−kˉ)!m!​

其中， k ˉ \bar{k} kˉ表示每种切法中切割出不同长度钢卷数目的平均值，可以看出当 m m m增大时，切法的数目会呈指数级增大。

虽然依然能用单纯形法求解上述模型 P 1 , P 2 P_1,P_2 P1​,P2​。先不谈单纯形法求解该模型的效率，仅从需要枚举所有的切割方案 x j x_j xj​这一操作，就能知道上述求解方法一定不是最优的。而且上述枚举法存在以下缺点：

  • 当原材料长度增加时，切割方案呈指数式增长，我们无法枚举所有的 x j x_j xj​；

  • 枚举出来的切割方案可能会剩余过多的木料，这种方案通常不会被选择。

单纯形法要求我们列出所有的变量（ x j x_j xj​）。回忆一下，单纯形法在每次迭代的过程中，只能用一个非基变量替换一个基变量，而最终的基变量个数仅与约束条件的数量有关。这么看来，有大部分的 x j x_j xj​是不会被用到的（也就是说这些方案都不够好）。因此，我们不需要，也没必要列出所有的 x j x_j xj​；再者，对于大规模的问题，我们无法列出满足条件的所有变量。基于此，有学者就提出了求解大规模线性规划的高效算法 —— 列生成法，主要适用于**变量数目（列）大于约束条件个数（行）**的线性规划问题。

2 Column Generation

2.1 列生成基本思想

列生成算法的主要思想是对那些无法列出所有变量的原问题（master problem，MP），先考虑有限变量的主问题（restricted master problem ，RMP)，其余变量在有需要的时候，再添加到RMP中。什么变量是需要被加入到RMP中的呢？在单纯形法中，我们根据非基变量检验数的正负来选择进基变量；在列生成中，通过子问题（subproblem）来寻找满足reduced cost条件的变量，如果找到的话，就将它加入到RMP中，直至找不到为止，那么MP就求得了最优解。

2.2 引例

2.2.1 Model

Master Problem

x j x_j xj​：将要使用模式 j j j的次数；

A i j A_{ij} Aij​：满足需求 d i d_i di​所需的每个模式 j j j的项数 i i i。
min ⁡ ∑ j x j s.t. ∑ i , j A i j x j ≥ d i ∀ i ∈ I x j ≥ 0 ∀ j ∈ J \min \sum_j x_j\\ \text{s.t.}\quad \sum_{i,j}A_{ij}x_j\ge d_i\quad\forall i\in I\\ x_j\ge0\quad\forall j\in J minj∑​xj​s.t.i,j∑​Aij​xj​≥di​∀i∈Ixj​≥0∀j∈J

Subproblem

π \pi π：master problem当前解的对偶变量向量， π i \pi_i πi​即master problem的对偶变量；

W W W：卷材宽度， W i W_i Wi​是从卷材上切下的模式中可使用的带材 i i i的宽度；

A i A_i Ai​：该subproblem的建模变量，即在此subproblem中尺寸 i i i的切割次数。
min ⁡ 1 − ∑ i π i A i s.t. ∑ i W i A i ≤ W A i ≥ 0 ∀ i ∈ I \min\quad 1-\sum_i\pi_iA_i\\ \text{s.t.}\quad \sum_iW_iA_i\le W\\ A_i\ge0\quad\forall i\in I min1−i∑​πi​Ai​s.t.i∑​Wi​Ai​≤WAi​≥0∀i∈I
Subproblem的目标函数是新列的Reduced Cost。

2.2.2 Codes

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

/*using name space for cplex program*/
ILOSTLBEGIN

#define RC_EPS 1.0e-6

static void readData(const char* filename, IloNum& rollWidth, IloNumArray& size, IloNumArray& amount){
    ifstream in(filename);
    if(in){
        in >> rollWidth;//input the standard width of the steel roll
        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;
        throw(1);
    }
}

/*report for the master problem*/
static void report1(IloCplex& cutSolver, IloNumVarArray Cut, IloRangeArray Fill){
    cout << endl;
    cout << "Using " << cutSolver.getObjective() << " rolls" << endl;//output the total number of used rolls
    cout << endl;

    /*output the number of each pattern, aka x_j*/
    for (IloInt j = 0; j < Cut.getSize();j++){
        cout << " Cut" << j << " = " << cutSolver.getValue(Cut[j]) << endl;
    }
    cout << endl;
    /*output the shadow price of each pattern*/
    for (IloInt i = 0;i<Fill.getSize();i++){
        cout << " Fill" << i << " = " << cutSolver.getDual(Fill[i]) << endl;
    }
    cout << endl;
}

/*report for the subproblems*/
static void report2(IloAlgorithm& patSolver, IloNumVarArray Use, IloObjective obj){
    cout << endl;
    cout << "Reduced cost is " << patSolver.getValue(obj) << endl;//output the reduced cost of the pattern
    cout << endl;

    /*output the cofficients of each size (aka A_i) in the pattern*/
    if(patSolver.getValue(obj)<=-RC_EPS){
        for (IloInt i = 0; i < Use.getSize(); i++) {
            cout << " Use " << i << " = " << patSolver.getValue(Use[i]) << endl;
        }
        cout << endl;
    }
}

/*report for the final master problem*/
static void report3(IloCplex& cutSolver, IloNumVarArray Cut){
    cout << endl;
    cout << "Best integer solution uses " << cutSolver.getObjValue() << " rolls" << endl;//output the final total number of used rolls
    cout << endl;

    /*output the number of each pattern, aka x_j*/
    for (IloInt j = 0; j < Cut.getSize();j++){
        cout << " Cut" << j << " = " << cutSolver.getValue(Cut[j]) << endl;
    }
}

/*main program*/

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

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

        IloNum rollWidth;//the standard width of steel roll (there is no need to add "(env)" for IloNum/IloInt)
        IloNumArray amount(env);//the wanted amount of cilents
        IloNumArray size(env);

        readData("cutstock.dat", rollWidth, size, amount);//input the data of the problem instance
        //system("pause");
        /*cutting-optimization problem, aka master problem*/
        IloModel cutOpt(env);//declare the model for the master problem

        IloObjective RollsUsed = IloAdd(cutOpt, IloMinimize(env));//add objective to model cutOpt
        IloRangeArray Fill = IloAdd(cutOpt, IloRangeArray(env, amount, IloInfinity));//add constraints to model cutOpt, amount is the lower bound and Ilofinity is the upper bound

        IloNumVarArray Cut(env); //declare decision variables for the master problem, aka x_j; Note that the variables are initialized with numVar

        IloInt nWdth = size.getSize();//the number of different size
        for (j = 0; j < nWdth;j++){
            /*
            add the coefficients of variables by column;
            "RollsUsed(1)"" means the coefficient of the varibale is 1;
            "Fill[j](int(rollWidth / size[j]))" means the coefficient of the variable in the j-th constraint is int(rollWidth / size[j]);
            this program cuts a standard roll with single size as initial solution. 
            */
            Cut.add(IloNumVar(RollsUsed(1) + Fill[j](int(rollWidth / size[j]))));
        }

        IloCplex cutSolver(cutOpt);//declare the solver for the master problem

        /*pattern-generation problem, aka subproblem*/
        IloModel patGen(env);//declare the model for subproblems

        IloObjective ReducedCost = IloAdd(patGen, IloMinimize(env, 1));//"IloMinimize(env, 1)" means the objective function initialize with a constant 1, 1 is the coefficient of the variable in the objective function
        IloNumVarArray Use(env, nWdth, 0.0, IloInfinity, ILOINT);//declare integer decision variables for the subproblems, aka A_i
        patGen.add(IloScalProd(size, Use) <= rollWidth);//(product) add constraints

        IloCplex patSolver(patGen);//declare the solver for subproblems

        /*column-generation procedure*/
        IloNumArray price(env, nWdth);//declare container for the shadow price, aka the cofficients of objective function for subproblems
        IloNumArray newPatt(env, nWdth);//declare the container for new column

        for (;;){
            /*optimize the current pattern, aka solve the current master problem*/
            cutSolver.solve();//solve the master proble
            cutSolver.exportModel("MP.lp");//export the model as lp file
            report1(cutSolver, Cut, Fill);
            /*find and add a new pattern*/
            for (i = 0; i < nWdth;i++){
                price[i] = -cutSolver.getDual(Fill[i]);
            }
            ReducedCost.setLinearCoefs(Use, price);//set the linear coefficients of objective function for subproblems

            patSolver.solve();//solve the subproblem
            patSolver.exportModel("SP.lp");
            report2(patSolver, Use, ReducedCost);

            /*if the redeced cost of the pattern is postive, stop searching for columns*/
            if(patSolver.getValue(ReducedCost)>-RC_EPS){
                break;
            }

            patSolver.getValues(newPatt, Use);//copy the values of "Use" to "newPatt", aka let newPattern be the new column
            Cut.add(IloNumVar(RollsUsed(1) + Fill(newPatt)));//add new column to the master problem
            system("pause");//pause to check out the lp files
        }

        cutOpt.add(IloConversion(env, Cut, ILOINT));//transfer numVar "Cut" into intVar

        cutSolver.solve();//solve the last master problem
        cout << "Solution status: " << cutSolver.getStatus() << endl;//output the solving status
        report3(cutSolver, Cut);
    }
    catch(IloException& ex){
        cerr << "Error: " << ex << endl;
    }
    catch(...){
        cerr << "Error" << endl;
    }

    env.end();
    return 0;
}

2.2.2.1 Analysis

代码来源于IBM CPLEX安装文件夹中的cutstock.cpp示例文件，根据自己的代码习惯和理解做了修改添加了注释。

该问题中，卷材只有一种标准长度，客户需求的长度可以有多种，数量不定。

2.2.2.2 Structure

Meaning code
卷材标准长度 rollWidth
需求数量 amount
需求尺寸 size
Master Problem Subproblem
元素 code 类型 code 类型
模型 cutOpt IloModel patGen IloModel
目标函数 RollsUsed IloObjective ReducedCost IloObjective
约束 Fill IloRangeArray 直接添加
变量 Cut IloNumVarArray
IloIntVarArray		 Use IloNumVarArray
求解算法 cutSolver IloCplex patSolver IloCplex

2.2.3 Problem Instance

rollWidth=115, size=[25, 40, 50, 55, 70], amount=[50, 36, 24, 8, 30].

Iteration 1

MP1
min ⁡ x 1 + x 2 + x 3 + x 4 + x 5 s.t. 4 x 1 > = 50 2 x 2 > = 36 2 x 3 > = 24 2 x 4 > = 8 x 5 > = 30 x i ≥ 0 i = { 1 , … , 5 } \min\quad x_1 + x_2 + x_3 + x_4 + x_5\\ \text{s.t.}\quad 4 x_1 >= 50\\ 2 x_2 >= 36\\ 2 x_3 >= 24\\ 2 x_4 >= 8\\ x_5 >= 30\\ x_i{\color{red}\ge0}\quad i=\{1,\dots,5\} minx1​+x2​+x3​+x4​+x5​s.t.4x1​>=502x2​>=362x3​>=242x4​>=8x5​>=30xi​≥0i={1,…,5}

Cut: [12.5, 18, 12, 4, 30]

dual(Fill): [0.25, 0.5, 0.5, 0.5, 1]

SP1
min ⁡ − 0.25 x 1 − 0.5 x 2 − 0.5 x 3 − 0.5 x 4 − x 5 + 1 s.t. 25 x 1 + 40 x 2 + 50 x 3 + 55 x 4 + 70 x 5 < = 115 x i ≥ 0 i = { 1 , … , 5 } \min\quad - 0.25 x_1 - 0.5 x_2 - 0.5 x_3 - 0.5 x_4 - x_5 + 1\\ \text{s.t.}\quad 25 x_1 + 40 x_2 + 50 x_3 + 55 x_4 + 70 x_5 <= 115\\ x_i\ge0 \quad i=\{1,\dots,5\} min−0.25x1​−0.5x2​−0.5x3​−0.5x4​−x5​+1s.t.25x1​+40x2​+50x3​+55x4​+70x5​<=115xi​≥0i={1,…,5}
Reduced cost = -0.5

Use: [0, 1, 0, 0, 1]

Iteration 2

MP2
min ⁡ x 1 + x 2 + x 3 + x 4 + x 5 + x 6 s.t. 4 x 1 > = 50 2 x 2 + x 6 > = 36 2 x 3 > = 24 2 x 4 > = 8 x 5 + x 6 > = 30 x i ≥ 0 i = { 1 , … , 6 } \min\quad x_1 + x_2 + x_3 + x_4 + x_5 + x_6\\ \text{s.t.}\quad 4 x_1 >= 50\\ 2 x_2 + x_6 >= 36\\ 2 x_3 >= 24\\ 2 x_4 >= 8\\ x_5 + x_6 >= 30\\ x_i\ge0\quad i=\{1,\dots,6\} minx1​+x2​+x3​+x4​+x5​+x6​s.t.4x1​>=502x2​+x6​>=362x3​>=242x4​>=8x5​+x6​>=30xi​≥0i={1,…,6}
Cut: [12.5, 3, 12, 4, 0, 30]

dual(Fill): [0.25, 0.5, 0.5, 0.5, 0.5]

SP2
min ⁡ − 0.25 x 1 − 0.5 x 2 − 0.5 x 3 − 0.5 x 4 − 0.5 x 5 + 1 s.t. 25 x 1 + 40 x 2 + 50 x 3 + 55 x 4 + 70 x 5 < = 115 x i ≥ 0 i = { 1 , … , 5 } \min\quad - 0.25 x_1 - 0.5 x_2 - 0.5 x_3 - 0.5 x_4 - 0.5 x_5 + 1\\ \text{s.t.}\quad 25 x_1 + 40 x_2 + 50 x_3 + 55 x_4 + 70 x_5 <= 115\\ x_i\ge0 \quad i=\{1,\dots,5\} min−0.25x1​−0.5x2​−0.5x3​−0.5x4​−0.5x5​+1s.t.25x1​+40x2​+50x3​+55x4​+70x5​<=115xi​≥0i={1,…,5}
Reduced cost = -0.25

Use: [1, 2, 0, 0, 0]

Iteration 3

MP3
min ⁡ x 1 + x 2 + x 3 + x 4 + x 5 + x 6 + x 7 s.t. 4 x 1 + x 7 > = 50 2 x 2 + x 6 + 2 x 7 > = 36 2 x 3 > = 24 2 x 4 > = 8 x 5 + x 6 > = 30 x i ≥ 0 i = { 1 , … , 7 } \min\quad x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7\\ \text{s.t.}\quad 4 x_1 +x_7 >= 50\\ 2 x_2 + x_6 + 2 x_7 >= 36\\ 2 x_3 >= 24\\ 2 x_4 >= 8\\ x_5 + x_6 >= 30\\ x_i\ge0\quad i=\{1,\dots,7\} minx1​+x2​+x3​+x4​+x5​+x6​+x7​s.t.4x1​+x7​>=502x2​+x6​+2x7​>=362x3​>=242x4​>=8x5​+x6​>=30xi​≥0i={1,…,7}
Cut: [11.75, 0, 12, 4, 0, 30, 3]

dual(Fill): [0.25, 0.375, 0.5, 0.5, 0.625]

SP3
min ⁡ − 0.25 x 1 − 0.375 x 2 − 0.5 x 3 − 0.5 x 4 − 0.625 x 5 + 1 s.t. 25 x 1 + 40 x 2 + 50 x 3 + 55 x 4 + 70 x 5 < = 115 x i ≥ 0 i = { 1 , … , 5 } \min\quad - 0.25 x_1 - 0.375 x_2 - 0.5 x_3 - 0.5 x_4 - 0.625 x_5 + 1\\ \text{s.t.}\quad 25 x_1 + 40 x_2 + 50 x_3 + 55 x_4 + 70 x_5 <= 115\\ x_i\ge0 \quad i=\{1,\dots,5\} min−0.25x1​−0.375x2​−0.5x3​−0.5x4​−0.625x5​+1s.t.25x1​+40x2​+50x3​+55x4​+70x5​<=115xi​≥0i={1,…,5}
Reduced cost = -0.125

Use: [3, 1, 0, 0, 0]

Iteration 4

MP4
min ⁡ x 1 + x 2 + x 3 + x 4 + x 5 + x 6 + x 7 + x 8 s.t. 4 x 1 + x 7 + 3 x 8 > = 50 2 x 2 + x 6 + 2 x 7 + x 8 > = 36 2 x 3 > = 24 2 x 4 > = 8 x 5 + x 6 > = 30 x i ≥ 0 i = { 1 , … , 8 } \min\quad x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7 + x_8\\ \text{s.t.}\quad 4 x_1 +x_7 + 3 x_8 >= 50\\ 2 x_2 + x_6 + 2 x_7 + x_8 >= 36\\ 2 x_3 >= 24\\ 2 x_4 >= 8\\ x_5 + x_6 >= 30\\ x_i\ge0\quad i=\{1,\dots,8\} minx1​+x2​+x3​+x4​+x5​+x6​+x7​+x8​s.t.4x1​+x7​+3x8​>=502x2​+x6​+2x7​+x8​>=362x3​>=242x4​>=8x5​+x6​>=30xi​≥0i={1,…,8}
Cut: [8, 0, 12, 4, 0, 30, 0, 6]

2.3 Multiple standard width

2.3.1 Codes

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

/*using name space for cplex program*/
ILOSTLBEGIN

typedef IloArray<IloModel> IloModelArray; // model array container
typedef IloArray<IloObjective> IloObjectiveArray; // objective array container
typedef IloArray<IloCplex> IloCplexArray; // Cplex array container
typedef IloArray<IloNumVarArray> NumVarMatrix2; // 2-d float variables container

#define RC_EPS 1.0e-6

static void readData(const char* filename, IloNumArray& rollWidth, IloNumArray& rollCost, IloNumArray& size, IloNumArray& amount){
    ifstream in(filename);
    if(in){
        in >> rollWidth;//input the standard width of the steel roll
        in >> rollCost;//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);
    }
}

/*report for the master problem*/
static void report_MP(IloCplex& costSolver, IloNumVarArray Cut, IloRangeArray Fill){
    cout << endl;
    cout << "The total costs of used rolls: " << costSolver.getObjective() << endl; //output the total costs of used rolls
    cout << endl;

    /*output the number of each pattern, aka x_j*/
    for (IloInt j = 0; j < Cut.getSize();j++){
        cout << " Cut" << j << " = " << costSolver.getValue(Cut[j]) << endl;
    }
    cout << endl;
    /*output the shadow price of each pattern*/
    for (IloInt i = 0;i<Fill.getSize();i++){
        cout << " Fill" << i << " = " << costSolver.getDual(Fill[i]) << endl;
    }
    cout << endl;
}

/*report for the subproblems*/
static void report_SP(IloAlgorithm& patSolver, IloNumVarArray Use, IloObjective obj){
    cout << endl;
    cout << "Reduced cost is " << patSolver.getValue(obj) << endl;//output the reduced cost of the pattern
    cout << endl;

    /*output the cofficients of each size (aka A_i) in the pattern*/
    if(patSolver.getValue(obj) >= -RC_EPS){
        for (IloInt i = 0; i < Use.getSize(); i++) {
            cout << " Use " << i << " = " << patSolver.getValue(Use[i]) << endl;
        }
        cout << endl;
    }
}

/*report for the final master problem*/
static void report_FP(IloCplex& cutSolver, IloNumVarArray Cut){
    cout << endl;
    cout << "Best integer solution costs: " << cutSolver.getObjValue() << endl;//output the final total number of used rolls
    cout << endl;

    /*output the number of each pattern, aka x_j*/
    for (IloInt j = 0; j < Cut.getSize();j++){
        cout << " Cut" << j << " = " << cutSolver.getValue(Cut[j]) << endl;
    }
}

/*main program*/
int main(){
    /*---------------------------------Timer Starter---------------------------------*/
    clock_t startTime, endTime;
    startTime = clock();
    IloEnv env;//declare the cplex environment

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

        IloNumArray rollWidth(env);//the standard width of rolls
        IloNumArray rollCost(env);//the standard cost of rolls
        IloNumArray amount(env);//the wanted amount of cilents
        IloNumArray size(env);//the wanted size of cilents

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

        /*cutting-optimization problem, aka master problem*/
        IloModel costOpt(env);//declare the model for the master problem
        IloObjective RollsUsed = IloAdd(costOpt, IloMinimize(env));//add objective to model costOpt
        IloRangeArray Fill = IloAdd(costOpt, IloRangeArray(env, amount, IloInfinity));//add constraints to model costOpt, amount is the lower bound and Ilofinity is the upper bound
        IloNumVarArray Cut(env); //declare decision variables for the master problem, aka x_j; Note that the variables are initialized with numVar
        /*initial cutting pattern*/
        for (j = 0; j < size.getSize();j++){
            /*
            add the coefficients of variables by column;
            "RollsUsed(rollCost[j])"" means the coefficient of the varibale is rollCost[i];
            "Fill[j](int(rollWidth / size[j]))" means the coefficient of the variable in the j-th constraint is int(rollWidth / size[j]);
            this program cuts a standard roll with single size as initial solution. 
            */
            Cut.add(IloNumVar(RollsUsed(rollCost[0]) + Fill[j](int(rollWidth[0] / size[j]))));
        }
        IloCplex costSolver(costOpt);//declare the solver for the master problem

        /*pattern-generation problems, aka subproblems*/
        IloModelArray patGen(env, rollWidth.getSize());
        IloObjectiveArray ReducedCost(env, rollWidth.getSize());
        IloCplexArray patSolver(env, rollWidth.getSize());
        NumVarMatrix2 Use(env, rollWidth.getSize());

        for (i = 0; i < rollWidth.getSize();i++){
            patGen[i] = IloModel(env);//declare the model for subproblems
            ReducedCost[i] = IloAdd(patGen[i], IloMaximize(env, -rollCost[i]));//"IloMinimize(env, -rollCost[i])" means the objective function initialize with a constant -rollCost[i], -rollCost[i] is the coefficient of the variable in the objective function
            Use[i] = IloNumVarArray(env, size.getSize());
            for (j = 0; j < size.getSize();j++){
                Use[i][j]=IloNumVar(env, 0.0, IloInfinity, ILOINT); //declare integer decision variables for the subproblems, aka A_i
            }
            patGen[i].add(IloScalProd(size, Use[i]) <= rollWidth[i]);//(product) add constraints
            patSolver[i] = IloCplex(patGen[i]);//declare the solver for subproblems
        }

        /*column-generation procedure*/
        IloNumArray price(env, size.getSize());//declare container for the shadow price, aka the cofficients of objective function for subproblems
        IloNumArray newPatt(env, size.getSize());//declare the container for new column
        IloBool Z;
        IloInt iter = 0;
        IloInt index;
        IloNum max;

        for (;;){
            Z = 0;
            iter++;
            cout << "Iteration " << iter << "......" << endl;
            /*optimize the current pattern, aka solve the current master problem*/
            cout << "Solving Master Problem..." << endl;
            costSolver.solve();//solve the master proble
            costSolver.exportModel("MP.lp");//export the model as lp file
            report_MP(costSolver, Cut, Fill);

            /*find and add a new pattern*/
            for (i = 0; i < size.getSize(); i++)
            {
                price[i] = costSolver.getDual(Fill[i]);
            }

            for (i = 0; i < rollWidth.getSize();i++){
            ReducedCost[i].setLinearCoefs(Use[i], price);//set the linear coefficients of objective function for subproblems

            cout << "Solving Subproblem " << i + 1 << "..." << endl;
            patSolver[i].solve();//solve the subproblem
            patSolver[i].exportModel("SP.lp");
            report_SP(patSolver[i], Use[i], ReducedCost[i]);
            //system("pause");//pause to check out the lp files
            }

            /*if no redeced cost of the patterns is postive, stop searching for columns*/
            for (i = 0; i < rollWidth.getSize();i++){
                if(patSolver[i].getValue(ReducedCost[i]) >= RC_EPS){
                    Z = 1;
                }
            }

            index = 0;
            max = patSolver[0].getValue(ReducedCost[0]);
            if(rollWidth.getSize() > 1){
            for (i = 1; i < rollWidth.getSize();i++){
                if(patSolver[i].getValue(ReducedCost[i])>max){
                    max = patSolver[i].getValue(ReducedCost[i]);
                    index = i;
                }
            }
            }

            if (!Z){
                break;
            }

            patSolver[index].getValues(newPatt, Use[index]);//copy the values of "Use" to "newPatt", aka let newPattern be the new column
            Cut.add(IloNumVar(RollsUsed(rollCost[index]) + Fill(newPatt)));//add new column to the master problem
            //system("pause");//pause to check out the lp files
        }

        costOpt.add(IloConversion(env, Cut, ILOINT));//transfer numVar "Cut" into intVar
        costSolver.solve();//solve the last master problem
        cout << "Solution status: " << costSolver.getStatus() << endl;//output the solving status
        report_FP(costSolver, Cut);

    }
    catch(IloException& ex){
        cerr << "Error: " << ex << endl;
    }
    catch(...){
        cerr << "Error" << endl;
    }

    env.end();

    /*---------------------------------Timer end---------------------------------*/
    endTime = clock();
    cout << endl;
    cout << "Totle Time : " << (double)(endTime - startTime) / CLOCKS_PER_SEC << "s" << endl;
    
    return 0;
}

2.3.2 Analysis

该问题中，卷材只有一种标准长度，客户需求的长度可以有多种，数量不定。目标函数为所有卷材的成本（卷材单位成本为rollCost，要计算所用卷材数量的话，可将所有卷材成本设为1）。

上一篇：Problem J: IP地址


下一篇：centos/redhat图转字，字转图