运筹优化在物流行业应用demo案例，Gurobi求解

作者 | xulu1352 目前在一家互联网公司从事推荐算法工作（知乎：xulu1352）

编辑 | auroral-L

（本篇文章已获得作者的独家授权，全文共4650字，预计阅读时间20分钟。）

缘起

最近在清理笔记本硬盘文件，看到一个压箱底文件夹（本人早已脱离靠运筹模型恰饭的行业，进了推荐算法的坑），此文件是当年本人在老东家（某快递）做的一个组内分享gurobi求解数学建模的资料。国内互联网上关于gurobi这种土豪级专业运筹求解器的学习资料并不多见，今天就把当年做的demo案例分享给有需要的人，相关资料已经做了脱敏处理。

问题描述

在航空货运领域，使用货机运输快件。与干线(陆运)最大的不同是，飞机装载货物是通过不同型号的板箱，且飞机有固定的航线。

已知的条件：

1.目前除了全勤的货航飞机（下称全货机）以外，还有一些市场上的客机行李舱（下称散航）可以租用；

2.全货机通常是物流公司的固定资产，因此成本远低于散航；

3.路由，这里指航线，其输入组成：全货机-全货机；全货机-散航；散航直达；全货机直达；

4.全货机的板箱有多种类型：AKE、FC、DQF、PAG等，受飞机内部结构影响做成的不同形状的板箱，其中PAG板是规则的板，可以在飞机转运过程中适用于各种机型，因此PAG板在人工经验下一般装载相同目的地的货物，以便于转运时装卸，从而提高时效；

5.名词解释：od指从始发地到目的地的一票件，为缩小计算规模，通常会按半小时汇总一次；

6.全货机与散航都有机型的区别，载重与装载率限制，货物可以按票计量也可以按重计量；

下述模型解决的问题：如何在货航运输中最有效地利用PAG板，使得整个航空运输网络的时效最大化，时效一致的情况下选择成本更低的航线；

Problem Description

目标函数：优先全网时效最大化，时效相同的情况下选择成本较低的航线

约束条件：

（1）：整板中转路由的板数限制

（2）：混板中转/直发路由的板数限制

（3）：货量百分比<=1

（4）：全货机板数的限制

（5）：散航可获取舱位的限制

（6）：转板最低装载率的限制

（7）：转运的板箱PAG数量相等的限制

Sets:

：OD集合(d)

：全货机集合（c）

：散航集合（b）

：路由集合（r）

：全货机整板中转路由集合（rt,RT⊆CR）

：全货机直发路由/混板中转集合（rs,RS⊆CR）

：满载率系数 0.95

：最低装载率系数 0.5

：惩罚系数

Parameter:

：OD pair d的货量

：OD pair d 的票量

：散航b的容量

：全货机c的容量

：全货机c可承载的PAG板箱个数（针对需要中转的路由）

：全货机c除了PAG外剩余舱位

：全货机PAG板箱的容量，可以容纳的OD重量

：路由r的派送达成率（包括次晨和次日）

：路由r的单公斤成本

：全货机直发路由/混板中转路由rs是否包含全货机c，0-1

：全货机整板中转路由rt是否包含全货机c，0-1

：路由r是否包含此散航b，0-1

Decision Variable:

： [0,1],d货量分配在路由r上的百分比,连续变量

：全货机c在全货机整板中转路由rt中分配的PAG个数,整数变量

：全货机c在全货机直发/混板中转路由rs中分配的PAG个数,整数变量

加载数据

