用Python和MATLAB搞定数学建模：从报童问题到轧钢浪费，手把手教你搭建概率模型

用Python和MATLAB搞定数学建模：从报童问题到轧钢浪费，手把手教你搭建概率模型

数学建模的魅力在于将现实世界的复杂问题转化为可计算的数学模型。对于初学者来说，最大的挑战往往不是理解概率论的概念，而是不知道如何将这些抽象的理论转化为实际可运行的代码。本文将聚焦两个经典案例——报童问题和轧钢浪费问题，通过Python和MATLAB双语言实现，带你体验从问题分析到代码落地的完整建模流程。

1. 环境准备与工具选择

在开始建模之前，我们需要准备好开发环境。Python和MATLAB各有优势：Python开源免费且生态丰富，MATLAB则在矩阵运算和工程建模方面有独特优势。

1.1 Python环境配置

推荐使用Anaconda管理Python环境，它集成了数据科学常用的库：

conda create -n math_modeling python=3.8 conda activate math_modeling conda install numpy scipy matplotlib pandas jupyter

关键库的作用：

• NumPy：高效数值计算
• SciPy：科学计算工具集
• Matplotlib：数据可视化
• Pandas：数据处理

1.2 MATLAB基础设置

MATLAB开箱即用，但有几个设置能提升效率：

% 在启动脚本(startup.m)中添加 format compact % 紧凑显示输出 set(0,'DefaultFigureWindowStyle','docked') % 停靠图形窗口

提示：对于大型计算任务，MATLAB的并行计算工具箱(Parallel Computing Toolbox)能显著提升性能。

2. 报童问题建模实战

报童问题是一个经典的随机优化问题，核心是确定最优订购量以平衡供不应求和供过于求的风险。

2.1 问题数学化

设：

• 零售价：a
• 购进价：b
• 退回价：c
• 每日需求r的概率密度函数：f(r)

目标函数为期望利润：

$$ G(n) = \sum_{r=0}^n [(a-b)r - (b-c)(n-r)]f(r) + \sum_{r=n+1}^\infty (a-b)nf(r) $$

2.2 Python实现

import numpy as np from scipy.stats import norm import matplotlib.pyplot as plt def newsboy_problem(a, b, c, mu, sigma): """报童问题求解""" critical_ratio = (a - b)/(a - c) # 关键比率 n_opt = norm.ppf(critical_ratio, mu, sigma) # 逆CDF求最优解 return n_opt # 参数设置 a, b, c = 1.0, 0.75, 0.6 # 价格参数 mu, sigma = 500, 50 # 需求分布参数 # 计算最优订购量 optimal_order = newsboy_problem(a, b, c, mu, sigma) print(f"最优订购量：{optimal_order:.0f}份") # 可视化分析 n_values = np.linspace(mu-3*sigma, mu+3*sigma, 100) profits = [] for n in n_values: r = np.linspace(0, mu+3*sigma, 1000) prob = norm.pdf(r, mu, sigma) profit = np.sum(((a-b)*np.minimum(r,n) - (b-c)*np.maximum(n-r,0))*prob) profits.append(profit) plt.plot(n_values, profits) plt.axvline(optimal_order, color='r', linestyle='--') plt.xlabel('订购量') plt.ylabel('期望利润') plt.title('报童问题利润分析') plt.grid(True) plt.show()

2.3 MATLAB实现对比

function optimal_order = newsboy_matlab(a, b, c, mu, sigma) % 报童问题MATLAB实现 critical_ratio = (a - b)/(a - c); optimal_order = norminv(critical_ratio, mu, sigma); % 可视化 n_values = linspace(mu-3*sigma, mu+3*sigma, 100); profits = zeros(size(n_values)); for i = 1:length(n_values) n = n_values(i); r = linspace(0, mu+3*sigma, 1000); prob = normpdf(r, mu, sigma); profit = sum(((a-b)*min(r,n) - (b-c)*max(n-r,0)).*prob); profits(i) = profit; end figure; plot(n_values, profits); hold on; xline(optimal_order, '--r'); xlabel('订购量'); ylabel('期望利润'); title('报童问题利润分析'); grid on; end

注意：实际应用中，需求分布可能不是正态分布。可以通过历史数据拟合更合适的分布，如泊松分布适用于离散需求场景。

3. 轧钢浪费问题建模

轧钢问题需要优化粗轧的均值设定，以最小化因长度不合格造成的材料浪费。

3.1 问题建模

设：

• 成品规定长度：l
• 粗轧长度均值：m（可调）
• 粗轧长度标准差：σ（固定）

目标是最小化每根成品材的平均浪费：

$$ J(m) = \frac{m}{P(m)} \quad \text{其中} \quad P(m) = P(X \geq l) = 1 - \Phi\left(\frac{l-m}{\sigma}\right) $$

3.2 Python实现

