数学建模竞赛：交通网络优化与可达率提升的线性规划实战

1. 项目概述：从赛题到实战的完整拆解

又到了一年一度的五一数学建模竞赛，今年的B题“未来新城背景下的交通需求规划与可达率问题”一出来，就在我们几个老建模人的小群里炸开了锅。这道题的核心，说白了就是给你一个未来城市的交通网络雏形，让你去规划如何分配有限的交通资源（比如道路容量、公交线路、共享单车点），才能让尽可能多的人，从他们住的地方，能够方便快捷地到达他们想去的地方（比如工作区、商业中心、学校）。这个“尽可能多”的比例，就是题目里反复提到的“可达率”。它不是一个简单的路径搜索问题，而是一个典型的、带有强约束条件的资源优化配置问题。你需要考虑不同交通方式（题目背景里暗示了可能存在多种出行模式）的成本、容量限制，以及不同区域之间动态变化的需求。这题非常适合用数学规划（尤其是线性规划、整数规划）或者网络流模型来啃，Matlab和Python都是绝佳的求解工具。接下来，我就结合自己多年打比赛和带队的经验，把这道题的解题思路、核心模型构建、以及具体的代码实现（Matlab和Python双版本）掰开揉碎了讲清楚，无论你是初次参赛的新手，还是想寻找优化思路的老手，都能从这里找到可以直接“抄作业”的干货。

2. 核心问题解析与建模思路确立

2.1 题目深层需求挖掘

首先，我们得跳出字面意思，理解出题人到底想考察什么。“未来新城”意味着交通网络可能不是现状，而是规划图，这给了我们设定变量和参数的自由度，但也要求模型具备合理性和前瞻性。“交通需求规划”是手段，“可达率”是核心目标。这里隐含了几个关键点：

1. 需求是动态且空间分布的：居民出行需求（O-D矩阵，即从起点Origin到终点Destination的流量）不是均匀的，通常早晚高峰特征明显，且居住区、就业区、商业区的分布会导致需求矩阵非常稀疏（即不是所有区域之间都有大量出行）。
2. 资源是有限且有成本的：无论是修建道路、开设公交线路还是投放共享单车，都有预算或物理限制。提高可达率不能靠无限堆资源，必须在成本和效益之间做权衡。
3. 可达性的定义：这是模型的灵魂。是单纯的地理上连通就算可达？还是必须在某个时间阈值内（比如30分钟）到达才算？或者是考虑出行成本（时间、费用）低于某个门槛？题目没有明说，这就需要我们根据“未来新城”的背景赋予其一个合理的、可量化的定义。通常，我们会定义一个“效用函数”或“广义成本”，当成本低于阈值时，认为该O-D对是可达的。

2.2 模型框架选择：为什么是线性/整数规划？

面对这种“在约束下最大化或最小化某个目标”的问题，数学规划是首选。具体到本题：

• 线性规划：如果你的模型假设所有变量都是连续的，并且目标函数和约束条件都是线性的，那么LP是最高效的。例如，将每条路径上的流量视为连续变量，目标是最小化总出行时间或最大化满足需求的流量（即可达率）。它的优势是求解速度快，理论成熟。但缺点也很明显：它无法处理“是否开设某条公交线路”这种0-1决策问题。
• 整数规划/混合整数线性规划：这才是更贴合本题复杂性的框架。你需要引入0-1决策变量来表示是否在某个候选走廊上建设轨道交通、是否在某个小区设置共享单车枢纽。同时，流量分配可以用连续变量。MILP能够完美整合网络设计（离散决策）和流量分配（连续决策）问题。虽然求解比LP困难，但对于竞赛规模的题目，现代求解器（如Gurobi, CPLEX，或开源的OR-Tools、SCIP）完全能应付。

