2023 第九届数维杯B题 建模解析，小鹿学长带队指引全代码文章与思路

我是小鹿学长，就读于上海交通大学，截至目前已经帮200+人完成了建模与思路的构建的处理了～
这回带大家体验一下第九届数维杯B题呀！

问题重述

问题一：通过数学建模，分析催化剂（脱硫灰）在棉秆（cotton stalk）和模型化合物（CE和LG）的热解过程中对产物（焦油、水、焦渣、合成气）产量的影响，特别关注不同混合比例下脱硫灰在催化热解中的作用。

问题二：利用实验数据，深入研究三种热解组合（脱硫灰-棉秆、脱硫灰-CE、脱硫灰-LG）中混合比例对热解气体产物（H2、CO、CO2、CH4等）产量的影响，通过图形化结果进行解释。

问题三：在相同脱硫灰催化比例下，探讨纤维素（CE）和木质素（LG）的热解产物产量及热解气体组分的异同，提供详细解释。

问题四：建立脱硫灰对模型化合物（CE和LG）的催化反应机理模型，进行反应动力学分析，并通过最小二乘法和统计检验验证模型参数的显著性。

问题五：采用机器学习方法（支持向量回归等），基于给定数据建立模型，预测催化剂在不同条件下对产物产量的影响，实现对热解产物的数量级预测。

## 建模思路 
## 问题一

问题一的建模思路采用多元线性回归，其中混合比例是自变量，而产物的产量是因变量。以下是具体步骤：

1. **数据准备：**
  - 将实验数据整理成适合多元线性回归的格式，确保数据集中包含混合比例和对应的产物产量。

2. **变量选择：**
  - 选择混合比例作为自变量，即 \(X\)。产物产量作为因变量，即 \(Y\)。

3. **多元线性回归模型：**
  - 构建多元线性回归模型：\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_n X_n + \varepsilon \]
  - 这里，\(Y\) 是产物产量，\(X_1, X_2, \ldots, X_n\) 是混合比例的不同组分，\(\beta_0, \beta_1, \ldots, \beta_n\) 是回归系数，\(\varepsilon\) 是误差。

4. **模型训练：**
  - 将数据集拆分为训练集和测试集。
  - 使用训练集对多元线性回归模型进行训练。

5. **模型评估：**
  - 使用测试集数据对模型进行评估，考察模型的预测性能。
  - 可以使用均方误差（MSE）等指标来评估模型的拟合程度。

6. **回归系数的显著性检验：**
  - 对回归系数进行显著性检验，以判断混合比例对产物产量是否具有显著影响。
  - 通常，p-value 小于设定的显著性水平（如0.05）表明相关系数显著。

7. **模型解释：**
  - 解释回归系数的物理意义，了解混合比例对产物产量的具体影响。可能需要考虑交互作用等因素。

8. **模型应用：**
  - 使用训练好的模型对新的混合比例进行预测，以估计产物产量。

```python
import pandas as pd
import statsmodels.api as sm
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

# 假设你的数据文件为 'data.xlsx'，包含混合比例和产物产量
data = pd.read_excel('data.xlsx')

# 选择自变量和因变量
X = data[['混合比例1', '混合比例2', '其他混合比例', ...]] # 根据你的实际情况替换列名
y = data['产物产量']

# 添加截距项
X = sm.add_constant(X)

# 分割数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 构建多元线性回归模型
model = sm.OLS(y_train, X_train).fit()

# 模型预测
y_pred = model.predict(X_test)
#见完整代码

```

## 问题二

问题二的建模思路采用主成分分析（PCA），通过降维将多个热解气体产物转换为几个主成分，并对这些主成分进行回归分析。

主成分分析（Principal Component Analysis，简称PCA）是一种常用的数据降维技术，用于发现数据中的主要结构并减少数据的维度。PCA的主要目标是通过线性变换将原始数据映射到一个新的坐标系，使得映射后的数据在新坐标系下的方差最大化。

具体来说，PCA的步骤如下：

1. **数据标准化：**
  - 对原始数据进行标准化处理，使得每个特征的均值为0，方差为1。

2. **构造协方差矩阵：**
  - 计算标准化后的数据的协方差矩阵。协方差矩阵反映了不同特征之间的关系。

3. **计算特征值和特征向量：**
  - 对协方差矩阵进行特征值分解，得到特征值和对应的特征向量。特征向量表示新坐标系的方向，而特征值表示数据在这个方向上的方差大小。

4. **选择主成分：**
  - 将特征值按大小排序，选择前k个特征值对应的特征向量作为主成分，其中k是希望降维后的维度。

5. **构造投影矩阵：**
  - 将选定的前k个特征向量按列构成投影矩阵。

6. **数据投影：**
  - 将原始数据投影到选定的前k个特征向量所构成的子空间上，得到降维后的数据。

PCA的应用包括但不限于：
- 数据降维：减少数据的维度，去除冗余信息。
- 数据可视化：将高维数据映射到二维或三维空间，方便可视化。
- 噪声过滤：去除数据中的噪声。
- 特征提取：从高维数据中提取最重要的特征。
------
结合问题二，用PCA方法来解决的详细步骤如下：