from itertools import groupby
from itertools import combinations
import pandas as pd
import numpy as np
import pickle
# from docplex.mp.model import Model
# from docplex.util.environment import get_environment
import gurobipy#od集合
with open('./data/od', 'rb') as f:od = pickle.load(f)
#路由集合
with open('./data/route', 'rb') as f:route = pickle.load(f)
#散航集合
with open('./data/cargo_set_data', 'rb') as f:cargo_set_data = pickle.load(f)
#全货机集合
with open('./data/cargo_route_data', 'rb') as f:cargo_route_data = pickle.load(f)
#全货机路由映射表
with open('./data/cargo_trans_data', 'rb') as f:cargo_trans_data = pickle.load(f)
#全货机直发/混板中转路由映射表
with open('./data/cargo_straight_data', 'rb') as f:cargo_straight_data = pickle.load(f)
# 全货机整板中转路由映射表
with open('./data/passenger_map_data', 'rb') as f:passenger_map_data = pickle.load(f)
# 散航映射表
with open('./data/passenger_data', 'rb') as f:passenger_data = pickle.load(f)#全货机直发/混板中转路由与OD
route_data = route[['od_id','route_id']].copy()
DCS_df = pd.merge(cargo_straight_data,route_data,on=['route_id'],how='left')[['od_id','route_id','cargo_straight_id']].copy()
DCS_df.rename(columns={'cargo_straight_id':'move_id'},inplace=True)
#全货机中转路由与OD
DCT_df = route[(route['aircraft_type'] == 'cargo_aircraft')&(route['is_trans'] == '1')][['od_id','route_id','move_id']].copy()
#散航路由与OD
DBS_df = route[(route['aircraft_type'] == 'passenger_aircraft')][['od_id','route_id','move_id']].copy()DCS = list(map(tuple, DCS_df.values))
DCT = list(map(tuple, DCT_df.values))
DBS = list(map(tuple, DBS_df.values))DBS = sorted(DBS, key=lambda x:x[2])
DBS_iter = groupby(DBS, lambda x:x[2])
DBS_dict = dict([(key,list(group)) for key, group in DBS_iter])DCT = sorted(DCT, key=lambda x:x[2])
DCT_iter = groupby(DCT, lambda x:x[2])
DCT_dict = dict([(key,list(group)) for key, group in DCT_iter])DCS = sorted(DCS, key=lambda x:x[2])
DCS_iter = groupby(DCS, lambda x:x[2])
DCS_dict = dict([(key,list(group)) for key, group in DCS_iter])

创建模型及变量

创建模型

model = gurobipy.Model()
#Compute Server job ID: aadf4435-ae7e-4a24-b4d5-917e0a5acdc5
#Capacity available on '*.*.*.*:80' - connecting...
#Established HTTP unencrypted connection

Z变量 ([0,1],d货量分配在路由r上的百分比)

DR = list(map(tuple, route[['od_id', 'route_id']].values))
DR = sorted(DR, key=lambda x:x[0])
DR_iter = groupby(DR, lambda x:x[0])
DR_dict = dict([(key,list(group)) for key, group in DR_iter])Z = model.addVars(DR, lb=0.0, ub=1, vtype=gurobipy.GRB.CONTINUOUS, name="Z")
# Z[('010R#020W#2019-09-17 18:30:00#T4', 29819)]

XT变量（全货机c在全货机整板中转路由rt中分配的PAG个数,整数变量）

#整板中转路由与全货机映射关系
RT_C = cargo_trans_data[['cargo_transportation_id', 'flight_id']].copy()
RT_C.drop_duplicates(inplace=True)
RT_C_MAP = list(map(tuple, RT_C.values)) #全货机中转路由和全货机映射关系
RT_C_MAP = sorted(RT_C_MAP, key=lambda x:x[0])
RT_C_iter = groupby(RT_C_MAP, lambda x:x[0])
RT_C_dict = dict([(key,list(group)) for key, group in RT_C_iter])
RT = RT_C.cargo_transportation_id.unique()XT = model.addVars(RT_C_MAP, lb=0, ub=25, vtype=gurobipy.GRB.INTEGER, name="XT")
model.update()
# XT[('CAN#O36979#WUX&WUX#O36902#CGO', 'CAN#O36979#WUX')]

XS变量（全货机c在全货机直发/混板中转路由rs中分配的PAG个数,整数变量）

RS_C = cargo_straight_data[['cargo_straight_id','flight_id']].copy()
RS_C.drop_duplicates(inplace=True)
RS_C_MAP = list(map(tuple, RS_C.values)) #全货机直发路由和全货机映射关系
# model.XS = model.integer_var_dict(RS_C_MAP,lb=0, ub=25, name=lambda xs: 'XS_{}_{}'.format(xs[1], xs[0]))
XS = model.addVars(RS_C_MAP, lb=0, ub=25, vtype=gurobipy.GRB.INTEGER, name="XS")
model.update()
# XS[('SZX#O36879#CGO', 'SZX#O36879#CGO')]

参数定义

集合

