CN104698976A

Movatterモバイル変換

Info

Publication number: CN104698976A
Application number: CN201410811191.XA
Authority: CN
Inventors: 李丽娟; 王凯; 张晓晓; 周梦迪
Original assignee: Nanjing Tech University
Current assignee: Nanjing Tech University
Priority date: 2014-12-23
Filing date: 2014-12-23
Publication date: 2015-06-10
Anticipated expiration: 2034-12-23
Also published as: CN104698976B

Abstract

Translated fromChinese

本发明公开了一种预测控制模型性能下降的深度诊断方法，利用生产运行数据计算过程扰动信号，然后利用生产数据和预测控制器设计阶段的阶跃响应系数计算模型预测误差，通过二者构造的模型性能指标判断模型整体性能优劣。对于性能恶化的模型，进一步用逐一去除输入变量的方法计算新的模型性能指标，通过性能指标的变化情况判断去除的输入对应的子模型性能，从而实现对每一个子模型性能进行监控。本发明仅利用生产过程数据和设计数据，不仅能够对多变量预测控制模型整体性能给出评估，更能对所有输入对应的子模型性能进行评估，为工程师进行控制系统维护给出建议，能够大幅降低预测控制器维护成本。

The invention discloses a method for in-depth diagnosis of performance degradation of a predictive control model. The process disturbance signal is calculated by using the production operation data, and then the model prediction error is calculated by using the production data and the step response coefficient in the design stage of the predictive controller. The model performance index judges the overall performance of the model. For models with deteriorating performance, the new model performance index is calculated by removing input variables one by one, and the performance of the sub-model corresponding to the removed input is judged by the change of the performance index, so as to monitor the performance of each sub-model. The present invention only uses the production process data and design data, not only can evaluate the overall performance of the multivariable predictive control model, but also can evaluate the performance of sub-models corresponding to all inputs, and give suggestions for engineers to maintain the control system, which can greatly Reduce predictive controller maintenance costs.

Description

Translated fromChinese

一种预测控制模型性能下降的深度诊断方法A Deep Diagnosis Method for Predictive Control Model Performance Degradation

技术领域technical field

本发明属于工业预测控制系统性能监控领域，涉及基于模型性能指标及移除输入变量深度检验的预测模型监控技术。The invention belongs to the field of performance monitoring of industrial predictive control systems, and relates to a predictive model monitoring technology based on model performance indicators and in-depth inspection of removed input variables.

背景技术Background technique

由于生产过程动态特性随时间发生变化，工业模型预测控制(MPC)系统通常在投运不到一年就需要进行维护，维护过程需要对象重新建模，周期长成本高，而且通常需要生产过程停车，造成很大经济损失。况且，在维护之前判断是否有必要对整个对象重新建模，或者重新建模是否确定会提高整个控制性能，目前缺乏有效的分析方法。为了避免重新建模和控制器维护的不必要投入，有必要研究一种方法来判断对象哪部分模型性能出现了下降，以及下降到了什么程度。As the dynamic characteristics of the production process change over time, the industrial model predictive control (MPC) system usually needs to be maintained less than one year after it is put into operation. The maintenance process requires object remodeling, which has a long cycle and high cost, and usually requires the shutdown of the production process. , causing great economic losses. Moreover, there is currently a lack of effective analysis methods to judge whether it is necessary to remodel the entire object before maintenance, or whether remodeling will definitely improve the overall control performance. In order to avoid unnecessary investment in remodeling and controller maintenance, it is necessary to study a method to judge which part of the object model performance has declined, and to what extent.

控制性能监控可以给企业工程师们对MPC是否需要维护给出建议。控制性能监控的一般方法是Harris提出的极小方差性能指标，把实际的输出方差和理论最小方差进行比较，从而判断出实际控制器性能。后续相继出现了多输入多输出系统，时变系统，有约束系统的极小方差性能评估指标。由于多输入多输出系统的极小方差评估技术需要计算关联矩阵，而关联矩阵的计算往往需要获取比较困难的时滞和脉冲响应系数。为避免关联矩阵计算，出现了一些更实用的评估指标，如线性二次型高斯基准，用户定义的基准，设计指标和实际指标的对比基准，历史最好性能基准，等等。Control performance monitoring can advise enterprise engineers on whether MPCs need maintenance. The general method of control performance monitoring is the minimum variance performance index proposed by Harris, which compares the actual output variance with the theoretical minimum variance to judge the actual controller performance. Subsequently, the minimum variance performance evaluation indicators for multi-input and multi-output systems, time-varying systems, and constrained systems appeared one after another. Because the minimum variance evaluation technique of the MIMO system needs to calculate the correlation matrix, and the calculation of the correlation matrix often needs to obtain the time delay and impulse response coefficients which are difficult. In order to avoid the calculation of the correlation matrix, some more practical evaluation indicators have emerged, such as linear quadratic Gaussian benchmarks, user-defined benchmarks, comparison benchmarks between design indicators and actual indicators, historical best performance benchmarks, and so on.