我的思路选择：我倾向于构建一个两阶段模型或集成MILP模型。第一阶段（或模型的一部分）用整数变量决定交通基础设施的布局（哪里建站，哪里开线）；第二阶段（或模型的另一部分）在确定的网络基础上，用线性规划或用户均衡模型分配流量，计算可达率。这样结构清晰，也便于分步调试。

2.3 关键参数与变量定义

在写代码之前，必须把模型“翻译”成数学语言。以下是一个核心的变量与参数集合，你可以根据题目具体数据调整：

• 集合：

  • N: 网络中的所有节点集合（如交通小区、交叉口）。
  • A: 有向弧集合（即路段，(i,j)表示从节点i到节点j）。
  • K: 出行模式集合（如步行、自行车、公交、地铁、自驾）。
  • W: O-D对集合（(r,s)表示从起点r到终点s的需求）。

• 参数：

  • d_rs: O-D对(r,s)之间的出行需求量（人/小时）。
  • c_ij^k: 使用模式k通过弧(i,j)的广义成本（可能是时间、费用或混合）。
  • u_ij: 弧(i,j)的容量限制（车辆/小时）。
  • B: 总预算限制。
  • cost_ij^k: 在弧(i,j)上提供模式k服务的单位建设或运营成本。
  • M: 一个足够大的常数（Big-M法中用）。

• 决策变量：

  • x_ij^k, rs: 连续变量，表示O-D对(r,s)之间使用模式k在弧(i,j)上的流量。
  • y_ij^k: 0-1变量，表示是否在弧(i,j)上提供交通模式k的服务（网络设计决策）。
  • z_rs: 0-1变量，表示O-D对(r,s)的需求是否被满足（即可达）。

• 目标函数：最大化总可达率，即Maximize Σ_(r,s) d_rs * z_rs。或者，最大化被满足的需求量。

• 核心约束：

  1. 流量守恒约束：对于每个O-D对(r,s)，每个节点i，每个模式k，流入等于流出（除非是起点或终点）。
  2. 容量约束：对于每个弧(i,j)，所有O-D对和模式的总流量不能超过该弧容量：Σ_(r,s,k) x_ij^k, rs <= u_ij。
  3. 模式可用性约束（逻辑约束）：如果y_ij^k = 0（不提供该模式服务），则对应所有O-D的流量必须为0。这需要用Big-M法线性化：x_ij^k, rs <= M * y_ij^k。
  4. 预算约束：Σ_(i,j,k) cost_ij^k * y_ij^k <= B。
  5. 可达性判定约束：这是连接流量分配和可达率的关键。一种常见方法是，定义一个O-D对(r,s)的可达成本C_rs（例如，其最短路成本或均衡成本），当C_rs <= Threshold时，z_rs = 1。这同样需要线性化处理，可能引入辅助变量和Big-M约束。

注意：这里的模型是一个高度简化的框架。实际比赛中，你需要根据题目给出的具体数据表（如节点坐标、需求矩阵、路段基础属性、模式特性成本等）来细化参数和约束。例如，可能还需要考虑不同模式间的换乘成本、需求的时间分布等。

3. 求解工具链与代码环境搭建

3.1 Matlab vs. Python：如何选择？

两者都能出色完成任务，选择取决于你的团队熟悉度和题目细微需求。

• Matlab：

  • 优势：优化工具箱（optimproblem）语法直观，特别适合表述线性规划和整数规划问题。内置的linprog,intlinprog求解器对于中小规模问题非常方便。矩阵运算和绘图功能强大，快速验证模型时很高效。
  • 劣势：处理超大规模整数规划可能不如专业商业求解器快。许可证问题（虽然比赛通常提供）。
  • 适合人群：习惯Matlab矩阵操作、追求快速原型验证、问题规模适中的队伍。

• Python：

  • 优势：生态丰富。你可以用PuLP、CVXPY（用于凸优化）进行建模，然后调用如CBC（开源）、Gurobi、CPLEX（商业，但学术免费）等强大求解器。Pandas处理表格数据（如O-D矩阵）比Matlab更灵活。NetworkX用于复杂网络分析。代码更易集成和版本管理。
  • 劣势：环境配置稍复杂，需要安装多个库。
  • 适合人群：熟悉Python数据处理、需要处理复杂数据或大规模问题、希望代码更具通用性的队伍。

