基于PPLSR的产品设计方法

基于PPLSR模型的工业过程产品设计

摘要

基于偏最小二乘回归（PLSR）模型的产品设计方法已被广泛用于获得期望产品质量。然而，由于PLSR模型是一种确定性模型，该方法无法应对实际工业过程中存在的显著不确定性。本研究提出了一种基于概率偏最小二乘回归（PPLSR）模型的产品设计方法学。从概率视角出发，利用已有数据建立操作条件与产品质量之间的回归关系。Then计算能够产生期望产品质量的操作条件的概率分布，并对操作条件的不确定性进行定量分析。此外，在给定操作条件下预测产品质量，进而计算并分析产品质量误差的数值及其成因。所提出的方法学在青霉素发酵生产过程中进行了测试。

关键词 —产品质量，产品设计，概率偏最小二乘回归，条件分布，模型反演

引言

在不同操作条件下运行的过程可以生产出不同的产品。例如，在青霉素发酵生产过程中，通过使用不同的底物进料速率、通气速率和搅拌功率，可以生产出不同青霉素浓度的产品。为了生产具有期望产品质量的所需产品，必须确定适当的操作条件，这正是产品设计的任务。随着对提高产品质量和成本比率的要求不断提高，产品设计在工业生产中变得越来越重要。[1-2]。早期的产品设计方法基于工业生产过程的精确机理模型以及期望产品质量来计算操作条件[3-4]。然而，在实际工业过程中建立准确的机理模型是困难的。

数据驱动方法不需要工业过程的精确模型，而是利用历史数据建立模型。将历史操作条件作为输入数据，相应的产品质量作为输出数据，通过回归分析建立输入与输出之间的关系。然后通过模型反演计算出适当的操作条件。典型示例包括主成分回归（PCR）方法和偏最小二乘回归（PLSR）方法。这为确定操作条件提供了新思路。1998年，Jaeckle和MacGregor首次将PLSR模型引入产品设计。该方法利用工业过程的输入数据（操作条件）和输出数据（产品质量）建立PLSR模型，并获得操作条件与产品质量之间的关系，然后通过模型反演方法设计操作条件，并取得了良好效果[7-8]。此后，基于PLSR模型的产品设计受到广泛关注，针对不同工艺提出了多种扩展方法，并应用于药物制剂和小麦加工的开发中[9-10]。

然而，由于实际工业过程中的随机噪声，输入和输出实际上是随机变量，而经典数据驱动方法建立的模型是确定性模型，无法反映变量的真实特性[11]。因此，概率模型更具优势，其可采用与经典确定性模型类似的方式建立，并将输入和输出数据视为概率性的。例如，Tipping基于主成分分析(PCA)提出了概率PCA[12]，并将其应用于工业过程监控[13],，取得了良好的应用效果；Li将概率分布引入PLSR模型[14]：假设所有变量、潜在变量和残差均满足正态分布，建立了概率偏最小二乘回归(PPLSR)模型。

由于PPLSR相比PLSR具有优越性，为了更好地描述操作条件的不确定性，本文将PPLSR模型引入到工业过程产品设计中。考虑到过程数据受随机噪声影响所导致的不确定性，采用PPLSR模型捕捉操作条件与产品质量之间的回归关系，并获取满足期望产品质量的操作条件的概率分布。该概率分布用于量化在期望产品质量下操作条件的不确定性。进一步地，基于此概率分布，在相应操作条件下计算产品质量，定量计算产品的误差，并分析误差产生的原因。以青霉素发酵过程的工业案例进行验证，结果表明了所提出方法论的可行性与优越性。

II. PLSR模型与产品设计方法

A. PLSR模型

偏最小二乘回归是一种常用的多元回归建模方法。通过提取潜在变量，将高维数据空间投影到低维数据空间，从而考虑输入和输出变量之间的相关性，有效处理数据共线性问题。然后建立输入和输出变量之间的回归关系[15]。给定一个过程操作条件矩阵 $ \mathbf{X} \in \mathbb{R}^{I \times D_x} $ 和相应的产品质量矩阵 $ \mathbf{Y} \in \mathbb{R}^{I \times D_y} $，其中 $ I $ 是不同操作条件下的生产批次数量。$ D_x $ 是测量的操作条件的数量。$ D_y $ 是不同操作条件下产品质量的数量。然后，PLSR模型可以表示为：

