自动化学报  2017, Vol. 43 Issue (2): 271-279   PDF    
基于PLS交叉积矩阵非相似度分析的MPC性能监控与诊断
尚林源, 田学民, 曹玉苹, 蔡连芳     
中国石油大学(华东) 信息与控制工程学院 青岛 266580
摘要: 针对传统基于输出协方差矩阵的性能监控方法未充分考虑过程变量与输出变量之间的相关性问题,提出一种基于偏最小二乘(Partial least squares,PLS)交叉积矩阵非相似度分析的性能监控与诊断方法,用于多变量模型预测控制(Model predictive control,MPC)系统.首先,考虑模型预测控制系统的控制结构,构造包含预测误差的增广过程变量与输出变量相关性的PLS交叉积矩阵,通过非相似度分析方法将交叉积矩阵的非相似度比较转化为转换矩阵特征值的比较.然后提取转换矩阵中表征最大非相似度的l个特征值构造实时性能指标,对MPC系统进行性能监控.检测到性能下降后,进一步利用转换矩阵的特征值诊断性能恶化源.Wood-Berry二元精馏塔上的仿真结果表明,所提方法能够有效地提高监控性能,并准确地定位性能恶化源.
关键词: 模型预测控制     性能监控与诊断     偏最小二乘     交叉积矩阵     非相似度分析    
MPC Performance Monitoring and Diagnosis Based on Dissimilarity Analysis of PLS Cross-product Matrix
SHANG Lin-Yuan, TIAN Xue-Min, CAO Yu-Ping, CAI Lian-Fang     
College of Information and Control Engineering, China University of Petroleum (East China), Qingdao 266580
Received: 2015-11-20, Accepted: 2016-06-06.
Foundation Item: Supported by National Natural Science Foundation of China (61273160, 61403418), the Fundamental Research Funds for the Central Universities (15CX06063A), and Natural Science Foundation of Shandong Province (ZR2014FL016, ZR2016FQ21)
Author brief: SHANG Lin-Yuan Ph. D. candidate at the College of Information and Control Engineering, China University of Petroleum. His main research interest is control system performance assessment;
CAO Yu-Ping Ph. D., lecturer at the College of Information and Control Engineering, China University of Petroleum. Her research interest covers fault diagnosis and prediction;
CAI Lian-Fang Ph. D. at the College of Information and Control Engineering, China University of Petroleum. His research interest covers fault detection and diagnosis in industrial processes
Corresponding author. TIAN Xue-Min Professor at the College of Information and Control Engineering, China University of Petroleum. His research interest covers modeling, advanced process control and optimization technology, intelligent monitoring, and fault diagnosis. Corresponding author of this paper
Recommended by Associate Editor ZHONG Mai-Ying
Abstract: Performance monitoring methods for control systems based on output covariance matrix can not sufficiently exploit the correlation between the process variables and output variables. To solve this problem, a performance monitoring and diagnosis method based on dissimilarity analysis of partial least squares (PLS) cross-product matrix is proposed for multivariate model predictive control (MPC) systems. Firstly, the PLS cross-product matrix, which contains the correlation information of augmented process variables and output variables, is constructed. And dissimilarity analysis is carried out to transform dissimilarity comparison of cross-product matrixes to eigenvalue comparison of transformed matrixes. Then, using the l eigenvalues, which include the maximum dissimilarity information, a new performance index is constructed to monitor the performance of MPC system. Finally, the index is further improved to meet the requirement of diagnosing the root cause of performance deterioration. Simulation results on the Wood-Berry binary distillation column demonstrate that the proposed method can effectively enhance the monitoring performance and accurately locate the source of performance deterioration.
Key words: Model predictive control (MPC)     performance monitoring and diagnosis     partial least squares (PLS)     cross-product matrix     dissimilarity analysis    

模型预测控制(Model predictive control, MPC) 作为最具代表性的先进控制策略, 已经在复杂工业过程控制中取得了广泛应用[1-3].然而, 在实际运行过程中, 在各种因素的影响下(如干扰特性变化、模型失配、约束饱和等), MPC系统的性能会逐渐下降, 从而影响产品质量甚至生产过程的安全[4].因此, 对MPC系统进行性能维护受到学术界和工程界的高度重视, 而性能的维护需要高效、可靠的控制性能监控技术.

Harris[5]于1989年提出的单变量最小方差控制(Minimum variance control, MVC) 性能基准, 开创了对控制性能评价技术研究的先例.此后, 许多更为实用的性能评价基准被相继提出, 如线性二次型高斯基准(Linear quadratic Gaussian, LQG)[6]、用户自定义基准[7]以及基于历史指标的性能基准[4]等.在现代过程工业中, 分布式系统(Distributed control systems, DCS) 和SCADA (Supervisory control and data acquisiton) 技术得到广泛应用, 实时采集和保存了大量过程数据, 从而促进了多变量统计过程控制(Multivariate statistical process control, MSPC) 方法在控制系统性能监控和诊断领域的应用和研究[8-20]. MSPC方法具有简单有效、非侵入等特点, 从而受到研究人员的关注[10]. Zhang等[11]使用基准数据建立PCA (Principal component analysis) 模型, 并根据改进的统计量评价MPC控制性能. AlGhazzawi等[12]使用静态和递归PCA/PLS (Partial least squares) 辨识MPC控制回路的异常工况并分离可能的恶化源. Chen等[13]将PCA与ARMA (Auto-regressive and moving average model) 模型结合解决变量的互相关和自相关问题, 用于评价系统的可达最小方差性能. Shang等[14]使用慢特征分析方法进行厂级性能监控,并用贡献图方法确定导致性能变化的过程变量. Das等[15]基于赫斯特指数和马氏距离进行多变量控制回路的性能评估. Yu等[16-17]提出了基于统计协方差行列式的性能评价方法, 该方法考虑了协方差的非对角线元素, 使用行列式表征协方差所张成的超椭球体体积, 更加符合实际应用.通过分析MPC模型预测误差, 田学民等[18-19]扩展了监控变量集, 提出了基于2-范数的协方差性能基准和基于协方差预测残差的性能监控方法.考虑到输出协方差超椭球体中的方向信息, Li等[20]提出了基于非相似度分析的协方差性能监控方法.上述基于协方差的方法在进行性能监控时只考虑了所有测量变量或输出变量的协方差, 没有考虑性能变化对过程变量与输出变量间相关性的影响.

