作者 | 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. Time
-
0
7.9463817e+06 0.000000e+00 6.256257e+07 16s
-
Concurrent
spin time: 0.04s
-
-
-
Solved
with dual simplex
-
-
-
Root
relaxation: objective 1.855799e+08, 3937 iterations, 3.93 seconds
-
-
-
Nodes
| Current Node | Objective Bounds | Work
-
Expl
Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
-
-
-
0
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% - 27s
-
0
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% - 31s
-
0
0 1.8542e+08 0 111 1.8532e+08 1.8542e+08 0.05% - 32s
-
0
0 1.8542e+08 0 111 1.8532e+08 1.8542e+08 0.05% - 33s
-
0
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% - 36s
-
0
0 1.8540e+08 0 69 1.8535e+08 1.8540e+08 0.03% - 39s
-
0
0 1.8540e+08 0 69 1.8535e+08 1.8540e+08 0.03% - 39s
-
0
0 1.8540e+08 0 69 1.8535e+08 1.8540e+08 0.03% - 40s
-
0
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% - 43s
-
0
0 1.8540e+08 0 49 1.8537e+08 1.8540e+08 0.02% - 44s
-
0
0 1.8540e+08 0 49 1.8537e+08 1.8540e+08 0.02% - 44s
-
0
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% - 47s
-
0
0 1.8540e+08 0 48 1.8538e+08 1.8540e+08 0.01% - 49s
-
0
0 1.8540e+08 0 48 1.8538e+08 1.8540e+08 0.01% - 50s
-
0
0 1.8540e+08 0 48 1.8538e+08 1.8540e+08 0.01% - 51s
-
0
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% - 55s
-
-
-
Cutting
planes:
-
Gomory:
60
-
Implied
bound: 810
-
Clique:
1
-
MIR:
255
-
Flow
cover: 13
-
-
-
Explored
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+07
-
-
-
Optimal
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'])
转载：https://blog.csdn.net/qq_43045873/article/details/114312658
查看评论