$$
\mathbf{X} = \mathbf{T}\mathbf{P}^T + \mathbf{E} \tag{1}
$$
$$
\mathbf{Y} = \mathbf{T}\mathbf{Q}^T + \mathbf{F} \tag{2}
$$

PLSR模型将样本数据空间划分为得分空间和残差子空间，其中 $ \mathbf{E} \in \mathbb{R}^{I \times D_x} $ 和 $ \mathbf{F} \in \mathbb{R}^{I \times D_y} $ 分别是X和Y的残差矩阵，$ \mathbf{P} \in \mathbb{R}^{D_x \times A} $ 和 $ \mathbf{Q} \in \mathbb{R}^{D_y \times A} $ 是载荷矩阵。$ \mathbf{T} \in \mathbb{R}^{I \times A} $ 是潜在变量，用于充分解释X和Y的数据信息，并满足 $ \mathbf{T} = \mathbf{XW} $，其中 $ \mathbf{W} \in \mathbb{R}^{D_x \times A} $ 是X的投影矩阵，是 $ \mathbf{X}^T\mathbf{Y}\mathbf{Y}^T\mathbf{X} $ 的特征向量，并满足 $ \mathbf{W}^T\mathbf{X}^T\mathbf{Y}\mathbf{Y}^T\mathbf{X}\mathbf{W} = \lambda $。A是潜在变量的维度，由交叉验证确定。

B. 基于PLSR模型的产品设计方法

偏最小二乘回归（PLSR）方法通过提取输入变量和输出变量的关键潜在变量信息，最大化其相关性，建立输入‐输出关系。基于PLSR模型的产品设计方法利用过程的输入和输出数据建立PLSR模型，并获得操作条件与产品质量之间的回归关系。通过模型反演推导出操作条件。模型反演的目标是估计可行的操作条件 $ \mathbf{x} {\text{new}} $，以实现期望产品质量 $ \mathbf{y} {\text{des}} $。

由于T通过载荷矩阵X和Y关联矩阵P和Q，并将高维数据空间投影到低维数据空间。为了解决共线性问题，可将潜在变量矩阵T缩放为 $ \mathbf{T} = \mathbf{US} $，其中U为标准正交矩阵，$ \mathbf{S} = (\mathbf{T}^T\mathbf{T})^{1/2} $ 是一个 $ (A \times A) $ 对角矩阵。则相应的PLS模型变为：

$$
\mathbf{X} = \mathbf{USP}^T \tag{3}
$$
$$
\mathbf{Y} = \mathbf{USQ}^T \tag{4}
$$

给定期望产品质量 $ \mathbf{y}_{\text{des}} $，它满足：

$$
\mathbf{y} {\text{des}} = \mathbf{u} {\text{new}}^T \mathbf{SQ}^T \tag{5}
$$

然后，通过模型反演得到的新操作条件将取决于潜在变量的维度和产品质量。当产品质量的维度 $ D_y $ 大于或等于过程操作空间 $ \mathbf{X} $ 的维度时，即 $ D_y \geq A $。在这种情况下，新操作条件 $ \mathbf{u}_{\text{new}} $ 的坐标由下式给出：

$$
\mathbf{u} {\text{new}} = (\mathbf{Q}^T\mathbf{Q})^{-1}\mathbf{Q}^T\mathbf{y} {\text{des}} \mathbf{S}^{-1} \tag{6}
$$

对应的 $ \mathbf{x}_{\text{new}} $ 值为：

$$
\mathbf{x} {\text{new}}^T = \mathbf{u} {\text{new}}^T \mathbf{SP}^T \tag{7}
$$
$$
\mathbf{x} {\text{new}}^T = \mathbf{y} {\text{des}}^T (\mathbf{Q}^T\mathbf{Q})^{-1}\mathbf{Q}^T\mathbf{P}^T \tag{8}
$$

此时， $ \mathbf{x}_{\text{new}} $ 的解是唯一的。

当 $ D_y < A $ 时，可通过求解伪逆矩阵得到 $ \mathbf{u}_{\text{new}} $ 的值：