1. **数据准备：**
  - 将实验数据整理成适合主成分分析的格式，确保数据集中包含热解气体产物的信息。

2. **主成分分析（PCA）：**
  - 使用主成分分析对热解气体产物进行降维，将多个相关的产物转换为几个主成分，以捕捉主要的变化。
  - 选择主成分的数量，通常选择能够解释大部分方差的主成分。

3. **回归模型：**
  - 使用主成分作为新的自变量进行回归分析。这里的因变量是混合比例。

4. **模型训练：**
  - 将数据集拆分为训练集和测试集。
  - 使用训练集对回归模型进行训练。

5. **模型评估：**
  - 使用测试集数据对模型进行评估，考察模型的预测性能。
  - 可以使用均方误差（MSE）等指标来评估模型的拟合程度。

6. **主成分的解释：**
  - 解释主成分的物理意义，了解它们对原始热解气体产物的贡献。

7. **模型解释：**
  - 解释回归系数的物理意义，了解主成分对混合比例的影响。

8. **模型应用：**
  - 使用训练好的模型对新的混合比例进行预测，以估计热解气体产物的变化趋势。

```python
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

# 读取数据
data = pd.read_excel('你的数据文件.xlsx')

# 选择自变量（混合比例）和因变量
X = data[['混合比例1', '混合比例2', '其他混合比例', ...]] # 根据你的实际情况替换列名
y = data['热解产物产量']

# 数据标准化（重要：PCA对数据的尺度敏感，需要先进行标准化）
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# PCA降维
pca = PCA(n_components=2) # 选择主成分数量
X_pca = pca.fit_transform(X_scaled)

# 构建新的数据框，包含主成分和因变量
df_pca = pd.DataFrame(X_pca, columns=['PC1', 'PC2'])
df_pca['热解产物产量'] = y

# 分割数据集
X_train, X_test, y_train, y_test = train_test_split(df_pca[['PC1', 'PC2']], y, test_size=0.2, random_state=42)

# 线性回归模型
model = LinearRegression()
model.fit(X_train, y_train)
```

## 问题三

1. **数据准备：**
  - 整理实验数据，按照时间序列组织，包括纤维素（CE）和木质素（LG）的热解产物产量以及热解气体组分的信息。确保数据集中包含时间戳，以便建立时间序列模型。

2. **数据标准化：**
  - 对数据进行标准化处理，确保每个时间点的数据在相同尺度上。

3. **构建RNN模型：**
  - 使用深度学习框架（如TensorFlow或PyTorch）构建循环神经网络模型。选择合适的网络结构，可以考虑使用LSTM（长短时记忆网络）或GRU（门控循环单元）等变种。
  - 输入序列为纤维素和木质素的热解产物产量以及热解气体组分的信息，输出为相应的产物产量。

4. **模型训练：**
  - 将数据集划分为训练集和测试集。
  - 使用训练集对RNN模型进行训练，调整模型参数。

5. **模型评估：**
  - 使用测试集对模型进行评估，考察模型的预测性能。
  - 可以使用均方误差（MSE）等指标来评估模型的拟合程度。

6. **结果解释：**
  - 分析模型输出，了解在相同脱硫灰催化比例下，纤维素和木质素的热解产物产量以及热解气体组分的差异。

7. **可视化：**
  - 使用可视化工具展示模型对纤维素和木质素热解过程的预测结果，以更直观地展示差异。

8. **调整和优化：**
  - 根据模型评估的结果，可能需要调整模型结构、超参数或采用其他技术手段来优化模型的性能。

```python
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense
from sklearn.metrics import mean_squared_error

# 读取数据
data = pd.read_excel('你的数据文件.xlsx')

# 选择特征和目标
features = data[['时间戳', 'CE产物产量', 'LG产物产量', 'CE热解气体组分', 'LG热解气体组分']]
target = data[['CE产物产量', 'LG产物产量', 'CE热解气体组分', 'LG热解气体组分']]

# 数据标准化
scaler = StandardScaler()
features_scaled = scaler.fit_transform(features)
target_scaled = scaler.fit_transform(target)

# 构建时间序列数据
time_steps = 10 # 设置时间步长，根据实际情况调整
X, y = [], []
for i in range(len(features_scaled) - time_steps):
  X.append(features_scaled[i:i + time_steps, :])
  y.append(target_scaled[i + time_steps, :])

X, y = np.array(X), np.array(y)

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

# 构建LSTM模型
model = Sequential()
model.add(LSTM(units=50, activation='relu', input_shape=(X_train.shape[1], X_train.shape[2])))
model.add(Dense(units=target.shape[1])) # 输出层的神经元数量等于目标的特征数量

# 编译模型
model.compile(optimizer='adam', loss='mse')

# 训练模型
model.fit(X_train, y_train, epochs=50, batch_size=32, validation_split=0.1)

# 模型预测
y_pred = model.predict(X_test)

# 反向转换标准化的预测结果
y_pred_rescaled = scaler.inverse_transform(y_pred)
y_test_rescaled = scaler.inverse_transform(y_test)

# 模型评估
mse = mean_squared_error(y_test_rescaled, y_pred_rescaled)
print(f'Mean Squared Error: {mse}')
```

