1 前言

本文以两道经典建模题为例, 进一步介绍 Gurobi 与 Python 的交互, 以及其在建模中的应用. 阅读本文前, 建议读者先配置好 Gurobi 环境, 并且对数学建模有一定的认识 (吹水, 不考虑绝对的严谨性)。

本文也可作为建模小白的“入门指南”, 全文都是按照我的思维过程进行书写, (笔者在建模中承担建模 + 编程 + 模型主框架书写任务, 会将自己的模型以朴素的语言告知队友) 小白可以放心食用, 感受完成一道建模题所需要进行的分析工作.

2 Problem 1 最佳组队安排问题 (出处不详)

2.1 问题描述

在一年一度的我国和美国大学生数学建模竞赛活动中, 任何一个参赛院校都会遇到如何选拔最优秀的队员和科学合理地组队问题, 这是一个最实际的、而且首先需要解决的数学模型问题.

现假设有 20 名队员准备参加竞争, 根据队员的能力和水平要选出 18 名优秀队员分别组成 6 个队, 每个队 3 名队员去参加比赛. 选择队员主要考虑的条件依次为有关学科成绩(平均成绩)、智力水平(反映思维能力、分析问题和解决问题的能力等)、动手能力(计算机的使用和其他方面实际操作能力)、写作能力、外语能力、协作能力(团结协作能力)和其他特长. 每个队员的基本条件量化后如下表所示:

假设所有队员接受了同样的培训, 外部环境相同, 竞赛中不考虑其他的随机因素, 竞赛水平的发挥只取决于表中所给的各项条件, 并且参赛队员都能正常发挥自己的水平. 现在的问题是:

(1) 如何在20名队员中选拔18名优秀队员参加竞赛;

(2) 给出由18名队员组成6个队的组队方案, 使整体竞赛技术水平最高,并给出每个队的竞赛技术水平;

(3) 如果有更多的学生报名参加, 你的模型是否仍然可行?

2.2 问题分析

建模组队问题是经典的团队分工与总效益评估问题. 一般来说, 数学建模竞赛中的队员应该至少包含三项技能: 擅长文字写作与论文排版、灵活运用的算法思想与编程能力、优秀的数学基础与建模能力, 一般的分配就是三位同学各有专长. 这样的分配对于校方而言更有利于针对不同工作的学生进行集中培训, 避免资源浪费; 对于建模团队而言, 每位同学发挥自己的长处, 合作完成一篇建模论文也是大有裨益的.笔者实践经验认为, 优秀的建模团队应该是每位队员都能独立地完成一道建模题目, 具备独立解决问题的基本技能, 在建模工作中相互渗透, 队员各有所长, 这有助于团队的协调分工、减少沟通成本.

但这样的队员可遇不可求, 退而求其次, 能做好自己的本职工作、具有较好的团队合作能力即可.

搜索资料容易查到武汉理工大学2019年美国大学生数学建模竞赛选拔, 可以作为参考. 建模队员的选拔, 往往根据不同的职位, 有不同的能力考察. 比如说编程的同学可能需要考虑算法功底、编程软件使用、动手实践能力、项目经验、智力水平等, 这也启发我们, 对于学生能力的评估, 不能用统一的评价体系一锅端. 我们要解决的第一个问题就是怎么选拔出适合从事单项工作的同学. 选拔出从事单项工作的队员后, 下一步要解决的问题就是如何进行组队.计算方差, 让总体实力均衡? 计算平均数, 让总体能力均值最高? 似乎都可行, 但是联想到高中的时候开设的火箭班 (尖子班), 校方的考量应该是让优秀的同学冲刺更高的奖项 (更好的大学), 中间的学生尽可能拿奖 (不那么差的学校). 因此就涉及到第三个问题, 如何衡量一只队伍获奖的潜力.

基于以上分析, 解决该问题的思路为:

评价各个学生适合分别从事三项工作的得分值 -> 将学生进行工作分队 -> 建立获奖潜力评估模型 -> 组队安排, 让获奖的总效益最大化当然, 注意到题目说的 “任何一个参赛院校都会遇到如何选拔最优秀的队员和科学合理地组队问题”, 以及第三问要求分析更多学生时的编排策略, 在建模中不应过早地将模型局限在 20 人选 6 只队伍上. 而应该建立一个普适性较强的模型. 再以20位同学的数据进行应用, 分析结论.

2.3 基本假设

(1) 学校通过对学生三种不同的能力评估 (定义为建模能力, 编程能力, 写作能力), 确定学生适合进行什么工作, 再进入各项工作的编制中;

(2) 学校根据高级别奖项优先原则, 在各个能力编制选出合适的队员组成一对;

(3) 假设所有学生在能力评估/竞赛中都能正常发挥;

(4) 假设每位同学都对参加选拔的同学的能力有一定的了解, 为了确保自己能被选上, 会选择对自己最有利的方向报名.

2.4 模型建立

2.4.1 学生单项能力综合评估

此部分较为简单, 能联想到需要根据不同的能力分别选取不同的指标, 建立三个评价模型即可. 当然, 作为建模团队而言, 建模同学能对写作、编程有一定的了解的话可能也是不错的选择, 即除了单独衡量一个方向外, 还可以结合上其他的能力值进行综合赋值.

我们假设三种方向主要考察以下指标:

建模能力评分: 科学成绩

、智力水平

、其他特长

编程能力评分: 智力水平

、动手能力

、其他特长

写作能力评分: 写作能力

、外语水平

、其他特长

TOPSIS法(优劣解距离法) 和 RSR(秩和比综合评价法) 都是简单方便的方法. 在评价建模能力的时候, 选择科学成绩、智力水平、其他特长作为评价指标, 使用 熵权法 + TOPSIS 方法进行评估, 设学生

从事

工作的能力评分为

, 则学生的单项能力得分矩阵为:

考虑到建模中团队协作能力的影响, 定义学生

的协作能力为

, 则学生的综合能力得分矩阵为:

也可以考虑某位同学选择建模方向, 则主要考虑建模能力, 编程能力和写作能力值折扣记分, (或其他的换算函数), 如:

当然, 也可以对 “其他特长” 进一步细化考虑. 毕竟... 学习成绩同样的人, 特长越多越好么~

此部分总体思路如下:

2.4.2 团队总能力评估

定义团队总能力为

. 设学生

组成同一只建模队伍

, 且学生

从事建模工作, 学生

从事编程工作, 学生

从事写作工作, 则团队的总能力为:

这里将团队各成员的单项能力得分作为总得分, 并人为规定顺序.

# 团队总能力评估

def team_score(model, code, write):

"""团队能力 = 建模能力 * 协作能力 + 编程能力 * 协作能力 + 写作能力 * 协作能力"""

return power_data.at[model, '建模'] + power_data.at[code, '编程'] + power_data.at[write, '写作']

2.4.3 获奖潜力评估

获奖概率难以量化, 但我们可以通过抽样仿真来进行估计.

假设该校所有报名的学生的能力得分分布与全体参加竞赛同学的得分分布近似. 这样我们就能通过对所有学生进行随机组合, 得到全体参加竞赛团队的成绩分布. 结合模型假设: 所有团队在竞赛中都能正常发挥, 那就可以按照团队实力的分布来代替他们比赛时的得分分布.

接着根据获奖比例, 设定三个获奖阈值, 只要团队总能力达到阈值水平就认为有可能获该奖. 达到

阈值认为有可能获得特等奖, 达到

阈值认为有可能获得一等奖, 达到

阈值认为有可能获得二等奖.实际论文写作的话, 这里需要插入一张流程图.

定义指标变量:

仿真部分具体实现请看代码: (这里对整体水平 + 20%, 实际上可以多设几组, 观察组队变化情况)