我的建议：优先使用Python + Gurobi。Gurobi的学术许可很容易申请，其求解MILP的性能是顶级的，而且Python的gurobipy接口建模非常自然。如果时间紧或对Python不熟，Matlab的intlinprog是可靠的保底选择。

3.2 环境配置要点

Python环境（以Anaconda为例）：

1. 创建独立环境：conda create -n math_modeling python=3.9
2. 激活环境：conda activate math_modeling
3. 安装核心库：

   pip install numpy pandas scipy matplotlib jupyter pip install pulp # 可选，轻量级建模 pip install gurobipy # 首选！记得去Gurobi官网申请学术许可证并安装 pip install networkx # 用于网络分析

Matlab环境：确保已安装Optimization Toolbox。在命令窗口输入ver查看。建模主要使用optimproblem对象。

实操心得：比赛开始前，务必在你们的电脑上完整跑通一个简单的线性规划示例，比如经典的运输问题。这能确保所有环境、许可证都工作正常，避免赛时浪费宝贵时间在配置上。把示例代码存成模板。

4. 数据预处理与网络构建实战

题目通常会提供节点表、边表、需求表。这部分是模型的基础，处理不好，后面全盘皆输。

4.1 数据读取与清洗

Python (Pandas) 示例：

import pandas as pd import numpy as np # 假设有CSV文件 nodes_df = pd.read_csv('nodes.csv') # 列可能包括：node_id, x, y, type(居住/就业/商业) edges_df = pd.read_csv('edges.csv') # 列可能包括：from_node, to_node, length, base_capacity, speed_limit demand_df = pd.read_csv('demand.csv') # 列可能包括：origin, destination, demand_peak, demand_offpeak # 数据清洗：检查缺失值、重复边、非数值数据 print(nodes_df.info()) print(edges_df.isnull().sum()) # 处理缺失值，例如用均值填充或删除 edges_df['capacity'].fillna(edges_df['capacity'].median(), inplace=True) # 构建节点映射字典，将节点ID映射到连续的索引（0,1,2,...），这对后续矩阵运算至关重要 node_id_to_index = {id: idx for idx, id in enumerate(nodes_df['node_id'].unique())} num_nodes = len(node_id_to_index)

Matlab 示例：

% 读取数据 nodes = readtable('nodes.csv'); edges = readtable('edges.csv'); demand = readtable('demand.csv'); % 获取基本维度 numNodes = height(nodes); numEdges = height(edges); % 同样建议构建映射 [uniqueNodeIDs, ~, nodeIdx] = unique(nodes.node_id); % nodeIdx 现在可以用于索引

4.2 网络拓扑与成本矩阵生成

我们需要计算任意两点间的最短路径成本（作为可达性判定的基础或初始解）。这里以时间成本为例，假设成本与长度成正比，与速度成反比。

Python (NetworkX) 示例：

import networkx as nx G = nx.DiGraph() # 创建有向图 # 添加带权重的边 for _, row in edges_df.iterrows(): from_idx = node_id_to_index[row['from_node']] to_idx = node_id_to_index[row['to_node']] # 计算时间成本（小时），假设速度单位是km/h，长度是km travel_time = row['length'] / row['speed_limit'] G.add_edge(from_idx, to_idx, weight=travel_time, capacity=row['base_capacity']) # 计算所有节点对的最短路径长度（成本） # 注意：对于大规模网络，所有节点对最短路径计算量很大(O(n^3))，谨慎使用。 # 通常只计算有需求的O-D对之间的最短路径。 predecessors, shortest_dist = nx.floyd_warshall_predecessor_and_distance(G, weight='weight') # shortest_dist[i][j] 就是节点i到j的最短时间成本 # 更实用的方法：只为存在需求的O-D对计算 od_pairs = list(zip(demand_df['origin'].map(node_id_to_index), demand_df['destination'].map(node_id_to_index))) od_costs = {} for o, d in od_pairs: try: length = nx.shortest_path_length(G, source=o, target=d, weight='weight') od_costs[(o, d)] = length except nx.NetworkXNoPath: od_costs[(o, d)] = np.inf # 不可达