## 问题四
对于问题四，使用简化的表观动力学模型（Apparent Kinetics Model）来描述脱硫灰对模型化合物（CE和LG）的催化反应。这个模型假设反应服从一级反应动力学，并通过Arrhenius方程考虑温度的影响。具体步骤如下：

1. **一级反应动力学模型：**
  - 表达式：$r = k \cdot [CE]$
  - $r$是反应速率，$CE$是CE的浓度，k是速率常数。

2. **Arrhenius方程：**
  - 将速率常数 k与温度联系起来：$k = A \cdot e^{-\frac{E}{RT}}$
  - A是预指数因子，E是活化能，R 是理想气体常数，T是温度。

3. **整体反应速率方程：**
  - 将上述两个方程结合，得到整体反应速率方程：$r = A \cdot e^{-\frac{E}{RT}} \cdot [CE]$

4. **参数估计：**
  - 利用实验数据，通过最小二乘法等方法，估计Arrhenius方程中的参数 A和E，以及一级反应动力学模型中的速率常数。

5. **模型验证：**
  - 使用估计的参数，将模型应用于其他实验条件，与实验数据进行比较，验证模型的预测性能。

6. **灵敏度分析：**
  - 进行参数的灵敏度分析，了解模型对参数变化的敏感性。这有助于确定哪些参数对于模型的输出影响较大。

7. **模型解释和物理意义：**
  - 解释模型中的参数的物理意义，理解模型对实验观察到的现象的贡献。验证模型是否与实验室观察到的反应机理相符。

8. **优化模型：**
  - 根据验证和解释的结果，对模型进行优化。可能需要调整动力学模型的形式，进一步细化反应机理，或者考虑其他因素的影响。

```python
import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 假设的一级反应动力学模型
def first_order_kinetics(C, A, E, T):
  R = 8.314 # 理想气体常数
  k = A * np.exp(-E / (R * T))
  return k * C

# 生成一些模拟数据
CE_data = np.array([concentration1, concentration2, ...]) # 你的CE浓度数据
rate_data = np.array([rate1, rate2, ...]) # 你的反应速率数据

# 参数估计的初始值
initial_guess = [1.0, 1.0, 300.0] # A, E, 初始温度

# 利用最小二乘法进行参数估计
params, covariance = curve_fit(first_order_kinetics, CE_data, rate_data, p0=initial_guess)

# 获取估计得到的参数
A_estimate, E_estimate, T_estimate = params

# 预测模型
predicted_rates = first_order_kinetics(CE_data, A_estimate, E_estimate, T_estimate)

# 绘制拟合结果
plt.scatter(CE_data, rate_data, label='实验数据')
```

## 问题五
此处使用支持向量机（SVM）解决问题五，采用回归问题的方法，即预测热解产物的产量或数量。以下是详细的建模思路：

1. **数据准备：**
  - 整理实验数据，包括各种条件下的热解产物的产量或数量，以及相关的实验条件。确保数据中包含特征（如温度、混合比例等）和目标变量（热解产物的产量或数量）。

2. **特征选择和标准化：**
  - 选择与热解产物产量或数量相关的特征，并对这些特征进行标准化。确保特征在相同尺度上，以提高支持向量机的性能。

3. **数据划分：**
  - 将数据集分为训练集和测试集，通常采用交叉验证的方式。

4. **支持向量机回归模型选择：**
  - 选择合适的支持向量机回归模型。对于回归问题，通常使用核函数为径向基函数（RBF）的支持向量机。根据数据特点选择合适的超参数。

5. **模型训练：**
  - 使用训练集对支持向量机回归模型进行训练。

6. **模型评估：**
  - 使用测试集对模型进行评估，通常使用均方误差（MSE）等指标来评估预测精度。

7. **调整和优化：**
  - 根据评估结果，调整支持向量机的超参数，如调整核函数的参数，以优化模型性能。

8. **预测：**
  - 在有限的数据条件下，使用训练好的支持向量机回归模型进行预测。输入实验条件，模型将给出相应的热解产物产量或数量的预测值。

9. **模型解释：**
  - 解释模型的预测结果，了解模型认为哪些条件对于热解产物的产量或数量有较大影响。支持向量机虽然是一种黑盒模型，但通过一些方法可以进行特征重要性分析。

10. **迭代和改进：**
  - 如果有新的实验数据可用，可以将其用于进一步改进和迭代模型，以提高预测的准确性。

以下是使用Python的`sklearn`库进行支持向量机回归的简化代码示例，假设数据集包括特征`X`和目标变量`y`：

```python
from sklearn.svm import SVR
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from sklearn.preprocessing import StandardScaler

# 特征标准化
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# 数据划分
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.2, random_state=42)

# 创建支持向量机回归模型
svm_model = SVR(kernel='rbf', C=1.0, epsilon=0.1)

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

# 模型预测
y_pred = svm_model.predict(X_test)

# 评估模型性能
mse = mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse}')
```