控制性能监控技术给出了评估整个控制系统性能的方法，但是，控制性能下降可能由于控制器性能下降，预测模型不匹配，新的干扰出现以及执行器非线性等多种因素。近年来很多研究工作集中在对控制器性能下降的原因进行诊断。如文献Xuemin Tian,GongquanChen,Sheng Chen.A data-based approach for multivariate model predictive control performance monitoring.Neurocomputing 74(2011)588–597，提出采用数据驱动的分类方法对控制性能下降的种类进行甄别，文献Keck Voon Ling,Weng Khuen Ho,Yong Feng,and Bingfang Wu.Integral-Square-Error Performance of Multiplexed Model Predictive Control.IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS,2011,7(2):196-203，提出了判定控制器本身性能下降的ISE指标，从而识别出控制系统性能下降是否由控制器引起。针对预测模型性能下降的诊断也有一些研究工作，如通过设计的控制误差基准判断模型性能(如Abhijit S.Badwe,Rohit S.Patwardhan,Sirish L.Shah,Sachin C.Patwardhan,Ravindra D.Gudi.Quantifying the impact of model-plant mismatch on controller performance.Journal of Process Control 20(2010)408–425)，通过测试的过程频率响应基准判定模型性能(G.Ji,K.Zhang,Y.Zhu.A method of MPC model error detection.Journal of Process Control,22(2012)635–642)，或通过预测模型余量和实际干扰量的对比判定模型性能等方法(Zhijie Sun,S.Joe Qin,Ashish Singhal,Larry Megan.Performance monitoring of model-predictive controllers via model residual assessment.Journal of Process Control,23(2013)473-482)。Control performance monitoring techniques give a way to evaluate the performance of the entire control system, however, control performance degradation may be due to various factors such as controller performance degradation, mismatch of predictive models, emergence of new disturbances, and actuator nonlinearity. Much research work in recent years has focused on diagnosing the causes of controller performance degradation. For example, Xuemin Tian, Gongquan Chen, Sheng Chen.A data-based approach for multivariate model predictive control performance monitoring.Neurocomputing 74(2011)588–597, proposed to use data-driven classification method to identify the types of control performance decline, the document Keck Voon Ling, Weng Khuen Ho, Yong Feng, and Bingfang Wu. Integral-Square-Error Performance of Multiplexed Model Predictive Control. IEEE TRANSACTIONS ON INDUSTRIAL INFORMATICS, 2011, 7(2): 196-203, proposed to determine the performance of the controller itself Decreased ISE metrics to identify if control system performance degradation is caused by the controller. There are also some research works on the diagnosis of prediction model performance degradation, such as judging model performance through designed control error benchmarks (such as Abhijit S.Badwe, Rohit S.Patwardhan, Sirish L.Shah, Sachin C.Patwardhan, Ravindra D.Gudi.Quantifying the impact of model-plant mismatch on controller performance. Journal of Process Control 20(2010) 408–425), through the test process frequency response benchmark to determine model performance (G.Ji, K.Zhang, Y.Zhu.A method of MPC model error detection.Journal of Process Control, 22(2012)635–642), or by comparing the predicted model margin and the actual interference to determine the performance of the model (Zhijie Sun, S. Joe Qin, Ashish Singhal, Larry Megan .Performance monitoring of model-predictive controllers via model residual assessment.Journal of Process Control,23(2013)473-482).

对于预测模型性能的评估，上述研究工作存在以下问题：(1)实际工业过程大都是多输入多输出(MIMO)的多变量系统，每个输出变量对应每个输入变量的模型都是一个单输入单输出(SISO)模型，目前的研究工作只给出了整个模型性能是否下降的评估结论，但工程师们更感兴趣的是具体哪个变量的哪个模型性能下降导致了控制性能下降，从而在控制系统维护时减少工作量，节省成本；(2)性能指标的定义往往需要难以获取的过程机理知识或要进行对象特性测试，难以向实际应用推广。因此，需要一种只利用生产过程信息而不影响生产过程，能够在线分析出模型整体性能是否下降以及诊断具体哪个子模型性能下降的方法，在系统维护之前为工程师们给出维护建议，减少不必要的维护成本。For the evaluation of the performance of the prediction model, the above research work has the following problems: (1) Most of the actual industrial processes are multi-variable systems with multiple input and multiple output (MIMO), and each output variable corresponds to each input variable. The model is a single input Single output (SISO) model, the current research work only gives the evaluation conclusion of whether the performance of the whole model is degraded, but engineers are more interested in which model performance degrades in which variable leads to the degraded control performance, so that in the control system Reduce workload and save costs during maintenance; (2) The definition of performance indicators often requires difficult-to-obtain process mechanism knowledge or object characteristic testing, which is difficult to promote to practical applications. Therefore, there is a need for a method that only uses the production process information without affecting the production process, can analyze online whether the overall performance of the model is degraded and diagnoses which sub-model performance degrades, and provides maintenance suggestions for engineers before system maintenance, reducing Unnecessary maintenance costs.

发明内容Contents of the invention

本发明的目的在于提供一种预测控制模型性能监控的方法，该方法只利用生产过程运行数据而不影响生产过程正常运行，能够给出模型整体性能评估指标，并且在整体性能下降时进一步诊断出具体哪个子模型性能下降以及下降程度，为预测控制系统的维护减少工作量，节省成本。The purpose of the present invention is to provide a method for monitoring the performance of predictive control models, which only uses the operating data of the production process without affecting the normal operation of the production process, can give the overall performance evaluation index of the model, and can further diagnose when the overall performance declines Specifically which sub-model performance declines and the degree of decline reduces the workload and costs for the maintenance of the predictive control system.

为解决上述技术问题，本发明提供一种预测控制模型性能下降的深度诊断方法，包括以下步骤：In order to solve the above technical problems, the present invention provides a method for in-depth diagnosis of predictive control model performance degradation, comprising the following steps:

步骤1：对于某多输入多输出过程，k时刻某输出变量y(k)对应于N_u个输入变量u_j(k)(j＝1,2，…，N_u)，e^o(k)为过程扰动信号，p为时间窗长度，m为输出变量阶次，n为输入变量阶次，定义时间窗p内输出向量y_p(k)，输入向量u_jp(k)(j＝1,2，…，N_u)，过程扰动信号向量以及由过去时刻输入输出变量构造的联合向量Step 1: For a certain MIMO process, an output variable y(k) at time k corresponds to N_u input variables u_j (k) (j=1,2,...,N_u ), e^o (k) is the process disturbance signal, p is the length of the time window, m is the order of the output variable, n is the order of the input variable, define the output vector y_p (k) in the time window p, and the input vector u_jp (k) (j=1, 2,...,N_u ), process disturbance signal vector and the joint vector constructed from the input and output variables at past moments

步骤2：由输出向量y_p(k)和输入输出联合向量计算过程扰动信号向量序列 $e_{p}^{o} (k) = [\begin{matrix} e^{o} (k) & e^{o} (k - 1) & . . . & e^{o} (k - p) \end{matrix}];$ Step 2: From the output vector y_p (k) and the input-output joint vector Computing Process Perturbation Signal Vector Sequence $e_{p}^{o} (k) = [\begin{matrix} e^{o} (k) & e^{o} (k - 1) & . . . & e^{o} (k - p) \end{matrix}];$

