2024辽宁省数学建模竞赛C题思路代码分析：基于响应面分析法

2024辽宁省数学建模竞赛ABC题完整成品论文和全部问题的解题代码更新
https://www.yuque.com/u42168770/qv6z0d/midyg0c9wypw9s5m

摘要

本文通过数学建模方法，系统地研究了2024辽宁省数学建模竞赛C题改性生物碳对水中砷离子[As(V)]和洛克沙胂（ROX）的吸附问题。研究围绕三个核心问题展开：探讨反应温度、溶液pH和吸附剂用量对As(V)和ROX去除率的影响；优化As(V)和ROX的总吸附量；设计额外的实验以进一步改进模型。通过建立多元非线性响应面回归模型、多目标粒子群优化模型和多目标自适应实验设计优化模型，我们成功地解析了吸附过程的复杂机制，优化了吸附条件，并提出了高效的实验设计方案。

针对问题一，我们构建了多元非线性响应面回归模型来分析各因素对As(V)和ROX去除率的影响。模型采用了二阶多项式形式，包含主效应、二次效应和交互效应。为提高模型的鲁棒性和解释性，引入了LASSO正则化。模型的求解采用了高阶迭代优化算法，结合了坐标下降法和近似牛顿法。求解结果显示，对于ROX，(求解结果略)。对于As(V)，(求解结果略)。模型的创新点在于综合考虑了多个因素的非线性和交互作用，并通过正则化技术提高了模型的泛化能力。

问题二的解决采用了响应面多目标粒子群优化模型（RSMPSOM）。该模型结合了响应面法的建模能力、多目标优化的灵活性和粒子群算法的全局搜索能力。我们首先建立了As(V)和ROX吸附量的二阶多项式模型，然后构建了一个加权和目标函数来平衡两种污染物的吸附量。优化过程使用了自适应多目标粒子群优化算法（AMPSOA），该算法引入了自适应参数调整和多目标权重自适应机制。求解结果得到的最优条件为：(求解结果略)。模型的创新点在于引入了自适应机制，能够在优化过程中动态平衡不同目标，提高了算法的鲁棒性和效率。

对于问题三，我们提出了一种多目标自适应实验设计优化模型（MOAEDOM）。该模型综合考虑了空间填充、模型不确定性、预测优化和多目标平衡等多个方面。模型的核心是一个自适应多准则实验设计算法（AMCEDA），该算法通过迭代优化来选择新的实验点。算法使用高斯过程回归来构建响应面模型，并引入了自适应权重调整机制。求解结果给出了5个新的实验点，(求解结果略)。这一结果暗示了在这些条件下可能存在更有趣或更重要的吸附行为。模型的创新点在于其自适应性和多准则决策机制，能够在探索新区域和优化已知区域之间取得平衡。
（后略）

关键词：多元非线性响应面回归、粒子群优化、高斯过程回归、自适应实验设计、As(V)和ROX吸附、水体污染处理

问题重述

目前全球水体砷污染问题严重，极大地危害了人类的健康。砷在水体中主要以砷离子[As(V)]和洛克沙胂（ROX）两种形式存在。改性生物碳具有优秀的吸附性能，可用于砷污染的处理。改性生物碳表面的官能团可以通过静电吸引、配位作用或络合作用与As(V)和ROX 结合，从而使As(V)和ROX 吸附在生物炭表面，达到降低环境中As(V)和ROX 浓度的目的。改性生物碳吸附净化污水示意图如下图所示。

在吸附过程中，反应温度，溶液pH，吸附剂用量对As(V)和ROX 的去除率和吸附量将产生影响（名词解释见附录）。因此，探索改性生物碳高效吸附有机无机砷的工艺条件，对于改善水体砷污染问题具有重要的意义和价值。
某化学实验室针对不同吸附条件进行了一系列实验，结果如附件1 所示。请通过数学建模完成下列问题：

