A data-driven LQG method for linear time-varying systems to monitor the controller performance in the batch process

Abstract

Purpose

The flexibility of batch process enables its wide application in fine-chemical, pharmaceutical and semi-conductor industries, whilst its complexity necessitates control performance monitoring to ensure high operation efficiency. This paper proposes a data-driven approach to carry out controller performance monitoring within batch based on linear quadratic Gaussian (LQG) method.

Design/methodology/approach

A linear time-varying LQG method is proposed to obtain the joint covariance benchmark for the stochastic part of batch process input/output. From historical golden operation batch, linear time-varying (LTV) system and noise models are identified based on generalized observer Markov parameters realization.

Findings

Open/closed loop input and output data are applied to identify the process model as well as the disturbance model, both in Markov parameter form. Then the optimal covariance of joint input and output can be obtained by the LQG method. The Hotelling's Tˆ2 control chart can be established to monitor the controller.

Originality/value

(1) An observer Markov parameter approach to identify the time-varying process and noise models from both open and closed loop data, (2) a linear time-varying LQG optimal control law to obtain the optimal benchmark covariance of joint input and output and (3) a joint input and output multivariate control chart based on Hotelling's T² statistic for controller performance monitoring.

Citation

Chen, M. and Xie, L. (2023), "A data-driven LQG method for linear time-varying systems to monitor the controller performance in the batch process", Journal of Intelligent Manufacturing and Special Equipment, Vol. 4 No. 1, pp. 24-46. https://doi.org/10.1108/JIMSE-09-2022-0016

1. Introduction

Batch chemical processes are widely used in the production of fine chemicals, pharmaceutical products, semiconductor, polymers and many other materials (Caccavale et al., 2010). In contrast to continuous process, the operation of batch process often contains multiple stages and cannot be regarded as a linear time invariant (LTI) case. In addition, batch chemical reaction mechanism is typically complex (Aumi et al., 2013); the control of its dynamic process and thermal conduction speed is characterized by nonlinear dynamics, whose effects are further emphasized by intrinsically unsteady operating conditions (Mhaskar et al., 2019).

The complex nonlinear dynamics of batch process are often difficult to obtain and can be approximated with a moderate set of local LTI models (Garg, 2018). Each model corresponds to a characteristic operation region and is usually described by a set of constraints. To operate a batch process in an optimal fashion, such local approximations are not always sufficient since the transitions between them are missing to provide a complete description. Furthermore, even if one specific set of constraints holds true for a long period, local LTI models still cannot describe the time variation due to change of hold-ups or compositions (Yin et al., 2013). Therefore, multiple LTI models fail to describe the batch process adequately, and the corresponding control strategy often leads to poor performance. To tackle the above deficiencies, linear time varying (LTV) model (Xu et al., 2013) is employed to provide more accurate description of batch process with restricted complexity (Liu and Alleyne, 2015).

Due to the increasing competition in chemical industry, control performance assessment and monitoring (CPA/CPM) has received much attention to ensure high efficiency of process operations. Several techniques using minimum variance control (MVC) (Harris et al., 1996; Harris, 1989; Fu et al., 2012; Ko and Edgar, 2004) have been proven useful in prioritizing the LTI control loop maintenance activities since it only requires process time delay and routine operating data to evaluate controlled variable variance. To perform CPA/CPM of batch process which consists a sequence of operations impacting the final product yield and quality, Xu and Huang (2006) extended the MVC benchmark to the LTV process based on the assumption that the disturbance is time varying and piece-wise. For the case of both LTV process and disturbance models, recent studies (Huang, 2002; Wang et al., 2017; Li and Evans, 1997) supposed the LTV system satisfies the condition of pseudo-commutation, then established LTV-MVC benchmark and LTV-GMVC (generalized minimum variance control) benchmark, respectively. On the other hand, the batch process can be characterized by its frequent repetition of batch runs, and a high degree of reproducibility is necessary to obtain successful batch production. Taking this into account, Chen and Kong (2009) developed an optimal iterative learning control (ILC) algorithm to establish the MVC tracking and regulating achievable benchmark bounds. More recently, Farasat and Huang (2013) elaborate the optimal MVC solution to make a trade-off between stochastic and deterministic CPA/CPM (Ge, 2017) of batch processes.

In real practice, many types of restrictions, for example, process input/output constraints, render the MVC benchmark too ideal to achieve, and it may give rise to conservative performance assessment conclusions. One more practical technique, data-driven CPA/CPM, determines satisfactory control performance from the historical golden operation period (Yu and Joe Qin, 2008). Control performance index can be identified by comparing the statistical property of real-time data with that of benchmark data. Another deficiency of MVC benchmark relates to its preference in the variation of system output. More precisely, MVC does not take the input variability (such as valve movement) into account, which may accelerate the equipment wear. In contrast, linear quadratic Gaussian (LQG) benchmark is more appropriate to balance the output/input variance. Kadali and Huang (2002) proposed a subspace matrix-based method to obtain the LQG benchmark variances from closed loop data for the LTI process. But for the batch process, little has been done regarding making a trade-off between both variance of input and output.

