Supercised Learning & Reinforesment Learning

决策树分类

前提

假设我们有一个包含 $N$ 个样本的训练集，每个样本包含 $M$ 个特征，用 $X_{i,j}$ 表示第 $i$ 个样本的第 $j$ 个特征的取值，用 $y_i$ 表示第 $i$ 个样本的目标值（或标签）。决策树模型的目标是学习从特征到目标值的映射关系，并用学习到的模型对新的样本进行预测。

流程

计算特征阈值和判断最优特征

通过基尼指数判断最优特征

基尼指数越小表示数据集越纯。

对于一个特征 $A$ 和其取值 $a_i$，假设 $D$ 表示当前节点的数据集，$D_i$ 表示特征 $A$ 取值为 $a_i$ 时的子数据集，则基尼指数计算公式为：

$$Gini(D,A)=\sum_{i=1}^{k}\frac{|D_i|}{|D|}Gini(D_i)$$

其中，$k$ 表示特征 $A$ 的取值数目，$|D_i|$ 表示子数据集 $D_i$ 中样本的数目，$|D|$ 表示当前节点的样本数目，$Gini(D_i)$ 表示子数据集 $D_i$ 的基尼指数，计算公式为：

$$Gini(D_i)=1-\sum_{j=1}^{c}p_{ij}^2$$

其中，$c$ 表示类别数目，$p_{ij}$ 表示子数据集 $D_i$ 中属于第 $j$ 类的样本占比。

基尼指数越小表示使用特征 $A$ 进行划分后，得到的子数据集的纯度越高。因此，决策树算法选择基尼指数最小的特征作为最优特征进行划分。

计算特征阈值(用信息增益表示特征阈值)

信息增益是评价特征重要性的一种指标，它表示特征 $x_j$ 对分类的贡献程度。信息增益越大，说明特征的重要性越高。可以通过以下公式计算信息增益：

$$ \begin{aligned} Gain(D,x_j) &= H(D) - \sum_{i=1}^m \frac{|D_i|}{|D|} H(D_i) \\ H &= -\sum p(x) \times log_2[p(x)] \end{aligned} $$

其中 $H(D)$ 表示数据集 $D$ 的熵，$m$ 表示特征 $x_j$的可能取值个数，$D_i$ 表示数据集 $D$ 在特征 $x_j$ 取值为 $i$ 时对应的子集，$|D_i|$ 表示子集 $D_i$ 的大小，$|D|$ 表示数据集 $D$ 的大小。

对特征 $x_j$ 的取值进行排序，得到一个有序的列表 $[x_1,x_2,\cdots,x_m]$
对于相邻的取值 $x_i$ 和 $x_{i+1}$，将它们作为特征阈值 $s$，将数据集划分为两个子集 $D_1={x\in D|x_j\leq s}$ 和 $D_2={x\in D|x_j>s}$
计算使用特征阈值 $s$ 进行划分的信息增益 $Gain(D,x_j,s)$

根据 $x_j \leq s$ 和 $x_j > s$ 将数据集 $D$ 划分为两个子集 $D_1$ 和 $D_2$ ，然后计算信息增益

$$Gain(D,x_j,s) = H(D) - \frac{|D_1|}{|D|} H(D_1) - \frac{|D_2|}{|D|} H(D_2)$$

重复步骤2和3，直到计算完所有相邻取值的信息增益为止。
选择信息增益最大的特征阈值作为最终的特征阈值。

筛选最优特征

计算训练集的经验熵（entropy）$H(D)$，表示训练集中样本的不确定性

对于每个特征$A$，计算其对训练集的条件熵（conditional entropy）$H(D|A)$，表示在特征$A$已知的条件下，训练集中样本的不确定性，公式为： $$H(D|A)=\sum_{i=1}^{|A|} \frac{|D_i|}{|D|} H(D_i)$$ 其中，$|A|$表示特征$A$可能的取值个数，$|D_i|$表示在特征$A$取值为第$i$个值时，训练集中样本的个数，$H(D_i)$表示在特征$A$取值为第$i$个值时，训练集中样本的经验熵。