（1）当其他反应条件一定时，分别探讨反应温度，溶液pH，吸附剂用量对As(V)和ROX 去除率的影响。

（2）如何选择反应温度，溶液pH，吸附剂用量让As(V)和ROX 的总吸附量尽可能大。

（3）如果允许再增加5 次实验，应如何设计，并给出详细理由。

问题分析

2024辽宁省数学建模竞赛C题这道题目探讨了改性生物碳对水中洛克沙胂(ROX)和砷离子[As(V)]的吸附问题。题目的核心是研究不同条件(反应温度、溶液pH、吸附剂用量)对吸附效果的影响,并优化这些条件以达到最佳吸附效果。这是一个典型的多变量优化问题,涉及数据分析、模型构建和优化算法的应用。

对于问题(1),我们需要分别分析反应温度、溶液pH和吸附剂用量对ROX和As(V)去除率的影响。这需要对给定的实验数据进行深入分析,可以考虑使用多元回归分析或响应面法来建立各个因素与去除率之间的关系模型。通过这些模型,我们可以量化每个因素对去除率的影响程度,并识别出最显著的影响因素。此外,还可以使用方差分析(ANOVA)来评估各因素的显著性,或者采用主成分分析(PCA)来降维并找出主要影响因素。在可视化方面,可以使用热力图或三维曲面图来直观展示各因素对去除率的影响。

问题(2)要求我们找出能够最大化As(V)和ROX总吸附量的最佳条件组合。这是一个典型的优化问题,可以考虑使用响应面法(RSM)来构建总吸附量与各因素之间的关系模型,然后利用优化算法(如梯度下降法、遗传算法或粒子群优化算法)来寻找最优解。另一种方法是使用多目标优化技术,将As(V)和ROX的吸附量作为两个目标函数,寻找帕累托最优解集。此外,也可以考虑使用机器学习方法,如支持向量回归(SVR)或随机森林回归,来建立预测模型,然后通过网格搜索或贝叶斯优化来找到最佳参数组合。

问题(3)涉及实验设计,目的是通过增加5次实验来获得更多有价值的信息。这里可以考虑使用正交实验设计或中心复合设计(CCD)来安排新的实验点,以覆盖更广泛的参数空间。另一种方法是使用基于模型的实验设计,如最大信息设计或D-最优设计,这些方法可以基于已有的模型预测,选择最具信息量的新实验点。此外,还可以考虑使用自适应实验设计方法,根据已有实验结果动态调整新实验的参数设置。在选择新的实验点时,应考虑填补现有数据的空白,探索边界条件,以及验证模型在特定区域的准确性。

模型假设

以下是对问题1-问题3的模型建立与求解过程中使用的模型假设的总结：

（由于篇幅问题，假设略）

符号说明

以下是问题1-问题3的模型建立与求解过程中使用的符号及其说明的总结表格：（由于篇幅问题，略）

模型的建立与求解

问题1模型的分析与建立

思路分析

2024辽宁省数学建模竞赛C问题1要求我们探讨反应温度、溶液pH和吸附剂用量对As(V)和ROX去除率的影响。这是一个典型的多因素影响分析问题，需要我们深入研究实验数据，建立合适的数学模型来描述各因素与去除率之间的关系。考虑到问题的复杂性和数据的特征，我们可以采用多元非线性回归分析结合响应面法来构建一个综合性的数学模型。这种方法不仅能够捕捉各因素对去除率的非线性影响，还能够揭示因素之间的交互作用，为后续的优化提供理论基础。

在模型构建过程中，我们首先需要对原始数据进行预处理，包括数据清洗、异常值检测和必要的标准化处理。然后，我们将使用多项式回归的方法来拟合各因素与去除率之间的关系，同时考虑因素间的交互项。为了避免过拟合问题，我们将采用正则化技术，特别是LASSO（Least Absolute Shrinkage and Selection Operator）正则化，以实现特征选择和模型简化。此外，我们还将利用响应面法来可视化模型结果，直观地展示各因素对去除率的影响。