In this paper, a linear time-varying LQG method is proposed to obtain the joint covariance benchmark for the stochastic part of batch process input/output. From historical golden operation batch, the LTV system and noise models are identified based on generalized observer Markov parameters realization (Majji et al., 2010a; Majji and Junkinst, 2008). In such models, the process output at each step depends on the linear convolution of the impulse response and the inputs applied until that instant. Based on the LTV models and benchmark covariance of joint input and output, the control charts based on Hotelling's T² statistics are built up to carry out CPA/CPM for batch process.

The contribution of the paper is composed of three parts. The first one is the initial identification of LTV process for either open loops or closed systems. The Markov parameters of state space are obtained from operation dataset which performed well, and the parameters are necessary for next calculation for LQG benchmark. The second is to compute the covariance of input and output variables. It is difficult to capture the real operating models under data driven method in LTV loops. If the system is time-varying, the time-invariant covariance benchmark might be not proper for monitoring. Hence, LQG optimal control law using identification to combine both models and dataset together is adopted to complete this task. The parameters to calculate LQG benchmark are from the first identification. After that, the covariance is used to get the Hotelling's T² statistic for further assessment. And the third one is to adopt the dynamic Hotelling's T² statistic to the process monitoring.

The rest of this article is arranged as follows: Section 2 derives the LQG method for the LTV batch process and establishes multivariate control charts based on Hotelling's T² index. Section 3 and 4 present the approach to obtain the system model Markov parameter matrices required for LQG analysis from historical open loop data and closed loop data, respectively. Section 5 summarizes the whole CPM procedure for batch process. Section 6 shows the application results from an industrial batch reactor process, which is followed by conclusions in Section 7.

2. LQG design and CPM for LTV system

2.1 LQG design for LTV system

Consider the following system with r inputs and m outputs represented by a LTV process with additive noise at the output in transfer function form:

(1)yt=∑i=1tHt,iui+∑i=1tNt,iai

where H_t,i and N_t,i are system Markov parameter matrices of the LTV process and the noise, respectively, a_i represents an uncorrelated random noise vector at discrete time i with zero-mean and covariance matrix Σ_a. Both y_i and u_i are expressed as deviation variables. We can first consider the influence of a single noise a₁ introduced at t = 1 and then apply the principle of superposition for the all noise case.

When a single random noise a₁ is introduced at t = 1, the whole tf step output evolution equations for the LTV batch process can be expressed as follows:

(2)y1y2y3⋮ytf−N1,1N2,1N3,1⋮Ntf,1a1=H1,100⋯0H2,1H2,20⋯0H3,1H3,2H3,3⋯0⋮⋮⋮⋱⋮Htf,1Htf,2Htf,3⋯Htf,tf.u1u2u3⋮utf

In order to match the state space form for convenient identification computations (details are in Section 3), some changes are introduced into the noise part in Equation (2). By defining

(3)e1=N1,1a1

(4)N¯i,1e1=Ni,1a1(i=2,3,⋯tf)

Equation (2) can be transformed into

(5)y1y2y3⋮ytf⏟y1|tf−ImN¯2,1N¯3,1⋮N¯tf,1⏟Le1|tf,1e1=H1,100⋯0H2,1H2,20⋯0H3,1H3,2H3,3⋯0⋮⋮⋮⋱⋮Htf,1Htf,2Htf,3⋯Htf,tf⏟Lu1|tf.u1u2u3⋮utf⏟u1|tf

and in compact form,

(6)y1|tf=Lu1|tfu1|tf+Le1|tf,1e1

with e₁ represents an uncorrelated zero-mean random noise vector with covariance matrix Σe=N1,1ΣaN1,1T∈Rm×m.

For batch process, the variable u_i, y_i can be separated into two parts, ui=uisp+uie, yi=yisp+yie, respectively. ytsp stands for the output deviation influenced by the deterministic setpoint while yie stands for that influenced by the stochastic noise. uisp is designed to track ytsp while uie is designed to counteract the effect of yie, so Equation (6) can be divided into two parts:

(7)y1|tfsp=Lu1|tfu1|tfsp

(8)y1|tfe=Lu1|tfu1|tfe+Le1|tf,1e1

From Equation (7), the ideal setpoint input part u^1|tfsp=(Lu1|tf)†y1|tfsp can be designed to counteract the effect of y1|tfsp.

In order to acquire optimal control law for the whole batch process, we consider the regulatory finite-horizon LQG control objective function for the stochastic noise effect in Equation (8) as follows:

(9)J=∑i=1tf(yi−yisp)Tyi−yisp+(ui−u^isp)T(λI)(ui−u^isp)=∑i=1tf[(yie)Tyie+(uie)T(λI)uie]=(y1|tfe)Ty1|tfe+(u1|tfe)T(λI)u1|tfe