计算信息增益（information gain）$g(D,A)$，表示特征$A$对训练集划分的贡献程度，公式为： $$g(D,A)=H(D)-H(D|A)$$

选择信息增益最大的特征作为当前节点的划分特征，将训练集划分成更加纯净的子集，并递归地重复上述步骤，直到所有子集都变得足够纯或达到预定的停止条件为止，生成一棵完整的决策树。

需要注意的是，信息增益存在一个缺陷，即它倾向于选择可能取值较多的特征，因为这些特征往往可以划分出更多的子集。为了克服这个缺陷，可以使用其他的划分指标，例如基尼不纯度（Gini impurity）或者增益率（gain ratio），来选择最优的特征。

信息增益比（Information gain ratio）

信息增益比是信息增益的一种改进，它考虑了特征取值数目对信息增益的影响。可以通过以下公式计算信息增益比：

$$Gain \ ratio \ (D,x_j) = \frac{Gain(D,x_j)}{IV(x_j)}$$

其中 $IV(x_j)$ 表示特征 $x_j$ 的固有值，表示特征 $x_j$ 取值的多样性，可以通过以下公式计算：

$$IV(x_j) = - \sum_{i=1}^m \frac{|D_i|}{|D|} \log_2 \frac{|D_i|}{|D|}$$

在特征 $x_j$ 的阈值为 $s$ 时，可以根据 $x_j \leq s$ 和 $x_j > s$ 将数据集 $D$ 划分为两个子集 $D_1$ 和 $D_2$，然后计算信息增益比。

总结

以上三种方法都可以用于计算特征阈值，其中基尼指数适用于连续和离散特征，而信息增益和信息增益比适用于离散特征。在实际应用中，我们可以根据具体问题选择不同的评价指标来计算特征阈值，并选择最优的阈值来进行样本划分。

结点判定规则

假设第 $t$ 个节点的判定规则为 $f_t(x)$，其中 $x$ 表示一个样本的特征向量，$f_t(x)$ 返回一个取值为 $0$ 或 $1$ 的标记，表示该样本应该进入该节点的左子树（$f_t(x)=0$）还是右子树（$f_t(x)=1$）。可以将判定规则表示为：

$$f_t(x) = \begin{cases} 0, & \text{if}\ x \in R_{t,\text{left}} \ 1, & \text{otherwise} \end{cases}$$

其中 $R_{t,\text{left}}$ 是第 $t$ 个节点左子树的样本所在区域，可以用一个特征阈值组合来表示：

$$R_{t,\text{left}} = {x | x_j \leq s_t}$$

其中 $x_j$ 表示样本 $x$ 的第 $j$ 个特征的取值，$s_t$ 是第 $t$ 个节点的特征阈值。右子树的样本所在区域可以表示为 $R_{t,\text{right}} = { x | x_j > s_t }$。

树的结构

决策树模型的树的结构可以用一个树形图或者一组规则来表示。假设决策树模型的根节点为 $t=1$，可以用以下的数学公式来表示该模型：

$$\hat{y}i = \sum{t=1}^T c_t \cdot I(x_i \in R_{t})$$

其中 $\hat{y}_i$ 表示用决策树模型预测样本 $i$ 的目标值，$T$ 表示树的深度或节点个数，$c_t$ 表示第 $t$ 个节点的预测值，可以是该节点所属的样本的平均值或者分类问题中该节点所属的样本中出现最多的类别值， $I(x_i \in R_t)$ 是一个指示函数，如果样本 $i$ 属于第 $t$ 个节点的子树，取值为 $1$，否则取值为 $0$

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19


import pandas as pd
from sklearn.tree import DecisionTreeClassifier, export_graphviz
from sklearn.preprocessing import OneHotEncoder
import graphviz

df = pd.read_excel('~/数模/原表格.xlsx',header=None)
X = df.iloc[1:98, [2,3,4,8,24]]
y = df.iloc[1:98, 26]
y = y.astype('str')

