Plant-model mismatch estimation from closed-loop data for state

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Plant−Model Mismatch Estimation from Closed-Loop Data for StateSpace Model Predictive Control¶ Jodie M. Simkoff,‡ Siyun Wang,‡ Michael Baldea,*,‡ Leo H. Chiang,§ Ivan Castillo,§ Rahul Bindlish,§ and David B. Stanley§ ‡

McKetta Department of Chemical Engineering, The University of Texas at Austin 200 East Dean Keeton Street, Stop C0400, Austin, Texas 78712, United States § The Dow Chemical Company, Freeport, Texas 77541, United States ABSTRACT: The area of controller performance monitoring, assessment, and diagnosis for model predictive control (MPC) has seen an increase in interest in recent years. A frequently identified cause of degraded performance is mismatch between the plant model used in the controller and the true dynamics of the plant. Most recent research focuses on locating plant− model mismatch in order to reduce the considerable effort required to re-identify the model. In this paper, we present a novel autocovariance-based plant−model mismatch estimation approach for unconstrained MPC based on linear state space models. We show that (additive) plant−model mismatch can be quantified as the solution of an optimization problem which minimizes the discrepancy between the sample autocovariance of plant outputs and the corresponding value obtained from theoretical predictions. We illustrate our theoretical results with two simulation case studies, demonstrating good performance in estimating parametric mismatch.

1. INTRODUCTION Model predictive control (MPC) has been widely applied in industry in the recent decades. It has gained particularly strong traction in the chemical industry, where it was shown that MPC can perform much better than multiloop feedback control, especially for multivariate systems with significant interactions. Additionally, the MPC framework facilitates incorporation of constraints on inputs, states, and outputs.1,2 The model-centric nature of the MPC paradigm inherently renders the controller vulnerable to plant−model mismatch (PMM). Mismatch may be structural, that is, unmodeled dynamics, or parametric, whereby the model structure is correct but the parameters are not. Such mismatch arises naturally over the lifespan of a physical system; in the chemical industry, it is due among others to phenomena such as corrosion, catalyst deactivation, fouling, and plugging of pipes and ducts. A realization of this fact has spurred research efforts toward defining the means to monitor, diagnose and repair the performance of MPC systems. Loosely stated, the goal of controller performance monitoring (CPM) approaches is to assess how well a controller is performing and to detect deviations from normal operation. CPM can be carried out on several levels. First, one must identify when control performance has degraded. Often this is done by way of a metric that

can be computed from routine operating data, where the current value is compared to the corresponding value computed using data from a “golden,” benchmark period of good operation. Another category of CPM benchmarks compares the current performance of the plant with the best achievable performance; the most famous of these are the Minimum Variance benchmark3−5 and the Linear Quadratic Gaussian benchmark,6−8 but user-defined benchmarks have also been proposed.9−11 Then, controller performance monitoring efforts can focus on diagnosing the cause of performance degradation. This is particularly challenging in the context of MPC, where performance can degrade as a result of multiple factors; improper controller tuning, unmeasured disturbances, and errors in the plant or disturbance models are chief among them.12,13 Such diagnoses can rely on benchmarking the value of the MPC objective function.14,15 Notably, Schafer and Cinar14 propose a method to discern between errors in tuning and the presence of unmeasured disturbances, and changes in measurable disturbances and/or process dynamics, as the potantial source of the control performance change. In a similar vein, a result by Harrison and Qin16 allows for discriminating between mismatch in the disturbance and plant models based Received: Revised: Accepted: Published:

¶

A preliminary version of this material was submitted for presentation at the 2017 IEEE Conference on Decision and Control © XXXX American Chemical Society

A

November 27, 2017 February 12, 2018 February 20, 2018 February 20, 2018 DOI: 10.1021/acs.iecr.7b04917 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research on the order of the autocorrelation in the Kalman Filter innovations. If errors in the process model are identified as the root cause of the degraded control, practitioners must decide how to proceed. System identification comprises a significant portion of the initial MPC commissioning (Sun et al. estimate that it accounts for as much as 80% of the implementation effort13). Therefore, frequent, comprehensive model reidentification and updating are undesirable from an economic point of view. To address this problem, researchers have developed methods for isolating mismatch to a set of submodels which should be reidentified. Badwe et al.17 used partial correlation analysis to identifying input−output channels where mismatch is present and, in a subsequent work,18 quantified the effect of these mismatches on control performance via the designed sensitivity function. Botelho et al.19 proposed a related method using a nominal sensitivity function, which is independent of set point. Bachnas et al.20 proposed a data-driven adaptive MPC scheme based on representing the model with a finite set of orthogonal basis functions and iteratively reidentifying the function coefficients. In a similar vein, Kheradmandi et al.21 have recently proposed a framework that pairs a model quality monitoring metric based on prediction errors with a closedloop reidentification step, which is undertaken as indicated by the quality metric. Our approach differs from those described above in that it both isolates and quantifies plant−model mismatch from existing closed-loop data sets (unlike Badwe et al.17,18 and Botelho et al.19), without performing a full model reidentification online (as in Bachnas et al.20 and Kheradmandi et al.21). In our previous work,22,23 we made initial steps toward simultaneously detecting the mismatched parameters and providing a quantitative estimate of parametric mismatch. Our approach was based on a closed-form expression for the theoretical value of the process output autocovariance matrices. The plant−model mismatch was estimated by solving an optimization problem aimed at minimizing the discrepancy between this theoretical value and the sample estimate of the autocovariance of process outputs computed from routine operating data. For simplicity, these initial developments considered controller formulations based on parametric (transfer function) or nonparametric, finite step response model structures. However, the most recent developments in the literature point toward state-space models becoming the plant representation of choice for future practical applications, given that such models can capture a wider range of dynamics and noise models in a more parsimonious format than convolution models.2 Motivated by this, in the present work, we develop an autocovariance-based approach for estimating MPC plant−model mismatch for controllers using stochastic state-space models, with a Kalman filter (KF) performing state estimation. We demonstrate that an asymptotically consistent estimate can be obtained when the size of the available data set tends to infinity; nevertheless, the two case studies we present show that very reasonable results are available from finitedimensional data sets.