Matlab 示例：Matlab没有现成的复杂图论库像NetworkX，但可以借助稀疏矩阵和最短路径算法。

% 构建邻接矩阵（成本矩阵） costMatrix = inf(numNodes, numNodes); for i = 1:numEdges fromIdx = find(uniqueNodeIDs == edges.from_node(i)); toIdx = find(uniqueNodeIDs == edges.to_node(i)); travelTime = edges.length(i) / edges.speed_limit(i); costMatrix(fromIdx, toIdx) = travelTime; end % 对角线设为0 for i = 1:numNodes costMatrix(i, i) = 0; end % 使用Floyd-Warshall算法计算所有点对最短路径（适用于节点数不多时） for k = 1:numNodes for i = 1:numNodes for j = 1:numNodes if costMatrix(i, k) + costMatrix(k, j) < costMatrix(i, j) costMatrix(i, j) = costMatrix(i, k) + costMatrix(k, j); end end end end % 现在 costMatrix(i,j) 就是节点i到j的最短时间成本

注意事项：网络规模很大时（节点>500），所有点对最短路径计算会成为瓶颈。务必根据题目数据规模选择算法。对于稀疏网络和特定O-D对，使用Dijkstra算法（nx.single_source_dijkstra_path_lengthin Python,graphshortestpathin Matlab）多次计算更高效。赛题数据通常会在合理范围内。

5. 数学规划模型构建与代码实现

这是最核心的部分。我将以一个简化的、但包含网络设计（0-1变量）和流量分配（连续变量）的MILP模型为例，分别用Python+Gurobi和Matlab实现。

5.1 模型假设与简化

为了代码清晰，我们做以下简化：

1. 只考虑一种出行模式（如小汽车）。
2. “可达”定义为：该O-D对间被分配到的流量至少达到其需求的某个比例（例如80%），或者其最短路时间成本低于阈值。
3. 决策变量y_ij表示是否提升该路段的容量（代表一种投资）。
4. 目标：在投资预算B下，最大化被满足的O-D需求总量。

5.2 Python + Gurobi 实现