$$
\mathbf{u} {\text{new}} = \mathbf{y} {\text{des}}^T \mathbf{S}^{-1} \mathbf{Q}^T (\mathbf{QS}^T\mathbf{Q})^{-1} \tag{9}
$$

则 $ \mathbf{x}_{\text{new}} $ 的解为：

$$
\mathbf{x} {\text{new}}^T = \mathbf{y} {\text{des}}^T \mathbf{S}^{-1} \mathbf{Q}^T (\mathbf{QS}^T\mathbf{Q})^{-1} \mathbf{P}^T \tag{10}
$$

然而，如果 $ D_y < A $，那么 $ \mathbf{x}_{\text{new}} $ 仅投影到一个 $ D_y $ 维空间，并且 $ (A - D_y) $ 维零空间的存在将弥补实际操作空间。

操作条件可以从方程(8)和方程(10)中确定。需要特别说明的是，当主维度小于输入维度时，操作条件的确定由 $ D_y $ 维空间的操作条件和 $ (A - D_y) $ 维的零空间共同组成。因此，在设计产品质量时，操作条件仍具有一定的自由度，需根据实际需求进行调整。

III. PPLSR模型与新产品设计方法

值得注意的是，基于PLSR模型的产品设计方法利用数据建模提取操作条件与产品质量之间的相关性。当给定期望产品质量时，通过模型反演重新计算操作条件，从而实现新产品设计。然而，PLSR模型作为确定性模型，无法定量反映由随机噪声等因素影响的工业过程中的不确定性。并且由于操作条件通常具有一定的自由度，对操作条件不确定性的描述显得尤为重要。但是，基于PLSR模型的产品设计在模型反演中未能定量描述操作条件的不确定性，也未将操作条件的不确定性与产品质量相关联。因此，本文提出了一种基于PPLSR模型的产品设计方法。

A. PPLSR模型

PPLSR模型是PLSR模型的概率扩展，其模型结构与PLSR模型相同。假设工业过程中的输入和输出数据为：

$$
\mathcal{X} = { \mathbf{x}_n | n = 1, \dots, I }, \quad \mathbf{x}_n \in \mathbb{R}^{D_x} \tag{11}
$$
$$
\mathcal{Y} = { \mathbf{y}_n | n = 1, \dots, I }, \quad \mathbf{y}_n \in \mathbb{R}^{D_y} \tag{12}
$$

然后，PPLSR模型可以表示为：

$$
\mathbf{x}_n = \mathbf{P}\mathbf{u}_n + \boldsymbol{\mu}_x + \boldsymbol{\xi}_n \tag{13}
$$
$$
\mathbf{y}_n = \mathbf{C}\mathbf{u}_n + \boldsymbol{\mu}_y + \boldsymbol{\epsilon}_n \tag{14}
$$

其中 $ \mathbf{u}_n $ 是主成分矩阵，且满足 $ \mathbf{u}_n \sim N(0, \mathbf{I}) $， $ \mathbf{P} \in \mathbb{R}^{D_x \times K} $， $ \mathbf{C} \in \mathbb{R}^{D_y \times K} $，K是 $ \mathbf{u}_n $ 的主成分维度。 $ \boldsymbol{\mu}_x $ 是X的中心， $ \boldsymbol{\mu}_y $ 是Y的中心， $ \boldsymbol{\xi}_n $ 和 $ \boldsymbol{\epsilon}_n $ 是各向同性的零均值高斯噪声，即 $ \boldsymbol{\Sigma}_x = \sigma_x^2 \mathbf{I} $， $ \boldsymbol{\Sigma}_y = \sigma_y^2 \mathbf{I} $。

PPLSR模型首先假设主成分的先验分布，然后根据观测样本的最大似然函数推断主成分的后验分布，并基于贝叶斯公式化利用期望最大化（EM）算法进行参数估计。随后得到输入和输出变量之间的回归关系矩阵为：

$$
\mathbf{B}_{\text{PPLSR}} = \mathbf{C}(\mathbf{P}^T\mathbf{P} + \boldsymbol{\Sigma}_x)^{-1}\mathbf{P}^T \tag{15}
$$

关于改进的偏最小二乘回归（PPLSR）的详细步骤可参见文献[14]。