在模型评估阶段，我们将使用交叉验证技术来评估模型的泛化能力，并通过计算均方根误差（RMSE）、决定系数（R²）等指标来量化模型的性能。最后，我们将基于建立的模型，对各因素的主效应和交互效应进行详细分析，得出As(V)和ROX去除率随各因素变化的规律。

多元非线性响应面回归模型建立

为了全面分析反应温度、溶液pH和吸附剂用量对As(V)和ROX去除率的影响，我们提出了一个多元非线性响应面回归模型。这个模型不仅考虑了各因素的主效应，还包含了它们之间的交互作用，以及可能存在的二次效应。模型的一般形式可以表示为：

$$Y={\beta }_0+\sum _{i=1}^k{\beta }_i{X}_i+\sum _{i=1}^k{\beta }_{ii}{X}_i^2+\sum _{i<j}^k{\beta }_{ij}{X}_i{X}_j+\epsilon$$

其中，$Y$ 表示响应变量（As(V)或ROX的去除率），$X_i$ 和 $X_j$ 表示自变量（反应温度、溶液pH和吸附剂用量），${\beta }_0$ 是截距项，${\beta }_i$ 是一次项系数，${\beta }_{ii}$ 是二次项系数，${\beta }_{ij}$ 是交互项系数，$\epsilon$ 是随机误差项，$k$ 是自变量的数量（在本问题中 $k=3$）。

具体到本问题，我们可以将模型展开为：

$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+{\beta }_3{X}_3+{\beta }_{11}{X}_1^2+{\beta }_{22}{X}_2^2+{\beta }_{33}{X}_3^2+{\beta }_{12}{X}_1{X}_2+{\beta }_{13}{X}_1{X}_3+{\beta }_{23}{X}_2{X}_3+\epsilon$$

其中，$X_1$ 表示反应温度，$X_2$ 表示溶液pH，$X_3$ 表示吸附剂用量。

为了提高模型的鲁棒性和解释性，我们引入LASSO正则化项，修改后的目标函数变为：
（略，见完整版）

高阶迭代优化算法步骤

为了有效地估计模型参数并实现模型的优化，我们提出了一种高阶迭代优化算法。该算法结合了坐标下降法和近似牛顿法，能够高效地处理带有L1正则化项的非线性优化问题。算法的主要步骤如下：

1. 数据预处理：对原始数据进行标准化处理，使所有变量的尺度一致。标准化公式为：

   $$X_{normalized}=\frac{X-\mu }{\sigma }$$

   其中，$\mu$ 是变量的均值，$\sigma$ 是标准差。

2. 初始化：随机初始化模型参数 $\beta$，设定正则化参数 $\lambda$ 和收敛阈值 $ϵ$。

3. 外层循环：(略)

4. 内层循环：（略）

5. 响应面分析：基于估计的模型参数，构建响应面，并分析各因素对去除率的影响。

算法通过结合坐标下降法和近似牛顿法，既能有效处理L1正则化带来的非平滑优化问题，又能快速收敛于非线性模型的最优解。坐标下降法在处理L1正则化项时特别有效，因为它可以直接利用软阈值函数进行参数更新。而近似牛顿法则可以加速非线性优化的收敛速度，特别是在处理交互项和二次项时。（后略）

问题一模型的求解

预处理2024辽宁省数学建模竞赛C题数据：变成二维表格，方便代码直接加载数据

以下是基于前面建立的多元非线性响应面回归模型，用于求解问题1的详细Python代码，包括数据处理、模型训练、结果可视化和分析：（完整代码见附件）

# 分离特征和目标变量
X = data[['吸附剂用量(g/L)', '温度', 'PH']]
y_ROX = data['ROX去除率(%)']
y_As = data['As(V)去除率(%)']

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

