Handling Model Plant Mismatch in State Estimation Using a Multiple

Subscriber access provided by Eastern Michigan University | Bruce T. Halle Library

Article

Handling model plant mismatch in state estimation using a multiple model based approach Kevin Arulmaran, and Jinfeng Liu Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b00234 • Publication Date (Web): 19 Apr 2017 Downloaded from http://pubs.acs.org on April 25, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Handling model plant mismatch in state estimation using a multiple model based approach Kevin Arulmaran and Jinfeng Liu∗ Department of Chemical & Materials Engineering, University of Alberta, Edmonton, AB T6G 2V4, Canada E-mail: [email protected] Phone: +1 780 492 1317. Fax: +1 780 492 2881 Abstract Accurate state estimates are important for the success of MPC. State estimates are obtained using a model but in real plants there will always be model plant mismatch (MPM) which affects these estimates. In this work, we present a multiple model based approach to obtain unbiased state estimates in the presence of MPM. Necessary assumptions on the source of mismatch and models used are presented. It is shown that unbiased output estimates do not guarantee unbiased state estimates. Our approach is shown to provide unbiased state estimates when all the assumptions are met using a froth flotation system. A model identification based control approach using our multiple model estimation approach with a conventional MPC was tested on the froth flotation system and was found to successfully provide offset free reference tracking when all the necessary assumptions for unbiased state estimation were met. A nonlinear offset free MPC was also tested on the froth flotation system but was not able to provide offset free reference tracking as some necessary conditions were not met.

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

State estimation plays a fundamental role in the success of model predictive control (MPC). 1 Kalman filtering and moving horizon state estimation (MHE) are commonly used state estimation methods in MPC applications (e.g. Refs. 2–4). However, these model based state estimation methods are sensitive to model-plant mismatch (MPM) which is always present in real systems due to factors such as external disturbances, parameter drift, sensor failures and others. 5 Due to the popularity of MPC, it becomes necessary to design state estimation methods that provide accurate results even in the presence of MPM. Numerous approaches have been proposed in the literature to address this issue. One approach, as described in Ref. 6, is to use a data driven model to capture the effects of the mismatch. However, while this method can provide unbiased state estimates it does not provide any information about the cause of the mismatch. This is because the data driven model has no physical significance while the states and parameters in chemical systems do. An alternative approach, as described in Ref. 7, is to reduce model plant mismatch by updating the system parameters at each time step by solving an optimization problem. This approach provides diagnostic information through the updated parameter estimates but has a high computational load as an optimization problem has to be solved at each time step. If a large number of parameters have to be optimized then computation time can be significant. Similar to offset free model predictive control, it is possible to augment the system model with constant step disturbances and use special observer designs to obtain unbiased estimates. 8,9 However these approaches suffer the same downside as that in Ref. 6 in that they do not provide any information about the potential cause of the mismatch. A common approach in tracking and fault detection applications for mechanical and electrical systems is the use of Multiple Model (MM) methods. 10–15 These methods work by having multiple models (each model representing one mode of the system) run in parallel. Model probabilities and the overall estimate are obtained according to the rules of the particular method being used. It is possible to use augmented models with MM methods to 2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


obtain updated parameter estimates which provide diagnostic information. 10 MM methods offer the advantage of low computational cost while still providing diagnostic information, however there are only limited examples of their use in chemical systems. 16,17 Furthermore, these works focus on unbiased output estimation and control which does not necessarily guarantee unbiased state estimation. The design of an appropriate model set to use with MM methods is also challenging as there needs to be enough separation between models based on output residuals 15 and enough models to capture the range of system dynamics but using too many models will decrease performance. 18 In this paper, we present an algorithm for the use of a MM method for state estimation in the presence of MPM caused by parameter mismatch. The model set used in our approach includes models augmented with different parameters (augmented models do not have the same states). Guidelines for model set design and assumptions on model properties are also presented. Simulations examples are included to illustrate the effectiveness of this approach. Although the proposed approach has been successfully applied to nonlinear systems, it has only been proved for linear systems. While the simulation results are obtained using the extended Kalman filter (EKF), the proposed approach does not restrict the choice of estimator thus MHE (or any other estimator) can be used as well. Finally, the MM estimation approach is combined with a conventional MPC and compared with a traditional offset free MPC.

2

Problem description

Consider a nonlinear plant whose dynamics are characterized as follows:

x˙ = f (x, θ∗ , u, w)

(1a)

y = h(x) + v

(1b)

3



Page 4 of 36

where θ∗ is the vector of actual plant parameters (θ∗ ∈ Rm ), x is the state vector (x ∈ Rn ), u is the vector of inputs (u ∈ Rp ) and y is the vector of outputs (y ∈ Rq ). w and v represent process noise and measurement noise respectively. We assume a model of the nonlinear plant is developed and has the form:

x˙ = f (x, θ, u, 0)

(2a)

y = h(x)

(2b)