x(k + 1) = Ax(k) + Bu(k) + w(k) y(k) = Cx(k) + v(k)

(1)

with x(k) ∈ R , u(k) ∈ R , and y(k) ∈ R . The noise signals are Gaussian random variables, where w(k) ∈ Rn ∼ N(0, Qw) and v(k) ∈ Rp ∼ N(0, Rv) The plant is under model predictive control. A Kalman filter (KF) estimates the process states at each time step. A schematic of the system is shown in Figure 1. n

m

p

Figure 1. Schematic of KF-MPC loop.

The model used in both the MPC and KF is xμ(k + 1) = A μxμ(k) + Bμu(k) yμ(k) = Cμxμ(k)

(2)

The objective of this work is to estimate parametric (not structural) mismatch in the matrices A and B of the state-space model: ΔA = A − A μ ΔB = B − Bμ

(3)

2.2. MPC Formulation. We first consider an unconstrained MPC controller which solves, at each time step k, the optimization problem min J = Ê (k + 1)T QÊ (k + 1) + ΔU(k)T RΔU(k)

ΔU(k)

(4)

where Ê (k + 1) is the predicted output error over prediction horizon P, ΔU(k) is the sequence of optimal control moves over control horizon M, and Q and R are weighting matrices. An overview of the MPC framework can be found in several published works.24,25 We note that this approach is not limited to unconstrained MPC, and its application to constrained MPC is discussed in Remark 1 in section 4. The unconstrained problem has the explicit solution: o

ΔU(k) = K cÊ (k + 1)

(5)

where the control law Kc is a function of the controller tuning and dynamics, and the unforced predicted error is defined as o Ê (k + 1) = Yr (k + 1) − Ŷ (k + 1|k)

(6)

where Yr(k + 1) is the reference trajectory and Ŷ (k + 1|k) is the (pP)-dimensional vector of future unforced (i.e., if no control moves were made) predicted outputs. A controller designed as above would not, however, guarantee offset-free control in the presence of either plant− model mismatch or unmeasured disturbances.26,27 In practical applications, integral action must therefore be incorporated into

2. PRELIMINARIES 2.1. Problem Statement. We consider a linear timeinvariant (LTI) discrete-time system, where the plant is represented as B

DOI: 10.1021/acs.iecr.7b04917 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

⎡ CB ⎤ ̅ ̅ ⎢ ⎥ ⎡ CA ⎤ ̅ ̅ ̅ ̅ ̅ CB ̅ ̅ ⎢ CAB ⎥ ⎢ ⎥ ⎢ ⋮ ⎥ 2 ... ⋱ ⎢ CA ̅ ̅ ⎥ ⎢ ⎥ = M̅ = ⎢ G , ̅ ⎥ ⎢ ⋮ ⎥ ⋮ ... ... CB ̅ ̅ ⎢ ⎥ ⎢ ⎥ P⎥ ⎢⎣ CA ⋮ ⋮ ... ... ̅ ̅ ⎦ ⎢ ⎥ ⎢⎣ CA P − 1B P−M−1 ⎥ CA B̅ ⎦ ̅ ̅ ̅ ̅ ̅

the controller design. A common way of doing so is via a constant output disturbance model. We note that the theory in the following sections applies equally regardless of the type of disturbance model used (state/input or output). The augmented process model becomes ̃ ̃(k) + Bu ̃ (k ) x̃(k + 1) = Ax ̃ ̃ (k ) y(k) = Cx

(7)

(14)

where,

The form of eq 13 leads, when substituted in eq 4, to the following result for the control law Kc introduced in eq 5

⎡ x(k ) ⎤ ⎥ x ̃ (k ) = ⎢ ⎢⎣ d(k)⎥⎦

K c = (G̅ TQG̅ + R)−1G̅ TQ

(8)

2.3. State Estimation. At each time step, a Kalman filter is used to estimate the process and disturbance states.

and, for a constant output disturbance, ⎡ Aμ 0 ⎤ ⎡ Bμ ⎤ Ã = ⎢ ⎥ , B̃ = ⎢ ⎥ , C̃ = [C I ] ⎣ 0 I⎦ ⎣0 ⎦

∼x̂(k + 1|k) = Ã ∼x̂(k|k) + B̃ u(k) μ μ ̃ ∼̂(k + 1|k) ŷ (k + 1|k) = Cx

(9)

26