from scipy.optimize import minimize_scalar def steel_waste(l, sigma): """轧钢浪费优化""" def objective(m): P = 1 - norm.cdf(l, m, sigma) # 合格概率 return m / P # 初始猜测为l + sigma result = minimize_scalar(objective, bounds=(l, l+3*sigma), method='bounded') return result.x # 参数设置 l = 2.0 # 米 sigma1 = 0.2 # 20厘米 sigma2 = 0.1 # 10厘米 # 优化计算 m_opt1 = steel_waste(l, sigma1) m_opt2 = steel_waste(l, sigma2) print(f"当σ=20cm时，最优均值：{m_opt1:.3f}米") print(f"当σ=10cm时，最优均值：{m_opt2:.3f}米") # 可视化分析 m_values = np.linspace(l, l+0.5, 100) J_values1 = [objective(m) for m in m_values] plt.plot(m_values, J_values1, label='σ=20cm') plt.axvline(m_opt1, color='r', linestyle='--') plt.xlabel('粗轧均值(m)') plt.ylabel('平均浪费J(m)') plt.title('轧钢浪费优化') plt.legend() plt.grid(True) plt.show()

3.3 MATLAB实现

function m_opt = steel_waste_matlab(l, sigma) % 轧钢问题MATLAB实现 objective = @(m) m ./ (1 - normcdf(l, m, sigma)); % 使用fminbnd优化 options = optimset('Display','off'); m_opt = fminbnd(objective, l, l+3*sigma, options); % 可视化 m_values = linspace(l, l+0.5, 100); J_values = arrayfun(objective, m_values); figure; plot(m_values, J_values); hold on; xline(m_opt, '--r'); xlabel('粗轧均值(m)'); ylabel('平均浪费J(m)'); title('轧钢浪费优化'); grid on; end

4. 模型扩展与实战技巧

4.1 蒙特卡洛模拟验证

对于复杂模型，蒙特卡洛模拟是验证解析解的有效方法。以报童问题为例：

def monte_carlo_validation(a, b, c, mu, sigma, n_opt, n_sim=10000): """蒙特卡洛验证""" demands = np.random.normal(mu, sigma, n_sim) profits_opt = a * np.minimum(demands, n_opt) - b * n_opt + c * np.maximum(n_opt - demands, 0) # 测试附近点 n_test = [int(n_opt*0.9), int(n_opt), int(n_opt*1.1)] avg_profits = [] for n in n_test: profits = a * np.minimum(demands, n) - b * n + c * np.maximum(n - demands, 0) avg_profits.append(np.mean(profits)) return n_test, avg_profits n_test, profits = monte_carlo_validation(a, b, c, mu, sigma, optimal_order) plt.bar([str(n) for n in n_test], profits) plt.xlabel('订购量') plt.ylabel('平均利润') plt.title('蒙特卡洛验证') plt.show()

4.2 敏感性分析

考察参数变化对最优解的影响：

# 报童问题价格敏感性分析 b_values = np.linspace(0.7, 0.8, 20) optimal_orders = [newsboy_problem(a, b, c, mu, sigma) for b in b_values] plt.plot(b_values, optimal_orders) plt.xlabel('购进价(b)') plt.ylabel('最优订购量') plt.title('购进价敏感性分析') plt.grid(True)

4.3 实用技巧

1. 分布选择：

   • 对于计数数据（如需求），考虑泊松分布
   • 对于连续正数，考虑伽马分布
   • 对于有界数据，考虑Beta分布

2. 优化加速：

   # 使用Numba加速Python代码 from numba import njit @njit def fast_profit_calc(demands, n, a, b, c): return np.sum(a * np.minimum(demands, n) - b * n + c * np.maximum(n - demands, 0))

3. MATLAB并行计算：

   parfor i = 1:numel(n_values) n = n_values(i); profits(i) = calculate_profit(n); end

5. 常见问题与调试技巧

在实际建模过程中，经常会遇到一些典型问题：

1. 模型不收敛：

   • 检查目标函数是否平滑
   • 尝试不同的初始值
   • 考虑约束条件是否合理

2. 结果不符合直觉：

   • 验证概率分布的选择
   • 检查参数单位是否一致
   • 进行量纲分析

3. 性能瓶颈：

   • 向量化操作替代循环
   • 使用更高效的算法
   • 考虑近似计算

# 诊断示例：检查概率分布拟合 from scipy.stats import probplot # 生成模拟需求数据 demand_samples = np.random.normal(mu, sigma, 1000) probplot(demand_samples, plot=plt) plt.title('正态性检验')

对于MATLAB用户，类似的诊断工具同样可用：

% MATLAB正态性检验 qqplot(demand_samples)

在项目实践中，保持代码的模块化和文档完整非常重要。每个关键函数都应该有清晰的文档字符串：

def critical_ratio(a, b, c): """计算报童问题的关键比率 参数： a: 零售价 b: 购进价 c: 退回价 返回： float: 关键比率，用于确定最优订购量 """ return (a - b) / (a - c)