数学建模的常用知识python

一、规划问题

1、线性规划

1. 定义：线性规划是研究线性约束条件下线性目标函数的极值问题的数学理论和方法。
2. 步骤：

• 根据影响所要达到目的的因素找到决策变量。
• 由决策变量和所在达到目的之间的函数关系确定目标函数。
• 由决策变量所受的限制条件确定决策变量所要满足的约束条件。

1. 特点：目标函数是决策变量的线性函数。根据具体问题可以是最大化或最小化，二者统称为最优化。约束条件也是决策变量的线性函数。

from scipy import optimize
import numpy as np
'''
max z = 2x + 3y - 5z

x + y + z = 7
2x - 5y + z >= 10
x + 3y + z <= 12
x , y , z > 0
'''

c = np.array([2, 3, -5])
A = np.array([[-2, 5, -1], [1, 3, 1]])
B = np.array([-10, 12])
Aeq = np.array([[1, 1, 1]])
beq = np.array([7])

# 求最小要在c前面加一个负号
res = optimize.linprog(-c, A, B, Aeq, beq)
print(res)
'''
con: array([1.807134e-09])
fun: -14.571428565645057
message: 'Optimization terminated successfully.'
nit: 5
slack: array([-2.24586572e-10,  3.85714286e+00])
status: 0
return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
success: True
x: array([6.42857143e+00, 5.71428571e-01, 2.35900788e-10])
'''

2、整数规划

1. 定义：整数规划是指规划中的变量（全部或部分）限制为整数，若在线性模型中，变量限制为整数，则称为整数线性规划。目前所流行的求解整数规划的方法往往只适用于整数线性规划。
2. 分类：在整数规划中，如果所有变量都限制为整数，则称为纯整数规划；如果仅一部分限制为整数，则称为混合整数规划。整数规划的一种特殊情形是0-1规划，它的变数仅限于0或1。
3. 0-1规z划：问题中许多量具有不可分割的性质（最优调度的车辆数、设置的销售网点数…），或者问题的解必须满足一些特殊的约束条件（满足逻辑条件、顺序…），需引入逻辑变量（0-1变量）以表示“是”与“非”。这类问题的模型均为整数规划。

#导入numpy
import numpy as np
#导入numpy
import cvxpy as cp

#设置目标函数中变量个数
n=3

#输入目标函数的系数
c=np.array([3,1,3])

#输入约束条件的系数矩阵（3×3）
a=np.array([[-1,2,1],[0,4,-3],[1,-3,2]])

#输入b值（3×1）
b=np.array([4,2,3])

#创建x，个数是3
x=cp.Variable(n,integer=True)

#明确目标函数（此时c是3×1，x是3×1,但python里面可以相乘）
objective=cp.Maximize(cp.sum(c*x))

#明确约束条件，其中a是3×3，x是3×1,a*x=b(b为3×1的矩阵)
constriants=[0<=x,a*x<=b]
#求解问题
prob=cp.Problem(objective,constriants)
#这里solver必须使用cp.CPLEX,否则计算不出来，而CPLEX需要pip intall CPLEX(建议使用清华镜像)
resluts=prob.solve(solver=cp.CPLEX)

#输入结果
print(prob.value)#目标函数的值
print(x.value)#各x的值

3、非线性规划

• 非线性规划可以简单分两种，目标函数为凸函数or非凸函数凸函数的非线性规划，比如fun = x^2 +y^2 +xy,有很多常用库完成，比如cvxpy
• 非凸函数的非线性规划（求极值），可以尝试以下方法:
  纯数学方法，求导求极值
• 神经网络、深度学习(反向传播算法中链式求导过程)scipy. optimize.minimize

from scipy.optimize import minimize
import numpy as np


# 计算(2+x1)/ (1+x2)- 3*x1+4*x3的最小值 x1,x2,x3的范围都在0.1到0.9 之间
def fun(args):
    a, b, c, d = args
    v = lambda x: (a + x[0]) / (b + x[1]) - c * x[0] + d * x[2]
    return v