Pannochia and Rawlings have shown that the number of disturbance states d(k) should be equal to the number of measured outputs in order to guarantee offset-free control. The added states provide this guarantee in the presence of either an external disturbance (entering in the outputs, states, or inputs), and/or an effective disturbance due to plant−model mismatch. The latter is the topic under investigation in this work. We thus begin by considering the scenario where plant−model mismatch is the only source of offset; that is, the plant− model mismatch is the sole cause of biased state estimates. Gaussian white noise is still affecting the system as described in eq 1, and the covariance of the noise signals is considered to be known correctly (see ref 16 for a method that can be used to distinguish between PMM and incorrect noise covariances). Since the decision variable in eq 4 is Δu(k) while the statespace model of the process is in terms of u(k), we must define a second augmented system in order to obtain the control law Kc in eq 5:

∼x̂(k + 1|k + 1) = ∼x̂(k + 1|k) + L(k + 1)[y(k + 1) − ŷ(k + 1|k)] (17)

where L is the optimal steady-state Kalman filter gain obtained by solution of the algebraic Riccati equation.28 In principle, a deterministic observer design, that is, a Luenberger observer, could also be used. As it will be shown in the following sections, the matrix L is simply a parameter in the mismatch estimation procedure.

3. AUTOCOVARIANCE-MISMATCH RELATION We now establish an expression for the autocovariance of the process outputs as a function of known parameters and plant− model mismatch. We first make a number of assumptions that enable this development. Assumption 1. The model in eq 2 is a minimal state-space realization. Assumption 2. There is constant, time-invariant mismatch in the state and input matrices A and B. The output matrix C is fixed. Since all minimal realizations of a given system are related by a similarity transformation, fixing C does not restrict the identification of plant−model mismatch. Assumption 3. Time series data reflecting periods of steadystate operation of the process are available. This also implies that, if a Kalman Filter is employed, the filter gain L has converged to a steady-state value. Assumption 4. The noise covariances matrices Qw and Rv are known accurately. Our approach then is to express the process outputs as an autoregressive moving average (ARMA) process driven by state and measurement white noise vectors. We define the Green’s function (impulse response) coefficients of this process, from which the autocovariance function of the outputs easily follows. These developments are presented below:

(10)

where ⎡ x(k ) ⎤ ⎢ ⎥ x̅ (k + 1) = ⎢ d(k) ⎥ ⎢ ⎥ ⎢⎣ u(k − 1)⎥⎦

(11)

and ⎡ A μ 0 Bμ ⎤ ⎡ Bμ ⎤ ⎢ ⎥ ⎢ ⎥ A̅ = ⎢ 0 I 0 ⎥ , B̅ = ⎢ 0 ⎥ , C̅ = [C I 0 ] ⎢⎣ 0 0 I ⎥⎦ ⎢⎣ 0 ⎥⎦

(12)

We can then write Ŷ (k + 1|k) = Mx ̅ ̅ (k) + G̅ ΔU(k)

(16)

The measurement update for the a posteriori state estimate is computed as follows:

x̅ (k + 1) = Ax ̅ ̅ (k) + B̅ Δu(k) y(k) = Cx ̅ ̅ (k )

(15)

(13)

where M̅ is the observability matrix predicted one-step ahead, and G̅ is a block-Toeplitz matrix: C


Article

Industrial & Engineering Chemistry Research ∞ ⎧ ⎪ ⎨ y(k) = ∑ CA B ∑ γw , j w(k − i − j) ⎪ i=1 ⎩ j=1

We begin by using eq 1 together with eqs 5, 16, and 17 to express u(k) as a linear combination of past inputs, and state and measurement noise:

N

k−i

+

which leads to a Green’s function representation of y(k) similar to eq 20: ∞

y(k) =

i=1

=

∑ θi w(z)z

Λw , i

(18)

Λv , i

+

i=1

∑ ωi v(z)z

∑ γv ,i v(z)z−i

i=1

i=0

∞

4. PLANT−MODEL MISMATCH DIAGNOSIS AND ESTIMATION In the previous section we established an explicit expression for the expected value of the autocovariance of the process outputs, as a function of known model parameters and the additive plant−model mismatch terms we would like to estimate. These developments allow us to introduce our main result: Theorem 1 Consider a system of the form (1), under the control of an unconstrained model predictive controller as in eq 4. Assume that a data set of length Ns collected during routine operation is available, and that Assumptions 1−4 are met. Then, the solution to the minimization problem

−i

i=0

Ly

(20)

min J =

ΔA, ΔB

∑ (R yy(π ) − R̂ yy(π ))T (R yy(π ) − R̂ yy(π )) π=0

(27)

is an asymptotically consistent estimator of the plant−model mismatch ΔA, ΔB. In eq 27, R yy (π) ∈ R p×p is the unbiased sample autocovariance estimator R yy(π ) =

k−i

1 Ns − π

Ns − π

∑ k=1

y(k)y(k + π ) (28)

are the entries of the sample estimate of the autocovariance matrix from the closed-loop data, and Ly is a user-defined parameter specifying the number of autocovariance lags considered in the estimation.

Bu(k − i) + Cw(k − 1) + v(k)

i=1

i=0

(26)

N

∑ CA

∞

i=1