# 仿真: 重复实验, 假设该校水平与整体水平相差 20%

team_scores = []

for i in range(21000):

team = list(power_data.sample(3, replace=True, random_state=i).index)

team_scores.append(team_score(*team) * 1.2)

# 排序, 设 .43% 特等奖, 8.88% 一等奖, 37.97% 二等奖

team_scores.sort(reverse=True)

g0 = team_scores[int(21000 * 0.0043)]

g1 = team_scores[int(21000 * 0.0932)]

g2 = team_scores[int(21000 * 0.3829)]

2.4.4 队员组队分配模型建立

假设某次建模比赛选拔中需要从

名同学中选拔

只队伍,

与

应满足:

每位学生最多只能从事1项工作, 这是 0-1 指派问题. 同时, 要满足问题 (1) 选拔出优秀的队员, 再满足问题 (2) 进行组队分配, 这是个多目标优化模型.

定义决策变量:

定义第一级目标函数为全体同学的总能力得分最大化, 实际意义也很好理解, 优秀的团队应该从优秀的同学中产生:

定义第二级目标函数为获奖效益最大化:

10 : 5 : 1 是我瞎掰的, 可以自行设定. 也可以拆分为第三级、第四级目标.

约束条件:

① 每个团队每项工组都要有人负责

② 每一个人至多在一个队伍里负责1项工作

③ 每项任务都必须有m个人

④ 每个团队都必须有3个人

⑤ 获奖情况指标变量约束

实际上，约束 ③ 已经包含了 约束 ④

约束 ⑤ 中, 优化目标是全体

,

,

尽可能大, 在可能的情况下 每一个都会尽量取1, 但取1就必须要满足团队的总能力得分大于各个奖项的得分阈值. 因此, 这样的约束就能用来转化为获奖队伍数量.

2.5 代码及求解20选18的问题太简单了, 类似的改一下代码就行了

使用的编程语言：python3.7.1 (Anaconda3)

使用的编辑器：Sublime Text 3

使用的模块：pandas、statsmodels、scipy

代码一: topsis.py 文件

import pandas as pd

import numpy as np

def topsis(data, weight=None):

"""TOPSIS评价法"""

# 归一化

data = data / np.sqrt((data ** 2).sum())

# 最优最劣方案

Z = pd.DataFrame([data.min(), data.max()], index=['负理想解', '正理想解'])

# 距离

weight = entropyWeight(data) if weight is None else np.array(weight)

Result = data.copy()

Result['正理想解'] = np.sqrt(

((data - Z.loc['正理想解']) ** 2 * weight).sum(axis=1))

Result['负理想解'] = np.sqrt(

((data - Z.loc['负理想解']) ** 2 * weight).sum(axis=1))

# 综合得分指数

Result['综合得分指数'] = Result['负理想解'] / (Result['负理想解'] + Result['正理想解'])

Result['排序'] = Result.rank(ascending=False)['综合得分指数']

return Result, Z, weight

def entropyWeight(data):

"""熵权法"""

data = np.array(data)

# 归一化

P = data / data.sum(axis=0)

# 计算熵值

E = np.nansum(-P * np.log(P) / np.log(len(data)), axis=0)

# 计算权系数

return (1 - E) / (1 - E).sum()

代码二: main.py文件

# -*- coding: utf-8 -*-

# @Author: suranyi

# @Date: 2019-03-30 11:32:17

# @Last Modified by: suranyi

# @Last Modified time: 2019-03-30 12:54:51

import numpy as np

import pandas as pd

import gurobipy

import topsis

# %% 1.创建测试数据

# 创建 600 名同学

classmates = [f'同学{i}' for i in range(1, 601)]

# 组队队伍数

m = 100

# 能力

power = ['建模', '编程', '写作']

# 随机生成他们的成绩

np.random.seed(0)

data = pd.DataFrame(7 + np.random.randn(600, 7), index=classmates, columns=['科学成绩', '智力水平', '动手能力', '写作能力', '外语水平', '协作能力', '其他特长'])

# %% 2.能力值评估

# 单项能力评估

""" 熵权法赋权 + TOPSIS建模能力: (科学成绩、智力水平、其他特长) * 协作能力编程能力: (智力水平、动手能力、其他特长) * 协作能力写作能力: (写作能力、外语水平、其他特长) * 协作能力(建模能力 + 编程能力 + 写作能力) * 总协作能力 属于二次规划问题, 对于规模较小的模型可以快速解, 对于大规模模型难以求解"""

power_data = pd.DataFrame(index=classmates, columns=power)

power_data['建模'] = topsis.topsis(data.loc[:, ['科学成绩', '智力水平', '其他特长']])[0]['综合得分指数'] * data.loc[:, '协作能力']

power_data['编程'] = topsis.topsis(data.loc[:, ['智力水平', '动手能力', '其他特长']])[0]['综合得分指数'] * data.loc[:, '协作能力']

power_data['写作'] = topsis.topsis(data.loc[:, ['写作能力', '外语水平', '其他特长']])[0]['综合得分指数'] * data.loc[:, '协作能力']

# 团队总能力评估

def team_score(model, code, write):

"""团队能力 = 建模能力 * 协作能力 + 编程能力 * 协作能力 + 写作能力 * 协作能力"""

return power_data.at[model, '建模'] + power_data.at[code, '编程'] + power_data.at[write, '写作']

# 仿真: 重复实验, 假设该校水平与整体水平相差 20%

team_scores = []

for i in range(21000):

team = list(power_data.sample(3, replace=True, random_state=i).index)

team_scores.append(team_score(*team) * 1.2)

# 排序, 设 .43% 特等奖, 8.88% 一等奖, 37.97% 二等奖

team_scores.sort(reverse=True)

g0 = team_scores[int(21000 * 0.0043)]

g1 = team_scores[int(21000 * 0.0932)]

g2 = team_scores[int(21000 * 0.3829)]

# %% 3.Gurobipy求解

# 创建模型

MODEL = gurobipy.Model()

# 创建变量

price = ["Outstanding Winner/Finalist", "Meritorious Winner", "Honorable Mention", "Successful Participant"]

x = MODEL.addVars(classmates, power, m, vtype=gurobipy.GRB.BINARY, name="X")

p = MODEL.addVars(m, price, vtype=gurobipy.GRB.BINARY, name='Price')

# 更新变量环境

MODEL.update()

# 创建目标函数

f = [sum(x[i, j, k] * power_data.at[i, j] for i in classmates for j in power) for k in range(m)] # 每只队伍的总能力值

MODEL.setObjectiveN(- sum(f), priority=3, index=0) # Obj1: 团队总得分最大

MODEL.setObjectiveN(- 10 * sum(p.select('*', price[0])) - 5 * sum(p.select('*', price[1])) - sum(p.select('*', price[2])), priority=2, index=1) # Obj2: 特等奖/一等奖人数尽可能多, 加权

# 创建约束条件

MODEL.addConstrs(sum(x.select('*', j, k)) == 1 for j in power for k in range(m)) # 每个团队每项工作都要有人负责

MODEL.addConstrs(sum(x.select(i, '*', '*')) <= 1 for i in classmates) # 每一个人至多在一个队伍里负责1项工作

MODEL.addConstrs(sum(x.select('*', j, '*')) == m for j in power) # 每项任务都必须有m个人

MODEL.addConstrs(sum(x.select('*', '*', k)) == 3 for k in range(m)) # 每个团队都必须有3个人

MODEL.addConstrs(p[k, price[0]] * g0 <= f[k] for k in range(m))

MODEL.addConstrs(p[k, price[1]] * g1 <= f[k] for k in range(m))

MODEL.addConstrs(p[k, price[2]] * g2 <= f[k] for k in range(m))

# 自动排序结果, 但该条件太苛刻了, 无法产生最优解

