System identification and control design of a nonlinear continuously stirred tank reactor

Journal Pre-proof System identification and control design of a nonlinear continuously stirred tank reactor Abolfazl Simorgh, Abolhassan Razmini, Vladimir I. Shiryaev

PII: DOI: Reference:

S0378-4754(20)30017-3 https://doi.org/10.1016/j.matcom.2020.01.010 MATCOM 4934

To appear in:

Mathematics and Computers in Simulation

Received date : 27 November 2018 Revised date : 24 December 2019 Accepted date : 20 January 2020 Please cite this article as: A. Simorgh, A. Razmini and V.I. Shiryaev, System identification and control design of a nonlinear continuously stirred tank reactor, Mathematics and Computers in Simulation (2020), doi: https://doi.org/10.1016/j.matcom.2020.01.010. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. c 2020 Published by Elsevier B.V. on behalf of International Association for Mathematics and ⃝ Computers in Simulation (IMACS).

Journal Pre-proof

System Identification and Control Design of a Nonlinear Continuously Stirred Tank Reactor

a

Abolfazl Simorgha

Abolhassan Razminiaa∗

Vladimir I. Shiryaevb

[email protected]

[email protected]

[email protected]

Dynamical Systems & Control (DSC) Research Lab., Department of Electrical Engineering, School of Engineering, Persian Gulf University, 75169, Bushehr, Iran Department of Automatic Control Systems, South Ural State University, Chelyabinsk, 454080, Russia

roo

f

b

Abstract

lP

rep

This study presents a black-box identification and a self-tuning control design of a nonlinear continuously stirred tank reactor (CSTR). Analyzing and controlling the CSTR is a challenging topic because of its highly nonlinear dynamics and its sensitivity, which, in turn, requires a high-performance control design. The data set used in the identification is collected from the nonlinear dynamical equations of the CSTR in the presence of white and colored noises for considering realistic conditions. linear model structures are selected, and by employing the online and offline estimation methods, the best representation of the CSTR is determined based on a fitness index. By employing the results of the online identification algorithms, a self-tuning proportional-integral-derivative (PID) controller is designed which provides an ability to control the CSTR in the presence of unwanted signals (i.e., noise and disturbance) and variations in system dynamics by considering several tuning parameters including forgetting factor and suitable initial values. The performance of the proposed approach is compared with several works in the literature, reflecting the superiority of the presented technique.

rna

Keywords: Continuous stirred tank reactor; System identification; PID controller; Online estimation; Forgetting factor.

1 Introduction

Jo u

Continuously stirred tank reactor (CSTR) is one of the main components in the industry with numerous applications [1, 2]. The CSTR plays a crucial role in maximizing the productivity in chemical processes, however, analyzing and controlling such systems is a challenging task due to the complex behavior referring to their highly nonlinear dynamics and sensitivity [1–5]. To predict the behavior of a system, a model description is required, which can be obtained using physical laws such as the second Newton law and mass balance equation or determining the model experimentally using system identification (SI) techniques [6–8]. The dynamical equations of a system derived from physical laws mostly result in highly nonlinear equations which are difficult to be employed to design a controller such as the CSTR [9–11]. Recently, numerous approaches have been suggested in the literature for controlling the CSTR, including the sliding mode control [1, 9, 12], the feedback linearization [13], the neural networks [5, 14–16], the PID control [10, 17], and the model predictive control [18]. An important issue in ∗ Corresponding

Author. Email: [email protected], Tel.: +98 773 1222164.

1

Journal Pre-proof

rep

roo

f

such control strategies is robust and adaptive performances in the presence of model uncertainties in the model-based controllers and varying environmental conditions for all techniques. Due to the presence of uncertainty in the dynamical equations of the CSTR [18–20], a datadriven identification and control design is proposed in this study to analyze and control the system. The data-set is collected in the presence of white and colored noises from the nonlinear dynamical equations of the CSTR. To predict the behavior of the CSTR with the assumption of no information of its structure, the black-box identification techniques are employed. In the black-box identification, the dynamical model of systems is assumed to be unknown, and only input/output data are available [21] (see [22–27] for the applications). It should be noted that the dynamical equations of CSTR are used only for collecting data. In this manuscript, linear model structures are selected, and by employing the online and offline estimation methods, the best representation of the CSTR is determined by means of a fitness index. The PID is known as one of the most efficient controllers in industrial applications [28–30]. For instance, the design of PID controllers for the CSTR using the meta-heuristic methods such as modified particle swarm optimization (PSO) and artificial bee colony (ABC) have been presented in [10, 17], respectively. These controllers are designed offline by the assumption that the system is time-invariant. However, in practice, owing to the presence of uncertainties, noises, and disturbances, the offline methods may fail to control the system, and can cause instability. In this paper, based on the online estimation methods, such as the recursive least-square and its modification with forgetting factor, a self-tuning PID controller using Ziegler-Nichols tuning algorithm is designed, providing an ability to control the CSTR in the presence of unwanted signals (i.e., noise and disturbance) and variations in system dynamics. The applicability of the proposed method is supported by comparing it with several works in the literature. This paper is organized as follows. The data-driven identification of the CSTR is presented in Section 2, where the simulation results are provided in Subsection 2.3. Section 3 deals with designing an adaptive PID controller based on the recursive identification algorithm, whose simulations are performed in a subsection. In the end, some concluding remarks are provided in Section 4.

2 Identification of CSTR

rna

lP

In this section, a black-box identification of a CSTR is studied by employing linear model structures and estimation methods. In the black-box identification, only a set of input/output data is available, and the structure of the system is unknown. The data used in this study is collected from the following nonlinear representation of CSTR [10, 17]

x˙ 1 (t) = −x1 (t) + Da (1 − x1 (t)) exp(

x2 (t) ) 1 + x2 (t)φ−1