import gurobipy as gp from gurobipy import GRB def solve_transportation_planning(num_nodes, edges_list, demand_list, budget, base_capacity, upgrade_cost, upgrade_capacity_gain, bigM=1e6): """ 求解交通规划MILP模型 edges_list: 列表，每个元素为 (from_idx, to_idx, travel_time, base_cap) demand_list: 列表，每个元素为 (origin_idx, dest_idx, demand_value) """ model = gp.Model("FutureCityTransportation") # --- 创建变量 --- # 流量变量 x_{ij}^{rs}， 这里简化了，先假设为总流量 x_{ij}， 实际需要更精细的守恒约束 # 更准确的建模需要为每个O-D对设置路径流量或弧流量，这里用弧流量近似 x = {} # 弧上的总流量 for (i, j, _, _) in edges_list: x[i, j] = model.addVar(lb=0.0, vtype=GRB.CONTINUOUS, name=f"x_{i}_{j}") # 网络设计变量 y_{ij}: 是否升级弧(i,j) y = {} for (i, j, _, _) in edges_list: y[i, j] = model.addVar(vtype=GRB.BINARY, name=f"y_{i}_{j}") # 可达性变量 z_{rs}: O-D对(r,s)的需求是否被满足（流量>=需求比例） z = {} for (r, s, d) in demand_list: z[r, s] = model.addVar(vtype=GRB.BINARY, name=f"z_{r}_{s}") # --- 设置目标函数：最大化满足的需求总量 --- obj_expr = gp.quicksum(d * z[r, s] for (r, s, d) in demand_list) model.setObjective(obj_expr, GRB.MAXIMIZE) # --- 添加约束 --- # 1. 流量守恒约束（简化版，实际应对每个节点，每个O-D对） # 这里用一个高度简化的约束：所有流入某个目的地的流量 >= 该目的地的总需求 * 可达率？ 这并不准确。 # 正确的做法是构建节点-弧关联矩阵，为每个O-D对设置流量平衡。篇幅所限，这里展示更复杂的模型框架。 # 我们假设已经通过最短路或其它方式得到了一个候选路径集P_rs，决策变量为路径流量 f_p。 # 以下代码块示意一个更接近真实但复杂的建模思路： print("构建精确的节点-弧流量平衡约束...") # 首先，为每个O-D对，创建从起点到终点的流变量，并满足平衡 x_rs = {} # 四级字典: x_rs[r][s][i][j] 或使用GRB的tuplelist # 此处省略详细的、每个O-D对每个节点的流量平衡约束代码，它需要大量循环。 # 取而代之，我们用一个聚合的、示意性的容量约束代替，以保持示例可运行。 # 2. 弧容量约束 (与升级决策相关) for (i, j, _, base_cap) in edges_list: # 升级后的容量 = 基础容量 + 升级增益 * 升级决策 model.addConstr(x.get((i,j), 0) <= base_cap + upgrade_capacity_gain * y[i, j], name=f"cap_{i}_{j}") # 3. 预算约束 model.addConstr(gp.quicksum(upgrade_cost * y[i, j] for (i, j, _, _) in edges_list) <= budget, name="budget") # 4. 可达性逻辑约束 (简化：如果分配到目的节点s的总流量 >= 总需求的某个比例，则z_rs=1) # 首先计算每个目的地s的总需求 total_demand_to_s = {} for (r, s, d) in demand_list: total_demand_to_s[s] = total_demand_to_s.get(s, 0) + d # 对于每个目的地s，计算流入它的总流量（简化，假设所有弧都指向它） # 这非常不精确，仅用于示意。实际必须用节点净流量=需求来约束。 inflow_to_s = {} for s in set([s for (_, s, _) in demand_list]): inflow_expr = gp.quicksum(x[i, j] for (i, j, _, _) in edges_list if j == s) # 假设j是终点 # 逻辑约束：如果 inflow_expr >= 0.8 * total_demand_to_s[s]，则 z_rs 可以取1（对于所有以s为终点的O-D对） # 这里需要为每个(r,s)对单独约束，再次简化 for (r, s_d, d) in demand_list: if s_d == s: # Big-M 线性化: inflow_expr - 0.8*total_demand_to_s[s] >= -M*(1-z[r,s]) # 同时 inflow_expr - 0.8*total_demand_to_s[s] <= M*z[r,s] - epsilon (迫使满足时z必须为1) # 篇幅原因，此处省略具体实现。这是一个关键且容易出错的点。 pass # 由于上述精确约束代码冗长，我们临时修改目标为：最大化总流量，同时受预算和容量限制 # 这变成了一个简单的网络流问题，以演示Gurobi求解。 print("注意：为示例可运行，模型已简化为最大化总流量，未体现精确的O-D对可达性约束。") model.setObjective(gp.quicksum(x[i,j] for (i,j) in x.keys()), GRB.MAXIMIZE) # --- 求解模型 --- model.optimize() # --- 输出结果 --- if model.status == GRB.OPTIMAL: print(f"最优目标值（总流量）: {model.objVal:.2f}") print("升级的弧:") for (i, j, _, _) in edges_list: if y[i, j].X > 0.5: print(f" 弧 ({i} -> {j})") # 输出流量较大的弧 print("\n流量较大的弧:") for (i, j) in x.keys(): if x[i, j].X > 1.0: print(f" 弧 ({i} -> {j}): 流量 = {x[i, j].X:.2f}") else: print("未找到最优解") return model, x, y # 假设一些数据 edges = [(0,1,1.2,100), (1,2,0.8,80), (0,2,2.0,150), (2,3,1.5,90), (1,3,1.0,120)] demands = [(0,3,50), (1,3,30)] BUDGET = 500 UPGRADE_COST = 100 CAPACITY_GAIN = 50 model, x_vars, y_vars = solve_transportation_planning( num_nodes=4, edges_list=edges, demand_list=demands, budget=BUDGET, base_capacity=None, # edges中已包含 upgrade_cost=UPGRADE_COST, upgrade_capacity_gain=CAPACITY_GAIN )