本文以充分挖掘系统过程变量与输出变量数据间相关关系及数据分布特征为出发点, 通过分析PLS与典型相关分析(Canonical correlation analysis, CCA) 之间的联系构造过程变量与输出变量数据的交叉积特征矩阵, 并以基准性能和监控性能下特征矩阵之间的整体非相似度为准则提取两个特征矩阵之间的最大非相似度信息, 用于建立可量化的MPC系统性能指标; 进一步利用非相似度分析后得到的特征值构造性能诊断指标, 定位性能恶化源.最后, 在Wood-Berry二元精馏塔上对所提方法的有效性进行验证.

1 PLS交叉积矩阵

假设两个数据集分别为${ { {X}}}\in{\rm{{\bf R}}}^{N \times p}$, ${ { {Y}}}\in{\rm{{\bf R}}}^{N \times m}$, 其中Npm分别为样本长度和数据集的变量个数.使用PLS算法建立两个数据集之间的关系模型

$\begin{array}{l} { X}= { {TP}}^{\rm T} + { F} \\ { Y}= { {UQ}}^{\rm T} + { G} \\ \end{array}$ (1)

其中, TU$N \times l$维得分向量矩阵, PQ分别为$p \times l$$m \times l$维载荷矩阵, FG为残差矩阵.经典的PLS算法是非线性迭代偏最小二乘(Nonlinear iterative partial least squares, NIPALS) 算法[21], 该算法中X的权重向量$ {\pmb {w}}$Y的权重向量$\bf{\pmb {c}}$满足

$\max {\rm{ }}\left[{{\rm cov}( {\pmb t}, {\pmb u})} \right]^2 = \mathop {\max }\limits_{_{\left|{ \pmb w} \right| = \left|{ \pmb c }\right| = 1} } \left[{ {\rm cov}( X{\pmb w}, Y{\pmb c})} \right]^2$ (2)

其中, ${{\rm cov}({\pmb t}, { \pmb u})}$为得分向量$ {\pmb {t}}$$ {\pmb {u}}$之间的样本协方差.数据集Y的第1个权重向量$ {\pmb {c}}_{\rm 1}$等于式(3) 特征值求解问题中最大特征值对应的单位特征向量[22]

$Y^{\rm T} XX^{\rm T} Y {\pmb c}=\lambda {\pmb c}$ (3)

CCA是一种基于相关系数最大化的降维技术, 常用来分析两个数据集之间的相关性.由式(4) 可以看出, PLS是一种带惩罚的CCA形式[23]

$\begin{aligned} &\max {\rm{ }}\left[{{\mathop{\rm cov}} ( {\pmb t}, {\pmb u})} \right]^2 = \\ &\mathop {\max }\limits_{_{\left| {\pmb w} \right|=\left| \pmb c \right|=1} } {\rm{ }}\left[{{\rm cov}( X{\pmb w}, Y{\pmb c)}} \right]^2 =\\ &\mathop {\max }\limits_{_{\left| \pmb w \right|=\left| {\pmb c} \right|=1} } {\rm{ }}\left\{ {{\rm var}( X {\pmb w})\left[{\rm corr}( X {\pmb w}, Y{\pmb c})\right]^2 {\rm var}( Y{\pmb c})} \right\} \\ \end{aligned}$ (4)

其中, ${\rm corr}$${\rm var}$分别表示样本相关性和方差.

在进行MPC性能评价研究时, 令数据集X为系统过程变量数据集, 数据集Y为系统输出变量数据集, 建立两个数据集之间的PLS关系模型.当系统性能发生变化时, 数据集XY以及两个数据集的相关性往往会发生变化.输出变量的方差直接体现系统性能的优劣, 因此, 为突出输出变量的重要性, 在评价系统性能时将过程变量数据集X的方差惩罚项${\rm var}(X {\pmb w})$从式(4) 中移除, 从而得到新的目标函数

$\begin{aligned} &\mathop {\max }\limits_{_{\left| {\pmb w} \right| = \left|{ \pmb c} \right| = 1} } {\rm{ }} \left\{ {\left[{\rm corr}( X{\pmb w}, Y{\pmb c})\right]^2 {\rm var}( Y{\pmb c})} \right\} =\\ &\mathop {\max }\limits_{_{\left| {\pmb w} \right| = \left| { \pmb c }\right| = 1} } \left\{ {\frac{{\left[{{\rm cov}( X{\pmb w}, Y{\pmb c})} \right]^2 }}{{{\rm var}( X{\pmb w})}}} \right\} \\ \end{aligned}%%%%%%%\label{dddss}, , , \ref{ddddsss}, \cite{ds}$ (5)

式中, 等式左边包含了两个数据集的典型相关分析和对输出数据集Y的PCA分解, 且可以转变为如下特征值求解问题[23-24]

$\begin{aligned} Y^{\rm T} { {X}}( X^{\rm T} X)^{{\rm{ - }}1} X^{\rm T} Y{\pmb c}=\,&\\ Y^{\rm T} \tilde { X} \tilde { X}{}^{\rm T} Y{\pmb c} =\,&\lambda { \pmb c} \\ \end{aligned}$ (6)

其中, $\tilde X=X{({X^{\rm{T}}}X)^{ -1/2}}$代表去相关和标准化后的过程变量数据集.定义Y$\tilde {X}$数据集的交叉积矩阵为

$ {\varPsi} = Y^{\rm T} \tilde { X}$ (7)

令特征矩阵$ {\varPhi}={\varPsi} {\varPsi} ^{\rm T} $, 则有

$ {\varPsi} {\varPsi} ^{\rm T} { \pmb c} = {\varPhi} {\pmb c}=\lambda {\pmb c}$ (8)

由式(5)~(8) 可知, 在新的目标函数下构造的交叉积矩阵${ {\varPsi}}$及其特征矩阵${ {\varPhi}}$既能够表征输出数据Y的方差特征, 也能够表征YX间的相关性特征.

考虑到模型预测控制系统的结构特点, 其内模结构描述如图 1所示[18]. $ \pmb G_{ p} $$\hat{ \pmb G}_ {p}$分别为实际过程和过程模型, ${ \pmb y} \in { {\bf R}}^{ m}$$\hat{ \pmb y} \in { {{\bf R}}}^m$分别为过程输出和模型输出, $ \pmb G_{ c} $${ \pmb u}_{ c} \in { {\bf R}}^{n}$分别为模型预测控制器和控制器输出, $ \pmb r$${ \pmb v}$分别为设定值和外部干扰, n为控制变量个数.则MPC控制系统的模型预测误差$ \pmb {ep}$