where tf is the final instant of the whole batch and λ is the user-defined nonnegative input weighting parameter. Partial differentiation of J and considering the model in Equation (8) yield the LQG optimal control law:

(10)u1|tfopt=−(Lu1|tfTLu1|tf+λI)−1Lu1|tfTLe1|tf,1e1

Substituting it in Equation (8) gives the corresponding optimal output as follows:

(11)y1|tfopt=I−Lu1|tf(Lu1|tfTLu1|tf+λI)−1Lu1|tfTLe1|tf,1e1

Define the following two matrices:

(12)Φ1(1)Φ2(1)⋮Φtf(1)=−(Lu1|tfTLu1|tf+λI)−1Lu1|tfTLe1|tf,1

(13)Γ1(1)Γ2(1)⋮Γtf(1)=I−Lu1|tf(Lu1|tfTLu1|tf+λI)−1Lu1|tfTLe1|tf,1

The coefficient matrix Φ_i(j) (j ≤ i) can be regarded as the impulse response coefficient from the jth step noise e_j to the ith step optimal input, so as the coefficient Γ_i(j) for the output. When all the random from t = 1 to tf are taken into account, we have

(14)Φt(t)Φt+1(t)⋮Φtf(t)=−(Lut|tfTLut|tf+λI)−1Lut|tfTLet|tf,1(k=1,2,3,⋯tf)

Applying the principle of superposition obtains the optimal sequence of control inputs as follows:

(15)u1opt=Φ1(1)e1u2opt=Φ2(2)e2+Φ2(1)e1u3opt=Φ3(3)e3+Φ3(2)e2+Φ3(1)e1⋮utopt=∑i=1tΦt(i)ei⋮utfopt=Φtf(tf)etf+Φtf(tf−1)etf−1+⋯+Φtf(1)e1

Equation (15) describes the relationship between utopt and all the noise occurred before the time instant of t. Similarly, the sequence of the corresponding optimal outputs is

(16)ytopt=∑i=1tΓt(i)ei

From the above equations, it can be seen that the optimal input and output at instant t under LQG control are both related to the noise before instant t. From Equation (15) and (16), the optimal input and output for the stochastic noise effect part is determined by coefficients Φ_t(i) and Γ_t(i). Section 3 and 4 will elaborate how to obtain these coefficients from historical golden batch data.

2.2 Hotelling's T² CPM index for LTV system

Based on the optimal sequence of inputs and outputs, a benchmark can be set up to monitor the control performance of the whole batch. First, the covariance between the optimal input and output are calculated as follows:

(17)Cov(ujopt,ukopt)=∑i=1min(j,k)Φj(i)RΦkT(i)

(18)Cov(yjopt,ykopt)=∑i=1min(j,k)Γj(i)RΓkT(i)

(19)Cov(ujopt,ykopt)=∑i=1min(j,k)Φj(i)RΓkT(i)

Since the optimal inputs and outputs are dependent with each other, they are put together to formulate an augmented vector, stopt=((u1opt)T,(y1opt)T,⋯(utopt)T,(ytopt)T)T(t=1,2,…,tf) with covariance defined as follows:

(20)Cov(stopt,stopt)=Cov(u1opt,u1opt)⋯Cov(u1opt,ytopt)Cov(y1opt,u1opt)⋯Cov(y1opt,ytopt)Cov(u2opt,u1opt)⋯Cov(u2opt,ytopt)Cov(ytopt,u1opt)⋯Cov(ytopt,ytopt)∈Rt(r+m)×t(r+m)

Since the covariance of stopt may be poorly conditioned, its Singular Value Decomposition (SVD) decomposition is performed first to obtain a reliable monitoring statistic:

(21)Cov(stopt,stopt)=VtoptDtopt(Vtopt)T

where Dtopt∈Rdt×dt is a diagonal matrix storing the first d_t dominant eigenvalues of Cov(stopt,stopt) and the corresponding eigenvectors are in Vtopt∈Rt(r+m)×d.

When the a new batch data obtained in closed loop, the data sequences ste=(u1eT,y1eT,u2eT……uteT,yteT)T until time instant t are obtained by subtracting the deterministic part, the following Hotelling's T² CPM index can be established,

(22)Tt2=steTVtopt(Dtopt)−1(Vtopt)Tste

The index Tt2 quantifies the comparison between the new batch data and benchmark data. And the batch operation is considered to have good performance for a given significance level α if

(23)Tt2≤Tα2=(N2−1)dtN(N−dt)Fα(dt,N−dt)

where N is the total number of golden batches and F_α(d_t, N−d_t) is the critic value of the Fisher–Snedecor distribution with N and N−d_t degrees of freedom.

