【2025妈妈杯MathorcupD题】详细解题思路代码分享

Better Rose

已于 2025-04-18 10:32:49 修改

阅读量1.8k

点赞数 34

分类专栏：数学建模文章标签：数学建模妈妈杯

于 2025-04-18 10:31:05 首次发布

本文链接：https://blog.csdn.net/Magnolia_He/article/details/147321867

版权

数学建模专栏收录该内容

108 篇文章

订阅专栏

2025妈妈杯Mathorcup D题

详细解题思路代码分享

在这里插入图片描述

一、问题一

1.1 问题分析

短途运输货量预测需综合考虑历史货量规律、预知货量数据的修正以及突发因素影响。根据附件2历史货量数据和附件3预知货量数据，需解决以下核心问题：

数据修正：预知数据存在发运异常、未下单货量及取消订单的偏差，需通过历史数据对比修正。
时间序列建模：短途运输货量具有日周期性（如发运节点6点和14点）和趋势性（如周末高峰）。
颗粒度拆解：将总货量拆分为10分钟颗粒度，需结合货量到达规律（如包裹集中到达分拣中心的时间分布）。

1.2 数学模型

货量修正模型

设预知货量为 $\ P_{\text{pre}}$ ，实际货量为 $\ P_{\text{real}}$ ，修正系数 $\ \alpha$ 通过历史数据回归得到：

$P_{\text{real}} = \alpha \cdot P_{\text{pre}} + \epsilon$

其中 $\ \epsilon$ 为随机误差项，通过正态分布建模。

时间序列预测模型（ARIMA）

采用季节性差分自回归移动平均模型 SARIMA(p, d, q)(P, D, Q)，参数通过网格搜索确定。模型表达式为：

$\left(1 - \sum_{i=1}^{p} \phi_i L^i\right)\left(1 - \sum_{s=1}^{P} \Phi_s L^{s}\right) \nabla^d \nabla_s^D P_t = \left(1 + \sum_{i=1}^{q} \theta_i L^i\right)\left(1 + \sum_{s=1}^{Q} \Theta_s L^{is}\right) \epsilon_t$

其中 $\ L$ 为滞后算子， $\ s=24$ （日周期性）， $\ \nabla$ 为差分操作。

颗粒度拆解模型

基于历史货量分布，定义时间权重系数 $W_t$ ：

$P_{10\text{ min}} = P_{\text{total}} \cdot W_t$

$W_t$ 通过核密度估计（KDE）拟合历史数据分布得到。

1.3 Python代码

import pandas as pd
import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.neighbors import KernelDensity
from scipy.stats import norm

# 模拟数据生成（示例）
date_rng = pd.date_range(start='2023-01-01', end='2023-01-30', freq='10T')
historical_data = pd.DataFrame({
    'timestamp': date_rng,
    '货量': np.random.poisson(50, len(date_rng)) + np.sin(np.linspace(0, 8*np.pi, len(date_rng)))*20
})

# SARIMA模型训练
def sarima_forecast(data):
    model = SARIMAX(data, 
                    order=(1,1,1), 
                    seasonal_order=(1,1,1,144))  # 144=24*6 (10分钟颗粒度)
    results = model.fit(disp=False)
    forecast = results.get_forecast(steps=6*24)  # 预测24小时
    return forecast.predicted_mean

# 颗粒度拆解
def time_decomposition(total_forecast, historical):
    # 提取小时和分钟作为特征
    historical['hour'] = historical['timestamp'].dt.hour + historical['timestamp'].dt.minute/60
    # KDE拟合
    kde = KernelDensity(bandwidth=0.5).fit(historical['hour'].values.reshape(-1,1))
    
    # 生成10分钟间隔的时间点
    time_points = np.linspace(0, 24, 24*6).reshape(-1,1)
    log_density = kde.score_samples(time_points)
    weights = np.exp(log_density)
    weights /= weights.sum()
    
    return total_forecast * weights

# 主程序
if __name__ == "__main__":
    # 训练预测模型
    ts_data = historical_data.set_index('timestamp')['货量']
    forecast_total = sarima_forecast(ts_data)
    
    # 时间分解
    decomposed = time_decomposition(forecast_total.sum(), historical_data)
    
    # 输出结果
    result_df = pd.DataFrame({
        '时间戳': pd.date_range(start='2023-01-31', periods=144, freq='10T'),
        '预测货量': decomposed
    })
    result_df.to_csv('problem1_forecast.csv', index=False)