model = DecisionTreeClassifier(criterion='gini', splitter='best',max_depth=3, max_features=5, min_samples_split=3)
model.fit(X, y)

# 将树可视化为图形
dot_data = export_graphviz(model, out_file=None, feature_names=['年龄','术前CEA','距肛距离','放疗与手术时间差(W)','阳性淋巴结初始数量估计'], class_names= ['LRG0','LRG1','LRG2','LRGX'], filled=True, rounded=True, special_characters=True)

graph = graphviz.Source(dot_data)
graph.render("决策树") 
graph.view()

随机森林分类

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28


import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# 读取数据集
df = pd.read_excel('', sheet_name='')
X = df.iloc[:, 1:6]
y = df.iloc[:, 0]

# 如果标签是字符
y = y.astype('str')

# 划分数据集为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 创建随机森林模型
model = RandomForestClassifier(n_estimators=10, criterion='gini', max_depth=3, max_features=3, min_samples_split=10)

# 在训练集上训练模型
model.fit(X_train, y_train)

# 在测试集上进行预测
y_pred = model.predict(X_test)

# 计算模型在测试集上的准确性
accuracy = accuracy_score(y_test, y_pred)
print("模型在测试集上的准确性：", accuracy)

XGBoost(eXtreme Gradient Boosting) + 网格搜索 + 交叉验证

原理

网格搜索

穷举所有给出的参数取值的组合
根据交叉验证计算模型质量
给出质量最高的模型

交叉验证

将数据集分成k个相等大小的子集（折）。
对于每个折，将其作为验证集，剩余的k-1个折作为训练集。
使用训练集进行模型的训练，并在验证集上进行性能评估。
计算模型在当前折上的评估指标（如准确率、F1 分数、均方误差等）。
重复步骤2-4，直到每个折都作为验证集并得到了一个评估指标。
统计k次评估指标的平均值作为模型的最终评估结果。

XGBoost

梯度提升树原理：XGBoost是基于梯度提升树（Gradient Boosting Tree）的算法。梯度提升树是一种迭代的集成学习算法，通过迭代地训练一系列决策树来逐步改进模型的预测能力。每个新的决策树都会尝试纠正前一棵树的预测误差，从而逐步减小模型的残差。
XGBoost的目标函数：XGBoost的目标函数由两部分组成：损失函数（Loss Function）和正则化项（Regularization Term）。损失函数衡量模型的预测误差，正则化项控制模型的复杂性。目标函数的目标是最小化损失函数和正则化项的和。
模型训练步骤：XGBoost的模型训练步骤如下：

初始化模型：将初始预测值设为全局平均值（或使用先验概率）。
迭代训练：对于每一轮迭代：
- 计算残差：计算当前模型的预测值与实际值之间的残差。
- 训练一个决策树：使用残差作为目标值，训练一个决策树模型。决策树的构建过程采用贪婪算法，通过优化目标函数来确定最佳的分裂点。
- 更新模型：根据学习率（步长）和决策树的输出，更新模型的预测值。
- 正则化：根据正则化项对树的结构进行修剪，以控制模型的复杂性。
重复步骤b，直到达到预先设定的迭代次数（即树的数量）或达到停止条件。

预测步骤：使用训练好的模型进行预测时，将待预测样本输入到每棵树中，根据树的结构和叶子节点的预测值得到最终的预测结果。最终的预测结果由所有树的预测值加权求和得到。

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42


import pandas as pd
import xgboost as xgb
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import accuracy_score

# 读取数据集
df = pd.read_excel('~/数模/2020C/附件/指标选择.xlsx', sheet_name='第一题')
X = df.iloc[:, 1:6]
y = df.iloc[:, 0]

# 划分数据集为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 将数据集转换为XGBoost专用的DMatrix格式
dtrain = xgb.DMatrix(X_train, label=y_train)

# 定义初始参数范围
param_grid = {
    'max_depth': [3, 4, 5], # 树的最大深度
    'eta': [0.1, 0.01, 0.001], # 学习率
    'subsample': [0.8, 0.9],
    'colsample_bytree': [0.8, 0.9],
}