The Green’s function coefficients, γw,i and γv,i, of this MA(∞) process can be obtained implicitly.29 Next, we express eq 1 in terms of the impulse response of the plant and current noise terms. Recall that we assume the states and outputs are only affected by Gaussian noise with known covariance, and that plant−model mismatch is the only persistent “disturbance” acting on the system. y(k) =

i = 1, ..., ∞

R yy(π ) = Q w ∑ (Λw , iΛw , i + π ) + R v ∑ (Λv , iΛv , i + π )

∞

∑ γw ,i w(z)z−i +

⎧I i = 1 ⎪ i =⎨ j−1 ⎪∑ C(A μ + ΔA) (Bμ + ΔB) × γv , i − j ⎩ j=1

Recalling that w(k) ∼ 5(0, Q w) and v(k) ∼ 5(0, R v), the autocovariance of y(k) is

and obtain the Green’s function coefficients via the partial fraction expansion of each of the two terms on the right-hand side, which correspond to the process and measurement noise vectors, respectively. ∞

i = 2, ..., ∞

(25)

N−1 −i

⎧C i = 0 ⎪ i−1 =⎨ j ⎪∑ C(A μ + ΔA) (Bμ + ΔB) × γw , i − j ⎩ j=1

(24)

(19)

u (z ) =

(23)

i=0

where

N

∑ ϕi u(z)z

∞

∑ Λw ,iw(k − i) + ∑ Λv ,iv(k − i) i=1

The coefficients are determined recursively using the MPC and Kalman filter equations, and are defined in Appendix A. We note that eq 18 is a vectorial autoregressive moving average (VARMA) process,29 and that we can define the coefficients up to a sufficient horizon N such that the effect of inputs and noise signals on u(k) is negligible for t < k − N. To arrive at a discrete-time transfer function representation of eq 18, we take the z-transform: u (z ) −

⎭

(22)

⎡ v(k) ⎤ ⎢ ⎥ ⎢ v(k − 1) ⎥ + [ω0 ω1 ... ωN − 1]⎢ ⎥ ⎢⋮ ⎥ ⎢ v(k − N + 1)⎥ ⎣ ⎦

−i

∑ γv ,i v(k − i − j)⎬⎪ + Cw(k − 1) + v(k) j=0

⎡ w(k − 1) ⎤ ⎢ ⎥ ⎢ w(k − 2) ⎥ ⎥ + [θ1θ2 ... θN ]⎢⋮ ⎢ ⎥ ⎢ w(k − N )⎥ ⎢⎣ ⎥⎦

N

⎫ ⎪

∞

⎡ u(k − 1) ⎤ ⎢ ⎥ ⎢ u(k − 2) ⎥ u(k) = [ϕ1ϕ2 ... ϕN ]⎢ ⎥ ⎢⋮ ⎥ ⎢ u(k − N )⎥ ⎣ ⎦

(21)

Substituting eq 20 in eq 21, we have D


Article

Industrial & Engineering Chemistry Research Proof. Since the sample autocovariance estimator in eq 28 is unbiased, as the number of samples approaches infinity it will tend to the theoretical value in eq 26. Then eq 27 will attain its lowest possible value at 0. Remark 1. The unconstrained linear MPC framework introduced in section 2.2 is used for clarity of exposition in the derivations that follow in section 3. However, the autocovariance-based approach can be easily extended to handle constrained linear MPC in the following way:23 First, the available data should be segmented into operating periods during which the active set of the controller remains the same, since the derivation of eq 26 requires an explicit control law. Then, for each data segment, the appropriate control law is used to compute the coefficients in eq 18. If no constraints are active, the control law is equal to the constrained control law. If a constraint is active, the corresponding column in the input matrix Bμ should be replaced with a vector of zeros before computing eqs 18−26. Then the modified PMM estimation optimization problem becomes S

min J =

ΔA, ΔB

Remark 4. In contrast to Theorem 1, practical implementation requires that the sums in eq 20 and 23−and as a result eq 26−be computed up to a finite number of terms. Thus, we utilize high-order moving average models rather than infinite order models, thereby introducing a nonzero bias. However, if the order of the moving average model is much larger than the dominant time constant of the discretized process, the effect on the PMM estimates is negligible.

5. SIMULATION CASE STUDIES 5.1. Square MIMO System under Unconstrained Linear MPC. We consider the following multi-input multioutput (MIMO) system where n = m = p = 3. The model used in the MPC is ⎡−0.4 0 0.6 ⎤ ⎡ 0.9 0.1 0 ⎤ ⎢ ⎥ ⎢ ⎥ A μ = ⎢ 0 0.8 0 ⎥ , Bμ = ⎢ 0.5 −2 0 ⎥ , C = I3 ⎢⎣ 0.1 0 0.5⎥⎦ ⎢⎣ 0.8 0 1 ⎥⎦ (30)

Ly

The weighting parameters in the MPC objective function are Q = I3 and R = 0.1I3, and the sample time is 1 min. The process and measurement noise covariances, which are assumed to be known accurately, are

∑ ∑ (R y ̅ y ̅ ,s(π ) − R̂ y ̅ y ̅ ,s(π ))T s=1 π=0

(R y ̅ y ̅ , s(π ) − R̂ y ̅ y ̅ , s(π ))

(29)

⎡ 0.02 0.004 0 ⎤ ⎢ ⎥ Q w = ⎢ 0.004 0.01 0 ⎥ , R v = 0.06I3 ⎢⎣ 0 0 0.03⎥⎦

