用 Python 做数学建模

2023-07-20 13:35| 来源: 网络整理| 查看: 265

本文由 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 = 的项对应的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 x_1=6.42857143, x_2=0.57142857, x_3=0 x1=6.42857143,x2=0.57142857,x3=0, 对应的最优值 z = 14.571428571428571 z=14.571428571428571 z=14.571428571428571.

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

整数规划

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

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

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

e.g. y = x 2 y=x^2 y=x2、 y = 12 − x y=12-x y=12−x 与 x 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] list: res = [] res.append(sum(x) - 400) res.append(x[0] + 2 * x[1] + 2 * x[2] + x[3] + 6 * x[4] - 800) res.append(2 * x[0] + x[1] + 6 * x[2] - 200) res.append(x[2] + x[3] + 5 * x[4] - 200) return res random.seed(time.time) pb = 0 xb = [] for i in range(10 ** 6): x = [random.randint(0, 99) for i in range(5)] # 产生一行五列的区间[0, 99] 上的随机整数 rf = f(x) rg = g(x) if all((a > import numpy as np >>> from scipy import optimize >>> c = [[3,8,2,10,3], [8,7,2,9,7], [6,4,2,7,5], [8,4,2,3,5], [9,10,6,9,10]] >>> c = np.array(c) >>> optimize.linear_sum_assignment(c) (array([0, 1, 2, 3, 4]), array([4, 2, 1, 3, 0])) # 对应 x15、x23、x32、x44、x51 >>> c[[0, 1, 2, 3, 4], [4, 2, 1, 3, 0]].sum() # 结果代入 cost_matrix，求得最优值 21 非线性规划非线性规划模型

第三方依赖库：numpy , scipy。

我们之前多次使用的 scipy.optimize 中集成了一系列用来求规划的函数，其中当然不乏解决非线性规划的方法。例如用其中的 minimize 函数就可以解决很多在 Matlab 中用 fmincon 解的非线性规划问题。

minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

常用的参数：

fun:：待求最小值的目标函数，fun(x, *args) -> float， x 是 shape (n,)的 1-D arrayx0：初始猜测值， shape (n,)的 1-D arraybounds：设置参数范围/约束条件，tuple，形如 ((0, None), (0, None))constraints：约束条件，放一系列 dict 的 tuple，({'type': TYPE, 'fun': FUN}, ...) TYPE：'eq'表示函数结果等于0 ； 'ineq' 表示表达式大于等于0FUN：约束函数