x˙ 2 (t) = −(1 + δ)x2 (t) + B · Da (1 − x2 (t)) exp(

x2 (t) ) + δu(t) 1 + x2 (t)φ−1

(1)

Jo u

where x1 (t) and x2 (t) are the dimensionless concentration of the reactant and temperature of the reactor, respectively, and u(t) is the dimensionless temperature of the cooling jacket. x1 (t) is the output of the system, and u(t) is the controlling input. There are some parameters in the Eq. (1) including φ, δ, Da , and B that refer to the activated energy, heat transfer coefficient, Damokhler number, and heat of reaction, respectively. As it is known, the model descriptions based on the physical laws such as the mass balance equations represent systems behavior satisfactorily, however, they mostly result in highly nonlinear representations, in turn, challenging to treat. Analyzing and controlling the considered CSTR is a challenging task due to the nonlinearity and unstable equilibrium points [10]. In this paper, a study on the identification of the considered CSTR is proposed, and based on the results, a controller will be designed. The nonlinear dynamical equation of the CSTR is employed here for collecting data in the presence of noises. To this end, a stochastic term v(t) is added to the left-hand side of Eq. (1) for considering realistic conditions such as measurement 2

Journal Pre-proof

Table 1: Some applicable linear model structures. Name

Structure

Autoregressive with exogenous input (ARX)