It should be noted that Hotelling's T² index for CPM in this paper is different from that often used for fault detection and diagnosis (FDD) in definition and sense. In definition, T² index for CPM in this paper comes from choosing all the components that have nonzero eigenvalues in the covariance of joint input and output because any covariance information on all input and output should not be left out, while T² index for FDD comes from choosing part of the principle components with the purpose of dimensionality reduction. In the sense, T² index for CPM in this paper reflects whether the real-time data are optimal under control or not, while T² index for FDD reflects whether the real-time data are normal or not, so T² index in CPM and in FDD have different meaning.

3. Computation of Markov parameter matrix from open loop data

This section elaborates how to compute the LTV process model Markov parameter matrix from open loop data. Unlike the LTI case, the identification of LTV model is not straightforward. For example, for the model in the state space form, the subspace identification regression analysis approach (Corbett and Mhaskar, 2016) directly from the Hankel matrix of process inputs and outputs for LTI case cannot be applied to the LTV case because of the variability of system parameters. On the other hand, although we can make full use of input–output data to identify the process model Markov parameter H_t,i in Equation (1) with classic least-square method, the noise model Markov parameter cannot be identified from input–output data directly.

An effective and simple method for LTV process model identification was developed by Majji (Majji et al., 2010b), using generalized observer Markov parameters which link the process model Markov parameter with noise model Markov parameter. In this section, the transfer function system model is transformed into state space form in order to obtain the generalized observer Markov parameters from open loop data. Then the relationship between the generalized observer/system Markov parameters is derived.

3.1 Input–output representation: observer Markov parameters

The system described by Equation (1) can be reformulated as

(24)xk+1=Akxk+Bkuk+Kkek

(25)yk=Ckxk+Dkuk+ek

with the state, output, input and noise dimensions xk∈Rn, yk∈Rm, uk∈Rr, ek∈Rm and the system matrices Ak∈Rn×n, Bk∈Rn×r, Ck∈Rm×n, Dk∈Rm×r, Kk∈Rn×m, respectively, where n is the state order, K_k is the Kalman filter gain at time k and e_k represents an uncorrelated zero-mean random noise vector with covariance matrix R∈Rm×m.

The whole tf step output evolution equations for the LTV batch process can be expressed as follows:

(26)y1y2y3⋮ytf=C1C2A1C3A2A1⋮CtfAtf−1⋯A1x1+D10⋯0C2B1D2⋯0C3A2B1C3B2⋯0⋮⋮⋯⋮CtfAtf−1⋯B1CtfAtf−1⋯B2⋯Dtf.u1u2u3⋮utf+Im0⋯0C2K1D2⋯0C3A2K1C3K2⋯0⋮⋮⋯⋮CtfAtf−1⋯K1CtfAtf−1⋯K2⋯Im.e1e2e3⋮etf

The transfer function form can be transformed into different state space form by choosing different state variable x_k, but the input–output relationship holds the same. By setting initial state x₁ = 0 and consider the same single noise case as Section 2, Equation (26) becomes

(27)y1y2y3⋮ytf−ImC2K1C3A2K1⋮CtfAtf−1⋯K1e1= D10⋯0C2B1D2⋯0C3A2B1C3B2⋯0⋮⋮⋯⋮CtfAtf−1⋯B1CtfAtf−1⋯B2⋯Dtf.u1u2u3⋮utf

From Equation (27) and (28), the relationship between the system parameter matrices A_k,B_k,C_k,D_k,K_k and the process model system Markov parameters H_i,k can be derived as

(28)Hi,k=CiAi−1⋯Ak+1Bk,∀i>k+1Ck+1Bk,i=k+1Dk,i=k0,∀i<k

And the noise model system Markov parameters as

(29)N¯i,k=CiAi−1⋯Ak+1Kk,∀i>k+1Ck+1Kk,i=k+1Im,i=k0,∀i<k

Further, the system described by Equation (24) and (25) can also be represented in the following forms:

(30)xk+1=A¯kxk+B¯kvk

(31)yk=Ckxk+Dkuk+ek

where

(32)A¯k=Ak−KkCkB¯k=Bk−KkDkKkvk=ukyk

Starting from the initial time step t = 1 and taking x₁ = 0, the input/output relation can be written as

(33)y1=D1u1+e1y2=D2u2+C2B¯1v1+e2y3=D3u3+C3B¯2v2+C3A¯2B¯1v1+e3⋮ytf=Dtfutf+∑j=1tf−1H¯tf,jvj+etf

where the process model observer Markov parameters H¯i,k is being defined as

(34)H¯i,k=CiA¯i−1⋯A¯k+1B¯k,∀i>k+1Ck+1B¯k,i=k+1Dk,i=k0,∀i<k

In order to identify H¯i,k, the datasets are collected from N(N > tf) batches, and the least-squares method is employed. Displaying Equation (33) in a vector matrix form for any time of k without loss of generality, we have

(35)yk=DkH¯k,k−1⋯H¯k,1ukvk−1vk−2⋮v1+ek

Repeating Equation (35) with N batches of recorded data yields