e.g. 求下列非线性规划： min f ( x ) = x 1 2 + x 2 2 + x 3 2 + 8 , s.t. { x 1 2 + x 2 2 + x 3 2 ≥ 0 , x 1 + x 2 2 + x 3 2 ≤ 20 , − x 1 − x 2 2 + 2 = 0 , x 2 + 2 x 3 2 = 3 , x 1 , x 2 , x 3 ≥ 0. \textrm{min } f(x)=x_1^2+x_2^2+x_3^2+8,\\ \textrm{s.t. } \left \{ \begin{array}{ll} x_1^2+x_2^2+x_3^2 \ge 0,\\ x_1+x_2^2+x_3^2 \le 20,\\ -x_1-x_2^2+2=0,\\ x_2+2x_3^2=3,\\ x_1,x_2,x_3 \ge 0. \end{array}\right. min f(x)=x12+x22+x32+8,s.t. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x12+x22+x32≥0,x1+x22+x32≤20,−x1−x22+2=0,x2+2x32=3,x1,x2,x3≥0.

import numpy as np from scipy import optimize f = lambda x: x[0] ** 2 + x[1] **2 + x[2] ** 2 + 8 # Notice：eq ==; ineq >= cons = ({'type': 'ineq', 'fun': lambda x: x[0]**2 - x[1] + x[2]**2}, {'type': 'ineq', 'fun': lambda x: -x[0] - x[1] - x[2]**3 + 20}, {'type': 'eq', 'fun': lambda x: -x[0] - x[1]**2 + 2}, {'type': 'eq', 'fun': lambda x: x[1] + 2 * x[2]**2 - 3}) res = optimize.minimize(f, (0, 0, 0), constraints=cons) print(res)

输出：

fun: 10.651091840572583 jac: array([1.10433471, 2.40651834, 1.89564812]) message: 'Optimization terminated successfully.' nfev: 86 nit: 15 njev: 15 status: 0 success: True x: array([0.55216734, 1.20325918, 0.94782404])

即，求得当 ( x 1 , x 2 , x 3 ) = ( 0.55216734 , 1.20325918 , 0.94782404 ) (x_1,x_2,x_3)=(0.55216734, 1.20325918, 0.94782404) (x1,x2,x3)=(0.55216734,1.20325918,0.94782404) 时，最小值 y = 10.651091840572583 y=10.651091840572583 y=10.651091840572583.

无约束问题的 Python 解法符号解

第三方依赖库：sympy。

e.g. 求多元函数 f ( x ， y ) = x 3 − y 3 + 3 x 2 + 3 y 2 − 9 x f(x，y) = x^3 - y^3 + 3x^2 + 3y^2 - 9x f(x，y)=x3−y3+3x2+3y2−9x 的极值

思路：求驻点，代入Hessian矩阵，正定则为极小值，负定极大，不定不是极值点。

解：

import sympy f, x, y = sympy.symbols("f x y") f = x**3 - y**3 + 3 * x**2 + 3 * y**2 - 9 * x funs = sympy.Matrix([f]) args = sympy.Matrix([x, y]) df = funs.jacobian(args) # 一阶偏导 d2f = df.jacobian(args) # Hessian 矩阵 stationaryPoints = sympy.solve(df) # 驻点 for i in stationaryPoints: a = d2f.subs(x, i[x]).subs(y, i[y]) # 驻点处的Hessian b = a.eigenvals(multiple=True) # 求Hessian矩阵的特征值 fv = f.subs(x, i[x]).subs(y, i[y]) # 驻点处的函数值 if all((j > 0 for j in b)): print('点({x}, {y})是极小值点，对应的极小值为: f({x}, {y}) = {f}'.format(x=i[x], y=i[y], f=fv)) elif all((j >> import sympy >>> x = sympy.Symbol('x') >>> s = sympy.solve(x**3 - x**2 + 2 * x -3) # 求符号解 >>> s [1/3 + (-1/2 - sqrt(3)*I/2)*(65/54 + 5*sqrt(21)/18)**(1/3) - 5/(9*(-1/2 - sqrt(3)*I/2)*(65/54 + 5*sqrt(21)/18)**(1/3)), 1/3 - 5/(9*(-1/2 + sqrt(3)*I/2)*(65/54 + 5*sqrt(21)/18)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(65/54 + 5*sqrt(21)/18)**(1/3), -5/(9*(65/54 + 5*sqrt(21)/18)**(1/3)) + 1/3 + (65/54 + 5*sqrt(21)/18)**(1/3)] >>> for i in s: ... print(i.evalf()) # 近似数值 ... -0.137841101825493 - 1.52731225088663*I -0.137841101825493 + 1.52731225088663*I 1.27568220365098 求如下方程组的解

{ x 2 + y − 6 = 0 y 2 + x − 6 = 0 \left\{\begin{array}{l} x^2+y-6=0\\ y^2+x-6=0\\ \end{array}\right. {x2+y−6=0y2+x−6=0

>>> import sympy >>> x, y = sympy.symbols('x y') >>> s = sympy.solve((x**2+y-6, y**2+x-6)) >>> s [{x: -3, y: -3}, {x: 2, y: 2}, {x: 6 - (1/2 - sqrt(21)/2)**2, y: 1/2 - sqrt(21)/2}, {x: 6 - (1/2 + sqrt(21)/2)**2, y: 1/2 + sqrt(21)/2}] >>> for i in s: ... print('({x}, {y})'.format(x=i[x].evalf(), y=i[y].evalf())) ... (-3.00000000000000, -3.00000000000000) (2.00000000000000, 2.00000000000000) (2.79128784747792, -1.79128784747792) (-1.79128784747792, 2.79128784747792) 约束极值问题

参考上文 “无约束问题的 Python 解法”。

附：numpy

推荐一份质量比较高的 numpy 中文文档：https://www.numpy.org.cn/index.html.

创建向量/矩阵

用 np.array()

>>> import numpy as np >>> np.array([1, 2, 3]) array([1, 2, 3]) >>> np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) 矩阵的拼接列合并/扩展：np.column_stack()行合并/扩展：np.row_stack() >>> import numpy as np >>> A = np.array([[1,1], [2,2]]) >>> B = np.array([[3, 3], [4, 4]]) >>> A array([[1, 1], [2, 2]]) >>> B array([[3, 3], [4, 4]]) >>> np.column_stack((A, B)) # 行 array([[1, 1, 3, 3], [2, 2, 4, 4]]) >>> np.row_stack((A, B)) # 列 array([[1, 1], [2, 2], [3, 3], [4, 4]]) 对角线元素赋值

比如说我们有一个 np.array X，我想对角线的所有值设置为0，可以用 np.fill_diagonal(X, [0, 0, 0, ...])。

>>> D = np.zeros([3, 3]) >>> np.fill_diagonal(D, [1, 2, 3]) >>> D array([[1., 0., 0.], [0., 2., 0.], [0., 0., 3.]])

…

未完但不用待续了，应该没有后续了，我对用 python 做数学建模以及数学建模都失去兴趣了。

【本文地址】

公司简介

联系我们