def con(args):
    # 约束条件 分为eq 和ineq
    # eq表示 函数结果等于0 ； ineq 表示 表达式大于等于0
    x1min, x1max, x2min, x2max, x3min, x3max = args
    cons = ({'type': 'ineq', 'fun': lambda x: x[0] - x1min},
            {'type': 'ineq', 'fun': lambda x: -x[0] + x1max},
            {'type': 'ineq', 'fun': lambda x: x[1] - x2min},
            {'type': 'ineq', 'fun': lambda x: -x[1] + x2max},
            {'type': 'ineq', 'fun': lambda x: x[2] - x3min},
            {'type': 'ineq', 'fun': lambda x: -x[2] + x3max})
    return cons


if __name__ == "__main__":
    # 定义常量值
    args = (2, 1, 3, 4)  # a,b,c,d
    # 设置参数范围/约束条件
    args1 = (0.1, 0.9, 0.1, 0.9, 0.1, 0.9)  # x1min, x1max, x2min,x2max
    cons = con(args1)
    # 设置初始猜测值
    x0 = np.asarray((0.5, 0.5, 0.5))
    res = minimize(fun(args), x0, method='SLSQP', constraints=cons)
    print(res.fun)
    print(res.success)
    print(res.x)

'''
-0.773684210526435
True
[0.9 0.9 0.1]
'''

二、数值逼近问题

1、一维插值

1. 定义：在离散数据的基础上补插连续函数，使得这条连续曲线通过全部给定的离散数据点。
2. 作用：插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值情况，估算出函数在其他点处的近似值。
3. 区别：从定义上看，插值和拟合有一定的相似度，但插值要求近似函数通过给定的所有离散数据，而拟合并不要求这样，只要近似函数能较好的反映数据变化的趋势即可（近似含义不同），当测量值是准确的，没有误差时，一般用插值；当测量值与真实值有误差时，一般用数据拟合。

2、二维插值

3、最小二乘拟合

拟合指的是已知某函数的若干离散函数值{f1,f2,…,fn}，通过调整该函 数中若干待定系数f(λ1, λ2,…,λn)，使得该函数与已知点集的差别(最小 二乘意义)最小。

三、微分方程问题

0、分类

• 常微分方程
• 偏微分方程

1、解析解

利用SymPy库进行求解

2、数值解

用odeint求其数值解

3、传染病

传染病模型主要包括：SI、SIS、SIR、SIRS、SEIR、 SEIRS共六个模型。

四、建模图论问题

1、Dijkstra

2、Floyd

五、差分方程问题

1、递推关系

差分方程建模的关键在于如何得到第n组数据与第n+1组数据之间的关系。

2、显式差分

3、马尔科夫链

六、蒙特卡洛算法

​ 当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种 “实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的 某些数字特征，并将其作为问题的解。

• 样例：经典的蒙特卡洛方法求圆周率
• 基本思想：在图中区域产生足够多的随机数点，然后计算落在圆内 的点的个数与总个数的比值再乘以4，就是圆周率。

M = input(‘请输入一个较大的整数’)
N = 0
import math
import random
for i in range(int(M)):
    x = random.random()
    y = random.random()
if math.sqrt(x**2+y**2)<1:
N = N+1
pi = 4*N/int(M)
print(pi)

七、时间序列

时间序列预测在日常分析中常会用到，是重要的时序数据处理方法。

0、误差计算

from sklearn.metrics import mean_squared_error
from math import sqrt
rms = sqrt(mean_squared_error(test['Count'], y_hat['naive']))
print(rms)

1、朴素法

如果数据集在一段时间内都很稳定，想要预测第二天的价格，可以取前面一天的价格，预测第二天的值。这种假设第一个预测点和上一个观测点相等的预测方法就叫朴素法。

• Hence y(t+1) = y(t)

y_hat[‘naive’]=dd[len(dd)-1]

朴素法并不适合变化很大的数据集，最适合稳定性很高的数据集。

2、简单平均法