B. 产品设计方法

与标准PLSR模型相比，PPLSR模型基于贝叶斯法则，从概率视角解释了模型中输入和输出变量之间的关系。

当获得期望产品质量（ $ \mathbf{y} {\text{des}} $ ）时，可通过条件分布得到新的操作条件（ $ \mathbf{x} {\text{new}} $ ）：

$$
p(\mathbf{x} {\text{new}} | \mathbf{y} {\text{des}}) = \int p(\mathbf{x} {\text{new}} | \mathbf{u}_n) p(\mathbf{u}_n | \mathbf{y} {\text{des}}) d\mathbf{u}_n \tag{16}
$$

由公式(13)和公式(14)可得，新操作条件( $ \mathbf{x} {\text{new}} $ )关于 $ \mathbf{u}_n $ 的条件分布以及期望产品质量( $ \mathbf{y} {\text{des}} $ )关于 $ \mathbf{u}_n $ 的条件分布为：

$$
p(\mathbf{x} {\text{new}} | \mathbf{u}_n) = N(\mathbf{P}\mathbf{u}_n + \boldsymbol{\mu}_x, \boldsymbol{\Sigma}_x) \tag{17}
$$
$$
p(\mathbf{y} {\text{des}} | \mathbf{u}_n) = N(\mathbf{C}\mathbf{u}_n + \boldsymbol{\mu}_y, \boldsymbol{\Sigma}_y) \tag{18}
$$

$ \mathbf{u}_n $ 的先验分布为：

$$
p(\mathbf{u}_n) = N(0, \mathbf{I}) \tag{19}
$$

根据贝叶斯法则 $ p(\mathbf{u} n | \mathbf{y} {\text{des}}) p(\mathbf{y} {\text{des}}) = p(\mathbf{y} {\text{des}} | \mathbf{u}_n) p(\mathbf{u}_n) $ 和公式(19)，我们可以得到新操作条件的条件分布为：

$$
p(\mathbf{u} n | \mathbf{y} {\text{des}}) = N(\mathbf{M}^{-1}\mathbf{C}^T(\mathbf{y}_{\text{des}} - \boldsymbol{\mu}_y), \boldsymbol{\Sigma}_y \mathbf{M}^{-1}) \tag{20}
$$

其中 $ \mathbf{M} = \mathbf{C}^T\mathbf{C} + \boldsymbol{\Sigma}_y $。

由公式(16)、公式(17)和公式(20)可得，在期望产品质量下新操作条件的概率分布为：

$$
p(\mathbf{x} {\text{new}} | \mathbf{y} {\text{des}}) = \int p(\mathbf{x} {\text{new}} | \mathbf{u}_n) p(\mathbf{u}_n | \mathbf{y} {\text{des}}) d\mathbf{u}_n = N(\mathbf{c}, \mathbf{V}) \tag{21}
$$

其中 $ \mathbf{c} = \mathbf{P}\mathbf{M}^{-1}\mathbf{C}^T(\mathbf{y}_{\text{des}} - \boldsymbol{\mu}_y) + \boldsymbol{\mu}_x $, $ \mathbf{V} = \mathbf{P}\mathbf{M}^{-1}\mathbf{P}^T + \boldsymbol{\Sigma}_x $。

当期望产品质量为 $ \mathbf{y} {\text{des}} $ 时，操作条件 $ \mathbf{x} {\text{new}} $ 的期望值为：

$$
\mathbf{x} {\text{new}} = \mathbf{c} = \mathbf{P}\mathbf{M}^{-1}\mathbf{C}^T(\mathbf{y} {\text{des}} - \boldsymbol{\mu}_y) + \boldsymbol{\mu}_x \tag{22}
$$

C. 产品设计方法误差计算与分析

为了确定新设计的操作条件是否满足期望产品质量，需要进行误差计算和分析。首先，根据产品设计值计算操作条件，然后将操作条件添加到生产过程中以获得产品质量。产品质量与期望值之间的差异定义为设计误差。

当计算出操作条件的概率分布后，重新计算操作条件分布下新产品质量分布为：