步骤3：利用输入变量u_j(k)、输出变量y(k)的常规运行数据和预测模型系数a_ij,分别计算k至k-p时刻的模型预测误差，得到模型预测误差序列e_p(k)＝[e(k) e(k-1) ... e(k-p)]；Step 3: Using the input variable u_j (k), the output variable y(k)’s regular operating data and the prediction model coefficient a_ij , respectively calculate the model prediction error from k to kp, and obtain the model prediction error sequence e_p (k) =[e(k) e(k-1) ... e(kp)];

步骤4：由步骤2得到过程扰动信号序列[e^o(k) e^o(k-1) ... e^o(k-p)]和步骤3得到的模型预测误差序列[e(k) e(k-1) ... e(k-p)]，计算模型性能指标η_R；Step 4: Get the process disturbance signal sequence [e^o (k) e^o (k-1) ... e^o (kp)] from step 2 and the model prediction error sequence [e(k) e(k -1)...e(kp)], calculation model performance index η_R ;

步骤5：比较模型性能指标η_R与预先确定的模型性能指标门槛值η_R0，如果模型性能指标η_R高于该门槛值，则认为该输出变量的模型性能良好，否则令l＝1，转到步骤6；Step 5: Compare the model performance index η_R with the predetermined model performance index threshold value η_R0 , if the model performance index η_R is higher than the threshold value, it is considered that the model performance of the output variable is good, otherwise set l=1, go to Go to step 6;

步骤6：移除第l个输入变量u_l，计算相应的第二过程扰动信号 $e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}];$ Step 6: Remove the lth input variable u_l and calculate the corresponding second process disturbance signal $e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}];$

步骤7：移除第l个输入变量u_l，分别计算k至k-p时刻的第二模型预测误差，得到第二模型预测误差序列e_pl＝[e_l(k) e_l(k-1) ... e_l(k-p)]；Step 7: Remove the lth input variable u_l , and calculate the second model prediction error at time k to kp respectively, and obtain the second model prediction error sequence e_pl =[e_l (k) e_l (k-1) . .. e_l (kp)];

步骤8：由步骤6得到的第二过程扰动信号 $e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}]$ 和步骤7得到的第二模型预测误差序列[e_l(k) e_l(k-1) ... e_l(k-p)]，计算第二模型性能指标η_Rl以及性能指标比κ_l；Step 8: The second process disturbance signal obtained from step 6 $e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}]$ And the second model prediction error sequence [e_l (k) e_l (k-1) ... e_l (kp)] obtained in step 7, calculate the second model performance index η_Rl and performance index ratio κ_l ;

步骤9：根据性能指标比κ_l判断输入变量u_l(k)对应的子模型性能，如果κ_l>1，表示输入变量u_l(k)对应的子模型性能恶化，需要维护，反之如果κ_l<1，表示相应的子模型性能良好；Step 9: Judging the performance of the sub-model corresponding to the input variable u_l (k) according to the performance index ratio κ_l , if κ_l > 1, it means that the performance of the sub-model corresponding to the input variable u_l (k) deteriorates and needs to be maintained; otherwise, if κ_l < 1, indicating that the corresponding sub-model has good performance;

步骤10：令l＝l+1，返回步骤6，直至l＝N_u，所有的输入处理完毕；Step 10: set l=l+1, return to step 6 until l=N_u , all input processing is completed;

步骤11：返回步骤1处理下一个输出变量y(k)，直至所有的输出变量评估完毕。Step 11: Return to step 1 to process the next output variable y(k), until all output variables are evaluated.

进一步，在所述步骤1中，构造k时刻时间窗p内向量过程如下：Further, in the step 1, the process of constructing the vector in the time window p at time k is as follows:

y_p(k)＝[y(k) y(k-1) ... y(k-p)] y_p (k)=[y(k) y(k-1) ... y(kp)]

$e_{p}^{o} (k) = [\begin{matrix} e^{o} (k) & e^{o} (k - 1) & . . . & e^{o} (k - p) \end{matrix}],$ 其中j＝1,...N_u $e_{p}^{o} (k) = [\begin{matrix} e^{o} (k) & e^{o} (k - 1) & . . . & e^{o} (k - p) \end{matrix}],$ where j=1,...N_u

u_jp(k)＝[u_j(k) u_j(k-1) ... u_j(k-p)] u_jp (k)＝[u_j (k) u_j (k-1) ... u_j (kp)]

${Y Y}_{m m} ((k k - - 11)) = = [\begin{matrix} {y the y}_{p p} ((k k - - 11)) \\ {y the y}_{p p} ((k k - - 22)) \\ . . \\ . . \\ . . \\ {y the y}_{p p} ((k k - - m m)) \end{matrix}],, {U u}_{n no} ((k k - - 11)) = = [\begin{matrix} {u u}_{11 p p} ((k k - - 11)) \\ . . \\ . . \\ . . \\ {u u}_{11 p p} ((k k - - n no)) \\ . . \\ . . \\ . . \\ {u u}_{{N N}_{u u} p p} ((k k - - 11)) \\ . . \\ . . \\ . . \\ {u u}_{{N N}_{u u} p p ((k k - - n no))} \end{matrix}],, {\overset{&OverBar; &OverBar;}{Z Z}}_{p p} ((k k)) = = [\begin{matrix} {Y Y}_{m m} ((k k - - 11)) \\ {U u}_{n no} ((k k - - 11)) \end{matrix}] - - - - - - ((11))$

其中，Y_m(k-1)由k-1时刻至k-m时刻的y_p(k)构造，U_n(k-1)由k-1时刻至k-n时刻的u_jp(k)构造，m、n分别输出变量和输入变量的阶次。Among them, Y_m (k-1) is constructed from y_p (k) from time k-1 to time km, U_n (k-1) is constructed from u_jp (k) from time k-1 to time kn, m, n is the order of the output variable and the input variable, respectively.

