从建模到使用Ortools代码求解

Ortools简介

Ortools是一个用于运筹优化的开源优软件套件，用于解决车辆路径，流程，整数和线性规划以及约束等世界上棘手问题。

模型构建

某一工厂生产三种产品P1，P2，P3。每种产品需经过无先后顺序的工步，每种产品对应的工步信息如表1所示。
每种工步只能在一种设备上进行，每种设备可进行多种工步。每件产品执行某一工步的时间，以及机器数量如表2所示。
现有一批订单信息，每个订单说明产品的需求数量以及订单的开始和开始和交付时间如表3所示。

表1：产品工步表

表2：机台产能表

表3：订单信息表

问题要求如下：

1.每天每种产品分配的数量必定为整数。
2.产品每道工序间没有前后关系。
3.一种工步只能有一种机型进行加工, 一种机台可以加工多种工步, 具体看机台产能信息表。
4.一件产品的所有工步都只能安排在同一天进行生产
5.如果出现机台数量不足的情况, 允许机台超负荷工作, 也就是负荷120%, 那么其产能也是120%
6. 机台每天的产能负荷尽可能的均匀(优先级最高)
7. 机台尽可能不出现超负荷现象(优先级最高)
8. 尽可能用最少数量的机台数量(优先级次高)
9. 同订单每天分配的数量尽可能均匀(优先级最低)
10.要求使用线性混合整数规划求解
11.要求使用python求解
12.要求程序运算时间小于1秒
13.如果无法找到最优解, 请给出找可行解或优秀解的算法

问题分析
1.每天每种产品分配的数量必定为整数，可通过设产品每天生产量的变量为整数实现。
$x$${}_{ijk}$: 订单$i$的第$j$产品在第$k$天生产的数量，为整数。
2.产品的每道工序间没有前后关系。表明每件产品只需要完成所有工步即可。
根据决策变量设计规则，只要满足每天产能的约束，该约束自动满足。
设备Mc约束：

${\sum }_{i=1}^5$${\sum }_{j=1}^3$$\frac{x_{ijk}}{m_{jc}}$$<=$$p$${}_c$ $\forall c\in 1,2,3,4,5,6,7,8,9,10$

3.一种工步只能有一种机型进行加工，一种机台可以加工多种工步, 具体看机台产能信息表。该约束可通过机台产能约束实现，和要求2相同。

4.一件产品的所有工步都只能安排在同一天进行生产。要求1中要求每天生产的产品数量为整数，通过决策变量设计，该约束可自动满足。

5.如果出现机台数量不足的情况, 允许机台超负荷工作, 也就是负荷120%, 那么其产能也是120%。该要求是一种额外放松约束，可通过构建让机台在不超负荷情况下，构建模型的约束条件，如果模型有可行解，则该约束要求自动满足。如果不存在可行解，可识别瓶颈设备，将该设备的产能扩大为原来的120%。

6.机台每天的产能负荷尽可能的均匀(优先级最高)
该要求希望机台每天产能平滑，为目标函数。可计算完成这些订单，设备每天的平均负荷$burd$${}_c$,依次计算每台设备的$MSE$即
$MSE$${}_c$$=$${\sum }_{k=1}^{30}$${\sum }_{i=1}^4$${\sum }_{j=1}^3$$($$x$${}_{ijk}$$-$$burd$${}_c$$)$${}^2$ $\forall c\in 1,2,3,4,5,6,7,8,9,10$

所以目标函数为：
$Min$${\sum }_{c=1}^{10}$$MSE$${}_c$

7.机台尽可能不要出现超负荷现象。该要求和要求5的解决方式使一样的。

8.尽可能用最少数量的机台数量(优先级次高)。该要求可在实现要求1和2的基础上调整，由于计算结果为订单每种产品的每天生产数量，因而可知，可知该要求和通过调整每天的生产数量自动得到满足。

9.同订单每天分配的数量尽可能均匀(优先级最低)。该要求首先需满足设备产能平滑要求。概要求同样可通过构建类似于平滑设备产能要求的式子构建约束。

综上所述分析，构建的模型为：

$Min$${\sum }_{c=1}^{10}$$MSE$${}_c$

