Operations Research (2): Optimization Algorithms课程作业代码及结果（调度｜选址）

来源：Coursera课程Operations Research (2): Optimization Algorithms
Operations Research课程是National Taiwan University开设的运筹优化课程，共分成三个系列：Models and Applications；Optimization Algorithms和Theory，其中Part2主要介绍了一些求解算法，比如单纯形法，梯度下降，分支定界和启发式等，本文记录了该课程的结课小作业

1.调度问题

给定15个作业和加工时间，这些作业可以在3台机器上加工，作业之间存在冲突，即某两个作业不能在同一台机器上加工，要求给出作业加工的调度方案来最小化最大完工时间，使用Gurobi求解

首先建立数学模型

import gurobipy as gp
from gurobipy import GRB

jobs = [i for i in range(1, 16)]
machines = [i for i in range(1, 4)]

conflicts = {1:None, 2:(5,8), 
             3:None, 4:None, 
             5:(2,8), 6:(9,),
             7:(10,), 8:(2,5), 
             9:(6,), 10:(7,), 
             11:(15,), 12:None, 
             13:None, 14:None, 15:(11,)}

process = [7, 4, 6, 9, 12, 8, 10, 11, 8, 7, 6, 8, 15, 14, 3]

model = gp.Model()

# 决策变量
decision = model.addVars(jobs, machines, vtype=GRB.BINARY)
downtime = model.addVars(machines, vtype=GRB.CONTINUOUS)
makespan = model.addVar(vtype=GRB.CONTINUOUS)
#  目标函数
model.setObjective(makespan, GRB.MINIMIZE)
# 约束条件
model.addConstrs(makespan >= downtime[m] for m in machines)
model.addConstrs(downtime[m] == gp.quicksum(decision[j, m]*process[j-1] 
                 for j in jobs) for m in machines) 
model.addConstrs(gp.quicksum(decision[j, m] for m in machines) == 1
                 for j in jobs)
for j in jobs:
    if conflicts[j] is not None:
        for c in conflicts[j]:
            for m in machines:
                model.addConstr(decision[j, m] * decision[c, m] == 0) 

model.optimize()

for m in machines:
    print(downtime[m].x)

Optimal solution found (tolerance 1.00e-04)
Best objective 4.300000000000e+01, best bound 4.300000000000e+01, gap 0.0000%

最优目标值为43，与官方给定结果一致

2.选址问题

2.1 模型求解

有n个街区和m台救护车，已知街区之间的距离和人口，要求选定救护车的地址，使得最大加权消防时间最小，每个街区的加权消防时间=与所有救护车最近的距离✖️该街区的人口，比如只有一台救护车安排在街区1，那么对于街区2，加权消防时间=3✖️30=90，目标是最小化所有街区的加权消防时间最大值。

首先建立数学模型


基于上述模型使用python+Gurobi求解

import gurobipy as gp
from gurobipy import GRB
from itertools import combinations

districts = [i for i in range(1, 9)]
populations = [40, 30, 35, 20, 15, 50, 45, 60]
distances = [[0, 3, 4, 6, 8, 9, 8, 10],
             [3, 0, 5, 4, 8, 6, 12, 9],
             [4, 5, 0, 2, 2, 3, 5, 7],
             [6, 4, 2, 0, 3, 2, 5, 4],
             [8, 8, 2, 3, 0, 2, 2, 4],
             [9, 6, 3, 2, 2, 0, 3, 2],
             [8, 12, 5, 5, 2, 3, 0, 2],
             [10, 9, 7, 4, 4, 2, 2, 0]]

ambulance_nb = 3
# edges = [(i, j) for i, j in combinations(districts, 2)]
model = gp.Model()

# 决策变量
decision = model.addVars(districts, vtype=GRB.BINARY)
weight_district = model.addVars(districts, vtype=GRB.CONTINUOUS)
weight_max = model.addVar(vtype=GRB.CONTINUOUS)
close = model.addVars(districts, districts, vtype=GRB.BINARY)

# 目标函数
model.setObjective(weight_max, GRB.MINIMIZE)

# 约束
model.addConstrs(weight_max >= weight_district[i] for i in districts)
model.addConstr(gp.quicksum(decision[i] for i in districts) == ambulance_nb)

# 对街区i,街区集合districts里最近的那个一定有救护车（筛选边）
model.addConstrs(close[i, j] <= decision[j] for i in districts 
                for j in districts)
# 对街区i,街区集合districts里只有一个是最近的（筛选出唯一边）
model.addConstrs(gp.quicksum(close[i, j] for j in districts) == 1 for i in districts)

model.addConstrs(weight_district[i] >= close[i, j]*distances[i-1][j-1]*populations[i-1] 
                                                   for j in districts for i in districts)

model.optimize()

Optimal solution found (tolerance 1.00e-04)
Best objective 1.350000000000e+02, best bound 1.350000000000e+02, gap 0.0000%

在n=8，m=2时，最优目标值为135，与官方给定结果一致

2.2 启发式求解

对于同样问题，用启发式方法求解，启发式规则如下：首先选择第一辆救护车，假设选中街区1，获取一个最大加权消防时间，同样假设选中街区2在获取一个最大加权消防时间，在所有假设条件中，选择最大加权消防时间最小的假设，确定为第一辆救护车选址。在第一轮选择基础上，开启第二轮选择，以同样的方式决定。设置n=8，m=3，要求计算Gurobi求解和启发式求解两种方法的性能差距。
启发式代码如下：

populations = [40, 30, 35, 20, 15, 50, 45, 60]
distances = [[0, 3, 4, 6, 8, 9, 8, 10],
             [3, 0, 5, 4, 8, 6, 12, 9],
             [4, 5, 0, 2, 2, 3, 5, 7],
             [6, 4, 2, 0, 3, 2, 5, 4],
             [8, 8, 2, 3, 0, 2, 2, 4],
             [9, 6, 3, 2, 2, 0, 3, 2],
             [8, 12, 5, 5, 2, 3, 0, 2],
             [10, 9, 7, 4, 4, 2, 2, 0]]
target = 3
district = 8

# 记录
select = []
rest = [i for i in range(district)]
count = 0
final_max_weight = None

def get_mindist(start, select):
    res = 10000
    for end in select:
        res = min(res, distances[start][end] * populations[start])
    return res

# 循环
while count < target:
    flag = None 
    final_max_weight = None
    for d in rest:
        max_weight = 0 
        
        temp_select = select.copy()
        temp_rest = rest.copy()
        temp_select.append(d)
        temp_rest.remove(d)
        
        for a in temp_rest:
            min_dist = get_mindist(a, temp_select)
            if min_dist >= max_weight:
                max_weight = min_dist
        if final_max_weight == None or final_max_weight > max_weight:
            final_max_weight = max_weight
            flag = d
        
    select.append(flag)
    rest.remove(flag)
    count += 1
    print(select, rest, final_max_weight, count)
    print('*****')

print(final_max_weight)

[3] [0, 1, 2, 4, 5, 6, 7] 240 1
*****
[3, 0] [1, 2, 4, 5, 6, 7] 240 2
*****
[3, 0, 7] [1, 2, 4, 5, 6] 100 3
*****
100

在n=8，m=2时，最优目标值为100，与Gurobi结果一致