用 Python 做数学建模

---

本文由 CDFMLR 原创，收录于个人主页 https://clownote.github.io，并同时发布到 CSDN。本人不保证 CSDN 排版正确，敬请访问 clownote 以获得良好的阅读体验。

---

用 Python 做数学建模

⚠️【注意】这是一篇没有完成也不会被完成的文章，这里我只写了几个简单问题的 Python 解法，我发布它只是希望说明 Python 做数学建模是有可行性的。

线性规划

第三方依赖库：scipy。

用 scipy.optimize.linprog 可以解线性规划问题：

linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='simplex', callback=None, options=None)

其规定的问题标准型式为：

        Minimize:     c^T * x
    
        Subject to:   A_ub * x <= b_ub
                      A_eq * x == b_eq

e.g. 求接下列线性规划问题：
$$\begin{gathered} \text{max }z=2x_1+3x_2-5x_3, \\ \text{s.t. }\{\begin{array}{l}{x}_1+x_2+x_3=7\\ 2x_1-5x_2+x_3\ge 10\\ x_1+3x_2+x_3\le 12\\ x_1,x_2,x_3\ge 0\end{array} \end{gathered}$$
解：

#! /usr/bin/python3

import numpy as np
from scipy import optimize


c = [-2, -3, 5]
a = [[-2, 5, -1], [1, 3, 1]]
b = [-10, 12]
aeq = [[1, 1, 1]]
beq = [7]
bounds = [[0, None], [0, None], [0, None]]	# (0, None) means non-negative, this is a default value
result = optimize.linprog(c, a, b, aeq, beq, bounds)
print(result)

⚠️【注意】这里题目是求max，标准型是min，所以在写c矩阵的时候把值都写成了题目中的负；类似地， >= 的项对应的a、b中值要取之负。然后最终结果也要取fun的负。

输出：

     fun: -14.571428571428571
 message: 'Optimization terminated successfully.'
     nit: 2
   slack: array([3.85714286, 0.        ])
  status: 0
 success: True
       x: array([6.42857143, 0.57142857, 0.        ])

最优解为$x_1=6.42857143,x_2=0.57142857,x_3=0$, 对应的最优值$z=14.571428571428571$.

若要取 fun、x的值，可以直接用 result.fun 和 result.x。

整数规划

线性整数规划问题可以转换为线性规划问题求解。

对于非线性的整数规划，在穷举不可为时，在一定计算量下可以考虑用 蒙特卡洛法 得到一个满意解。

蒙特卡洛法（随机取样法）

e.g. $y=x^2$、$y=12-x$ 与 $x$ 轴 在第一象限围成一个曲边三角形。设计一个随机试验，求该图像面积的近似值。

解： 设计的随机试验思想如下：在矩形区域 [0, 12] * [0, 9] 上产生服从均匀分布的 10^7 个随机点，统计随机点落在曲边三角形的频数，则曲边三角形的面积近似为上述矩形面积乘于频率。

import random


x = [random.random() * 12 for i in range(0, 10000000)]
y = [random.random() * 9 for i in range(0, 10000000)]

p = 0
for i in range(0, 10000000):
    if x[i] <= 3 and y[i] < x[i] ** 2:
            p += 1
    elif x[i] > 3 and y[i] < 12 - x[i]:
            p += 1
 
area_appr = 12 * 9 * p / 10 ** 7
print(area_appr)

结果在 49.5 附近。

e.g. 已知非线性整数规划为：
max  z = x 1 2 + x 2 2 + 3 x 3 2 + 4 x 4 2 + 2 x 5 2 − 8 x 1 − 2 x 2 − 3 x 3 − x 4 − 2 x 5 s.t.  { 0 ≤ x