$$
p(\mathbf{y} {\text{new}} | \mathbf{x} {\text{new}}) = \int p(\mathbf{y} {\text{new}} | \mathbf{u}_n) p(\mathbf{u}_n | \mathbf{x} {\text{new}}) d\mathbf{u}_n \tag{23}
$$

From Eq.(17), we can get the distribution by Bayesian as:

$$
p(\mathbf{u} n | \mathbf{x} {\text{new}}) = N(\mathbf{M} x^{-1}\mathbf{P}^T(\mathbf{x} {\text{new}} - \boldsymbol{\mu}_x), \boldsymbol{\Sigma}_x \mathbf{M}_x^{-1}) \tag{24}
$$

其中 $ \mathbf{M}_x = \mathbf{P}^T\mathbf{P} + \boldsymbol{\Sigma}_x $。

由公式(18)、公式(23)和公式(24)可得，在新的操作条件下，新产品质量的条件分布为：

$$
p(\mathbf{y} {\text{new}} | \mathbf{x} {\text{new}}) = N(\mathbf{m}, \mathbf{S}) \tag{25}
$$

根据新的操作条件，获得新产品质量期望如下：

$$
\mathbf{y} {\text{new}} = \mathbf{m} = \mathbf{B} {\text{PPLSR}}(\mathbf{x} {\text{new}} - \boldsymbol{\mu}_x) + \boldsymbol{\mu}_y = \mathbf{C}(\mathbf{P}^T\mathbf{P} + \boldsymbol{\Sigma}_x)^{-1}\mathbf{P}^T(\mathbf{x} {\text{new}} - \boldsymbol{\mu}_x) + \boldsymbol{\mu}_y \tag{26}
$$
$$
\mathbf{S} = \sigma^2 \mathbf{C}(\mathbf{P}^T\mathbf{P} + \boldsymbol{\Sigma}_x)^{-1}\mathbf{C}^T + \boldsymbol{\Sigma}_y \tag{27}
$$

协方差矩阵 $ \mathbf{S} $ 表示在新的操作条件下获得的新产品质量的不确定性。

利用得到的新产品质量与预期的产品质量计算误差，该误差如下所示：

$$
\mathbf{y} {\text{new}} - \mathbf{y} {\text{des}} = (\mathbf{C}(\mathbf{P}^T\mathbf{P} + \boldsymbol{\Sigma} x)^{-1}\mathbf{P}^T + \mathbf{C}^T\boldsymbol{\Sigma}_y^{-1}\mathbf{C} - \mathbf{I})(\mathbf{y} {\text{des}} - \boldsymbol{\mu}_y) \tag{28}
$$

From Eq.(28), when $ \boldsymbol{\Sigma}_x = \boldsymbol{\Sigma}_y = 0 $, the error is zero. But $ \boldsymbol{\Sigma}_x $ 并且X和Y的残差 $ \boldsymbol{\Sigma}_y $ 通常不为零，因此误差是存在的，且由设计操作条件引起的误差与X和Y的残差相关，PPLSR模型残差表示实际观测值与回归估计值之间的差异。

IV. 实验

在本节中，以青霉素发酵生产为实验对象，验证所提出方法的优越性。

A. 过程与数据

青霉素发酵过程是一种复杂的生化反应过程[18-19]，需要在适当的操作条件下促进青霉素的合成。该发酵过程主要分为两个阶段：第一阶段是青霉素发酵过程中的微生物培养阶段，主要产生青霉素发酵所需的菌体；第二阶段是通过进料促进青霉素的合成。

本文中的数据来自西纳尔开发的青霉素发酵过程模拟软件Pensim2.0[20]。该软件通过设置不同的操作条件和控制参数，对青霉素发酵过程进行仿真，并生成不同的青霉素发酵过程数据。青霉素发酵过程的具体原理和介绍参见[18]。

本章的基本思路是在Pensim2.0模拟软件中，模拟正常条件下的青霉素发酵过程，操作条件和输入控制参数设置为：

变量	默认值	值	集合范围
底物浓度(g/L)	15	1-35
溶解氧浓度(/Lmmole)	1.16	1.0‐1.2
生物量浓度(g/L)	0.1	0‐3.0
penicillin浓度(g/L)	0	0
培养液体积(L)	100.0	50.0‐1000.0
PH	5.0	4.5‐6.5
温度(开尔文)	298	297‐300