进一步，在所述步骤2中，过程扰动信号向量计算过程如下：Further, in the step 2, the process disturbance signal vector The calculation process is as follows:

定义行空间的正交补投影：definition Orthogonal complement projection to row space:

${Π Π}_{{\overset{&OverBar; &OverBar;}{Z Z}}_{p p} ((k k))}^{&perp; &perp;} = = I I - - {\overset{&OverBar; &OverBar;}{Z Z}}_{p p} {((k k))}^{T T} {[[{\overset{&OverBar; &OverBar;}{Z Z}}_{p p} ((k k)) {\overset{&OverBar; &OverBar;}{Z Z}}_{p p} {((k k))}^{T T}]]}^{- - 11} {\overset{&OverBar; &OverBar;}{Z Z}}_{p p} ((k k)) - - - - - - ((22))$

其中I为(m+n*N_u)维的单位矩阵。Wherein I is a (m+n*N_u )-dimensional identity matrix.

当p足够大时，通过正交补投影计算过程扰动信号向量When p is large enough, the process perturbation signal vector is calculated by orthogonal complementary projection

${e e}_{p p}^{o o} ((k k)) = = {y the y}_{p p} ((k k)) {Π Π}_{{\overset{&OverBar; &OverBar;}{z z}}_{p p}}^{&perp; &perp;} - - - - - - ((33))$

为简化计算，对 ${[\begin{matrix} {\overset{&OverBar;}{Z}}_{p} (k) & y_{p} (k) \end{matrix}]}^{T}$ 实施LQ分解To simplify the calculation, the ${[\begin{matrix} {\overset{&OverBar;}{Z}}_{p} (k) & {the y}_{p} (k) \end{matrix}]}^{T}$ Implementing LQ Decomposition

$[\begin{matrix} {\overset{&OverBar; &OverBar;}{Z Z}}_{p p} ((k k)) \\ {y the y}_{p p} ((k k)) \end{matrix}] = = [\begin{matrix} {L L}_{1111} \\ {L L}_{21 twenty one} & {L L}_{22 twenty two} \end{matrix}] [\begin{matrix} {Q Q}_{11} \\ {Q Q}_{22} \end{matrix}] - - - - - - ((44))$

其中向量Q₁，Q₂正交， $[\begin{matrix} L_{11} \\ L_{21} & L_{22} \end{matrix}]$ 为LQ分解后得到的下三角阵。where the vectors Q₁ , Q₂ are orthogonal, $[\begin{matrix} L_{11} \\ L_{twenty one} & L_{twenty two} \end{matrix}]$ It is the lower triangular matrix obtained after LQ decomposition.

然后计算过程扰动信号向量Then calculate the process disturbance signal vector

其中符号其中表示广义逆运算。where symbol where Represents the generalized inverse operation.

进一步，在所述步骤3中，计算k时刻模型预测误差e(k)的过程如下： Further, in said step 3, the process of calculating the model prediction error e(k) at time k is as follows:

计算模型预测输出y_m(k)：Calculate the model predicted output y_m (k):

${y the y}_{m m} ((k k)) = = {Σ Σ}_{j j = = 11}^{{N N}_{u u}} {Σ Σ}_{i i = = 11}^{{N N}_{00} - - 11} {a a}_{ji the ji} {Δu Δu}_{j j} ((k k - - i i)) + + {Σ Σ}_{j j = = 11}^{{N N}_{u u}} {a a}_{{jN n}_{00} {u u}_{j j} ((k k - - {N N}_{00})) - - - - - - ((66))}$

其中N₀为预测模型长度,N_u为对应于输出变量y的输入变量数,a_ji(i＝1,...N_u,j＝1,...,N₀)为预测控制器设计阶段确定的针对第j个输入变量的第i个阶跃响应系数，Δu_j(k-i)为第j个输入变量两个相邻采样时刻间的变化量Δu_j(k-i)＝u_j(k-i)-u_j(k-i-1)；where N₀ is the length of the predictive model, N_u is the number of input variables corresponding to the output variable y, a_ji (i=1,...N_u ,j=1,...,N₀ ) is the predictive controller design The i-th step response coefficient for the j-th input variable determined in the stage, Δu_j (ki) is the variation of the j-th input variable between two adjacent sampling moments Δu_j (ki)=u_j (ki) -u_j (ki-1);

根据输出变量y(k)、模型预测输出y_m(k)计算模型预测误差Calculate the model prediction error based on the output variable y(k) and the model prediction output y_m (k)

e(k)＝(1-q^-1)(y(k)-y_m(k)) (7)e(k)=(1-q^-1 )(y(k)-y_m (k)) (7)

其中q^-1为一步延迟算子。Where q^-1 is a one-step delay operator.

进一步，在所述步骤4中，模型性能指标η_R计算过程如下：Further, in said step 4, the model performance index η_R calculation process is as follows:

${η η}_{R R} = = \frac{{Σ Σ}_{i i = = 00}^{p p} {[[{e e}^{o o} ((k k - - i i))]]}^{22}}{{Σ Σ}_{i i = = 00}^{p p} {[[e e ((k k - - i i))]]}^{22}} - - - - - - ((88))$

进一步，步骤6中计算移除输入变量u_l后相应的第二过程扰动信号过程如下：Further, in step 6, calculate the corresponding second process disturbance signal after removing the input variable u_l The process is as follows:

移除输入变量u_l，类似式(1)重新构造时间窗内向量：Remove the input variable u_l and reconstruct the vector in the time window similar to formula (1):

y_p(k)＝[y(k) y(k-1) ... y(k-p)] y_p (k)=[y(k) y(k-1) ... y(kp)]

$e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}],$ 其中j＝1,...,l-1,l+1,...N_u $e_{pl}^{o} (k) = [\begin{matrix} e_{l}^{o} (k) & e_{l}^{o} (k - 1) & . . . & e_{l}^{o} (k - p) \end{matrix}],$ where j=1,...,l-1,l+1,...N_u