where S is the number of data segments, the subscript y̅ indicates that the data have been mean-centered, and all other terms are defined as in eq 27. Remark 2. Measurable disturbances and reference tracking can also be handled by performing the expectation operation on inputs and outputs in eqs 18 and 21, and then subtracting deterministic components that arise due to measurable disturbances or reference trajectory terms.22 Remark 3. The presence of unmeasured disturbances besides Gaussian noise during the operating period of the data collection will bias plant−model mismatch estimates obtained using this data set and the proposed autocovariance-based approach. The use of sufficient disturbance states means that the controller will still provide offset-free tracking of feasible setpoints. This problem of unmeasured disturbances may be encountered during any modeling or model reidentification procedure. In industrial implementations, it is sometimes assumed that these unmeasured disturbances are short in duration, so that if the data set collected is large relative to the length of the disturbance, their effect will not be significant. On the other hand, it is possible that disturbances persist throughout the relevant operating period. In this case, it is very difficult in practice to distinguish between plant−model mismatch and sustained, unmeasured disturbances. However, since the plant−model mismatch estimate will incorporate the effect of unmeasured disturbances on the resultant behavior of the system, the estimated plant may provide improved performance as long as the disturbance continues to be present at the same level. A recent work by Harrison and Qin16 proposed a method for discriminating between process plant−model mismatch and an incorrect noise covariance matrix by analysis of the Kalman filter innovations; this technique is applicable to the class of systems considered in this paper and could be employed to perform additional performance diagnostics. However, the more general problem of discriminating between process plant−model mismatch and disturbance plant−model mismatch remains an open research topic.30

(31)

Since there are three outputs, it is necessary to augment the state-space with three disturbance states. In this system, they are treated as constant output disturbances so that the augmented system is of the form in eq 9. The noise covariance matrix for the disturbance states is 0.1I3, and these noise signals are assumed uncorrelated with the process state noise. Then, by solving the algebraic Riccati equation we obtain the optimal steady-state Kalman filter gain: ⎡ 0.062 ⎢ ⎢ 0.020 ⎢ 0.0027 L=⎢ ⎢ 0.67 ⎢−0.0195 ⎢ ⎣ 0.0083

0.0017 0.027 − 0.0020 0.0071 0.6847 0.0035

0.0094 ⎤ ⎥ 0.0023 ⎥ 0.051 ⎥ ⎥ − 0.0095 ⎥ − 0.0029 ⎥ ⎥ 0.6475 ⎦

(32)

We now introduce mismatch into some parameters of the discrete-time state space model, such that ⎡ 0 −0.06 0 ⎤ ⎡0 0 0⎤ ⎢ ⎥ ⎢ ⎥ ΔA = ⎢ 0 0.01 0 ⎥ , ΔB = ⎢ 0 0.8 0 ⎥ ⎣0 0 ⎣0 0 0⎦ 0⎦

(33)

To explore the effect of data set size on PMM estimation performance, 100 Monte Carlo simulations were performed in MATLAB for each of five data set sizes: N = 1000, N = 5000, N = 10000, N = 20000, and N = 50000. For each trial, the optimization problem in eq 27 was formulated using Ly = 20 lags and solved using IPOPT31 via the MATLAB Opti Toolbox.32 To assess the quality of the estimates, one approach is to establish, for example, 95% confidence intervals on each of the parameter mismatch estimates (i.e., each element of ΔA and E


Article


Table 1. Results of Mismatch Estimation in State Space Model Parameters as a Function of Data Size (Case Study 5.1) F-statistic (Fαp,n−p = 0.05 = 1.73) 2-norm of mismatch estimate covariance matrix sum of squared errors of mean step response

N = 1000

N = 5000

N = 10000

N = 20000

N = 50000

2.84 6.54 × 10−3 26.0

2.07 6.61 × 10−3 8.10

1.50 5.59 × 10−3 3.68

1.19 3.25 × 10−3 0.552

1.10 1.14 × 10−3 0.116

ΔB). However, by computing individual, scalar confidence intervals on each parameter mismatch estimate we are not accounting for the significant covariance between the elements of the parameter mismatch vector. Therefore, we proceed to establish a multivariate confidence region for this vector by performing an F-test, based on the following fact:33 p(n − 1) −1 Fp , n − p (x̅ − μ)T Σ̂ x̅ (x̅ − μ) ∼ n−p

N = 5000 data set is 8.10, while the SSE between the original MPC model and the true plant is 2458. Figures 2 and 3 show the step responses of the model used in the MPC, the “true plant” and the mean value of the new

(34)