$\begin{aligned} { \pmb {ep}}(z) &= { \pmb y}(z) - {\hat{ \pmb y}}(z)= \\ &({ \pmb I} + \Delta { \pmb G}_{ c} )^{ - 1} \Delta { \pmb G}_{ c}{ \pmb r}(z) + ({ \pmb I} + \Delta { \pmb G}_{ c} )^{ - 1}{ \pmb v}(z) \\ \end{aligned}$ (9)

其中, $ \pmb I$为单位阵, ${\Delta}=\pmb G_{ p}-\hat{ \pmb G}_{ p} $为模型失配程度.由式(9) 可知, 模型预测误差受设定值$ \pmb r$、外部干扰$ \pmb v$和模型失配${\Delta}$三方面因素的影响.在稳态条件下, 设定值$ \pmb r$为常数, 此时, 预测误差的变化能够反映外部干扰和模型失配的变化.因此, 为了更为有效地监控MPC系统的性能, 将预测误差纳入到过程变量中得到增广的过程变量集X, 并计算矩阵$ {\varPhi}$.

图 1 模型预测控制系统内模结构 Figure 1 Schematic diagram of the internal model control structure for model predictive control

为了能够充分利用矩阵$ {\varPhi}$的特征信息分析当前工况相对于基准工况的性能变化, 采用非相似度分析方法分析该矩阵的特征变化, 对系统性能进行监控和诊断.

2 基于交叉积矩阵非相似度分析的性能监控及诊断 2.1 基于非相似度分析的性能监控

在性能监控时, 由工程师选取一段运行状态良好的工况数据作为基准工况, 假设矩阵${ {\varPsi} _{\rm ben}}$${ {\varPsi} _{\rm act} }$分别为控制系统基准工况和实际工况下的交叉积矩阵, 则容易看出${ {\varPhi} _{\rm ben} }$${ {\varPhi} _{\rm act} }$都是正定的, 定义它们的联合矩阵如下式

$ {\varPhi} _{\rm s} = {\varPhi} _{\rm ben} + {\varPhi} _{\rm act}$ (10)

则矩阵${ {\varPhi} _{\rm s}}$也是正定矩阵, 对${ {\varPhi} _{\rm s}}$进行特征值分解如下式

$ {\varPhi} _{\rm s} ={ M}_{\rm 0} {\varLambda} { M}_{\rm 0} ^{ - 1}$ (11)

其中, ${ M}_{\rm 0}$为正交矩阵, 其列向量为特征向量, ${ {\varLambda}}$为矩阵${ {\varPhi} _{\rm s}}$的特征值构成的对角阵.定义转换矩阵M如下式

$ M = M_{\rm 0} {\varLambda} ^{ - \frac{1}{2}}$ (12)

则矩阵M满足

${ M}^{\rm T} {\varPhi} _{\rm s} M = I$ (13)

利用M对矩阵${ {\varPhi} _{\rm ben} }$${ {\varPhi} _{\rm act} }$进行变换, 可得到各自的转换矩阵如下:

$\begin{array}{l} {\varPhi} ^\prime _{\rm ben} = M^{\rm T} {\varPhi} _{\rm ben} M, \quad {\varPhi} ^\prime _{\rm act} = M^{\rm T} {\varPhi} _{\rm act} M{\rm{ }}\\ \end{array}$ (14)

则由式(13) 和(14) 可推导出转换后的矩阵满足

$ M^{\rm T} {\varPhi} _{\rm s} M = {\varPhi} ^\prime _{\rm ben} + {\varPhi} ^\prime _{\rm act} = I$ (15)

对转换矩阵$ {\varPhi} ^\prime _{\rm ben}$$ {\varPhi} ^\prime _{\rm act}$进行特征值分解如下式

$\begin{array}{l} {\varPhi}^\prime _{\rm ben} \boldsymbol{\upsilon} _i^{\textrm{ben}} = \lambda _i^{\rm ben} \boldsymbol{\upsilon} _i^{\rm ben} {\rm{, }}\quad {\varPhi} ^\prime _{\rm act} \boldsymbol{\upsilon} _i^{\rm act} = \lambda _i^{\rm act} \boldsymbol{\upsilon} _i^{\rm act}\\ \end{array}$ (16)

其中, $\boldsymbol{\upsilon} _i^{\rm ben}$$\boldsymbol{\upsilon} _i^{\rm act}$分别为转换矩阵$\boldsymbol{\varPhi} ^\prime _{\rm ben}$$\boldsymbol{\varPhi} ^\prime _{\rm act}$的特征向量, $\lambda_i^{\rm ben}$$\lambda_i^{\rm act}$分别为转换矩阵$\boldsymbol{\varPhi} ^\prime _{\rm ben}$$\boldsymbol{\varPhi} ^\prime _{\rm act}$的特征值, $i=1, \cdots, m$.由式(15) 和(16) 可以得到:

$ {\varPhi} ^\prime _{\rm ben} \boldsymbol{\upsilon} _i^{\rm act} = (1 - \lambda _i^{\rm act} )\boldsymbol{\upsilon} _i^{\rm act}$ (17)

进而由式(16) 和(17) 得到[25]:

$\begin{array}{l} \boldsymbol{\upsilon} _i^{\rm ben} = \boldsymbol{\upsilon} _i^{\rm act} {\rm{, }} \quad \lambda _i^{\rm ben} + \lambda _i^{\rm act} = 1 \\ \end{array}$ (18)

由式(18) 可知, 经过矩阵M投影后的转换矩阵$ {\varPhi} ^\prime _{\rm ben}$$ {\varPhi} ^\prime _{\rm act}$有相同的特征向量, 即具有相同的主元方向.特征值$\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$关于0.5对称且$\lambda _i^{\rm ben}\in (0, 1)$, $\lambda _i^{\rm act}\in (0, 1)$, 因此, 特征值$\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$分布于图 2所示${\lambda _i^{\rm ben} + \lambda _i^{\rm act}=1 }$的线段上.通过比较转换矩阵$ {\varPhi} ^\prime _{\rm ben}$$ {\varPhi} ^\prime _{\rm act}$特征值的接近程度可以衡量矩阵$ {\varPhi}_{\rm ben}$$ {\varPhi}_{\rm act}$的非相似度.