#系数
alpha = 0.95 #满载
beta = 0.5 #最低装载
theta = 100 #时效激励系数
gamma = 1 #成本惩罚系数
omega = 80 #无路由惩罚系数
# Parameter
#路由集合
R = route[['route_id', 'rate','route_cost']].set_index('route_id').to_dict()
# R_id = list(route['route_id'].unique())
Tr = R['rate'] #路由时效
Lr = R['route_cost'] #路由成本
#od集合
D = od[['od_id', 'package_qty','real_weight']].set_index('od_id').to_dict()
D_id = list(od['od_id'].unique())
Nd = D['package_qty'] #od票数
Wd = D['real_weight'] #od重量
#全货机集合
C = cargo_set_data[['flight_id','capacity','PAG_qty','non_PAG_capacity','PAG_capacity']].set_index('flight_id').to_dict()
C_id = list(cargo_set_data['flight_id'].values) #全货机id
Mc = C['capacity'] #全货机业载
Pc = C['PAG_qty'] #全货机c可承载的PAG板箱个数
Qc = C['non_PAG_capacity'] #全货机c除了PAG外剩余舱位
# V = C['PAG_capacity'] #全货机PAG板箱的容量
#散航集合
B = passenger_data[['move_id','capacity']].set_index('move_id').to_dict()
B_id = list(passenger_data['move_id'].values) #散航id
Mb = B['capacity'] #散航可获取舱位# 路有单个PAG容量
V = cargo_route_data[['cargo_route_id','PAG_capacity']].drop_duplicates().set_index('cargo_route_id').to_dict()['PAG_capacity']

约束条件

#整板中转路由 全货机分配的PAG
model.addConstrs(gurobipy.quicksum(Z[(dr[0],dr[1])] * Wd[dr[0]] for dr in DCT_dict[rt_c[0]]) <= alpha * XT[rt_c] * V[rt_c[0]] for rt_c in RT_C_MAP)# 直发路由 全货机分配的货量
model.addConstrs(gurobipy.quicksum(Z[(dr[0],dr[1])] * Wd[dr[0]] for dr in DCS_dict[rs_c[0]]) <= alpha * XS[rs_c] * V[rs_c[0]] + Qc[rs_c[1]] for rs_c in RS_C_MAP)# Z变量约束
model.addConstrs(gurobipy.quicksum(Z[dr] for dr in DR_dict[d]) <= 1 for d in D_id)# 全货机PAG板约束
model.addConstrs(gurobipy.quicksum(XT[rt_c] for rt_c in RT_C_MAP if rt_c[1] == c) + gurobipy.quicksum(XS[rs_c] for rs_c in RS_C_MAP if rs_c[1] == c) <= Pc[c] for c in C_id)#散航可获取舱位约束
model.addConstrs(gurobipy.quicksum(Z[(dr[0],dr[1])] * Wd[dr[0]] for dr in DBS_dict[b])<= alpha * Mb[b] for b in B_id)#转板最低装载率的限制
model.addConstrs(gurobipy.quicksum(Z[(dr[0],dr[1])] * Wd[dr[0]] for dr in DCT_dict[rt_c[0]]) >= alpha * (XT[rt_c] - 1) * V[rt_c[0]] + beta * V[rt_c[0]] for rt_c in RT_C_MAP)# 整板中转路由PAG数相等的限制
for rt in RT:rtc = RT_C_dict[rt]rtc_comb = list(combinations(rtc, 2))for rtc_c in rtc_comb:model.addConstr(XT[rtc_c[0]] == XT[rtc_c[1]])model.update()

目标函数

model.setObjective(gurobipy.quicksum(theta * Z[dr] * Tr[dr[1]] * Nd[dr[0]] - gamma * Lr[dr[1]] * Z[dr] * Wd[dr[0]] for dr in DR) -gurobipy.quicksum(omega * (1 - gurobipy.quicksum(Z[dr] for dr in DR_dict[d])) * Nd[d] for d in D_id),gurobipy.GRB.MAXIMIZE)# model.Params.LogToConsole=True # 显示求解过程
# model.Params.outputFlag = 0
logfile = open('cb.log', 'w')
model._logfile = logfile
# model.Params.MIPGap=0.0000001 # 百分比界差
model.optimize()
logfile.close()

求解日志