where p in eq 34 is the number of estimated parameters, n is the number of trials, x ̅ is the sample mean, Σ̂x̅ is the sample covariance of the mean, and μ is the true mean. If the F-statistic computed using eq 34 with μ set to be the true value of the plant−model mismatch vector, is less than the critical value at a given significance level, then we should fail to reject the null hypothesis. The latter postulates that the mean value computed using our proposed approach is the same as the true value of the plant−model mismatch. The results of applying this assessment to the five data sets described above are presented in Table 1. In the first row, the F-statistic for each of the data set sizes is shown. For 100 trials, the critical F-value at the 0.05 significance level is 1.73. For the N = 1000 and N = 5000 experiments, the F-statistic fails the Ftest by a small margin, indicating bias in the estimates. The three larger data sets have an F-statistic lower than the critical value, so the true parametric mismatch vector is contained in the multivariate 95% confidence region. This trend with respect to data set length is in accordance with Theorem (1). Table 1 also presents the matrix norm of the sample covariance matrix computed from the mismatch estimate vectors, as a measure of the magnitude of the variance of the estimates. In addition to the reduction in bias, as the data set size increases, the variance of the estimates also decreases. This is expected, since the covariance of the unbiased sample autocovariance estimator decreases with sample size as follows:34 ⎛ 1 ⎞ cov(R yy, π) = O⎜ ⎟ ⎝ Ns − π ⎠

Figure 2. Open loop step response comparison: mean value of estimates for N = 5000 sample case.

model obtained by adding the PMM estimate to the original model. Since, as seen in Table 1, the bias in parameter estimates is small, the updated step response in Figure 2 is quite a good representation of the plant for the 5000 sample data set. A marginal improvement in the SSE of the updated model for the 20000 sample data set can be seen in 3. To conclude, although the shorter data sets in these simulation studies may fail the Ftest at the 5% significance level, the magnitude of bias in the estimates is relatively small. 5.2. Nonsquare MIMO System under Constrained Linear MPC. We now consider a nonsquare system with three states, two inputs, and three outputs:

(35)

where Ns in eq 35 is the size of the data set, and π is the lag considered. Finally, despite the indication of some bias in parametric mismatch estimates from the shorter data sets, the mean value of the parametric mismatch estimates is quite close to the true value for these data sets. One way to quantify this is to consider the sum of squared errors (SSE) between the true plant response and the updated model resulting from adding the mean parametric mismatch estimate to the original model used in the MPC. This is presented in the last row of Table 1. It is evident that the SSE value decreases as the data set size increases and thus the difference between the expected value and true value (the bias) is reduced. However, we can see that this bias is small even for the shorter data sets; for example the SSE between the updated and true plant step responses for the

⎡ 0.9 −0.1 0 ⎤ ⎡ 0.6 0.2 ⎤ ⎢ ⎥ ⎢ ⎥ A μ = ⎢ 0.1 0.8 −0.1⎥ , Bμ = ⎢ 0.3 −1⎥ , C = I3 ⎢⎣ 0 0.8 ⎥⎦ ⎢⎣−0.1 −0.1 0.6 ⎥⎦ (36)

Outputs y1 and y3 are selected as the controlled variables. The weighting parameters in the MPC objective function are Q = I2 and R = 10I2, and the sample time is 1 min. The presence of three measured outputs requires that three disturbance states be added to the model. The same noise model in eq 31 is used. For this model, the steady-state Kalman filter gain obtained by solving the algebraic Riccati equation is F


Article

Industrial & Engineering Chemistry Research ⎡ 0 −0.06 0 ⎤ ⎢ ⎥ ΔA = ⎢ 0 0.01 0 ⎥ , ⎣0 0 0⎦

⎡ 0 −0.2 ⎤ ⎢ ⎥ ΔB = ⎢ 0.2 0 ⎥ ⎣ 0 −0.4 ⎦

(38)

Simulations were carried out as follows: for each data set size (N = 5400, N = 10800, and N = 21600), 30 Monte Carlo trials were performed. In each, two set point changes are made with respect to the controlled variables, dividing the data into three segments of equal length. The first set point change is achieved by the controller without changing its active set, while the second set point change drives u1 to its lower bound, and the controlled variables are maintained at their best achievable levels (the set points are not feasible). Figure 4 shows the outputs and inputs for one such simulation data set (N = 10800). The results are presented in Table 2. The shortest data set has an F-statistic that is slightly above the 5% significance Table 2. Results of Mismatch Estimation in State Space Model Parameters as a Function of Data Size (Case Study 5.2) (Fαp,n−p = 0.05

= 2.40) F-statistic 2-norm of mismatch estimate covariance matrix sum of squared errors of mean step response

Figure 3. Open loop step response comparison: mean value of estimates for N = 20000 sample case.

⎡ 0.1893 ⎢ ⎢ 0.0017 ⎢−0.0026 L=⎢ ⎢ 0.2848 ⎢ 0.0487 ⎢ ⎣−0.0324

− 0.0372 ⎤ ⎥ 0.1157 − 0.3068 ⎥ −0.0391 0.1518 ⎥ ⎥ −0.0657 0.0360 ⎥ 0.3005 − 0.0058 ⎥ ⎥ 0.0196 0.2791 ⎦ 0.0716

N = 5400

N = 10800

N = 21600

2.98 13.9 × 10−3

1.19 6.77 × 10−3

2.36 4.22 × 10−3

8.96

1.34

1.22

critical value, while the two longer data sets do contain the true parametric mismatch estimates within the 95% multivariate confidence region. Again we note that the parameter variance− covariance matrix norm and the SSE of the mean step response both decrease with data size. Thus, the approach performs well for this nonsquare, constrained system. It should be noted that this performance depends upon the availability of data sets where all of the manipulated variables are free (not at their bounds) for at least part of the data collection period, as well as upon the correct specification of the MPC for nonsquare systems such that the inputs and outputs are stationary processes at steady-state.