(36)Yk=yk(1)yk(2)yk(3)⋯yk(N)+Ek=MkVk+EkMk= DkH¯k,k−1H¯k,k−2⋯H¯k,1Vk= uk(1)vk−1(1)vk−2(1)⋮v1(1)uk(2)vk−1(2)vk−2(2)⋮v1(2)⋯uk(N)vk−1(N)vk−2(N)⋮v1(N)

The least squares solution for the process model observer Markov parameters and the estimated residual (innovation) e^k at time of k are given by

(37)M^k=YkVk†

(38)E^k=e^(1)ke^(2)ke^(3)k⋯e^(N)k=Yk−M^kVk

where (.)^† denotes the pseudo inverse of a matrix. So the noise covariance matrix R can be estimated as

(39)R^=Ee^ke^kT=CovE^kT

3.2 Computation of process model system Markov parameters from observer Markov parameters

A recursive method can be developed to calculate the process model system Markov parameters from the observer Markov parameters. From the definition of the observer Markov parameters, we have

(40)H¯k+1,k=Ck+1B¯k=Ck+1Bk−KkDkKk=H¯k+1,k(1)H¯k+1,k(2)

So

(41)H¯k+1,k(1)+H¯k+1,k(2)Dk=Ck+1Bk=Hk+1,k

A similar expression for Markov parameters of H¯k+2,k can be obtained as follows:

(42)H¯k+2,k=Ck+2A¯k+1B¯k=Ck+2A¯k+1Bk−KkDkKk=H¯k+2,k(1)H¯k+2,k(2)

(43)H¯k+2,k(1)+H¯k+2,k(2)Dk=Ck+2A¯k+1Bk=Ck+2(Ak+1−Kk+1Ck+1)Bk=Ck+2Ak+1Bk−H¯k+2,k+1(2)Ck+1Bk=Hk+2,k−H¯k+2,k+1(2)Hk+1,k

In general formulation, we have

(44)Hk+p,k(1)+H¯k+p,k(2)Dk=Ck+pA¯k+p−1A¯k+p−2⋯A¯k+1Bk=Ck+pA¯k+p−1⋯A¯k+2(Ak+1−Kk+1Ck+1)Bk=Ck+pA¯k+p−1⋯A¯k+2Ak+1Bk−H¯k+p,k+1(2)Hk+1,k=Ck+pA¯k+p−1⋯A¯k+3(Ak+2−Kk+2Ck+2)Ak+1Bk−H¯k+p,k+1(2)Hk+1,k=Ck+pA¯k+p−1⋯A¯k+3Ak+2Ak+1Bk−H¯k+p,k+2(2)Hk+2,k−H¯k+p,k+1(2)Hk+1,k⋯=Ck+pAk+p−1⋯Ak+1Bk−H¯k+p,k+p−1(2)Hk+p−1,k−H¯k+p,k+p−2(2)Hk+p−2,k−⋯=Hk+p,k−H¯k+p,k+p−1(2)H

Writing the derived recursive relationships between H_i,k and H¯i,k together yields a series of equations as follows:

(45)Hk+p,k+p−1=H¯k+p,k+p−1(1)+H¯k+p,k+p−1(2)Dk+p−1Hk+p,k+p−2−H¯k+p,k+p−1(2)Hk+p−1,k+p−2=H¯k+p,k+p−2(1)+H¯k+p,k+p−2(2)Dk+p−2⋮Hk+p,k−H¯k+p,k+p−1(2)Hk+p−1,k−H¯k+p,k+p−2(2)Hk+p−2,k+p−⋯−H¯k+p,k+1(2)Hk+1,k=H¯k+p,k(1)+H¯k+p,k(2)Dk

Denoting H̃i,j=H¯i,j(1)+H¯i,j(2)Dj, we obtain the matrix form description of linear equations relating the process model system/observer Markov parameters as follows:

(46)Im00⋱0−H¯k+2,k+1(2)Im0⋱0⋮⋮⋱⋱⋮−H¯k+p−1,k+1(2)−H¯k+p−1,k+2(2)⋯Im0−H¯k+p,k+1(2)−H¯k+p,k+1(2)⋯−H¯k+p,k+p−1(2)Im×Hk+1,k0⋯00Hk+2,kHk+2,k+1⋯00⋮⋮⋱⋮⋮Hk+p−1,kHk+p−1,k+1⋯Hk+p−1,k+p−20Hk+p,kHk+p,k+1⋯Hk+p,k+p−2Hk+p,k+p−1=H̃k+1,k0⋯00H̃k+2,kH̃k+2,k+1⋯00⋮⋮⋱⋮⋮H̃k+p−1,kH̃k+p−1,k+1⋯H̃k+p−1,k+p−20H̃k+p,kH̃k+p,k+1⋯H̃k+p,k+p−2H̃k+p,k+p−1