5.3 Matlab 优化工具箱实现

Matlab中使用optimproblem对象来构建模型，更贴近数学表达。

% 假设数据 edges = [0 1 1.2 100; 1 2 0.8 80; 0 2 2.0 150; 2 3 1.5 90; 1 3 1.0 120]; % 列：from, to, time, base_cap demands = [0 3 50; 1 3 30]; % 列：origin, dest, demand budget = 500; upgradeCost = 100; capacityGain = 50; bigM = 1e6; numEdges = size(edges, 1); numDemands = size(demands, 1); nodeList = unique([edges(:,1); edges(:,2)]); numNodes = length(nodeList); % 创建优化问题 prob = optimproblem('Description', 'Future City Transport MILP'); % 1. 定义变量 % 流量变量 (连续，非负) x = optimvar('x', numEdges, 'LowerBound', 0); % 升级决策变量 (0-1) y = optimvar('y', numEdges, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); % 可达性变量 (0-1) - 这里简化，假设一个总体的“需求满足率”变量，实际应对应每个O-D对 % z = optimvar('z', numDemands, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); % 2. 目标函数：最大化总流量 (简化版) prob.Objective = sum(x); % 更真实的目标：最大化满足的需求 sum(demands(:,3) .* z) % 3. 约束条件 % 容量约束 (每条边的流量 <= 基础容量 + 升级增益 * 升级决策) capacityConstraints = optimconstr(numEdges); for e = 1:numEdges capacityConstraints(e) = x(e) <= edges(e, 4) + capacityGain * y(e); end prob.Constraints.capacityConstraints = capacityConstraints; % 预算约束 prob.Constraints.budgetConstraint = sum(upgradeCost * y) <= budget; % 流量平衡约束 (高度简化版) % 实际应为：对于每个节点，流入量 - 流出量 = 净需求（起点为+需求，终点为-需求，中间点为0） % 这里用一个简单的“总流出>=总需求”约束代替，仅作演示 totalDemand = sum(demands(:,3)); % 找出所有以需求终点为结束的边 destNode = 3; % 假设所有需求都去节点3 edgesToDest = find(edges(:,2) == destNode); if ~isempty(edgesToDest) prob.Constraints.demandSatisfaction = sum(x(edgesToDest)) >= totalDemand * 0.8; % 至少满足80%需求 end % 4. 求解 [sol, fval, exitflag, output] = solve(prob); % 5. 输出结果 if exitflag > 0 fprintf('最优总流量: %.2f\n', fval); fprintf('升级的边:\n'); for e = 1:numEdges if sol.y(e) > 0.5 fprintf(' 边 %d (%d -> %d)\n', e, edges(e,1), edges(e,2)); end end fprintf('\n各边流量:\n'); for e = 1:numEdges if sol.x(e) > 0 fprintf(' 边 %d (%d -> %d): 流量 = %.2f\n', e, edges(e,1), edges(e,2), sol.x(e)); end end else fprintf('求解失败。\n'); disp(output); end