Optimize a model with 128960 rows, 605647 columns and 1379343 nonzeros
Variable types: 604159 continuous, 1488 integer (0 binary)
Coefficient statistics:Matrix range     [4e-01, 6e+03]Objective range  [5e-02, 7e+05]Bounds range     [1e+00, 2e+01]RHS range        [1e+00, 1e+04]
Found heuristic solution: objective 6.374111e+07
Presolve removed 86364 rows and 394040 columns (presolve time = 6s) ...
Presolve removed 103872 rows and 467534 columns (presolve time = 12s) ...
Presolve removed 103872 rows and 467534 columns
Presolve time: 11.97s
Presolved: 25088 rows, 138113 columns, 294625 nonzeros
Found heuristic solution: objective 8.552618e+07
Variable types: 137781 continuous, 332 integer (125 binary)Deterministic concurrent LP optimizer: primal and dual simplex
Showing first log only...Root simplex log...Iteration    Objective       Primal Inf.    Dual Inf.      Time0    7.9463817e+06   0.000000e+00   6.256257e+07     16s
Concurrent spin time: 0.04sSolved with dual simplexRoot relaxation: objective 1.855799e+08, 3937 iterations, 3.93 secondsNodes    |    Current Node    |     Objective Bounds      |     WorkExpl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time0     0 1.8558e+08    0  176 8.5526e+07 1.8558e+08   117%     -   21s
H    0     0                    1.848550e+08 1.8558e+08  0.39%     -   22s
H    0     0                    1.849884e+08 1.8558e+08  0.32%     -   27s0     0 1.8543e+08    0  123 1.8499e+08 1.8543e+08  0.24%     -   29s
H    0     0                    1.853229e+08 1.8543e+08  0.06%     -   31s0     0 1.8542e+08    0  111 1.8532e+08 1.8542e+08  0.05%     -   32s0     0 1.8542e+08    0  111 1.8532e+08 1.8542e+08  0.05%     -   33s0     0 1.8541e+08    0   69 1.8532e+08 1.8541e+08  0.04%     -   35s
H    0     0                    1.853492e+08 1.8541e+08  0.03%     -   36s0     0 1.8540e+08    0   69 1.8535e+08 1.8540e+08  0.03%     -   39s0     0 1.8540e+08    0   69 1.8535e+08 1.8540e+08  0.03%     -   39s0     0 1.8540e+08    0   69 1.8535e+08 1.8540e+08  0.03%     -   40s0     0 1.8540e+08    0   52 1.8535e+08 1.8540e+08  0.03%     -   41s
H    0     0                    1.853721e+08 1.8540e+08  0.02%     -   43s0     0 1.8540e+08    0   49 1.8537e+08 1.8540e+08  0.02%     -   44s0     0 1.8540e+08    0   49 1.8537e+08 1.8540e+08  0.02%     -   44s0     0 1.8540e+08    0   50 1.8537e+08 1.8540e+08  0.02%     -   46s
H    0     0                    1.853783e+08 1.8540e+08  0.01%     -   47s0     0 1.8540e+08    0   48 1.8538e+08 1.8540e+08  0.01%     -   49s0     0 1.8540e+08    0   48 1.8538e+08 1.8540e+08  0.01%     -   50s0     0 1.8540e+08    0   48 1.8538e+08 1.8540e+08  0.01%     -   51s0     0 1.8540e+08    0   41 1.8538e+08 1.8540e+08  0.01%     -   53s
H    0     0                    1.853855e+08 1.8540e+08  0.01%     -   55sCutting planes:Gomory: 60Implied bound: 810Clique: 1MIR: 255Flow cover: 13Explored 1 nodes (7317 simplex iterations) in 55.71 seconds
Thread count was 8 (of 8 available processors)Solution count 9: 1.85385e+08 1.85378e+08 1.85372e+08 ... 8.51746e+07Optimal solution found (tolerance 1.00e-04)
Best objective 1.853854947260e+08, best bound 1.854012555127e+08, gap 0.0085%

解析结果

df_xt = pd.DataFrame([[rt_c[0], rt_c[1], XT[rt_c].X] for rt_c in RT_C_MAP], columns=['cargo_transportation_id', 'flight_id', 'PAG'])
df_xs = pd.DataFrame([[rs_c[0], rs_c[1], XS[rs_c].X] for rs_c in RS_C_MAP], columns=['cargo_straight_id','flight_id', 'PAG'])
df_z = pd.DataFrame([[dr[0], dr[1], Z[dr].X] for dr in DR], columns=['od_id', 'route_id', 'rate'])