【优化算法】大M法在软时间窗的VRP问题中的应用

在CVRPTW问题中需要设定时间窗约束，以最小化等待成本和惩罚成本。设想中的目标函数是设置的：

m.setObjective(
    quicksum(x[i, j, k] * data.distanceMatrix[i, j] for i, j, k in X_set_tplst)
    + penaltyCost * quicksum((s[i, k] - data.dueTime[i]) * (s[i, k] > data.dueTime[i]) for i in data.customerId for k in data.vehicleId)
    + penaltyWait * quicksum((data.readyTime[i] - s[i, k]) * (s[i, k] < data.readyTime[i]) for i in data.customerId for k in data.vehicleId),
    sense=GRB.MINIMIZE)

目标函数由距离成本、提前到达的等待成本、迟到的惩罚成本构成。其中x[i, j, k]是二元变量，代表车辆k从客户点i到客户点j，s[i, k]为时间戳辅助变量（类似于MTZ辅助变量u）。

该模型中软时间窗可以使用一下方法设置约束：

# 时间窗约束（1）
m.addConstrs((s[i, k] + data.distanceMatrix[i][j] - M * (1-x[i, j, k]) <= s[j, k] for i, j, k in X_set_tplst),'-')
# 时间窗约束（2）
m.addConstrs((s[i, k] >= data.readyTime[i] for i, k in S_set_tplst), 'LowerBound')  # 保证时间戳不小于时间窗下界
m.addConstrs((s[i, k] <= data.dueTime[i] for i, k in S_set_tplst), 'UpperBound')  # 保证时间戳不大于时间窗上界

在约束中，如果 x[i, j, k] 为 1（表示车辆 k 从节点 i 移动到节点 j），那么 s[j, k] 必须大于等于 s[i, k] + data.distanceMatrix[i][j] 以符合逻辑。 如果 x[i, j, k] 为 0（表示该路径不被使用），则为了避免不等式限制 s[j, k] 的值，需要通过 -M * (1 - x[i, j, k]) 将这一部分“释放”。

但是目标函数中由于s[i, k] 是一个Gurobi的决策变量，它的值在优化过程中会被求解并赋值。s[i, k] > data.dueTime[i] 这部分试图比较一个决策变量与一个浮点数，这将导致 TypeError。所以我们添加二元辅助决策变量y1[i, k]和y2[i, k]来表示是否违法时间窗约束，其中：

• y1[i, k] 用于指示是否违反时间窗，当 s[i, k] 超过 dueTime[i] 时，y1[i, k] 为 1，表示迟到
• y2[i, k] 用于指示是否违反时间窗，当 s[i, k] 小于 readyTime[i] 时，y2[i, k] 为 1，表示早到

此时，软时间窗约束修改如下：

# 时间窗约束（1）
m.addConstrs((s[i, k] + data.distanceMatrix[i][j] - M * (1-x[i, j, k]) <= s[j, k] for i, j, k in X_set_tplst),'-')  # 保证变量s可行性，这里注意将t_ij用c_ij进行了代替
# 时间窗约束（2）
m.addConstrs((s[i, k] >= data.readyTime[i] for i, k in S_set_tplst), 'LowerBound')  # 保证时间戳不小于时间窗的下界
m.addConstrs((s[i, k] <= data.dueTime[i] for i, k in S_set_tplst), 'UpperBound')  # 保证时间戳不大于时间窗的上界
# 时间窗约束（3），用于定义变量y以便将惩罚项添加至目标
for i, k in S_set_tplst:
    m.addConstr(s[i, k] - data.dueTime[i] <= M * y1[i, k])
    m.addConstr(s[i, k] - data.dueTime[i] >= -M * (1 - y1[i, k]))
    m.addConstr(s[i, k] - data.readyTime[i] <= M * y2[i, k])
    m.addConstr(s[i, k] - data.readyTime[i] >= -M * (1 - y2[i, k]))
    m.addConstr(y1[i, k] + y2[i, k] <= 1)

此时的目标函数设置为：

m.setObjective(
    quicksum(x[i, j, k] * data.distanceMatrix[i, j] for i, j, k in X_set_tplst)
    + quicksum(penaltyCost * y1[i, k] for i, k in S_set_tplst)
    + quicksum(penaltyWait * y2[i, k] for i, k in S_set_tplst),
    sense=GRB.MINIMIZE)