关键提示：以上两个代码示例中的流量平衡约束都被大幅简化了。在实际的交通分配或网络流问题中，必须为每个节点和每个O-D对建立精确的流量守恒方程。这是建模中最容易出错的地方。完整的实现需要构建节点-弧关联矩阵，并为每个O-D对定义一组流变量。由于代码非常冗长，上述示例旨在展示模型框架和求解器调用方法。你需要根据题目具体的数据结构，补全这部分核心逻辑。

6. 结果可视化与方案分析

模型求解后，输出一堆数字是不够的，需要用图表让结果说话。

6.1 网络拓扑与流量热力图

Python (Matplotlib + NetworkX) 示例：

import matplotlib.pyplot as plt plt.figure(figsize=(12, 8)) # 1. 绘制基础网络 pos = nx.spring_layout(G, seed=42) # 节点位置 nx.draw_networkx_nodes(G, pos, node_color='lightblue', node_size=300) nx.draw_networkx_labels(G, pos, font_size=10) # 2. 绘制边，宽度和颜色代表流量 edge_widths = [sol_x.get((i, j), 0) / 10 for (i, j) in G.edges()] # 缩放流量便于显示 edge_colors = [sol_x.get((i, j), 0) for (i, j) in G.edges()] edges = nx.draw_networkx_edges(G, pos, width=edge_widths, edge_color=edge_colors, edge_cmap=plt.cm.Blues, alpha=0.7) # 添加颜色条 pc = plt.cm.ScalarMappable(cmap=plt.cm.Blues, norm=plt.Normalize(vmin=min(edge_colors), vmax=max(edge_colors))) pc.set_array(edge_colors) plt.colorbar(pc, label='Traffic Flow') # 3. 高亮显示被升级的边 upgraded_edges = [(i, j) for (i, j) in G.edges() if sol_y.get((i, j), 0) > 0.5] nx.draw_networkx_edges(G, pos, edgelist=upgraded_edges, width=3, edge_color='red', style='dashed', label='Upgraded Links') plt.title("Future City Traffic Flow & Upgrade Plan") plt.axis('off') plt.legend() plt.tight_layout() plt.show() # 4. 绘制可达率分析柱状图 # 假设我们计算了每个O-D对的可达状态 od_pairs = list(od_costs.keys()) od_labels = [f'{o}-{d}' for (o, d) in od_pairs] reachable = [1 if od_costs[pair] <= TIME_THRESHOLD else 0 for pair in od_pairs] # 假设有个时间阈值 plt.figure(figsize=(10, 6)) plt.bar(od_labels, reachable, color=['green' if r else 'red' for r in reachable]) plt.xlabel('O-D Pair') plt.ylabel('Reachable (1) / Unreachable (0)') plt.title('Accessibility Analysis for Each O-D Pair') plt.xticks(rotation=45, ha='right') plt.tight_layout() plt.show()

6.2 灵敏度分析与方案对比

求解一次得到的是“最优解”，但作为优秀论文，还需要分析这个解的稳健性和参数敏感性。

• 预算灵敏度：逐渐增加或减少预算B，观察可达率的变化。绘制“预算-可达率”曲线。这能告诉决策者，投入的边际效益在哪里开始递减。

  budget_range = np.linspace(100, 1000, 10) reachability_rates = [] for b in budget_range: # 修改模型中的预算约束，重新求解 # 注意：每次求解需重新构建模型或修改约束，耗时。对于MILP，可以考虑用参数化求解或多次调用。 # 这里示意流程 model.getConstrByName('budget').RHS = b model.optimize() if model.status == GRB.OPTIMAL: # 计算当前解的可达率 total_demand = sum(d for _,_,d in demands) satisfied_demand = sum(d * sol_z[r,s].X for (r,s,d) in demands) # 假设有z变量 reachability_rates.append(satisfied_demand / total_demand) else: reachability_rates.append(0) plt.plot(budget_range, reachability_rates, 'bo-') plt.xlabel('Total Budget') plt.ylabel('Overall Reachability Rate') plt.grid(True) plt.show()