From Equation (46), it can be seen that the lower triangular component of Lu_1|tf in Equation (6) can be identified directly from the process model observer Markov parameters.

3.3 Computation of noise model system Markov parameters from observer Markov parameters

The last section describes the method to identify the process model. This section considers how to obtain noise model system Markov parameter matrix N¯i,k from observer Markov parameters H¯i,k.

Note the definition of Le_k|k + p,1 as follows:

(47)Lek|k+p,1=ImCk+1KkCk+2Ak+1Kk⋮Ck+pAk+p−1⋯Kk=ImOk+1(p)Kk

where

(48)Ok(p)=CkCk+1Ak⋮Ck+p−1Ak+p−2⋯Ak

Ok(p) is the system observability matrix, which can be obtained from the SVD decomposition of the system Markov parameter Hankel matrix Hk(p,q),

(49)Hk,k−1Hk,k−2⋯Hk,k−qHk+1,k−1Hk+1,k−2⋯Hk,k−q⋮⋮⋮⋮Hk+p−1,k−1Hk+p−1,k−2⋯Hk+p−1,k−q=Hk(p,q)=O^k(p)R^k−1(q)=CkCk+1AkCk+2Ak+1Ak⋮Ck+p−1Ak+p−2⋯AkBk−1Ak−1Bk−2Ak−1Ak−2Bk−3⋮Ak−1Ak−2⋯Bk−qT=UkΣk1/2Σk1/2Vk

where the parameters q are chosen such that the system Hankel matrix retains the state dimension n.

Equation (46) at consecutive time instants can be used to calculate the process model system Markov parameters, from which Hk(p,q) can be constructed. The observability matrix Ok(p) can be obtained from the SVD decomposition as shown in Equation (49). Then the parameter C_k can be identified from Ok(p) directly as follows:

(50)C^k=O^k(p)(1:m,:)

where the notation 1: m denotes the first m rows of the matrix. Then from Equation (40), the following equation holds

(51)K^k=C^k+1†H¯k+1,k(2)

And Le_k|k + p,1 can be estimated by

(52)Lek|k+p,1=ImO^k+1(p)K^k=ImO^k+1(p)O^k+1(p)(1:m,:)†H¯k+1,k(2)

4. Computation of Markov parameter matrix from closed loop data

Open loop operation of the process may not always be possible for safety reasons and other practical limitations. Estimation of the system model matrix from closed loop data is desirable in such cases.

Extend the process model Equation (5) to the case including all the noise effects as follows:

(53)y1y2y3⋮ytf⏟y1|tf=H1,100⋯0H2,1H2,20⋯0H3,1H3,2H3,3⋯0⋮⋮⋮⋱⋮Htf,1Htf,2Htf,3⋯Htf,tf⏟Lu1|tfu1u2u3⋮utf⏟u1|tf+Im00⋯0N¯2,1Im0⋯0N¯3,1N¯3,2Im⋯0⋮⋮⋮⋱⋮N¯tf,1N¯tf,2N¯tf,3⋯Im⏟Le1|tfe1e2e3⋮etf⏟e1|tf

or in compact form,

(54)y1|tf=Lu1|tfu1|tf+Le1|tfe1|tf

Consider the process system described by Equation (54) under closed loop with an LTI feedback controller L_C, which is expressed in an impulse response matrix form as follows:

(55)u1u2u3⋮utf⏟u1|tf=c100⋯0c2c10⋯0c3c2c1⋯0⋮⋮⋮⋱⋮ctfctf−1ctf−2⋯c1⏟LC.r1−y1r2−y2r3−y3⋮rtf−ytf⏟r1|tf−y1|tf

In compact form,

(56)u1|tf=LC(r1|tf−y1|tf)

where r_1|tf is the setpoint trajectory for the process output. Assume that the controller does not cancel any plant dynamics, substituting Equation (56) into Equation (54) yields

(57)y1|tf=Lyrr1|tf+Lyee1|tf

(58)u1|tf=Lurr1|tf+Luee1|tf

where

(59)Lyr=(I+Lu1|tfLC)−1Lu1|tfLCLye=(I+Lu1|tfLC)−1Le1|tfLur=LCI−(I+Lu1|tfLC)−1Lu1|tfLCLue=−LC(I+Lu1|tfLC)−1Le1|tf

From the above equations, and note that the setpoint is uncorrelated with the noise, least squares method can be employed to estimate the related matrix from N batch sets of data as

(60)y1|tf(1)y1|tf(2)⋯y1|tf(N)⏟YCL=Lyrr1|tf(1)r1|tf(2)⋯r1|tf(N)⏟RCL+Lyee1|tf(1)e1|tf(2)⋯e1|tf(N)⏟ECL

(61)Lyr=YCLRCLYCL=y1|tf(1)y1|tf(2)⋯y1|tf(N)RCL=r1|tf(1)r1|tf(2)⋯r1|tf(N)

L_ur is estimated in the same way. With the estimation of the subspace matrices, the predictions Y^CL,U^CL can be written as

