Cplex解决FSP问题

第一次写博客。

流水调度问题（FSP：flow shop scheduling problem）：研究n个工件在m台机器上的流水加工过程，每个工件在各机器上加工顺序相同，每个工件在每台机器上只加工一次，同一时刻一台机器上只能加工一个工件。已知每个工件在每台机器上的加工时间，求调度方案使得某项指标最优。更加详细的概念请看----[1]王凌. 车间调度及其遗传算法[M]. 清华大学出版社, 2003.

此处研究的问题：n个工件，m台机器，每个工件都需要依次经过机器1、机器2、...、机器m的加工才算完成。每个工件在每台机器上的加工时间已知，安排工件的加工顺序，使得最大完工时间（makespan）最小。

Cplex建模：

数据文件(后缀名为.dat)：

job=11;
machine=5;
order=11;
pTimes=[ [375   12   142  245  412],
         [632      452  758  278  398],
         [12      876  124  534  765],
         [460      542  523  120  499],
         [528      101  789  124  999],
         [796      245  632  375  123],
         [532      230  543  896  452],
         [14      124  214  543  785],
         [257      527  753  210  463],
         [896      896  214  258  259],
         [532      302  501  765  988]];

11个工件，5台机器，pTime表示各工件在各机器上的加工时间（横坐标表示工件，纵坐标表示机器）。

模型文件(后缀名为.mod)：

首先从.dat导入数据，并声明变量：

int job=...;
int machine=...;
int order=...;

range j=1..job;//a-->j
range m=1..machine;//b-->m
range k=1..order;//c-->k

float pTimes[j][m]=...;

dvar boolean z[j][k];
dvar float+ startT[m][k];
dvar float+ endT[m][k];

...表示从数据文件中导入，注意命名一致。range表示范围，dvar表示声明的变量。z用来表示工件的加工顺序，startT表示可开工时间，endT表示结束时间。

优化目标：

minimize endT[machine][order];

约束条件：

subject to{
    forall(a in j)
        eachOrderHaveOneOrder://每个订单只有一个加工顺序（每一行的和为1）
        sum(b in k)z[a][b]==1;
    forall(a in j)
        eachOrderIsWhichOrder://每个订单是哪个加工顺序（每一列的和为1）
        sum(b in k)z[b][a]==1;
    forall(c in k)
        forall(b in 1..machine-1)
            nowMachineEndTimeGreatThanPreMach://现在机器的可开工时间大于等于上一机器的完工时间（同一加工阶段）
            endT[b][c]<=startT[b+1][c];
    forall(b in m)
        forall(c in 1..order-1)
            nextProcessStartTimeGreatThanNowProcessEndTime://下一阶段的可开工时间大于等于现在阶段的完工时间（同一机器）
            endT[b][c]<=startT[b][c+1];
    forall(b in m,c in k)
        timeRestriction://完工时间大于等于开工时间加上加工时间
        endT[b][c]>=startT[b][c]+sum(a in j)z[a][c]*pTimes[a][b];
}

总共有5个约束，前两个约束的含义如下：

假设有4个工件，z如下面表格所示。

0	0	1	0
1	0	0	0
0	1	0	0
0	0	0	1

横坐标表示工件，纵坐标表示顺序，必须保证每一行都只有一个为1(1为true，0为false)，每一列都只有一个为1。第一列表示第一个加工的，第二列表示第二个加工的，以此类推，那么上面的z表示的加工顺序为工件2-->3-->1-->4。

后面三个约束有注释，不多赘述。

另外，为了显示出最后的加工顺序，还需要来个后处理：

execute display{
    writeln("工件加工的顺序为：")

    for(var c in k){
        for(var a in j){
            if(z[a][c]==1){
                write(a+"-->");
                break;            
            }        
        }    
    }
}

 

总的模型文件为：

int job=...;
int machine=...;
int order=...;

range j=1..job;//a-->j
range m=1..machine;//b-->m
range k=1..order;//c-->k

float pTimes[j][m]=...;

dvar boolean z[j][k];
dvar float+ startT[m][k];
dvar float+ endT[m][k];

minimize endT[machine][order];
subject to{
    forall(a in j)
          eachOrderHaveOneOrder://每个订单只有一个加工顺序（每一行的和为1）
        sum(b in k)z[a][b]==1;
    forall(a in j)
          eachOrderIsWhichOrder://每个订单是哪个加工顺序（每一列的和为1）
        sum(b in k)z[b][a]==1;
    forall(c in k)
        forall(b in 1..machine-1)
              nowMachineEndTimeGreatThanPreMach://现在机器的可开工时间大于等于上一机器的完工时间（同一加工阶段）
            endT[b][c]<=startT[b+1][c];
    forall(b in m)
        forall(c in 1..order-1)
              nextProcessStartTimeGreatThanNowProcessEndTime://下一阶段的可开工时间大于等于现在阶段的完工时间（同一机器）
            endT[b][c]<=startT[b][c+1];
    forall(b in m,c in k)
          timeRestriction://完工时间大于等于开工时间加上加工时间
        endT[b][c]>=startT[b][c]+sum(a in j)z[a][c]*pTimes[a][b];
}
execute display{
    writeln("工件加工的顺序为：")

    for(var c in k){
        for(var a in j){
            if(z[a][c]==1){
                write(a+"-->");
                break;            
            }        
        }    
    }
}

结果：

// solution (optimal) with objective 7038
// Quality Incumbent solution:
// MILP objective                                7.0380000000e+003
// MILP solution norm |x| (Total, Max)           3.61200e+005 7.03800e+003
// MILP solution error (Ax=b) (Total, Max)       0.00000e+000 0.00000e+000
// MILP x bound error (Total, Max)               0.00000e+000 0.00000e+000
// MILP x integrality error (Total, Max)         0.00000e+000 0.00000e+000
// MILP slack bound error (Total, Max)           0.00000e+000 0.00000e+000
// 

endT = [[14
             389 917 929 1838 2095 2627 3087 3719 4515 5411]
             [138 753 1018 1894 2196 2723 3177 3719 4171 5347 6307]
             [352 895 1968 2196 2697 3523 4066 4589 5347 5979 6521]
             [895 1140 2092 3091 3856 4066 4962 5082 6146 6521 6779]
             [1680 2092 3091 3856 4844 5307 5759 6258 6656 6779 7038]];
z = [[0 1 0 0 0 0 0 0 0 0 0]
             [0 0 0 0 0 0 0 0 1 0 0]
             [0 0 0 1 0 0 0 0 0 0 0]
             [0 0 0 0 0 0 0 1 0 0 0]
             [0 0 1 0 0 0 0 0 0 0 0]
             [0 0 0 0 0 0 0 0 0 1 0]
             [0 0 0 0 0 0 1 0 0 0 0]
             [1 0 0 0 0 0 0 0 0 0 0]
             [0 0 0 0 0 1 0 0 0 0 0]
             [0 0 0 0 0 0 0 0 0 0 1]
             [0 0 0 0 1 0 0 0 0 0 0]];
startT = [[0 14 389 917 1306 1838 2095 2627 3087 3719 4515]
             [14 741 917 1018 1894 2196 2947 3177 3719 5102 5411]
             [138 753 1179 2072 2196 2770 3523 4066 4589 5347 6307]
             [352 895 1968 2557 3091 3856 4066 4962 5868 6146 6521]
             [895 1680 2092 3091 3856 4844 5307 5759 6258 6656 6779]];

由上述结果可知，最优解为7038。

工件加工的顺序为：

8-->1-->5-->3-->11-->9-->7-->4-->2-->6-->10-->