# 创建模型
model = xgb.XGBClassifier(objective='multi:softmax', num_class=4) # 多分类问题，类别数

# 定义网格搜索对象
grid_search = GridSearchCV(estimator=model, param_grid=param_grid, cv=5)

# 在训练集上进行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数组合
print("最佳参数组合：", grid_search.best_params_)

# 使用最佳参数的模型进行预测
y_pred = grid_search.predict(X_test)

# 计算模型在测试集上的准确性
accuracy = accuracy_score(y_test, y_pred)
print("模型在测试集上的准确性：", accuracy)

Adaboost (Adaptive Boosting)

原理

AdaBoost是一种集成学习算法，用于提高分类器的准确性。它通过迭代地训练一系列弱分类器（即在某个特定任务上表现略好于随机猜测的分类器），并将它们组合成一个强分类器。

下面是AdaBoost算法的基本原理：

初始化样本权重：对于包含N个样本的训练集，初始时每个样本的权重相等，即1/N。
迭代训练弱分类器：对于每个迭代轮次t（从1到T），进行以下步骤：

使用当前样本权重训练一个弱分类器（例如决策树桩），得到分类器ht。
计算分类器ht在训练集上的误差率εt，即被错误分类的样本的权重之和。
计算分类器ht的权重系数αt，αt = 0.5 * ln((1-εt)/εt)。这个权重系数表示分类器ht在最终的强分类器中的重要程度，当分类器的误差率较低时，其权重系数会更大。
更新样本权重：根据样本是否被分类器ht正确分类来更新样本的权重，被错误分类的样本的权重会增加，而被正确分类的样本的权重会减少。

构建强分类器：将每个弱分类器的预测结果乘以其权重系数αt，并将它们线性组合得到最终的强分类器。
预测：使用强分类器对新样本进行预测。对于二分类问题，预测结果是通过对弱分类器的预测结果进行加权投票得到的。

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53


import pandas as pd
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.ensemble import AdaBoostClassifier
from sklearn.metrics import accuracy_score

# 读取数据集
df = pd.read_excel('~/数模/2020C/附件/指标选择.xlsx', sheet_name='第一题')
X = df.iloc[:, 1:6]
y = df.iloc[:, 0]

# 划分数据集为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=36)

# 创建AdaBoost模型
ada_model = AdaBoostClassifier()

# 定义参数网格
param_grid = {
    'n_estimators': [50, 100, 150],  # 弱分类器数量
    'learning_rate': [0.1, 0.01, 0.001]  # 学习率
}

# 创建GridSearchCV对象
grid_search = GridSearchCV(estimator=ada_model, param_grid=param_grid, cv=5)

# 在训练集上进行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数组合
print("最佳参数组合：", grid_search.best_params_)

# 使用最佳参数的模型进行预测
y_pred = grid_search.predict(X_test)

# 计算模型在测试集上的准确性
accuracy = accuracy_score(y_test, y_pred)
print("模型在测试集上的准确性：", accuracy)

# 针对数据量较小的模型
'''
# 创建AdaBoost模型
ada_model = AdaBoostClassifier(n_estimators=50, learning_rate=1.0)

# 在训练集上拟合AdaBoost模型
ada_model.fit(X_train, y_train)

# 使用训练好的模型进行预测
y_pred = ada_model.predict(X_test)

# 计算模型在测试集上的准确性
accuracy = accuracy_score(y_test, y_pred)
print("AdaBoost模型在测试集上的准确性：", accuracy)
'''

模拟退火 (Simulated Annealing, SA)

数学原理和过程 (默认求最小值)

令$T = T_0$ 表示开始退火的初始温度，随机产生一个初始解$x_t$,并计算对应的目标函数值
令$T= kT, k \in (0,1)$, 为温度下降速率
对当前解施加随机扰动，在其邻域内产生一个新解$x_{t+1}$，通过以下公式(MetroPolis 准则)判断是否接受这个解的概率