# for k in range(m - 1):

# MODEL.addConstr(f[k] >= f[k + 1])

# 执行最优化, 由于 k 取值上是等价的, 会有很多组最优解长得一样, 只是 k 不同, 控制求解时长为300s

MODEL.Params.TimeLimit = 300

MODEL.optimize()

# 输出组队情况

Team = pd.DataFrame(columns=power, index=range(m))

for k in range(m):

for i in classmates:

for j in power:

if x[i, j, k].X:

Team.at[k, j] = i

Team['团队总得分'] = Team.apply(lambda item: team_score(*item), axis=1)

Team['奖项'] = Team.apply(lambda item: price[0] if item['团队总得分'] > g0 else price[1] if item['团队总得分'] > g1 else price[2] if item['团队总得分'] > g2 else price[3], axis=1)

Team.sort_values(by='团队总得分', inplace=True, ascending=False)

Team.index = range(1, m + 1)

Team.to_excel('组队结果.xlsx')

pd.set_option('display.expand_frame_repr', False)

pd.set_option('display.max_columns', None)

pd.set_option('display.max_rows', None)

print(Team)仿真结果

Gurobi 计算速度极快, 在使用pandas.DataFrame访问代价矩阵元素的时候, 一定要使用 .at 方法 (.at 方法的速度是.loc 的千倍以上).

2.6 模型推广

实际在分配中, 可能会出现某些同学希望能在同一只队伍进行比赛, 某些同学不希望在同一只队伍比赛, 为便于模型描述, 我们假设这些选择是双向的, 即

与

要么共同偏好与对方组队, 要么都厌恶与对方组队. 并且, 每位同学至多建立2个偏好关系。定义同学

的偏好集合

和厌恶集合

.

则有:

这样的组队偏好也有一定风险: 如果大佬要带2只菜鸡, 团队总能力/获奖期望达不到最低要求, 可能会集体淘汰.

2.7 总结

以上就是本题的完整思维过程 + 代码实现, 不是很难吧? 从选拔队员、组队上都合理地分析, 仿真试验, 评估建模, 基本上能解决所有提出的问题, 在选拔规模 <1000 的情况下都能快速求解, > 1000 的情况我没具体分析. 另外, 不管是 gurobi 还是 lingo, 对线性模型对求解速度是非常快的, 二次规划问题就不是那么好做了. 比如将团队能力评估模型设为:

该式的含义为, 完成一篇论文需要三个人共同的努力, 三个人的协作能力才是一个团队总能力能发挥的关键. 看似很有道理, 但这个模型是二次规划问题, 数据量大了, 就根本跑不出来了. (可以考虑做一个模型对比, 或者求解效率分析)

最重要的是, 建模的出发点不要太低, 太低的话你的模型推广、模型的优缺点就没法“吹”了. 我们这里写的是

个人选拔

只队伍, 对于其他团队比赛是不是也适用呢?

3 Problem 2 参会安排问题 —— 同济大学校内数模竞赛2018年 A题

3.1 问题描述

3.2 问题分析

最小代价出行安排问题本质上是算法设计问题, 如中国邮递员问题. 但在实际应用/建模中, 不仅要考虑时间上冲突、经费限制等硬性约束, 会议影响力权衡、教师出行交通工具选择、住宿方案也是需要我们考虑的。

出行安排问题是NP难问题, 在规模量较大的时候难以快速求解. 本题中, 需要安排 13 个地点, 18位教师, 如果一次性遍历所有方案再挑出最优解, 程序跑不完, 如果采用贪心算法, 不一定能取得最优解. 考虑到该问题中教师去不去某个地方开会, 可以构成 0-1 变量, 而各种条件就是模型的约束, 由此想到可以使用 0-1 整数规划求解. 但直接使用 0-1 整数规划, 约束太多, 很多条件也不能很好的表示, 是不是还存在着其他的方法呢? (如果对0-1规划解法感兴趣, 可以参考该校的优秀论文).

会议时间限制是这道题区别于常规TSP问题的关键. 有会议时间限制, 教师的出行安排地点方案其实是由时间约束好了的. 比如说去了北京就一定不可能去上海, 最长的路径是 "北京-兰州-昆明-南京-杭州-济南-大连", 如果再考虑其他的模型假设 (如: 参加会议的教师需要提前一天报道), 那么实际上教师可以选择的路径数量是很少的 (可能只有 1, 2 百来条路径). 这样想, 那就可以把教师去哪个地点以及地点间的相互组合, 转换为教师可以选择哪一条既定的路径方案的问题. 模型本身也将“时间不能冲突”约束转化为“可以取的解”问题.

到这里我们已经有了初步的设想, 先把所有可能的路径方案列出来, 转为路径分配问题. 当然, 只要路径能完整的确定下来, 它的代价就是完全确定的 (而常规 0-1 变量没办法解决这个问题, 或者说能解决但需要考虑的约束式太多. 即使相邻地点确定下来, 乘坐什么交通工具, 是先去目标地多付目标地点住宿费, 还是留在现在的地方, 多付这个地方的住宿费? ). 只要满足各种路径中去某一些地点的教师职称满足参会最低要求即可. 再设 0-1 变量, 表示教师执不执行这条路径方案, 而约束就是满足最低参会要求和个别教师参会数量限制.

接着我们要考虑,怎么样列出可能的路径, 以及结合实际来思考一些隐含的约束. 我们先举两个例子, 便于分析规律:如果某位教师的路径是 “北京 - 咸阳”, 这位教师是北京开完会后留在北京, 还是提前去咸阳呆在咸阳? 此外, 对于时间间隔较长的地方, 是否还存在着"回校待命"的选择? 即: "北京 - 上海待命 - 咸阳". (上海不需要支付住宿费, 能节省一大笔开销)

还是以 "北京 - 咸阳" 为例, 该教师相邻两次会议的时间间隔足够长, 那么该教师是乘坐飞机还是乘坐火车呢? (这个问题非常重要)

这两个问题其实也反映了建模中需要考虑的问题: 列出的路径是指所有可能的路径, 还是最经济实惠的路径? 从目标函数为最小化总费用的角度出发, 如果同一条路径能满足参会需求, 但是有不同的交通出行方案、住宿方案, 那么我们自然希望选择开销最小的一条. 也就是说, 对于 "北京-咸阳" 这样的路径, 甭管你怎么走, 只要开销最小就行. 可能的路径是指可能的最小费用路径.笔者认为, 全面地思考一些极端的例子, 有助于简化模型和更细致地分析问题.

这样一考虑, 就清晰多了, 但似乎还有点不对劲? 如果我们考虑两点间的球面距离作为城市的距离, 火车的费用可是一直是比飞机还要低的! 按照我们上面的设想, 应该会导致所有教师都是坐火车出行. 要是真的这么安排, 估计该校教师会说“这学校也太抠门了”吧. 另外, 上海有虹桥机场、浦东国际机场, 我坐飞机是回哪个机场, 距离又要怎么算呢? 如果是北京坐火车去咸阳, 那还可以省下好几天的住宿费呢... 且慢, 建模可不是让你考虑这么多无关紧要的问题的, 要不然也不会把会场交通费写上去了! 这个地方就需要合理的假设. 比如, 教师不能连续乘坐12小时火车(别说太娇气了, 老师们可是去参加大型国际学术会议的, 坐火车两三天不洗澡吃泡面, 像话吗! 另外, 如果想考虑的很仔细的话, 你还得考虑什么时候发车什么时候到达的问题呢! ), 教师开会前一天必须到会场进行签到, 教师在会议结束次日才能离开该地 (比如北京 7.20-7.25, 则必须在 7.19 到达北京, 7.26 才能离开北京去下一个城市, 最早参加 7.27 的会议. 欧吼, 这样以来去了北京的根本就不可能去兰州、成都, 减少了好多条路径!)