(62)Y^CLU^CL=LyrLurRCL

4.1 Estimation of the process model

From Equation (59)

(63)Lur=LC(I−Lyr)

(64)Lu1|tfLC=(I+Lu1|tfLC)Lyr=Lyr+Lu1|tfLCLyr

So

(65)Lyr=Lu1|tfLC(I−Lyr)=Lu1|tfLur

Thus, the process model is obtained as follows:

(66)Lu1|tf=Lyr(Lur)−1

4.2 Estimation of the noise model

Note that there is always at least one step delay in the real closed-loop process, that is to say, Lu_1|tf, L_C, Le_1|tf are all low triangle matrices and

(67)Lye=(I+Lu1|tfLC)−1Le1|tf=Im00⋯0b2,1Im0⋯0b3,1b3,2Im⋯0⋮⋮⋮⋱⋮btf,1btf,2btf,3⋯Im

where b_i,j in L_ye are coefficients to be identified. Considering Equation (60) and noting the form of L_ye in Equation (67), we can obtain

(68)YCL(1:m,:)−Y^CL(1:m,:)=ECL(1:m,:)=e1(1)e1(2)⋯e1(N)

And the following equation also holds

(69)YCL(m+1:2m,:)−Y^CL(m+1:2m,:)=b2,1e1(1)e1(2)⋯e1(N)+e2(1)e2(2)⋯e2(N)

Substituting e1(1)e1(2)⋯e1(N) obtained from Equation (68), the least square method can be applied to Equation (69) to estimate b_2,1 and residual e2(1)e2(2)⋯e2(N).

Repeat the above least square method in every block row [YCL((i−1)m+1:im,:)−Y^CL((i−1)m+1:im,:)](i=1,2,⋯tf) by substituting et(1)et(2)⋯et(N)(t=1,2,⋯i−1) estimated from the previous steps, and we can get all the coefficients b_i,j in L_ye.

Then from Equations (57)–(59), note that

(70)LC=−(Luee1|tf)(Lyee1|tf)†=−(UCL−U^CL)(YCL−Y^CL)†

Finally we get

(71)Le1|tf=(I+Lu1|tfLC)Lye=I−Lyr(Lur)−1(UCL−U^CL)(YCL−Y^CL)†Lye

5. Procedures of CPM for the whole batch process

In order to carry out CPM for the whole batch process, the benchmark covariance of joint input and output should be available for every instant k. The procedures of CPM consist of the following key steps:

1. Collect N(N > tf) historical batches with golden operation, which achieve satisfactory control performance, as benchmark datasets and scale them for identifying the observer Markov parameters.

2. Compute process model Markov parameters Lu_t|tf and noise model Markov parameters Le_t|tf at every instant t from observer Markov parameters.

3. Compute Φ_t(k) and Γ_t(k) from Lu_t|tf and Le_t|tf to establish benchmark covariance Cov(stopt,stopt) of joint input and output under LQG control at every time instant t.

4. Extract the data ste=(u1eT,y1eT,u2eT……uteT,yteT)T of real-time batch at every instant, scale it and calculate its Tt2 index by Equation (22). Draw the Hotelling's T² control chart to check if the current Tt2 exceeds the upper limit, from which we can know whether the batch process at sampling time is in optimal condition.

6. Case study

In this section, the proposed framework is applied to a batch reactor example in Luyben (Luyben, 1989) (Section 5.7) in which the differential equations and the parameter values describing the reaction process are detailed.

The reaction system involves two consecutive first-order reactions: A → B → C and two stages are run in the system. In the first stage, the steam in the jacket initially heats up the reactor content until the exothermic heat of reaction generated is significantly enough. In the second stage, the cooling water in the jacket is used to remove the exothermic heats of reaction. Two split-ranged control valves, i.e. a steam valve and a water valve, control the reactor temperature. The instrumentation is all pneumatic and the controller output, P_c varies from 3 to 15 psig. Feedback control is used to eliminate the disturbance effect on the manipulated variables and keep the heat requirement with the batch run. The duration of each batch is 80 min. The sampling time of each batch is 1 min. The process is affected by a persistent disturbance (N(0, 1)) in the two control valves of the nominal upstream pressures. The measured jacked and the measured reactor temperatures with noise (N(0, 0.01)) are included. Also, the initial concentration C_A with noise (N(0,0.0032)) is introduced. To get the desired qualities at the end of the batch run, the temperature reference profile or setpoint (T_s) is shown in Figure 1.

The proportional plus integral controller (PI) is tuned for the above system and the control parameters as chosen as [3, 0.1/s]. The controlled variable, temperature output (T), with the corresponding controller output (P_c) of one typical batch are also illustrated in Figure 1. A total of 50 batch runs of the data under golden operation condition (GOC) with above PI controller parameters are collected as benchmark data to establish benchmark covariance under LQG control at every time instant.