$$ \begin{aligned} P &= \begin{cases} 1 & E_{t+1} \le E_t \\ e^{\dfrac{-(E_{t+1} - E_t)}{kT}} & E_{t+1} \gt E_t \end{cases} \end{aligned} $$

流程图

应用

用于VLSI(超大集成电路设计)

双退火代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13


from scipy.optimize import dual_annealing

'''
scipy.optimize.dual_annealing(func, bounds, args=(), maxiter=1000, minimizer_kwargs=None, initial_temp=5230.0, restart_temp_ratio=2e-05, visit=2.62, accept=-5.0, maxfun=10000000.0, seed=None, no_local_search=False, callback=None, x0=None)

func(x, *args): args是Optional

bounds: (min,max)

args: optional

maxiter: 最大迭代次数, optional
'''

Regression

多元逻辑回归(Multinomial Logistic Regression)

原理

在多元逻辑回归中，我们有一个因变量（多类别）和多个自变量。假设有 k 个类别，我们需要建立 k-1 个二元逻辑回归模型来预测每个类别相对于某个基准类别的概率。
数学上，多元逻辑回归使用 softmax 函数将线性组合的结果转化为概率分布。对于第 i 个类别（其中 i = 1, 2, …, k-1），我们定义一个线性组合：

$$ z_i = b_{i0} + b_{i1}x_1 + b_{i2}x_2 + … + b_{ip}x_p $$

通过应用 softmax 函数，我们可以计算每个类别的概率。对于第$i$个类别，其概率为：

$$ P(Y=i|X) = exp(z_i) / (exp(z_1) + exp(z_2) + … + exp(z_{k-1})) $$

意义

训练多元逻辑回归模型的目标是最大化似然函数或最小化交叉熵损失函数。通常使用梯度下降等优化算法来寻找最优的系数。

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35


import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression

# 提取数据
data1 = pd.read_excel('~/2020C/附件/指标选择.xlsx', sheet_name='第一题')
data2 = pd.read_excel('~/2020C/附件/指标选择.xlsx', sheet_name='第二题')

# 获取 训练集 预测集
X_train = data1.iloc[:, 1:6].to_numpy()
y_train = data1.iloc[:, 0].to_numpy()  # 将目标变量转换为整数类型
X_test = data2.iloc[:, 0:5].to_numpy()

# 创建逻辑回归模型
model = LogisticRegression()

# 训练模型
model.fit(X_train, y_train)

# 使用模型进行预测
predictions = model.predict(X_test)

# 企业代号
code = ['']*302

for i in range(1, 303):
    code[i-1] += 'E' + str(i)

# 导出数据
data_out = pd.DataFrame({
    '企业代号': code,
    '企业信誉评级': predictions
    })

data_out.to_excel('第二题信誉评级.xlsx', index = False)

Ridge Regression

原理

减小Overfitting

LASSO Regression

Classification

Bayes Algorithms

$$ \begin{gather} \begin{aligned} P(c|X) &= \frac{ P(x|c) P(c)}{ P(X) } \\

\end{aligned} \end{gather} $$

Other Algorithms

贪心算法

一种可能的贪心算法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26


def chosen_suppliers_II(suppliers_A_names, suppliers_A, suppliers_B_names, suppliers_B, suppliers_C_names, suppliers_C):
    
    # 将A、B、C的供应商名称和供货量合并到一个列表中
    suppliers = list(zip(suppliers_A_names, suppliers_A, ['A'] * len(suppliers_A))) + list(zip(suppliers_B_names, suppliers_B, ['B'] * len(suppliers_B))) + list(zip(suppliers_C_names, suppliers_C, ['C'] * len(suppliers_C)))

    # 先按照A、B、C, 再按照供货量从大到小对suppliers列表进行排序
    suppliers = sorted(suppliers, key=lambda x: (x[2], -x[1]))

    # 按照A、B、C的顺序选取供应商，并购买其提供的原材料，直到所需量达到要求为止
    total_amount = 0
    total_productivity = 0
    purchased = []
    for name, supply, type in suppliers:
        if supply == 0: # 跳过供应量为0的供应商
            continue
        if total_productivity >= 28200:
            break
        purchased_amount = supply
        total_amount += purchased_amount
        if type == 'A' :
            total_productivity += purchased_amount/0.6
        elif type == 'B':
            total_productivity += purchased_amount/0.66
        else:
            total_productivity += purchased_amount/0.72
        purchased.append((name, purchased_amount))