# 添加二次项和交互项
X_poly = np.column_stack((X_scaled, X_scaled**2, 
                          X_scaled[:, 0] * X_scaled[:, 1],
                          X_scaled[:, 0] * X_scaled[:, 2],
                          X_scaled[:, 1] * X_scaled[:, 2]))

# 定义模型函数
def model(X, beta):
    return np.dot(X, beta)

# 定义目标函数（包含L1正则化）
return np.mean((y - model(X, beta))**2) + lambda_ * np.sum(np.abs(beta))

# 定义梯度函数
def gradient(beta, X, y, lambda_):
    n = X.shape[0]
    return -2/n * np.dot(X.T, (y - model(X, beta))) + lambda_ * np.sign(beta)

# 定义优化函数
def optimize_model(X, y, lambda_=0.1, max_iter=1000, tol=1e-6):
    beta_init = np.zeros(X.shape[1])
    result = minimize(objective, beta_init, args=(X, y, lambda_),
                      method='L-BFGS-B', jac=gradient, options={'maxiter': max_iter, 'gtol': tol})
    return result.x

# 训练模型
（略）
# 计算R²和RMSE
y_pred_ROX = model(X_poly, beta_ROX)
y_pred_As = model(X_poly, beta_As)

r2_ROX = r2_score(y_ROX, y_pred_ROX)
rmse_ROX = np.sqrt(mean_squared_error(y_ROX, y_pred_ROX))
r2_As = r2_score(y_As, y_pred_As)
rmse_As = np.sqrt(mean_squared_error(y_As, y_pred_As))

# 可视化：响应面图
def plot_response_surface(X, y, beta, title, filename):
    fig = plt.figure(figsize=(18, 5))
    
    for i in range(3):
        ax = fig.add_subplot(1, 3, i+1, projection='3d')
        x1, x2 = np.meshgrid(np.linspace(-2, 2, 100), np.linspace(-2, 2, 100))
        X_plot = np.column_stack((x1.ravel(), x2.ravel(), np.zeros_like(x1.ravel())))
        if i == 0:
            X_plot[:, 2] = X[:, 2].mean()
            xlabel, ylabel = '吸附剂用量', '温度'
        elif i == 1:
            X_plot[:, 1] = X[:, 1].mean()
            xlabel, ylabel = '吸附剂用量', 'PH'
        else:
            X_plot[:, 0] = X[:, 0].mean()
            xlabel, ylabel = '温度', 'PH'
        
        X_plot_poly = np.column_stack((X_plot, X_plot**2, 
                                       X_plot[:, 0] * X_plot[:, 1],
                                       X_plot[:, 0] * X_plot[:, 2],
                                       X_plot[:, 1] * X_plot[:, 2]))
        z = model(X_plot_poly, beta).reshape(x1.shape)
        
        surf = ax.plot_surface(x1, x2, z, cmap='viridis')
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_zlabel('去除率')
        fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)
    
    plt.suptitle(title, fontsize=16)
    plt.tight_layout()
    plt.savefig(f'问题1_{filename}.png', dpi=300, bbox_inches='tight')
    plt.show()

模型性能：

ROX模型摘要:
                            OLS Regression Results                            
==============================================================================
Dep. Variable:              ROX去除率(%)   R-squared:                       0.776
Model:                            OLS   Adj. R-squared:                  0.758
Method:                 Least Squares   F-statistic:                     44.14
Date:                Sat, 13 Jul 2024   Prob (F-statistic):           2.92e-33
Time:                        13:08:36   Log-Likelihood:                -498.26
No. Observations:                 125   AIC:                             1017.
Df Residuals:                     115   BIC:                             1045.
Df Model:                           9                                         
Covariance Type:            nonrobust                                         
==============================================================================
  （后略）

从打印的R²和RMSE值可以看出，两个模型都具有较好的拟合效果。