以上就能解决第1问的问题, 接着第2问就要对会议影响力进行评估. 这个时候没有参会人数要求的限制了, 可不可以采用攒星星的模式呢? (去了这个会议就拿到对应的星星, 直接相加). 再往后看第3问暗示我们, 不同职称的教师对会议影响力的量化评价是不同的, 人数也有 "至少" 的限制.为了便于分析, 我们可以理解为: 会议没有最低参会人数要求, 但是只有满足 "最低参会人数要求" 的教师团队才能进行大会报告, 对应的会议也才能攒到会议影响力.这样思考的好处就是, 问题直接转化为在5万经费下, 去不去参加某个会议的问题, 以及对应的影响力评估.

总结以上,本问题可选取的模型: 0-1 整数规划 + 大量约束; 路径分配 + 随机漫步遗传算法; 路径分配 + 0-1 整数规约 + 少量约束. 需要解决如下问题:

距离转换 -> 构造最小代价候选路径列表 -> 影响力评估、代价评估 -> 0-1整数模型 -> 求解算法设计以上只是笔者的看法, 充分结合实际简化模型. 也没有刻意忽略题干的某些条件.

你也可以通过查找比较详细的距离, 如果花费时间超过12小时, 则教师需要提前一天, 以此类推, 通过时间来体现交通工具的区别. (不过我感觉这样太麻烦了, 本文的重点不在于怎么分析距离)思维过程

3.3 基本假设

(1) 教师在会议时间内必须留在会场地点, 不能中途离开. 因此, 若参加 start-end 的会议, 教师在 start - 1 天的时候到达会议地点, end + 1 天的时候才能离开会议地点.

(2) 教师到达会议地点当天就产生了住宿费用, 离开会议地点的时候不产生住宿费用

(3) 若教师需要参加的两场会议时间间隔较长, 允许教师回到上海待命, 并且教师在上海不产生注册费用、住宿费用、会场交通费, 只有回到上海的交通费.

(4) 交通费用包括火车/飞机交通费和会场交通费. 假设城市间的距离为两城市中心的球面距离, 从火车站/机场到达会场的费用都包含在会场交通费中. (此外, 会场交通费包括来回的费用, 即一个会议地点的交通费只收1次钱)

(5) 不考虑自然灾害导致交通工具晚点的问题

(6) 考虑到教师的休息时间即会议的各项准备工作, 教师乘坐交通工具的时间不能超过12小时.

(7) 查阅资料, 全国列车的平均时速为 70.18 公里/小时 (参考文献[1]), 飞机的平均时速在 800公里/时 ~ 1000公里/时 之间 (参考文献[2]). 我国领土东西距离约5200公里,南北距离约为5500公里, 因此乘坐飞机一定能在 12 小时内到达.

(8) 每个学校只能有一位教师的学术报告被选为大会报告,(6) 还可以细化, 若教师结束会议返回上海, 则一律乘坐火车. 本文不考虑这一点, 统一标准.

3.4 模型建立先将数据转化为该表形式, 便于程序编写该文件记为 Base_Data.xlsx

定义几个全文通用的符号: :

为会议地点,

为会议地点集合.

: 出发点、终点. 本题中出发点和终点都是上海, 即:

: 会议的开始时间、结束时间

:

会议的住宿费用

:

会议的注册费用

:

会议的会场交通费用

:

会议的影响力 (星级)

:

会议要求的参会教授人数

:

会议要求的参会副教授或教授人数

: 火车的平均速度, 飞机的平均速度

:

会议城市到

会议城市的球面距离

:

会议城市到

会议城市的交通费 (火车或飞机费用)

: 教师, 其中定义 1~5 教师为教授, 6 ~ 13 教师为副教授, 14~18 教师为讲师, 且 1 为教研室主任, 2 为副主任

3.4.1 地理编码与距离计算

使用百度地理编码工具, 获取这13个城市的经纬度坐标, 再计算两点的球面距离(经纬度跨度较大, 最好不要直接用直线距离). 设火车的平均速度为

, 飞机的平均速度为

,

到

的球面距离为

, 则

到

的交通费为

有关地理编码与距离计算可参考:https://zhuanlan.zhihu.com/p/60945589​zhuanlan.zhihu.com

代码如下:

import pandas as pd

from Geocoders import Geocoders

from Geodistance import geo_distance

# 导入数据

data = pd.read_excel('../table/Base_Data.xlsx', index_col='地点')

# 查询经纬度

geo = Geocoders(data.index, api_key=input("Your Baidu's api_key: "))

latitude = geo.export2series().str.split(',').apply(lambda item: item[0]).astype(float)

longitude = geo.export2series().str.split(',').apply(lambda item: item[1]).astype(float)

latitude.name = '纬度'

longitude.name = '经度'

data = data.join(latitude).join(longitude)

# 获取价格矩阵

distance_cost = pd.DataFrame(index=data.index, columns=data.index)

transportation = distance_cost.copy()

for i1 in data.index:

for i2 in data.index:

distance = geo_distance((data.at[i1, '纬度'], data.at[i1, '经度']), (data.at[i2, '纬度'], data.at[i2, '经度']))

if distance / 70.18 > 12:

print(f'{i1}到{i2}太远啦! 坐飞机吧 ~')

transportation.at[i1, i2] = '飞机'

distance_cost.at[i1, i2] = distance * 0.8

else:

transportation.at[i1, i2] = '火车'

distance_cost.at[i1, i2] = distance * 0.5

# 导出数据

distance_cost.to_excel("../table/distance_cost.xlsx")

transportation.to_excel("../table/transportation.xlsx")

怕模型不够准确? 看看下面俩张图.飞机距离差不多, 火车距离偏小, 无碍, 反正设置了不超过12小时, 并且教师要提前到达, 有一定的容错

3.4.2 单次会议基本开销

记会议地点为

,每一个地点对应的会议时间为

, 住宿费用为

, 注册费用为

, 会场交通费为

, 会议影响力

.

单次会议基本开销指参加某场会议所需的基本费用. 包含会议的住宿费用、注册费用、会场交通费. 根据模型假设(1), 会议的住宿费用包括注册当天的住宿费用及会议全程的住宿费用. 因此, 参加会议

的基本开销

为:

data['基本费用'] = data['住宿费用'] * (data['结束时间'] - data['开始时间'] + 2) + data['注册费用'] + data['会场交通费']

data.to_excel("../table/data.xlsx")

引入基本开销之后, 我们不再考虑会场交通费用. 下文所说的交通费指的是 "火车/飞机" 的费用.

3.4.3 邻接路径方案及其代价

邻接路径指从地点

到地点

的住宿、交通方案. 这里的住宿是非会议时间住宿, 我们称之为"滞留代价". 先定义"可达"概念, 若

和

满足:

那就称

是"可达"的, 即是一条时间上不冲突的可行路径. 邻接路径是从可达路径中产生的, 否则生成的路径是没有意义的.

定义符号

表示

地会议与

地会议的净时间间隔.

从

有如下几种情形:若

, 即两场会议时间时间间隔2天, 则该路径方案为

飞机/火车直达

. 交通费用为

, 滞留代价为

若

, 即两场会议间隔2天以上, 则该路径方案可以为以下三种之一:在

地留宿

天, 则该路径方案为

(滞留

天)

飞机/火车直达

. 交通费用为

, 滞留代价为

在

地留宿

天, 则该路径方案为

飞机/火车直达

(滞留

天)

. 交通费用为

, 滞留代价为

从

地飞回上海 (返校待命), 在校滞留

天, 再前往

地. 则该路径方案为

上海中转

天

. 交通费用为

, 滞留代价为

定义邻接路径代价为:

# 构造可达路径