动态规划算法

0-1背包问题

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29


def knapsack(weights, values, capacity):
    n = len(weights)
    # dp[i][j]: 在前 i 个物品中，背包容量为 j 时的最大价值
    dp = [[0] * (capacity + 1) for _ in range(n + 1)] 

    for i in range(1, n + 1):
        for j in range(1, capacity + 1):
            if weights[i - 1] <= j:
                dp[i][j] = max(dp[i - 1][j], values[i - 1] + dp[i - 1][j - weights[i - 1]])
            else:
                dp[i][j] = dp[i - 1][j]

    selected_items = []
    i, j = n, capacity
    while i > 0 and j > 0:
        if dp[i][j] != dp[i - 1][j]:
            selected_items.append(i)
            j -= weights[i - 1]
        i -= 1

    selected_items.reverse()
    return dp[n][capacity], selected_items

# 示例用法
weights = [2, 3, 4, 5, 10, 23, 1, 2, 34, 242, 23, 3]
values = [3, 4, 5, 6, 21, 23, 23, 51, 14, 41, 13, 14]
capacity = 100
max_value, selected_items = knapsack(weights, values, capacity)
print(f"选择第{selected_items}个物体可以获得最大价值为{max_value}")

调度

时间片(轮转)

需要的已知数据

time_slice: 时间片的长度
time_remaining: 该任务还需要多长时间完成
time_available: 当前时间片内的剩余时间

代码

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16


# 时间片的分配
while time_remaining > 0
    # 订单可以在当前时间切片内完成 (更新当前时间和剩下时间)
    if time_remaining <= time_available:
        current_time += time_remaining
        time_remaining = 0
        time_available -= time_remaining
        # 计算超时时间
        order_timeout = max(current_time - task['deadline'], 0)
        line_timeout += order_timeout
        
    # 订单无法在当前时间切片内完成，继续下一个时间切片
    else:
        current_time += time_available
        time_remaining -= time_available
        time_available = time_slice

EDF + Dynamic Programming + Linear Programming

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200


import pandas as pd
import numpy as np
import pulp

class ProductLine:
    def __init__(self, line_id=0, num_workstations=0, total_time=0, orders = []):
        # 属性: 产线编号、订单信息、工位数、总时间
        self.line_id = line_id
        self.orders = orders
        self.num_workstations = num_workstations
        self.total_time = total_time

    def add_orders(self, item):
        self.orders.append(item)

    def __str__(self):
        orders_str = "\n".join(str(order) for order in self.orders)
        return f"L{self.line_id}, 工位: {self.num_workstations}, 总耗时: {self.total_time}, 订单:\n{orders_str}"

class Order:
    def __init__(self, order_id, product_name, product_number, ddl, line_id, start_time = 0, finish_time = 0, overtime = 0, total_time = 0, standard_time = []):
        # 属性: 订单号 产品名称 产品数目 截止时间 对应产线 开始时间 结束时间 超时时间 总用时 产品标准工时
        self.order_id = order_id
        self.product_name= product_name
        self.product_number = product_number
        self.standard_time = standard_time
        self.ddl = ddl
        self.line_id =line_id
        self.start_time = start_time
        self.finish_time = finish_time
        self.overtime = overtime
        self.total_time = total_time
    
    def add_standard_time(self, item):
        self.standard_time.append(item)

    def __str__(self):
        return f"订单号: {self.order_id}, 总用时: {self.total_time}, 超时时间:{self.overtime}"