两个矩阵整体的非相似度指标定义为[25]

$D = \frac{{4\sum\limits_{i = 1}^m {(\lambda _i^{\rm act} - 0.5)^2 } }}{m}$ (19)

i转换特征值对整体非相似度的贡献率为

$Con(i) = \frac{{(\lambda _i^{\rm act} - 0.5)^2 }}{{\sum\limits_{i = 1}^m {(\lambda _i^{\rm act} - 0.5)^2 } }} \times 100\%$ (20)

由式(20) 可知, 转换矩阵的特征值$\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$的值越接近0或1, 该特征值对两个矩阵整体非相似度贡献率越高, 相应主元描述了更多的非相似度信息. $\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$越接近0.5, 对非相似度贡献率越低, 相应主元描述了更少的非相似度信息, 则该主元描述了更多的相似度信息.转换矩阵特征值的分布位置与非相似度的关系如图 2所示.

图 2 转换矩阵特征值分布示意图 Figure 2 Distribution diagram of the eigenvalues of transformed matrixes

在进行性能评价时, 实际控制系统的性能越接近基准工况的性能, 实际工况数据的分布特征与基准工况数据的分布特征差异越小, $ {\varPhi} ^\prime _{\rm ben}$$ {\varPhi} ^\prime _{\rm act}$的特征值$\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$越靠近0.5;而另一方面, 实际控制系统的性能越偏离基准工况的性能, 实际工况数据与基准工况数据的分布特征差异越大, $\lambda _i^{\rm ben}$$\lambda _i^{\rm act}$的值越接近0或1[20].根据上述分析和图 2可知, 距离0.5最远的l个特征值对非相似度的贡献率最大, 包含了两个工况数据中最多的非相似度信息量, 即差异信息量, 从而能够更敏锐地反映出由于性能变化而导致的数据分布特征的变化, 因此可以选取距离0.5最远的l个特征值来评价系统性能的变化.利用两个矩阵的整体非相似度D作为目标准则来评估特征值对非相似度的重要程度, 如下式所示:

$4(\lambda _i^{\rm act} -0.5)^2 \ge D$ (21)

由式(21) 可知, $ \lambda _i^{\rm act} \in (0, \frac{{1 -\sqrt D }}{2}] \cup [\frac{{1 + \sqrt D }}{2}, 1)$.在特征值$\{ \lambda _1^{\rm act}, \cdots, \lambda _m^{\rm act} \}$中, 将满足式(21) 的特征值按照其对非相似度的贡献率由大到小的顺序重新排列, 得到相应的l个表征最大非相似度信息的特征值${\rm{\{ }}\tilde \lambda _1^{\rm act}, \cdots, \tilde \lambda _l^{\rm act} {\rm{\} }}$, $l < m$.则上述特征值所对应的主元为非相似度主元, 对应特征向量${\rm{\{ }}\tilde {\bf\pmb p}_{1}^{\rm act}, \cdots, \tilde {\bf\pmb p}_{l}^{\rm act} {\rm{\} }}$组成的子空间为非相似度子空间.

因此基于上述分析, 利用选取的l个特征值定义新的性能指标$\eta$.

$\eta = 1- \frac{{2\sum\limits_{i = 1}^l {\left| {\tilde \lambda _i^{\rm act}- 0.5} \right|} }}{l}$ (22)

传统的基于输出协方差行列式基准的性能指标为[16]

$\eta _{\det } = \frac{{\det ( {C} _{\rm ben} )}}{{\det ( {C} _{\rm act} )}} = \frac{{\prod\limits_{j = 1}^m {\varepsilon _j^{\rm ben} } }}{{\prod\limits_{j = 1}^m {\varepsilon _j^{\rm act} } }}$ (23)

其中, $ {C} _{\rm ben}$$ {C} _{\rm act}$分别为控制系统基准工况下和实际工况下的输出协方差矩阵, $\varepsilon _j^{\rm ben}$$\varepsilon _j^{\rm act}$分别为协方差矩阵$ {C} _{\rm ben}$$ {C} _{\rm act}$进行特征值分解后的特征值.

$\begin{array}{l} {C} _{\rm ben} = {\pmb W}_{\rm ben} { E}^{\rm ben} {\pmb W}_{\rm ben} ^{ - 1}{\rm{, }}\quad {C} _{\rm act} = {\pmb W}_{\rm act} { E}^{\rm act} {\pmb W}_{\rm act} ^{ - 1} \\ \end{array}$ (24)

其中, $\boldsymbol{W} _{\rm ben}$$\boldsymbol{W} _{\rm act}$分别为协方差矩阵$ {C} _{\rm ben}$$ {C} _{\rm act}$的特征向量, ${ E}^{\rm ben}={\rm diag}\{\varepsilon _1^{\rm ben}, \cdots, \varepsilon _{m}^{\rm ben} \}$${ E}^{\rm act}={\rm diag}\{\varepsilon _1^{\rm act}, \cdots, \varepsilon _{m}^{\rm act} \}$分别为$ {C} _{\rm ben}$$ {C} _{\rm act}$的特征值组成的对角矩阵.

通过将式(18)、(22) 与式(23)、(24) 对比可知, 上述基于交叉积矩阵非相似度分析的性能评价指标同时考虑了特征矩阵$ {\varPhi}$的特征值和特征向量信息, 而基于输出协方差行列式基准的性能评价指标只考虑了输出协方差矩阵的特征值信息而忽略了特征向量信息.此外, 所提性能指标更注重实际工况数据与基准工况数据的最大非相似度信息, 因而对系统性能的变化更为灵敏.在实时监控系统性能时, 用于获得性能指标的矩阵$ {\varPhi}$为统计量, 需一定长度的样本数据计算获得, 因此, 采用滑动时间窗口的形式实时计算该性能指标.将待监控的变量数据集转化为实时更新的形式:

$\left\{\!\!\!\begin{array}{l} \!{ \pmb Y}_{ r }(k)\! =\! \left[{ Y^{\rm T} (k-d + 1), \cdots, Y^{\rm T} (k-1), Y^{\rm T} (k)} \right]^{\rm T} \\ \!{ \pmb X}_{ r }(k)\! =\! \left[{ X^{\rm T} (k-d + 1), \cdots, X^{\rm T} (k-1), X^{\rm T} (k)} \right]^{\rm T} \\ \end{array} \right.$ (25)

其中, d为滑动时间窗口的大小, k为采样时刻.可以得到实时交叉积矩阵$ {\varPsi} _{ r} (k)$和对应的$ {\varPhi} _{ r} (k)$, 通过对$ {\varPhi} _{ r} (k)$$ {\varPhi} _{\rm ben}$的非相似度分析得到实时更新的MPC性能指标:

$\eta _{ r}(k) = 1 -\frac{{2\sum\limits_{i = 1}^{l_{ r}(k) } {\left| {\tilde \lambda _i^{\rm act} (k) -0.5} \right|} }}{{l_{ r} (k)}}$ (26)

其中, $\tilde \lambda _i^{\rm act}(k) \; (i=1, \cdots, l_{ r} (k))$, 为实时计算的$l_{ r}(k)$个表征最大非相似度的特征值, $l_{ r}(k)$为每个时刻实时提取的特征值数目, $\eta _{ r} \in (0, 1]$. $\eta _{ r}$越接近于1, 表明控制系统性能越接近于基准工况; $\eta _{ r}$越接近于0, 则表明控制系统性能越偏离基准工况.滑动窗口d的选取应保证指标$\eta _{ r}$具有适度平滑的变化趋势[26].如果d较小, 指标$\eta _{ r}$的变化趋势中波动较大, 可能因为不准确或敏感而导致误报.如果d较大, 指标$\eta _{ r}$的变化趋势较平缓, 对过程的动态变化响应较慢. d的大小可利用工况良好的基准数据离线确定.

2.2 性能诊断

在模型预测控制系统中, 引起MPC性能恶化的常见因素有模型失配、约束饱和、干扰特性变化等.这些因素都会影响过程数据的分布特征, 且不同性能恶化源会导致过程数据呈现不同的分布特征.

现代过程工业中, DCS系统和SCADA技术的广泛应用使得过程数据的实时采集和保存成为可能, 因此可以通过分析工业历史数据或实时采集的数据及其对应的实际工况, 建立已知的性能恶化模式库.假设$ {\varPsi} _f$为已知模式库中第f个性能恶化类别下控制系统的交叉积矩阵, $ {\varPhi} _f={\varPsi} _f {\varPsi} _f^{\rm T}$, 则该工况下的数据分布特征可由$ {\varPhi} _f$的特征向量与特征值表征.由第2.1节分析可知, 非相似度分析可将两个矩阵投影为具有相同特征向量的转换矩阵, 则可以利用两个转换矩阵特征值之间的差异来度量两个数据集的相似度[19].假设待诊断工况下的交叉积矩阵为$ {\varPhi} _x$, $ {\varPhi} _x={\varPsi} _x {\varPsi} _x^{\rm T}$, 因此矩阵$ {\varPhi} _f$$ {\varPhi} _x$之间的相似度指标定义为

${{SI}}_{xf} = 1 - \frac{{\sum\limits_{i = 1}^m {(\lambda _i^x - \lambda _i^f )^2 } }}{m}$ (27)

其中, $\lambda _i^x$$\lambda _i^f$分别为矩阵$ {\varPhi} ^\prime _x$$ {\varPhi} ^\prime _f$的特征值, $f=1, 2, \cdots, g $, g为性能恶化类别总数. ${{SI}}_{xf}$越接近1, 两个数据集的相似性越大; ${{SI}}_{xf}$越接近0, 两个数据集的相似性越小.相似度指标${{SI}}_{xf}$可以定量地描述两个数据集的相似程度, 可以用于性能诊断.

在进行性能诊断时, 首先根据先验知识对第f个性能恶化工况${\rm {CL}}_{f}$下的数据进行标准化并计算交叉积矩阵$ {\varPsi} _f$$ {\varPhi} _f$, 建立控制系统的性能模式库$ \{ {\rm {CL}}_{f}, { {\varPhi}} _{f}\}$.当监控到性能下降时, 计算$ {\varPhi} _x$与性能模式库中$ {\varPhi} _f$的相似度指标${{SI}}_{xf}$, 根据式(28) 所示判别策略确定性能恶化类别.

$\begin{align*} f_0&= \\ &\begin{cases} \arg \max ({ {SI}}_{xf} ), &\!\!\exists { {SI}}_{xf} \ge \alpha, f \in \{ 1, 2, \cdots, g\} \\ \mbox{未知故障}, &\!\!\mbox{否则}\\ \end{cases} \end{align*}$ (28)

其中, $f_0$为判别得到的性能恶化类别, $\alpha$为相似度阈值.较大的阈值意味着待诊断的性能数据需要与性能模式库中的数据在分布上有更高的相似性, 从而易将已知性能故障判别为未知故障; 较小的阈值意味着判别条件放松, 则易将未知性能故障判别为已知故障, 因此需根据具体过程的经验知识选取适当的阈值[27-28].当诊断结果为未知的性能恶化类别时, 需通过现场工程师和专家分析确定性能恶化原因, 建立该性能模式类别, 对模式库进行及时更新和补充.

由式(26) 和(27) 可知, 本文提出的性能监控和诊断方法均基于非相似度分析方法.由于在对矩阵$ { \varPhi}$的非相似度分析过程中充分挖掘了矩阵的特征值和特征向量信息, 同时提取了描述最大非相似度的l个特征值用于性能监控, 有助于进一步提高性能监控效果.

3 仿真研究

Wood-Berry二元精馏塔是甲烷-水的精馏塔模型, 为典型的具有较大纯滞后的多输入多输出系统[20], 输出$X_D$$X_B$分别为塔顶馏出物浓度和塔底液相浓度, 控制量RS分别为塔顶回流量和塔底再沸蒸汽量, 进料流量F为不可测干扰变量, 该过程模型为

$\begin{aligned}%{array}{l} \!\!\left[\!\! \begin{array}{l} X_D (s) \\ X_B (s) \\ \end{array}\!\! \right] \!= \! &\left[{\begin{aligned} \frac{{12.8\textrm{e}^{-s} }}{{16.7s + 1}}&\quad \frac{{-18.9\textrm{e}^{-3s} }}{{21.0s + 1}}\\ \frac{{6.6\textrm{e}^{ - 7s} }}{{10.9s + 1}}&\quad \frac{{ - 19.4\textrm{e}^{ - 3s} }}{{14.4s + 1}}\\ \end{aligned}} \right]\!\!\left[\!\!\begin{array}{l} R(s)\\ S(s)\\ \end{array} \!\!\right]\!+\! \\ &\left[\begin{aligned} \frac{{3.8\textrm{e}^{-8s} }}{{14.9s + 1}} \\ \frac{{4.9\textrm{e}^{-3s} }}{{13.2s + 1}} \\ \end{aligned} \right]F(s) \\ \end{aligned}$ (29)

假设干扰为满足$F\sim \textrm{N}(0, 0.1^2)$的白噪声, 采样步长为1 min, 设计动态矩阵控制器, 控制器的二次型优化目标函数为

$\begin{align} \!\!J\!(k)\!\!=\!\!&\sum\limits_{i = 1}^P\!{[{\bf\pmb r}(k\!+\!i)\!-\!\hat {\bf\pmb y}(k\!+\!i)]^{\rm T} } { H}[{\bf\pmb r}(k\!+\!i)\!-\!\hat {\bf\pmb y}(k\!+\!i)]+\nonumber\\ &\sum\limits_{i = 1}^M \!{\Delta {\bf\pmb u}_c (k\!+\!i\!-\!1)^{\rm T} } { L}\Delta {\bf\pmb u}_c (k\!+\!i\!-\!1) \end{align}$ (30)

其中, $\Delta {\bf\pmb u}_c$为控制输出增量, PM分别为预测时域和控制时域, HL分别为预测输出误差和控制器输出增量的加权系数矩阵.控制器参数为$P=10$, $M=1$, $ { H}=I $, $ { L}=I $, 系统输出设定值为$r_{X_D }=1$, $r_{X_B }=0.4$, 控制器预测误差为$ep_D, ep_B$.则控制系统的增广过程变量数据集为$ {\bf\pmb X}{\rm{=[}}R, S, ep_D, ep_B]$, 输出变量数据集为$ {\bf\pmb Y}=[X_D, X_B]$.仿真控制系统的良好运行工况, 采集1 000个观测样本作为基准数据, 分别模拟过程模型失配($\rm {CP}_1$)、干扰标准差变化($\rm {CP}_2$) 和输出约束饱和($\rm {CP}_3$) 3种MPC系统性能恶化工况, 具体参数设置见表 1.对每种恶化工况分别采集1 000个样本作为观测数据建立Wood-Berry精馏塔的MPC性能恶化模式库$\{ \rm {CL}_1, {\varPhi} _1 \}$$\{ \rm {CL}_2, {\varPhi} _2 \}$$\{ \rm {CL}_3, {\varPhi} _3 \}$.同时, 选取滑动窗口$d=1\, 000$, 采用滑窗的形式实时计算基于输出协方差行列式的性能指标$\eta _{\det }$ [16]和基于输出协方差非相似度的性能指标$\eta _{\rm dissim}$[20], 并比较所提方法与上述两种方法的性能监控效果.

表 1 恶化性能类别及参数设置 Table 1 Classes and parameters of performance deterioration
3.1 性能评价

按照表 1仿真3种性能恶化工况, $t=2\, 000\min$时, 设置MPC系统性能恶化, $t=4\, 000\min$时, 使系统恢复正常.为直观对比所提方法下的矩阵和输出协方差矩阵在性能恶化下的变化差异, 令$ \bar { {\varPhi}}={ \varPhi} /(N -1)={ \varPsi \varPsi} ^{\rm T} /(N -1)$, $ {C}=Y^{\rm{T}} Y{{/(N -1)}}$, 可分别画出矩阵$\bar { {\varPhi}}$和矩阵C在基准工况和不同恶化工况下所构成的椭圆(高维数据时为超椭球体).

当模型失配时, 矩阵C$\bar { {\varPhi}}$构成的椭圆如图 3所示, 由图 3可知, 矩阵$\bar { {\varPhi}}$在基准工况下和该性能恶化工况下的差异更显著, 因此所提方法对性能变化的检测更为灵敏.该工况下的实时性能监控曲线如图 4所示, 由图 4可知, 在$ t=0 \sim 2\, 000\min$的正常工况阶段, 性能指标$\eta _{ r}$$\eta _{\rm dissim}$曲线均在$0.9\sim1$的区域波动, $\eta _{\det } $在1附近波动, 表明此时MPC系统具有良好的控制性能.当失配发生后, 性能指标$\eta _{ r}$$\eta _{\rm dissim}$均明显小于1, 检测到MPC系统性能的恶化, 但$\eta _{ r}$下降更为显著, 表明该指标对性能恶化的监控效果更好, 而指标$\eta _{\det }$却仍在1附近波动, 没有监控到性能的下降.这是因为$\eta _{\det }$只考虑了输出协方差矩阵的特征值而忽略了特征向量, 从而对某些情况导致的性能变化给出了错误的评价结果[20], 而基于非相似度分析的性能指标同时考虑了特征值和特征向量, 从而很好地解决了这一问题.

图 3 模型失配时矩阵C与矩阵$\bar { {\varPhi}}$所构成的椭圆 Figure 3 Ellipses figure formed by matrices C and $\bar{\varPhi}$ with model mismatch
图 4 模型失配时的实时监控曲线 Figure 4 Real-time monitoring curves of model mismatch

表 2设置不同的传递函数矩阵首行增益, 取性能恶化工况下的1 000个样本分别计算性能指标$\eta _{\det }$$\eta _{\rm dissim}$$\eta _{ r}$, 结果如表 2所示.由表 2可知, 当模型失配时, 指标$\eta _{ r}$$\eta _{\rm dissim}$的值明显小于1, 表明都检测到系统性能的下降, 且模型失配越严重, 两个指标值越小, 表明性能恶化趋于严重.而基于协方差行列式的指标$\eta _{\det }$随着失配程度的增大反而大于1, 即表明系统的性能提升, 与实际工况不符.

表 2 不同程度模型失配下的性能指标 Table 2 Performance under different degree of model mismatch

当干扰特征变化时, 矩阵C$\bar { {\varPhi}}$所构成的椭圆如图 5所示.由图 5可知, 由于干扰标准差变化没有改变两数据集之间的相关性, 矩阵$\bar { {\varPhi}}$在基准工况下和该恶化工况下的差异与$ {C}$的差异相当.该工况下的实时性能监控曲线如图 6所示, 由图 6可知, 3种性能指标均检测到了性能下降, 但指标$\eta _{\det }$的性能监控效果最优, $\eta _{ r}$的监控效果次优, $\eta _{\rm dissim}$的最差.按表 3所示设置不同干扰标准差, 由表 3可知, 三个指标均检测到系统性能的恶化, 且随着标准差的增大, 各指标值变小, 表明性能恶化趋于严重.

图 5 干扰特征变化时矩阵C与矩阵$\bar { {\varPhi}}$所构成的椭圆 Figure 5 Ellipses figure formed by matrices C and $\bar { {\varPhi}}$ with changing disturbance
图 6 干扰特征变化时的实时监控曲线 Figure 6 Real-time monitoring curves of changing disturbance
表 3 不同干扰标准差下的性能指标 Table 3 Performance under different standard deviation of disturbance

当输出约束饱和时, 矩阵C$\bar { {\varPhi}}$所构成的椭圆如图 7所示.由图 7可知, 矩阵$\bar { {\varPhi}}$在基准工况下和该恶化工况下的差异更显著.该工况下的实时性能监控曲线如图 8所示, 由图 8可知, 3种性能评价指标均检测到了性能下降.但对于指标$\eta _{\det }$, 在$t=2\, 000\min$性能下降时, 由于在长度为$d=1\, 000$的滑窗中包含了性能恶化瞬间过程输出的动态变化, 即滑窗中的样本数据既有正常工况数据也有恶化工况数据, 因而导致指标$\eta _{\det }$迅速减小.在$t=3\, 000\min$后; 滑窗中的样本数据全部为恶化数据, 因此指标$\eta _{\det }$有所上升并趋于平稳, 且此时$\eta _{\det }$明显小于1表明系统性能恶化; 在$t=4\, 000\min$时, 系统恢复正常工况, 同样由于在滑窗中包含了过程输出的动态变化导致指标$\eta _{\det }$迅速减小, 并在$t=5\, 000\min$后迅速上升并趋于平稳.可知指标$\eta _{\det }$对输出约束饱和导致的性能恶化的监控效果会受到滑窗长度的影响, 在性能恶化和恢复时均会产生较大的下降波动, 从而在性能恢复的初期易造成性能误报, 而指标$\eta_{ r}$$\eta_{\rm dissim}$则没有受到影响.同时, 指标$\eta_{ r}$比指标$\eta_{\rm dissim}$的值更小, 监控效果最好.按表 4所示设置不同的输出约束, 由表 4仿真结果可以看出, 三个指标都检测到系统性能的恶化, 且随着约束饱和程度增强, 各指标值变小, 表明MPC系统性能越差.

图 7 约束饱和时矩阵C与矩阵$\bar { {\varPhi}}$所构成的椭圆 Figure 7 Ellipses figure formed by matrices C and $\bar { {\varPhi}}$ with output constraint saturation
图 8 输出约束饱和时的实时监控曲线 Figure 8 Real-time monitoring curves of output constraint saturation
表 4 不同程度输出约束饱和下的性能指标 Table 4 Performance under different degree of output constraint saturation

为进一步比较指标$\eta_{ r}$$\eta_{\rm dissim}$对性能监控的灵敏度.设置性能恶化源为较小程度的模型失配,假设过程传递函数首行增益在$ t=2\, 000\min$时由$(12.8, {\rm{ \ -}}18.9)$变为$(15.36, {\rm{ \ -}}22.68)$, $t=4\, 000\min$时系统恢复正常, 此时的实时性能监控曲线如图 9所示.由图 9可知, 当发生较小程度的模型失配时, 在发生失配的时间区间内, 指标$\eta_{ r}$能够显著地监控到性能的变化, 而指标$\eta_{\rm dissim}$却没有监控到性能的改变.结果表明本文所提方法相比于传统基于非相似度的性能指标对性能变化的监控更为灵敏.

图 9 较小模型失配程度下的实时监控曲线 Figure 9 Real-time monitoring curves under small extent of model mismatch

综上可知, 在三种性能指标中, 指标$\eta_{ r}$对模型失配和输出约束饱和导致的性能下降监控效果最优, 对于噪声标准差变化导致的性能下降监控效果次优; $\eta _{\det }$对于噪声标准差变化导致的性能恶化监控效果最优, 但其没有监控到模型失配导致的性能恶化, 且当约束饱和导致性能恶化时, 在性能恢复初期易造成性能误报; $\eta_{\rm dissim}$没有监控到较小程度的模型失配导致的性能恶化.因此, 本文所提性能指标$\eta_{ r}$在性能监控效果上整体优于传统性能指标$\eta _{\det }$$\eta_{\rm dissim}$.

3.2 性能诊断

表 2表 4不同类型不同程度性能恶化下的数据作为待诊断数据, 计算待诊断数据与模式库$ \{{\rm {CL}}_f, {\varPhi} _f \}$的相似度指标${{SI}}_{xf}$, 取$\alpha=0.85$, 性能诊断结果见表 5.由表 5可知, 对于模型失配$\rm {CP}_1$, 相似度指标${{SI}}_{x1}$最接近于1, 表明其与$\rm {CL}_1$性能模式相似度最高, 从而诊断出当前性能恶化源为模型失配.同样, 对于干扰标准差变化$\rm {CP}_2$和输出约束饱和$\rm {CP}_3$, 也能够正确诊断出性能恶化源.因此, 该方法能够正确识别MPC系统的性能下降源.

表 5 各类测试数据的性能诊断结果 Table 5 Performance diagnosis results corresponding various test data
4 结论

本文提出了基于偏最小二乘交叉积矩阵非相似度分析的性能监控与诊断方法.利用PLS交叉积矩阵考虑了系统过程变量与输出变量相关性以及MPC控制器的预测误差信息, 通过非相似度分析的方法充分利用了矩阵中特征值和特征向量信息.同时, 定义了非相似度主元的选取方法, 并通过所选取的l个包含最大非相似度信息的特征值构造性能指标, 更有效地利用非相似度信息监控系统性能的变化.所提方法比较有效地解决了MPC系统的性能评价和诊断问题.在Wood-Berry塔上的仿真结果表明, 所提方法能够更加显著地检测到MPC系统的性能下降, 有效诊断出性能恶化原因.

参考文献
1 Xi Yu-Geng, Li De-Wei, Lin Shu. Model predictive control-status and challenges. Acta Automatica Sinica, 2013, 39 (3): 222–236.
( 席裕庚, 李德伟, 林姝. 模型预测控制-现状与挑战. 自动化学报, 2013, 39 (3): 222–236. DOI:10.1016/S1874-1029(13)60024-5 )
2 Ellis M, Christofides P D. Economic model predictive control of nonlinear time-delay systems:closed-loop stability and delay compensation. AIChE Journal, 2015, 61 (12): 4152–4165. DOI:10.1002/aic.v61.12
3 Müller M A, Angeli D, Allgöwer F. On the performance of economic model predictive control with self-tuning terminal cost. Journal of Process Control, 2014, 24 (8): 1179–1186. DOI:10.1016/j.jprocont.2014.05.009
4 Wang L, Li N, Li S Y. Performance monitoring of the data-driven subspace predictive control systems based on historical objective function benchmark. Acta Automatica Sinica, 2013, 39 (5): 542–547.
5 Harris T J. Assessment of control loop performance. The Canadian Journal of Chemical Engineering, 1989, 67 (5): 856–861. DOI:10.1002/cjce.v67:5
6 Tian Xue-Min, Luo Zhi-Fen, Wang Ping. LQG-based sensitivity analysis and tuning guidelines in economic performance assessment of predictive controller. Acta Automatica Sinica, 2013, 39 (10): 1735–1740.
( 田学民, 罗芝芬, 王平. 基于LQG基准的预测控制器经济敏感度分析及调节准则. 自动化学报, 2013, 39 (10): 1735–1740. DOI:10.3724/SP.J.1004.2013.01735 )
7 Liu C Y, Huang B, Wang Q L. Control performance assessment subject to multi-objective user-specified performance characteristics. IEEE Transactions on Control Systems Technology, 2011, 19 (3): 682–691. DOI:10.1109/TCST.2010.2051669
8 Ge Z Q, Song Z H. An overview of conventional MSPC methods. Multivariate Statistical Process Control. London: Springer, 2013: 5-11.
9 Yan Z B, Chan C L, Yao Y. Multivariate control performance assessment and control system monitoring via hypothesis test on output covariance matrices. Industrial & Engineering Chemistry Research, 2015, 54 (19): 5261–5271.
10 Jelali M. Statistical process control. Control Performance Management in Industrial Automation. London: Springer, 2013: 209-217.
11 Zhang Q, Li S Y. Performance monitoring and diagnosis of multivariable model predictive control using statistical analysis. Chinese Journal of Chemical Engineering, 2006, 14 (2): 207–215. DOI:10.1016/S1004-9541(06)60060-8
12 AlGhazzawi A, Lennox B. Model predictive control monitoring using multivariate statistics. Journal of Process Control, 2009, 19 (2): 314–327. DOI:10.1016/j.jprocont.2008.03.007
13 Chen J H, Wang W Y. PCA-ARMA-based control charts for performance monitoring of multivariable feedback control. Industrial & Engineering Chemistry Research, 2010, 49 (5): 2228–2241.
14 Shang C, Huang B, Fan Y, Huang D X. Slow feature analysis for monitoring and diagnosis of control performance. Journal of Process Control, 2016, 39 : 21–34. DOI:10.1016/j.jprocont.2015.12.004
15 Das L, Srinivasan B, Rengaswamy R. Multivariate control loop performance assessment with Hurst exponent and Mahalanobis distance. IEEE Transactions on Control Systems Technology, 2016, 24 (3): 1067–1074. DOI:10.1109/TCST.2015.2468087
16 Yu J, Qin S J. Statistical MIMO controller performance monitoring. Part I:data-driven covariance benchmark. Journal of Process Control, 208, 18 (3-4): 277–296.
17 Yu J, Qin S J. Statistical MIMO controller performance monitoring. Part II:performance diagnosis. Journal of Process Control, 2008, 18 (3-4): 297–319. DOI:10.1016/j.jprocont.2007.09.003
18 Tian X M, Chen G Q, Chen S. A data-based approach for multivariate model predictive control performance monitoring. Neurocomputing, 2011, 74 (4): 588–597. DOI:10.1016/j.neucom.2010.09.018
19 Tian Xue-Min, Shi Ya-Jie, Cao Yu-Ping. Real-time performance monitoring of MPC based on covariance index prediction. Acta Automatica Sinica, 2013, 39 (5): 658–663.
( 田学民, 史亚杰, 曹玉苹. 基于协方差指标预测的MPC实时性能监控. 自动化学报, 2013, 39 (5): 658–663. )
20 Li C, Huang B, Zheng D, Qian F. Multi-input-multi-output (MIMO) control system performance monitoring based on dissimilarity analysis. Industrial & Engineering Chemistry Research, 2014, 53 (47): 18226–18235.
21 Wold H. Soft modeling by latent variables:the nonlinear iterative partial least squares approach. Perspectives in Probability and Statistics:Papers in Honour of M.S. Bartlett. New York: Academic Press, 1975: 520-540.
22 Höskuldsson A. PLS regression methods. Journal of Chemometrics, 1988, 2 (3): 211–228. DOI:10.1002/(ISSN)1099-128X
23 Rosipal R, Trejo L J, Matthews B. Kernel PLS-SVC for linear and nonlinear classification. In:Proceedings of the 20th International Conference on Machine Learning. Washington D.C., USA, 2003. 640-647
24 Barker M, Rayens W. Partial least squares for discrimination. Journal of Chemometrics, 2003, 17 (3): 166–173. DOI:10.1002/(ISSN)1099-128X
25 Kano M, Hasebe S, Hashimoto I, Ohno H. Statistical process monitoring based on dissimilarity of process data. AIChE Journal, 2002, 48 (6): 1231–1240. DOI:10.1002/(ISSN)1547-5905
26 Yuan Q L, Lennox B, McEwan M. Analysis of multivariable control performance assessment techniques. Journal of Process Control, 2009, 19 (5): 751–760. DOI:10.1016/j.jprocont.2008.10.001
27 Tan S, Wang F L, Peng J, Chang Y Q, Wang S. Multimode process monitoring based on mode identification. Industrial & Engineering Chemistry Research, 2012, 51 (1): 374–388.
28 Zheng Y, Qin S J, Wang F L. PLS-based similarity analysis for mode identification in multimode manufacturing processes. In:Proceedings of the 9th IFAC Symposium on Advanced Control of Chemical Processes. Whistler, British Columbia, Canada:Elsevier, 2015, 48(8):777-782