reachability_roads = []

for start in data.index:

for end in data.index:

if data.at[end, '开始时间'] - data.at[start, '结束时间'] >= 2 and start != end:

reachability_roads.append([start, end])

3.4.4 最小费用邻接路径

最小费用候选路径列表是由一系列最小费用候选路径构成. 而最小代价路径指完成一整条不冲突的地点(基本开销是固定的)及其最小交通、滞留方案构成(由3.4.3知道, 根据不同的方案会有不同的开销).

若

是"可达"的, 参加

和

会议的最小交通-滞留代价是:

这两个地点构成的路径的最小费用为: 基本费用和最小交通-滞留代价之和, 即

定义标记符号为

, 用来表示从

到

的最小费用路径方案,

使用箭头表示参加会议的顺序, 引入记号

. : 会议地点

: 提前到达

的天数, 标记为 "提前

天"

: 停留在

的天数, 标记为 "滞留

天". 如果

包含了回到

待命, 则将

所在的分枝改为"中转

天"

: 交通工具, 火车或飞机

例如:

对应的开销为:

# 构造最小费用邻接路径

shortest_nearest_roads = []

for road in reachability_roads:

start, end = road

if data.at[end, '开始时间'] - data.at[start, '结束时间'] == 2:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, end]}]->{end}", data.at[start, '基本费用'] + distance_cost.at[start, end] + data.at[end, '基本费用']])

else:

# 滞留天数

wait_day = data.at[end, '开始时间'] - data.at[start, '结束时间'] - 2

# start 地驻留

way1_cost = data.at[start, '基本费用'] + data.at[start, '住宿费用'] * wait_day + distance_cost.at[start, end] + data.at[end, '基本费用']

# end 地驻留

way2_cost = data.at[start, '基本费用'] + data.at[end, '住宿费用'] * wait_day + distance_cost.at[start, end] + data.at[end, '基本费用']

# 中转

way3_cost = data.at[start, '基本费用'] + distance_cost.at[start, '上海'] + distance_cost.at['上海', end] + data.at[end, '基本费用']

# 最小费用邻接路径方案代价

min_cost = min(way1_cost, way2_cost, way3_cost)

# 储存最小费用邻接路径及代价

if min_cost == way1_cost:

shortest_nearest_roads.append([f"{start}(提前{wait_day}天)[{transportation.at[start, end]}]->{end}", way1_cost])

elif min_cost == way2_cost:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, end]}]->(滞留{wait_day}天){end}", way2_cost])

else:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, '上海']}]->上海(中转{wait_day}天)[{transportation.at['上海', end]}]->{end}", way3_cost])

3.4.5 最小费用候选路径列表

由于时间先后及会议的相对独立性, 最小费用候选路径在每一个独立的会议阶段都是最小费用邻接路径 (如 从北京到昆明到咸阳这一条路径的子路径中, 北京到昆明的滞留方案、出行方式并不会影响昆明到咸阳的方案, 整条路径最优必然有子路径 北京到昆明 是最优的, 昆明到咸阳也是最优的). 因此, 我们可以先构造出所有的2个地点会议的最小费用邻接路径, 并按照首位组合的方式更新3个会议地点、4个会议地点、5个会议地点, ... 的情况, 由此得到所有路径的出行安排情况. 由于每个会议至少占用 4 天时间 (提前一天到达, 迟一天离开), 22 天内最多参与进行 5 场会议.

该部分具体流程如下:

Step1先根据所有的2个会议地点的最小费用邻接路径, 判断是否可首位相接, 即如果有

和

, 那么:

由此生成所有具有3个会议地点的最小费用路径

Step2再根据所有的3个会议地点的最小费用邻接路径与2个会议地点的最小费用邻接路径, 判断是否可首位相接, 即如果有

和

, 那么:

由此生成所有具有4个会议地点的最小费用路径

.

Step3根据所有的4个会议地点的最小费用邻接路径与2个会议地点的最小费用邻接路径, 判断是否可首位相接, 即如果有

和

, 那么:

由此生成所有具有5个会议地点的最小费用路径

.

Step4将上述得到的所有2个会议地点的最小费用路径集合、3个会议地点的最小费用路径集合、4个地点的最小费用路径集合、5个会议地点的最小费用路径集合合并, 得到了 2~5 个会议的最小费用路径集合.

# 拼接最小费用邻接路径, 以下可以用while True, 迭代创建, 为便于读者理解, 此处拆开

shortest_nearest_roads_2 = shortest_nearest_roads.copy()

shortest_nearest_roads_3 = []

for road1 in shortest_nearest_roads_2:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_3.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads_4 = []

for road1 in shortest_nearest_roads_3:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_4.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads_5 = []

for road1 in shortest_nearest_roads_4:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_5.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads.extend(shortest_nearest_roads_3 + shortest_nearest_roads_4 + shortest_nearest_roads_5)总共是136条路径!至此, 再回顾我们的工作: 根据时间不冲突原则构造可达路径, 根据可达路径对应的多种交通工具、滞留方案并搜寻最小费用邻接路径, 再将最小费用邻接路径首尾拼接为最小费用候选路径

当然, 这样的路径费用并不是完整的. 我们忽略了只去1个地点的路径, 也没有将出发程和返程的路径代价加上. 下面就要进行这样的工作, 获取完整的最小费用候选路径列表及其代价.从

出发, 完成最小费用候选路径再回到

的总费用按照如下方式进行计算:

# 获取最小费用邻接路径列表

# 先加入所有只去1个地点的路径

shortest_roads_list = []

for location in data.index:

shortest_roads_list.append([f'{location}', data.at[location, '基本费用']])

shortest_roads_list += shortest_nearest_roads

# 再为所有的路径添加出发程、返程路径代价

for road in shortest_roads_list:

road[1] += distance_cost.at['上海', road[0][:2]] + distance_cost.at[road[0][-2:], '上海']

road[0] = f"[出发]上海[{transportation.at['上海', road[0][:2]]}]->{road[0]}[{transportation.at[road[0][-2:], '上海']}]->上海[返回]"

3.4.6 路径分配模型建立

在建立模型前, 我们先将上述路径拆分成包含地点与否的独热编码矩阵. 这样处理的好处是可以直接从矩阵上读取约束, 便于代码编写, 此外还便于统计该路径的影响力、总费用. 记该矩阵为

, 矩阵的每一行都是一条 "最小费用候选路径", 第1列元素为该路径对应的代价 (我们简记

表示第

条候选路径), 第2列~第14列对应路径

是否包含着地点

的0-1示性函数

.

roads = pd.DataFrame({'价格': dict(shortest_roads_list)})

for location in data.index:

# 注意, 这里要先去掉首尾

roads[location] = roads.index.str.split('->').map(lambda items: int(location in ''.join([item for item in items[1:-1] if "中转" not in item])))

roads.index.name = '路径'

roads.to_excel("../table/shortest_roads_list.xlsx")

(1) 构建决策变量

定义路径分配的 0-1 变量

:

(2) 构建多水平目标函数

1) 总费用水平

2) 总影响力水平

教师参与会议

的影响力评估因子为:

总会议影响力评估因子为:

这样考虑主要是想拉开五星级会议和四星级的差距、四星级和三星级的差距, 只要参加了会议就认为能获得相应的会议影响力

这里可以引入最高职称教师对会议影响力的影响, 进行不同的换算. 具体可以参考下面第三问的模型.

(3) 全局约束条件

定义学校是否组织教师参与

地会议的 0-1 变量

:

每位教师最多只能安排一条参会路径:

主任和副主任至多参加3个会议:

若学校组织参加某次会议, 则参会教师及职称需要满足会议的最低要求:

这里的

是一个很大的数, 这里可以取18, 表示同一个地点最多只能去18个人. 约束 (24) 表明, 如果该学校组织教师参加