def read_data(file_path):
    df = pd.read_excel(file_path, sheet_name = "Sheet1")
    orders = df["生产订单"].to_numpy()
    products = df["产品名称"].to_numpy()
    standard_time = df["标准工时"].to_numpy()
    ddls = df["完工截止时间"].to_numpy()
    product_lines = df["产线代码"].to_numpy()

    tasks = [] # 定义字典，合并、存储所有订单
    current_order = None

    for i in range(len(orders)):

        # 新订单
        if orders[i] != orders[i-1]:
            # 添加上一个订单到字典
            if i != 0:
                tasks.append(new_order)
            # 添加新的订单项
            current_order = orders[i]
            new_order = Order(order_id=orders[i], product_name=products[i], product_number=1, ddl = ddls[i], line_id=product_lines[i], start_time = 0, finish_time= 0, overtime= 0, total_time= 0, standard_time=[])
            new_order.add_standard_time(standard_time[i])
        # 合并相同订单的项
        else:
            new_order.product_number += 1
            new_order.add_standard_time(standard_time[i])

    return tasks

def save_data(file_path, lines):
    data_out = pd.DataFrame({
        "生产订单": [order.order_id for line in lines for order in line.orders],
        "开始时间": [order.start_time for line in lines for order in line.orders],
        "结束时间": [order.finish_time for line in lines for order in line.orders],
        "超时时间": [order.overtime for line in lines for order in line.orders]
        })
    
    data_out.to_excel(file_path, sheet_name="Sheet1", index=False)

def min_processing_time(products, num_stations):
    if num_stations >= len(products):
        return max(products)
    elif num_stations == 1:
        return sum(products)

    m = len(products)
    n = num_stations
    dp = [[0] * (n + 1) for _ in range(m + 1)]

    # 初始化边界状态
    for i in range(1, m + 1):
        dp[i][1] = sum(products[:i])

    # 通过状态转移方程求解问题
    for i in range(1, m + 1):
        for j in range(2, n + 1):
            dp[i][j] = float('inf')
            for k in range(j - 1, i):
                dp[i][j] = min(dp[i][j], max(dp[k][j - 1], sum(products[k:i])))

    return dp[m][n]

def Linear_Programming(num_workers, line, station_number): # 员工数 生产线数 工位分布
    skill_matrix = pd.read_excel('~/杉树杯/附件/附件4.xlsx', usecols='B:M', skiprows=2, header=None).value
    
    # 将数据转换为NumPy数组并进行值转换
    skill_matrix[skill_matrix == 'N'] = 0
    skill_matrix[skill_matrix == 'O'] = 0.8
    skill_matrix[skill_matrix == 'E'] = 1

    # 构建线性规划模型
    prob = pulp.LpProblem('Product_Allocation', pulp.LpMaximize)
    
    # 定义决策变量 - 每个员工被分配到每个工位的binary变量
    assignment = [(i,j) for i in range(num_workers) for j in range(line)]
    vars = pulp.LpVariable.dicts('assignment', assignment, cat='Binary')
    
    # 定义目标函数 - 最大化效率
    prob += pulp.lpSum([vars[(i,j)] * skill_matrix[i][j] for i in range(num_workers) for j in range(line)]), 'Efficiency'
    
    # 添加约束条件
    for i in range(num_workers):
        prob += pulp.lpSum([vars[(i,j)] for j in range(line)]) <= 1  # 每个员工只能被分配到一个工位
        
    for j in range(line):
        prob += pulp.lpSum([vars[(i,j)] for i in range(num_workers)]) <= station_number[j] # 每个工位只能分配一个员工
    # 求解模型
    prob.solve() 
    
    # 提取解并存储在矩阵中
    coordination = []
    for v in prob.variables():
        if v.varValue == 1:
            # 提取坐标
            coord = v.name.replace('assignment_', '')
            coord = coord.replace(',_', ',')
            coordination.append(coord)
    
    # (员工索引，生产线索引) 工作人数
    coordination = [tuple(int(num) for num in element.strip('()').split(',')) for element in coordination]
    working_number = len(coordination)
    
    # 通过员工分配矩阵计算效率矩阵
    arrangment = np.zeros((num_workers, line))
    efficiency_matrix = np.zeros((num_workers,line))
    
    for i in range(working_number):
        arrangment[coordination[i][0]][coordination[i][1]] = 1

    efficiency_matrix = arrangment*skill_matrix

    # 根据效率矩阵计算效率 根据员工分配矩阵获得员工分配列表
    efficiency = [[] for _ in range(line)]
    arrangment_list = [[] for _ in range(line)]
    
    for x in range(num_workers): # 遍历员工
        for y in range(line): # 遍历生产线
    
            if efficiency_matrix[x][y] != 0: # 若是非零元素
                efficiency[y].append(efficiency_matrix[x][y]) # 加入对应生产线的效率
                arrangment_list[y].append('PE' + "{:03d}".format(x)) # 加入员工分配列表
    
    for i in range(line):
        if (len(efficiency[i]) == 0):
            efficiency[i].append(0)

    # 返回效率列表，员工列表
    return efficiency, arrangment_list