${Y Y}_{m m} ((k k - - 11)) = = [\begin{matrix} {y the y}_{p p} ((k k - - 11)) \\ {y the y}_{p p} ((k k - - 22)) \\ . . \\ . . \\ . . \\ {y the y}_{p p} ((k k - - m m)) \end{matrix}],,$

${U u}_{nl nl} ((k k - - 11)) = = {[\begin{matrix} {u u}_{11 p p} ((k k - - 11)) & . . . . . . & {u u}_{11 p p} ((k k - - n no)) & . . . . . . & {u u}_{((l l - - 11)) p p} ((k k - - 11)) & . . . . . . & {u u}_{((l l - - 11)) p p} ((k k - - n no)) & {u u}_{((l l + + 11)) p p} ((k k - - 11)) & . . . . . . & {u u}_{((l l + + 11)) p p} ((k k + + 11)) & . . . . . . & {u u}_{{N N}_{u u} p p} ((k k - - 11)) & . . . . . . & {u u}_{{N N}_{u u} p p} ((k k - - n no)) \end{matrix}]}^{T T},,$

${\overset{&OverBar; &OverBar;}{Z Z}}_{pl pl} ((k k)) = = [\begin{matrix} {Y Y}_{m m} ((k k - - 11)) \\ {U u}_{nl nl} ((k k - - 11)) \end{matrix}] - - - - - - ((99))$

其中T符号表示转置。然后实施LQ分解where T symbol means transpose. Then implement the LQ decomposition

$[\begin{matrix} {\overset{&OverBar; &OverBar;}{Z Z}}_{pl pl} ((k k)) \\ {y the y}_{p p} ((k k)) \end{matrix}] = = [\begin{matrix} {L L}_{1111}^{l l} \\ {L L}_{21 twenty one}^{l l} & {L L}_{22 twenty two}^{l l} \end{matrix}] [\begin{matrix} {Q Q}_{11}^{l l} \\ {Q Q}_{22}^{l l} \end{matrix}] - - - - - - ((1010))$

其中向量与正交， $[\begin{matrix} L_{11}^{l} \\ L_{21}^{l} & L_{22}^{l} \end{matrix}]$ 为LQ分解后得到的下三角阵。where the vector and Orthogonal, $[\begin{matrix} L_{11}^{l} \\ L_{twenty one}^{l} & L_{twenty two}^{l} \end{matrix}]$ It is the lower triangular matrix obtained after LQ decomposition.

计算移除输入变量ul后的过程扰动信号Calculate the process disturbance signal after removing the input variable ul

进一步，在所述步骤7中,移除输入变量u_l后，计算k时刻的模型预测误差e_l(k)过程如下：Further, in the step 7, after removing the input variable u_l , the process of calculating the model prediction error e_l (k) at time k is as follows:

移除输入变量u_l，式(6)第二项中去除a_liΔu_l(k-i)项计算y_ml(k)Remove the input variable u_l , remove the item a_li Δu_l (ki) from the second item in formula (6) to calculate y_ml (k)

${y the y}_{ml ml} ((k k)) = = {Σ Σ}_{j j = = 11}^{l l - - 11} {Σ Σ}_{i i = = 11}^{{N N}_{00} - - 11} {a a}_{ji the ji} {Δu Δu}_{j j} ((k k - - i i)) + + {Σ Σ}_{j j = = l l + + 11}^{{N N}_{u u}} {Σ Σ}_{i i = = 11}^{{N N}_{00} - - 11} {a a}_{ji the ji} {Δu Δ u}_{j j} ((k k - - i i)) + + {Σ Σ}_{j j = = 11}^{l l - - 11} {a a}_{j j {N N}_{00}} {u u}_{j j} ((k k - - {N N}_{00})) + + {Σ Σ}_{j j = = l l + + 11}^{{N N}_{00}} {a a}_{{jN n}_{00}} {u u}_{j j} ((k k - - {N N}_{00})) - - - - - - ((1212))$

按式(13)计算第二模型预测误差e_l(k)：Calculate the prediction error e_l (k) of the second model according to formula (13):

e_l(k)＝(1-q^-1)(y(k)-y_ml(k)) (13) e_l (k)=(1-q^-1 )(y(k)-y_ml (k)) (13)

进一步，在所述步骤8中，第二模型性能指标η_Rl以及性能指标比κ_l计算过程如下：Further, in said step 8, the calculation process of the second model performance index η_R1 and performance index ratio κ₁ is as follows:

${η η}_{Rl Rl} = = \frac{{Σ Σ}_{i i = = 00}^{p p} {(({e e}_{l l}^{o o} ((k k - - i i))))}^{22}}{{Σ Σ}_{i i = = 00}^{p p} {(({e e}_{l l} ((k k - - i i))))}^{22}} - - - - - - ((1414))$

${κ κ}_{l l} = = \frac{{η η}_{Rl Rl}}{{η η}_{R R}} - - - - - - ((1515))$

本发明所达到的有益效果：本发明为工程师们提供了一种预测控制模型深度监控的方法，该方法无需中断生产过程，也不需要工艺先验知识，只利用生产过程数据对预测控制整体模型和各子模型的性能进行评估，从而判断对象哪部分模型性能出现了下降，以及下降到了什么程度，使工程师们在预测控制系统维护前预知模型是否需要维护，以及哪部分需要维护，节省了不必要的再建模成本。Beneficial effects achieved by the present invention: the present invention provides engineers with a method for in-depth monitoring of the predictive control model. The method does not need to interrupt the production process and does not require prior knowledge of the process, and only uses the production process data to predict the overall model of control. And the performance of each sub-model is evaluated to judge which part of the model performance of the object has declined and to what extent, so that engineers can predict whether the model needs to be maintained and which part needs to be maintained before predicting the control system maintenance, saving a lot of money. Necessary remodeling costs.

附图说明Description of drawings

图1是模型性能指标方法原理图；Figure 1 is a schematic diagram of the model performance index method;

图2是深度诊断子模型性能方法示意图；Fig. 2 is a schematic diagram of the performance method of the deep diagnosis sub-model;

图3是Wood-Berry蒸馏塔实验结果图。Figure 3 is a diagram of the experimental results of the Wood-Berry distillation column.