The controlled output consists of two parts, temperature and pressure, respectively. The temperature is the setpoint variables and also belongs to the monitored outputs. The temperature setpoint consists of four sections. The first one is a linear increase from 80 to 200 Fahrenheit degree, and the second is a maintaining process at the top degree. And the third one is the decline section while the fourth one is a new balance in the top picture in Figure 1.

There are three lines in the chart, Ts for temperature setpoint trajectory in red as well as the golden data (GOC) in blue and Case I in actual process in black. It is easy to find that the temperature was tracked by and large and that the GOC data had a lower fluctuation than Case I.

However, the condition of controlled pressure state is not as good as temperature in the below chart in Figure 1. The golden data in blue color has a smaller variance or a kind of much slighter fluctuation than that of real process in Case I. And the mean value of pressure is around 8 psig. The blue line has an undulation between 6 to 10, which shows that the maximum deviation of intensity of pressure is 2 psig. In the meantime, the upper and lower bounds of black line are 15 and 5, respectively, which implies that the pressure condition is not ideal. Therefore, Case I is in bad condition and should be monitored and detected.

Two more methods including similar benchmark covariance are compared with the proposed method (Method I). Method in reference (Kadali and Huang, 2002) (Method II) considers only the autocovariance of input and that of output at the present moment, respectively, and gives the benchmark covariance as

Cov(utopt,utopt)00Cov(ytopt,ytopt)

and does not include cross-covariance of input and output. Based on Method II, we can define the following benchmark covariance (Method III):

Cov(utopt,utopt)Cov(utopt,ytopt)Cov(utopt,ytopt)TCov(ytopt,ytopt).

Both of them do not take past covariance information into consideration. The controller performance assessment results of these three methods are compared in the following two cases.

Figure 2 is the method we introduce for the batch process monitoring. The excess of threshold from black points towards the red curve signifies that the condition is deteriorating in data and some measure in next control should be considered.

But Figure 3 with Method II and Figure 4 with Method III only detected a slight difference out of gauge. In other word, two more methods show an insensitive detection of a light worsening condition.

6.1 Case II: scrap value impulse noise effect

In industrial environment, the scrap value impulse noise often occurs, resulting from temporary failure of the sensor and data sampling system. It does not mean that the behavior of the batch operation itself is not good. Suppose the scrap value impulse noise affects the batch process at the instant 20 and 30, the controlled temperature output (noise effect) is plotted in Figure 5. The corresponding Hotelling's T² control charts of this batch by three methods are plotted in Figure 6, Figure 7 and Figure 8.

Figure 5 exhibits a normal correct condition for batch process no matter in temperature tracking or in pressure fluctuations. The responses to the correct data in Figure 7 and Figure 8 return two wrong assessment results. Method II in Figure 7 and Method III in Figure 8 show two excess of the index line in red color by the black points. However, Method I we proposed has no outstripping evaluation towards the red line.

Method II as well as Method III is formed at time “t-th”. Therefore, their covariance benchmarks (in red) are fixed at that moment. Our method to compute the covariance is based on LQG, which adopts the vectors form time “1-st”(Group 1) to “t-th”(Group t) in Equation (20). Therefore, our benchmark has a dynamic red curve to monitor the process, which is not as rigid as the next two methods.

It can be seen that the Hotelling's T² of our method do not behave so sensitive because it contains all past covariance information of input and output, and the inertia is much greater. While the Hotelling's T² in Method I and II only focus on the input and output at the present instant. So it reacts at once, exceeds the T² limit. But in fact, the batch process operation is still normal after the impulse noise appears. So Method I and II are much sensitive to the noise and trigger a false alarm which tells the process behavior is become worse. In summary, the method of this paper is more adequate for the latent and slow variation and is not so sensitive to the sudden impulse noise.

7. Conclusion

In this paper, a data-driven controller performance monitoring method based on LQG approaches and Hotelling's T² index for LTV batch process is proposed. Open/closed loop input and output data are applied to identify the process model as well as the disturbance model, both in Markov parameter form. Then the optimal covariance of joint input and output can be obtained by the LQG method. The Hotelling's T² control chart can be established to monitor the controller performance of a new batch. An industrial batch reactor example is studied and clearly demonstrated the effectiveness and merits of the proposed method.

Figures

Figure 1

Controlled outputs of Operation I for batch reactor

Figure 2

The Hotelling's T² control chart of Case I in the Method I of this paper

Figure 3

The Hotelling's T² control chart of Case I in Method II

Figure 4

The Hotelling's T² control chart of Case I in Method III

Figure 5

Controlled outputs in Operation II for batch reactor

Figure 6

The Hotelling's T² control chart of Case II in the Method I of this paper

Figure 7

The Hotelling's T² control chart of Case II in Method II

Figure 8

The Hotelling's T² control chart of Case II in Method III

Funding: This research was funded by the National Natural Science Foundation of P.R. China (NSFC: 62073286).