约束条件为：
1.产品交期约束：即在从开始到交期这段时间，某一订单某种产品生产的数量应等于订单数量。
$T$${}_i$:表示订单$i$的交期。
${\sum }_{k=1}^{T_i}$$x$${}_{ijk}$$=$$o$${}_{ij}$ $\forall i\in 1,2,3,4$
$\forall j\in 1,2,3$
2.设备产能约束：
${\sum }_{i=1}^5$${\sum }_{j=1}^3$$\frac{x_{ijk}}{m_{jc}}$$<=$$p$${}_c$ $\forall c\in 1,2,3,4,5,6,7,8,9,10$
$\forall k\in 1,2,3,...,30$

3.整数约束：
$\forall x_{ijk}\in \mathbb{Z}$

Ortool代码实现

from ortools.sat.python import cp_model
import numpy as np
import pandas as pd

def main():

    # 首先是创建模型
    model = cp_model.CpModel()
    var_upper_bound = 100
    #定义变量
    x = {}
    for i in range(0,5,1):
        for j in range(0,3,1):
            for t in range(0,30,1):
                a = 'x%s%s%s'%(i,j,t)
                b = str(a)
                x[i, j, t] = model.NewIntVar(0, var_upper_bound,'%s'%(b))

    a = {}
    b = {}
    for i in range(0,10,1):
        for t in range(0,30,1):
            f = 'a%s%s'%(i,t)
            h = 'b%s%s'%(i,t)
            a[i,t] = model.NewIntVar(0,10000,'%s'%(f))
            b[i,t] = model.NewIntVar(0,10000,'%s'%(h))

    sum1 = 0
    for t in range(0,15,1):
        sum1 += x[0,0,t]
    model.Add(sum1 == 120)

    sum2 = 0
    for t in range(5,20,1):
        sum2 += x[1,1,t]
    model.Add(sum2 == 140)

    sum3 = 0
    for t in range(10,25,1):
        sum3 += x[2,2,t]
    model.Add(sum3 == 130)

    sum4 = 0
    for t in range(5,20,1):
        sum4 += x[3,1,t]
    model.Add(sum4 == 160)

    sum5 = 0
    for t in range(10,30,1):
        sum5 += x[4,2,t]
    model.Add(sum5 == 200)

    for t in range(0,5,1):
        model.Add(x[0,0,t] <= 15)
    #
    for t in range(5,10,1):
        model.Add(4 * x[0,0,t] + 5 * x[1,1,t] + 5 * x[3,1,t] <= 180)

    for t in range(5,10,1):
        model.Add(x[1,1,t] + x[3,1,t] <= 40)

    for t in range(5,10,1):
        model.Add(x[0,0,t] <= 15)

    for t in range(5,10,1):
        model.Add(2 * x[0,0,t] + x[1,1,t] + x[3,1,t] <= 60)

    for t in range(5,10,1):
        model.Add(2 * x[0,0,t] + 3 * x[1,1,t] + 3 * x[3,1,t] <= 150)

    for t in range(5,10,1):
        model.Add(3 * x[0,0,t] + 2 * x[1,1,t] + 2 * x[3,1,t] <= 90)

    for t in range(10,15,1):
        model.Add(4 * x[0,0,t] + 5 * x[1,1,t] + 3 * x[2,2,t] + 5 * x[3,1,t] + 3 * x[4,2,t] <= 180)

    for t in range(10,15,1):
        model.Add(x[1,1,t] + 2 * x[2,2,t] + x[3,1,t] + 2 * x[4,2,t] <= 40)

    for t in range(10,15,1):
        model.Add(x[0,0,t] <= 15)

    for t in range(10,15,1):
        model.Add(2 * x[0,0,t] + x[1,1,t] + 2 * x[2,2,t] + x[3,1,t] + 2 * x[4,2,t] <= 60)

    for t in range(10,15,1):
        model.Add(4 * x[0,0,t] + 6 * x[1,1,t] + 3 * x[2,2,t] +6 * x[3,1,t] + 3 * x[4,2,t] <= 300)

    for t in range(10,15,1):
        model.Add(2 * x[1,1,t] + x[2,2,t] +2 * x[3,1,t] + x[4,2,t] <= 100)

    for t in range(10,15,1):
        model.Add(3 * x[0,0,t] + 2 * x[1,1,t] + 2 * x[3,1,t] <= 90)

    for t in range(10,15,1):
        model.Add(3 * x[0,0,t] + 2 * x[2,2,t] + 2 * x[4,2,t] <= 60)

    for t in range(10,15,1):
        model.Add(3 * x[0,0,t] + 2 * x[1,1,t] + 3 * x[2,2,t] + 2 * x[3,1,t] + 3 * x[4,2,t] <= 90)

    for t in range(15,20,1):
        model.Add(5 * x[1,1,t] + 3 * x[2,2,t] + 5 * x[3,1,t] + 3 * x[4,2,t] <= 180)

    for t in range(15,20,1):
        model.Add(x[1,1,t] + 2 * x[2,2,t] + x[3,1,t] + 2 * x[4,2,t] <= 48)

    for t in range(15, 20,1):
        model.Add(2 * x[1,1,t] + x[2,2,t] + 2 * x[3,1,t] + x[4,2,t] <= 100)

    for t in range(15,20,1):
        model.Add(x[2,2,t] + x[4,2,t] <= 30)

    for t in range(15,20,1):
        model.Add(2 * x[1,1,t] + 3 * x[2,2,t] + 2 * x[3,1,t] + 3 * x[4,2,t] <=90)

    for t in range(20,25,1):
        model.Add(x[2,2,t]+x[4,2,t] <= 24)

    for t in range(25,30,1):
        model.Add(x[4,2,t] <= 24)

    # 定义目标函数
    for t in range(0, 30, 1):
        model.Add(a[0,t] - b[0,t] == 4 * x[0,0,t] + 5 * x[1,1,t] + 3 * x[2,2,t] + 5 * x[3,1,t] + 3 * x[4,2,t] - 99)
    for t in range(0,30,1):
        model.Add(a[1,t]-b[1,t] == x[1,1,t] + 2 * x[2,2,t] + x[3,1,t] + 2 * x[4,2,t] - 32 )

    for t in range(0,30,1):
        model.Add(a[2,t] - b[2,t] == 2* x[0,0,t] - 8)

    for t in range(0,30,1):
        model.Add(a[3,t]-b[3,t] == 2*x[0,0,t] + x[1,1,t] + 2*x[2,2,t] + x[3,1,t] + 2*x[4,2,t]-40)

    for t in range(0,30,1):
        model.Add(a[4,t]-b[4,t] == 4*x[0,0,t] + 6*x[1,1,t] + 3*x[2,2,t] + 6*x[3,1,t] + 3*x[4,2,t]-109)

    for t in range(0,30,1):
        model.Add(a[5,t]-b[5,t] == 3 * x[0,0,t] + 2 * x[2,2,t] + 2*x[4,2,t] - 34)

    for t in range(0,30,1):
        model.Add(a[6,t]-b[6,t] == 2* x[1,1,t] + x[2,2,t] + 2*x[3,1,t] + x[4,2,t] - 31)

    for t in range(0,30,1):
        model.Add(a[7,t] - b[7,t] == 3*x[0,0,t] + 2*x[1,1,t] + 2*x[3,1,t]-32)

    for t in range(0,30,1):
        model.Add(a[8,t] - b[8,t] == 3*x[0,0,t] + 2 * x[2,2,t] + 2*x[4,2,t]-34)
    for t in range(0,30,1):
        model.Add(a[9,t] - b[9,t] == 3*x[0,0,t] + 2*x[1,1,t] +3*x[2,2,t]+2*x[3,1,t]+3*x[4,2,t]-65)
    M = 0
    for i in range(10):
        for t in range(30):
            M +=a[i,t]+b[i,t]

    model.Minimize(M)

    # 求解并打印结果
    solver = cp_model.CpSolver()
    status = solver.Solve(model)

    if status == cp_model.OPTIMAL:
        for i in range(0,5,1):
            for j in range(0,3,1):
                for t in range(0,30,1):

                    print(x[i,j,t],solver.Value(x[i, j, t]))
    else:
        print('无可行解')

if __name__ == '__main__':
    main()