具体实施方式Detailed ways

下面结合附图和具体实施方式对本发明进行详细说明。The present invention will be described in detail below in conjunction with the accompanying drawings and specific embodiments.

本发明所采用的技术方案原理如图1所示。The principle of the technical solution adopted by the present invention is shown in Fig. 1 .

考虑图1所示的预测控制系统，其中G^o(q)和H^o(q)分别表示过程模型和过程干扰模型，G_m(q)和H(q)分别表示预测控制模型和预测干扰模型，G_c(q)是预测控制器，e^o(k)和d^o(k)是过程干扰信号和扰动量，e(k)和d(k)是预测误差和控制模型误差，r(k)为参考轨迹，u(k)为控制输入(操作变量)，y(k)为输出变量(被控变量)，为预测输出。Consider the predictive control system shown in Fig. 1, where G^o (q) and H^o (q) denote the process model and process disturbance model respectively, G_m (q) and H(q) denote the predictive control model and predictive disturbance model , G_c (q) is a predictive controller, e^o (k) and d^o (k) are process disturbance signals and disturbances, e(k) and d(k) are prediction errors and control model errors, r(k ) is the reference trajectory, u(k) is the control input (operated variable), y(k) is the output variable (controlled variable), for the predicted output.

假定控制模型和扰动模型均匹配，显然，预测误差e(k)应等于过程扰动信号e^o(k)。反之，如果存在模型失配，由于预测误差e(k)包含模型失配信息，应大于过程扰动信号e^o(k)，因此定义模型性能指标Assuming that the control model and the disturbance model are matched, obviously, the prediction error e(k) should be equal to the process disturbance signal e^o (k). Conversely, if there is a model mismatch, since the prediction error e(k) contains model mismatch information, it should be greater than the process disturbance signal e^o (k), so define the model performance index

${η η}_{R R} = = \frac{{Σ Σ}_{k k = = 11}^{N N} {(({e e}^{o o} ((k k))))}^{22}}{{Σ Σ}_{k k = = 11}^{N N} {((e e ((k k))))}^{22}}$

其中N为评估时域数据长度，e^o(k)可通过输出变量y(k)和输入变量u(k)的常规运行数据得到，而e(k)可由模型阶跃响应系数和运行数据得到。η_R的范围为(0,1]，如果η_R接近于1表明模型性能良好，η_R接近于零则表明模型性能恶化。where N is the length of the evaluation time-domain data, e^o (k) can be obtained from the conventional operating data of the output variable y(k) and input variable u(k), and e(k) can be obtained from the model step response coefficient and operating data . The range of η_R is (0,1], if η_R is close to 1, it indicates that the model performance is good, and if η_R is close to zero, it indicates that the model performance is deteriorated.

性能指标η_R能对整体模型性能给出评估，但对于多输入多输出过程，预测模型是由多个单输入单输出的子模型构成的，模型整体性能下降并不一定意味着所有的子模型性能都下降。因此，对于性能下降的模型，本发明提出了进一步诊断各子模型性能的方法，如图2所示。The performance index η_R can evaluate the overall model performance, but for the multi-input multi-output process, the prediction model is composed of multiple single-input and single-output sub-models, and the decline in the overall performance of the model does not necessarily mean that all sub-models Performance drops. Therefore, for models with degraded performance, the present invention proposes a method for further diagnosing the performance of each sub-model, as shown in FIG. 2 .

假定判断针对第l个输入u_l的子模型性能，由于难以从输出变量y(k)中获取对应于u_l的分量，因此针对u_l的子模型性能无法通过η_R的计算式直接计算得到。本发明提出的方法为，逐一去除输入变量u_l后计算新的模型性能指标，如果去除后模型性能指标变好，说明移除的输入变量对应的子模型性能较差，需要进行维护。Assuming that the performance of the sub-model for the l-th input u_l is judged, since it is difficult to obtain the component corresponding to u_l from the output variable y(k), the performance of the sub-model for u_l cannot be directly calculated by the formula η_R . The method proposed by the present invention is to calculate new model performance indicators after removing the input variables u_l one by one. If the model performance indicators become better after removal, it means that the sub-models corresponding to the removed input variables have poor performance and need to be maintained.

由于输入变量u_l对输出的作用无法从输出变量y(k)中去除，因此把u_l的作用从输入通道转移到扰动通道，这样所有输入通道和扰动通道对输出的总作用不变，y(k)仍可以用来计算e^o(k)，这样计算得到的e^o(k)包含所移除的u_l的作用。计算e(k)时也从预测模型中去除u_l相关项，计算结果e(k)中也包含所去除的第l个子模型的作用。Since the effect of the input variable u_l on the output cannot be removed from the output variable y(k), the effect of u_l is transferred from the input channel to the perturbation channel, so that the total effect of all input channels and perturbation channels on the output remains unchanged, y (k) can still be used to calculate e^o (k), so that the calculated e^o (k) contains the effect of the removed u_l . When calculating e(k), u_l related items are also removed from the prediction model, and the calculated result e(k) also includes the effect of the removed lth sub-model.

具体实现步骤如下：The specific implementation steps are as follows:

y_p(k)＝[y(k) y(k-1) ... y(k-p)] y_p (k)=[y(k) y(k-1) ... y(kp)]

其中，Y_m(k-1)由k-1时刻至k-m时刻的y_p(k)构造，U_n(k-1)由k-1时刻至k-n时刻的u_jp(k)构造。Among them, Y_m (k-1) is constructed from y_p (k) from time k-1 to time km, and U_n (k-1) is constructed from u_jp (k) from time k-1 to time kn.

计算模型预测输出y_m(k)：Calculate the model predicted output y_m (k):

${y the y}_{m m} ((k k)) = = {Σ Σ}_{j j = = 11}^{{N N}_{u u}} {Σ Σ}_{i i = = 11}^{{N N}_{00} - - 11} {a a}_{ji the ji} {Δu Δ u}_{j j} ((k k - - i i)) + + {Σ Σ}_{j j = = 11}^{{N N}_{u u}} {a a}_{{jN n}_{00} {u u}_{j j} ((k k - - {N N}_{00})) - - - - - - ((66))}$