• 需求波动分析：将需求数据d_rs按一定比例（如±20%）随机扰动，多次运行模型，观察最优解（哪些边被升级）是否稳定。如果解变化很大，说明方案对需求预测很敏感，需要在论文中提出弹性规划建议。
• 方案对比：设计不同的基准方案，例如：
  1. 无投资方案：仅使用现有网络。
  2. 均匀投资方案：预算平均分配到所有边。
  3. 启发式方案：优先升级最拥堵的边（基于初始流量分配）。 将你的优化方案与这些基准方案在可达率、平均出行时间等指标上进行对比，突出优化模型的价值。

7. 论文写作要点与避坑指南

模型和代码搞定了，最后要把你的工作清晰地展现在论文里。

7.1 论文结构建议

1. 摘要：用一段话概括问题、你的方法、模型、算法、主要结果和结论。避免细节，突出亮点。
2. 问题重述与分析：用自己的话梳理题目，明确假设，定义关键指标（可达率）。
3. 模型假设与符号说明：列出所有假设（如需求固定、单模式等），用表格清晰列出所有符号。
4. 模型建立：这是核心。
   • 先给出整体建模思路框图。
   • 然后分小节详细阐述：网络构建、可达性定义、目标函数、约束条件（流量守恒、容量、预算、逻辑约束）。务必解释每个约束的物理或经济意义。
   • 给出完整的数学模型（公式）。
5. 模型求解与算法设计：
   • 说明你将复杂的MILP模型如何转化为求解器可解的形式（如线性化技巧）。
   • 描述求解流程（数据预处理→模型构建→调用求解器→结果提取）。
   • 如果是大规模问题，你设计了什么启发式或分解算法？要说明。
6. 结果分析与可视化：
   • 展示关键结果（最优可达率、投资方案、流量分布）。
   • 用图表说话（网络图、柱状图、灵敏度分析图）。
   • 对结果进行深入分析：为什么这些边被选中升级？方案有何特点？
7. 灵敏度分析与模型检验：展示预算、需求等参数变化时结果如何变化，检验模型稳定性。
8. 模型评价与推广：客观评价模型的优点（如综合考虑成本与效益）和缺点（如假设需求固定），并提出改进方向（如加入动态交通分配、多模式竞争）。

7.2 常见踩坑点与应对策略

• 坑1：模型求解时间爆炸。MILP问题规模稍大就可能算很久。
  • 对策：合理简化网络（聚合小区）；使用启发式算法先求一个较好解作为MILP的初始解；设置求解时间限制；对于超大规模问题，考虑分解算法（如Benders分解）。
• 坑2：流量不守恒，出现“幽灵流量”。这是建模中最常见的错误。
  • 对策：在添加流量平衡约束后，用一个小型网络（3-4个节点）手动计算验证。确保对于每个O-D对，起点只有流出，终点只有流入，中间点净流量为零。
• 坑3：结果不直观或违反常识。比如投资全部集中在一条偏僻的路上。
  • 对策：检查目标函数和约束是否准确反映了现实。例如，是否忽略了出行时间成本？可达性定义是否合理？回去审视你的模型假设。
• 坑4：代码调试困难。
  • 对策：模块化编程。将数据读取、网络构建、模型定义、求解、后处理分成独立函数。先在小数据集上测试通过。多用print或disp语句输出中间变量值。利用求解器的输出信息（如Gurobi的model.write('model.lp')可以输出模型文件检查）。
• 坑5：论文像代码说明书。
  • 对策：论文的重点是思路、模型和结论，不是代码。代码可以放附录。在正文中用流程图、公式和文字描述你的算法。图表要有标题和注释，做到“不言自明”。

最后，记住数学建模竞赛是“建模”而不是“编程”比赛。清晰的逻辑、合理的假设、创新的视角、严谨的论证和美观的表达，比单纯的代码技巧更重要。把上面提供的思路和代码作为你的脚手架，结合题目具体数据，构建出属于你们队伍的、有说服力的解决方案。祝大家在五一赛中取得好成绩！