会议, 则在足够大的

下该约束无效; 若该学校不组织教师参加

会议, 则

, 此时不论是哪位教师, 都不能选择包含地点

的路径

.式 21 代表若参加该会议, 则教授人数要达到要求; 式 22 代表若参加该会议, 则副教授人数要达到要求, 式 23 代表若参加该会议, 则总教师人数要达到要求.

3.5 针对性模型扩展

针对性问题求解的约束是全局约束+对应题意约束构成. 这样的写法可提高模型复用性、易迁移性, 并且条理清晰, 避免大量重复工作(如果你想凑字数的话就多复制几遍吧).

3.5.1 第一问题目没有说一定要全部会议都去. 但如果认为可以有某些会议不去, 那有可能所有教师都去同一个地点(最省钱的路径), 不符合 “提升学校影响力” 的要求. 因此, 这里应该要体现出所有会议都要参加.

此外题目中的 "科研室主任和副主任两人" 应这么断句 "科研室主任和副主任 | 两人".

每位教师至少参加两个会议:

当然, 我们可以在构建决策变量的时候, 直接将只去了一个地点开会的路径去掉.

所有会议都要参加:

目标函数:

模型求解结果, 总费用 125630.20 元

import pandas as pd

import gurobipy

# %% 1. 数据准备

# 导入文件

shortest_roads = pd.read_excel("../table/shortest_roads_list.xlsx")

data = pd.read_excel("../table/data.xlsx", index_col="地点")

# 教师

T = ['教授1[主任]', '教授2[副主任]', '教授3', '教授4', '教授5',

'副教授1', '副教授2', '副教授3', '副教授4', '副教授5', '副教授6', '副教授7', '副教授8',

'讲师1', '讲师2', '讲师3', '讲师4', '讲师5']

# 会议地点

I = data.index

# 路线编号

R = shortest_roads.index

# %% 2. gurobipy 求解

MODEL = gurobipy.Model()

# 创建变量, 分开写便于后文代码编写

X = MODEL.addVars(T, R, vtype=gurobipy.GRB.BINARY, name="Teacher")

Z = MODEL.addVars(I, vtype=gurobipy.GRB.BINARY, name='Z')

# 更新变量环境

MODEL.update()

# 创建目标函数 —— 总费用水平最小化

MODEL.setObjective(sum(X[t, j] * shortest_roads.at[j, '价格'] for t in T for j in R), gurobipy.GRB.MINIMIZE)

# 创建约束条件

# 全局约束

MODEL.addConstrs(sum(X[t, j] for j in R) <= 1 for t in T) # 每位教师至多只能安排一条参会路径

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for i in I for j in R) <= 3 for t in T[:2]) # 主任和副主任至多参加3个会议

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= Z[i] * data.at[i, '教授要求'] for i in I) # 最低教授参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:13] for j in R) >= Z[i] * data.at[i, '副教授及以上要求'] for i in I) # 最低副教授及教授参会人数要

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) >= Z[i] * data.at[i, '最低人数要求'] for i in I) # 最低参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) <= Z[i] * 18 for i in I) # 若不参加该会议, 则全部教师都不去

# 扩展模型约束

MODEL.addConstrs(Z[i] == 1 for i in I)

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for i in I for j in R) >= 2 for t in T)

# 执行线性规划模型

MODEL.optimize()

# 输出模型结果

result = pd.DataFrame(index=T, columns=I)

for t in T:

for r in R:

if X[t, r].X:

result.at[t, '方案'] = shortest_roads.at[r, '路径']

result.at[t, '费用'] = shortest_roads.at[r, '价格']

for i in I:

if shortest_roads.at[r, i]:

result.at[t, i] = 1

result.to_excel('../table/Question1.xlsx')

3.5.2 第二问所有最小费用候选路径中, 总费用最小的路径是

对应的费用是:

, 若所有教师都按照第一问的约束(至少参加两个会议), 则总开销已经达到35280元. 而其他路径都明显大于 980, 因此, 这一问是允许有教师不参加会议的.

总费用不超过5万元:

笔者认为, 学校参会最低要求是写在题干里的, 应该是全局约束. 而 "每位教师至少要参加2个会议" 是第1问的条件.

这里但是这里保留了教研室主任、副主任的出行次数约束, 笔者认为即使是题目没有涉及, 也可以往这方面进行推广, 读者可以自行甄别是否应该删除 (小细节, 不case啦)

目标函数:

模型求解结果, 剩余经费 4244.24 元, 会议总影响力 107

import pandas as pd

import gurobipy

# %% 1. 数据准备

# 导入文件

shortest_roads = pd.read_excel("../table/shortest_roads_list.xlsx")

data = pd.read_excel("../table/data.xlsx", index_col="地点")

# 教师

T = ['教授1[主任]', '教授2[副主任]', '教授3', '教授4', '教授5',

'副教授1', '副教授2', '副教授3', '副教授4', '副教授5', '副教授6', '副教授7', '副教授8',

'讲师1', '讲师2', '讲师3', '讲师4', '讲师5']

# 会议地点

I = data.index

# 路线编号

R = shortest_roads.index

# %% 2. gurobipy 求解

MODEL = gurobipy.Model()

# 创建变量, 分开写便于后文代码编写

X = MODEL.addVars(T, R, vtype=gurobipy.GRB.BINARY, name="Teacher")

Z = MODEL.addVars(I, vtype=gurobipy.GRB.BINARY, name='Z')

# 更新变量环境

MODEL.update()

# 创建目标函数 —— 总费用水平最小化

MODEL.setObjective(sum(data.at[i, '会议影响力'] ** 2 * Z[i] for i in I), gurobipy.GRB.MAXIMIZE)

# 创建约束条件

# 全局约束

MODEL.addConstrs(sum(X[t, j] for j in R) <= 1 for t in T) # 每位教师至多只能安排一条参会路径

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for i in I for j in R) <= 3 for t in T[:2]) # 主任和副主任至多参加3个会议

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= Z[i] * data.at[i, '教授要求'] for i in I) # 最低教授参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:13] for j in R) >= Z[i] * data.at[i, '副教授及以上要求'] for i in I) # 最低副教授及教授参会人数要

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) >= Z[i] * data.at[i, '最低人数要求'] for i in I) # 最低参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) <= Z[i] * 18 for i in I) # 若不参加该会议, 则全部教师都不去

# 扩展模型约束

MODEL.addConstr(sum(X[t, j] * shortest_roads.at[j, '价格'] for t in T for j in R) <= 50000) # 费用不超过5万元

# 执行线性规划模型

MODEL.optimize()

# 输出模型结果

result = pd.DataFrame(index=T, columns=I)

for t in T:

for r in R:

if X[t, r].X:

result.at[t, '方案'] = shortest_roads.at[r, '路径']

result.at[t, '费用'] = shortest_roads.at[r, '价格']

for i in I:

if shortest_roads.at[r, i]:

result.at[t, i] = 1

print(f"会议影响力 ={MODEL.ObjVal}")

print(f"剩余经费={50000 - result.loc[:, '费用'].sum()}")

result.to_excel('../table/Question2.xlsx')

3.5.3 第三问设教授的学术报告被选为大会报告的概率为

, 副教授为

, 讲师为

, 则由题意有:

, 无解. 说明不要想太多, 这不是概率论的问题. 直接假设每个学校只能有一位教师的学术报告被选为大会报告, 算期望就行了.

若被选为大会报告, 则可以获得对应的影响力评估因子；若在分会场进行报告, 则影响力评估因子减半. 以影响因子的期望替代:

仿造参数

, 设定以下几组 0-1 变量:75% -> 参加

会议有2位教授的指示变量为

, 它的含义为至少包含2位教授