物品价格会随机上涨和下跌，平均价格会保持一致。我们经常会遇到一些数据集，虽然在一定时期内出现小幅变动，但每个时间段的平均值确实保持不变。这种情况下，我们可以预测出第二天的价格大致和过去天数的价格平均值一致。这种将预期值等同于之前所有观测点的平均值的预测方法就叫简单平均法。

y_hat_avg = test.copy()
y_hat_avg['avg_forecast'] = train['Count'].mean()

3、 移动平均法

明显这里的逻辑是只有最近的值最要紧，用某些窗口期计算平均值的预测方法就叫移动平均法。
针算移动平均值涉及到一个有时被称为“滑动窗口”的大小值p。使用简单的移动平均模型，我们可以根据之前数值的固定有限数p的平均值预测某个时序中的下一个值。

y_hat_avg = test.copy()
y_hat_avg['moving_avg_forecast'] = train['Count'].rolling(60).mean().iloc[-1]

4、简单指数平滑法

相比更早时期内的观测值，给近期的观测值赋予更大的权重。按照这种原则工作的方法叫做简单指数平滑法。通过加权平均值计算出预测值，其中权重随着观测值从早期到晚期的变化呈指数级下降，最小的权重和最早的观测值相关:

• Hence y(t+1) = αy(t) + α(1-α)y(t-1) + α(1-α)^2 y(t-2) + ······ (0<=α<=1 平滑参数)

rom statsmodels.tsa.api import SimpleExpSmoothing
y_hat_avg = test.copy()
fit=SimpleExpSmoothing(np.asarray(train['Count'])).fit(smoothinglevel=0.6,optimized=False)
y hat_avg ['SES']= fit.forecast(len(test))

5、霍尔特线性趋势法

如果物品的价格是不断上涨的，上面的方法并没有考虑这种趋势，即我们在一段时间内观察到的价格的总体模式。在上图例子中，我们可以看到物品的价格呈上涨趋势。虽然上面这些方法都可以应用于这种趋势，但我们仍需要一种方法可以在无需假设的情况下，准确预测出价格趋势。这种考虑到数据集变化趋势的方法就叫做霍尔特线性趋势法。

每个时序数据集可以分解为相应的几个部分：趋势（Trend），季节性(Seasonal)和残(Residual)。 任何呈现某种趋势的数据集都可以用霍尔特线性趋势法用于预测。

from statsmodels.tsa.api import Holt
y_hat_avg = test.copy()
fit=Holt(np.asarray(train['Count'])).fit(smoothing_level=0.3,smoothing_slope=0.1)
y_hat_avg['Holt_linear'] = fit.forecast(len(test))

6、 Holt-Winters季节性预测模型

在应用这种算法前，我们先介绍一个新术语。假如有家酒店坐落在半山腰上，夏季的时候生意很好，顾客很多，但每年其余时间顾客很少。因此，每年夏季的收入会远高于其它季节，而且每年都是这样，那么这种重复现象叫做“季节性”(Seasonality)。如果数据集在一定时间段内的固定区间内呈现相似的模式，那么该数据集就具有季节性。

from statsmodels.tsa.api import ExponentialSmoothing
y_hat_avg = test.copy()
fit1 = ExponentialSmoothing(np.asarray(train['Count']),seasonal_periods=7,trend='add',seasonal='add',).fit()
y_hat_avg['Holt_Winter'] = fit1.forecast(len(test))

7、自回归移动平均模型（ARIMA）

另一个场景的时序模型是自回归移动平均模型（ARIMA）。指数平滑模型都是基于数据中的趋势和季节 性的描述，而自回归移动平均模型的目标是描述数据中彼此之间的关系。ARIMA的一个优化版就是季节 性ARIMA。它像Holt-Winters季节性预测模型一样，也把数据集的季节性考虑在内

import statsmodels.api as sm
y_hat_avg = test.copy()
fit1 = sm.tsa.statespace.SARIMAX(train.Count, order=(2, 1, 4),seasonal_order=(0, 1, 1, 7)).fit()
y_hat_avg['SARIMA']=fit1.predict(start="2013-11-1",end="2013-12-31",dynamic=True)

八、支持向量机SVM