(37)

In this case study, there are constraints on the manipulated variables: u1 and u2 are constrained to within ±0.5 of their nominal operating values. We now introduce the following mismatch in the state and input matrices:

Figure 4. Input−output data for case study 5.2 (N = 10800). G


Article


Figure 5. Open-loop step response comparison: N = 5400, 30 trials.

Figure 6. Open loop step response comparison: N = 10 800, 30 trials.

In Figures 5−7, we plot the step responses of the MPC model, the true plant and the updated models obtained in all 30 trials (where each trial differs only by the random noise added to the states and outputs). These plots demonstrate the variability of the resulting updated models that may be obtained and show that as the data set size increases, the updated models are becoming more closely centered around the true plant value. Moreover, it is apparent that, in closed-loop, the output

autocovariance is more sensitive to the accuracy of the openloop response models of some input−output channels and less sensitive to others, especially the u2−y3 channel in this simulation case study. In a practical setting, of course, one has only a single realization of this stochastic process with which to perform the plant−model mismatch estimation. However, since simulations like those presented above can easily be performed offline, one H


Article


Figure 7. Open loop step response comparison: N = 21 600, 30 trials.

■

APPENDIX A In this appendix we define the coefficients vectors ϕ, θ, and ω which were introduced in eq 18. These coefficients were obtained using the control law in eq 5 and Kalman filter equations 16 and 17 together with the expression for the evolution of the true plant in eq1, assuming plant−model mismatch is the only persistent disturbance affecting the system.

approach is to generate, via simulation, a confidence region around the nominal model case (all elements of the mismatch vector are zero). Then, this confidence region can be used to assess the statistical significance of a parametric mismatch estimate obtained by applying eq 29 to a closed-loop data set.

6. CONCLUSIONS Performance monitoring and assessment for model predictive controllers has been the subject of much research activity in recent years. Although MPC has been widely used for at least three decades, many implementations in industry report poor performance; plant−model mismatch is an oft-cited culprit. A comprehensive framework for detecting and correcting plant− model mismatch in these control loops is thus desirable for practitioners, who currently must weigh the negative effects of poor quality control (e.g., reduced productivity, increased emissions, damage to equipment) against the prospect of a costly system reidentification effort to improve the model. Although the recent literature offers several techniques for identifying submodels affected by mismatch or general sources of error in MPC, our approach is the firstto our knowledgeto both locate and quantify plant−model mismatch from routine closed loop data without the requirement of external excitation. In this paper we have extended our formulation to include model predictive controllers using statespace models and online state estimation, since such representations are increasingly prevalent in industrial applications. Future work will focus on applying the proposed method to industrial data sets and providing a comprehensive evaluation of performance compared with other available methods recently proposed in the literature, as well as on extending the theory to the case of nonlinear systems.

⎧ ̃ ̃ ⎪ I m − K c ,1M1Fϕ , m − K c ,1M 2 m = 1 ϕm = ⎨ ⎪ m = 2, ..., N ⎩−K c ,1M̃ 1Fϕ , m

(39)

where Fϕ,m is defined as

Fϕ , m

⎧(B̃ − LCB̃ ) + LC(B + ΔB) μ ⎪ ⎪ m=1 ⎪ ⎨ = (Ã − LCÃ )Fϕ , m − 1 + LC(A μ ⎪ m−1 ⎪ + ΔA) (Bμ + ΔB) ⎪ m = 2, ..., N ⎩ (40)

We also have

θm = −K c ,1M̃ 1Fθ , m

(41)

where ⎧ LC m=1 ⎪ = ⎨(Ã − LCÃ )Fθ , m − 1 + LC(A μ m = 2, ..., N ⎪ m−1 ⎩ + ΔA)

(42)