响应面图分析：响应面图直观地展示了去除率随两个因素变化的趋势。从图中可以观察到，ROX去除率对pH和吸附剂用量的变化较为敏感，呈现明显的非线性关系。As(V)去除率则对温度和pH的变化更为敏感。这些信息有助于我们理解各因素对去除率的影响程度和方式。

热力图分析：热力图更清晰地展示了去除率在不同因素组合下的分布。对于ROX，我们可以看到在高pH和高吸附剂用量时，去除率达到最高。而对于As(V)，中等pH和较高温度时去除率较高。这些信息可以指导我们选择最佳的操作条件。

偏依赖图分析：（略）

在具体的系数分析中，我们发现几个变量对ROX去除率有显著影响….略

现在，让我们转向As(V)去除率模型。在As(V)模型的具体系数分析中，我们发现….略

问题二模型的建立与求解

思路分析

问题二要求我们确定反应温度、溶液pH和吸附剂用量的最佳组合，以使As(V)和ROX的总吸附量达到最大。这是一个典型的多目标优化问题，我们需要同时考虑两种不同污染物的吸附量，并在三个变量的空间中寻找最优解。考虑到问题的复杂性和非线性特征，我们可以采用一种综合的方法来解决这个问题。首先，我们需要建立As(V)和ROX吸附量与三个变量之间的关系模型。基于问题一的分析结果，我们知道这些关系可能是非线性的，并且存在交互效应。因此，我们可以考虑使用响应面法（Response Surface Methodology，RSM）来构建这些模型。RSM不仅可以捕捉变量之间的非线性关系，还能反映它们之间的交互作用。

考虑到优化问题的复杂性和可能存在的多个局部最优解，我们需要一个强大的全局优化算法。在这种情况下，我们可以选择使用粒子群优化算法（Particle Swarm Optimization，PSO）。PSO是一种受生物群体行为启发的元启发式算法，它能够有效地在高维空间中搜索全局最优解，并且对于非凸优化问题表现良好。

响应面多目标粒子群优化模型建立

基于上述思路，我们提出一个"响应面多目标粒子群优化模型"（Response Surface Multi-objective Particle Swarm Optimization Model，RSMPSOM）来解决这个问题。这个模型结合了响应面法的建模能力、多目标优化的灵活性和粒子群算法的全局搜索能力。

首先，我们使用响应面法为As(V)和ROX的吸附量建立二阶多项式模型。对于每种污染物，模型可以表示为：

$$Q_i={\beta }_0+\sum _{j=1}^3{\beta }_j{x}_j+\sum _{j=1}^3{\beta }_{jj}{x}_j^2+\sum _{j<k}^3{\beta }_{jk}{x}_j{x}_k+\epsilon$$

其中，$Q_i$ 表示污染物 $i$ (As(V)或ROX)的吸附量，$x_j$ 代表第 $j$ 个自变量（温度、pH或吸附剂用量），$\beta$ 是回归系数，$\epsilon$ 是随机误差项。

然后，我们构建一个加权和目标函数：

$$F=w_1{Q}_{As(V)}+w_2{Q}_{ROX}$$

其中，$w_1$ 和 $w_2$ 是权重系数，满足 $w_1+w_2=1$。

最后，我们使用PSO算法来最大化这个目标函数。PSO算法中的每个粒子代表一个可能的解（即温度、pH和吸附剂用量的组合），通过迭代更新粒子的位置和速度来搜索最优解。

自适应多目标粒子群优化算法步骤

为了提高算法的性能和适应性，我们提出一种"自适应多目标粒子群优化算法"（Adaptive Multi-objective Particle Swarm Optimization Algorithm，AMPSOA）。该算法的主要步骤如下：（略）

问题二模型的求解

基于上述模型和算法，我们使用Python实现求解过程。（部分，完整代码见附件）