def __main__():
    # 基本参数
    num_lines = 12
    num_workstations = [1,2,1,3,1,1,2,1,1,3,1,1]
    lines = [ProductLine(line_id = i+1, num_workstations=num_workstations[i], total_time= 0, orders= []) for i in range(num_lines)]

    # 读取数据
    tasks = read_data("~/数模/杉树杯/附件/附件1.xlsx")

    # 给每个生产线分配订单
    for task in tasks:
        lines[task.line_id - 1].add_orders(task)

    # 计算订单的时间信息
    line_total_overtime = [0]*num_lines
    for line in lines:
        for order in line.orders:
            order.start_time = line.total_time
            order.total_time = min_processing_time(order.standard_time, num_workstations[line.line_id - 1])
            order.finish_time = order.start_time + order.total_time
            line.total_time = order.finish_time
            order.overtime = max(0, order.finish_time -order.ddl)
            line_total_overtime[line.line_id - 1] += order.overtime
        print(f"{line.line_id}号生产线的超时时间为{line_total_overtime[line.line_id-1]}")

    print(f"第一题的超时时间为: {sum(line_total_overtime)} 分钟")

    # 存储数据
    save_data("第一题分配方案.xlsx", lines)

if __name__ == "__main__":
   __main__()

SRIMAX(滑动平均的改良版)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39


import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from matplotlib.font_manager import FontProperties

# 读取数据(可以修改)
data = pd.read_excel('Data2.xlsx', sheet_name='订货量', usecols='C:IH', skiprows=51, nrows=3, header=None)
data = data.T
data.index = pd.to_datetime(data.index)

# 设置中文字体
font = FontProperties(fname='/Users/downeyflyfan/Library/Fonts/NotoSansCJKsc-Regular.otf')

plt.figure(figsize=(10, 5))
for i, product in enumerate(['A', 'B', 'C']):

    # 进行差分处理
    diff = np.diff(data[i].values)
    
    # 拟合SARIMA模型
    model = SARIMAX(data[i], order=(4,1,4), seasonal_order=(4,1,4,24))  # (p, d, q)是非季节性部分的阶数，(P, D, Q)是季节性部分的阶数，m是季节性周期。
    results = model.fit()
    forecast = results.predict(start=len(data[i]), end=len(data)+23, dynamic=True).round(0)

    # 将预测值限制在大于0的范围内
    forecast = np.clip(forecast, a_min=0, a_max=None)
    
    # 绘制原始数据和预测数据
    plt.plot(data[i], label='原材料{}原始数据'.format(product))
    plt.plot(forecast, label='原材料{}预测数据'.format(product))

plt.xlabel('时间(W)', fontproperties=font, fontsize=20)
plt.ylabel('订货量($m^3$)', fontproperties=font, fontsize=20)
plt.title('订货量预测图像', fontproperties=font, fontsize=24)

plt.tight_layout()
plt.legend(prop=font)
plt.show()