50% -> 参加

会议有1名教授和一名副教授的指示变量为

, 它的含义为有且仅有1位教授和至少包含1位副教授

35% -> 参加

会议有1名教授和其余均为讲师的指示变量为

, 它的含义为有且仅有1位教授和不能有副教授

35% -> 参加

会议有2名副教授和其余均为讲师的指示变量为

, 它的含义为至少包含2名副教授和不能有教授

10% -> 其余情况的指示变量为

, 它的含义为没有教授和至多包含1位副教授

则任意一场会议, 都会满足:

表示需要安排参加该会议, 此时只有一个

;

表示不需要安排参加该会议, 此时所有的

则有以下约束 (这些约束都是根据上面5小点得到的, 可以根据

上角标的标记找到对应语句, 看看是啥意思~ 这也是很重要的约束技巧):

将式 (17) 定义的

进一步拆分为五个分量, 表示在五个指标变量

下的取值.

本题中已经提供了不同职称教师对会议影响力的影响, 我们就不再使用自定义的职称影响力, 将会议影响力更新为:

各

是表达式，不是变量。

总会议影响力为:

优化目标是使用尽可能少的钱获得尽可能多的会议影响力. 这样的表述其实有点问题, 因为没有指明钱究竟多少才算 "少". 一种处理方法是加权组合, 第二种处理方法是使单位影响力的代价尽可能小, 即

, 但这个优化目标不是线性的, 不好求解, 需要进行线性化处理. 如

, 则实际优化目标是

. 第三种处理方法就是使用分层序列法求解. 我们设定第一级优化目标是使得总会议影响力最大, 第二级优化目标是使得会议总开销最小. 这样就合理多了! 具体表达式如下:

模型求解结果, 总费用 113912.00 元, 会议总影响力 149.95

import pandas as pd

import gurobipy

# %% 1. 数据准备

# 导入文件

shortest_roads = pd.read_excel("../table/shortest_roads_list.xlsx")

data = pd.read_excel("../table/data.xlsx", index_col="地点")

# 教师

T = ['教授1[主任]', '教授2[副主任]', '教授3', '教授4', '教授5',

'副教授1', '副教授2', '副教授3', '副教授4', '副教授5', '副教授6', '副教授7', '副教授8',

'讲师1', '讲师2', '讲师3', '讲师4', '讲师5']

# 会议地点

I = data.index

# 路线编号

R = shortest_roads.index

# %% 2. gurobipy 求解

MODEL = gurobipy.Model()

# 创建变量, 分开写便于后文代码编写

X = MODEL.addVars(T, R, vtype=gurobipy.GRB.BINARY, name="Teacher")

Z = MODEL.addVars(I, vtype=gurobipy.GRB.BINARY, name='Z')

beta = MODEL.addVars(I, range(1, 6), vtype=gurobipy.GRB.BINARY, name='beta')

# 更新变量环境

MODEL.update()

# 创建目标函数 —— 总费用水平最小化

F_1 = [beta[i, 1] * (0.75 + 0.25 / 2) * data.at[i, '会议影响力'] ** 2 for i in I]

F_2 = [beta[i, 2] * (0.5 + 0.5 / 2) * data.at[i, '会议影响力'] ** 2 for i in I]

F_3 = [beta[i, 3] * (0.35 + 0.65 / 2) * data.at[i, '会议影响力'] ** 2 for i in I]

F_4 = [beta[i, 4] * (0.35 + 0.65 / 2) * data.at[i, '会议影响力'] ** 2 for i in I]

F_5 = [beta[i, 5] * (0.1 + 0.9 / 2) * data.at[i, '会议影响力'] ** 2 for i in I]

MODEL.setObjectiveN(- sum(F_1 + F_2 + F_3 + F_4 + F_5), priority=1, index=0)

MODEL.setObjectiveN(sum(X[t, j] * shortest_roads.at[j, '价格'] for t in T for j in R), priority=0, index=1)

# 创建约束条件

# 全局约束

MODEL.addConstrs(sum(X[t, j] for j in R) <= 1 for t in T) # 每位教师至多只能安排一条参会路径

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for i in I for j in R) <= 3 for t in T[:2]) # 主任和副主任至多参加3个会议

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= Z[i] * data.at[i, '教授要求'] for i in I) # 最低教授参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:13] for j in R) >= Z[i] * data.at[i, '副教授及以上要求'] for i in I) # 最低副教授及教授参会人数要

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) >= Z[i] * data.at[i, '最低人数要求'] for i in I) # 最低参会人数要求

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T for j in R) <= Z[i] * 18 for i in I) # 若不参加该会议, 则全部教师都不去

# 扩展模型约束

MODEL.addConstrs((1 - Z[i]) + sum(beta.select(i, '*')) == 1 for i in I) # 约束30

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= 2 * beta[i, 1] for i in I) # 约束31

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= beta[i, 2] for i in I) # 约束32

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) <= 5 - 4 * beta[i, 2] for i in I) # 约束32

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[5:13] for j in R) >= beta[i, 2] for i in I) # 约束33

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) >= beta[i, 3] for i in I) # 约束34

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) <= 5 - 4 * beta[i, 3] for i in I) # 约束34

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[5:13] for j in R) <= 8 * (1 - beta[i, 3]) for i in I) # 约束35

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) <= 5 * (1 - beta[i, 4]) for i in I) # 约束36

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[5:13] for j in R) >= 2 * beta[i, 4] for i in I) # 约束37

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[:5] for j in R) <= 5 * (1 - beta[i, 5]) for i in I) # 约束38

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[5:13] for j in R) <= 8 - 7 * beta[i, 5] for i in I) # 约束39

MODEL.addConstrs(sum(X[t, j] * shortest_roads.at[j, i] for t in T[13:] for j in R) >= 2 * beta[i, 5] for i in I) # 约束40

# 执行线性规划模型

MODEL.optimize()

# 输出模型结果

result = pd.DataFrame(index=T, columns=I)

for t in T:

for r in R:

if X[t, r].X:

result.at[t, '方案'] = shortest_roads.at[r, '路径']

result.at[t, '费用'] = shortest_roads.at[r, '价格']

for i in I:

if shortest_roads.at[r, i]:

result.at[t, i] = 1

for index, i in enumerate(I):

result.at['预估影响力', i] = (F_1[index] + F_2[index] + F_3[index] + F_4[index] + F_5[index]).getValue()

print(f"会议影响力 = {-MODEL.ObjVal}")

print(f"经费={result.loc[:, '费用'].sum()}")

result.to_excel('../table/Question3.xlsx')

目录结构如下图所示:最佳出行安排codeGeocoders.py

Geodistance.py

data_preprocessing.py

shortest_roads.py

Question1.py

Question2.py

Question3.py

tableBase_Data.xlsx

代码一: Geocoders.py, 地理编码与反编码工具.

import pandas as pd

from geopy.geocoders import Bing, Baidu

class Geocoders(object):

def __init__(self, addresses: "[str/(latitude, longitude)]", api_key: str, timeout: int=5, api: 'Baidu / Bing'='Baidu'):

self.values = dict.fromkeys(addresses, "N/A")

self.timeout = timeout

self.gps = self.get_gps(api, api_key)

self.shape = len(str(len(self.values)))

self.translate()

def __encode(self, address: str) -> "latitude, longitude":

location = self.gps.geocode(address, timeout=self.timeout)

return "%.2f,%.2f" % (location.latitude, location.longitude)

def __decode(self, address: str) -> "address":

location = self.gps.reverse(address, timeout=self.timeout)

return location.address if type(location.address) is str else location.address.decode()

def get_gps(self, api, api_key: str) -> "Geopy Server":

if api.lower() == "baidu":

return Baidu(api_key)

elif api.lower() == "bing":

