ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
Research Article
A fast closed-loop process dynamics characterization Miroslav R. Mataušek n, Tomislav B. Šekara Faculty of Electrical Engineering, University of Belgrade, Belgrade 11120, Serbia
art ic l e i nf o
a b s t r a c t
Article history: Received 17 January 2013 Received in revised form 9 October 2013 Accepted 3 December 2013
Stable, integrating and unstable processes, including dead-time, are analyzed in the loop with a known PI/PID controller. The ultimate gain and frequency of an unknown process Gp(s), and the angle of tangent to the Nyquist curve Gp(iω) at the ultimate frequency, are determined from the estimated Laplace transform of the set-point step response of amplitude r0. Gain Gp(0) is determined from the measurements of the control variable and known r0. These estimates define a control relevant model Gm(s), making possible the use of the previously determined and memorized look-up tables to obtain PID controller guaranteeing desired maximum sensitivity and desired sensitivity to measurement noise. Simulation and experimental results, from a laboratory thermal plant, are used to demonstrate the effectiveness and merits of the proposed method. & 2013 ISA. Published by Elsevier Ltd. All rights reserved.
Keywords: PID control Closed-loop identification Measurement noise Robustness Tuning
1. Introduction The importance of control relevant process dynamics characterization is discussed in detail in an overview paper [1], and reexamined recently in [2]. One of the basic conclusions from [1,2] follows the idea from [3], where it is proposed to estimate a highorder ARX model first and then to perform a model reduction in the frequency domain, to obtain a reduced-order model used for controller tuning. Three problems are related to the adequate control relevant process dynamics characterization: the model structure, the control relevant region of frequencies and the dilemma open-loop versus closed-loop process identification. Closed-loop system performance/robustness tradeoff strongly depends on a priory knowledge defined by the structure of the model used for the PID controller tuning. The two parameter model, obtained from the Ziegler–Nichols time domain tuning [4], can be used to tune PID controller for lag-dominated processes Gp(s) [5]. This model is known at the present time as integrator plus dead-time (IPDT) model. Its relationship with the Ziegler– Nichols frequency domain tuning is discussed in detail in [5]. The second one, used for the controller tuning, is the three parameter model. It is represented by the first-order plus dead-time (FOPDT) model, introduced first by Cohen–Coon [6], and by the integrating first-order plus dead-time (IFOPDT) model. FOPDT model can be used to approximate balanced and dead-time dominated processes, while IFOPDT model defines a better approximation of
n
Corresponding author. Tel.: þ 381 11 337 01 64; fax: þ381 11 324 86 81. E-mail addresses: [email protected] (M.R. Mataušek), [email protected] (T.B. Šekara).
the lag-dominated processes, than the IPDT model [7]. The four parameter second-order plus dead-time (SOPDT) model, and other reduced-order models, used for the PI/PID controller tuning, are summarized in [8]. Mathematical models are always an approximation of reality. Appearance of a dead-time in a model might be a consequence of low-order modeling [9,10], for example based on the FOPDT, IFOPDT or SOPDT models, or might be a consequence of the real time-delay caused by some physical phenomena [11]. The control relevant region of frequencies is the region around the ultimate frequency ωu of a process Gp(s). This is confirmed by the PI/PID controller optimization [12–15], based on the frequency response of the process Gp(iω), under constraints on the robustness. The frequency ω0, where the sensitivity function has its maximum, occurs in the region around the ultimate frequency ωu, as demonstrated in [14]. The importance of the process dynamics characterization based on the ultimate frequency estimation is firstly recognized by Ziegler and Nichols [4] and further developed by Åström and Hägglund [16,17]. The extension of the Ziegler– Nichols process dynamics characterization proposed in [14], by introducing an additional parameter in the frequency domain, improved considerably the possibility of better process modeling in the wider region of frequencies around the ultimate frequency ωu of a large class of stable processes, processes with oscillatory dynamics, integrating processes and unstable processes Gp(s), including dead-time. Dilemma open-loop versus closed-loop process identification in academia is treated mainly as a problem of ensuring the best statistical accuracy [1,3]. However, in industry, breaking of control loops in operation is mainly ignored by plant operators. When some initially tuned controller is in operation, a procedure for fine tuning can be easily accepted if it can be activated/deactivated
0019-0578/$ - see front matter & 2013 ISA. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.isatra.2013.12.006
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
2
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
without breaking control loops in operation [18]. A special attention should be devoted when estimating parameters of a continuous model including dead-time. For example, this can be done by applying the extended polynomial closed-loop identification [9], or to obtain model parameters, including dead-time, from the estimated frequency response of the process, obtained from the closed-loop step response tests [10,11]. However, in these methods, the dead-time free part of the model transfer function must be specified, based on some a priory knowledge, and different model structure is used for stable, integrating and unstable processes. In this paper, a new insight into the problem of control relevant process dynamics characterization and an effective solution of this problem is presented. It is based on a continuous model Gm(s) with the unified structure for a large class of stable processes, processes with oscillatory dynamics, integrating processes and unstable processes Gp(s), including dead-time. Proposed in [14] as an effective extension of the Ziegler–Nichols process dynamics characterization in the frequency domain, this model is defined by the quadruplet {ku, ωu, φ, A}. Parameter ku is the ultimate gain and φ ¼ argð∂Gp ðiωÞ=∂ωÞjω ¼ ωu is the angle of the tangent to the Nyquist curve Gp(iω) at the ultimate frequency ωu of a process Gp(s). In this model an equivalent dead-time τ is defined by τ¼ φ/ωu. Two procedures [15,18] are proposed for estimating the quadruplet {kuest, ωuest, φest, Aest}. Parameter Aest ¼A0 is defined in the frequency domain by A0 ¼ 2j∂Gp ðiωÞ=∂ωjω¼1 ωu =ku [15]. The new PLL (phase-locked-loop) estimator [15], further improved in [5], requires only that the controller in operation is a linear controller, while the relay SheMa estimator [18] does not require this preliminary information. Estimation methods [5,15,18] are performed without breaking the loop containing a controller in operation. The model Gm(s), defined by the quadruplet {kuest, ωuest, φest, Aest}, approximates the Nyquist curve of a large class of stable processes, processes with oscillatory dynamics, integrating and unstable processes Gp(s), including dead-time, in a large region around ωu, and can be effectively used in PID controller constrained optimization [5]. However, from the industry viewpoint, disadvantage of estimation methods [5,15,18] can be a longer period of time required for determining the quadruplet {kuest, ωuest, φest, Aest}. Namely, determination of φest and Aest by applying PLL and SheMa estimators requires two additional experiments, as explained in Appendix A. Besides, and this is of essential importance, the basic definition of parameter A¼ ωukuGp(0)/(1 þkuGp(0)) [14] offers the possibility to classify a large class of stable processes, processes with oscillatory dynamics, integrating and unstable processes, including deadtime, in the ρ φ plane [19], where ρ ¼κ/(1þ κ), κ ¼kuGp(0). Stable and unstable processes are classifiedpas processes inside and ffiffiffiffiffiffiffiffiffiffi outside the region 0 oρo1, 0 o φ oπ= ρ þ1, while the integratpffiffiffi ing processes are classified as processes with ρ¼ 1, 0 o φ o π= 2, since Gp(0)¼ 71. The possibility to classify stable processes, processes with oscillatory dynamics, integrating, and unstable processes, including deadtime, in a two parameter plane, is important from the PID controller tuning viewpoint. The desired performance/robustness tradeoff can be obtained by applying the previously memorized process indepen-
Fig. 1. Process Gp(s), with a two-degree-of-freedom controller. The set-point r(t), controlled variable y(t), control variable u(t), load and output disturbances, d(t) and n(t), represent variations around their values in the nominal regime.
dent look-up tables in this ρ φ plane [19]. However, to be effective, such a possibility must be supported by a fast estimation of the quadruplet {ku, ωu, φ, Gp(0)}. Then, from the estimated {ku, ωu, φ, Gp(0)} and the previously memorized look-up tables one obtains directly the gains and the noise filter time constant of a real parallel PID controller, as presented in Appendix B, for stable processes and illustrated by experimental results. In the present paper a new simple and effective closed-loop procedure is proposed for estimating quadruplet {ku, ωu, φ, Gp(0)} in a short time interval, without braking the loop with the controller in operation. It is assumed that the controller in operation is a known linear controller. The estimation procedure is defined in Section 2. In Section 3, a test batch, consisting of stable processes, process with oscillatory dynamics, integrating and unstable processes, including dead-time, in the loop with PI and PID controllers, is used to demonstrate the properties of the proposed estimation method. Finally, in Section 4, the experimental verification of the proposed method is presented. The quadruplet {kuest, ωuest, φest, Gpest(0)} is determined by applying a PI controller to a laboratory thermal process [20], with noisy measurements, and used for the PID controller tuning by applying the previously memorized look-up tables in the ρ φ plane, presented in Appendix B.
2. Determination of the quadruplet {ku, ωu, φ, Gp(0)} from the closed-loop system set-point step response Model Gm(s) of an unknown stable process, process with oscillatory dynamics, integrating process and unstable process, 7 6 5 4 3 2 1 0
0
100
200 300 Time [sec]
400
500
Fig. 2. The set-point (dashed), the measured closed-loop response y(mTs) of the real plant (solid-blue), the data yj (circles-black) and yid(t), defined by (solid-red) lines for 0r tr 400 s, obtained by linear interpolation of yj. For t4400 s the response is approximated with its value defined by the known set-point r0 ¼ 5 and presented also by solid-red line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 1 Parameters of the controllers used to obtain closed-loop set-point step responses for processes Gpl(s), l ¼1,2,…,7. Process/controller
k
ki
kd
Tf
b
Gp1(s)/PI Gp2(s)/PI Gp3(s)/PI Gp4(s)/PID1 Gp5(s)/PI Gp6(s)/PI Gp7(s)/PI
2.1631 0.7903 0.1204 0.5620 0.1010 1.3273 0.3553
4.7980 0.0654 0.0946 0.0830 0.00255 0.01354 0.0030
0 0 0 1.0620 0 0 0
0 0 0 0.1770 0 0 0
1 1 1 0 1 0.15 0.25
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
From Fig. 1, for d ¼0 and n ¼ 0, one obtains that process transfer function Gp(s) is defined by
including dead-time, proposed in [14], is defined by relations Gm ðsÞ ¼
WðsÞ 1 Aωu expð τsÞ ; WðsÞ ¼ ; τ ¼ φ=ωu ; 1 WðsÞ ku s2 þ ω2u
ð1Þ Gp ðsÞ ¼
A ¼ ρωu ; ρ ¼ ku Gp ð0Þ=ð1 þku Gp ð0ÞÞ;
3
ð2Þ
HðsÞ YðsÞ ; ; HðsÞ ¼ C ff ðsÞ HðsÞCðsÞ RðsÞ
ð3Þ
where H(s) is the closed-loop system transfer function from r to y in Fig. 1. In the present paper, the step set-point R(s) ¼r0/s, of known amplitude r0, is used to collect measurements of the controlled variable applied to obtain an estimate Yest(s). Then, for
and by the quadruplet {ku, ωu, φ, Gp(0)}. Unknown process Gp(s), in the loop with known linear controllers C(s) and Cff(s), is presented in Fig. 1.
Table 2 True values ku, ωu, φ, Gp(0) compared with the estimated kuest, ωuest, φest, Gpest(0) of processes Gpl(s), l ¼1,2,…,7, obtained for the three pairs T1J1 ¼ T2J2 ¼ T3J3, T1 o T2 o T3. Time period T is in seconds. Process
ku
ωu
φ
Gp1(s)
11.5919
9.8696
0.7854
1
Gp2(s)
4.6626
0.2144
0.8271
1
Gp3(s)
2.5443
0.5884
0.5648
1
Gp4(s)
0.5640
0.4080
0.7200
1
Gp5(s)
0.2371
0.2291
0.9716
1
Gp6(s)
3.1865
0.4287
0.3903
1
Gp7(s)
0.6341
0.5828
0.7603
4
Gp(0)
kuest
ωuest
φest
Gpest(0)
T
J
11.660 11.791 12.335 4.7021 4.7675 4.8814 2.5448 2.6287 2.9567 0.5641 0.5673 0.5731 0.2380 0.2405 0.2514 3.1901 3.2015 3.2307 0.6358 0.6415 0.6460
9.8817 9.8796 9.8628 0.2143 0.2142 0.2140 0.5870 0.5881 0.5850 0.4085 0.4100 0.4127 0.2291 0.2289 0.2285 0.4280 0.4260 0.4226 0.5836 0.5835 0.5835
0.8055 0.8032 0.7659 0.7828 0.7695 0.7480 0.6074 0.6244 0.5900 0.7167 0.7125 0.7047 0.9628 0.9528 0.9083 0.3920 0.3937 0.5632 0.8289 0.8193 0.8105
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1 1 1 1 1 1 1.000 1.000 1.000 4.005 4.005 4.005
0.02 0.04 0.08 1 2 3 0.5 1 2 0.5 1 1.5 1 2 4 0.8 1 2 0.4 0.8 1
200 100 50 150 75 50 200 100 50 120 60 40 200 100 50 200 160 80 200 100 80
0.15
0.03 0.02
0.1
0.01 0.05
0 -0.01
0
-0.02
-0.05
-0.03
-0.1
Gm1(s) Gp1 (s)
-0.04 -0.05 -0.06 -0.12
-0.11
-0.1
-0.09
-0.08
-0.07
-0.06
Gm3 (s) Gp3 (s)
-0.15
-0.05
4
-0.2 -0.6
-0.5
-0.4
-0.3
-0.2
0.015 0.01
2
0.005 0
0 -0.005
-2
-0.01
-4
-8
-0.015
Gm5 (s) Gp5 (s)
-6
-7
-6
-5
-4
-3
-2
-1
0
Gp6 (s) Gm6 (s)
-0.02
1
-0.025 -0.38
-0.36
-0.34
-0.32
-0.3
-0.28
-0.26
Fig. 3. Nyquist curves, around the critical point ( 1/ku,i0): (a) stable process Gp1(s); (b) process with oscillatory dynamics Gp3(s); (c) integrating process Gp5(s); and (d) unstable process Gp6(s). In all cases: exact (solid-blue), estimated (dashed-red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
H(s) ¼sYest(s)/r0 and for the known linear controllers C(s) and Cff(s), from (3) one obtains an estimate Gpest(s), used to determine parameters kuest, ωuest and φest. In industry applications PI controller predominates [21]. Thus, in the present paper it is assumed that the two-degree-of-freedom PI controller, or parallel PID controller is used to determine the quadruplet {ku, ωu, φ, Gp(0)}. In (3) and Fig. 1, these controllers are given by C(s)¼ (kds2 þks þ ki)/(s(Tfsþ 1)) and Cff(s)¼(bksþki)/s. They are implemented as defined by relations k UðsÞ ¼ ðbkRðsÞ kY f ðsÞÞ þ i ðRðsÞ Y f ðsÞÞ kd sY f ðsÞ; 0 rb; Y f ðsÞ s YðsÞ ; ¼ Tf sþ1
ð4Þ
used also in [13,15]. For the PI controller kd ¼0, and it is adopted that Tf ¼0 as in [12,15]. The proposed estimation of quadruplet {kuest, ωuest, φest, Gp(0)} is performed in three steps. In the first step, a small sampling period Ts is used to collect the closed-loop system set-point step response y(mTs) and the control signal u(mTs) m ¼ 0,1,2,…,M. To reduce the effect of measurement noise, the averaging procedure is applied to obtain the signal yj defined by yj ¼
1 jn 1 yðmT s Þ; j ¼ 1; 2; …; J; y0 ¼ 0; yJ þ 1 ¼ r 0 : ∑ n m ¼ ðj 1Þn
ð5Þ
In (5) parameter n, significantly greater than 1, is defined by period T ¼nTs. The time period T is defined as Tp ¼pT, where Tp is the time-to-go from 10% of the set-point step response first maximum to the first maximum. Alternatively, if the well damped closed-loop system set-point step response is obtained by the implemented PI/PID controller, it is recommended to define Tp as the rise time. For integer p it is recommended to use p¼ 5–20. Then, signal yid(t), used to determine the Laplace transform Yest(s), is obtained by applying the linear interpolation of signal yj, j¼1,2, …,J. For j ZJþ 1 it is adopted yj ¼r0. The reason is obvious. The relation lim yid ðtÞ ¼ r 0 holds true for the closed-loop system t-1 set-point step response of a process in the loop with the PI/PID controller. This first step is illustrated in Fig. 2. The measured closed-loop set-point step response y(mTs), obtained experimentally for Ts ¼0.5 s on the real thermal plant from Section 4, is presented in Fig. 2 with data yj and signal yid(t), obtained by linear interpolation of yj. For the time time-to-go to the first maximum equal to Tp E140 s and for p¼ 7 one obtains T¼ 20 s. A small number of points J¼ 20 is used to calculate Yest(s). The response is truncated at t¼ 400 s to demonstrate that the extremely short time interval JTþ T/2 is necessary for obtaining adequate estimates
kuest, ωuest, φest and Gpest(0), guaranteeing an excellent tuning of the PID controller, as demonstrated in Section 4. In the second step the Laplace transform Yest(s) is determined, by applying data yid(t) prepared in the first step of the proposed estimation procedures. Taking into account that for j ZJ þ1 it is adopted yj ¼r0, the Laplace transform Yest(s) is defined by the relation Z Y est ðsÞ ¼
JT þ T=2 0
yid ðtÞe st dt þ
r 0 sðJT þ T=2Þ e ; s
ð6Þ
where the time JTþ T/2 corresponds to the settling time estimated from the memorized closed-loop system set-point step response. Data yj, j¼ 1,2,…,J, used to define signal yid(t) by the linear interpolation, appears at points jT T/2. Thus, from (6) one obtains 2y Y est ðsÞ ¼ 1 T þ
Z
T=2
te
st
J
dt þ ∑
j¼1
0
Z
jT þ T=2 jT T=2
! ðaj þ bj tÞe
st
dt
r 0 sðJT þ T=2Þ e ; s
ð7Þ
where aj ¼yj (yj þ 1 yj)(j 1/2) and bj ¼ (yj þ 1 yj)/T. By solving integrals in (7), the analytical solution of the Laplace transform Yest(s) is obtained, defined by the measured data and following relations: Y est ðsÞ ¼ 2y1
1 ð1 þsT=2Þe sT=2 r0 þ Y 1est ðsÞ Y 2est ðsÞ þ e sðJT þ T=2Þ ; s s2 T ð8Þ
ðsTyj þyj þ 1 yj Þe sðjT T=2Þ ; s2 T j¼1 J
ð9Þ
Y 1est ðsÞ ¼ ∑
Fig. 5. The laboratory thermal plant. Positions of temperature sensors are denoted by 2 and 3, while position of the heater is denoted by 1.
0.1 0 Load disturbance response
Set-point response
1 0.8 0.6 0.4 0.2 0
0
10
20
30 Time [sec]
40
50
-0.1 -0.2 -0.3 -0.4 -0.5 -0.6
0
10
20
30
40
50
Time [sec]
Fig. 4. Process Gp4(s) and its model Gm4(s), defined by kuest ¼0.5731, ωuest ¼0.4127, φest ¼0.7047 and Gpest(0)¼ 1, in the loop with PIDtun controller. Responses of the controlled variable for b¼ 0: (a) set-point step response, R(s) ¼ 1/s; (b) load disturbance step response, D(s) ¼ 0.2/s.
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
ðsTyj þ 1 þ yj þ 1 yj Þe sðjT þ T=2Þ : s2 T j¼1 J
ð10Þ
Y 2est ðsÞ ¼ ∑
It can be easily confirmed that the condition limHðsÞ ¼ 1 is s-0 satisfied, for the closed-loop system transfer function defined by H(s) ¼sYest(s)/r0, with Yest(s) given by (8)–(10). A time scaled set-point response is used to perform calculations, with the time scale factor defined by JT. If higher values of Ts are used for memorizing the set-point step response, the following improvement is proposed. Instead from (5), the signal yj, j¼1,2,…, R jT J, can be obtained from the integral yj ¼ ð ðj 1ÞT ym ðtÞdtÞ=T, where ym(t) is defined by the linear interpolation of y(mTs). In the third step, from H(s)¼sYest(s)/r0 and (3), the process transfer function Gpest(s) is determined and used to estimate the ultimate frequency ωuest from argðGpest ðiωuest ÞÞ ¼ π. Then, the ultimate gain kuest is obtained from kuest ¼ jGpest ðiωuest Þj 1 and the angle φest is obtained from φest ¼ argð∂Gpest ðiωÞ=∂ωÞjω ¼ ωuest , as defined in [14]. The gain Gp(0) is estimated from 8 > Z JT þ T=2 < r 0 =u0 ; for u0 a 0; 1 for u0 -0; kuest 4 0; um ðtÞdt; Gpest ð0Þ ¼ 1; u0 ¼ > T JT T=2 : 1; for u -0; k o 0: uest
0
ð11Þ by using the known amplitude r0 of the set-point response and linear interpolation of the memorized control signal u(mTs), m ¼0,1,2,…,M, used to calculate um(t) in (11).
3. Simulation analysis of the proposed estimation method Effectiveness of the proposed method is demonstrated first by the following test batch: 1 ð2s þ 1Þe 4s pffiffiffiffiffi; Gp2 ðsÞ ¼ ; Gp1 ðsÞ ¼ ð10s þ 1Þð7s þ 1Þð3s þ1Þ cosh 2s Gp3 ðsÞ ¼
es e 0:5s ; Gp4 ðsÞ ¼ 9s2 þ 2:4s þ 1 sð1:2s þ 1Þ3
Gp5 ðsÞ ¼
e 5s ; sðs þ 1Þð0:5s þ 1Þð0:25s þ 1Þð0:125s þ 1Þ
Gp6 ðsÞ ¼
e 0:5s 4e 2s ; Gp7 ðsÞ ¼ ð5s 1Þð2s þ 1Þð0:5s þ 1Þ 4s 1
5
The impact of reducing the number of points J on estimation is analyzed for the three pairs T1J1 ¼T2J2 ¼ T3J3, such that T1 oT2 oT3. As demonstrated in Table 2, the impact of reducing considerably number of points J is small in all examples, even for unstable processes. In Fig. 3, Nyquist curves of processes Gpl(s), l ¼1,3,5,6, are compared for T1 and J1 with the corresponding Nyquist curves of the estimated Gpestl(s). Excellent estimation of the Nyquist curves, presented in Fig. 3, is obtained within a large region around the critical point ( 1/ku,i0). This is important characteristic of the proposed identification method, since only this result makes possible to determine parameter φ. An important quality of the proposed identification method is demonstrated in Fig. 4. The closed-loop set-point and load disturbance step responses of Gp4(s) and Gm4(s), defined by kuest, ωuest, φest and Gpest(0) in Table 2 (Gp4, T¼1.5, J¼ 40), in the loop with the PIDtun controller, are almost exactly equal. Parameters of the PIDtun controller in (4) are obtained by applying kuest, ωuest, φest and Gpest(0) from Table 2 (Gp4, T¼1.5, J¼ 40), and tuning formulae defined by equations (30)–(37) in [15], to obtain critically damped responses for: k¼0.3787, ki ¼ 0.0324, kd ¼ 0.8453, Tf ¼0.2423, b¼ 0. Results presented in Fig. 4 confirm that the performance of the considered loop can be analyzed by applying the proposed identification method and simulation of the tuned controller in the loop with the model Gm(s), defined by (1) and (2), instead of the real process.
Table 3 Estimated values kuest, ωuest, φest, Gpest(0), obtained for two values of J from the real plant closed-loop set-point step response of amplitude r0 ¼ 5 1C, presented in Fig. 6b. Plant
kuest
ωuest
φest
Gpest(0)
ρest
Ts [s]
T [s]
J
(a) (b)
28.6582 27.0675
0.04458 0.04461
0.6377 0.6814
0.4104 0.4082
0.9216 0.9170
0.5 0.5
20 10
20 40
Table 4 Parameters of the PID controllers obtained, for the laboratory plant quadruplet {kuest, ωuest, φest, Gpest(0)} of Table 3, by applying tuning defined by look-up Tables B1–B3 of Appendix B.
of processes in the loop with the PI/PID controllers (4), defined by Table 1. The estimates kuest, ωuest, φest and Gpest(0) are compared in Table 2 with the true values of ku, ωu, φ and Gp(0).
PIDplant
K
ki
kd
Tf
J
Ms
Mn
(a) (b)
14.3795 13.5302
0.1774 0.1668
378.7222 340.0592
6.6076 6.2817
20 40
2.00 2.00
57.32 54.13
65
52
60
51
55
50
Controlled variable
Control variable
Set-point
50
49
45 48 40 47
35
46
30
45
25 20
44 1000
1200
1400
1600
Time [sec]
1800
1000
1200
1400
1600
1800
Time [sec]
Fig. 6. Experimental results: (a) and (b) set-point step responses of the plant with PIplant controller. Response in (b) is used to calculate Yest(s), as defined by Fig. 2. Control variable in [%], controlled variable in [1C].
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
of temperature at x¼ 0 is used to keep the temperature T(0,t) r70 1C [20]. The closed-loop set-point step responses in Fig. 6a–b are obtained on the laboratory thermal plant with the PIplant controller (4), defined by: k ¼6.3810, ki ¼0.055 and b¼1. The step response shown in Fig. 6b is used (as defined by Fig. 2) to obtain Yest(s) from (8)–(10) and parameters kuest, ωuest, φest in Table 3. The static gain Gpest(0) in Table 3 is obtained from (11) and the response shown Fig. 6a. By visual inspection one obtains Gp(0)E 0.4 1C/(%) from Fig. 6, which is in accordance with the estimate Gpest(0), given in Table 3.
4. Experimental verification by a laboratory thermal plant The laboratory thermal plant, presented in Fig. 5, consists of a thin plate made of aluminum, long L ¼0.1 m and wide h ¼0.03 m [20]. The temperature T(x,t) is distributed along the plate, from x ¼0 to x ¼L, and measured by precision sensors LM35 (TO92), at positions x ¼0 and x ¼L. The plate is heated by terminal adjustable regulator LM317 (TO 220) at position x ¼0. The manipulated variable is dissipated power of the heater at x¼ 0. The input to heater is control variable u(t) [%], defined by the output of the controller. The controlled variable is y(t)¼T(L,t) [1C]. Measurement
52
70
Controlled variable PI
51
60 Control variable
PID
Controlled variable PID
50
50
49 Set-point
48
40 PI
47
30
46 20
45
10
44 1000
1200
1400
1600
1800
1000
2000
1200
Control variable
Time [sec]
55
90
54
80
53 PID
51
50
50
40
49
30
48
20
47
10
46 400
500
600
1800
2000
Controlled variable PID
52
60
0 300
1600
Time [sec]
100
70
1400
700
800
900
45 300
Set-point
400
500
600
700
800
900
Time [sec]
Time [sec]
Fig. 7. Experimental results: (a) and (b) set-point and load step ( 20% change of the controller output at t¼ 1500 s) responses of the plant with the PIplant and PIDplant controllers; (c) and (d) disturbance responses of the plant with the PIDplant controller, for disturbance induced by activating/deactivating the fan. Control variable in [%], controlled variable in [1C].
Table B1 Stable processes (Ms ¼ 2, mn ¼ 2), kn [19]. kn (φ1, ρ) φ¼ 10 φ¼ 20 φ¼ 30 φ¼ 40 φ¼ 50 φ¼ 60 φ¼ 70 φ¼ 80 φ¼ 90 φ¼ 100 φ¼ 110 φ¼ 115 φ¼ 120 φ¼ 125 φ¼ 130
ρ ¼ 0.5
0.4977 0.5084 0.5198 0.5324 0.5393 0.5467 0.5547 0.5634
ρ ¼ 0.55
0.4968 0.5060 0.5159 0.5269 0.5330 0.5396 0.5468
ρ¼ 0.6
0.4891 0.4964 0.5043 0.5129 0.5227 0.5282 0.5341
ρ ¼0.65
0.4848 0.4903 0.4963 0.5031 0.5106 0.5194 0.5243 0.5297
ρ ¼0.7
0.4871 0.4914 0.4964 0.5022 0.5089 0.5167 0.5212
ρ¼ 0.75
0.4871 0.4893 0.4925 0.4967 0.5016 0.5075 0.5145 0.5186
ρ¼ 0.8
0.4910 0.4902 0.4913 0.4936 0.4969 0.5012 0.5064 0.5128
ρ ¼0.85
0.5000 0.4948 0.4929 0.4931 0.4996 0.4973 0.5008 0.5054 0.5113
ρ ¼ 0.9
ρ ¼ 0.95
0.5165 0.5046 0.4982 0.4954 0.4947 0.4956 0.4976 0.5006 0.5047
0.5446 0.5218 0.5086 0.5013 0.4976 0.4962 0.4965 0.4979 0.5004 0.5041
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Estimates in Table 3 are used to design the PIDplant controllers with parameters in Table 4. Parameters are calculated as presented in Appendix B, by using data from the three points in the look-up Tables B1–B3, determined previously for stable processes in [19]. These three points are defined by {ρ2,1 ¼0.90, φ1,1 ¼ 301}, {ρ2,1 ¼ 0.90, φ2,1 ¼401} and {ρ2,2 ¼0.95, φ2,2 ¼ 401}, since ρest and φest, presented in Table 3, correctly classify the laboratory plant in this stable region of the ρ φ plane. Points used to obtain the gains kni,j, kini,j, kdni,j, i,j ¼1,2, in Appendix B, are bolded in Tables B1–B3. The impact of reducing the number of points J on estimates kuest, ωuest, φest, Gpest(0) and on the PID controller tuning is illustrated in Tables 3 and 4. In Table 4, maximum sensitivity Ms is defined by M s ¼ max j1=ð1 þ CðsÞGm ðsÞÞjs ¼ iω , where Gm(s) is the ω
model (1) and (2), defined by parameters in Table 3, and C(s) is the PIDplant controller (4), defined by parameters in Table 4. In both cases, a) and b) in Table 4, the obtained maximum sensitivity is exactly equal to Ms ¼2, used to calculate Tables B1–B3 in [19]. Also, the obtained sensitivity to the high-frequency measurement noise Mn ¼|kd|/Tf [15] is exactly equal to the value mn|ku|, for mn ¼2 used to calculate Tables B1–B3 in [19]. Besides, in Table 4, insignificantly different values of the sensitivity to the high-frequency measurement noise Mn are obtained for different values of J. The anti-reset windup implementation of the PID controller (4) is applied, given by U C ðsÞ ¼ T aw
bks þ ki kd s2 þ ks þ ki 1 RðsÞ YðsÞ þ UðsÞ; T aw s þ 1 T aw s þ 1 ðT aw s þ1ÞðT f s þ 1Þ ð12Þ
7
as in [20], where the saturation element is defined by the input uC(t) and output u(t): 8 uC ðtÞ rllow l ; > < low u ðtÞ; l C ð13Þ uðtÞ ¼ low ouC ðtÞ o lhigh : > :l ; u C ðtÞ Zlhigh high In the linear region of the saturation element, for uC(t) u(t) one obtains (4) from (12). Closed-loop responses presented in Fig. 7a–d, were obtained with the PIDplant controller (12) and (13), having the gains k, ki, kd and noise filter time constant Tf presented in Table 4 for smaller value of J ¼20, and for b ¼0.4, llow ¼0 and lhigh ¼100%. Parameter Taw ¼37 s is obtained from Taw ¼ aTi þ(1 a)Td, for a ¼0.2 satisfying the condition Td oTaw oTi, for Ti ¼ k/ki and Td ¼kd/k. Advantages of the PIDplant controller, compared to the PIplant controller used to determine kuest, ωuest, φest and Gpest(0) in Table 3, are demonstrated in Fig. 7a–b. This experiment starts from the temperature T(L,t) E45 1C, for to1000 s. Then, at t ¼1000 s the set-point is changed for r0 ¼5 1C. At t¼ 1500 s a load disturbance is inserted as a 20% step change of the controller output. As presented in Fig. 7b, improvement of the performance is obtained by the PIDplant controller. Closed-loop experiment in Fig. 7c–d starts from the temperature T(L,t)E50 1C for to300 s, by activating a fan at t¼300 s. Then, at t¼ 600 s the fan is switched-off. Action of the fan induced an intensive cooling of the plate, resulting into a strong disturbance as seen from the control signal u(t) in Fig. 7c. Anti-reset windup action is activated two times, around 320 s and 620 s.
Table B2 Stable processes (Ms ¼ 2, mn ¼ 2), kin [19]. kin (φ1, ρ) φ ¼10 φ ¼20 φ ¼30 φ ¼40 φ ¼50 φ ¼60 φ ¼70 φ ¼80 φ ¼90 φ ¼100 φ ¼110 φ ¼115 φ ¼120 φ ¼125 φ ¼130
ρ¼ 0.5
0.3393 0.3250 0.3119 0.3004 0.2952 0.2905 0.2861 0.2822
ρ ¼0.55
0.2931 0.2803 0.2688 0.2587 0.2542 0.2500 0.2462
ρ¼ 0.6
ρ ¼ 0.65
ρ ¼ 0.7
ρ¼ 0.75
0.2666 0.2543 0.2430 0.2328 0.2239 0.2199 0.2163
0.2441 0.2324 0.2214 0.2113 0.2023 0.1945 0.1910 0.1879
0.2137 0.2032 0.1933 0.1842 0.1763 0.1694 0.1664
0.1976 0.1875 0.1778 0.1689 0.1609 0.1538 0.1478 0.1451
ρ ¼0.6
ρ ¼ 0.65
ρ ¼0.7
ρ ¼ 0.75
ρ ¼ 0.8
ρ ¼ 0.85
0.1835 0.1738 0.1645 0.1557 0.1477 0.1406 0.1343 0.1290
0.1716 0.1620 0.1530 0.1443 0.1364 0.1292 0.1228 0.1173 0.1126
ρ ¼0.8
ρ¼ 0.85
ρ ¼0.9
ρ ¼ 0.95
0.1622 0.1522 0.1430 0.1345 0.1266 0.1194 0.1129 0.1072 0.1023
0.1558 0.1446 0.1349 0.1261 0.1181 0.1108 0.1043 0.0985 0.0935 0.0892
ρ ¼ 0.9
ρ ¼ 0.95
0.8286 0.6798 0.5566 0.4515 0.3591 0.2758 0.1987 0.1257 0.0552
0.9886 0.8040 0.6569 0.5332 0.4314 0.3403 0.2582 0.1823 0.1104 0.0409
Table B3 Stable processes (Ms ¼ 2, mn ¼ 2), kdn [19]. kdn (φ1, ρ) φ ¼10 φ ¼20 φ ¼30 φ ¼40 φ ¼50 φ ¼60 φ ¼70 φ ¼80 φ ¼90 φ ¼100 φ ¼110 φ ¼115 φ ¼120 φ ¼125 φ ¼130
ρ ¼0.5
0.4412 0.3556 0.2722 0.1899 0.1490 0.1082 0.0675 0.0270
ρ ¼ 0.55
0.3934 0.3099 0.2287 0.1488 0.1090 0.0694 0.0299
0.4383 0.3527 0.2712 0.1920 0.1141 0.0754 0.0368
0.4919 0.4018 0.3178 0.2380 0.1607 0.0846 0.0468 0.0091
0.4586 0.3700 0.2876 0.2093 0.1336 0.0591 0.0221
0.5256 0.4293 0.3420 0.2611 0.1843 0.1100 0.0370 0.0008
0.6061 0.4982 0.4032 0.3173 0.2377 0.1623 0.0894 0.0178
0.7047 0.5801 0.4736 0.3800 0.2954 0.2171 0.1429 0.0712 0.0009
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
8
The analysis of the same laboratory thermal plant in [20], and PID controller tuning, was based on the 5th-order model, obtained by the open-loop 100th-order ARX identification, in the time interval Texp ¼1000 s, and model reduction. The estimates kuest, ωuest, φest and Gpest(0) in (Table 3, row a), are obtained within the shorter time interval Texp ¼410 s in Fig. 2, defined by (6) as Texp ¼(TJþ T/2). Finally, it is important to observe that the estimates kuest, ωuest, φest and Gpest(0) are obtained for Texp E3Tu, since the ultimate period Tu ¼2π/ωuest is equal to Tu ¼141 s for ωuest in (Table 3, row a). 5. Conclusions The proposed procedure for determining the model Gm(s) of stable processes, processes with oscillatory dynamics, integrating and unstable processes, including dead-time, by estimating the quadruplet {ku, ωu, φ, Gp(0)} has been performed in an extremely short time interval. The proposed fast process dynamics characterization, performed without breaking the loop containing a controller in operation, offers the following: (i) first, as demonstrated by experimental results, the controller tuning for the desired performance/robustness tradeoff is obtained by using the estimated kuest, ωuest, φest and Gpest(0), and previously calculated and memorized process-independent look-up tables of parameters kn, kin and kdn, as demonstrated in Fig. 7a–d; (ii) a reliable simulation performance analysis of the loop considered is possible, as demonstrated in Fig. 4. Acknowledgements Authors gratefully acknowledge Dr. Ribić his help in implementing anti-windup PID controller on the laboratory plant. T.B. Šekara gratefully acknowledges the financial support from the Ministry of Education, Science and Technological Development of Serbia (Project TR33020). Appendix A. Determination of φest and Aest by the PLL and SheMa estimators When applying PLL and SheMa estimators, kuest and ωuest are determined by the first experiment. As demonstrated in [5], two additional experiments must be performed, to obtain for ϕ 7 ¼ 7 π=N 0 ; N 0 36: ωu7 , defined by argðGp ðiωu7 ÞÞ ¼ π þ ϕ 7 , 7 7 ρ 7 ¼ ku =ku and ku ¼ jGp ðiωu7 j 1 . Then, φest and Aest are determined from relations [14,15]: φþ þφ ρ 7 sin ϕ 7 A þ þA0 ; Aest ¼ 0 ; φ 7 ¼ arctan 7 ; φest ¼ 7 2 2 ρ cos ϕ 1 2j ωu ωu7 j A07 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 1 2ρ 7 cos ϕ 7 þðρ 7 Þ2
ðA1Þ
Appendix B. PID controller tuning by applying look-up tables in the ρ φ plane The process independent normalized gains kn ¼kn(ρ,φ), kin ¼kin (ρ,φ), kdn ¼kdn(ρ,φ) of a virtual PIDn controller, calculated in [19, Tables A1–A3], are presented here as look-up Tables B1–B3. The gains k, ki, kd and the noise filter time constant Tf, in the real
PID controller, are defined by [19] k ¼ kuest kn ; ki ¼ kuest ωuest kin ; kd ¼ kuest kdn =ωuest ; T f ¼ T fn =ωuest ;
ðB1Þ
where the normalized noise filter time constant is defined as Tfn ¼| kdn(ρ,φ)|/mn. From (B1) it follows that the real sensitivity to the high-frequency measurement noise Mn ¼|kd|/Tf, used in [15], is given by Mn ¼ mn|ku|. Results in Tables B1–B3 were obtained in [19] by optimization under constraints on the maximum sensitivity Ms ¼2 and for the normalized sensitivity to measurement noise mn ¼2. However, optimization can be performed under desired constraints on the Ms, mn and for the desired closed-loop system damping ratio ζ ¼1, to obtain the critically damped set-point and load disturbance step responses [15], as done in [19] for processes with oscillatory dynamics. Parameters kn, kin and kdn of the virtual PIDn controller are obtained by using the gains kni,j, kini,j, kdni,j, i,j ¼1,2, from the lookup Tables B1–B3, and by applying interpolation from the three points [5], as defined by kn ¼ ð1 α βÞkn2;1 þ αkn2;2 þ βkn1;1 , kin ¼ ð1 α βÞkin2;1 þ αkin2;2 þ βkin1;1 , kdn ¼ ð1 α βÞkdn2;1 þ αkdn2;2 þ βkdn1;1 , α¼(ρest–ρ2,1)/(ρ2,2 ρ2,1) and β¼(φ2,1 φest)/(φ2,1 φ1,1). Angles in Tables B1–B3 are in degrees and φ1,1 rφ2,1, φ1,1 rφest rφ2,1, ρ1,1 rρest r ρ2,1. For details see [5]. References [1] Hjalmarsson H. From experiment design to closed-loop control. Automatica 2005;41:393–438. [2] Álvarez JD, Guzmán JL, Rivera DE, Berenguel M, Dormido S. Perspectives on control relevant identification through the use of interactive tools. Control Eng Pract 2013;21:171–83. [3] Wahlberg B. Model reductions of high-order estimated models: the asymptotic ML approach. Int. J Control 1989;49:169–92. [4] Ziegler JG, Nichols NB. Optimum settings for automatic controllers. Trans ASME 1942;64:759–68. [5] Šekara TB, Mataušek MR. PID controller tuning based on the classification of stable, integrating and unstable processes in a parameter plane. In: Serra. GLO, editor. Frontiers in advanced control systems. InTechOpen; 2012. p. 117–42. [6] Cohen GH, Coon GA. Theoretical consideration of retarded control. Trans ASME 1953;75:827–34. [7] Mataušek MR, Ribić AI. Design and robust tuning of control scheme based on the PD controller plus Disturbance Observer and low-order integrating firstorder plus dead-time model. ISA Trans 2009;48:410–6. [8] O’Dwyer A. Handbook of PI and PID controller tuning rules. 3rd ed.Imperial College Press; 2009. [9] Lo WL, Rad AB, Li CK. Self-tuning control of systems with unknown time delay via extended polynomial identification. ISA Trans 2003;42:259–72. [10] Mei H, Li S. Decentralized identification for multivariable integrating processes with time delays from closed-loop step tests. ISA Trans 2007;46:189–98. [11] Hashimoto T, Ishida Y. An adaptive I-PD controller based on frequency domain system identification. ISA Trans 2000;39:71–8. [12] Åström KJ, Panagopoulos H, Hägglund T. Design of PI controllers based on non-convex optimization. Automatica 1998;34:585–601. [13] Panagopoulos H, Åström KJ, Hägglund T. Design of PID controllers based on constrained optimization. IEE Proc. Control Theory Appl 2002;34:585–601. [14] Šekara TB, Mataušek MR. Revisiting the Ziegler–Nichols process dynamics characterization. J Process Control 2010;20:360–3. [15] Mataušek MR, Šekara TB. PID controller frequency-domain tuning for stable, integrating and unstable processes, including dead-time. J Process Control 2011;21:17–27. [16] Åström KJ, Hägglund T. Automatic tuning of simple regulators with specifications on phase and amplitude margins. Automatica 1984;20:645–51. [17] Åström KJ, Hägglund T. PID controllers: theory, design and tuning. 2nd ed.. Research Triangle Park, NC: ISA; 1995; 27709. [18] Šekara TB, Mataušek MR. Relay-based critical point estimation of a process with the PID controller in the loop. Automatica 2011;47:1084–8. [19] Šekara TB, Mataušek MR. Classification of dynamic processes and PID controller tuning in a parameter plane. J Process Control 2011;21:620–6. [20] Ribić AI, Mataušek MR. A dead-time compensating PID controller structure and robust tuning. J Process Control 2012;22:1340–9. [21] S Yamamoto, I Hashimoto. Present status and future needs: the view from Japanese industry. In: Proceedings of the IV conference on chemical process control, Texas, USA, 1991.
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
Contents lists available at ScienceDirect
ISA Transactions journal homepage: www.elsevier.com/locate/isatrans
Research Article
A fast closed-loop process dynamics characterization Miroslav R. Mataušek n, Tomislav B. Šekara Faculty of Electrical Engineering, University of Belgrade, Belgrade 11120, Serbia
art ic l e i nf o
a b s t r a c t
Article history: Received 17 January 2013 Received in revised form 9 October 2013 Accepted 3 December 2013
Stable, integrating and unstable processes, including dead-time, are analyzed in the loop with a known PI/PID controller. The ultimate gain and frequency of an unknown process Gp(s), and the angle of tangent to the Nyquist curve Gp(iω) at the ultimate frequency, are determined from the estimated Laplace transform of the set-point step response of amplitude r0. Gain Gp(0) is determined from the measurements of the control variable and known r0. These estimates define a control relevant model Gm(s), making possible the use of the previously determined and memorized look-up tables to obtain PID controller guaranteeing desired maximum sensitivity and desired sensitivity to measurement noise. Simulation and experimental results, from a laboratory thermal plant, are used to demonstrate the effectiveness and merits of the proposed method. & 2013 ISA. Published by Elsevier Ltd. All rights reserved.
Keywords: PID control Closed-loop identification Measurement noise Robustness Tuning
1. Introduction The importance of control relevant process dynamics characterization is discussed in detail in an overview paper [1], and reexamined recently in [2]. One of the basic conclusions from [1,2] follows the idea from [3], where it is proposed to estimate a highorder ARX model first and then to perform a model reduction in the frequency domain, to obtain a reduced-order model used for controller tuning. Three problems are related to the adequate control relevant process dynamics characterization: the model structure, the control relevant region of frequencies and the dilemma open-loop versus closed-loop process identification. Closed-loop system performance/robustness tradeoff strongly depends on a priory knowledge defined by the structure of the model used for the PID controller tuning. The two parameter model, obtained from the Ziegler–Nichols time domain tuning [4], can be used to tune PID controller for lag-dominated processes Gp(s) [5]. This model is known at the present time as integrator plus dead-time (IPDT) model. Its relationship with the Ziegler– Nichols frequency domain tuning is discussed in detail in [5]. The second one, used for the controller tuning, is the three parameter model. It is represented by the first-order plus dead-time (FOPDT) model, introduced first by Cohen–Coon [6], and by the integrating first-order plus dead-time (IFOPDT) model. FOPDT model can be used to approximate balanced and dead-time dominated processes, while IFOPDT model defines a better approximation of
n
Corresponding author. Tel.: þ 381 11 337 01 64; fax: þ381 11 324 86 81. E-mail addresses: [email protected] (M.R. Mataušek), [email protected] (T.B. Šekara).
the lag-dominated processes, than the IPDT model [7]. The four parameter second-order plus dead-time (SOPDT) model, and other reduced-order models, used for the PI/PID controller tuning, are summarized in [8]. Mathematical models are always an approximation of reality. Appearance of a dead-time in a model might be a consequence of low-order modeling [9,10], for example based on the FOPDT, IFOPDT or SOPDT models, or might be a consequence of the real time-delay caused by some physical phenomena [11]. The control relevant region of frequencies is the region around the ultimate frequency ωu of a process Gp(s). This is confirmed by the PI/PID controller optimization [12–15], based on the frequency response of the process Gp(iω), under constraints on the robustness. The frequency ω0, where the sensitivity function has its maximum, occurs in the region around the ultimate frequency ωu, as demonstrated in [14]. The importance of the process dynamics characterization based on the ultimate frequency estimation is firstly recognized by Ziegler and Nichols [4] and further developed by Åström and Hägglund [16,17]. The extension of the Ziegler– Nichols process dynamics characterization proposed in [14], by introducing an additional parameter in the frequency domain, improved considerably the possibility of better process modeling in the wider region of frequencies around the ultimate frequency ωu of a large class of stable processes, processes with oscillatory dynamics, integrating processes and unstable processes Gp(s), including dead-time. Dilemma open-loop versus closed-loop process identification in academia is treated mainly as a problem of ensuring the best statistical accuracy [1,3]. However, in industry, breaking of control loops in operation is mainly ignored by plant operators. When some initially tuned controller is in operation, a procedure for fine tuning can be easily accepted if it can be activated/deactivated
0019-0578/$ - see front matter & 2013 ISA. Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.isatra.2013.12.006
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
2
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
without breaking control loops in operation [18]. A special attention should be devoted when estimating parameters of a continuous model including dead-time. For example, this can be done by applying the extended polynomial closed-loop identification [9], or to obtain model parameters, including dead-time, from the estimated frequency response of the process, obtained from the closed-loop step response tests [10,11]. However, in these methods, the dead-time free part of the model transfer function must be specified, based on some a priory knowledge, and different model structure is used for stable, integrating and unstable processes. In this paper, a new insight into the problem of control relevant process dynamics characterization and an effective solution of this problem is presented. It is based on a continuous model Gm(s) with the unified structure for a large class of stable processes, processes with oscillatory dynamics, integrating processes and unstable processes Gp(s), including dead-time. Proposed in [14] as an effective extension of the Ziegler–Nichols process dynamics characterization in the frequency domain, this model is defined by the quadruplet {ku, ωu, φ, A}. Parameter ku is the ultimate gain and φ ¼ argð∂Gp ðiωÞ=∂ωÞjω ¼ ωu is the angle of the tangent to the Nyquist curve Gp(iω) at the ultimate frequency ωu of a process Gp(s). In this model an equivalent dead-time τ is defined by τ¼ φ/ωu. Two procedures [15,18] are proposed for estimating the quadruplet {kuest, ωuest, φest, Aest}. Parameter Aest ¼A0 is defined in the frequency domain by A0 ¼ 2j∂Gp ðiωÞ=∂ωjω¼1 ωu =ku [15]. The new PLL (phase-locked-loop) estimator [15], further improved in [5], requires only that the controller in operation is a linear controller, while the relay SheMa estimator [18] does not require this preliminary information. Estimation methods [5,15,18] are performed without breaking the loop containing a controller in operation. The model Gm(s), defined by the quadruplet {kuest, ωuest, φest, Aest}, approximates the Nyquist curve of a large class of stable processes, processes with oscillatory dynamics, integrating and unstable processes Gp(s), including dead-time, in a large region around ωu, and can be effectively used in PID controller constrained optimization [5]. However, from the industry viewpoint, disadvantage of estimation methods [5,15,18] can be a longer period of time required for determining the quadruplet {kuest, ωuest, φest, Aest}. Namely, determination of φest and Aest by applying PLL and SheMa estimators requires two additional experiments, as explained in Appendix A. Besides, and this is of essential importance, the basic definition of parameter A¼ ωukuGp(0)/(1 þkuGp(0)) [14] offers the possibility to classify a large class of stable processes, processes with oscillatory dynamics, integrating and unstable processes, including deadtime, in the ρ φ plane [19], where ρ ¼κ/(1þ κ), κ ¼kuGp(0). Stable and unstable processes are classifiedpas processes inside and ffiffiffiffiffiffiffiffiffiffi outside the region 0 oρo1, 0 o φ oπ= ρ þ1, while the integratpffiffiffi ing processes are classified as processes with ρ¼ 1, 0 o φ o π= 2, since Gp(0)¼ 71. The possibility to classify stable processes, processes with oscillatory dynamics, integrating, and unstable processes, including deadtime, in a two parameter plane, is important from the PID controller tuning viewpoint. The desired performance/robustness tradeoff can be obtained by applying the previously memorized process indepen-
Fig. 1. Process Gp(s), with a two-degree-of-freedom controller. The set-point r(t), controlled variable y(t), control variable u(t), load and output disturbances, d(t) and n(t), represent variations around their values in the nominal regime.
dent look-up tables in this ρ φ plane [19]. However, to be effective, such a possibility must be supported by a fast estimation of the quadruplet {ku, ωu, φ, Gp(0)}. Then, from the estimated {ku, ωu, φ, Gp(0)} and the previously memorized look-up tables one obtains directly the gains and the noise filter time constant of a real parallel PID controller, as presented in Appendix B, for stable processes and illustrated by experimental results. In the present paper a new simple and effective closed-loop procedure is proposed for estimating quadruplet {ku, ωu, φ, Gp(0)} in a short time interval, without braking the loop with the controller in operation. It is assumed that the controller in operation is a known linear controller. The estimation procedure is defined in Section 2. In Section 3, a test batch, consisting of stable processes, process with oscillatory dynamics, integrating and unstable processes, including dead-time, in the loop with PI and PID controllers, is used to demonstrate the properties of the proposed estimation method. Finally, in Section 4, the experimental verification of the proposed method is presented. The quadruplet {kuest, ωuest, φest, Gpest(0)} is determined by applying a PI controller to a laboratory thermal process [20], with noisy measurements, and used for the PID controller tuning by applying the previously memorized look-up tables in the ρ φ plane, presented in Appendix B.
2. Determination of the quadruplet {ku, ωu, φ, Gp(0)} from the closed-loop system set-point step response Model Gm(s) of an unknown stable process, process with oscillatory dynamics, integrating process and unstable process, 7 6 5 4 3 2 1 0
0
100
200 300 Time [sec]
400
500
Fig. 2. The set-point (dashed), the measured closed-loop response y(mTs) of the real plant (solid-blue), the data yj (circles-black) and yid(t), defined by (solid-red) lines for 0r tr 400 s, obtained by linear interpolation of yj. For t4400 s the response is approximated with its value defined by the known set-point r0 ¼ 5 and presented also by solid-red line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Table 1 Parameters of the controllers used to obtain closed-loop set-point step responses for processes Gpl(s), l ¼1,2,…,7. Process/controller
k
ki
kd
Tf
b
Gp1(s)/PI Gp2(s)/PI Gp3(s)/PI Gp4(s)/PID1 Gp5(s)/PI Gp6(s)/PI Gp7(s)/PI
2.1631 0.7903 0.1204 0.5620 0.1010 1.3273 0.3553
4.7980 0.0654 0.0946 0.0830 0.00255 0.01354 0.0030
0 0 0 1.0620 0 0 0
0 0 0 0.1770 0 0 0
1 1 1 0 1 0.15 0.25
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
From Fig. 1, for d ¼0 and n ¼ 0, one obtains that process transfer function Gp(s) is defined by
including dead-time, proposed in [14], is defined by relations Gm ðsÞ ¼
WðsÞ 1 Aωu expð τsÞ ; WðsÞ ¼ ; τ ¼ φ=ωu ; 1 WðsÞ ku s2 þ ω2u
ð1Þ Gp ðsÞ ¼
A ¼ ρωu ; ρ ¼ ku Gp ð0Þ=ð1 þku Gp ð0ÞÞ;
3
ð2Þ
HðsÞ YðsÞ ; ; HðsÞ ¼ C ff ðsÞ HðsÞCðsÞ RðsÞ
ð3Þ
where H(s) is the closed-loop system transfer function from r to y in Fig. 1. In the present paper, the step set-point R(s) ¼r0/s, of known amplitude r0, is used to collect measurements of the controlled variable applied to obtain an estimate Yest(s). Then, for
and by the quadruplet {ku, ωu, φ, Gp(0)}. Unknown process Gp(s), in the loop with known linear controllers C(s) and Cff(s), is presented in Fig. 1.
Table 2 True values ku, ωu, φ, Gp(0) compared with the estimated kuest, ωuest, φest, Gpest(0) of processes Gpl(s), l ¼1,2,…,7, obtained for the three pairs T1J1 ¼ T2J2 ¼ T3J3, T1 o T2 o T3. Time period T is in seconds. Process
ku
ωu
φ
Gp1(s)
11.5919
9.8696
0.7854
1
Gp2(s)
4.6626
0.2144
0.8271
1
Gp3(s)
2.5443
0.5884
0.5648
1
Gp4(s)
0.5640
0.4080
0.7200
1
Gp5(s)
0.2371
0.2291
0.9716
1
Gp6(s)
3.1865
0.4287
0.3903
1
Gp7(s)
0.6341
0.5828
0.7603
4
Gp(0)
kuest
ωuest
φest
Gpest(0)
T
J
11.660 11.791 12.335 4.7021 4.7675 4.8814 2.5448 2.6287 2.9567 0.5641 0.5673 0.5731 0.2380 0.2405 0.2514 3.1901 3.2015 3.2307 0.6358 0.6415 0.6460
9.8817 9.8796 9.8628 0.2143 0.2142 0.2140 0.5870 0.5881 0.5850 0.4085 0.4100 0.4127 0.2291 0.2289 0.2285 0.4280 0.4260 0.4226 0.5836 0.5835 0.5835
0.8055 0.8032 0.7659 0.7828 0.7695 0.7480 0.6074 0.6244 0.5900 0.7167 0.7125 0.7047 0.9628 0.9528 0.9083 0.3920 0.3937 0.5632 0.8289 0.8193 0.8105
1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1 1 1 1 1 1 1.000 1.000 1.000 4.005 4.005 4.005
0.02 0.04 0.08 1 2 3 0.5 1 2 0.5 1 1.5 1 2 4 0.8 1 2 0.4 0.8 1
200 100 50 150 75 50 200 100 50 120 60 40 200 100 50 200 160 80 200 100 80
0.15
0.03 0.02
0.1
0.01 0.05
0 -0.01
0
-0.02
-0.05
-0.03
-0.1
Gm1(s) Gp1 (s)
-0.04 -0.05 -0.06 -0.12
-0.11
-0.1
-0.09
-0.08
-0.07
-0.06
Gm3 (s) Gp3 (s)
-0.15
-0.05
4
-0.2 -0.6
-0.5
-0.4
-0.3
-0.2
0.015 0.01
2
0.005 0
0 -0.005
-2
-0.01
-4
-8
-0.015
Gm5 (s) Gp5 (s)
-6
-7
-6
-5
-4
-3
-2
-1
0
Gp6 (s) Gm6 (s)
-0.02
1
-0.025 -0.38
-0.36
-0.34
-0.32
-0.3
-0.28
-0.26
Fig. 3. Nyquist curves, around the critical point ( 1/ku,i0): (a) stable process Gp1(s); (b) process with oscillatory dynamics Gp3(s); (c) integrating process Gp5(s); and (d) unstable process Gp6(s). In all cases: exact (solid-blue), estimated (dashed-red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
4
H(s) ¼sYest(s)/r0 and for the known linear controllers C(s) and Cff(s), from (3) one obtains an estimate Gpest(s), used to determine parameters kuest, ωuest and φest. In industry applications PI controller predominates [21]. Thus, in the present paper it is assumed that the two-degree-of-freedom PI controller, or parallel PID controller is used to determine the quadruplet {ku, ωu, φ, Gp(0)}. In (3) and Fig. 1, these controllers are given by C(s)¼ (kds2 þks þ ki)/(s(Tfsþ 1)) and Cff(s)¼(bksþki)/s. They are implemented as defined by relations k UðsÞ ¼ ðbkRðsÞ kY f ðsÞÞ þ i ðRðsÞ Y f ðsÞÞ kd sY f ðsÞ; 0 rb; Y f ðsÞ s YðsÞ ; ¼ Tf sþ1
ð4Þ
used also in [13,15]. For the PI controller kd ¼0, and it is adopted that Tf ¼0 as in [12,15]. The proposed estimation of quadruplet {kuest, ωuest, φest, Gp(0)} is performed in three steps. In the first step, a small sampling period Ts is used to collect the closed-loop system set-point step response y(mTs) and the control signal u(mTs) m ¼ 0,1,2,…,M. To reduce the effect of measurement noise, the averaging procedure is applied to obtain the signal yj defined by yj ¼
1 jn 1 yðmT s Þ; j ¼ 1; 2; …; J; y0 ¼ 0; yJ þ 1 ¼ r 0 : ∑ n m ¼ ðj 1Þn
ð5Þ
In (5) parameter n, significantly greater than 1, is defined by period T ¼nTs. The time period T is defined as Tp ¼pT, where Tp is the time-to-go from 10% of the set-point step response first maximum to the first maximum. Alternatively, if the well damped closed-loop system set-point step response is obtained by the implemented PI/PID controller, it is recommended to define Tp as the rise time. For integer p it is recommended to use p¼ 5–20. Then, signal yid(t), used to determine the Laplace transform Yest(s), is obtained by applying the linear interpolation of signal yj, j¼1,2, …,J. For j ZJþ 1 it is adopted yj ¼r0. The reason is obvious. The relation lim yid ðtÞ ¼ r 0 holds true for the closed-loop system t-1 set-point step response of a process in the loop with the PI/PID controller. This first step is illustrated in Fig. 2. The measured closed-loop set-point step response y(mTs), obtained experimentally for Ts ¼0.5 s on the real thermal plant from Section 4, is presented in Fig. 2 with data yj and signal yid(t), obtained by linear interpolation of yj. For the time time-to-go to the first maximum equal to Tp E140 s and for p¼ 7 one obtains T¼ 20 s. A small number of points J¼ 20 is used to calculate Yest(s). The response is truncated at t¼ 400 s to demonstrate that the extremely short time interval JTþ T/2 is necessary for obtaining adequate estimates
kuest, ωuest, φest and Gpest(0), guaranteeing an excellent tuning of the PID controller, as demonstrated in Section 4. In the second step the Laplace transform Yest(s) is determined, by applying data yid(t) prepared in the first step of the proposed estimation procedures. Taking into account that for j ZJ þ1 it is adopted yj ¼r0, the Laplace transform Yest(s) is defined by the relation Z Y est ðsÞ ¼
JT þ T=2 0
yid ðtÞe st dt þ
r 0 sðJT þ T=2Þ e ; s
ð6Þ
where the time JTþ T/2 corresponds to the settling time estimated from the memorized closed-loop system set-point step response. Data yj, j¼ 1,2,…,J, used to define signal yid(t) by the linear interpolation, appears at points jT T/2. Thus, from (6) one obtains 2y Y est ðsÞ ¼ 1 T þ
Z
T=2
te
st
J
dt þ ∑
j¼1
0
Z
jT þ T=2 jT T=2
! ðaj þ bj tÞe
st
dt
r 0 sðJT þ T=2Þ e ; s
ð7Þ
where aj ¼yj (yj þ 1 yj)(j 1/2) and bj ¼ (yj þ 1 yj)/T. By solving integrals in (7), the analytical solution of the Laplace transform Yest(s) is obtained, defined by the measured data and following relations: Y est ðsÞ ¼ 2y1
1 ð1 þsT=2Þe sT=2 r0 þ Y 1est ðsÞ Y 2est ðsÞ þ e sðJT þ T=2Þ ; s s2 T ð8Þ
ðsTyj þyj þ 1 yj Þe sðjT T=2Þ ; s2 T j¼1 J
ð9Þ
Y 1est ðsÞ ¼ ∑
Fig. 5. The laboratory thermal plant. Positions of temperature sensors are denoted by 2 and 3, while position of the heater is denoted by 1.
0.1 0 Load disturbance response
Set-point response
1 0.8 0.6 0.4 0.2 0
0
10
20
30 Time [sec]
40
50
-0.1 -0.2 -0.3 -0.4 -0.5 -0.6
0
10
20
30
40
50
Time [sec]
Fig. 4. Process Gp4(s) and its model Gm4(s), defined by kuest ¼0.5731, ωuest ¼0.4127, φest ¼0.7047 and Gpest(0)¼ 1, in the loop with PIDtun controller. Responses of the controlled variable for b¼ 0: (a) set-point step response, R(s) ¼ 1/s; (b) load disturbance step response, D(s) ¼ 0.2/s.
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
ðsTyj þ 1 þ yj þ 1 yj Þe sðjT þ T=2Þ : s2 T j¼1 J
ð10Þ
Y 2est ðsÞ ¼ ∑
It can be easily confirmed that the condition limHðsÞ ¼ 1 is s-0 satisfied, for the closed-loop system transfer function defined by H(s) ¼sYest(s)/r0, with Yest(s) given by (8)–(10). A time scaled set-point response is used to perform calculations, with the time scale factor defined by JT. If higher values of Ts are used for memorizing the set-point step response, the following improvement is proposed. Instead from (5), the signal yj, j¼1,2,…, R jT J, can be obtained from the integral yj ¼ ð ðj 1ÞT ym ðtÞdtÞ=T, where ym(t) is defined by the linear interpolation of y(mTs). In the third step, from H(s)¼sYest(s)/r0 and (3), the process transfer function Gpest(s) is determined and used to estimate the ultimate frequency ωuest from argðGpest ðiωuest ÞÞ ¼ π. Then, the ultimate gain kuest is obtained from kuest ¼ jGpest ðiωuest Þj 1 and the angle φest is obtained from φest ¼ argð∂Gpest ðiωÞ=∂ωÞjω ¼ ωuest , as defined in [14]. The gain Gp(0) is estimated from 8 > Z JT þ T=2 < r 0 =u0 ; for u0 a 0; 1 for u0 -0; kuest 4 0; um ðtÞdt; Gpest ð0Þ ¼ 1; u0 ¼ > T JT T=2 : 1; for u -0; k o 0: uest
0
ð11Þ by using the known amplitude r0 of the set-point response and linear interpolation of the memorized control signal u(mTs), m ¼0,1,2,…,M, used to calculate um(t) in (11).
3. Simulation analysis of the proposed estimation method Effectiveness of the proposed method is demonstrated first by the following test batch: 1 ð2s þ 1Þe 4s pffiffiffiffiffi; Gp2 ðsÞ ¼ ; Gp1 ðsÞ ¼ ð10s þ 1Þð7s þ 1Þð3s þ1Þ cosh 2s Gp3 ðsÞ ¼
es e 0:5s ; Gp4 ðsÞ ¼ 9s2 þ 2:4s þ 1 sð1:2s þ 1Þ3
Gp5 ðsÞ ¼
e 5s ; sðs þ 1Þð0:5s þ 1Þð0:25s þ 1Þð0:125s þ 1Þ
Gp6 ðsÞ ¼
e 0:5s 4e 2s ; Gp7 ðsÞ ¼ ð5s 1Þð2s þ 1Þð0:5s þ 1Þ 4s 1
5
The impact of reducing the number of points J on estimation is analyzed for the three pairs T1J1 ¼T2J2 ¼ T3J3, such that T1 oT2 oT3. As demonstrated in Table 2, the impact of reducing considerably number of points J is small in all examples, even for unstable processes. In Fig. 3, Nyquist curves of processes Gpl(s), l ¼1,3,5,6, are compared for T1 and J1 with the corresponding Nyquist curves of the estimated Gpestl(s). Excellent estimation of the Nyquist curves, presented in Fig. 3, is obtained within a large region around the critical point ( 1/ku,i0). This is important characteristic of the proposed identification method, since only this result makes possible to determine parameter φ. An important quality of the proposed identification method is demonstrated in Fig. 4. The closed-loop set-point and load disturbance step responses of Gp4(s) and Gm4(s), defined by kuest, ωuest, φest and Gpest(0) in Table 2 (Gp4, T¼1.5, J¼ 40), in the loop with the PIDtun controller, are almost exactly equal. Parameters of the PIDtun controller in (4) are obtained by applying kuest, ωuest, φest and Gpest(0) from Table 2 (Gp4, T¼1.5, J¼ 40), and tuning formulae defined by equations (30)–(37) in [15], to obtain critically damped responses for: k¼0.3787, ki ¼ 0.0324, kd ¼ 0.8453, Tf ¼0.2423, b¼ 0. Results presented in Fig. 4 confirm that the performance of the considered loop can be analyzed by applying the proposed identification method and simulation of the tuned controller in the loop with the model Gm(s), defined by (1) and (2), instead of the real process.
Table 3 Estimated values kuest, ωuest, φest, Gpest(0), obtained for two values of J from the real plant closed-loop set-point step response of amplitude r0 ¼ 5 1C, presented in Fig. 6b. Plant
kuest
ωuest
φest
Gpest(0)
ρest
Ts [s]
T [s]
J
(a) (b)
28.6582 27.0675
0.04458 0.04461
0.6377 0.6814
0.4104 0.4082
0.9216 0.9170
0.5 0.5
20 10
20 40
Table 4 Parameters of the PID controllers obtained, for the laboratory plant quadruplet {kuest, ωuest, φest, Gpest(0)} of Table 3, by applying tuning defined by look-up Tables B1–B3 of Appendix B.
of processes in the loop with the PI/PID controllers (4), defined by Table 1. The estimates kuest, ωuest, φest and Gpest(0) are compared in Table 2 with the true values of ku, ωu, φ and Gp(0).
PIDplant
K
ki
kd
Tf
J
Ms
Mn
(a) (b)
14.3795 13.5302
0.1774 0.1668
378.7222 340.0592
6.6076 6.2817
20 40
2.00 2.00
57.32 54.13
65
52
60
51
55
50
Controlled variable
Control variable
Set-point
50
49
45 48 40 47
35
46
30
45
25 20
44 1000
1200
1400
1600
Time [sec]
1800
1000
1200
1400
1600
1800
Time [sec]
Fig. 6. Experimental results: (a) and (b) set-point step responses of the plant with PIplant controller. Response in (b) is used to calculate Yest(s), as defined by Fig. 2. Control variable in [%], controlled variable in [1C].
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
6
of temperature at x¼ 0 is used to keep the temperature T(0,t) r70 1C [20]. The closed-loop set-point step responses in Fig. 6a–b are obtained on the laboratory thermal plant with the PIplant controller (4), defined by: k ¼6.3810, ki ¼0.055 and b¼1. The step response shown in Fig. 6b is used (as defined by Fig. 2) to obtain Yest(s) from (8)–(10) and parameters kuest, ωuest, φest in Table 3. The static gain Gpest(0) in Table 3 is obtained from (11) and the response shown Fig. 6a. By visual inspection one obtains Gp(0)E 0.4 1C/(%) from Fig. 6, which is in accordance with the estimate Gpest(0), given in Table 3.
4. Experimental verification by a laboratory thermal plant The laboratory thermal plant, presented in Fig. 5, consists of a thin plate made of aluminum, long L ¼0.1 m and wide h ¼0.03 m [20]. The temperature T(x,t) is distributed along the plate, from x ¼0 to x ¼L, and measured by precision sensors LM35 (TO92), at positions x ¼0 and x ¼L. The plate is heated by terminal adjustable regulator LM317 (TO 220) at position x ¼0. The manipulated variable is dissipated power of the heater at x¼ 0. The input to heater is control variable u(t) [%], defined by the output of the controller. The controlled variable is y(t)¼T(L,t) [1C]. Measurement
52
70
Controlled variable PI
51
60 Control variable
PID
Controlled variable PID
50
50
49 Set-point
48
40 PI
47
30
46 20
45
10
44 1000
1200
1400
1600
1800
1000
2000
1200
Control variable
Time [sec]
55
90
54
80
53 PID
51
50
50
40
49
30
48
20
47
10
46 400
500
600
1800
2000
Controlled variable PID
52
60
0 300
1600
Time [sec]
100
70
1400
700
800
900
45 300
Set-point
400
500
600
700
800
900
Time [sec]
Time [sec]
Fig. 7. Experimental results: (a) and (b) set-point and load step ( 20% change of the controller output at t¼ 1500 s) responses of the plant with the PIplant and PIDplant controllers; (c) and (d) disturbance responses of the plant with the PIDplant controller, for disturbance induced by activating/deactivating the fan. Control variable in [%], controlled variable in [1C].
Table B1 Stable processes (Ms ¼ 2, mn ¼ 2), kn [19]. kn (φ1, ρ) φ¼ 10 φ¼ 20 φ¼ 30 φ¼ 40 φ¼ 50 φ¼ 60 φ¼ 70 φ¼ 80 φ¼ 90 φ¼ 100 φ¼ 110 φ¼ 115 φ¼ 120 φ¼ 125 φ¼ 130
ρ ¼ 0.5
0.4977 0.5084 0.5198 0.5324 0.5393 0.5467 0.5547 0.5634
ρ ¼ 0.55
0.4968 0.5060 0.5159 0.5269 0.5330 0.5396 0.5468
ρ¼ 0.6
0.4891 0.4964 0.5043 0.5129 0.5227 0.5282 0.5341
ρ ¼0.65
0.4848 0.4903 0.4963 0.5031 0.5106 0.5194 0.5243 0.5297
ρ ¼0.7
0.4871 0.4914 0.4964 0.5022 0.5089 0.5167 0.5212
ρ¼ 0.75
0.4871 0.4893 0.4925 0.4967 0.5016 0.5075 0.5145 0.5186
ρ¼ 0.8
0.4910 0.4902 0.4913 0.4936 0.4969 0.5012 0.5064 0.5128
ρ ¼0.85
0.5000 0.4948 0.4929 0.4931 0.4996 0.4973 0.5008 0.5054 0.5113
ρ ¼ 0.9
ρ ¼ 0.95
0.5165 0.5046 0.4982 0.4954 0.4947 0.4956 0.4976 0.5006 0.5047
0.5446 0.5218 0.5086 0.5013 0.4976 0.4962 0.4965 0.4979 0.5004 0.5041
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Estimates in Table 3 are used to design the PIDplant controllers with parameters in Table 4. Parameters are calculated as presented in Appendix B, by using data from the three points in the look-up Tables B1–B3, determined previously for stable processes in [19]. These three points are defined by {ρ2,1 ¼0.90, φ1,1 ¼ 301}, {ρ2,1 ¼ 0.90, φ2,1 ¼401} and {ρ2,2 ¼0.95, φ2,2 ¼ 401}, since ρest and φest, presented in Table 3, correctly classify the laboratory plant in this stable region of the ρ φ plane. Points used to obtain the gains kni,j, kini,j, kdni,j, i,j ¼1,2, in Appendix B, are bolded in Tables B1–B3. The impact of reducing the number of points J on estimates kuest, ωuest, φest, Gpest(0) and on the PID controller tuning is illustrated in Tables 3 and 4. In Table 4, maximum sensitivity Ms is defined by M s ¼ max j1=ð1 þ CðsÞGm ðsÞÞjs ¼ iω , where Gm(s) is the ω
model (1) and (2), defined by parameters in Table 3, and C(s) is the PIDplant controller (4), defined by parameters in Table 4. In both cases, a) and b) in Table 4, the obtained maximum sensitivity is exactly equal to Ms ¼2, used to calculate Tables B1–B3 in [19]. Also, the obtained sensitivity to the high-frequency measurement noise Mn ¼|kd|/Tf [15] is exactly equal to the value mn|ku|, for mn ¼2 used to calculate Tables B1–B3 in [19]. Besides, in Table 4, insignificantly different values of the sensitivity to the high-frequency measurement noise Mn are obtained for different values of J. The anti-reset windup implementation of the PID controller (4) is applied, given by U C ðsÞ ¼ T aw
bks þ ki kd s2 þ ks þ ki 1 RðsÞ YðsÞ þ UðsÞ; T aw s þ 1 T aw s þ 1 ðT aw s þ1ÞðT f s þ 1Þ ð12Þ
7
as in [20], where the saturation element is defined by the input uC(t) and output u(t): 8 uC ðtÞ rllow l ; > < low u ðtÞ; l C ð13Þ uðtÞ ¼ low ouC ðtÞ o lhigh : > :l ; u C ðtÞ Zlhigh high In the linear region of the saturation element, for uC(t) u(t) one obtains (4) from (12). Closed-loop responses presented in Fig. 7a–d, were obtained with the PIDplant controller (12) and (13), having the gains k, ki, kd and noise filter time constant Tf presented in Table 4 for smaller value of J ¼20, and for b ¼0.4, llow ¼0 and lhigh ¼100%. Parameter Taw ¼37 s is obtained from Taw ¼ aTi þ(1 a)Td, for a ¼0.2 satisfying the condition Td oTaw oTi, for Ti ¼ k/ki and Td ¼kd/k. Advantages of the PIDplant controller, compared to the PIplant controller used to determine kuest, ωuest, φest and Gpest(0) in Table 3, are demonstrated in Fig. 7a–b. This experiment starts from the temperature T(L,t) E45 1C, for to1000 s. Then, at t ¼1000 s the set-point is changed for r0 ¼5 1C. At t¼ 1500 s a load disturbance is inserted as a 20% step change of the controller output. As presented in Fig. 7b, improvement of the performance is obtained by the PIDplant controller. Closed-loop experiment in Fig. 7c–d starts from the temperature T(L,t)E50 1C for to300 s, by activating a fan at t¼300 s. Then, at t¼ 600 s the fan is switched-off. Action of the fan induced an intensive cooling of the plate, resulting into a strong disturbance as seen from the control signal u(t) in Fig. 7c. Anti-reset windup action is activated two times, around 320 s and 620 s.
Table B2 Stable processes (Ms ¼ 2, mn ¼ 2), kin [19]. kin (φ1, ρ) φ ¼10 φ ¼20 φ ¼30 φ ¼40 φ ¼50 φ ¼60 φ ¼70 φ ¼80 φ ¼90 φ ¼100 φ ¼110 φ ¼115 φ ¼120 φ ¼125 φ ¼130
ρ¼ 0.5
0.3393 0.3250 0.3119 0.3004 0.2952 0.2905 0.2861 0.2822
ρ ¼0.55
0.2931 0.2803 0.2688 0.2587 0.2542 0.2500 0.2462
ρ¼ 0.6
ρ ¼ 0.65
ρ ¼ 0.7
ρ¼ 0.75
0.2666 0.2543 0.2430 0.2328 0.2239 0.2199 0.2163
0.2441 0.2324 0.2214 0.2113 0.2023 0.1945 0.1910 0.1879
0.2137 0.2032 0.1933 0.1842 0.1763 0.1694 0.1664
0.1976 0.1875 0.1778 0.1689 0.1609 0.1538 0.1478 0.1451
ρ ¼0.6
ρ ¼ 0.65
ρ ¼0.7
ρ ¼ 0.75
ρ ¼ 0.8
ρ ¼ 0.85
0.1835 0.1738 0.1645 0.1557 0.1477 0.1406 0.1343 0.1290
0.1716 0.1620 0.1530 0.1443 0.1364 0.1292 0.1228 0.1173 0.1126
ρ ¼0.8
ρ¼ 0.85
ρ ¼0.9
ρ ¼ 0.95
0.1622 0.1522 0.1430 0.1345 0.1266 0.1194 0.1129 0.1072 0.1023
0.1558 0.1446 0.1349 0.1261 0.1181 0.1108 0.1043 0.0985 0.0935 0.0892
ρ ¼ 0.9
ρ ¼ 0.95
0.8286 0.6798 0.5566 0.4515 0.3591 0.2758 0.1987 0.1257 0.0552
0.9886 0.8040 0.6569 0.5332 0.4314 0.3403 0.2582 0.1823 0.1104 0.0409
Table B3 Stable processes (Ms ¼ 2, mn ¼ 2), kdn [19]. kdn (φ1, ρ) φ ¼10 φ ¼20 φ ¼30 φ ¼40 φ ¼50 φ ¼60 φ ¼70 φ ¼80 φ ¼90 φ ¼100 φ ¼110 φ ¼115 φ ¼120 φ ¼125 φ ¼130
ρ ¼0.5
0.4412 0.3556 0.2722 0.1899 0.1490 0.1082 0.0675 0.0270
ρ ¼ 0.55
0.3934 0.3099 0.2287 0.1488 0.1090 0.0694 0.0299
0.4383 0.3527 0.2712 0.1920 0.1141 0.0754 0.0368
0.4919 0.4018 0.3178 0.2380 0.1607 0.0846 0.0468 0.0091
0.4586 0.3700 0.2876 0.2093 0.1336 0.0591 0.0221
0.5256 0.4293 0.3420 0.2611 0.1843 0.1100 0.0370 0.0008
0.6061 0.4982 0.4032 0.3173 0.2377 0.1623 0.0894 0.0178
0.7047 0.5801 0.4736 0.3800 0.2954 0.2171 0.1429 0.0712 0.0009
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
M.R. Mataušek, T.B. Šekara / ISA Transactions ∎ (∎∎∎∎) ∎∎∎–∎∎∎
8
The analysis of the same laboratory thermal plant in [20], and PID controller tuning, was based on the 5th-order model, obtained by the open-loop 100th-order ARX identification, in the time interval Texp ¼1000 s, and model reduction. The estimates kuest, ωuest, φest and Gpest(0) in (Table 3, row a), are obtained within the shorter time interval Texp ¼410 s in Fig. 2, defined by (6) as Texp ¼(TJþ T/2). Finally, it is important to observe that the estimates kuest, ωuest, φest and Gpest(0) are obtained for Texp E3Tu, since the ultimate period Tu ¼2π/ωuest is equal to Tu ¼141 s for ωuest in (Table 3, row a). 5. Conclusions The proposed procedure for determining the model Gm(s) of stable processes, processes with oscillatory dynamics, integrating and unstable processes, including dead-time, by estimating the quadruplet {ku, ωu, φ, Gp(0)} has been performed in an extremely short time interval. The proposed fast process dynamics characterization, performed without breaking the loop containing a controller in operation, offers the following: (i) first, as demonstrated by experimental results, the controller tuning for the desired performance/robustness tradeoff is obtained by using the estimated kuest, ωuest, φest and Gpest(0), and previously calculated and memorized process-independent look-up tables of parameters kn, kin and kdn, as demonstrated in Fig. 7a–d; (ii) a reliable simulation performance analysis of the loop considered is possible, as demonstrated in Fig. 4. Acknowledgements Authors gratefully acknowledge Dr. Ribić his help in implementing anti-windup PID controller on the laboratory plant. T.B. Šekara gratefully acknowledges the financial support from the Ministry of Education, Science and Technological Development of Serbia (Project TR33020). Appendix A. Determination of φest and Aest by the PLL and SheMa estimators When applying PLL and SheMa estimators, kuest and ωuest are determined by the first experiment. As demonstrated in [5], two additional experiments must be performed, to obtain for ϕ 7 ¼ 7 π=N 0 ; N 0 36: ωu7 , defined by argðGp ðiωu7 ÞÞ ¼ π þ ϕ 7 , 7 7 ρ 7 ¼ ku =ku and ku ¼ jGp ðiωu7 j 1 . Then, φest and Aest are determined from relations [14,15]: φþ þφ ρ 7 sin ϕ 7 A þ þA0 ; Aest ¼ 0 ; φ 7 ¼ arctan 7 ; φest ¼ 7 2 2 ρ cos ϕ 1 2j ωu ωu7 j A07 ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : 1 2ρ 7 cos ϕ 7 þðρ 7 Þ2
ðA1Þ
Appendix B. PID controller tuning by applying look-up tables in the ρ φ plane The process independent normalized gains kn ¼kn(ρ,φ), kin ¼kin (ρ,φ), kdn ¼kdn(ρ,φ) of a virtual PIDn controller, calculated in [19, Tables A1–A3], are presented here as look-up Tables B1–B3. The gains k, ki, kd and the noise filter time constant Tf, in the real
PID controller, are defined by [19] k ¼ kuest kn ; ki ¼ kuest ωuest kin ; kd ¼ kuest kdn =ωuest ; T f ¼ T fn =ωuest ;
ðB1Þ
where the normalized noise filter time constant is defined as Tfn ¼| kdn(ρ,φ)|/mn. From (B1) it follows that the real sensitivity to the high-frequency measurement noise Mn ¼|kd|/Tf, used in [15], is given by Mn ¼ mn|ku|. Results in Tables B1–B3 were obtained in [19] by optimization under constraints on the maximum sensitivity Ms ¼2 and for the normalized sensitivity to measurement noise mn ¼2. However, optimization can be performed under desired constraints on the Ms, mn and for the desired closed-loop system damping ratio ζ ¼1, to obtain the critically damped set-point and load disturbance step responses [15], as done in [19] for processes with oscillatory dynamics. Parameters kn, kin and kdn of the virtual PIDn controller are obtained by using the gains kni,j, kini,j, kdni,j, i,j ¼1,2, from the lookup Tables B1–B3, and by applying interpolation from the three points [5], as defined by kn ¼ ð1 α βÞkn2;1 þ αkn2;2 þ βkn1;1 , kin ¼ ð1 α βÞkin2;1 þ αkin2;2 þ βkin1;1 , kdn ¼ ð1 α βÞkdn2;1 þ αkdn2;2 þ βkdn1;1 , α¼(ρest–ρ2,1)/(ρ2,2 ρ2,1) and β¼(φ2,1 φest)/(φ2,1 φ1,1). Angles in Tables B1–B3 are in degrees and φ1,1 rφ2,1, φ1,1 rφest rφ2,1, ρ1,1 rρest r ρ2,1. For details see [5]. References [1] Hjalmarsson H. From experiment design to closed-loop control. Automatica 2005;41:393–438. [2] Álvarez JD, Guzmán JL, Rivera DE, Berenguel M, Dormido S. Perspectives on control relevant identification through the use of interactive tools. Control Eng Pract 2013;21:171–83. [3] Wahlberg B. Model reductions of high-order estimated models: the asymptotic ML approach. Int. J Control 1989;49:169–92. [4] Ziegler JG, Nichols NB. Optimum settings for automatic controllers. Trans ASME 1942;64:759–68. [5] Šekara TB, Mataušek MR. PID controller tuning based on the classification of stable, integrating and unstable processes in a parameter plane. In: Serra. GLO, editor. Frontiers in advanced control systems. InTechOpen; 2012. p. 117–42. [6] Cohen GH, Coon GA. Theoretical consideration of retarded control. Trans ASME 1953;75:827–34. [7] Mataušek MR, Ribić AI. Design and robust tuning of control scheme based on the PD controller plus Disturbance Observer and low-order integrating firstorder plus dead-time model. ISA Trans 2009;48:410–6. [8] O’Dwyer A. Handbook of PI and PID controller tuning rules. 3rd ed.Imperial College Press; 2009. [9] Lo WL, Rad AB, Li CK. Self-tuning control of systems with unknown time delay via extended polynomial identification. ISA Trans 2003;42:259–72. [10] Mei H, Li S. Decentralized identification for multivariable integrating processes with time delays from closed-loop step tests. ISA Trans 2007;46:189–98. [11] Hashimoto T, Ishida Y. An adaptive I-PD controller based on frequency domain system identification. ISA Trans 2000;39:71–8. [12] Åström KJ, Panagopoulos H, Hägglund T. Design of PI controllers based on non-convex optimization. Automatica 1998;34:585–601. [13] Panagopoulos H, Åström KJ, Hägglund T. Design of PID controllers based on constrained optimization. IEE Proc. Control Theory Appl 2002;34:585–601. [14] Šekara TB, Mataušek MR. Revisiting the Ziegler–Nichols process dynamics characterization. J Process Control 2010;20:360–3. [15] Mataušek MR, Šekara TB. PID controller frequency-domain tuning for stable, integrating and unstable processes, including dead-time. J Process Control 2011;21:17–27. [16] Åström KJ, Hägglund T. Automatic tuning of simple regulators with specifications on phase and amplitude margins. Automatica 1984;20:645–51. [17] Åström KJ, Hägglund T. PID controllers: theory, design and tuning. 2nd ed.. Research Triangle Park, NC: ISA; 1995; 27709. [18] Šekara TB, Mataušek MR. Relay-based critical point estimation of a process with the PID controller in the loop. Automatica 2011;47:1084–8. [19] Šekara TB, Mataušek MR. Classification of dynamic processes and PID controller tuning in a parameter plane. J Process Control 2011;21:620–6. [20] Ribić AI, Mataušek MR. A dead-time compensating PID controller structure and robust tuning. J Process Control 2012;22:1340–9. [21] S Yamamoto, I Hashimoto. Present status and future needs: the view from Japanese industry. In: Proceedings of the IV conference on chemical process control, Texas, USA, 1991.
Please cite this article as: Mataušek MR, Šekara TB. A fast closed-loop process dynamics characterization. ISA Transactions (2014), http: //dx.doi.org/10.1016/j.isatra.2013.12.006i