B(q −1 ) 1 A(q −1 ) u(t) + A(q −1 ) e(t) −1 ) C(q −1 ) y(t) = B(q A(q −1 ) u(t) + A(q −1 ) e(t) −1 ) y(t) = B(q F (q −1 ) u(t) + e(t) −1 D(q −1 ) ) y(t) = B(q F (q −1 ) u(t) + C(q −1 ) e(t)

y(t) =

Autoregressive moving average exogenous (ARMAX) Output error (OE) Box-Jenkins (BJ)

noises. To have a reliable identification, the CSTR must be excited perfectly. According to the literature, the pseudo-random binary sequence (PRBS) is a suitable signal to provide an informative data set [6], that is defined as follows { } ( ) ( ) ( ) DN = u(1), y(1) , u(2), y(2) , · · · , u(N ), y(N )

roo

f

where N is the number of samples. In the next step, a model that fits the gathered data well should be selected. In this manuscript, the attention is restricted to the standard linear models, taking the following form (2)

y(t) = P (q, θ)u(t) + Π(q, θ)e(t)

P (q, θ) =

B(q −1 ) , A(q −1 )F (q −1 )

Π(q, θ) =

rep

where

D(q −1 ) , A(q −1 )C(q −1 )

q −1 is the shift operator (i.e., q −k x(t) = x(t − k)), θ is the parameters of the model, e(t) is a zero mean white noise and A(q −1 ) = 1 +

na ∑

ai q −i ,

lP

i=1

D(q −1 ) = 1 +

nd ∑

di q −i ,

B(q −1 ) =

nb ∑

bi q −i ,

C(q −1 ) = 1 +

i=1

F (q −1 ) = 1 +

i=1

nc ∑

ci q −i

i=1

nf ∑

fi q −i .

i=1

rna

Some applicable special cases of Eq. (2) are given in Table 1. After choosing a model structure, its parameters must be determined such that the data set be represented by the selected model acceptably. To this end, an algorithm based on the optimization needs to be conducted. The one-step-ahead predictor using Eq. (2) can be calculated as follows

Jo u

yˆ(t|θ) = [1 − Π−1 (q, θ)]y(t) + Π−1 (q, θ)P (q, θ)u(t)

(3)

where yˆ(t|θ) is the output of the selected model structure calculated from data set DN and the type of the model. The difference between the actual and predicted output is then defined as the prediction error, that is ε(t, θ) = y(t) − yˆ(t|θ), t = 1, 2, · · · , N. The model represents the system perfectly if the prediction error gets minimized for all samples. Prediction error method (PEM) is an identification technique that tries to find the parameters of the model such that a criterion based on the prediction error (ε(t)) gets minimized. The criterion considered in the PEM has the following general form VN (θ) =

N 1 ∑ Γ(ε(t, θ)) N t=1

3

Journal Pre-proof

where Γ(·) is a scalar-valued function. The optimal parameters θˆ are calculated from θˆP EM = arg min{VN (θ)}. θ

by means of offline and online algorithms.

2.1

Offline identification

In this case, the data has to be gathered first, and then the identification algorithm is implementable. A particular case of the PEM is the least square (LS) with the following criterion N

VN (θ) =

1∑ 2 ε (t, θ) 2 t=1

(4)

where ε(t) = y(t)− yˆ(t). If the ARX structure be chosen as the model, the one-step-ahead predictor yˆ(t, θ) can be written in the linear regression form as yˆ(t|θ) = ϕT (t)θ.

(5)

θ

roo

f

In this case, the estimated parameters will be received as [6] ] ]−1 [ ∑ [∑ N N T ˆ ϕ(t)y(t) ϕ(t)ϕ (t) . θLS (t) = arg min{VN (θ)} = t=1

(6)

t=1

rep

where θ ∈ Rn×1 is the vector of parameters to be estimated and ϕ(t) ∈ Rn×1 is a vector that contains the previous values of the output and input signals. Let assume that the data set has been collected from y(t) = ϕT (t)θ0 + v(t) (7)

lP

where θ0 is the exact parameters, and v(t) is a stochastic term that can represent the white and colored noises. Let v(t) be a white noise (i.e., have zero mean), then by substituting Eq. (7) in Eq. (6), one can yield N 1 ∑ −1 θˆLS(N ) − θ0 = KN · ϕ(t)v(t) (8) N t=1 ∑N where KN = N1 t=1 ϕ(t)ϕT (t). By increasing the collected data (N → ∞), we have E{θˆLS(N ) } − θ0 = E{ϕ(t)ϕT (t)}−1 · E{ϕ(t)v(t)}

(9)

Jo u

rna

where E{·} is the expectation operator. The stochastic term in Eq. (7) is considered as a white noise, which implies that ϕ(t) and v(t) are uncorrelated because v(t) does not rely on what has occurred to time t − 1, thus E{ϕ(t)v(t)} = 0. Consequently, to have an unbiased estimation i.e., E{θˆLS } = θ0 , E{ϕ(t)ϕT (t)} must be non-singular. For the colored noise, to have an unbiased estimation, the instrumental variable (IV) technique has been suggested. Since the ARX model structure is used, the regressor vector ϕ(t) contains the delayed outputs and inputs as [ ]T ϕ(t) = −y(t − 1) · · · −y(t − na ) u(t − 1) · · · u(t − nb ) which yields

E{ϕ(t)v(t)} ∝ E{y(t − 1)v(t)} ∝ E{v(t − 1)v(t)} ̸= 0

thus, the presence of colored noise leads to biased estimation. To handle this drawback, the delayed outputs in the regressor vector must be excluded and replaced with suitable variables known as IV, resulting in the following estimation ]−1 [ ∑ ] [ ∑ N N 1 1 · ζ(t)ϕT (t) ζ(t)v(t) (10) θˆIV (t) − θ0 = N t=1 N t=1 4

Journal Pre-proof

where ζ(t) is an instrumental variable. In this paper, by filtering the input signal ∆(q −1 )g(t) = Ω(q −1 )u(t) and taking the na delayed values of the resulting filtered input, the following IV is received [ ]T ζ(t) = −g(t − 1) −g(t − 2) · · · −g(t − na ) u(t − 1) · · · u(t − nb ) , yielding E{ζ(t)v(t)} = 0. The identification using the other linear models given in Table (1) results in a biased estimation when the LS method is employed because the one-step predictor in these cases can not be written in the linear regression form [6]. Two alternatives to identify such models are to use the extended LS or employing the statistical framework for identification, such as the maximum likelihood (ML) that has been presented in [6, 7].

2.2 Online identification

f

In the online identification, the parameters are estimated recursively as the data is received, which significantly reduces the computational time and can be employed by self-tuning controllers for controlling the uncertain dynamics [31]. Let the data set be generated from Eq. (1), and the criterion be considered as N 1 ∑ N −t 2 λ ε (t, θ), N t=1

roo

VN (θ) =

λ ∈ (0, 1]

(11)

rep

where λ is the forgetting factor parameter. If λ is set to one, Eq. (11) has the form as the classical least square, however, including such parameter in the criterion, leads to some valuable features ∑t T that will be discussed in follows. By defining P (t) := ϕ(i)ϕ (i) and using Eq. (6), the i=1 following recursive estimation of parameters at step t is obtained θˆt = θˆt−1 + P −1 (t)ζ(t).ε(t)

(12)

lP

where the recursion for P −1 using the Sherman–Morrison formula can be written as follows [ ] P −1 (t − 1)ζ(t)ϕT (t)P −1 (t − 1) 1 −1 −1 P (t − 1) − , P (t) = (13) λ λ + ϕT (t)P −1 (t − 1)ζ(t)

rna

handling the drawbacks of calculating the inverse of P (t) at each step, in turn, reducing the computational time. The following applicable special cases are considered and summarized as follows. • Recursive LS (RLS): Set λ = 1 and replace ζ(t) with ϕ(t). This is the same as offline LS, however, it is solved recursively, allowing to identify timevarying systems, and can be employed by self-tuning controllers.

Jo u

• Recursive IV (RIV): Set λ = 1. This case is the modification of RLS that provides an ability to have an unbiased estimation in the presence of colored noises. • Recursive LS with exponential forgetting factor (RLS-FF): Eq. (11). The RLS-FF uses the generalized LS criterion Eq. (11) that can efficiently identify the time-varying systems by means of forgetting factor. The selection of some parameters in the recursive identification algorithm, including the initial value of parameters θ(0), the initial value of P −1 (0), and the forgetting factor coefficient affects the result, which will be discussed in the next part.

5

Journal Pre-proof

2.3

Numerical simulations

In this part, the identification of the CSTR is carried out based on the online and offline estimation methods, and the discussion on results is given as well. The data is gathered from the nonlinear dynamical equation of CSTR Eq. (1) with the initial condition (x1 (0), x2 (0)) = (0.3, 0.3). Notice that the nominal value of parameters used in Eq. (1) are considered as follows [10]: Da = 0.072, φ = 20, B = 8, and δ = 0.3. As the first step, standard linear models provided in Table (1) are employed by offline estimators, then, the comparison between them is conducted based on the following fitness index

) (

y − y ˆ 2

Fitness = 100 1 −

y − y˜Υ 2

rep

roo

f

ˆ ∈ RN , y ∈ RN contain the samples of the predicted and actual outputs, where the vectors y respectively, y˜ is the average of actual output’s samples, and Υ ∈ RN is a vector of ones. To validate the proposed identification, the data set is divided into two subsets: train and test. The train set is used for identifying parameters while the test set validates the performance of the identified system. In this paper, N = 1000 samples are collected in 50 seconds, where the first 35 seconds are labeled as the train data, and the rest samples are labeled as the test. In the first step, the ARX model is selected with different orders, and the results are depicted in Fig. (1), and Table (2). In this case, the best and the worst fitness have been received 94.4920%, 68.12% for the train set and 86.1367%, 5.42% for the test set, respectively which have been obtained by selecting na = 3, nb = 2 and na = 2, nb = 1. As results show, the ARX model with higher-orders have provided better fitness than the other lower-order ones. In fact, the higher-order models are more flexible to track the CSTR output.

0.45 0.4 0.35 0.3 0.25

0.1 0

5

10

15

20

25

30

35

40

45

50

15

20

25

30

35

40

45

50

rna

0.45

lP

0.2 0.15

0.4 0.35 0.3 0.25 0.2

Jo u

0.15

0.1

0

5

10

Figure 1: Offline identification of CSTR using LS for the ARX model structure.

In the next cases, the ARMAX, BJ, and OE models with different orders are employed, and the results using the ML method to minimize the criterion Eq. (4) are depicted in Fig. (2), and Table (2). All these models, with the orders of two or more have successfully tracked the CSTR output, however, the BJ has the best fitness when the order nb = 3, nf = 3, nc = 3, and nd = 3 are selected, which is 97.52% for the train set and 96.6% for the test set. According to Table (1), 6

Journal Pre-proof

the structure of the BJ model is more flexible due to the presence of polynomials B(q −1 ), F (q −1 ), C(q −1 ), and D(q −1 ). Consequently, it can result in better fitness than the other linear models as has occurred.

0.32

0.4

0.31 0.3

0.3

0.2

0.29 0.28

0.1 0

5

10

15

20

25

30

35

35

40

45

50

40

45

50

40

45

50

0.32

0.4

0.31 0.3

0.3

0.2

0.29 0.28

0.1 0

5

10

15

20

25

30

35

35

0.32

f

0.4

roo

0.31

0.3

0.3

0.2

0.29 0.28

0.1 5

10

15

20

25

30

35

35

rep

0

Figure 2: Offline identification of CSTR using ML for the ARMAX, OE, and BJ model structures.

rna

lP

The previous simulations with suitable structures and orders have precisely fitted (i.e., more than 94%) the CSTR. In fact, the data are collected in a specific operating condition, and the obtained results may not be valid in the other operating conditions (e.g., when the amplitude of the control signal changes). In such a situation, the recursive algorithms are helpful since they can adapt themselves to new conditions. In this case, the second-order ARX model is employed, and the results using RLS are provided in Fig. (3). There are some tuning parameters in the RLS algorithm, including P −1 (0) and θ0 that can affect the quality of the identification. As can be observed in Fig. (3), when P −1 (0) = 105 is selected, the best fitness has been obtained and significant variations in the parameters can be seen compared to the two other values of P −1 (0). The reason can be explained as follows. The update law of parameters in the RLS method is

Jo u

θˆt = θˆt−1 + P −1 (t)ϕ(t).(y(t) − ϕT θ(t − 1))

(14)

where at the beginning of the run, ε(t) = y(t) − ϕT θ(t − 1) is large because the value of parameters are unknown (i.e., θ0 is chosen randomly). Therefore, when a significant value is assigned to P −1 (0), the value of the parameters in the next iteration will drift away from the previous values and tries to reduce the prediction error as soon as possible to decrease the effect of large value of P −1 (t) in P −1 (t)ϕ(t)ε(t). Since the data set is the same as the one used in the offline estimation algorithm, the parameters a1 , a2 , b1 , and b2 have converged to constant values due to the constant operating condition. As can be observed, the choice of P −1 (0) affects the convergence rate of the parameters in which, the largest one has resulted in the fastest convergence (especially b1 , b2 ), because the optimizer tries to seek the parameter that reduces the prediction error quickly to neutralize the effect of large assigned weight, i.e. P −1 (0).

7

Journal Pre-proof

0.4 0.3 0.33

0.2

0.325

0.1

0.32 10

0 0

2

4

6

8

10

12

10.5 14

11 16

11.5

12 18

20

1

0 -0.5

0.5

-1 0 -1.5 -2

-0.5 0

5 10

10

10

15

20

0

-3

10

10

5

5

10

15

20

5

10

15

20

-3

5

0 5

10

15

20

0

roo

0

f

0

Figure 3: Recursive identification of CSTR using RLS for the ARX model structure.

lP

rep

In the next case, the dynamic of the CSTR is considered to be time-varying, which is more realistic and practical. The data is collected from Eq. (1) in the presence of white noise and varying parameters as   B = 1, a = 0.3, φ = 20, Da = 0.072 0 ≤ t < 10 B = 4, a = 0.3, φ = 20, Da = 0.072 10 ≤ t < 20  B = 50, a = 1, φ = 50, Da = 0.2 t ≥ 20

Jo u

rna

where the operating conditions change during such variations in parameters. The RLS can be employed to identify time-varying systems, however, It will be shown that the modification of RLS with the forgetting factor can efficiently handle such systems. Consider the generalized LS criterion Eq. (11), where the forgetting coefficient should be selected λ ∈ (0, 1]. If the λ is set to a value less than one, the weight on the prediction error in previous samples in Eq. (11) becomes lower, which in turn results in forgetting previous information more quickly than the RLS. Moreover, in the recursive algorithm presented for calculating P −1 Eq. (13), λ−1 rises approximately exponentially that causes the P −1 to take significant values and as a result, reduces the prediction error quickly, referring to reason discussed in the previous case. The results of including the forgetting factor in the RLS are given in Fig. (4). When λ = 0.9, the update algorithm tries to neglect the previous info more quickly than the case when λ = 0.98 or 1 is selected. Fig. (4) shows the oscillatory behavior in the parameters’ trajectories for λ = 0.9 because, in this case, the model tries to capture the dynamics of the CSTR at each iterate based on the received data, not the previous ones. Therefore, it is a suitable alternative for identifying time-varying systems since it can capture even a small variation.

8

Journal Pre-proof

0.4 0.3 0.32 0.2 0.31 0.1 0.3 9.5

0 0

5

10

15

10

20

10.5 25

30

1

0 -0.5

0.5

-1 0 -1.5 -2

-0.5 0

5 10

10

10

15

20

25

30

0

-3

10

10

5

5

10

15

20

25

30

5

10

15

20

25

30

-3

5

0

0 5

10

15

20

25

30

0

roo

f

0

Figure 4: Recursive identification of CSTR using RLS with exponential forgetting factor for the ARX model structure.

lP

rep

In the last case, the data is gathered in the presence of colored noise, where the results of identification are depicted in Fig. (5). The comparison between the RLS and RIV with the same orders shows that the RIV has performed better than the RLS due to the unbiased identification property that it has in the presence of colored noises. Moreover, By applying the RIV to the third-order model, perfect identification has been received, as can be seen in Fig. (5).

0.4

0.3

0.1

0 0

5

10

0.36

0.34

0.32 20 15

20

25

20.5 30

21 35

21.5

22 40

45

50

45

50

0.02

Jo u

0.2

rna

0.2

0.15

0.01

0.1

0

0.05

-0.01 20

20.5

21

21.5

22

0

-0.05

0

5

10

15

20

25

30

35

40

Figure 5: Recursive identification of CSTR using RIV for the ARX model structure.

9

Journal Pre-proof

3 Control design This paper proposes a control design for the CSTR by employing the results of the previous section. As was previously shown, the ARX, ARMAX, OE, and BJ with suitable orders represented the CSTR successfully in the considered operating points. Therefore, the nonlinear CSTR can be approximated by linear models with slowly time-varying parameters referring to their operating conditions. By using this idea, a controller will be designed in this part. Among the presented models, the ARX is selected to represent the CSTR, which provides an ability to use most of the estimation algorithms, and can be employed to design an adaptive controller efficiently [32]. The controller that is considered to be designed adaptively is the PID type, which is the most commonly used controller in the industry [10, 17]. The PID controllers take the difference between the desired and the actual outputs (r(t) and y(t), respectively) as an input and based on the proportional, integral, and derivative terms try to shape a control policy that minimizes the error continuously. The parallel PID controller has the following structure ∫ [ dε(t) ] 1 t u(t) = Kp ε(t) + ε(τ )dτ + Td (15) Ti 0 dt

lP

rep

roo

f

where u(t) is the controller output, ε(t) = r(t)−y(t) and Kp , Ti , Td are tuning parameters. Despite broad applicability, the PID controllers suffer from not being tuned well. Through numerous techniques in tuning the PID controllers, the Ziegler-Nichols (ZN) methods have been received more attention due to the applicability and simplicity in implementation [33]. The tuning rule in ZN methods is summarized in some tables concerning the system’s characteristics such as the settling time, delay in response, and critical points. If the system is linear and completely known, the PID parameters can be simplicity found according to two suggested ZN methods known as the s-shape and ultimate gain techniques. In the ultimate gain method, the PID parameters are calculated as follows [32] Tu Tu 3Ku Td = , Ti = , Kp = 8 2 5 where Ku and Tu are the critical gain and period, respectively. Let the integral and derivative terms of PID controllers be turned off, then the value of Kp that transfers the closed-loop system to the stability boundary (i.e., when the closed-loop poles are located at imaginary axis), is the critical gain and the period of permanent oscillations in such situation is the critical period. The results in the previous section revealed that the third-order ARX model could successfully capture the CSTR by using the recursive estimation methods. Therefore, it is employed here to design the controller. Let the system be given as

rna

P (q −1 ) =

B(q −1 ) b1 q −1 + b2 q −2 + b3 q −3 = −1 A(q ) 1 + a1 q −1 + a2 q −2 + a3 q −3

(16)

Jo u

where the parameters are time-varying considering the operating condition. To design a PID controller for the system Eq. (16), its contiuous form Eq. (15) must be discretized. There exist various methods for discretizing the continuous-time functions, including zero-order hold and firstorder hold. In this paper, a Tustin transformation is adopted that provides the discretized form of Eq. (15) by replacing s operator in the frequency-domain representation of Eq. (15) with −1 2 1−q T0 1+q −1 [6, 32]. The Laplace transform of Eq. (15) can be simply calculated as follows E(s) 1 = Kp [1 + + Td s]. U (s) Ti s

(17)

The derivative component in Eq. (17) amplifies the process or measurement noises (i.e., highfrequency), resulting in undesirable behaviors in the output. Therefore, a modification of the derivative term with a low-pass filter is desired that is Td s ≈

Td s , τf s + 1 10

τf = ξ −1 Td

Journal Pre-proof

where ξ regulates the time constant of the filter. By using the Tustin transformation, the following digital PID will be received [32] u(k) = γ1 u(k − 1) + γ2 u(k − 2) + β0 ε(t)(k) + β1 ε(t)(k − 1) + β2 ε(t)(k − 2) where 4κf 2κf − 1 1 + 2(κf + κd ) + 0.5κi (1 + 2κf ) , γ2 = , β0 = Kp 1 + 2κf 1 + 2κf 1 + 2κf 0.5κi − 4(κf + κd ) κf (2 − κi ) + 2κd + 0.5κd − 1 β1 = Kp , β2 = Kp 1 + 2κf 1 + 2κf Td τf T0 κi = , κd = , κf = . Ti T0 T0 γ1 =

f

To tune the digital PID using the ZN method, the Ku and Tu must be determined. As mentioned, these parameters are obtained when the system operates in the stability boundary. Such a situation occurs If at least one of the closed-loop poles of the system is located on the unit circle, and the others are inside it. The characteristic polynomials of the closed-loop system with the proportional controller is

roo

C(q −1 , Kp ) = 1 + a1 q −1 + a2 q −2 + a3 q −3 + Kp (b1 q −1 + b2 q −2 + b3 q −3 ) where to find the critical points, two situations can occur: having complex or real poles at -1 [32]. For the first case, we have q 3 C(q −1 , Ku ) = q(q − α − j)(q − α + j)(1 − ζq −1 )

(18)

rep

where the parameters Ku , ζ, and α are unknown and can be determined by equating the coefficients of q, q 2 , and q 3 in both sides Eq. (18), which yields √ −µ1 ± δ α = 0.5[a3 − a1 + Ku (b3 − b1 )], Ku1,2 = 2µ2

lP

where

δ = µ21 − 4µ2 µ0 ;

µ2 =

b23

− b1 b3 ;

µ0 = (a3 − a1 )a3 − 1 + a2 ,

µ1 = b2 − (a1 − 2a3 )b3 − a3 b1 ,

rna

and ζ can be easily calculated using

ζ1 = a3 + b3 Ku1 ,

ζ2 = a3 + b3 Ku2

Jo u

from equating the coefficients of Eq. (18). The obtained parameters must satisfy one of the following conditions concerning the idealized location of closed-loop poles (i.e., critical condition) ζ1 < 1, Ku1 > 0 or ζ2 < 1, Ku2 > 0 (19) where the gain Kui for i = 1, 2 that satisfies the conditions will be selected as the critical gain. In this case, the critical period that depends on the location of poles in the z-plane is calculated as Tu = 2πT0 (cos−1 α)−1 . However, if the conditions Eq. (19) do not satisfy, the real poles must be considered which can be written as q 3 C(q −1 , Ku ) = q 2 (q + 1)(1 + ζ¯1 q −1 + ζ¯2 q −2 )

that by equating the coefficients, Ku can be calculated as follows Ku =

a2 − a3 − a1 + 1 b3 − b2 + b1

and in this case, the critical period is Tu = 2T0 .

11

Journal Pre-proof

Numerical simulations

f

In this part, based on the presented methodology, a self-tuning PID controller is designed for several scenarios. The objective is to find a control policy that causes the outputs of the nonlinear π CSTR Eq. (1) tracks a sinusoidal wave r(t) = 2π 15 sin( 25 t)+0.5 in the presence of various conditions. In all cases, the results are compared with two recently suggested methods for controlling CSTR known as the MPSO [17] and the ABC [10]. In these two approaches, the PID controller has been employed as the controller, and the PID parameters have been determined using meta-heuristics methods in an offline mode. Let consider the CSTR Eq. (1) to be controlled by the proposed PID controller using the RLS method for identifying the model’s parameters Eq. (16) and ZN tuning algorithm. The simulations are carried out and depicted in Fig. (6). As can be seen, the proposed controller has performed better than the two other methods except for the first 5 seconds, which refers to the uncertainty in the initial value of model parameters. However, it has adapted itself quickly and tracked the sinusoidal wave perfectly, thereafter. As was discussed in the identification section, when P −1 (0) is selected large enough, the model will fit the system more accurately and quickly, in turn, the controller will be designed in a more reliable way (see Fig. (6)).

1

roo

0.088 0.086

0.8

0.084 0.082

0.6

0.92

0.08 0.4

36

37

38

39

0.915

0 0

10

20

30

10 2

0.4

0

-3

lP

0.6

rep

0.2

0.2

40

0.91 61

62 60

50

63

64 70

80

90

100

70

80

90

100

-2

20

0

-0.4 0

rna

-0.2

10

20

40

30

60

40

80

100

50

60

Jo u

Figure 6: Online PID control design for tracking a sinusoidal wave.

In the next case, a constant disturbance d = 0.5 is added to the system in the interval 50 ≤ t ≤ 60. The results of THE simulation in Fig. (7) show that the proposed controller has successfully damped the effects of disturbance in comparison with MPSO and ABC method. The reason is that those controllers have been tuned offline in a specific condition and may not keep their performance in the other conditions as the one considered here. Whereas, due to the adaptive properties of the proposed PID design, the controller can efficiently handle various situations.

12

Journal Pre-proof

1 0.9 0.8 0.8 0.6 0.7 0.4 0.6

0.2 0 0

20

0.5

40

60

80

0.5 50

100

55

10 -3

0.4

60

65

0.1

1

0.3

0

0 0.2 -1

0.1

-0.1 50

20

0

40

60

80

52

54

56

100

-0.2 10

20

30

40

50

60

roo

0

f

-0.1

70

80

90

100

Figure 7: Online PID control design for tracking a sinusoidal wave in the presence of disturbance.

rep

The next simulations are performed when the parameters   B = 8, a = 0.3, φ = 20, Da = 0.072 B = 10, a = 0.5, φ = 21, Da = 0.076  B = 11, a = 0.6, φ = 19, Da = 0.076

of Eq. (1) vary with time as 0 ≤ t < 30 30 ≤ t < 60 t ≥ 60

Jo u

rna

lP

where the results of tracking errors are illustrated in Fig. (8). Since the system is time-varying in this case, the RLS-FF factor with λ = 0.95 has been used to determine the varying model’s parameters. Controlling a nonlinear system with varying parameters is a challenging task, however, in such a condition, the proposed controller design has successfully tracked the sinusoidal wave with superior performance than MPSO and ABC Fig. (8).

13

Journal Pre-proof

10

4

-3

2 0 -2 5

10 10

3

15

20

25

30

-3

2 1 0 -1 -2 30

35 10

2

40

45

50

55

60

-3

1

65

70

75

80

85

90

95

100

roo

-1 60

f

0

Figure 8: Online PID control design for tracking a sinusoidal wave in the presence of varying system

rep

dynamics.

lP

The considered CSTR has three equilibrium points [10]: two stable as p1 : (x1 , x2 ) = (0.144, 0.886), p2 : (x1 , x2 ) = (0.7665, 4.705), and one unstable p3 : (x1 , x2 ) = (0.445, 2.74). Transferring the states of a system from a stable point to an unstable point is demanding and needs a highperformance controller to be designed. In this paper, the proposed controller is employed to transfer the states from p1 and p2 to p3 . Consequently, the reference signal is set to r(t) = 0.445, and the initial values of states are selected as p1 and p2 , where the phase-plane trajectories of the obtained results are provided in Fig. (9). As can be observed, the stable equilibrium points have been successfully transferred to the unstable one.

rna

5 4.5 4 3.5

Jo u

3 2.5 2

1.5 1

0.5 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Figure 9: Phase-plane trajectories of controlling equilibrium points using the proposed controller.

14

Journal Pre-proof

4 Conclusion The data-driven identification and control design of a highly-nonlinear CSTR were presented in this paper. The data set used in the black-box identification was gathered in the presence of white and colored noises to consider realistic conditions. Then, linear model structures were selected, and based on the offline and online estimation methods, the best representations of the data by models were obtained. By employing the results of the online identification algorithms, a selftuning PID controller was designed. The main achievements of the paper can be summarized as follows: (i) It was validated that linear models such as ARX and ARMAX can efficiently capture the CSTR dynamics around a specific operating point. (ii) To represent the nonlinear CSTR entirely, the recursive algorithms, such as the RLS, must be employed to cover all the operating points by adapting the linear models to the new conditions. (iii) If there is an uncertainty in the initial parameters, the parameter P −1 (0) should be large enough to decrease the prediction error more efficiently.

roo

f

(iv) The RLS with exponential forgetting factor efficiently dealt with time-varying situations in which, the smaller the forgetting factor λ was, the faster new condition was captured. (v) Since the nonlinear CSTR can be modeled as a linear system with time-varying parameters, the developed controllers for linear systems such as the PID controller can be applied to the CSTR with time-varying coefficients.

rep

(vi) By using the RLS with exponential forgetting factor as the identifier in the controller design, the ability to control the nonlinear CSTR in the presence of unwanted signals (i.e., noise and disturbance) and variation in system dynamics were provided. The comparison of the proposed method with two reported works in the literature validated the applicability and the effectiveness of the presented approach.

Jo u

rna

lP

Appendix

15

Model ARX(2,1) ARX(2,2) ARX(3,2) ARX(3,3) ARMAX(2,2,1) ARMAX(3,3,3) BJ(2,3,2,3) BJ(3,3,3,3) OE(2,2) OE(3,3)

lP

A(q) 1 − 1.961q −1 + 0.9617q −2 1 − 1.918q −1 + 0.9211q −2 1 − 2.348q −1 + 1.796q −2 − 0.4454q −3 1 − 1.869q −1 − 0.8921q −2 − 0.8921q −3 1 − 1.906q −1 + 0.9092q −2 1 − 1.832q −1 + 0.8261q −2 − 0.0115q −3 -

rep

B(q) 9.806 × 10−6 q −2 −2.719 × 10−6 q −1 + 9.214 × 10−5 q −2 −2.907 × 10−6 q −1 + 9.19 × 10−5 q −2 −2.819 × 10−6 q −1 + 8.636 × 10−5 q −2 + 7.601 × 10−5 q −3 −3.787 × 10−6 q −1 + 1.1 × 10−4 q −2 −1.263 × 10−6 q −1 + 8.771 × 10−5 q −2 + 8.11 × 10−5 q −3 4.856 × 10−6 q −1 + 1.1 × 10−4 q −2 8.857 × 10−5 q −1 − 4.17 × 10−6 q −2 − 8.25 × 10−5 q −3 5.4743 × 10−5 q −1 + 1.1 × 10−4 q −2 8.407 × 10−5 q −1 + 4.748 × 10−6 q −2 + 4.97 × 10−5 q −3

f

roo

C(q) 1 − 0.2874q −1 −1 1 + 0.2191q + 0.2956q −2 + 0.2954q −3 1 − 0.7498q −1 − 0.02954q −2 + 0.2285q −3 1 − 0.2117q −1 − 0.167q −2 − 0.145q −3 -

D(q) 1 − 1.828q −1 + 0.8347q −2 1 − 1.79q −1 + 0.7q −2 + 0.903q −3 -

Table 2: The detailed results of identifying CSTR using offline estimation techniques.

rna

Jo u F (q) 1 − 2.271q −1 − 1.656q −2 − 0.3817q −3 1 − 2.822q −1 + 2.654q −2 − 0.8313q −3 1 − 1.848q −1 + 0.8535q −2 1 − 2.068q −1 + 1.266q −2 − 0.1934q −3

Fit. (train) 68.12 86.1911 89.0281 94.4920 88.63 94.5694 96 97.52 95.0849 95.3105

Fit. (test) 5.42 60.37 82.5394 86.1367 68.38 86.657 90.9525 96.60 85.73 85.64

Journal Pre-proof

16

Journal Pre-proof

References [1] D. Zhao, Q. Zhu, and J. Dubbeldam, “Terminal sliding mode control for continuous stirred tank reactor,” Chemical engineering research and design, vol. 94, pp. 266–274, 2015. [2] B. W. Bequette, Process Control: Modeling, Design and Simulation. Prentice Hall, 1 ed., 2003. [3] A. Favache and D. Dochain, “Power-shaping control of reaction systems: The cstr case,” Automatica, vol. 46, no. 11, pp. 1877–1883, 2010. [4] C.-T. Chen and C.-S. Dai, “Robust controller design for a class of nonlinear uncertain chemical processes,” Journal of Process Control, vol. 11, no. 5, pp. 469–482, 2001. [5] C.-T. Chen and S.-T. Peng, “Intelligent process control using neural fuzzy techniques,” Journal of process control, vol. 9, no. 6, pp. 493–503, 1999. [6] L. Ljung, “System identification,” Wiley Encyclopedia of Electrical and Electronics Engineering, 2001.

roo

f

[7] V. Pascu, H. Garnier, L. Ljung, and A. Janot, “Benchmark problems for continuous-time model identification: Design aspects, results and perspectives,” Automatica, vol. 107, pp. 511– 517, 2019. [8] R. S. Risuleo, F. Lindsten, and H. Hjalmarsson, “Bayesian nonparametric identification of wiener systems,” Automatica, vol. 108, p. 108480, 2019.

rep

[9] A. Sinha and R. K. Mishra, “Temperature regulation in a continuous stirred tank reactor using event triggered sliding mode control,” IFAC-PapersOnLine, vol. 51, no. 1, pp. 401 – 406, 2018. 5th IFAC Conference on Advances in Control and Optimization of Dynamical Systems ACODS 2018. [10] W.-D. Chang, “Nonlinear CSTR control system design using an artificial bee colony algorithm,” Simulation Modelling Practice and Theory, vol. 31, pp. 1–9, 2013.

lP

[11] A. Simorgh, A. Razminia, and J. T. Machado, “Optimal control of nonlinear fed-batch process using direct transcription method,” Computers & Chemical Engineering, vol. 130, p. 106561, 2019.

rna

[12] A. Sinha and R. K. Mishra, “Control of a nonlinear continuous stirred tank reactor via event triggered sliding modes,” Chemical Engineering Science, vol. 187, pp. 52–59, 2018. [13] B. Daaou, A. Mansouri, M. Bouhamida, and M. Chenafa, “Development of linearizing feedback control with a variable structure observer for continuous stirred tank reactors,” Chinese Journal of Chemical Engineering, vol. 20, no. 3, pp. 567 – 571, 2012.

Jo u

[14] S. Ge, C. Hang, and T. Zhang, “Nonlinear adaptive control using neural networks and its application to cstr systems,” Journal of process control, vol. 9, no. 4, pp. 313–323, 1999. [15] L. Dong-Juan, “Adaptive neural network control for continuous stirred tank reactor process,” IFAC Proceedings Volumes, vol. 46, no. 20, pp. 171 – 175, 2013. 3rd IFAC Conference on Intelligent Control and Automation Science ICONS 2013. [16] D.-J. Li, “Adaptive neural network control for a two continuously stirred tank reactor with output constraints,” Neurocomputing, vol. 167, pp. 451 – 458, 2015. [17] E. Salahshour, M. Malekzadeh, F. Gordillo, and J. Ghasemi, “Quantum neural network-based intelligent controller design for CSTR using modified particle swarm optimization algorithm,” Transactions of the Institute of Measurement and Control, vol. 41, no. 2, pp. 392–404, 2019.

17

Journal Pre-proof

[18] V. Ghaffari, S. V. Naghavi, and A. Safavi, “Robust model predictive control of a class of uncertain nonlinear systems with application to typical CSTR problems,” Journal of Process Control, vol. 23, no. 4, pp. 493–499, 2013. [19] O. Paladino and M. Ratto, “Robust stability and sensitivity of real controlled CSTRs,” Chemical Engineering Science, vol. 55, no. 2, pp. 321–330, 2000. [20] J. Gerhard, M. Mönnigmann, and W. Marquardt, “Robust stable nonlinear control and design of a cstr in a large operating range,” IFAC Proceedings Volumes, vol. 37, no. 9, pp. 209–214, 2004. [21] J. Sjöberg, Q. Zhang, L. Ljung, A. Benveniste, B. Delyon, P.-Y. Glorennec, H. Hjalmarsson, and A. Juditsky, “Nonlinear black-box modeling in system identification: a unified overview,” Automatica, vol. 31, no. 12, pp. 1691–1724, 1995. [22] K. Chen and R. Fung, “A mechatronic motor-table system identification based on an energetics fitness function,” IEEE/ASME Transactions on Mechatronics, vol. 22, pp. 2288–2295, Oct 2017.

roo

f

[23] J. Pan, S. W. Or, Y. Zou, and N. C. Cheung, “Sliding-mode position control of medium-stroke voice coil motor based on system identification observer,” IET Electric Power Applications, vol. 9, no. 9, pp. 620–627, 2015. [24] Z. Yang, D. Zhang, X. Sun, and X. Ye, “Adaptive exponential sliding mode control for a bearingless induction motor based on a disturbance observer,” IEEE Access, vol. 6, pp. 35425– 35434, 2018.

rep

[25] S. Gupta, R. Gupta, and S. Padhee, “Parametric system identification and robust controller design for liquid–liquid heat exchanger system,” IET Control Theory Applications, vol. 12, no. 10, pp. 1474–1482, 2018.

lP

[26] M. Espinoza, J. Rojas, R. Vilanova, and O. Arrieta, “Identification and control of chemical processes using the anisochronic modeling paradigm,” IFAC-PapersOnLine, vol. 48, no. 8, pp. 361 – 366, 2015. 9th IFAC Symposium on Advanced Control of Chemical Processes ADCHEM 2015.

rna

[27] A. Simorgh, A. Razminia, and V. I. Shiryaev, “Data-driven identification of a continuous type bioreactor,” Energy Sources, Part A: Recovery, Utilization, and Environmental Effects, pp. 1–29, 2019. [28] O. Khan, C. M. R. Madhuranthakam, P. Douglas, H. Lau, J. Sun, and P. Farrell, “Optimized PID controller for an industrial biological fermentation process,” Journal of Process Control, vol. 71, pp. 75–89, 2018.

Jo u

[29] Y. Wang, R. Rajamani, and A. Zemouche, “A quadratic matrix inequality based PID controller design for lpv systems,” Systems & Control Letters, vol. 126, pp. 67–76, 2019. [30] S. M. H. Mousakazemi and N. Ayoobian, “Robust tuned PID controller with PSO based on two-point kinetic model and adaptive disturbance rejection for a PWR-type reactor,” Progress in Nuclear Energy, vol. 111, pp. 183 – 194, 2019. [31] L. Ljung and T. Söderström, Theory and practice of recursive identification. MIT press, 1983. [32] V. Bobál, J. Böhm, J. Fessl, and J. Machácek, Digital self-tuning controllers: algorithms, implementation and applications. Springer Science & Business Media, 2006. [33] D. Valério and J. S. Da Costa, “Tuning of fractional PID controllers with ziegler–nichols-type rules,” Signal processing, vol. 86, no. 10, pp. 2771–2784, 2006.

18