return Bing(api_key)

def translate(self):

for n, address in enumerate(self.values, 1):

if self.values[address] == "N/A":

try:

self.values[address] = self.__encode(address) if isinstance(address, str) else self.__decode(address)

except Exception:

print(f"Warming: {n:{self.shape}}{address}查询失败")

finally:

print(f"{n:{self.shape}}{address}查询结果为{self.values[address]}")

else:

print(f"{n:{self.shape}}{address}查询结果为{self.values[address]}")

def export2excel(self, filename=None) -> "filename.xlsx":

filename = "查询结果.xlsx" if filename is None else filename

pd.Series(self.values, name="查询结果").to_excel(filename)

def export2series(self) -> pd.Series:

return pd.Series(self.values, name="查询结果")

def __getitem__(self, item): return self.values[item]

代码二: Geodistance.py 球面距离计算

from geopy.distance import great_circle, vincenty

# 注意，vincenty 是直线距离，可以直接替换为欧式距离，在 py3.7 中该函数已不可用

def geo_distance(add1: "(latitude, longitude)", add2: "(latitude, longitude)", default=True):

"""提供球面距离(默认)、经纬度距离"""

result = great_circle(add1, add2).kilometers if default else vincenty(add1, add2).kilometers

return result

代码三: data_preprocessing.py 数据预处理

import pandas as pd

from Geocoders import Geocoders

from Geodistance import geo_distance

# 导入数据

data = pd.read_excel('../table/Base_Data.xlsx', index_col='地点')

# 查询经纬度

geo = Geocoders(data.index, api_key=input("Your Baidu's api_key: "))

latitude = geo.export2series().str.split(',').apply(lambda item: item[0]).astype(float)

longitude = geo.export2series().str.split(',').apply(lambda item: item[1]).astype(float)

latitude.name = '纬度'

longitude.name = '经度'

data = data.join(latitude).join(longitude)

# 获取价格矩阵

distance_cost = pd.DataFrame(index=data.index, columns=data.index)

transportation = distance_cost.copy()

for i1 in data.index:

for i2 in data.index:

distance = geo_distance((data.at[i1, '纬度'], data.at[i1, '经度']), (data.at[i2, '纬度'], data.at[i2, '经度']))

if distance / 70.18 > 12:

print(f'{i1}到{i2}太远啦! 坐飞机吧 ~')

transportation.at[i1, i2] = '飞机'

distance_cost.at[i1, i2] = distance * 0.8

else:

transportation.at[i1, i2] = '火车'

distance_cost.at[i1, i2] = distance * 0.5

# 导出数据

distance_cost.to_excel("../table/distance_cost.xlsx")

# 单次会议基本开销

data['基本费用'] = data['住宿费用'] * (data['结束时间'] - data['开始时间'] + 2) + data['注册费用'] + data['会场交通费']

data.to_excel("../table/data.xlsx")

transportation.to_excel("../table/transportation.xlsx")

代码四: shortest_roads.py 获取最小费用候选路径列表, 并导出

矩阵文件

import pandas as pd

data = pd.read_excel('../table/data.xlsx', index_col='地点')

distance_cost = pd.read_excel('../table/distance_cost.xlsx', index_col='地点')

transportation = pd.read_excel('../table/transportation.xlsx', index_col='地点')

# 可达路径

reachability_roads = []

for start in data.index:

for end in data.index:

if data.at[end, '开始时间'] - data.at[start, '结束时间'] >= 2 and start != end:

reachability_roads.append([start, end])

# 构造最小费用邻接路径

shortest_nearest_roads = []

for road in reachability_roads:

start, end = road

if data.at[end, '开始时间'] - data.at[start, '结束时间'] == 2:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, end]}]->{end}", data.at[start, '基本费用'] + distance_cost.at[start, end] + data.at[end, '基本费用']])

else:

# 滞留天数

wait_day = data.at[end, '开始时间'] - data.at[start, '结束时间'] - 2

# start 地驻留

way1_cost = data.at[start, '基本费用'] + data.at[start, '住宿费用'] * wait_day + distance_cost.at[start, end] + data.at[end, '基本费用']

# end 地驻留

way2_cost = data.at[start, '基本费用'] + data.at[end, '住宿费用'] * wait_day + distance_cost.at[start, end] + data.at[end, '基本费用']

# 中转

way3_cost = data.at[start, '基本费用'] + distance_cost.at[start, '上海'] + distance_cost.at['上海', end] + data.at[end, '基本费用']

# 最短邻接路径方案代价

min_cost = min(way1_cost, way2_cost, way3_cost)

# 储存最短邻接路径及代价

if min_cost == way1_cost:

shortest_nearest_roads.append([f"{start}(提前{wait_day}天)[{transportation.at[start, end]}]->{end}", way1_cost])

elif min_cost == way2_cost:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, end]}]->(滞留{wait_day}天){end}", way2_cost])

else:

shortest_nearest_roads.append([f"{start}[{transportation.at[start, '上海']}]->上海(中转{wait_day}天)[{transportation.at['上海', end]}]->{end}", way3_cost])

# 拼接最小费用邻接路径, 以下可以用while True, 迭代创建, 为便于读者理解, 此处拆开

shortest_nearest_roads_2 = shortest_nearest_roads.copy()

shortest_nearest_roads_3 = []

for road1 in shortest_nearest_roads_2:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_3.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads_4 = []

for road1 in shortest_nearest_roads_3:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_4.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads_5 = []

for road1 in shortest_nearest_roads_4:

for road2 in shortest_nearest_roads:

if road1[0][-2:] == road2[0][:2]:

shortest_nearest_roads_5.append([road1[0][:-2] + road2[0], road1[1] + road2[1] - data.at[road1[0][-2:], '基本费用']])

shortest_nearest_roads.extend(shortest_nearest_roads_3 + shortest_nearest_roads_4 + shortest_nearest_roads_5)

# 获取最小费用邻接路径列表

# 先加入所有只去1个地点的路径

shortest_roads_list = []

for location in data.index:

shortest_roads_list.append([f'{location}', data.at[location, '基本费用']])

shortest_roads_list += shortest_nearest_roads

# 再为所有的路径添加出发程、返程路径代价

for road in shortest_roads_list:

road[1] += distance_cost.at['上海', road[0][:2]] + distance_cost.at[road[0][-2:], '上海']

road[0] = f"[出发]上海[{transportation.at['上海', road[0][:2]]}]->{road[0]}[{transportation.at[road[0][-2:], '上海']}]->上海[返回]"

roads = pd.DataFrame({'价格': dict(shortest_roads_list)})

for location in data.index:

# 注意, 这里要先去掉首尾

roads[location] = roads.index.str.split('->').map(lambda items: int(location in ''.join([item for item in items[1:-1] if "中转" not in item])))

roads.index.name = '路径'

roads.to_excel("../table/shortest_roads_list.xlsx")

代码五: Question1.py

代码六: Question2.py

代码七: Question3.py

[超字数了, 翻上面看吧]

3.8 总结

这道题有很多地方不够严谨, 但作为思维训练的题还是很不错的. 从拿到题目后仔细审题, 指出题目中的问题, 并加以分析, 最后化繁就简, 求解模型. 在我的电脑上每一个问题的求解程序都能在2s内求出解果, 也足以说明 Gurobi 的强大. 此外, 使用 Python 与 Gurobi 交互, 可以控制输出结果的显示效果.

以上就是 Gurobi 应用的2个经典案例 (我认为经典), 其中包含了指示变量、哑变量(独热编码)、含变量的数组表达式、多目标优化等等内容, 题目也较为简单, 适合作为小白的赛前训练使用.

参考文献

作者：张柳彬

如有疑问，请联系QQ：965579168

不接论文代做, 转载请声明出处