where θ is the vector of model parameters (θ ∈ Rm ). A common approach to deal with model parameter uncertainty is to augment the uncertain parameters as states and use the augmented model in estimators. 10 Parameters can be added to the state vector to obtain augmented systems of the form:     ˜ ˜ x, θ, u, 0) x˙ x  f (˜ x˜˙ =   =   0 x˙ θ

(3a)

˜ x) y = h(˜

(3b)

where x˜ represents the augmented state vector (˜ x ∈ Rn˜ =n+r ) consisting of the original state vector (xx ∈ Rn ) and the added parameters (xθ ∈ Rr ). θ˜ is the vector of remaining parameters (θ˜ ∈ Rm−r ). The objective of this work is to outline a Multiple Model (MM) approach that uses augmented models of the form (3) and is capable of providing unbiased state estimates when the following assumptions are satisfied: A1 The model and plant have the same structure with mismatch caused by model parameters not having the same values as plant parameters, i.e. θ 6= θ∗ . Furthermore, the plant parameters (θ∗ ) are time invariant or change slowly compared to the process dynamics.

4


Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


A2 Each parameter has a unique effect on the output, i.e. each unique set of parameter values results in a unique y trajectory for a given constant input. A3 The effects of the parameters are distinguishable from noise. In A1, parameters are allowed to vary provided their values change slowly compared to the process dynamics such that this will not affect estimator convergence. A2 and A3 are needed to ensure that the parameters can be distinguished from each other based on noisy output measurements which is a necessary condition for the success of any MM method. 15 As A2 and A3 are needed for parameter estimation in general, they are not restrictive. Some additional assumptions are also required and these are presented in Section 4.2. Remark 1 While the proposed approach is applied to cases of parameter mismatch, it can also be applied to cases of strucutral mismatch. The main challenge in this case is determining an appropriate model set. If the system has known failure modes with associated models then these can be used to create the model set provided that the models have adequate separation based on output trajectories.

3

Illustrative example - Two CSTRs

In this section, we present a plant which will be used to illustrate each step of our algorithm. The plant model is obtained from Ref. 19 and a schematic diagram is provided in Figure 1. The plant consists of two well-mixed, non-isothermal continuous stirred tank reactors (CSTRs) where three parallel, irreversible, elementary, exothermic reactions take place: k

k

k

1 2 3 A− → B, A − → U and A − → R. A is the reactant, B is the desired product and U and R are

5



Page 6 of 36

undesired byproducts. Under standard modeling assumptions, the plant model is: 3

X Fr F0 Q1 (T0 − T1 ) + (T2 − T1 ) + T˙1 = Gi (T1 )cA1 + V1 V1 ρcp V1 i=1

(4a)

3

c˙A1 =

X F0 Fr (cA0 − cA1 ) + (cA2 − cA1 ) − Ri (T1 )cA1 V1 V1 i=1

(4b)

3

X F3 Q2 F1 (T1 − T2 ) + (T03 − T2 ) + Gi (T2 )cA2 + T˙2 = V2 V2 ρcp V2 i=1

(4c)

3

c˙A2

X F3 F1 (cA1 − cA2 ) + (cA03 − cA2 ) − Ri (T2 )cA2 = V2 V2 i=1

(4d)

where Ri (Tj ) = ki exp(−Ei /RTj ), Gi (Tj ) = (−∆Hi /(ρcp ))Ri (Tj ) for j = 1, 2. Tj , cAj , Qj and Vj represent the temperature of the reactor, the concentration of A, the rate of heat input to the reactor and the reactor volume respectively with the subscript representing the CSTR number. ∆Hi , ki , Ei , i = 1, 2, 3 represent the enthalpies, pre-exponential constants and activation energies of the three reactions respectively. cp and ρ are the heat capacity and density of fluid in the reactor. The state vector is x = [T1 , cA1 , T2 , cA2 ]T and the input vector is u = [Q1 , Q2 ]T . There are two measured outputs, y = [x1 , x3 ]T . Nominal parameter values, steady states and associated inputs are provided in Table 1. Table 1: Nominal parameters, steady states and steady state inputs for two CSTRs F0 F1 F3 Fr V1 V2 R T0 T03 cA0 cA03 ∆H1 ∆H2 ∆H3

= = = = = = = = = = = = = =

4.998 m3 /h k1 3 39.996 m /h k2 3 30.0 m /h k3 34.998 m3 /h E1 1.0 m3 E2 3 3.0 m E3 8.314 kJ/kmol K cp 300.0 K ρ 300.0 K T1s csA1 4.0 kmol/m3 3 T2s 2.0 kmol/m −5.0 × 104 kJ/kmol csA2 −5.2 × 104 kJ/kmol Qs1 −5.4 × 104 kJ/kmol Qs2

6

= = = = = = = = = = = = = =

3.0 × 106 h−1 3.0 × 105 h−1 3.0 × 105 h−1 5.0 × 104 kJ/kmol 7.53 × 104 kJ/kmol 7.53 × 104 kJ/kmol 0.231 kJ/kg K 1000.0 kg/m3 303.7 K 2.5 kmol/m3 302.9 K 2.3 kmol/m3 1.0 × 105 kJ/h 1.0 × 105 kJ/h


Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


F0, T0, CA0

Fr, T2, CA2 F3, T03, CA03

F1, T1, CA1 F2, T2, CA2

CSTR 1 CSTR 2

Q1 Q2

Figure 1: Schematic diagram of 2 CSTR example The objective is to estimate the four states based on the two measured outputs in the presence of parameter uncertainty. In this example, we show that an estimator based on augmented models may give unbiased output estimates but that this does not guarantee unbiased state estimates. To emphasize the effects of the parameters, process and measurement noise are not considered in this example. Let us consider a scenario where the actual F3 = 50 m3 /h (for an unidentified reason) while the value used in the model is F3 = 30 m3 /h (the nominal value). Figure 2 shows the trajectories of the estimates given by an EKF using a model augmented with T0 and T03 and a non-augmented model (where parameter uncertainty is not considered). The initial state of the EKF is the nominal steady state (when F3 = 30 m3 /h). It can be seen from the figure that there is parameter uncertainty as the actual outputs move away from their nominal values. It can also be seen that if parameter uncertainty is not considered (EKF with a nonaugmented model), the EKF gives biased output estimates as well as biased state estimates. Further, it can be seen that the EKF with the augmented model gives unbiased output estimates even though the augmented parameters are not the ones that are mismatched and the state estimates are biased. Thus it can be concluded that unbiased output estimates do

7



not always guarantee unbiased state estimates. From the analysis presented in Section 4.2 for linear systems, it can be seen that this is caused by the choice of augmented models and not by the choice of estimator. This motivates our present work to develop a procedure for state estimation in the presence of parameter mismatch where unbiased output estimates guarantee unbiased state estimates. Our work will include guidelines on augmented model design, estimator selection and model selection.

303 303.8 302

302.8

x 3 (K)

x 1 (K)

303.6 302.78

303.4 303.2

302.76 3.6

303

3.7

302.5

3.8

301.9 301.8 3.6

302

3.7

3.8

302.8 0

1

2

3

4

0

1

Time (h)

2

3

4

3

4

Time (h)

2.5

2.3 Actual T 0 T 03 model

x 4 (kmol/m 3 )

x 2 (kmol/m 3 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 36

Non-augmented model

2.45

2.4

2.25

2.2

2.15 0

1

2

3

4

0

Time (h)

1

2

Time (h)

Figure 2: Trajectory of the actual states and estimates given by an EKF using a model augmented with T0 and T03 and a non-augmented model with model plant mismatch in F3 .

4

Proposed Solution

This section provides an overview of our procedure for state estimation using a MM approach. A flowchart of the propsed procedure is presented in Figure 3. The different steps are explained in the following subsections.

8


Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ĞƚĞƌŵŝŶĞƉĂƌĂŵĞƚĞƌƐĞƚŵŽƐƚůŝŬĞůǇƚŽ ĐŽŶƚĂŝŶƵŶĐĞƌƚĂŝŶƚǇͬŵŝƐŵĂƚĐŚ

ƌĞƚŚĞƌĞƉĂƌĂŵĞƚĞƌƐĨŽƌ ǁŚŝĐŚĂůůŽƵƚƉƵƚƐŚĂǀĞůŽǁ ƐĞŶƐŝƚŝǀŝƚǇ͍

zĞƐ

ZĞŵŽǀĞƚŚĞ ƉĂƌĂŵĞƚĞƌ;ƐͿ ĨƌŽŵƚŚĞ ƉĂƌĂŵĞƚĞƌƐĞƚ

EŽ

ƌĞĂƚĞĂƵŐŵĞŶƚĞĚŵŽĚĞůƐ

EŽ

ƌĞĂůůŵŽĚĞůƐ ŽďƐĞƌǀĂďůĞ͍

ZĞŵŽǀĞ ƵŶŽďƐĞƌǀĂďůĞ ŵŽĚĞůƐ

zĞƐ ĂƌƌǇŽƵƚƐƚĂƚĞĞƐƚŝŵĂƚŝŽŶĂŶĚŵŽĚĞů ƉƌŽďĂďŝůŝƚǇĐĂůĐƵůĂƚŝŽŶ

KďƚĂŝŶŽǀĞƌĂůůĞƐƚŝŵĂƚĞ

Figure 3: A flowchart of the proposed procedure for obtaining unbiased state estimates in the presence of parameter uncertainty.

9



4.1

Parameter set selection

The next step is to determine the sensitivity of the outputs to each parameter in S to ensure that the parameters satisfy assumption A2. If a parameter has a negligible effect on all the outputs it will not be possible to detect mismatch in that parameter using the output data even if there is no noise. There are four possible results for parameter sensitivity: 1. Outputs and unmeasured states are sensitive to the parameters. 2. Outputs are sensitive to the parameters but unmeasured states are not sensitive. 3. Outputs are not sensitive to the parameters but the unmeasured states are sensitive. 4. Neither outputs nor unmeasured states are sensitive to the parameters. In the first two cases the parameters should be left in S. In the last case the parameters can be removed from S without significantly affecting estimation performance. In the third case, state estimates will be affected by the parameters however it will not be possible to detect mismatch in these parameters based on the outputs. Thus these parameters can be removed from S as well. Additional measurements can be taken if accurate estimates of the affected states are critical (such as for safety reasons). Sensitivity is quantitative so while a sufficiently large mismatch in the parameters from cases three and four could be detected, a sufficiently large mismatch may be unlikely or even impossible based on the physical constraints of the system. As a result, leaving these parameters in S will increase computation time while not improving estimation performance as the model set will be larger due to the inclusion of models augmented with these insensitive parameters. Two methods to calculate parameter sensitivity are provided in the following subsections.

4.1.1

Step test

The step test only considers steady state information and can only handle variations in a single parameter. However, unlike the gramian approach discussed in the next section, it 10


Page 10 of 36

Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


can identify the sensitivity of unmeasured states to the parameters and is faster to calculate than a gramian. As such, it can be used as a preliminary test to reduce the size of S before carrying out the gramian test. The procedure for the step test is as follows: 1. Carry out simulations using the system model (2) and nominal parameter values (θnom ) s to find the nominal steady state (ynom ).

2. For a given non-zero input, increase (or decrease) a parameter (θi ) by a fixed percentage while holding the other parameters unchanged. Run the simulation to find the new steady state. s s (yjs − yj,nom )/yj,nom for j = 1, . . . , q 3. Calculate the normalized steady state gain Kij = (θi − θi,nom )/θi,nom (each output).

4. Repeat steps 2 and 3 for i = 1, . . . , s (each parameter in S). For a given θi , if the sum of the absolute normalized steady state gains Ki =

Pq

j=1

|Kij |

is smaller than a predetermined threshold (i.e., Ki < ), then the outputs are said to be insensitive to that parameter and the parameter should be removed from S. As step test results depend on the operating region, the parameter sensitivity should be checked over a typical range of manipulated input values. The reason non-zero inputs are used in step tests is based on a property of linear systems. In general, stable non-singular linear systems of the form:

x˙ = Ax + Bu

(5)

will converge to a zero steady state if u = 0 regardless of A. If u = 0 parameters will have no effect on the steady state and it will not be possible to distinguish between them based on step tests. Thus nonzero inputs are imposed for the sake of generality.

11



4.1.2

Page 12 of 36

Empirical observability gramian approach

The observability gramian provides a quantitative measure of observability however is difficult to calculate analytically for nonlinear systems. As a result, empirical observability gramians have been developed as local approximations to the analytical gramians for nonlinear systems. 20 Using these gramians for sensitivity calculations provides two major advantages over the use of step tests: transient information is used and it is possible to handle variations in multiple parameters. However gramian calculation is more computationally intensive than the step test. The procedure for calculating empirical gramians is as follows: 20 1. Pick a nominal operating point (states and parameters). Add the parameters under consideration to the state vector to obtain an augmented state vector (˜ x ∈ Rn˜ ). 2. Define step change directions for the augmented state vector. One option is to use two level factorial design but partial factorial designs may be used to reduce computation time. These directions are listed in a semi-orthogonal matrix T with each column being one direction. 3. Define nd step magnitudes (positive values) on a percentage scale denoting cd as the magnitude for d = 1, . . . , nd and define a scaling matrix S to convert the steps from percentage to actual values. 4. Define a number of experiments with one experiment for each combination of step direction and magnitude. The initial condition for each experiment is given by x˜id (0) = cd ST ei + xñom where ei is a standard unit vector in Rn˜ and the superscript id denotes the experiment. 5. Simulate the trajectory of each experiment until it reaches steady state. This data can be used to calculate the observability gramian (WO ) as follows:

WO =

nd X d=1

1 nd c2d 12

Z

tid f

T Ψd (t) T T dt

0


(6)

Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


d where tid ˜×n ˜ matrix f is the time to steady state for experiment id and Ψ (t) is an n

with the ij element defined as:

T jd jd jd Ψdij (t) = (y id (t) − y id (tid f )) (y (t) − y (tf ))

(7)

The observability gramian has the structure 



n×n n×m WXθ  WX  WOn˜ ×ñ =   m×n m×m WθX Wθ

(8)

where Wθ is the identifiability gramian of the parameters. The eigenvalues of Wθ provide information about the identifiability of the parameters with smaller eigenvalues indicating lower sensitivity. Parameters with significant components along eigenvectors of eigenvalues below a certain threshold will not be identifiable. This threshold will be zero in an ideal case but in reality will always be non-zero due to numerical errors. 20 The threshold can be approximated as the Frobenius-norm of E where E is found by linearizing the system about its nominal point and taking the difference between the empirical observability gramian and analytical observability gramian of this linearized system. 4.1.3

Two CSTR example

The results of the sensitivity tests for the two CSTR example are summarized in Table 2. Step tests were carried out by increasing parameters by 20% from their nominal values one at a time. Gramian tests were conducted by changing parameters one at a time by ±10% and ±20% of their nominal values. Both empirical gramian and step test indicate that F0 , Fr , F3 , k1 , k2 , k3 , V1 and V2 are not sensitive and should be removed from S. Removing these parameters results in S = {T0 , T03 }. Remark 2 In Table 2, it shows that k1 (the main reaction constant) has low sensitivity. It is important to note that this only means that changes in k1 do not greatly affect the outputs. 13



Page 14 of 36

Table 2: Results of parameter sensitivity analysis using step tests and empirical observability gramians for two CSTRs

F0 Fr F3 k1 k2 k3 V1 V2 T0 T03

Step test Empirical Gramian (Sum of absolute normalized gains) (Normalized eigenvalues) 0.004 7.150 × 10−9 0.002 7.507 × 10−8 0.014 5.810 × 10−6 2.578 × 10−6 6.032 × 10−11 −11 1.127 × 10 1.156 × 10−16 1.127 × 10−11 2.751 × 10−15 0.001 4.443 × 10−8 0.002 2.278 × 10−9 0.404 0.002 1.964 1.000

It does not necessarily mean that changes in k1 will not significantly affect unmeasured states (for example if the unmeasured states are only weakly observable). However for this system and operating region, changes in k1 do not significantly affect the unmeasured states. It is also important to note that these sensitivity results are only valid in a local region near the operating point and that changes in k1 may affect the states/outputs more significantly in different operating regions.

4.2

Augmented model creation

All augmented models used in the model set must satisfy certain conditions to ensure that unbiased output estimates will guarantee unbiased state estimates. These conditions are derived below for deterministic linear systems operating at steady state with mismatch in only the A matrix. Consider a stable (only negative eigenvalues) linear plant denoted by Aplant . This plant is modelled by Amodel . As the plant and model have the same structure, the model plant mismatch can be denoted by ∆A where:

Aplant = Amodel + ∆A

14


(9)

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with Aplant , Amodel and ∆A all having the same dimensions. In a similar manner, the state mismatch can be denoted by ∆x where:

xplant = xmodel + ∆x

(10)

with xplant , xmodel and ∆x all having the same dimensions. As there is no mismatch in the B matrix, Bplant = Bmodel = B. The plant dynamics can then be written as:

x˙ plant = Aplant xplant + Bu = (Amodel + ∆A)(xmodel + ∆x) + Bu yplant = Cxplant

(11a) (11b) (11c)

The model dynamics are given by:

x˙ model = Amodel xmodel + Bu

(12a)

ymodel = Cxmodel

(12b)

For a given constant nonzero input us , the system will converge to a nonzero steady state. At this steady state, we can check if output estimates are unbiased by directly comparing s s the values of yplant and ymodel . If we have unbiased output estimates, the following holds:

s s = ymodel yplant

(13a)

Cxsplant = Cxsmodel

(13b)

C(xsmodel + ∆xs ) = Cxsmodel

(13c)

C∆xs = 0

15


(13d)


Page 16 of 36

As B and u are known exactly, the following is also true at this steady state:

Aplant xsplant = −Bus

(14)

Amodel xsmodel = −Bus

(15)

From this we can derive the following:

Aplant xsplant = Amodel xsmodel

(16a)

(Amodel + ∆A)(xsmodel + ∆xs ) = Amodel xsmodel

(16b)

Amodel ∆xs + ∆Axsmodel + ∆A∆xs = 0

(16c)

For any row where the model matches the plant perfectly (i.e. any row where Aplant = Amodel ) the corresponding row of ∆A will be zero. The minimum number of zero rows in ∆A will occur when the mismatched parameters and augmented parameters are all on different rows. This minimum (nmin ) is equal to the number of rows (n) minus the number of mismatched parameters and augmented parameters (nmis + r), i.e. nmin = n − (nmis + r). Let Amodel,i represent row i of Amodel for i = 1, . . . , n. For any row where ∆A = 0, (16c) simplifies to Amodel,i ∆xs = 0. This can be combined with (13d) to obtain: 



Amodel,i  s   ∆x = 0 C | {z }

(17)

T

for i = {Rows of Amodel corresponding to zero rows of ∆A}. If T is full column rank, there exists a unique solution ∆xs = 0 which implies unbiased state estimation. T can only be full column rank if the number of rows (n − (nmis + r) + q) is greater than or equal to the number of columns (n).

16


Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


This can be summarized in the following additional assumptions used to guarantee unbiased state estimates if unbiased output estimates are obtained (regardless of the estimator used): A4 The number of outputs must be greater than or equal to the number of mismatched parameters plus the number of augmented parameters i.e. q ≥ nmis + r. A5 As the location of the mismatch is not known, T must be full column rank for ALL possible combinations of nmin rows of Amodel . Assumption A4 can be directly checked for nonlinear systems while assumption A5 can be checked using linearized system models. Although not proven, it can be seen from simulation results that the mismatch can be correctly identified if these two assumptions (and the ones presented in Section 2) hold. Once the number of parameters to augment has been decided, augmented models can be created by augmenting xmodel with elements of ∆A (the resulting models will be nonlinear). An estimator can then be used to obtain estimates of the augmented state and thus the mismatched parameters. For the two CSTR example, there are two outputs and we assume only one parameter is mismatched thus we can augment one parameter. This results in two augmented models: one augmented with T0 and the other augmented with T03 . Remark 3 Note that A4 and A5 do not guarantee that the mismatched parameter has been identified correctly but only that state estimates are unbiased. Note also that the proposed procedure is capable of handling cases where the mismatch is caused by the simultaneous effects of multiple parameters. A4 only limits the number of mismatched parameters that can be handled. It does not matter if the mismatched parameters are independent or described by a joint probability distribution.

17



4.3

Page 18 of 36

Model Observability

Given an observable linear system (A, C), the number of constant step disturbance states added must be less than or equal to the number of outputs to ensure that the resulting linear ˜ C) ˜ is observable. 1 In our algorithm, the added states are constant augmented system (A, parameters thus the augmented model will be nonlinear even if the original system is linear. As a result, the above condition does not guarantee observability of the augmented model and the observability of each augmented model must be checked before it can be added to the model set. A nonlinear model is observable if the nonlinear observability matrix Q is full rank, 21 where Q is defined as:    Q=  



˜ x) dh(˜ .. . ˜ x)) d(Lnf˜˜−1 h(˜

    

(18)

˜ (1 ≤ k ≤ n ˜ x) = [h ˜ 1 (˜ ˜ q (˜ ˜ − 1) the k-th with h(˜ x), . . . , h x)]T , n ˜ the size of the system and Lkf˜h ˜ with respect to f˜. This method requires the calculation of high order Lie derivative of h Lie derivatives which may be difficult to obtain. An alternative approach is to linearize the ˜ C) ˜ where (A˜ ,C) ˜ system along a typical trajectory to obtain linearized system models (A, are found by taking the Jacobian of (3a) and (3b) respectively at a given x˜ and u as shown below: 

∂ f˜1 ∂x ˜1

  . . A˜ =   . 

∂ fñ ∂x ˜1

··· .. .

∂ f˜1 ∂x ñ ˜

···

∂ fñ ∂x ñ ˜





 ..  .   

˜1 ∂h ∂x ˜1

  . . C˜ =   .  ˜

∂ hq ∂x ˜1

x ˜=˜ x(t),u=u(t)

18


··· .. .

˜1 ∂h ∂x ñ ˜

···

∂ hq ∂x ñ ˜



 ..  .    ˜ x ˜=˜ x(t),u=u(t)

(19)

Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The linearized models are observable if the following condition is satisfied:   ˜     C      C˜ A˜  rank  ..    .         C˜ Añ˜ −1

          ˜  =n        

(20)

where n ˜ is the size of the A˜ matrix. Any unobservable models should be removed from the model set.

4.4

State Estimation

After designing the model set, state estimation should be carried out using all the models in the model set simultaneously. As the proposed approach does not restrict the choice of estimator, we are free to select the best estimator for each application. If system constraints need to be handled then a common choice is the moving horizon estimator.

4.5

Overall Estimate

At each time step, the state estimates from the model with the highest probability are used. If the estimator chosen provides innovation covariance (e.g. EKF), model probability is found according to the formula used in the Autonomous Multiple Model (AMM) algorithm which is based off Bayes theorem: 22

Li,k pi,k

−1 exp(−0.5 · Ti,k · Si,k · i,k ) = |2πSi,k |1/2 pi,k−1 Li,k = t P pj,k−1 Lj,k−1

(21a) (21b)

j=1

where Li,k represents the likelihood of model i at time k, pi,k represents the probability of model i at time k and t is the number of models in the model set. i,k represents the 19



Page 20 of 36

innovation of model i at time k with Si,k representing the innovation covariance. An artificial lower limit δ is set on model probabilities with any probability that drops below δ set to δ. This prevents models from becoming inactive as the recursive nature of the probability calculation means that if a model probability drops to zero then it cannot be nonzero in future time steps. 17 If innovation covariance information Si,k is not available due to the choice of estimator (e.g. MHE), an alternate option for calculating model probabilities is to use the following equation for model likelihood: 17

Li,k = exp(−0.5 · Ti,k · ∆ · i,k )

(22)

where ∆ is a tuning parameter that is the same for all models in the model set. The AMM algorithm is used over newer MM algorithms such as Generalized Pseudo Bayesian (GPB) and Interacting Multiple Models (IMM) for the following reasons: 1. AMM assumes the true system does not change over time while GPB/IMM assume that the true system changes over time according to a Markov or semi-Markov chain. In our case the true system parameters and structure are constant with respect to time. 2. GPB/IMM fuse the estimates obtained from the different models which is only possible if all the models share the same states. If assumptions A1-A5 are satisfied there are two possible options when comparing model probabilities: 1. There is ONE clear winner and no output bias. This indicates that the mismatch has been identified correctly and occurs when the number of mismatched parameters (nmis ) is equal to the number of added states (r). 2. There are MANY winners and no output bias. This occurs when the number of mismatched parameters (nmis ) is less than the number of added states (r). All models 20


Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


which contain the mismatched parameters will have similar (or equal) probabilities. If any of the assumptions are violated, there are two possible options regardless of the number of winning models: 1. There is output bias in the winning model(s). This indicates that there is no augmented model that fully captures the mismatch and there will be bias in the state estimates. 2. There is no output bias in the winning model(s). This does not provide any meaningful information as there is no guarantee that state estimates are not biased or that the mismatched parameter(s) have been identified correctly. For our two CSTR example, let us consider the case where the actual T0 = 320 K while the value used in the model is T0 = 300 K. Figure 4 shows the trajectory of the actual states and the estimates given by EKFs using a model augmented with T0 and a model augmented with T03 . The initial state for both EKFs is the nominal steady state and nominal parameter values (T0 = 300 K). Assumptions A1-A5 are satisfied thus it can be seen that only the correct model (T0 ) removes both output and estimation bias. Figure 5 shows the trajectory of the T0 and T03 model probabilities obtained using the proposed algorithm. It can be seen that the correct model is picked. Remark 4 Note that in this work, we consider systems that can be described by pure ordinary differential equations (ODEs). This is also one of the fundamental requirements of classical EKF, MHE estimators. For many systems described by differential and algebraic equations (DAEs), they can be converted to pure ordinary differential equations. For example, by taking time derivatives of the algebraic equations. If the DAEs can be converted to ODEs, the proposed approach can still be used. For DAEs that cannot be converted to ODEs, it is not clear whether the proposed approach works or not. This remains as one of future research topics.

21



306 304 308.8 308.6 308.4 308.2 308

290

x 3 (K)

x 1 (K)

300

0

1

306.1

302

306 300

3.65

280

3.7

2

3.75

305.9 3.65

298

3

4

0

1

2

Time (h)

3.7

3.75 3

4

Time (h) 3

x 4 (kmol/m 3 )

2.8

x 2 (kmol/m 3 )

2.6 2.4 2.2

Actual T 0 model

2.8

T 03 model

2.6 2.4

2 0

1

2

3

4

0

1

2

Time (h)

3

4

Time (h)

Figure 4: Trajectory of the actual states and estimates given by an EKF using a model augmented with T0 and a model augmented with T03 with model plant mismatch in T0 .

0.9

T 0 model

0.8

T 03 model

0.7

Probabilities

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 36

0.6 0.5 0.4 0.3 0.2 0.1 0

0.5

1

1.5

2

2.5

3

3.5

4

Time (h)

Figure 5: Trajectory of T0 and T03 model probabilities obtained using the proposed algorithm with model plant mismatch in T0 .

22


Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


5

Application to a froth flotation system

The proposed algorithm is tested using a froth flotation system commonly used in Coal Handling and Preparation Plants (CHPPs). A schematic diagram of a froth flotation system consisting of five tanks is provided in Figure 6. This diagram was created based on the description of the froth flotation process provided in Ref. 23. Feed slurry enters the first tank where it is mixed with the input reagents. Froth is removed from the top of each tank while tailings are removed from the last tank. For simulation purposes, the actual plant is assumed to have the form: 23 V˙ ufi−1 V˙ ufi dcsti = (csti−1 ) − (csti ) − ri dt Vi Vi V˙ ufi−1 dclti V˙ ufi βi = (clti−1 ) − (clti ) − dt Vi Vi Vi ˙ ˙ Vufi−1 Vufi A˙ i dcati = (cati−1 ) − (cati ) − dt Vi Vi Vi

(23a) (23b) (23c)

where the subscript i = 1, . . . , 5 indicates the tank number, csti is the solids concentration (kg/m3 ), clti is the liquids concentration (kg/m3 ), cati is the ash concentration (kg/m3 ), V˙ ufi is the volumetric rate of the underflow (m3 /min) and Vi is the slurry volume (m3 ). ri is the rate of solid removal as defined below, βi is the mass flow rate of liquid to overflow (kg/min) and A˙ i is the mass flow rate of ash to overflow (kg/min). cst0 , clt0 , cat0 are the solids, liquids and ash concentration of the feed respectively. V˙ uf0 is the volumetric flow rate of the feed.

Figure 6: Schematic diagram of a froth floatation unit

23



Page 24 of 36

The rate of solid removal ri is given by: 23

ri = fr k(csti − c∞ )

(24)

where k is the rate constant (l/min), c∞ is the equilibrium solids concentration (kg/m3 ) and fr is the correction factor for industrial scale reactions. The values of k and c∞ depend on the input reagent concentrations. Additional equations are defined for each stage i, i = 1, . . . , 5, as follows: 23

V˙ ufi−1 = V˙ ufi + V˙ ofi

(25)

β ri Vi V˙ ofi = + ρl ρc ! i−1 i X X A˙ i = xAi M˙ sofj − A˙ j j=1

Ri =

(26) (27)

j=1

csofi V˙ ofi csofi V˙ ofi + csti V˙ ufi

! (100 − Ri−1 ) + Ri−1

xAi = g(Ri )

(28) (29)

where ρl is the liquid density (kg/m3 ), ρc is the coal density (kg/m3 ), Ri is the cumulative solid recovery (%) at stage i, xAi is the cumulative mass fraction of ash in the overflow solids at stage i, V˙ ofi is the volumetric flow rate of the overflow (m3 ), M˙ sofi is the mass flow rate of solids in the overflow (kg/min) and csofi is the concentration of solids in the overflow (kg/m3 ). g(Ri ) is an empirical function of Ri for a given input loading that is obtained from Ref. 23. The nominal model parameters are shown in Table 3. Parameter values are obtained from Ref. 23. The volumetric flow of feed liquid to the first tank is constant at 20 m3 /min while the volumetric flow of feed solids (and total feed) depends on the mass fraction of solids in the feed (cst0 ). 24


Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 3: Nominal parameters for the froth flotation system. k

= c∞ = fr = ρl = ρs = β = V = cst0 = clt0 = cat0 =

3.81 l/min 17.91 kg/m3 0.2 1000 kg/m3 1299 kg/m3 290.3 kg/min 18 m3 10 wt% of feed 90 wt% of feed 24.5 wt% of solids

The process model can be written in the following compact form:

x(t) ˙ = f (x(t), u(t), w(t))

(30)

y = Cx

(31)

T where the state vector is x(t) = x1 (t)T , x2 (t)T , x3 (t)T , x4 (t)T , x5 (t)T with xi (t)T = [ csti (t), clti (t), cati (t)] for i = 1, . . . , 5 and w(t) denotes random process noise. The liquid and ash concentrations in the last tank (clt5 and cat5 ) are measured resulting in: 



 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0  C=  0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

(32)

We select S = {solids fraction in feed (cst0 ), β, k, c∞ , ash fraction of solids in feed (cat0 )}. Solids fraction in the feed is selected because it is not controlled and can vary between 0 − 30%. 23 β is selected arbitrarily and mismatch in β could arise as a results of problems with the motors driving the froth removal paddles. k and c∞ are selected because they are calculated by extrapolating from batch data and depend on the coal type thus are highly likely to contain mismatch. Ash fraction of solids in the feed is selected because it depends on the coal type and can vary between samples of the same coal type. Sensitivity analysis was carried out on all the identified parameters and the results are summarized in Table 4.

25



Step tests were carried out by increasing parameters by 20% from their nominal values one at a time. Gramian tests were conducted by changing parameters one at a time by ±10% and ±20% of their nominal values. Table 4: Results of parameter sensitivity analysis using step tests and empirical observability gramians for froth flotation

cst0 β k c∞ cat0

Step test Empirical Gramian (Sum of absolute normalized gains) (Normalized eigenvalues) 0.619 0.010 0.059 1.105 × 10−5 0.348 0.001 0.324 0.003 1.478 1

Based on these results, we removed β from S resulting in S = {solids fraction in feed (cst0 ), k, c∞ , ash fraction in feed (cat0 )}. Models were augmented with one parameter resulting in four augmented models: cst0 model, k model, c∞ model and cat0 model. Simulations were carried out to test the efficacy of the algorithm using these models with model plant mismatch introduced by setting the actual cat0 to 31.85% (1.3 times the nominal value of 24.5%). Gaussian measurement noise was generated as v v [N (0, 1), N (0, 1)]T . Gaussian process noise was generated for each tank i, i = 1, . . . , 5, as wi (t) v [N (0, 1), N (0, 1), N (0, 1)]T . Small noise values were used to avoid violating assumption A3. The probability threshold δ was set as δ = 0.001. Figure 7 shows the trajectory of the actual states and the estimates given by an EKF using a k model, c∞ model and cat0 model for the last tank (tank 5) with model plant mismatch in cat0 . Estimates from the cst0 model are not shown because the cst0 model is not observable along the entire trajectory. Figure 8 shows the trajectory of the cst0 model, k model, c∞ model and cat0 model probabilities using the proposed algorithm with model plant mismatch in cat0 . Although the wrong model (c∞ model) is picked initially, the algorithm eventually picks the correct model (cat0 model). This is due to the initial lack of separation between the c∞ and cat0 models based on noisy output residuals. The separation between

26


Page 26 of 36

Page 27 of 36

these two models increases as the plant approaches steady state thus the correct model is picked when the plant is closer to steady state. It can also be seen that only the correct model (cat0 ) model is capable of removing bias in the cst estimates.

c st (kg/m3)

30 29 28 27 26 25 0

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

Time (min)

c lt (kg/m3)

981 980 979 978 977 0

1

2

3

4

5

Time (min)

Actual k model c model

28

c at (kg/m3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


26 24

c at model 0

22 20 0

1

2

3

4

5

Time (min)

Figure 7: Trajectory of the actual states and the estimates given by an EKF using a k model, c∞ model and cat0 model for the last tank (tank 5) with model plant mismatch in cat0 . Simulations were also carried out to test the case where some of assumptions A1-A5 were violated. The augmented models, measurement noise and process noise were the same as in the previous case but model plant mismatch was introduced by setting the actual cst0 and actual cat0 to be 1.3 times their nominal values. In this situation assumptions A4 and A5 are violated. Figure 9 shows the trajectory of the actual states and the estimates given by an EKF using a k model, c∞ model and cat0 model for the last tank (tank 5) in this case. Estimates obtained using the cst0 model are not shown as the cst0 model does not converge. Figure 10 shows the trajectory of the cst0 model, k model, c∞ model and cat0 model probabilities using the proposed algorithm in this case. It can be seen that 27



1 0.9 0.8 0.7

Probabilities

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 36

c st model 0

0.6

k model c model

0.5 0.4

c at model 0

0.3 0.2 0.1 0 0

1

2

3

4

5

6

7

8

9

10

Time (min)

Figure 8: Trajectory of the cst0 model, k model, c∞ model and cat0 model probabilities obtained using the proposed algorithm with model plant mismatch in cat0 . although one model (cat0 ) is the clear winner it does not remove bias in the cst estimates. As assumptions A4 and A5 are violated, unbiased output estimates do not guarantee unbiased state estimates.

6

Combined estimation and control of the forth flotation system

The end goal of our MM estimation approach is to use the state estimates obtained with a model based controller such as a MPC. In an ideal case (where assumptions A1-A5 are met), our approach will correctly identify the model which perfectly captures the plant dynamics (removes MPM) and provide unbiased state estimates. This model and associated state estimates can then be used with a MPC to achieve offset free reference tracking. This approach can be considered a model identification (ID) based approach. To test this identification based approach, we carried out state estimation using the proposed MM approach and used the highest probability model at each time step (and

28


Page 29 of 36

c st (kg/m3)

35

30

25 0

1

2

3

4

5

6

7

8

9

10

6

7

8

9

10

6

7

8

9

10

Time (min)

c lt (kg/m3)

980 979 978 977 976 975 0

1

2

3

4

5

Time (min)

c at (kg/m3)

35

Actual k model c model

30

c at model

25

0

20 0

1

2

3

4

5

Time (min)

Figure 9: Trajectory of the actual states and the estimates given by an EKF using a k model, c∞ model and cat0 model for the last tank (tank 5) with model plant mismatch in cat0 and cst0 .

1 0.9 0.8

c st model 0

0.7

Probabilities

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


k model c model

0.6 0.5

c at model

0.4

0

0.3 0.2 0.1 0 0

1

2

3

4

5

6

7

8

9

10

Time (min)

Figure 10: Trajectory of the cst0 model, k model, c∞ model and cat0 model probabilities obtained using the proposed algorithm with model plant mismatch in cat0 and cst0 .

29



Page 30 of 36

associated state estimates) with a MPC of the form: k+N X

min u(tk ),...,u(tk+N )

s.t

(ˇ y (ti ) − yref )T Q(ˇ y (ti ) − yref )

(33a)

i=k

˜ u(t), 0) xˇ˙ (t) = f˜(ˇ x(t), θ,

(33b)

yˇ(t) = g˜(ˇ x)

(33c)

xˇ(tk ) = x˜(tk )

(33d)

˜ u(t) ∈ U xˇ(t) ∈ X,

(33e)

˜ is the where N is the horizon length, Q is a weighting matrix on the output deviations, X set of all possible values of x˜ and U is the set of all possible values of u. yref is the reference signal. xˇ is the predicted state trajectory based on the initial state x˜(tk ) and model (33b). An alternate approach is to use an offset free MPC as presented in Ref. 24. In this approach, an arbitrary disturbance model (subject to certain conditions) is used to capture the effects of the mismatch. State and disturbance estimates are obtained at each time step with the disturbance estimates used to calculate new MPC targets based on the reference signal. The MPC then calculates an input sequence using these new targets. This method does not require the disturbance model to capture the mismatch dynamics perfectly (assumptions A1-A5 do not need to be met) to achieve offset free reference tracking but may result in worse performance than the identification based approach. This is because the parameter (mismatch) estimation is structured in the identification based approach while the offset free MPC treats the parameter mismatch as an arbitrary disturbance. As transient performance is better when the model is closer to the actual plant, 24 the identification based approach may perform better than the offset free MPC. Another disadvantage of the offset free MPC is that the state estimates obtained are not necessarily unbiased. Model parameters used are listed in Table 3. MPM was introduced by setting the plant cst0 = 0.2 instead of the model value of cst0 = 0.1. No process or measurement noise was

30


Page 31 of 36

included (w = v = 0). The model set only contained the cst0 model and cat0 model. The k model and c∞ model were not used as k and c∞ depend on u thus violate assumption A1 in closed loop operation. In order to ensure feasibility of the reference signal, yref was calculated based on the true plant and was set to yref = [968.5, 30.5]T which corresponds to a plant input of u = [0.08, 0.20]T . Figure 11 shows the reference signal and the outputs obtained using an offset free MPC and using the proposed model ID based approach. Figure 12 shows the associated inputs. The model ID based approach achieves offset free reference tracking as assumptions A1-A5 are met. The offset free MPC approach does not achieve offset free reference tracking as some of the conditions for the disturbance model (from Ref. 24) are not met.

982

Reference Offset free MPC Model ID based approach

y1 (kg/m3)

980 978 976 974 972 970 0

5

10

15

20

25

30

20

25

30

Time (min)

50 45

y2 (kg/m3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


40 35 30 25 20 15 0

5

10

15

Time (min)

Figure 11: Trajectory of the outputs using an offset free MPC and using a model ID based approach in the presence of mismatch in cst0

31


Page 32 of 36

0.2 0.15 0.1

0

5

10

15

20

25

30

25

30

Time (min) u2 (g/kg dry feed)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

u1 (g/kg dry feed)


0.4

Offset free MPC Model ID based approach

0.3 0.2 0.1 0

5

10

15

20

Time (min)

Figure 12: Trajectory of the inputs using an offset free MPC and using a model ID based approach in the presence of mismatch in cst0

7

Conclusion

In this paper, a method for obtaining unbiased state estimates in the presence of parameter mismatch was detailed. Necessary assumptions for the success of the method were provided. Proof was given for the simple case of deterministic linear systems operating at steady state. The algorithm was shown to work on a two CSTR system and the froth flotation system when the necessary assumptions were met. It was also shown that offset free reference tracking could be achieved by using the highest probability model (and associated estimates) identified by our proposed estimation approach with a conventional MPC.

References (1) Muske, K. R.; Badgwell, T. A. Disturbance modeling for offset-free linear model predictive control. Journal of Process Control 2002, 12, 617–632. (2) Haverbeke, N.; Van Herpe, T.; Diehl, M.; Van den Berghe, G.; De Moor, B. Nonlinear model predictive control with moving horizon esitmation - Application to the normal32


Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ization of blood glucose in the critically ill. Proceedings of the 17th World Congress, The International Federation of Automatic Control, Seoul, Korea, 2008. (3) Jacob, N. C.; Dhib, R. Unscented Kalman filter based nonlinear model predictive control of a LDPE autoclave reactor. Journal of Process Control, 2011, 21, 1332-1344. (4) Lee, J. H.; Ricker, N. L. Extended Kalman filter based nonlinear model predictive control. Industrial & Engineering Chemistry Research, 1994, 33, 1530-1541. (5) Botelho, V.; Trierweiler, J. O.; Farenzena, M.; Duraiski, R. Methodology for Detecting Model Plant Mismatches Affecting Model Predictive Control Performance. Industrial & Engineering Chemistry Research 2015, 54, 12072–12085. (6) Parlos, A. G.; Menon, S. K.; Atiya, A. F. An Adaptive State Filtering Algorithm for Systems With Partially Known Dynamics. Journal of Dynamic Systems, Measurement, and Control 2002, 124, 364. (7) Salhi, H.; Bouani, F. Nonlinear Parameters and State Estimation for Adaptive Nonlinear Model Predictive Control Design. Journal of Dynamic Systems, Measurement, and Control 2016, 138, 044502–044502. (8) Lee, D. J.; Park, Y.; Park, Y. s. Robust H∞ Sliding Mode Descriptor Observer for Fault and Output Disturbance Estimation of Uncertain Systems. IEEE Transactions on Automatic Control 2012, 57, 2928–2934. (9) Yang, Q.; Liu, Y. Adaptive state estimation of multi-input and multi-output non-linear systems with general uncertainties both in the state and output equations. IET Control Theory Applications 2016, 10, 354–362. (10) Semerdjiev, E.; Mihaylova, L.; Li, X. R. Variable- and fixed-structure augmented ZMM algorithms using coordinated turn model. Proceedings of the Third International Conference on Information Fusion FUSION 2000. 2000; pp MOD2/25–MOD2/32. 33



(11) Jiao, Z.; Li, W.; Zhang, Q.; Wang, P. Tracking for the maneuvering target based on multiple model and moving horizon estimation. 2015 IEEE International Conference on Information and Automation. 2015; pp 1936–1941 (12) Pitre, R. R.; Jilkov, V. P.; Li, X. R. A comparative study of multiple-model algorithms for maneuvering target tracking. Proceedings of the SPIE. 2005; pp 549–560. (13) Amirzadeh, A.; Karimpour, A.; Moeini, A. An IMM algorithm based on augmented Kalman filter for maneuvering target tracking. Scientific Research and Essays 2011, 6, 6787–6797. (14) Huang, D.; Leung, H. An expectation-maximization-based interacting multiple model approach for cooperative driving systems. IEEE Transactions on Intelligent Transportation Systems 2005, 6, 206–228. (15) Zhang, Y.; Li, X. R. Detection and diagnosis of sensor and actuator failures using IMM estimator. IEEE Transactions on Aerospace and Electronic Systems 1998, 34, 1293–1313. (16) Chetouani, Y. Design of a multi-model observer-based estimator for Fault Detection and Isolation (FDI) strategy: application to a chemical reactor. Brazilian Journal of Chemical Engineering 2008, 25, 777–788. (17) Kuure-Kinsey, M.; Bequette, B. W. Multiple Model Predictive Control Strategy for Disturbance Rejection. Industrial & Engineering Chemistry Research 2010, 49, 7983– 7989. (18) Li, X. R. Model-set design for multiple-model method. Part I. Proceedings of the Fifth International Conference on Information Fusion. 2002; pp 26–33. (19) Sun, Y.; El-Farra, N. H. Quasi-decentralized model-based networked control of process systems. Computers & Chemical Engineering 2008, 32, 2016–2029. 34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(20) Geffen, D. Parameter identifiability of biochemical reaction networks in systems biology. M.Sc. thesis, Queen’s University, Kingston, Ontario, Canada, 2008. (21) Marino, R.; Tomei, P. Nonlinear Control Design: Geometric, Adaptive, and Robust; Prentice-Hall, 1995. (22) Pitre, R. A Comparison of Multiple-Model Target Tracking Algorithms. M.Sc. thesis, University of New Orleans, New Orleans, Louisiana, United States, 2004. (23) Canright, G. S.; Brown, C. H.; Allgood, G. O.; Hamel, W. R. Dynamic Modeling and Control Analysis of Froth Floatation and Clean Coal Filtration as applied to Coal Benefictiation; 1981. (24) Morari, M.; Maeder, U. Nonlinear offset-free model predictive control. Automatica 2012, 48, 2059–2067.

35



Graphical TOC Entry

36


Page 36 of 36

Handling Model Plant Mismatch in State Estimation Using a Multiple

Recommend Documents