e(k)＝(1-q^-1)(y(k)-y_m(k)) (7)e(k)=(1-q^-1 )(y(k)-y_m (k)) (7)

其中q^-1为一步延迟算子。Where q^-1 is a one-step delay operator.

进一步，步骤6中计算移除输入变量u_l后相应的第二过程扰动信号过程如下： Further, in step 6, calculate the corresponding second process disturbance signal after removing the input variable u_l The process is as follows:

y_p(k)＝[y(k) y(k-1) ... y(k-p)] y_p (k)=[y(k) y(k-1) ... y(kp)]

计算移除输入变量u_l后的过程扰动信号Calculate the process disturbance signal after removing the input variable u_l

e_l(k)＝(1-q^-1)(y(k)-y_ml(k)) (13) e_l (k)=(1-q^-1 )(y(k)-y_ml (k)) (13)

${κ κ}_{l l} = = \frac{{η η}_{Rl Rl}}{{η η}_{R R}} - - - - - - ((1515))$

实验仿真与分析：Experimental simulation and analysis:

本发明提出的预测模型深度诊断方法在Wood-Berry蒸馏塔进行了仿真实验。该过程的传递函数矩阵G(s)来自参考文献R.K.Wood,M.W.Berry,Terminal composition control of a binary distillation column,Chemical Engineering Science28(1973)1707–1717中。The prediction model depth diagnosis method proposed by the present invention has been simulated in a Wood-Berry distillation tower. The transfer function matrix G(s) of this process comes from the references R.K.Wood, M.W.Berry, Terminal composition control of a binary distillation column, Chemical Engineering Science 28(1973) 1707-1717.

$G G ((s the s)) = = [\begin{matrix} \frac{{12.8 12.8 e e}^{- - s the s}}{16.7 16.7 s the s + + 11} & \frac{{- - 18.9 18.9 e e}^{- - 33 s the s}}{21.0 21.0 s the s + + 11} \\ \frac{{6.6 6.6 e e}^{- - 77 s the s}}{10.9 10.9 s the s + + 11} & \frac{{- - 19.4 19.4 e e}^{- - 33 s the s}}{1414 . . s the s + + 11} \end{matrix}] . . - - - - - - ((1616))$

回流量与蒸汽流量为该过程的输入变量(操作变量)，分别记为u₁和u₂，单位lb/min，塔顶和塔底产品组分为两个被控变量，分别记为y₁和y₂，单位mol％。假定采样周期为1分钟，离散化后过程传递函数矩阵G^o(q)为The reflux flow rate and the steam flow rate are the input variables (operating variables) of the process, which are denoted as u₁ and u₂ respectively, and the unit is lb/min. The product components at the top and bottom of the tower are two controlled variables, denoted as y₁ respectively and y₂ are in mol%. Assuming that the sampling period is 1 minute, the process transfer function matrix G^o (q) after discretization is

${G G}^{o o} ((q q)) = = [\begin{matrix} {q q}^{- - 22} \frac{0.744 0.744}{11 - - {0.9419 0.9419 q q}^{- - 11}} & {q q}^{- - 44} \frac{- - 0.8789 0.8789}{11 - - {0.9535 0.9535 q q}^{- - 11}} \\ {q q}^{- - 88} \frac{0.5786 0.5786}{11 - - {0.9123 0.9123 q q}^{- - 11}} & {q q}^{- - 44} \frac{- - 1.302 1.302}{11 - - {0.9329 0.9329 q q}^{- - 11}} \end{matrix}] - - - - - - ((1717))$

假定实际干扰过程模型H^o(k)为Suppose the actual disturbance process model H^o (k) is

${H h}^{o o} ((k k)) = = [\begin{matrix} \frac{11 - - {0.5 0.5 q q}^{- - 11}}{11 - - {q q}^{- - 11}} \\ \frac{11 - - {0.7 0.7 q q}^{- - 11}}{11 - - {q q}^{- - 11}} \end{matrix}] - - - - - - ((1818))$

施加的扰动信号e^o(k)是方差为diag{0.744²,0.13²}的独立白噪声。The applied perturbation signal e^o (k) is independent white noise with variance diag{0.744² ,0.13² }.

MPC预测时域和控制时域分别选为100和10。权矩阵选为Q＝diag{1,10}，S＝diag{1,10}。两个输出变量设定值分别为The MPC prediction time domain and control time domain are selected as 100 and 10, respectively. The weight matrix is selected as Q=diag{1,10}, S=diag{1,10}. The set values of the two output variables are

为了验证所研究方法的有效性，采用不同对象模型和扰动模型，按照前述方法步骤进行了试验，模型输出阶次m和输入阶次n均取30，时间窗长度p取3940。In order to verify the effectiveness of the researched method, different object models and disturbance models were used to carry out experiments according to the steps of the previous method. The model output order m and input order n were both set to 30, and the time window length p was set to 3940.

情况1：精确的对象模型和扰动模型；Case 1: exact object model and perturbation model;

情况2：精确的对象模型，扰动模型为DMC常用的模型形式Case 2: Accurate object model, disturbance model is a commonly used model form of DMC

$H h ((k k)) = = [\begin{matrix} \frac{11}{11 - - {q q}^{- - 11}} \\ \frac{11}{11 - - {q q}^{- - 11}} \end{matrix}] - - - - - - ((1919))$

情况3：对象模型和扰动模型均失配，过程模型选为Case 3: Both the object model and the disturbance model are mismatched, and the process model is chosen as

${G G}^{o o} ((q q)) = = [\begin{matrix} 00 . . 88 * * {q q}^{- - 22} \frac{0.744 0.744}{11 - - {0.9419 0.9419 q q}^{- - 11}} & 0.6 0.6 * * {q q}^{- - 44} \frac{- - 0.8789 0.8789}{11 - - {0.9535 0.9535 q q}^{- - 11}} \\ {q q}^{- - 88} \frac{0.5786 0.5786}{11 - - {0.9123 0.9123 q q}^{- - 11}} & 1.2 1.2 * * {q q}^{- - 44} \frac{- - 1.302 1.302}{11 - - {0.9329 0.9329 q q}^{- - 11}} \end{matrix}] - - - - - - ((2020))$