1.4 Matlab代码

% 数据生成（示例）
t = datetime(2023,1,1):minutes(10):datetime(2023,1,30);
y = poissrnd(50,1,length(t)) + 20*sin(linspace(0,8*pi,length(t)))';
data = timetable(t', y');

% SARIMA模型
model = arima('Constant',0, 'D',1, 'Seasonality',144,...
    'MALags',1, 'SMALags',1);
fit = estimate(model, data.y);
[forecast, ~] = forecast(fit, 144); % 预测24小时

% 时间分解
function weights = kde_decomposition(data)
    % 提取小时信息
    hours = hour(data.t) + minute(data.t)/60;
    [f, xi] = ksdensity(hours, linspace(0,24,144), 'Bandwidth',0.5);
    weights = f / sum(f);
end

% 主程序
weights = kde_decomposition(data);
decomposed = sum(forecast) * weights';

% 结果输出
timestamps = datetime(2023,1,31):minutes(10):datetime(2023,1,31,23,50);
result = table(timestamps', decomposed, 'VariableNames', {'Time','Forecast'});
writetable(result, 'problem1_forecast.csv');

二、问题二

2.1 问题分析及数学模型

需将货量转化为车次需求，并优化自有车辆调度。

车次需求计算

车次数 $\ N_{\text{truck}}$ 由累计货量 $\ P_{\text{cum}}$ 和单车货量C决定：

$N_{\text{truck}} = \left\lceil \frac{P_{\text{cum}}}{C} \right\rceil$

车辆调度优化模型（混合整数规划）

目标函数：

$\min(\lambda_1 \cdot \text{外部成本} + \lambda_2 \cdot (1 / \text{周转率}))$

2.2 Python代码

from ortools.sat.python import cp_model
import pandas as pd

class Task:
    def __init__(self, task_id, start, end, demand):
        self.id = task_id
        self.start_min = start  # 分钟数表示的时间
        self.end_min = end
        self.demand = demand

def vehicle_scheduling(tasks, vehicle_capacity=1000, max_merge=3):
    model = cp_model.CpModel()
    horizon = 24 * 60  # 全天分钟数
    
    # 创建任务变量
    intervals = []
    for task in tasks:
        start = model.NewIntVar(0, horizon, f'start_{task.id}')
        duration = model.NewIntVar(1, horizon, f'duration_{task.id}')
        end = model.NewIntVar(0, horizon, f'end_{task.id}')
        model.Add(start >= task.start_min)
        model.Add(end <= task.end_min)
        model.Add(end == start + duration)
        intervals.append((start, duration, end))
    
    # 车辆数量约束
    vehicle_used = []
    for i in range(len(tasks)):
        vehicle = model.NewBoolVar(f'vehicle_{i}')
        vehicle_used.append(vehicle)
    
    # 任务分配和合并约束
    for i in range(len(tasks)):
        for j in range(i+1, len(tasks)):
            # 合并任务数不超过max_merge
            merge = model.NewBoolVar(f'merge_{i}_{j}')
            model.Add(intervals[i][2] <= intervals[j][0]).OnlyEnforceIf(merge)
            model.Add(intervals[j][2] - intervals[i][0] <= 60 * max_merge).OnlyEnforceIf(merge)
    
    # 目标：最小化外部车辆使用
    model.Minimize(sum(vehicle_used))
    
    solver = cp_model.CpSolver()
    status = solver.Solve(model)
    
    if status == cp_model.OPTIMAL:
        return [(solver.Value(s), solver.Value(e)) for s,_,e in intervals]
    else:
        return None

# 示例任务生成
tasks = [
    Task(0, 360, 600, 800),   # 6:00-10:00
    Task(1, 720, 960, 1200),  # 12:00-16:00
    Task(2, 1080, 1320, 600)  # 18:00-22:00
]

schedule = vehicle_scheduling(tasks)
print(schedule)

2.3 Matlab代码

% 车辆调度问题（需要Optimization Toolbox）
function schedule = vehicle_scheduling(tasks)
    % tasks结构体数组：start_min, end_min, demand
    num_tasks = length(tasks);
    horizon = 24*60;
    
    prob = optimproblem('ObjectiveSense','min');
    
    % 变量定义
    start = optimvar('start', num_tasks, 'Type','integer','LowerBound',0,'UpperBound',horizon);
    duration = optimvar('duration', num_tasks, 'Type','integer','LowerBound',1,'UpperBound',horizon);
    vehicle = optimvar('vehicle', num_tasks, 'Type','integer','LowerBound',0,'UpperBound',1);
    
    % 约束
    for i = 1:num_tasks
        prob.Constraints.(['task_' num2str(i)]) = start(i) >= tasks(i).start_min;
        prob.Constraints.(['end_' num2str(i)]) = start(i) + duration(i) <= tasks(i).end_min;
    end
    
    % 车辆使用约束
    prob.Objective = sum(vehicle);
    
    % 求解
    options = optimoptions('intlinprog','Display','off');
    [sol,~,exitflag] = solve(prob,'Options',options);
    
    if exitflag == 1
        schedule = [sol.start, sol.start + sol.duration];
    else
        schedule = [];
    end
end

% 示例任务
tasks = struct(...
    'start_min', {360, 720, 1080},...
    'end_min', {600, 960, 1320},...
    'demand', {800, 1200, 600});
result = vehicle_scheduling(tasks);
disp(result);

三、问题三

3.1 问题分析及数学模型

调整模型：

装载量调整：( C = 800 )，装卸时间缩短至 10 分钟。
目标函数更新：增加容器使用决策变量 $\ u_i \in \{0,1\}$ 。

调整编码：

编码增加是否选择标准化容器的二元决策变量。

3.2 Python代码

# 在vehicle_scheduling函数中添加容器选择变量
def container_scheduling(tasks):
    model = cp_model.CpModel()
    # 新增容器选择变量
    container_used = [model.NewBoolVar(f'container_{i}') for i in range(len(tasks))]
    
    for i, task in enumerate(tasks):
        # 容器装载时间缩短
        load_time = model.NewIntVar(10, 45, f'load_{i}')
        model.Add(load_time == 10).OnlyEnforceIf(container_used[i])
        model.Add(load_time == 45).OnlyEnforceIf(container_used[i].Not())
        
        # 容量变化
        effective_cap = model.NewIntVar(800, 1000, f'cap_{i}')
        model.Add(effective_cap == 800).OnlyEnforceIf(container_used[i])
        model.Add(effective_cap == 1000).OnlyEnforceIf(container_used[i].Not())
        
        # 更新任务持续时间计算
        model.Add(duration[i] == load_time + task.demand // effective_cap * 60)
    
    # 目标函数增加容器成本
    model.Minimize(sum(container_used) * 500 + ...)  # 假设容器使用成本500/次

3.3 Matlab代码

function schedule = container_scheduling(tasks)
    % 添加容器选择变量
    container = optimvar('container', num_tasks, 'Type','integer','LowerBound',0,'UpperBound',1);
    
    % 装载时间和容量
    load_time = 10*container + 45*(1-container);
    effective_cap = 800*container + 1000*(1-container);
    
    % 更新任务持续时间
    duration = load_time + ceil(tasks.demand ./ effective_cap) * 60;
    
    % 目标函数
    prob.Objective = sum(vehicle) + 500*sum(container);  % 容器成本
end

四、问题四

4.1 问题分析及数学模型

敏感性分析
1、蒙特卡洛模拟：对货量预测值施加±20%随机扰动，统计成本及周转率变化。
2、鲁棒性指标：定义偏差容忍度δ，当δ<10%时调度方案仍有效。

4.2 Python代码

import numpy as np

def monte_carlo_simulation(base_forecast, num_sim=1000):
    cost_results = []
    for _ in range(num_sim):
        # 添加±20%扰动
        perturbed = base_forecast * np.random.uniform(0.8, 1.2, len(base_forecast))
        # 运行调度模型
        cost = run_scheduling(perturbed)  # 假设已实现调度函数
        cost_results.append(cost)
    
    # 分析统计量
    return np.mean(cost_results), np.std(cost_results)

# 示例：基于问题1预测结果
base = pd.read_csv('problem1_forecast.csv')['预测货量']
mean_cost, std_cost = monte_carlo_simulation(base)
print(f"平均成本: {mean_cost}, 标准差: {std_cost}")