ωm = {−K c ,1M̃ 1(Ã − LCÃ )m L m = 0, ..., N − 1

(43)

Fθ , m

I


Article

Industrial & Engineering Chemistry Research In eqs 39−43, Kc,1 is the (m × nP) submatrix of the controller gain matrix in (5) corresponding to the vector of inputs that are actually implemented at time k.

■

(19) Botelho, V.; Trierweiler, J. O.; Farenzena, M.; Duraiski, R. Methodology for Detecting Model Plant Mismatches Affecting Model Predictive Control Performance. Ind. Eng. Chem. Res. 2015, 54, 12072−12085. (20) Bachnas, A. A.; Weiland, S.; Tóth, R. Data driven predictive control based on model structures with orthonormal basis functions. IEEE Conference on Decision and Control 2015, 3026−3031. (21) Kheradmandi, M.; Mhaskar, P. Model predictive control with closed-loop re-identification. Comput. Chem. Eng. 2018, 109, 249−260. (22) Wang, S.; Simkoff, J. M.; Baldea, M.; Chiang, L. H.; Castillo, I.; Bindlish, R.; Stanley, D. B. Autocovariance-based MPC model mismatch estimation for systems with measurable disturbances. J. Process Control 2017, 55, 42−54. (23) Wang, S.; Simkoff, J. M.; Baldea, M.; Chiang, L. H.; Castillo, I.; Bindlish, R.; Stanley, D. B. Autocovariance-based plant−model mismatch estimation for linear model predictive control. Systems & Control Letters 2017, 104, 5−14. (24) Seborg, D. E.; Edgar, T. F.; Mellichamp, D. A.; Doyle, F. J., III Dynamics and Control, 3rd ed.; Wiley, 2011. (25) Rawlings, J. B.; Mayne, D. Q. Model Predictive Control: Theory and Design; Nob Hill Publishing, 2009. (26) Pannocchia, G.; Rawlings, J. B. Disturbance Models for OffsetFree Model-Predictive Control. AIChE J. 2003, 49, 426−437. (27) Muske, K. R.; Badgwell, T. A. Disturbance modeling for offsetfree linear model predictive control. J. Process Control 2002, 12, 617− 632. (28) Kalman, R. E.; Bucy, R. S. New results in linear filtering and prediction theory. J. Basic Eng. 1961, 83, 95−108. (29) Box, G.; Jenkins, G.; Reinsel, G. Time Series Analysis: Forecasting and Control, 3rd ed.; Prentice-Hall, Inc., 1994. (30) Sun, Z.; Qin, S.; Singhal, A.; Megan, L. Performance monitoring of model-predictive controllers via model residual assessment. J. Process Control 2013, 23, 473−482. (31) Wachter, A.; Biegler, L. T. On the Implementation of a PrimalDual Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming. Mathematical Programming 2006, 106, 25− 57. (32) Currie, J.; Wilson, D. I. OPTI: Lowering the Barrier Between Open Source Optimizers and the Industrial MATLAB User. Conference: Foundations of Computer-Aided Process Operations. Savannah, Georgia, 2012. (33) Anderson, T. W. An Introduction to Multivariate Statistical Analysis; Wiley: New York, 1984; Vol. 2; pp 163−165. (34) Odelson, B. J.; Rajamani, M. R.; Rawlings, J. B. A new autocovariance least-squares method for estimating noise covariances. Automatica 2006, 42, 303−308.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Michael Baldea: 0000-0001-6400-0315 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Partial financial support from The Dow Chemical Company is acknowledged with gratitude. Jodie M. Simkoff was partially supported by the NSF Graduate Research Fellowship Program through Grant DGE-1610403.

■

REFERENCES

(1) Qin, S. J.; Badgwell, T. A. A survey of industrial model predictive control technology. Chem. Eng. Prog. 2003, 11, 733−764. (2) Darby, M. L.; Nikolaou, M. MPC: Current practice and challenges. Contr. Eng. Prac. 2012, 20, 328−342. (3) Harris, T. J. Assessment of Closed Loop Performance. Can. J. Chem. Eng. 1989, 67, 856−861. (4) Ko, B.-S.; Edgar, T. F. Performance assessment of constrained model predictive control systems. AIChE J. 2001, 47, 1363−1371. (5) Yu, J.; Qin, S. J. MIMO control performance monitoring using left/right diagonal interactors. J. Process Control 2009, 19, 1267−1276. (6) Shah, S.; Patwardhan, R.; Huang, B. Multivariate controller performance analysis: methods, applications and challenges. AIChE Symp. Ser. 2002, 190−207. (7) Kadali, R.; Huang, B. Controller performance analysis with LQG benchmark obtianed under closed loop conditions. ISA Trans. 2002, 41, 521−537. (8) Harris, T.; Boudreau, F.; Macgregor, J. Performance assessment of multivariable feedback controlers. Automatica 1996, 32, 1505−1518. (9) Gao, J.; Patwardhan, R.; Akamatsu, K.; Hashimoto, Y.; Emoto, G.; Shah, S.; Huang, B. Performance evaluation of two industrial MPC controllers. Control Engineering Practice 2003, 11, 1371−1387. (10) Julien, R.; Foley, M.; Cluett, W. Performance assessment using a model predictive control benchmark. J. Process Control 2004, 14, 441− 456. (11) Yu, J.; Qin, S. J. Statistical MIMO controller performance monitoring. Part I: Data-driven covariance benchmark. J. Process Control 2008, 18, 277−296. (12) Patwardhan, R. S.; Shah, S. L.; Qi, K. Z. Assessing the performance of model predictive controllers. Can. J. Chem. Eng. 2002, 80, 954−966. (13) Sun, Z.; Qin, S. J.; Singhal, A.; Megan, L. Performance monitoring of model-predictive controllers via model residual assessment. J. Process Control 2013, 23, 473−482. (14) Schafer, J.; Cinar, A. MPC system performance assessment, monitoring, and diagnosis. J. Process Control 2004, 14, 113−129. (15) Zagrobelny, M.; Ji, L.; Rawlings, J. B. Quis custodiet ipsos custodes? Annu. Rev. Control 2013, 37, 260−270. (16) Harrison, C. A.; Qin, S. J. Discriminating between disturbance and process model mismatch in model predictive control. J. Process Control 2009, 19, 1610−1616. (17) Badwe, A. S.; Gudi, R. D.; Patwardhan, R. S.; Shah, S. L.; Patwardhan, S. Detection of model-plant mismatch in MPC applications. J. Process Control 2009, 19, 1305−1313. (18) Badwe, A.; Patwardhan, R.; Shah, S.; Patwardhan, S.; Gudi, R. Quantifying the impact of model-plant mismatch on controller performance. J. Process Control 2010, 20, 408−425. J


Plant-model mismatch estimation from closed-loop data for state

Recommend Documents