扰动模型为(19)式；The disturbance model is (19);

情况4：对象模型和扰动模型均失配，过程模型选为Case 4: Both the object model and the disturbance model are mismatched, and the process model is chosen as

${G G}^{o o} ((q q)) = = [\begin{matrix} {q q}^{- - 22} \frac{0.744 0.744}{11 - - {0.9419 0.9419 q q}^{- - 11}} & {q q}^{- - 44} \frac{- - 0.8789 0.8789}{11 - - {0.9535 0.9535 q q}^{- - 11}} \\ {q q}^{- - 88} \frac{0.5786 0.5786}{11 - - {0.9123 0.9123 q q}^{- - 11}} & {3.0 3.0 * * q q}^{- - 44} \frac{- - 1.302 1.302}{11 - - {0.65 0.65 q q}^{- - 11}} \end{matrix}] - - - - - - ((21 twenty one))$

扰动模型为(19)式；The disturbance model is (19);

情况5：情况4中移除第一个输入变量u₁，诊断y₂模型；Case 5: remove the first input variable u₁ in case 4, and diagnose the y₂ model;

情况6：情况4中移除第二个输入变量u₂，诊断y₂模型；Case 6: remove the second input variable u₂ in case 4, and diagnose the y₂ model;

相应的MQI值见表1，相应的直方图比较见图2。The corresponding MQI values are shown in Table 1, and the corresponding histogram comparisons are shown in Figure 2.

表1 Wood-Berry蒸馏塔模型评估结果Table 1 Wood-Berry distillation column model evaluation results

从表1和图2可以看出： It can be seen from Table 1 and Figure 2 that:

(1)情况1中模型匹配时，两个被控变量的MQI指标均接近于1，表明MQI指标能够精确反映模型性能。(1) When the model is matched in case 1, the MQI indicators of the two controlled variables are close to 1, indicating that the MQI indicators can accurately reflect the model performance.

如果过程模型或扰动模型出现性能下降(如情况2,情况3，情况4)，MQI指标即出现下降，而且下降程序与模型不匹配程度相对应。如情况3和4的MQI2，情况3的增益为1.2倍实际对象增益，情况4的增益则为3倍实际对象增益，相应的MQI2情况3为0.6090，而情况4为0.1188，因此MQI的数值大小能够反映模型的不匹配程度。If the performance of the process model or the perturbation model declines (such as case 2, case 3, case 4), the MQI index will decline, and the decline procedure corresponds to the degree of model mismatch. For example, MQI2 in cases 3 and 4, the gain in case 3 is 1.2 times the actual object gain, and the gain in case 4 is 3 times the actual object gain, the corresponding MQI2 case 3 is 0.6090, and case 4 is 0.1188, so the value of MQI It can reflect the degree of mismatch of the model.

(2)情况4中第4个子模型增益为3倍于实际增益，因此模型严重不匹配。计算得到的MQI2为0.1188，远小于1，因此可判断第二个输出变量模型性能严重下降，这与实际情况相符。这属于本发明方法的第一个阶段，即判断整体模型性能，尚不能判断y₂模型性能下降是由u₁对应子模型还是u₂对应子模型引起的。(2) In case 4, the gain of the fourth sub-model is 3 times of the actual gain, so the model is seriously mismatched. The calculated MQI2 is 0.1188, which is far less than 1. Therefore, it can be judged that the performance of the second output variable model is seriously degraded, which is consistent with the actual situation. This belongs to the first stage of the method of the present invention, that is, to judge the performance of the overall model, and it is not yet possible to judge whether the decline in the performance of the_y2 model is caused by the corresponding submodel of_u1 or the corresponding submodel_{of u2} .

(3)情况5和情况6属于本发明的第二个阶段，即采用逐一移除法计算移除输入变量后的MQI指标，根据其指标提升程度来判断移除的变量对应的子模型性能。由于情况4中判定y₂对应的模型性能恶化，情况5和情况6分别计算去除u₁和去除u₂后的MQI指标。(3) Case 5 and Case 6 belong to the second stage of the present invention, that is, the MQI index after removing the input variables is calculated by one-by-one removal method, and the performance of the sub-model corresponding to the removed variable is judged according to the degree of improvement of the index. Since the performance of the model corresponding to the decision y₂ in case 4 deteriorates, the MQI indicators after removing u₁ and u₂ are calculated in case 5 and case 6, respectively.

与情况4相比，情况5中MQI2从0.1188下降到0.0664，相应的κ值0.5589<1，表明移除了一个性能良好的子模型。实际上，情况4中与第一个输入u₁相关的子模型均为精确子模型，因此，判断结果与实际情况相吻合。Compared with case 4, the MQI2 in case 5 decreased from 0.1188 to 0.0664, corresponding to a κ value of 0.5589 < 1, indicating that a well-performing submodel was removed. In fact, the sub-models related to the first input u₁ in case 4 are all exact sub-models, so the judgment result is consistent with the actual situation.

相反，与情况4相比，情况6中MQI2大幅上升，从0.118上升到0.3727，相应的κ值3.1372>1，表明移除了一个性能恶化的子模型。实际上，情况4中与第二个输入u₂相关的子模型确实存在严重不匹配，因此，判断结果与实际情况相吻合。On the contrary, compared with case 4, the MQI2 in case 6 rises sharply, from 0.118 to 0.3727, with a corresponding κ value of 3.1372 > 1, indicating that a submodel with degraded performance is removed. In fact, there is indeed a serious mismatch in the sub-models related to the second input_u2 in case 4, so the judgment result is consistent with the actual situation.

以上所述仅是对本发明的较佳实施例而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所做的任何简单修改，等同变化与修饰，均属于本发明技术方案的范围内。The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention in any form. Any simple modifications made to the above embodiments according to the technical essence of the present invention, equivalent changes and modifications, all belong to this invention. within the scope of the technical solution of the invention.