# 自定义粒子群优化算法
def pso(objective, lb, ub, n_particles=30, max_iter=100):
    dim = len(lb)
    # 初始化粒子群
    particles = np.random.rand(n_particles, dim) * (np.array(ub) - np.array(lb)) + np.array(lb)
    velocities = np.zeros((n_particles, dim))
    
    # 初始化个体最优和全局最优
    pbest = particles.copy()
    pbest_obj = np.array([objective(p) for p in pbest])
    gbest = pbest[pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()
    
    # PSO参数
    w = 0.7  # 惯性权重
    c1 = 1.5  # 个体学习因子
    c2 = 1.5  # 社会学习因子
    
    for _ in range(max_iter):
        # 更新速度和位置
        r1, r2 = np.random.rand(2)
        velocities = (w * velocities 
                      + c1 * r1 * (pbest - particles) 
                      + c2 * r2 * (gbest - particles))
        particles += velocities
        
        # 边界处理
        particles = np.clip(particles, lb, ub)
        
        # 更新个体最优和全局最优
        obj_values = np.array([objective(p) for p in particles])
        improve = obj_values < pbest_obj
        pbest[improve] = particles[improve]
        pbest_obj[improve] = obj_values[improve]
        
        if obj_values.min() < gbest_obj:
            gbest = particles[obj_values.argmin()]
            gbest_obj = obj_values.min()
    
    return gbest, gbest_obj

# 设置PSO参数
lb = [0.2, 15, 3]  # 下界
ub = [1.0, 45, 9]  # 上界

# 运行PSO优化
xopt, fopt = pso(objective_function, lb, ub)

# 打印最优解
print(f"最优解: 吸附剂用量 = {xopt[0]:.4f} g/L, 温度 = {xopt[1]:.4f} °C, pH = {xopt[2]:.4f}")
print(f"最大总吸附量: {-fopt:.4f} mg/g")
后略

问题二模型求解结果可视化与分析

求解得到的最优解如下：

在获得最优解后，我们可以通过多种方式可视化结果并进行深入分析：

根据问题2的求解结果，我们可以对ROX和As(V)的吸附模型以及最优解进行深入的分析和解释：（略）

首先，让我们分析ROX模型的系数：（略）

吸附剂用量（第一个系数-8.10129411）对ROX吸附量有很强的负面影响。这意味着随着吸附剂用量的增加，ROX的吸附量反而会显著降低。这可能是由于高浓度的吸附剂导致了吸附位点之间的相互干扰或阻碍。

问题三模型的建立与求解

思路分析

问题三要求我们设计额外的5次实验，以进一步优化和验证As(V)和ROX的吸附过程。这是一个典型的实验设计问题，涉及到如何在有限的实验次数内最大化信息获取。考虑到我们已经有了初步的模型和优化结果，新的实验设计应该既能够验证现有模型的预测，又能够探索模型可能忽略的区域或关系。此外，由于As(V)和ROX的吸附行为存在差异，我们需要平衡这两种污染物的研究需求。基于这些考虑，我们可以采用一种结合了多准则决策、空间填充和模型不确定性分析的高级实验设计方法。

我们的设计思路应该包含以下几个关键点：首先，我们需要在现有实验点的基础上，选择能够最大程度填充实验空间的新点，以确保我们对整个参数空间有全面的了解。其次，我们应该关注模型预测的高不确定性区域，因为这些区域可能包含重要的信息，有助于改进模型的准确性。再次，我们应该在模型预测的最优区域附近进行一些实验，以验证和细化最优条件。最后，考虑到As(V)和ROX的行为差异，我们需要在实验设计中平衡这两种污染物的需求。

基于以上思考，我们提出一种"多目标自适应实验设计优化模型"（Multi-Objective Adaptive Experimental Design Optimization Model，MOAEDOM）。这个模型将综合考虑空间填充、模型不确定性、预测优化和多目标平衡等多个方面，通过一个迭代的过程来选择最优的新实验点。