在进行线性分类时,将分类面取在离两类样本距离较大的地方;进行非线性分类时通过高维空间变换,将非线性分类变成高维空间的线性分类问题。一般支持向量机分类采用的是线性可分 SVM。

个人觉得搭配pca降维简单又快捷。

九、主成分分析

降维，找到数据中的主成分，并利用这些主成分表征原始数据，从而达到降维的目的。

1. 对样本数据进行中心化处理；
2. 求样本协方差矩阵；
3. 对协方差矩阵进行特征值分解，将特征值从大到小排列；
4. 取特征值前 n 个最大的对应的特征向量 W1, W2, …, Wn ，这样将原来 m 维的样本降低到 n 维。

通过 PCA ，就可以将方差较小的特征给抛弃，这里，特征向量可以理解为坐标转换中新坐标轴的方向，特征值表示在对应特征向量上的方差，特征值越大，方差越大，信息量也就越大。这也是为什么选择前 n 个最大的特征值对应的特征向量，因为这些特征包含更多重要的信息。

PCA 是一种线性降维方法，这也是它的一个局限性。不过也有很多解决方法，比如采用核映射对 PCA 进行拓展得到核主成分分析(KPCA)，或者是采用流形映射的降维方法，比如等距映射、局部线性嵌入、拉普拉斯特征映射等，对一些 PCA 效果不好的复杂数据集进行非线性降维操作。

十、马尔可夫预测

对事件的全面预测，不仅要能够指出事件发生的各种可能结果，而且还必须给出每一种结果出现的概率，说明被预测的事件在预测期内出现每一种结果的可能性程度。这就是关于事件发生的概率预测。

马尔可夫（Markov)预测法，就是一种关于事件发生的概率预测方法。它是根据事件的目前状况来预测其将来各个时刻（或时期)变动状况的一种预测方法。马尔可夫预测法是地理预测研究中重要的预测方法之一。

十一、随机森林：

集成思想，涉及到决策树和集成学习，将若干个弱分类器的分类结果进行投票选择，从而组成一个强分类器。

随机森林的既可以用于回归也可以用于分类任务，并且很容易查看模型的输入特征的相对重要性。随机森林算法被认为是一种非常方便且易于使用的算法，因为它是默认的超参数通常会产生一个很好的预测结果。超参数的数量也不是那么多，而且它们所代表的含义直观易懂。随机森林有足够多的树，分类器就不会产生过度拟合模型。由于使用大量的树会使算法变得很慢，并且无法做到实时预测。一般而言，这些算法训练速度很快，预测十分缓慢。越准确的预测需要越多的树，这将导致模型越慢。在大多数现实世界的应用中，随机森林算法已经足够快，但肯定会遇到实时性要求很高的情况，那就只能首选其他方法。当然，随机森林是一种预测性建模工具，而不是一种描述性工具。也就是说，如果您正在寻找关于数据中关系的描述，那建议首选其他方法。

十二、模拟退火算法

在解空间随机寻找目标函数的全局最优解

十三、遗传算法

最优解，将方程求解问题转化为生存问题。

十四、熵值法

熵值法是一种客观赋权法,其根据各项指标观测值所提供的信息的大小来确定指标权重。设有 m个待评方案，n 项评价指标，形成原始指标数据矩阵 ，对于某项指标 ，指标值 的差距越大，则该指标在综合评价中所起的作用越大；如果某项指标的指标值全部相等，则该指标在综合评价中不起作用。

在信息论中，熵是对不确定性的一种度量。信息量越大，不确定性就越小，熵也就越小；信息量越小，不确定性就越大，熵也越大.根据熵的特性，我们可以通过计算熵值来判断一个方案的随机性及无序程度，也可以用熵值来判断某个指标的离散程度，指标的离散程度越大，该指标对综合评价的影响越大。因此，可根据各项指标的变异程度，利用信息熵这个工具，计算出各个指标的权重，为多指标综合评价提供依据。

十五、参考文章

1. 数学建模重要算法简介及算法实现
2. 数学建模常用知识与方法简介
3. 使用Python求解整数规划问题
4. Python 在数学建模中的应用