控制参数	默认值	值	设定范围
通气速率(升/小时)	8.6	8.0‐9.0
搅拌器功率(瓦)	30.0	29.0‐31.0
底物进料速率(升/小时)	0.042	0.03‐0.045
底物 feed 温度(开尔文)	296	293‐300
PH	5.0	4.5‐5.3
温度(开尔文)	298	297‐298

提供40组青霉素发酵过程数据，然后基于这些数据建立了青霉素发酵生产的数学模型。由此，针对给定的目标青霉素浓度反推出设计的操作条件。最后，将设计的操作条件代入Pensim2.0中，得到青霉素产物浓度，并与期望浓度进行比较，计算新方法的误差产品设计方法被计算以评估所提出的产品设计方法。

B. 实验结果

本实验获得的建模数据设置了不同的操作条件，包括通气速率、搅拌功率、底物进料速率、生物量浓度和底物浓度。共生成了40批次的过程数据，$\mathbf{X}$表示40批次的操作条件参数值，$\mathbf{Y}$表示每批次结束时的青霉素浓度值，其中30批次用作训练数据集，最后10批次用作测试数据集，以检验所提出方法的预测精度。预测精度如图1所示：

从图1可以看出，与PLSR模型相比，PPLSR模型对十个批次的预测值更接近真实值，表明其具有更高的预测精度。同时，表III对均方根误差 (RMSE) 进行了比较，也可以看出PPLSR模型的预测精度高于PLSR模型。PPLSR模型对输入‐输出模型回归关系的解释更为准确。

)

模型	RMSE
PLSR	0.0599
改进的偏最小二乘回归	0.0090

均方根误差（RMSE）的计算公式为：

$$
\text{RMSE} = \sqrt{\frac{1}{I} \sum_{i=0}^{I} (y_0 - \hat{y})^2} \tag{29}
$$

其中，$i = 1, 2, \dots, I$，$y_0$为真实值，$\hat{y}$为预测值。$I$为批次数量。RMSE越小，预测精度越高。

基于PPLSR模型更高的预测精度，设计实验以验证新产品设计方法的有效性和优越性。

首先，给出三组期望青霉素浓度：1.2 g/L、1.25 g/L、1.3 g/L，利用本文方法进行产品设计，在期望最终产品质量 $\mathbf{y} {\text{des}}$ 下设计新的操作条件 $\mathbf{x} {\text{new}}$，使用公式 (27) 和公式 (29) 计算新产品质量并计算误差，结果如下所示：

期望的最终产品质量 (g/L) $\mathbf{y}_{\text{des}}$	新产品 质量 (g/L) $\mathbf{y}_{\text{new}}$	误差 (g/L)
1.2	1.2206	0.0206
1.25	1.2495	0.0005
1.3	1.2784	0.0216

然后将设计出的新操作条件的期望值代入 Pensim2.0中，通过模拟生产得到三组新的青霉素浓度 $(g/L)$。参数值和均方根误差如表 V. 所示。

期望的最终产品质量 (g/L) $\mathbf{y}_{\text{des}}$	质量产 品 (g/L) $\mathbf{y}$	RMSE
1.2	1.203	0.0044
1.25	1.253	0.0044
1.3	1.299	0.0044

工业产品

可以看出，新的青霉素浓度接近预期的青霉素浓度，且得到的均方根误差较小，表明本文提出的方法能够实现良好的产品设计效果。

五、结论

本文将PPLSR模型应用于工业生产过程中的产品设计，提出了基于PPLSR模型的产品设计方法学。通过利用工业过程中不同操作条件和产品质量数据建立PPLSR模型，然后从概率角度进行模型反演。在期望产品质量下设计新的操作条件，并获得新操作条件的概率分布。新操作条件的协方差矩阵定量描述了操作条件的不确定性，并通过条件概率将产品质量与操作条件的不确定性相关联。计算了在新操作条件下的产品质量误差，并分析了产品质量误差的成因。最后，结合青霉素发酵的仿真，所提出的方法获得了更准确的输入‐输出回归关系。在产品设计过程中，通过设计实验以实现期望的产品设计结果，获得了期望产品质量，体现了基于PPLSR模型的产品设计方法学的优越性。