- Email: [email protected]
Rapid estimation of PID minimum variance
Accepted Manuscript Rapid estimation of PID minimum variance Farzam Shahni, Wei Yu, Brent Young
PII: DOI: Reference:
S0019-0578(18)30426-9 https://doi.org/10.1016/j.isatra.2018.10.047 ISATRA 2948
To appear in:
ISA Transactions
Received date : 1 March 2018 Revised date : 10 August 2018 Accepted date : 31 October 2018 Please cite this article as: Shahni F., Yu W. and Young B. Rapid estimation of PID minimum variance. ISA Transactions (2018), https://doi.org/10.1016/j.isatra.2018.10.047 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
Rapid estimation of PID minimum variance Farzam Shahni, Wei Yu, Brent Young
Chemical and Materials Engineering Department, The University of Auckland, New Zealand.
Corresponding author: Brent R Young Professor of Food & Process Systems Engineering Head of Chemical & Materials Engineering Chair of the Food & Health Programme University of Auckland T +64 9 923 5606 Email: [email protected] A Room G61A, 2-6 Park Avenue Grafton, Auckland, NZ 1023
*Highlights (for review)
The highlights: 1. Assessment of PID achievable performance for linear, stationary SISO processes. 2. Guidelines developed for choosing the appropriate size of the finite impulse response for the system. 3. Fast non-iterative method to estimate the PID benchmark.
*Blinded Manuscript - without Author Details Click here to view linked References
1
2
Abstract
3
PID controllers are the most common type of controllers used in industry. However, to find
4
the best PID performance is a challenging task. Assessment of the best PID performance
5
leads to the minimization of the deviation from set point as indicated by closed loop transfer
6
function. This is a non-convex optimisation problem and there is no direct solution for
7
finding this benchmark. Solving this problem iteratively as proposed in literature requires
8
long calculation times. Furthermore, iterative methods are not guaranteed to find the
9
minimum value. Using as many as possible infinite impulse response coefficients, or
10
increasing the number of finite impulse coefficients gradually, causes undesirably large
11
calculation times. We propose a fast method to evaluate the minimum variance by a fixed
12
number of finite impulse response coefficients. This finite impulse model is double the size
13
of the stable system impulse response (it contains the transient response and the same
14
duration of steady state). Thus it avoids undesired iteration, and consequently, the PID
15
minimum variance is quickly evaluated. Time is a critical factor as the purpose to assess this
16
index is for use in online monitoring. The proposed method is tested on benchmark
17
simulation
examples
from
1
literature.
18
1
Introduction
19
Surveys show that almost 60% of industrial control loops are working with poor performance
20
(60% in general, could be due to controller problem, final control element problems or
21
sensors problem, linear or non-linear, feedback or feedforward) and only 40% of processes
22
are working under desired optimal performance [1, 2]. Assessment of control loop
23
performance has been a topic of interest for many researchers [3-10]. Performance
24
assessment is mainly used to verify the working performance of one controller with
25
consideration of input and output variance. This assessment compares the current controller
26
performance and the best achievable performance. In the assessment of control loop
27
performance, factors such as time delay, disturbance and process models, and structural
28
constraints of the controller affect the system output performance. For a linear, discrete-time
29
system with additive disturbances, the lowest output variance is achieved when a minimum
30
variance (MV) controller is implemented. This MV controller becomes the standard
31
benchmark to assess control loop performance [11]. Even though using MV control as a
32
performance benchmark can provide reasonable performance, this controller is not practical
33
due to its demand for excessive control effort and poor robustness [6, 12]. Also, MV is not
34
the only and a complete objective in the reality of process control.
35 36
There is no doubt that more than 90% of the applied controllers in industry are implemented
37
with PID algorithms [7, 13-15]. Some even believe 98% of controllers in industry are PID
38
type [16, 17]. It is reported that 80% of PID controllers in industry are poorly tuned [16, 17]
39
(major reason of the 60% poor performance in general loops is due to the lack of the
40
resonable PID controllers tuning).
41
From the above-noted reports, it can be concluded that more than 70% of all controller bad
42
performance in industry is due to bad tuning of PID parameters. When the desired controller 2
43
is a PID type, performance indices based on MV controllers as presented in, e.g. [18, 19] are
44
not suitable to tell whether the PID type controller can provide good controller performance.
45
PID controllers have both order and structural constraints which result in the possibility of
46
not achieving MV performance, and thus the best performance of the practical controller is
47
unknown. As Kozub reported in [20], only 20% of loops in a refinery with PID controller
48
structure may achieve the MV control performance lower bound. Eriksson and Isaksson in
49
[21], proposed an index with an optimal PID controller as a benchmark to replace MV.
50
Hence, it is important to assess the minimum variance with PID control (MVPID). Although
51
the benchmark controller can provide reasonable performance, benchmark controllers are not
52
practical due to their demand for excessive control effort and poor robustness and in practice
53
a less aggressive controller will be selected.
54
The problem of assessing the MVPID produces a non-convex problem, and a globally optimal
55
solution cannot be guaranteed [22]. There is no direct and clear method to solve this complex
56
non-convex problem. Researchers have tried to solve the non-convex problem by various
57
convex methods. A gradient-based method was proposed by Ko and Edgar [23, 24]. An
58
analytical lower bound method was proposed by Kariwala [25], using the sum of the squares
59
approach proposed by Sendjaja and Kariwala [26]. Other similar approaches have also been
60
proposed [27-29].
61
There are two main drawbacks with the current methods of solving the MVPID assessment
62
problem: The first is calculation time (avoiding iteration), and the second is an accurate
63
estimation. Time is a critical factor as the purpose to evaluate this index is for use in online
64
monitoring. Most of the current literature uses iteration to solve this problem which causes
65
the calculation time to rise. Even if the iteration scheme is designed to increase the chance of
66
finding the minimum value, it may provide the minimum value or it may not. In these studies,
67
there were not any theoretical methods identified for finding the suitable finite length. A large
3
68
pre-defined finite size [27, 30] or increasing the size of finite length (impulse response
69
coefficients) from a small (e.g., 2d-1, where d is delay) to a large value (e.g., 4d, the size of
70
iteration is 2d here) was proposed by [25, 28]. These methods find the minimum from the
71
finite set and from the results (controller values) it calculates the variance of the closed-loop
72
system in each step. The results are saved and compared with the next sized data set. This
73
procedure causes calculation time to increase unnecessarily. The unclear, random selection
74
and insufficiently truncated the finite impulse response is the reason why these methods miss
75
the optimum time and optimum value.
76
This paper proposes a fast method to estimate MVPID which is capable of being used in online
77
monitoring. This technique reduces calculation times by using a fixed finite length and
78
removing the iteration. It produces reliable results as it is compared with the existing
79
methods. Due to the fast computation time, the method could be implemented for online
80
performance assessment. A full theoretical analysis is presented, and a simulation study is
81
performed for the assessments using the Matlab optimization toolbox. In addition to the fast
82
method proposed, we introduce a weighting method for more accurate results and as an
83
alternative solution for these nonconvex problems. This weighting method changes the
84
weight of the transient and the steady state impulse response and balances these two parts in
85
the cost function.
86
It is important to note that this is a simulation study, and constraints like manipulated variable
87
limitations are not considered. This is common practise in such studies and not usually a
88
major concern for process control as the plant dynamics are usually much slower than the
89
actuator dynamics. At this stage use of this method is also limited due to its demand for
90
information regarding the process and disturbance model. Thus it is not a feasible method for
91
loops where appropriate models do not exist or can be obtained and we would recommend to
92
use the proposed method only on important loops where models exist or are able to be
4
93
identified. The scope of this work is for linear plants and models, or any plants that have mild
94
nonlinearity which is well handled by PI and PID feedback (the majority of industrial process
95
plants).
96
This paper is organized as follows: Section 2 introduces the problem of assessing the
97
minimum variance of a PID controller. Section 3 consists of the analytical formulation of the
98
PID minimum variance calculation. Section 4 consists of the proposed methods. Section 5
99
and its subsections involve simulation study, examples and comparison. In this section,
100
fifteen examples have been used to verify the usefulness of the proposed method. Finally, this
101
paper is then summarized along with some concluding remarks in the last section.
102
2
103
A block diagram of a typical single-input-single-output (SISO) feedback control system is
104
shown in Figure 1 with output y (t ) , manipulated variable u (t ) , and zero mean white noise
105
a (t ) , where t is the sampling interval of the one-time unit. The process output can be written
106
as:
Best Achievable PID performance
y(t ) G(q 1 )u(t ) H (q 1 )a(t )
107
(1)
G(q 1 ) and H (q 1 ) represent the process and disturbance transfer functions,
108
where
109
respectively. For simplicity and brevity, the backward shift term q , q
110
omitted in the subsequent sections unless circumstances necessitate its presence.
1
1
y(t ) y(t 1) , is
111 112
When there is no set point change, the output y (t ) under closed loop conditions (with set point
113
zero), can be expressed by:
5
y (t )
114
H a(t ) Gcl a(t ) 1 G Gc
(2)
115
116
Gc is digital PID controller with the structure of: Gc
117
k1 k2 q 1 k3q 2 1 q 1
(3)
k1 = kP +kI +kD , k2 (kP 2kD ) , k 3 k D . k P , kI and k D denote proportional, integral
118
Where
119
and derivative gain, respectively. Gcl is defined as the closed loop transfer function between
120
y (t ) and a (t ) . The output variance under the PID controller is:
Var ( y ) Gcl 2 . a2 2
121
(4)
122
a2 represents the variance of the random noise
123
all its poles and zeroes inside the unit circle thus G contains no zeroes and poles on or
124
outside the unit circle except at infinity (time delay) and H contains no zeroes on or outside
125
the unit circle [31]. In this work, the variance of the white noise sequence ( a ) is fixed to be
126
unity.
127
The system output in Eq.(1) could be written in its impulse response form as follows:
a (t ) . For solving Eq. (4) we suppose Gcl has
2
i 1
i 1
y (t ) gi u (t i ) hi a (t i )
128
(5)
gi and hi are the impulse response coefficients of the process transfer function and
129
Where,
130
disturbance model, respectively. By substituting for the PID controller given in Eq. (3), the
131
output from the Eq. (5) becomes:
6
i 1
i 1
y (t ) gi Gc . y (t i ) hi a (t i )
132
133
By considering the definition of the backward shift operator ( q
134
rewrite Eq. (6) as follows:
1
(6)
y(t ) y(t 1) ), we can
y (t )
135
h a(t i) i 0
i
1 2 i 1 (k1 k2 q k3q ) si (q ) i 1
136
Where si represents the step response coefficient of the process transfer function. We can
137
write the impulse response of disturbance
139
1
0 D 1 0
1
0 0
(9)
The denominator of Eq. (7) can then be written in matrix form as follows:
S I (k1 k2 D 1 k3 D 2 )( s1D 1 s2 D 2 s3 D 3 ...)
142 143
(8)
In matrix form the backward shift operator ( q ) can be written as:
140
141
Hi (h0 , h1 , h2 , h3 ,...) in matrix form as follows:
h0 h h0 h1q 1 h2 q 2 h3q 3 ... 1 h2
138
(7)
(10)
Where I is the identity matrix. We then rewrite Eq. (7) as follows:
y (t ) i a (t i )
144
i 0
145
and
7
(11)
0 h 1 2 S
146
(12)
147
Where is the impulse response vector of the closed loop transfer function. The i ’s are
148
h0 h the closed loop impulse response coefficients of the SISO system and h 1 is the impulse h2
149
response vector of the disturbance transfer function. The output variance is defined as:
y2 i 2 a2
150
(13)
i 0
151
The minimum output variance with an unrestricted controller structure is known as the MV
152
benchmark [11] and is calculated as: d 1
2 MV hi2 a2
153
(14)
i 0
154
where d is the process time delay. If the controller structure must be PID, the achievable
155
minimum variance low bound MVPID is not equal to MV and cannot easily be found. To
156
minimize the output variance with a PID controller structure, one needs to find the
157
k1 , k2 and k3 values in Eq. (15), so that y2 reaches the minimum value. This leads to an
158
optimization problem which can be expressed as:
159
2
min y2 min
k1 , k2 , k3
k1 , k2 , k3
i 1
8
2 i
a2 min Gcl 2 a2 2
k1 , k2 , k3
(15)
160
Where . 2 denotes the 2 norm . Eq. (15) yields a non-convex optimization problem. There is
161
no direct and simple solution for such an equation [22]. Be noted that here it is considered
162
that a set point is constant in value and purpose is to minimize the effect of noise.
163
3
164
The H 2 norm of a stable transfer matrix S is:
Problem formulation 2
165
Gcl
2 2
( T )
(16)
166
From Eq.(2), due to the controller parameters being unknown as stated in Eq.(3), is
167
unidentified. In general is obtained from a finite impulse response [32]. However, it is
168
essential that the finite impulse response contains all non-zero coefficients. Here, to obtain
169
up to the m th element, the finite samples of the closed loop transfer function Gc l are
170
expressed by the following,
171
0 m 1 Sm 1 h m m1
(17)
172
Where hm is defined as the impulse response of the disturbance model H up to the m th
173
element. The term hm can be represented as,
174 175
hm (h0 , h1 , h2 ,..., hm1 )1Tm From Eq.(10) and (17), the truncated matrix S, is defined as:
9
(18)
1 s 1 S m s2 sm 1
176
177
If
m
is taken to be very large then
178
179
0
0 0
s1
0
sm 2
sm 3
0 0 0 1 mm
(19)
m and then from Eq.(16) and Eq.(17) we have: Gcl
2 2
( T ) mT m
(20)
In other words, we have the equation below:
0
180
mT m 1 ( T )
(21)
181
When sufficient impulse length is selected, the Eq. (21) becomes equal to unity, and if m is
182
very short Eq. (21) becomes zero. Directly solving Eq.(15) is impossible, instead Eq. (20) is
183
used to find the approximate minimum output variance for PID controllers [23-29, 32]: min y2 min mT m a2
184
PID
m
PID
(22)
is limited to the range 4d m 8d due to computational problems (the
185
In [28] the size of
186
iteration is 4d). Another way is to increase
187
to reach the convergence of the left-hand side (upper bound) and right-hand side (lower
188
bound) in Eq. (20) [26]. In this case the size of iteration could be anything between 2 to
189
infinity, however, in the program one should consider a limitation due to lack of computer
190
memories like 30 or 40. A pre-specified constant,
191
these problems. The authors here intentionally avoid iteration and instead increase the size of
192
the matrix which again causes the calculation time to increase.
m
starting from the analytical lower bound 2d -1
10
m (m=500) is also used in [27] to address
193
4 The proposed method
194
is separated into two parts, m , thus Eq. (21) can be rewritten as: r
mT m mT m 1 T T T m m m m r r r r
195
(23)
196
From Eq. (23) it can be concluded r r 0 . This term has led previous research to neglect
197
the system impulse response after
198
found by solving the optimisation problem in Eq.(22). For validation of the solution the
199
controller parameters which are found through Eq. (22) are used in Eq. (16). At this stage, if
200
the size of
201
Then the value found in the previous step is saved and is used to compare with the result of
202
the next iteration. In this case, some previous research gradually increased the matrix size, m
203
in Eq.(17). Gradually increase the size of
204
calculation time to increase and most results cause Eq. (16) (the upper bound), to become
205
infinitely large or far from the minimum amount. Some tried to choose a large
206
output variance as reported in [23-29, 32]. However, a pre-specified constant, large
207
similarly can lead to an unnecessary rise in calculation time.
208
4.1 Fast method
209
By increasing the size of the finite horizon the cost function Eq. (16) increases. However
210
from some point (appropriate m ) this increase in the cost function should stop (when it hits
211
the lower bound or the exact optimum) or becomes negligible (finds a local optimum). Now
212
we define the denominator of Eq. (17) as:
213
T
m is
m samples.
In the literature, the minimum variance is
not appropriate, the result will be far from the actual minimum variance.
m from a small value via trial and error, causes the
A A1 S m C 11
0 A mm
m to calculate
m
(24)
A1 consists of step response coefficients.
214
As is shown by Eq. (10) and Eq.(19) block matrix
215
Here m is selected to be whenever the step response of a process is within 99% of the final
216
value (to see the reason behind the selection of 99%, Section 5.1). The reason behind the
217
above selection of m is to have approximately constant coefficients in C. It is mentioned
218
before that the 2-norm is equal to the root-mean-square of the impulse response of the
219
system, and the variance is the squared of the 2-norm. Thus the first
220
contain the final minimum variance and the second
221
optimum variance. We consider the block matrix A as the transient part of the step response
222
and C as it’s approximately steady state part (the elements of being defined as the steady
223
state step response coefficients multiplied by the controller’s gains). The following
224
illustration Figure 2 is an example of the system model, G
225
of the matrix, m, is selected in this paper.
2
m impulse coefficients 2
m impulse coefficients contain the 2
q 5 , to clarify how the size q 0.8
226 227
The inverse of the lower triangular matrix
A1 A11 ( Sm ) 1 1 1 A CA
228
229
230
A1 is given below:
The matrices Aand
0 A1 mm
(25)
A1 are lower triangular matrices. Now matrix m is calculated as:
A1 m 1 A11 h m A1CA1 2
12
A1 h m2 0 nm A1 (hc CA1 h m ) A1 2 mm
(26)
231
h m Where h m2 and h c are defined as h m 2 , the impulse response samples of the disturbance hc
232
hm 2 h m 1 before and after term h m2 , respectively ( hc 2 ). If we increase the size of the matrix in h m1
233
Eq. (26) to 2m , 2m becomes:
2 m m S 2 m 1 h 2 m 2
234
235
1
h m2 hc hc hc
(27)
C C . Then 2m becomes: C C mm
2m
h m 1 A11nm 2 2 0 hc m h (28) 1 1 1 h c A11 hc A1 C1 A1 2 A1 3 hc hc 4 h c
A11 A 1C A 1 1 1 1
Then we have 3 , 4 as: m hc 3 1 1 h 2 A1 C1 A1 4 hc hc
238 239
0 0 0 A
We define the matrix C1
236
237
A 0 0 C A 0 C C A C C C
(29)
where, 3 A1 (CA1 I )(hc CA1 h ) m 2
240
and
4 A1CA1 (CA1 I )(h c CA1 h ) A1 (CA1 I )(nc CA1 h ) m 2
m 2
A1 (CA1 I ) 2 (h c CA1 h m2 )
13
(30)
241
We can continue this calculation and see the effect of the term (CA1 I )(nc CA1nm ) on future 2
242
prediction. 4m could be written as follows:
4m
243
244
A1 0 0 C A1 0 1 C1 C1 A1 C1 C1 C1
A1 h m2 h m2 1 1 1 A (hc CA h m2 ) hc 1 2 1 1 1 m ) A ( CA I )( h CA h c 3 2 0 hc A1 (CA1 I ) 2 (h CA1 h m ) 0 hc c 2 4 0 hc A1 (CA1 I )3 (h c CA1 h m ) 5 2 A1 hc A1 (CA1 I ) 4 (h CA1 h m ) 6 c 2 hc A1 (CA1 I )5 (hc CA1 h m ) 7 2 hc 8 A1 (CA1 I )6 (h c CA1 h m2 )
(31)
Eq. (16) can be rewritten as below: 2
245
Gcl ( T ) 1T 1 T2 2 T3 3 T4 4 T5 5 T6 6 T7 7 T8 8 ...
246
(32)
247
Thus the equation below is repeated over time (all rows from Eq. (31) include the term
248
below);
2
249
(33)
(CA1 I )(hc CA1 h m2 )
250
From the above equation (33), when the disturbance is assumed to be a step in the output
251
1 n ( the common case for a disturbance [33]), the impulse response of the closed-loop 1
252
system in each period
m is: 2
1 i A (CA I ) h m2 A (CA I ) 1 1
253
1
i 1
1
1
i 1
(34)
254
The impact of the
255
solution of Eq. (22). Thus for solving Eq. (15), we consider the first two terms of Eq. (32):
(hc CA1 h m2 )
term in Eq. (31) means we should consider the term for the
14
min S 2 a2 min( 1T 1 2T 2 ) a2 2
256
PID
m is
(35)
PID
selected as “large enough” and by chance may cover the
257
In previous research when
258
terms 1 , 2 , 3 ,... , the control parameters may be tuned well and achieve comparable results,
259
but purely by chance. An unnecessary large m causes extra calculation and extra time is
260
required for finding reliable optimum results. A small
261
To find the solution of Eq. (15) by solving Eq. (22), previous research [26-29] only focused
262
on the first
263
consider the important effect of the remaining samples. Thus the selection of impulse length
264
according to the criteria proposed in the first part of this section is appropriate, and it avoids
265
unnecessary iteration. In Figure 3 and Figure 4 it is shown that what will happen if we just
266
consider a short length of the finite impulse response and do not consider its future effects. In
267
these figures, the size of the matrix, m, was not appropriate. Regardless, the prior approach
268
tried to minimize the variance based on length m. It is depicted in these figures that selecting
269
finite length impulse response could minimize the right-hand side of Eq. (22). In these cases
270
in which the truncated length of the impulse response is not appropriate, the Eq. (21) is 0.
271
Example 1: By way of example we considered a system with the process model
272
G
273
impulse response; we compared 1 1 , 2 2 in Eq. (35) and saw its effects on Eq. (23),
274
based on optimum and non-optimum results, Table 1.
275
Table 1 shows how the selection of the wrong sized matrix (m) could affect the minimization
276
problem - the result is far away from the optimum value. In these cases, the short length of
277
the impulse response is considered without knowledge of its future.
m samples,
m
also leads to non-optimal results.
where they used iteration to find a suitable m , and they did not
0.1 0.1q 5 , and disturbance model H . Figure 3 is the 1 1 (1 q )(1 0.3q 1 )(1 0.6q 1 ) 1 0.8q T
T
15
278 279
Example2: In this example (Table 2) we consider a system with the process model
280
G
281
response of the system. Here we want to show where Eq. (35) is valid based on a comparison
282
of 1 1 , 2 2 . We also compare the left-hand side of the equation (the upper bound or
283
actual norm which is based on infinite impulse response) and right-hand side (lower bound
284
which is based on finite impulse response). These two bounds should be merged to validate
285
Eq. (35). Table 2 shows how bad tuning results due to wrong selection of impulse response
286
length from prior methods affects the Eq. (35).
287
As reported in literature for the non-convex optimization problems in general there is no
288
guaranteed method [22]. However, we have proved throughout Eq. (24) to Eq. (35) in Section
289
4 that if we select the correct length of finite horizon (which we proposed in section 4.1) the
290
method can guarantee a feasible solution and optimum value just for this problem. However
291
for more accuracy the weighting method is proposed.
292
4.2 Weighting method
293
When our objective is minimization, where one is concerned with obtaining more accurate
294
results, the impact or weight of the
295
greater than A1 h m2 , the first term in Eq. (22). Therefore, for more accuracy, the weight of the
296
first impulse response is reduced when solving Eq. (15). Hence, the optimization problem is
297
as follows,
298
(0.5285 0.2707q 1 )q 2 1
1
(1 0.375q )(1 0.3606q )
T
, and disturbance model H
1 . Figure 4 is the impulse 1 q 1
T
(hc CA1 h m2 )
term in a solution of Eq. (31), must be
min y2 min S 2 a2 min( 1T 1 2T 2 ) a2 2
PID
PID
PID
16
(36)
299
µ is a weighting coefficient which should be any value less than 1. When µ is small, the result
300
is far away from the minimum value. When µ is close to 1, we can tune to achieve an
301
optimum result. The scale in Figure 5 is enlarged to show how the weighting parameter
302
affects the optimization problem in Eq.(36). However, when µ is close to 1, the optimization
303
result is not too far from the minimum possible result. In previous studies, in some iterations,
304
the controller parameters found through Eq.(22), caused the variance of the system to become
305
infinite (Eq. (21) became zero) due to improper
306
small size and is gradually increased).The difference between Eq. (22) and Eq. (36) could be
307
shown as follows:
1
308
309
The right-hand side of the above equation,
Eq.(15) .
1T 1 2T 2 Gcl
2
(as mostly the size of
m
starts off being of
1T 1 2T 2 mT m 2 2 Gcl Gcl 2 2 mT m Gcl
310
m
2
(37)
, shows Eq. (22) is an under-estimation for
2
shows Eq. (36) is an unbiased estimation. It should be noted that
2
311
small leads to the result of the optimization problem in Eq.(15) to reach a large value for
312
1T 1 and smaller values for 2T 2 . The result should always be a balance between these
313
two parts of Eq. (36) to reach to an optimum value.
314 315
A detailed simulation study has been performed for assessment of the proposed method using
316
the MATLAB package with the optimization toolbox. We used large scale algorithm based
317
on interior-reflective Newton method. Any optimization function could be called to find this
318
local minimum of constrained nonlinear multivariable polynomials like the nonlinear
319
fmincon
function
in
MATLAB
17
with
using
option
settings
320
'TolCon ',1e 6,'TolX ',1e 6,'TolFun ',1e 6 . This optimization function finds a minimum
321
of a constrained nonlinear multivariable scalar function, but it is highly dependent on initial
322
values. To ensure for positive values of the controller parameters, inequalities k P
323
and
324
included in the analysis. For solving the optimization problem, the first step is to start the
325
optimization with an initial value of x0
326
order
327
calculation time drops from minutes to a few seconds.
328
This part is added as an alternative solution when the size of the finite horizon is not selected
329
properly. For better explanation please see Figures 6 and 7 below.
330 331
5 Simulation examples
332
models. However, if the system model is not known, system identification methods may be
333
used to obtain the open-loop process model, and by using output samples based on the
334
installed controller, the closed-loop disturbance model can be identified [7, 23]. We justify
335
this approach as follows. Even though no exact model can ever be obtained, for most
336
chemical processes system identification can find useful and appropriate models that are
337
sufficiently accurate and capture the essential dynamics that are necessary for this approach.
338
The method has been tested for more than 300 different scenarios and processes, for which
339
this method worked properly. From these 300 cases, 15 representative examples were
340
selected for the purpose of comparison with the previous methods. Computational
341
performance on various tests show the advantages of our approach.
0 , kI 0
kD 0 are considered. Thus inequality constraints such as k1 0 , k2 0 and k3 0 are
m
[0;0;0] . In this research, by selecting a fixed matrix
, which is described in Section 4, and by removing the undesired iteration, the
The proposed method is based on complete knowledge of the process and disturbance
18
342
Table 3 shows examples of process and disturbance models. These first order plus time delay
343
(FOPTD) and second order plus time delay (SOPTD) examples demonstrate the efficiency of
344
the proposed algorithm for general application.
345
The MV, results of other researchers, MVPID and optimal PID parameters obtained by the
346
proposed methods are given in In Table 2. In total 15 case studies were compared, and it was
347
found that our method either equalled or outperformed all the other previous methods [26-
348
30]. The Table 4 shows the weighting method is not case sensitive. It can be concluded that
349
the method is reliable in general and it can be used as a standard solution reference or
350
benchmark. In no cases did the proposed method perform inferiorly to previous studies, albeit
351
in some cases it yielded the same results as other approaches. It can be deduced that the
352
proposed method is more reliable as it always produced better results or performed the same
353
as other studies. The cases in which the method performed like others have been highlighted
354
with green colour in Table 4. The low accuracy results are highlighted with a red colour. As
355
is shown in Table 4 the previous methods do not perform well for higher order systems.
356
As expected, mostly there is a gap between MVPID and MV (as explained before in the
357
introduction section, this gap is due to limitations based on the fixed order and structure of
358
PID controllers). It is precisely indicated how the optimal behaviour of industrial controllers
359
(PI/PIDs) differ from MV. Examples 4-7 in Table 3 have the same process models, but they
360
have different disturbance transfer functions. The results of these examples in Table 2, show
361
the expected gap between MVPID and MV increases when disturbance changes from auto-
362
regressive moving average (ARMA) to auto-regressive integrated moving average (ARIMA).
363
The simulation study has been performed using an Intel Core i7 3.40 GHz CPU, 16 GB Ram
364
PC. Calculation time for the proposed method is two seconds as shown in Table 5 which is
365
much faster than prior methods. The second column of Table 5 is the reported calculation
366
times of the other methods [26-28]. The reported times in Table 5 have been updated by 19
367
repeating the simulation for all 15 cases using the same machine used for the proposed
368
method to have a fair comparison. The last column in this table is the calculation times of
369
second-order systems which show all the previous methods are time-consuming compared to
370
the proposed method. For second-order systems, as is shown there is almost a 600 seconds
371
gap between the proposed method and the previous methods. These updated times are shown
372
in the third column of Table5.
373
Table 5 shows the proposed method is fast and could be used in the online assessment. Based
374
on the comparison in Table 5, the method without weighting (the fast method) shows
375
reasonable results with just less than 0.005% deviation from the best possible results as
376
compared with weighting method (the standard solution) and the current methods proposed in
377
the literature. The simulation results verify the capability of the proposed method which
378
could be implemented and acceptable for online PID control performance assessment.
379
According to the obtained results, this method is the best method to achieve MVPID as shown
380
in Table 4.
381
Because of the simplicity of the method and the speed, the method could be potentially used
382
for general purposes and for higher order systems (Examples 12-15 in Table 3 are second-
383
order systems). Application of the method does not need high-performance computers for
384
complicated systems. It is worth mentioning that this method is also applicable to find PI
385
minimum variance. To find the PI minimum variance, in Eq.(3) we must set k 3 k D 0 , so
386
Eq.(3) becomes K PI
387
may be worthwhile to say here that in the presence of non-stationary disturbances, PID
388
controllers exhibit better performance than PI controllers, due to the prediction capability of a
389
derivative action in PID controllers. Our experience in this work study indicates that our
390
proposed method represents an efficient way to find minimum variance problem by PIDs.
k1 k2 q 1 and the rest of the procedure remains unchanged. However, it 1 q 1
20
391
5.1 Size of finite impulse
392
In Table 6 we show how the selection of matrix size (finite impulse length) affects the optimum value
393
of Eq.(36). As is shown in this table, there is almost no difference in the result of the optimization for
394
the
395
However, if m is selected to be only 50% of the final value, the optimization result is far
396
away from the optimum, and in some cases it becomes infinity. Figure 8 shows how the result
397
will be affected by the size of m, based on the percentage of final value. In Figure 8, the 15
398
different symbols belong to the optimization results of the 15 different case studies. To scale
399
values and compare all 15 examples, we divided all optimization results for each case by
400
numerical results (the matrix size m is approach to infinity). Figure 8 shows why we selected
401
the size of the finite impulse response to be based upon the approximation of the final value
402
as being within 99% to 90% of the steady-state step response in section 4.1. In this figure, 1
403
indicates the optimum value for the all 15 cases.
404
We investigate if there is almost the same result obtained if we select finite impulse response
405
coefficient (m) to be within 99% to 50% of the final value of the steady state step response.
406
The first four points are very similar however, if we take a closer look at the Figure8 we can
407
see the scores at the 99% level are the smallest. The proposed method is tested for more than
408
300 different scenarios and always worked properly.
409
5.2 Effect of delay and disturbance in size of the finite horizon:
410
Here we examine the probable impact of the delay and disturbance of system into the
411
appropriate finite impulse. Some researchers like [24, 26], suggested the start point for
412
optimization is related to delay of the system. They start the length of the finite horizon from
413
the double size of delay (2d-1) optimization. Table 7 shows there is no relation between the
414
proper size of m and the delay of the system. Also, a pre-specified large m (m=500 or
m
selected within 99% to 90% of final value.
21
415
m=1000) as proposed by other researchers [27, 30] could cause unnecessary long calculation
416
time. For example, for case 9, size of the finite impulse response coefficient should be 26
417
while [27, 30] considered the size of m be 500 or 1000. Thus the above mentioned approved
418
we cannot predefine m by chance, and it is case base. In Table 7 we also consider the
419
possibility any relation between the size of m and disturbance dynamic. It is shown there is
420
no relation between the size of the matrix, m, and disturbance dynamics.
22
421
6
Conclusions
422
This paper has proposed a new method to find the performance lower bound for PID type
423
controllers and single-input-single-output (SISO) systems to be used for online control
424
performance assessment. The calculation time for this method dropped to less than two
425
seconds which is capable it to be used for online performance assessment. The method
426
consists two parts. The first part selects the length of impulse response coefficients of the
427
closed loop transfer function between the output of and a disturbance to the system, which
428
causes the method, be fast and simple. The second part is optimization with or without a
429
weighted polynomial function arising from the impulse response. The purpose of the second
430
part is producing precise results. The simplicity of the method enables this method to be used
431
for general problems of higher order. It has been shown from theoretical and simulation
432
studies that the proposed algorithm is accurate and fast. The method is able to be extended to
433
find the achievable PID performance for the general class of all applications. This research
434
covers the most currently used controllers in the industry (PID and PI) and could identify
435
how a current controller is far from an ideal controller with the same performance or action.
436
As shown, the method is fast, so it potentially can be used for MIMO systems. However, it is
437
worthwhile noting here that finding the MV benchmark is not the only step in implementing
438
real process control. The next phase of this work is focused on extending the method to
439
MIMO processes. The most fundamental constraint in process control is the time delay for
440
SISO systems which are time invariant. Since the actual action of time invariance for MIMO
441
systems is done by an interactor matrix, the two main challenges here are large dimensional
442
matrixes which cause numerical errors and computational times that need to be short.
443
23
References
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
[1] [2]
[3] [4]
[5]
[6]
[7] [8] [9] [10] [11] [12] [13]
[14]
[15] [16]
[17] [18] [19] [20]
T. Samad and A. Annaswamy. (2014). The Impact of Control Technology (2nd ed.). Available: www.ieeecss.org T. Harris, C. Seppala, and L. Desborough, "A review of performance monitoring and assessment techniques for univariate and multivariate control systems," Journal of Process Control, vol. 9, no. 1, pp. 1-17, 1999. Y. Alipouri and J. Poshtan, "Minimum variance lower bound estimation and realization for desired structures," ISA Transactions, vol. 53, no. 3, pp. 787-792, 2014/05/01/ 2014. A. Silveira, R. Trentini, A. Coelho, R. Kutzner, and L. Hofmann, "Generalized minimum variance control under long-range prediction horizon setups," ISA Transactions, vol. 62, pp. 325-332, 2016/05/01/ 2016. Z. Yu and J. Wang, "Performance assessment of static lead-lag feedforward controllers for disturbance rejection in PID control loops," ISA Transactions, vol. 64, pp. 67-76, 2016/09/01/ 2016. B. Huang, "A pragmatic approach towards assessment of control loop performance," International Journal of Adaptive Control and Signal Processing, vol. 17, no. 7‐9, pp. 589-608, 2003. M. Jelali, "Control performance management in industrial automation," Advances in Industrial Control, DOI, vol. 10, pp. 978-1, 2013. B. Huang and S. L. Shah, Performance assessment of control loops: theory and applications. Springer Science & Business Media, 1999. M. Jelali, "An overview of control performance assessment technology and industrial applications," Control Engineering Practice, vol. 14, no. 5, pp. 441-466, 2006. S. Huang and C. Huang, "Optimal PID tuning for discrete-time systems with stochastic disturbance," in Control Conference (CCC), 2014 33rd Chinese, 2014, pp. 9002-9007: IEEE. T. J. Harris, "Assessment of control loop performance," The Canadian Journal of Chemical Engineering, vol. 67, no. 5, pp. 856-861, 1989. N. F. Thornhill and B. Huang, "Management and monitoring of process assets," Advanced control of industrial processes ADCONIP, pp. 4-7, 2008. L. Desborough and R. Miller, "Increasing customer value of industrial control performance monitoring-Honeywell's experience," in AIChE symposium series, 2002, pp. 169-189: New York; American Institute of Chemical Engineers; 1998. L. Desborough and R. Miller, "Increasing customer value of industrial control performance monitoring-Honeywell's experience " presented at the Chemical Process Control-VI, Tucson, Arizona, 2001. K. J. Åström and T. Hägglund, Advanced PID Control. ISA - The Instrumentation, Systems and Automation Society, 2006. P. Agrawal and S. Lakshminarayanan, "Tuning proportional-integral-derivative controllers using achievable performance indices," Industrial & engineering chemistry research, vol. 42, no. 22, pp. 5576-5582, 2003. P. Van Overschee and B. De Moor, "RAPID: The End of Heuristic PID Tuning," IFAC Proceedings Volumes, vol. 33, no. 4, pp. 595-600, 2000/04/01/ 2000. L. Desborough and T. Harris, "Performance assessment measures for univariate feedback control," The Canadian Journal of Chemical Engineering, vol. 70, no. 6, pp. 1186-1197, 1992. D. J. Kozub and C. E. Garcia, "Monitoring and diagnosis of automated controllers in the chemical process industries," in AIChE annual meeting, 1993, vol. 9. D. J. Kozub, "Controller performance monitoring and diagnosis: experiences and challenges," in Chemical process control conference, Lake Tahoe, USA, 1996, pp. 83-96.
24
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
[21]
[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]
[34]
[35]
[36]
P.-G. Eriksson and A. J. Isaksson, "Some aspects of control loop performance monitoring," in Control Applications, 1994., Proceedings of the Third IEEE Conference on, 1994, pp. 10291034: IEEE. S. P. Boyd, C. H. Barratt, S. P. Boyd, and S. P. Boyd, Linear controller design: limits of performance. Prentice Hall Englewood Cliffs, NJ, 1991. B. S. Ko and T. F. Edgar, "PID control performance assessment: The single‐loop case," AIChE Journal, vol. 50, no. 6, pp. 1211-1218, 2004. A. Visioli, Practical PID control. Springer Science & Business Media, 2006. V. Kariwala, "Fundamental limitation on achievable decentralized performance," Automatica, vol. 43, no. 10, pp. 1849-1854, 2007. A. Y. Sendjaja and V. Kariwala, "Achievable PID performance using sums of squares programming," Journal of Process Control, vol. 19, no. 6, pp. 1061-1065, 2009. R. Fu, L. Xie, Z. Song, and Y. Cheng, "PID control performance assessment using iterative convex programming," Journal of Process Control, vol. 22, no. 9, pp. 1793-1799, 2012. F. Shahni and G. Malwatkar, "Assessment minimum output variance with PID controllers," Journal of Process Control, vol. 21, no. 4, pp. 678-681, 2011. M. Veronesi and A. Visioli, "Global minimum-variance PID control," in Proceedings IFAC World Congress, 2011, pp. 7891-7896. F. Shahni, W. Yu, and B. Young, "Monitoring and diagnosis of PI controllers," in Control Conference (AUCC), 2015 5th Australian, 2015, pp. 1-5: IEEE. K. Åström, "Introduction to stochastic control theory. 1970," ed: Academic Press, New York. S. J. Qin, "Control performance monitoring—a review and assessment," Computers & Chemical Engineering, vol. 23, no. 2, pp. 173-186, 1998. L. G. Bergh and J. F. MacGregor, "Constrained minimum variance controllers: internal model structure and robustness properties," Industrial & engineering chemistry research, vol. 26, no. 8, pp. 1558-1564, 1987. J. Hanna, S. Upreti, A. Lohi, and F. Ein-Mozaffari, "Constrained minimum variance control using hybrid genetic algorithm–An industrial experience," Journal of Process Control, vol. 18, no. 1, pp. 36-44, 2008. K. Byung-Su and T. F. Edgar, "Assessment of achievable PI control performance for linear processes with dead time," in American Control Conference, 1998. Proceedings of the 1998, 1998, vol. 3, pp. 1548-1552 vol.3. N. Stanfelj, T. E. Marlin, and J. F. MacGregor, "Monitoring and diagnosing process control performance: the single-loop case," Industrial & Engineering Chemistry Research, vol. 32, no. 2, pp. 301-314, 1993.
527
25
Figure
H (q 1)
Gc (q 1)
Figure 1. Block diagram of a typical single-loop feedback control system
Figure 2. An illustration of how m is selected for the model system, G
1
q 5 . q 0.8
26
Impulse Response
x 10 8 6
Amplitude
4 2 0 -2 -4 -6 100
200
300
400
500
600
700
Time (seconds)
Figure 3. Example 1, where insufficient size of impulse length causes the system to become unstable. 25
4
Impulse Response
x 10
3 2
Amplitude
1 0 -1 -2 -3 -4
0
200
400
600
800
1000
1200
1400
1600
1800
Time (seconds)
Figure 4. Example 2, impulse response due to insufficient m.
2
2000
2.45
2.4496
2.4492
2.4488
2PID
2.4484
2.448
2.4476
2.4472
2.4468
2.4464
2.446 0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
Figure 5. MOV versus µ variation simulated for case study 15 from Table 3.
Figure 6. Even a poor selection of finite horizon gives a reasonable result.
3
Figure 7. The benefits of the weighting method with a poor selection of finite horizon.
Figure 8. Scaled optimization results based on the size of the impulse response.
4
Table
Table 1. Good and bad tuning based on the good and bad selection of finite impulse response size.
Methods Optimum tuning from our method Tuning from prior methods [24-28]
1T 1 0.4234
T2 2 0.0012
Tr r 0
Eq. (23) 1
0.4225
0.1015
infinity
0
Table 2. Best versus weak tuning.
Methods Optimum tuning from this proposed method tuning from prior methods
1T 1 2T 2
S
2.2980
2.2986
2.0971
1
2 2
No 1
Table 3. Some process and disturbance case studies taken from literature for comparison. G N Reference 0.2q
5
1 0.8q
2
1
4 5 6 7
8 9 10 11 12 13 14
0.08919
[16]
1 0.8669q 1
[34]
0.5108q 28 1 0.9604q 1
0.5108 1 0.9604q 1
[34]
q 6 1 0.8q 1
1 0.2q 1 (1 0.3q 1 )(1 0.4q 1 )(1 0.5q 1 )
[35]
q 6 1 0.8q 1
1 0.6q 1 (1 0.5q 1 )(1 0.6q 1 )(1 0.7q 1 )
[35]
q 6 1 0.8q 1
1 0.2q 1 (1 q 1 )(1 0.3q 1 )(1 0.4q 1 )(1 0.5q 1 )
[35]
q 6 1 0.8q 1
1 0.6q 1 (1 q )(1 0.5q 1 )(1 0.6q 1 )(1 0.7q 1 )
0.1q 5 1 0.8q 1
0.1 (1 q )(1 0.3q 1 )(1 0.6q 1 )
[35]
0.1q 3 1 0.8q 1
1 1 q 1
[23]
0.1q 6 1 0.8q 1
0.1 (1 q )(1 0.7q 1 )
[23]
0.1q 3 1 0.8q 1
0.001 (1 q )(1 0.2q 1 )
[36]
0.004679(1 0.9355q 1 )q 7 (1 .923q 1 )(1 .887q 1 )
0.07889(1 0.946q 1 )q 1 (1 q 1 )(1 0.8465q 1 ) 1 1 q 1
1
1
1
1
1
(0.2644 0.05q 1 )q 2 (1 0.8423q 1 )(1 0.003q 1 ) (0.5285 0.2707q 1 )q 2
[35]
[29] --
(1 0.375q )(1 0.3606q )
1 1 q 1
--
(0.6806 0.4834q 1 )q 2 (1 0.7859q 1 0.3679q 2 )
1 1 q 1
--
1
15
(1 q )(1 0.4q 1 )
0.08919q 12 1 0.8669q
3
1
1
1
2
Table 4. Comparison of PMR (Proposed Method’s Results) with previous studies for MVPID Prediction.
Case
Ref.
Ref.
Ref.
Ref.
Ref.
Ref.
[23]
[26]
[27]
[28]
[29]
[32]
MV
PMR
PMR
Without weighting
With weighting
Optimal controller
3.0728
3.0728
[2.8455;4.4166;1.7547]
0.0312
0.0310
[2.0266;3.7866;1.7600]
3.0256
3.0239
[0.5077;0.9839;0.4762]
1.1121
1.1121
1.1121
[0.0144;-0.0253 ;0.0109]
1
2.9427
3.0730
3.0728
3.0728
2
0.0306
0.0310
0.0311
0.0310
3
3.0112
3.0495
3.0442
3.0238
4
1.1121
5
3.4004
3.4065
3.4081
3.4088
3.4400
3.4106
3.4065
[0.1359;0.2527;0.1167]
6
11.9528
13.8243
13.8077
13.8077
17.75
13.8076
13.8076
[0.7241;1.2057;0.5178]
7
58.3406
89.6983
87.738
87.7519
123.54
87.7380
87.7377
[0.8306;1.3963;0.6074]
8
0.2978
0.4278
0.4247
0.4247
0.4247
0.4246
[8.0887;13.1802;5.5890]
9
3
3.2093
3.2032
3.2032
3.2032
3.2032
[6.5347;9.2397;3.3593]
10
0.3144
0.4288
0.4268
0.4268
0.4268
0.4268
[8.2159;13.7482;5.9545]
11
0.0023
0.0025
0.0024
0.0024
0.0024
0.0024
[6.1493;8.5376;3.0087]
12
1.0026
Infinity
4.7264
4.7270
4.7256
4.7256
[23.5491;44.6910;21.2945]
13
2
16.8678
2.1740
2.1740
[2.4471;3.2740;1.0699]
14
2
6.0002
2.2985
2.2985
[0.9815;1.2276;0.4788]
15
2
5.9892
2.4463
2.4463
[0.5950;0.6983;0.3244]
0.4564
*
*
*
*
4.2610
6.1283
*
*
5.0012
*
2.5804
*
5.0798
*
3.5791
3.1821
0.0312
*
*
*. Indicates we calculate the result by the method proposed in the literature.
3
Table 5. Calculation times of different methods. Updated Reported calculation time Case study calculation For examples 1-11 12-15 time <2s <2s <30s <40s Less than (< )100 s <90s ** <668s 4-514 s <476s** <564s 5-25 min <22min** <23min
Methods Our proposed method Method with weighing parameter [26] [27] [28]
** indicates without consideration of examples 12-15, second order systems.
Table 6. Effect of the selection of matrix size (m) by the percentage of system stability on optimum value. Case study
1%
2%
5%
10%
20%
30%
40%
50%
1
3.0728
3.0728
3.0728
3.0728
3.0728
3.0733
3.0797
3.0854
2
0.0312
0.0314
0.0314
0.0314
0.0314
0.0314
0.0314
0.0314
3
3.0239
3.0297
3.0314
3.0366
3.0386
3.0511
3.0686
3.0963
4
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
5
3.4106
3.4112
3.4137
3.4171
3.4239
3.4353
3.4482
3.4811
6
13.8076
13.8076
13.8078
13.8088
13.8096
13.9858
14.4394
∞
7
87.7380
87.7391
87.7582
87.8001
116.1792
99.9348
∞
∞
8
0.4247
0.4247
0.4247
0.4248
0.4250
0.4285
0.4698
∞
9
3.2032
3.2032
3.2032
3.2032
3.1798
3.2034
3.2036
3.2093
10
0.4268
0.4268
0.4270
0.4270
0.4274
0.4542
∞
∞
11
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
12
4.7256
4.7256
4.7256
4.7257
4.7258
4.7258
4.7258
4.7385
13
2.1740
2.1740
2.1740
2.1740
2.1741
2.1741
2.1741
2.1749
14
2.2985
2.2986
2.2986
2.2990
2.3642
2.3642
2.3642
∞
15
2.4463
2.4465
2.4466
2.4478
∞
∞
∞
∞
4
Table 7. Appropriate matrix order, m (at 99% and 90% of the process finite impulse response steady state) versus the delay of the system and disturbance steady state. Case
Delay
99%
90%
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
5 12 28 6 6 6 6 5 3 6 3 7 2 2 2
50 88 282 52 52 52 52 50 46 52 46 154 58 16 20
30 56 168 32 32 32 32 30 26 32 26 94 30 10 8
Time for the disturbance to achieve steady state 6 30 100 8 13 8 16 10 0* 15 4 30 0* 9 10
Note: Zero here is due to it being a step disturbance.
5
PII: DOI: Reference:
S0019-0578(18)30426-9 https://doi.org/10.1016/j.isatra.2018.10.047 ISATRA 2948
To appear in:
ISA Transactions
Received date : 1 March 2018 Revised date : 10 August 2018 Accepted date : 31 October 2018 Please cite this article as: Shahni F., Yu W. and Young B. Rapid estimation of PID minimum variance. ISA Transactions (2018), https://doi.org/10.1016/j.isatra.2018.10.047 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.
Rapid estimation of PID minimum variance Farzam Shahni, Wei Yu, Brent Young
Chemical and Materials Engineering Department, The University of Auckland, New Zealand.
Corresponding author: Brent R Young Professor of Food & Process Systems Engineering Head of Chemical & Materials Engineering Chair of the Food & Health Programme University of Auckland T +64 9 923 5606 Email: [email protected] A Room G61A, 2-6 Park Avenue Grafton, Auckland, NZ 1023
*Highlights (for review)
The highlights: 1. Assessment of PID achievable performance for linear, stationary SISO processes. 2. Guidelines developed for choosing the appropriate size of the finite impulse response for the system. 3. Fast non-iterative method to estimate the PID benchmark.
*Blinded Manuscript - without Author Details Click here to view linked References
1
2
Abstract
3
PID controllers are the most common type of controllers used in industry. However, to find
4
the best PID performance is a challenging task. Assessment of the best PID performance
5
leads to the minimization of the deviation from set point as indicated by closed loop transfer
6
function. This is a non-convex optimisation problem and there is no direct solution for
7
finding this benchmark. Solving this problem iteratively as proposed in literature requires
8
long calculation times. Furthermore, iterative methods are not guaranteed to find the
9
minimum value. Using as many as possible infinite impulse response coefficients, or
10
increasing the number of finite impulse coefficients gradually, causes undesirably large
11
calculation times. We propose a fast method to evaluate the minimum variance by a fixed
12
number of finite impulse response coefficients. This finite impulse model is double the size
13
of the stable system impulse response (it contains the transient response and the same
14
duration of steady state). Thus it avoids undesired iteration, and consequently, the PID
15
minimum variance is quickly evaluated. Time is a critical factor as the purpose to assess this
16
index is for use in online monitoring. The proposed method is tested on benchmark
17
simulation
examples
from
1
literature.
18
1
Introduction
19
Surveys show that almost 60% of industrial control loops are working with poor performance
20
(60% in general, could be due to controller problem, final control element problems or
21
sensors problem, linear or non-linear, feedback or feedforward) and only 40% of processes
22
are working under desired optimal performance [1, 2]. Assessment of control loop
23
performance has been a topic of interest for many researchers [3-10]. Performance
24
assessment is mainly used to verify the working performance of one controller with
25
consideration of input and output variance. This assessment compares the current controller
26
performance and the best achievable performance. In the assessment of control loop
27
performance, factors such as time delay, disturbance and process models, and structural
28
constraints of the controller affect the system output performance. For a linear, discrete-time
29
system with additive disturbances, the lowest output variance is achieved when a minimum
30
variance (MV) controller is implemented. This MV controller becomes the standard
31
benchmark to assess control loop performance [11]. Even though using MV control as a
32
performance benchmark can provide reasonable performance, this controller is not practical
33
due to its demand for excessive control effort and poor robustness [6, 12]. Also, MV is not
34
the only and a complete objective in the reality of process control.
35 36
There is no doubt that more than 90% of the applied controllers in industry are implemented
37
with PID algorithms [7, 13-15]. Some even believe 98% of controllers in industry are PID
38
type [16, 17]. It is reported that 80% of PID controllers in industry are poorly tuned [16, 17]
39
(major reason of the 60% poor performance in general loops is due to the lack of the
40
resonable PID controllers tuning).
41
From the above-noted reports, it can be concluded that more than 70% of all controller bad
42
performance in industry is due to bad tuning of PID parameters. When the desired controller 2
43
is a PID type, performance indices based on MV controllers as presented in, e.g. [18, 19] are
44
not suitable to tell whether the PID type controller can provide good controller performance.
45
PID controllers have both order and structural constraints which result in the possibility of
46
not achieving MV performance, and thus the best performance of the practical controller is
47
unknown. As Kozub reported in [20], only 20% of loops in a refinery with PID controller
48
structure may achieve the MV control performance lower bound. Eriksson and Isaksson in
49
[21], proposed an index with an optimal PID controller as a benchmark to replace MV.
50
Hence, it is important to assess the minimum variance with PID control (MVPID). Although
51
the benchmark controller can provide reasonable performance, benchmark controllers are not
52
practical due to their demand for excessive control effort and poor robustness and in practice
53
a less aggressive controller will be selected.
54
The problem of assessing the MVPID produces a non-convex problem, and a globally optimal
55
solution cannot be guaranteed [22]. There is no direct and clear method to solve this complex
56
non-convex problem. Researchers have tried to solve the non-convex problem by various
57
convex methods. A gradient-based method was proposed by Ko and Edgar [23, 24]. An
58
analytical lower bound method was proposed by Kariwala [25], using the sum of the squares
59
approach proposed by Sendjaja and Kariwala [26]. Other similar approaches have also been
60
proposed [27-29].
61
There are two main drawbacks with the current methods of solving the MVPID assessment
62
problem: The first is calculation time (avoiding iteration), and the second is an accurate
63
estimation. Time is a critical factor as the purpose to evaluate this index is for use in online
64
monitoring. Most of the current literature uses iteration to solve this problem which causes
65
the calculation time to rise. Even if the iteration scheme is designed to increase the chance of
66
finding the minimum value, it may provide the minimum value or it may not. In these studies,
67
there were not any theoretical methods identified for finding the suitable finite length. A large
3
68
pre-defined finite size [27, 30] or increasing the size of finite length (impulse response
69
coefficients) from a small (e.g., 2d-1, where d is delay) to a large value (e.g., 4d, the size of
70
iteration is 2d here) was proposed by [25, 28]. These methods find the minimum from the
71
finite set and from the results (controller values) it calculates the variance of the closed-loop
72
system in each step. The results are saved and compared with the next sized data set. This
73
procedure causes calculation time to increase unnecessarily. The unclear, random selection
74
and insufficiently truncated the finite impulse response is the reason why these methods miss
75
the optimum time and optimum value.
76
This paper proposes a fast method to estimate MVPID which is capable of being used in online
77
monitoring. This technique reduces calculation times by using a fixed finite length and
78
removing the iteration. It produces reliable results as it is compared with the existing
79
methods. Due to the fast computation time, the method could be implemented for online
80
performance assessment. A full theoretical analysis is presented, and a simulation study is
81
performed for the assessments using the Matlab optimization toolbox. In addition to the fast
82
method proposed, we introduce a weighting method for more accurate results and as an
83
alternative solution for these nonconvex problems. This weighting method changes the
84
weight of the transient and the steady state impulse response and balances these two parts in
85
the cost function.
86
It is important to note that this is a simulation study, and constraints like manipulated variable
87
limitations are not considered. This is common practise in such studies and not usually a
88
major concern for process control as the plant dynamics are usually much slower than the
89
actuator dynamics. At this stage use of this method is also limited due to its demand for
90
information regarding the process and disturbance model. Thus it is not a feasible method for
91
loops where appropriate models do not exist or can be obtained and we would recommend to
92
use the proposed method only on important loops where models exist or are able to be
4
93
identified. The scope of this work is for linear plants and models, or any plants that have mild
94
nonlinearity which is well handled by PI and PID feedback (the majority of industrial process
95
plants).
96
This paper is organized as follows: Section 2 introduces the problem of assessing the
97
minimum variance of a PID controller. Section 3 consists of the analytical formulation of the
98
PID minimum variance calculation. Section 4 consists of the proposed methods. Section 5
99
and its subsections involve simulation study, examples and comparison. In this section,
100
fifteen examples have been used to verify the usefulness of the proposed method. Finally, this
101
paper is then summarized along with some concluding remarks in the last section.
102
2
103
A block diagram of a typical single-input-single-output (SISO) feedback control system is
104
shown in Figure 1 with output y (t ) , manipulated variable u (t ) , and zero mean white noise
105
a (t ) , where t is the sampling interval of the one-time unit. The process output can be written
106
as:
Best Achievable PID performance
y(t ) G(q 1 )u(t ) H (q 1 )a(t )
107
(1)
G(q 1 ) and H (q 1 ) represent the process and disturbance transfer functions,
108
where
109
respectively. For simplicity and brevity, the backward shift term q , q
110
omitted in the subsequent sections unless circumstances necessitate its presence.
1
1
y(t ) y(t 1) , is
111 112
When there is no set point change, the output y (t ) under closed loop conditions (with set point
113
zero), can be expressed by:
5
y (t )
114
H a(t ) Gcl a(t ) 1 G Gc
(2)
115
116
Gc is digital PID controller with the structure of: Gc
117
k1 k2 q 1 k3q 2 1 q 1
(3)
k1 = kP +kI +kD , k2 (kP 2kD ) , k 3 k D . k P , kI and k D denote proportional, integral
118
Where
119
and derivative gain, respectively. Gcl is defined as the closed loop transfer function between
120
y (t ) and a (t ) . The output variance under the PID controller is:
Var ( y ) Gcl 2 . a2 2
121
(4)
122
a2 represents the variance of the random noise
123
all its poles and zeroes inside the unit circle thus G contains no zeroes and poles on or
124
outside the unit circle except at infinity (time delay) and H contains no zeroes on or outside
125
the unit circle [31]. In this work, the variance of the white noise sequence ( a ) is fixed to be
126
unity.
127
The system output in Eq.(1) could be written in its impulse response form as follows:
a (t ) . For solving Eq. (4) we suppose Gcl has
2
i 1
i 1
y (t ) gi u (t i ) hi a (t i )
128
(5)
gi and hi are the impulse response coefficients of the process transfer function and
129
Where,
130
disturbance model, respectively. By substituting for the PID controller given in Eq. (3), the
131
output from the Eq. (5) becomes:
6
i 1
i 1
y (t ) gi Gc . y (t i ) hi a (t i )
132
133
By considering the definition of the backward shift operator ( q
134
rewrite Eq. (6) as follows:
1
(6)
y(t ) y(t 1) ), we can
y (t )
135
h a(t i) i 0
i
1 2 i 1 (k1 k2 q k3q ) si (q ) i 1
136
Where si represents the step response coefficient of the process transfer function. We can
137
write the impulse response of disturbance
139
1
0 D 1 0
1
0 0
(9)
The denominator of Eq. (7) can then be written in matrix form as follows:
S I (k1 k2 D 1 k3 D 2 )( s1D 1 s2 D 2 s3 D 3 ...)
142 143
(8)
In matrix form the backward shift operator ( q ) can be written as:
140
141
Hi (h0 , h1 , h2 , h3 ,...) in matrix form as follows:
h0 h h0 h1q 1 h2 q 2 h3q 3 ... 1 h2
138
(7)
(10)
Where I is the identity matrix. We then rewrite Eq. (7) as follows:
y (t ) i a (t i )
144
i 0
145
and
7
(11)
0 h 1 2 S
146
(12)
147
Where is the impulse response vector of the closed loop transfer function. The i ’s are
148
h0 h the closed loop impulse response coefficients of the SISO system and h 1 is the impulse h2
149
response vector of the disturbance transfer function. The output variance is defined as:
y2 i 2 a2
150
(13)
i 0
151
The minimum output variance with an unrestricted controller structure is known as the MV
152
benchmark [11] and is calculated as: d 1
2 MV hi2 a2
153
(14)
i 0
154
where d is the process time delay. If the controller structure must be PID, the achievable
155
minimum variance low bound MVPID is not equal to MV and cannot easily be found. To
156
minimize the output variance with a PID controller structure, one needs to find the
157
k1 , k2 and k3 values in Eq. (15), so that y2 reaches the minimum value. This leads to an
158
optimization problem which can be expressed as:
159
2
min y2 min
k1 , k2 , k3
k1 , k2 , k3
i 1
8
2 i
a2 min Gcl 2 a2 2
k1 , k2 , k3
(15)
160
Where . 2 denotes the 2 norm . Eq. (15) yields a non-convex optimization problem. There is
161
no direct and simple solution for such an equation [22]. Be noted that here it is considered
162
that a set point is constant in value and purpose is to minimize the effect of noise.
163
3
164
The H 2 norm of a stable transfer matrix S is:
Problem formulation 2
165
Gcl
2 2
( T )
(16)
166
From Eq.(2), due to the controller parameters being unknown as stated in Eq.(3), is
167
unidentified. In general is obtained from a finite impulse response [32]. However, it is
168
essential that the finite impulse response contains all non-zero coefficients. Here, to obtain
169
up to the m th element, the finite samples of the closed loop transfer function Gc l are
170
expressed by the following,
171
0 m 1 Sm 1 h m m1
(17)
172
Where hm is defined as the impulse response of the disturbance model H up to the m th
173
element. The term hm can be represented as,
174 175
hm (h0 , h1 , h2 ,..., hm1 )1Tm From Eq.(10) and (17), the truncated matrix S, is defined as:
9
(18)
1 s 1 S m s2 sm 1
176
177
If
m
is taken to be very large then
178
179
0
0 0
s1
0
sm 2
sm 3
0 0 0 1 mm
(19)
m and then from Eq.(16) and Eq.(17) we have: Gcl
2 2
( T ) mT m
(20)
In other words, we have the equation below:
0
180
mT m 1 ( T )
(21)
181
When sufficient impulse length is selected, the Eq. (21) becomes equal to unity, and if m is
182
very short Eq. (21) becomes zero. Directly solving Eq.(15) is impossible, instead Eq. (20) is
183
used to find the approximate minimum output variance for PID controllers [23-29, 32]: min y2 min mT m a2
184
PID
m
PID
(22)
is limited to the range 4d m 8d due to computational problems (the
185
In [28] the size of
186
iteration is 4d). Another way is to increase
187
to reach the convergence of the left-hand side (upper bound) and right-hand side (lower
188
bound) in Eq. (20) [26]. In this case the size of iteration could be anything between 2 to
189
infinity, however, in the program one should consider a limitation due to lack of computer
190
memories like 30 or 40. A pre-specified constant,
191
these problems. The authors here intentionally avoid iteration and instead increase the size of
192
the matrix which again causes the calculation time to increase.
m
starting from the analytical lower bound 2d -1
10
m (m=500) is also used in [27] to address
193
4 The proposed method
194
is separated into two parts, m , thus Eq. (21) can be rewritten as: r
mT m mT m 1 T T T m m m m r r r r
195
(23)
196
From Eq. (23) it can be concluded r r 0 . This term has led previous research to neglect
197
the system impulse response after
198
found by solving the optimisation problem in Eq.(22). For validation of the solution the
199
controller parameters which are found through Eq. (22) are used in Eq. (16). At this stage, if
200
the size of
201
Then the value found in the previous step is saved and is used to compare with the result of
202
the next iteration. In this case, some previous research gradually increased the matrix size, m
203
in Eq.(17). Gradually increase the size of
204
calculation time to increase and most results cause Eq. (16) (the upper bound), to become
205
infinitely large or far from the minimum amount. Some tried to choose a large
206
output variance as reported in [23-29, 32]. However, a pre-specified constant, large
207
similarly can lead to an unnecessary rise in calculation time.
208
4.1 Fast method
209
By increasing the size of the finite horizon the cost function Eq. (16) increases. However
210
from some point (appropriate m ) this increase in the cost function should stop (when it hits
211
the lower bound or the exact optimum) or becomes negligible (finds a local optimum). Now
212
we define the denominator of Eq. (17) as:
213
T
m is
m samples.
In the literature, the minimum variance is
not appropriate, the result will be far from the actual minimum variance.
m from a small value via trial and error, causes the
A A1 S m C 11
0 A mm
m to calculate
m
(24)
A1 consists of step response coefficients.
214
As is shown by Eq. (10) and Eq.(19) block matrix
215
Here m is selected to be whenever the step response of a process is within 99% of the final
216
value (to see the reason behind the selection of 99%, Section 5.1). The reason behind the
217
above selection of m is to have approximately constant coefficients in C. It is mentioned
218
before that the 2-norm is equal to the root-mean-square of the impulse response of the
219
system, and the variance is the squared of the 2-norm. Thus the first
220
contain the final minimum variance and the second
221
optimum variance. We consider the block matrix A as the transient part of the step response
222
and C as it’s approximately steady state part (the elements of being defined as the steady
223
state step response coefficients multiplied by the controller’s gains). The following
224
illustration Figure 2 is an example of the system model, G
225
of the matrix, m, is selected in this paper.
2
m impulse coefficients 2
m impulse coefficients contain the 2
q 5 , to clarify how the size q 0.8
226 227
The inverse of the lower triangular matrix
A1 A11 ( Sm ) 1 1 1 A CA
228
229
230
A1 is given below:
The matrices Aand
0 A1 mm
(25)
A1 are lower triangular matrices. Now matrix m is calculated as:
A1 m 1 A11 h m A1CA1 2
12
A1 h m2 0 nm A1 (hc CA1 h m ) A1 2 mm
(26)
231
h m Where h m2 and h c are defined as h m 2 , the impulse response samples of the disturbance hc
232
hm 2 h m 1 before and after term h m2 , respectively ( hc 2 ). If we increase the size of the matrix in h m1
233
Eq. (26) to 2m , 2m becomes:
2 m m S 2 m 1 h 2 m 2
234
235
1
h m2 hc hc hc
(27)
C C . Then 2m becomes: C C mm
2m
h m 1 A11nm 2 2 0 hc m h (28) 1 1 1 h c A11 hc A1 C1 A1 2 A1 3 hc hc 4 h c
A11 A 1C A 1 1 1 1
Then we have 3 , 4 as: m hc 3 1 1 h 2 A1 C1 A1 4 hc hc
238 239
0 0 0 A
We define the matrix C1
236
237
A 0 0 C A 0 C C A C C C
(29)
where, 3 A1 (CA1 I )(hc CA1 h ) m 2
240
and
4 A1CA1 (CA1 I )(h c CA1 h ) A1 (CA1 I )(nc CA1 h ) m 2
m 2
A1 (CA1 I ) 2 (h c CA1 h m2 )
13
(30)
241
We can continue this calculation and see the effect of the term (CA1 I )(nc CA1nm ) on future 2
242
prediction. 4m could be written as follows:
4m
243
244
A1 0 0 C A1 0 1 C1 C1 A1 C1 C1 C1
A1 h m2 h m2 1 1 1 A (hc CA h m2 ) hc 1 2 1 1 1 m ) A ( CA I )( h CA h c 3 2 0 hc A1 (CA1 I ) 2 (h CA1 h m ) 0 hc c 2 4 0 hc A1 (CA1 I )3 (h c CA1 h m ) 5 2 A1 hc A1 (CA1 I ) 4 (h CA1 h m ) 6 c 2 hc A1 (CA1 I )5 (hc CA1 h m ) 7 2 hc 8 A1 (CA1 I )6 (h c CA1 h m2 )
(31)
Eq. (16) can be rewritten as below: 2
245
Gcl ( T ) 1T 1 T2 2 T3 3 T4 4 T5 5 T6 6 T7 7 T8 8 ...
246
(32)
247
Thus the equation below is repeated over time (all rows from Eq. (31) include the term
248
below);
2
249
(33)
(CA1 I )(hc CA1 h m2 )
250
From the above equation (33), when the disturbance is assumed to be a step in the output
251
1 n ( the common case for a disturbance [33]), the impulse response of the closed-loop 1
252
system in each period
m is: 2
1 i A (CA I ) h m2 A (CA I ) 1 1
253
1
i 1
1
1
i 1
(34)
254
The impact of the
255
solution of Eq. (22). Thus for solving Eq. (15), we consider the first two terms of Eq. (32):
(hc CA1 h m2 )
term in Eq. (31) means we should consider the term for the
14
min S 2 a2 min( 1T 1 2T 2 ) a2 2
256
PID
m is
(35)
PID
selected as “large enough” and by chance may cover the
257
In previous research when
258
terms 1 , 2 , 3 ,... , the control parameters may be tuned well and achieve comparable results,
259
but purely by chance. An unnecessary large m causes extra calculation and extra time is
260
required for finding reliable optimum results. A small
261
To find the solution of Eq. (15) by solving Eq. (22), previous research [26-29] only focused
262
on the first
263
consider the important effect of the remaining samples. Thus the selection of impulse length
264
according to the criteria proposed in the first part of this section is appropriate, and it avoids
265
unnecessary iteration. In Figure 3 and Figure 4 it is shown that what will happen if we just
266
consider a short length of the finite impulse response and do not consider its future effects. In
267
these figures, the size of the matrix, m, was not appropriate. Regardless, the prior approach
268
tried to minimize the variance based on length m. It is depicted in these figures that selecting
269
finite length impulse response could minimize the right-hand side of Eq. (22). In these cases
270
in which the truncated length of the impulse response is not appropriate, the Eq. (21) is 0.
271
Example 1: By way of example we considered a system with the process model
272
G
273
impulse response; we compared 1 1 , 2 2 in Eq. (35) and saw its effects on Eq. (23),
274
based on optimum and non-optimum results, Table 1.
275
Table 1 shows how the selection of the wrong sized matrix (m) could affect the minimization
276
problem - the result is far away from the optimum value. In these cases, the short length of
277
the impulse response is considered without knowledge of its future.
m samples,
m
also leads to non-optimal results.
where they used iteration to find a suitable m , and they did not
0.1 0.1q 5 , and disturbance model H . Figure 3 is the 1 1 (1 q )(1 0.3q 1 )(1 0.6q 1 ) 1 0.8q T
T
15
278 279
Example2: In this example (Table 2) we consider a system with the process model
280
G
281
response of the system. Here we want to show where Eq. (35) is valid based on a comparison
282
of 1 1 , 2 2 . We also compare the left-hand side of the equation (the upper bound or
283
actual norm which is based on infinite impulse response) and right-hand side (lower bound
284
which is based on finite impulse response). These two bounds should be merged to validate
285
Eq. (35). Table 2 shows how bad tuning results due to wrong selection of impulse response
286
length from prior methods affects the Eq. (35).
287
As reported in literature for the non-convex optimization problems in general there is no
288
guaranteed method [22]. However, we have proved throughout Eq. (24) to Eq. (35) in Section
289
4 that if we select the correct length of finite horizon (which we proposed in section 4.1) the
290
method can guarantee a feasible solution and optimum value just for this problem. However
291
for more accuracy the weighting method is proposed.
292
4.2 Weighting method
293
When our objective is minimization, where one is concerned with obtaining more accurate
294
results, the impact or weight of the
295
greater than A1 h m2 , the first term in Eq. (22). Therefore, for more accuracy, the weight of the
296
first impulse response is reduced when solving Eq. (15). Hence, the optimization problem is
297
as follows,
298
(0.5285 0.2707q 1 )q 2 1
1
(1 0.375q )(1 0.3606q )
T
, and disturbance model H
1 . Figure 4 is the impulse 1 q 1
T
(hc CA1 h m2 )
term in a solution of Eq. (31), must be
min y2 min S 2 a2 min( 1T 1 2T 2 ) a2 2
PID
PID
PID
16
(36)
299
µ is a weighting coefficient which should be any value less than 1. When µ is small, the result
300
is far away from the minimum value. When µ is close to 1, we can tune to achieve an
301
optimum result. The scale in Figure 5 is enlarged to show how the weighting parameter
302
affects the optimization problem in Eq.(36). However, when µ is close to 1, the optimization
303
result is not too far from the minimum possible result. In previous studies, in some iterations,
304
the controller parameters found through Eq.(22), caused the variance of the system to become
305
infinite (Eq. (21) became zero) due to improper
306
small size and is gradually increased).The difference between Eq. (22) and Eq. (36) could be
307
shown as follows:
1
308
309
The right-hand side of the above equation,
Eq.(15) .
1T 1 2T 2 Gcl
2
(as mostly the size of
m
starts off being of
1T 1 2T 2 mT m 2 2 Gcl Gcl 2 2 mT m Gcl
310
m
2
(37)
, shows Eq. (22) is an under-estimation for
2
shows Eq. (36) is an unbiased estimation. It should be noted that
2
311
small leads to the result of the optimization problem in Eq.(15) to reach a large value for
312
1T 1 and smaller values for 2T 2 . The result should always be a balance between these
313
two parts of Eq. (36) to reach to an optimum value.
314 315
A detailed simulation study has been performed for assessment of the proposed method using
316
the MATLAB package with the optimization toolbox. We used large scale algorithm based
317
on interior-reflective Newton method. Any optimization function could be called to find this
318
local minimum of constrained nonlinear multivariable polynomials like the nonlinear
319
fmincon
function
in
MATLAB
17
with
using
option
settings
320
'TolCon ',1e 6,'TolX ',1e 6,'TolFun ',1e 6 . This optimization function finds a minimum
321
of a constrained nonlinear multivariable scalar function, but it is highly dependent on initial
322
values. To ensure for positive values of the controller parameters, inequalities k P
323
and
324
included in the analysis. For solving the optimization problem, the first step is to start the
325
optimization with an initial value of x0
326
order
327
calculation time drops from minutes to a few seconds.
328
This part is added as an alternative solution when the size of the finite horizon is not selected
329
properly. For better explanation please see Figures 6 and 7 below.
330 331
5 Simulation examples
332
models. However, if the system model is not known, system identification methods may be
333
used to obtain the open-loop process model, and by using output samples based on the
334
installed controller, the closed-loop disturbance model can be identified [7, 23]. We justify
335
this approach as follows. Even though no exact model can ever be obtained, for most
336
chemical processes system identification can find useful and appropriate models that are
337
sufficiently accurate and capture the essential dynamics that are necessary for this approach.
338
The method has been tested for more than 300 different scenarios and processes, for which
339
this method worked properly. From these 300 cases, 15 representative examples were
340
selected for the purpose of comparison with the previous methods. Computational
341
performance on various tests show the advantages of our approach.
0 , kI 0
kD 0 are considered. Thus inequality constraints such as k1 0 , k2 0 and k3 0 are
m
[0;0;0] . In this research, by selecting a fixed matrix
, which is described in Section 4, and by removing the undesired iteration, the
The proposed method is based on complete knowledge of the process and disturbance
18
342
Table 3 shows examples of process and disturbance models. These first order plus time delay
343
(FOPTD) and second order plus time delay (SOPTD) examples demonstrate the efficiency of
344
the proposed algorithm for general application.
345
The MV, results of other researchers, MVPID and optimal PID parameters obtained by the
346
proposed methods are given in In Table 2. In total 15 case studies were compared, and it was
347
found that our method either equalled or outperformed all the other previous methods [26-
348
30]. The Table 4 shows the weighting method is not case sensitive. It can be concluded that
349
the method is reliable in general and it can be used as a standard solution reference or
350
benchmark. In no cases did the proposed method perform inferiorly to previous studies, albeit
351
in some cases it yielded the same results as other approaches. It can be deduced that the
352
proposed method is more reliable as it always produced better results or performed the same
353
as other studies. The cases in which the method performed like others have been highlighted
354
with green colour in Table 4. The low accuracy results are highlighted with a red colour. As
355
is shown in Table 4 the previous methods do not perform well for higher order systems.
356
As expected, mostly there is a gap between MVPID and MV (as explained before in the
357
introduction section, this gap is due to limitations based on the fixed order and structure of
358
PID controllers). It is precisely indicated how the optimal behaviour of industrial controllers
359
(PI/PIDs) differ from MV. Examples 4-7 in Table 3 have the same process models, but they
360
have different disturbance transfer functions. The results of these examples in Table 2, show
361
the expected gap between MVPID and MV increases when disturbance changes from auto-
362
regressive moving average (ARMA) to auto-regressive integrated moving average (ARIMA).
363
The simulation study has been performed using an Intel Core i7 3.40 GHz CPU, 16 GB Ram
364
PC. Calculation time for the proposed method is two seconds as shown in Table 5 which is
365
much faster than prior methods. The second column of Table 5 is the reported calculation
366
times of the other methods [26-28]. The reported times in Table 5 have been updated by 19
367
repeating the simulation for all 15 cases using the same machine used for the proposed
368
method to have a fair comparison. The last column in this table is the calculation times of
369
second-order systems which show all the previous methods are time-consuming compared to
370
the proposed method. For second-order systems, as is shown there is almost a 600 seconds
371
gap between the proposed method and the previous methods. These updated times are shown
372
in the third column of Table5.
373
Table 5 shows the proposed method is fast and could be used in the online assessment. Based
374
on the comparison in Table 5, the method without weighting (the fast method) shows
375
reasonable results with just less than 0.005% deviation from the best possible results as
376
compared with weighting method (the standard solution) and the current methods proposed in
377
the literature. The simulation results verify the capability of the proposed method which
378
could be implemented and acceptable for online PID control performance assessment.
379
According to the obtained results, this method is the best method to achieve MVPID as shown
380
in Table 4.
381
Because of the simplicity of the method and the speed, the method could be potentially used
382
for general purposes and for higher order systems (Examples 12-15 in Table 3 are second-
383
order systems). Application of the method does not need high-performance computers for
384
complicated systems. It is worth mentioning that this method is also applicable to find PI
385
minimum variance. To find the PI minimum variance, in Eq.(3) we must set k 3 k D 0 , so
386
Eq.(3) becomes K PI
387
may be worthwhile to say here that in the presence of non-stationary disturbances, PID
388
controllers exhibit better performance than PI controllers, due to the prediction capability of a
389
derivative action in PID controllers. Our experience in this work study indicates that our
390
proposed method represents an efficient way to find minimum variance problem by PIDs.
k1 k2 q 1 and the rest of the procedure remains unchanged. However, it 1 q 1
20
391
5.1 Size of finite impulse
392
In Table 6 we show how the selection of matrix size (finite impulse length) affects the optimum value
393
of Eq.(36). As is shown in this table, there is almost no difference in the result of the optimization for
394
the
395
However, if m is selected to be only 50% of the final value, the optimization result is far
396
away from the optimum, and in some cases it becomes infinity. Figure 8 shows how the result
397
will be affected by the size of m, based on the percentage of final value. In Figure 8, the 15
398
different symbols belong to the optimization results of the 15 different case studies. To scale
399
values and compare all 15 examples, we divided all optimization results for each case by
400
numerical results (the matrix size m is approach to infinity). Figure 8 shows why we selected
401
the size of the finite impulse response to be based upon the approximation of the final value
402
as being within 99% to 90% of the steady-state step response in section 4.1. In this figure, 1
403
indicates the optimum value for the all 15 cases.
404
We investigate if there is almost the same result obtained if we select finite impulse response
405
coefficient (m) to be within 99% to 50% of the final value of the steady state step response.
406
The first four points are very similar however, if we take a closer look at the Figure8 we can
407
see the scores at the 99% level are the smallest. The proposed method is tested for more than
408
300 different scenarios and always worked properly.
409
5.2 Effect of delay and disturbance in size of the finite horizon:
410
Here we examine the probable impact of the delay and disturbance of system into the
411
appropriate finite impulse. Some researchers like [24, 26], suggested the start point for
412
optimization is related to delay of the system. They start the length of the finite horizon from
413
the double size of delay (2d-1) optimization. Table 7 shows there is no relation between the
414
proper size of m and the delay of the system. Also, a pre-specified large m (m=500 or
m
selected within 99% to 90% of final value.
21
415
m=1000) as proposed by other researchers [27, 30] could cause unnecessary long calculation
416
time. For example, for case 9, size of the finite impulse response coefficient should be 26
417
while [27, 30] considered the size of m be 500 or 1000. Thus the above mentioned approved
418
we cannot predefine m by chance, and it is case base. In Table 7 we also consider the
419
possibility any relation between the size of m and disturbance dynamic. It is shown there is
420
no relation between the size of the matrix, m, and disturbance dynamics.
22
421
6
Conclusions
422
This paper has proposed a new method to find the performance lower bound for PID type
423
controllers and single-input-single-output (SISO) systems to be used for online control
424
performance assessment. The calculation time for this method dropped to less than two
425
seconds which is capable it to be used for online performance assessment. The method
426
consists two parts. The first part selects the length of impulse response coefficients of the
427
closed loop transfer function between the output of and a disturbance to the system, which
428
causes the method, be fast and simple. The second part is optimization with or without a
429
weighted polynomial function arising from the impulse response. The purpose of the second
430
part is producing precise results. The simplicity of the method enables this method to be used
431
for general problems of higher order. It has been shown from theoretical and simulation
432
studies that the proposed algorithm is accurate and fast. The method is able to be extended to
433
find the achievable PID performance for the general class of all applications. This research
434
covers the most currently used controllers in the industry (PID and PI) and could identify
435
how a current controller is far from an ideal controller with the same performance or action.
436
As shown, the method is fast, so it potentially can be used for MIMO systems. However, it is
437
worthwhile noting here that finding the MV benchmark is not the only step in implementing
438
real process control. The next phase of this work is focused on extending the method to
439
MIMO processes. The most fundamental constraint in process control is the time delay for
440
SISO systems which are time invariant. Since the actual action of time invariance for MIMO
441
systems is done by an interactor matrix, the two main challenges here are large dimensional
442
matrixes which cause numerical errors and computational times that need to be short.
443
23
References
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491
[1] [2]
[3] [4]
[5]
[6]
[7] [8] [9] [10] [11] [12] [13]
[14]
[15] [16]
[17] [18] [19] [20]
T. Samad and A. Annaswamy. (2014). The Impact of Control Technology (2nd ed.). Available: www.ieeecss.org T. Harris, C. Seppala, and L. Desborough, "A review of performance monitoring and assessment techniques for univariate and multivariate control systems," Journal of Process Control, vol. 9, no. 1, pp. 1-17, 1999. Y. Alipouri and J. Poshtan, "Minimum variance lower bound estimation and realization for desired structures," ISA Transactions, vol. 53, no. 3, pp. 787-792, 2014/05/01/ 2014. A. Silveira, R. Trentini, A. Coelho, R. Kutzner, and L. Hofmann, "Generalized minimum variance control under long-range prediction horizon setups," ISA Transactions, vol. 62, pp. 325-332, 2016/05/01/ 2016. Z. Yu and J. Wang, "Performance assessment of static lead-lag feedforward controllers for disturbance rejection in PID control loops," ISA Transactions, vol. 64, pp. 67-76, 2016/09/01/ 2016. B. Huang, "A pragmatic approach towards assessment of control loop performance," International Journal of Adaptive Control and Signal Processing, vol. 17, no. 7‐9, pp. 589-608, 2003. M. Jelali, "Control performance management in industrial automation," Advances in Industrial Control, DOI, vol. 10, pp. 978-1, 2013. B. Huang and S. L. Shah, Performance assessment of control loops: theory and applications. Springer Science & Business Media, 1999. M. Jelali, "An overview of control performance assessment technology and industrial applications," Control Engineering Practice, vol. 14, no. 5, pp. 441-466, 2006. S. Huang and C. Huang, "Optimal PID tuning for discrete-time systems with stochastic disturbance," in Control Conference (CCC), 2014 33rd Chinese, 2014, pp. 9002-9007: IEEE. T. J. Harris, "Assessment of control loop performance," The Canadian Journal of Chemical Engineering, vol. 67, no. 5, pp. 856-861, 1989. N. F. Thornhill and B. Huang, "Management and monitoring of process assets," Advanced control of industrial processes ADCONIP, pp. 4-7, 2008. L. Desborough and R. Miller, "Increasing customer value of industrial control performance monitoring-Honeywell's experience," in AIChE symposium series, 2002, pp. 169-189: New York; American Institute of Chemical Engineers; 1998. L. Desborough and R. Miller, "Increasing customer value of industrial control performance monitoring-Honeywell's experience " presented at the Chemical Process Control-VI, Tucson, Arizona, 2001. K. J. Åström and T. Hägglund, Advanced PID Control. ISA - The Instrumentation, Systems and Automation Society, 2006. P. Agrawal and S. Lakshminarayanan, "Tuning proportional-integral-derivative controllers using achievable performance indices," Industrial & engineering chemistry research, vol. 42, no. 22, pp. 5576-5582, 2003. P. Van Overschee and B. De Moor, "RAPID: The End of Heuristic PID Tuning," IFAC Proceedings Volumes, vol. 33, no. 4, pp. 595-600, 2000/04/01/ 2000. L. Desborough and T. Harris, "Performance assessment measures for univariate feedback control," The Canadian Journal of Chemical Engineering, vol. 70, no. 6, pp. 1186-1197, 1992. D. J. Kozub and C. E. Garcia, "Monitoring and diagnosis of automated controllers in the chemical process industries," in AIChE annual meeting, 1993, vol. 9. D. J. Kozub, "Controller performance monitoring and diagnosis: experiences and challenges," in Chemical process control conference, Lake Tahoe, USA, 1996, pp. 83-96.
24
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526
[21]
[22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]
[34]
[35]
[36]
P.-G. Eriksson and A. J. Isaksson, "Some aspects of control loop performance monitoring," in Control Applications, 1994., Proceedings of the Third IEEE Conference on, 1994, pp. 10291034: IEEE. S. P. Boyd, C. H. Barratt, S. P. Boyd, and S. P. Boyd, Linear controller design: limits of performance. Prentice Hall Englewood Cliffs, NJ, 1991. B. S. Ko and T. F. Edgar, "PID control performance assessment: The single‐loop case," AIChE Journal, vol. 50, no. 6, pp. 1211-1218, 2004. A. Visioli, Practical PID control. Springer Science & Business Media, 2006. V. Kariwala, "Fundamental limitation on achievable decentralized performance," Automatica, vol. 43, no. 10, pp. 1849-1854, 2007. A. Y. Sendjaja and V. Kariwala, "Achievable PID performance using sums of squares programming," Journal of Process Control, vol. 19, no. 6, pp. 1061-1065, 2009. R. Fu, L. Xie, Z. Song, and Y. Cheng, "PID control performance assessment using iterative convex programming," Journal of Process Control, vol. 22, no. 9, pp. 1793-1799, 2012. F. Shahni and G. Malwatkar, "Assessment minimum output variance with PID controllers," Journal of Process Control, vol. 21, no. 4, pp. 678-681, 2011. M. Veronesi and A. Visioli, "Global minimum-variance PID control," in Proceedings IFAC World Congress, 2011, pp. 7891-7896. F. Shahni, W. Yu, and B. Young, "Monitoring and diagnosis of PI controllers," in Control Conference (AUCC), 2015 5th Australian, 2015, pp. 1-5: IEEE. K. Åström, "Introduction to stochastic control theory. 1970," ed: Academic Press, New York. S. J. Qin, "Control performance monitoring—a review and assessment," Computers & Chemical Engineering, vol. 23, no. 2, pp. 173-186, 1998. L. G. Bergh and J. F. MacGregor, "Constrained minimum variance controllers: internal model structure and robustness properties," Industrial & engineering chemistry research, vol. 26, no. 8, pp. 1558-1564, 1987. J. Hanna, S. Upreti, A. Lohi, and F. Ein-Mozaffari, "Constrained minimum variance control using hybrid genetic algorithm–An industrial experience," Journal of Process Control, vol. 18, no. 1, pp. 36-44, 2008. K. Byung-Su and T. F. Edgar, "Assessment of achievable PI control performance for linear processes with dead time," in American Control Conference, 1998. Proceedings of the 1998, 1998, vol. 3, pp. 1548-1552 vol.3. N. Stanfelj, T. E. Marlin, and J. F. MacGregor, "Monitoring and diagnosing process control performance: the single-loop case," Industrial & Engineering Chemistry Research, vol. 32, no. 2, pp. 301-314, 1993.
527
25
Figure
H (q 1)
Gc (q 1)
Figure 1. Block diagram of a typical single-loop feedback control system
Figure 2. An illustration of how m is selected for the model system, G
1
q 5 . q 0.8
26
Impulse Response
x 10 8 6
Amplitude
4 2 0 -2 -4 -6 100
200
300
400
500
600
700
Time (seconds)
Figure 3. Example 1, where insufficient size of impulse length causes the system to become unstable. 25
4
Impulse Response
x 10
3 2
Amplitude
1 0 -1 -2 -3 -4
0
200
400
600
800
1000
1200
1400
1600
1800
Time (seconds)
Figure 4. Example 2, impulse response due to insufficient m.
2
2000
2.45
2.4496
2.4492
2.4488
2PID
2.4484
2.448
2.4476
2.4472
2.4468
2.4464
2.446 0.6
0.65
0.7
0.75
0.8
0.85
0.9
0.95
1
Figure 5. MOV versus µ variation simulated for case study 15 from Table 3.
Figure 6. Even a poor selection of finite horizon gives a reasonable result.
3
Figure 7. The benefits of the weighting method with a poor selection of finite horizon.
Figure 8. Scaled optimization results based on the size of the impulse response.
4
Table
Table 1. Good and bad tuning based on the good and bad selection of finite impulse response size.
Methods Optimum tuning from our method Tuning from prior methods [24-28]
1T 1 0.4234
T2 2 0.0012
Tr r 0
Eq. (23) 1
0.4225
0.1015
infinity
0
Table 2. Best versus weak tuning.
Methods Optimum tuning from this proposed method tuning from prior methods
1T 1 2T 2
S
2.2980
2.2986
2.0971
1
2 2
No 1
Table 3. Some process and disturbance case studies taken from literature for comparison. G N Reference 0.2q
5
1 0.8q
2
1
4 5 6 7
8 9 10 11 12 13 14
0.08919
[16]
1 0.8669q 1
[34]
0.5108q 28 1 0.9604q 1
0.5108 1 0.9604q 1
[34]
q 6 1 0.8q 1
1 0.2q 1 (1 0.3q 1 )(1 0.4q 1 )(1 0.5q 1 )
[35]
q 6 1 0.8q 1
1 0.6q 1 (1 0.5q 1 )(1 0.6q 1 )(1 0.7q 1 )
[35]
q 6 1 0.8q 1
1 0.2q 1 (1 q 1 )(1 0.3q 1 )(1 0.4q 1 )(1 0.5q 1 )
[35]
q 6 1 0.8q 1
1 0.6q 1 (1 q )(1 0.5q 1 )(1 0.6q 1 )(1 0.7q 1 )
0.1q 5 1 0.8q 1
0.1 (1 q )(1 0.3q 1 )(1 0.6q 1 )
[35]
0.1q 3 1 0.8q 1
1 1 q 1
[23]
0.1q 6 1 0.8q 1
0.1 (1 q )(1 0.7q 1 )
[23]
0.1q 3 1 0.8q 1
0.001 (1 q )(1 0.2q 1 )
[36]
0.004679(1 0.9355q 1 )q 7 (1 .923q 1 )(1 .887q 1 )
0.07889(1 0.946q 1 )q 1 (1 q 1 )(1 0.8465q 1 ) 1 1 q 1
1
1
1
1
1
(0.2644 0.05q 1 )q 2 (1 0.8423q 1 )(1 0.003q 1 ) (0.5285 0.2707q 1 )q 2
[35]
[29] --
(1 0.375q )(1 0.3606q )
1 1 q 1
--
(0.6806 0.4834q 1 )q 2 (1 0.7859q 1 0.3679q 2 )
1 1 q 1
--
1
15
(1 q )(1 0.4q 1 )
0.08919q 12 1 0.8669q
3
1
1
1
2
Table 4. Comparison of PMR (Proposed Method’s Results) with previous studies for MVPID Prediction.
Case
Ref.
Ref.
Ref.
Ref.
Ref.
Ref.
[23]
[26]
[27]
[28]
[29]
[32]
MV
PMR
PMR
Without weighting
With weighting
Optimal controller
3.0728
3.0728
[2.8455;4.4166;1.7547]
0.0312
0.0310
[2.0266;3.7866;1.7600]
3.0256
3.0239
[0.5077;0.9839;0.4762]
1.1121
1.1121
1.1121
[0.0144;-0.0253 ;0.0109]
1
2.9427
3.0730
3.0728
3.0728
2
0.0306
0.0310
0.0311
0.0310
3
3.0112
3.0495
3.0442
3.0238
4
1.1121
5
3.4004
3.4065
3.4081
3.4088
3.4400
3.4106
3.4065
[0.1359;0.2527;0.1167]
6
11.9528
13.8243
13.8077
13.8077
17.75
13.8076
13.8076
[0.7241;1.2057;0.5178]
7
58.3406
89.6983
87.738
87.7519
123.54
87.7380
87.7377
[0.8306;1.3963;0.6074]
8
0.2978
0.4278
0.4247
0.4247
0.4247
0.4246
[8.0887;13.1802;5.5890]
9
3
3.2093
3.2032
3.2032
3.2032
3.2032
[6.5347;9.2397;3.3593]
10
0.3144
0.4288
0.4268
0.4268
0.4268
0.4268
[8.2159;13.7482;5.9545]
11
0.0023
0.0025
0.0024
0.0024
0.0024
0.0024
[6.1493;8.5376;3.0087]
12
1.0026
Infinity
4.7264
4.7270
4.7256
4.7256
[23.5491;44.6910;21.2945]
13
2
16.8678
2.1740
2.1740
[2.4471;3.2740;1.0699]
14
2
6.0002
2.2985
2.2985
[0.9815;1.2276;0.4788]
15
2
5.9892
2.4463
2.4463
[0.5950;0.6983;0.3244]
0.4564
*
*
*
*
4.2610
6.1283
*
*
5.0012
*
2.5804
*
5.0798
*
3.5791
3.1821
0.0312
*
*
*. Indicates we calculate the result by the method proposed in the literature.
3
Table 5. Calculation times of different methods. Updated Reported calculation time Case study calculation For examples 1-11 12-15 time <2s <2s <30s <40s Less than (< )100 s <90s ** <668s 4-514 s <476s** <564s 5-25 min <22min** <23min
Methods Our proposed method Method with weighing parameter [26] [27] [28]
** indicates without consideration of examples 12-15, second order systems.
Table 6. Effect of the selection of matrix size (m) by the percentage of system stability on optimum value. Case study
1%
2%
5%
10%
20%
30%
40%
50%
1
3.0728
3.0728
3.0728
3.0728
3.0728
3.0733
3.0797
3.0854
2
0.0312
0.0314
0.0314
0.0314
0.0314
0.0314
0.0314
0.0314
3
3.0239
3.0297
3.0314
3.0366
3.0386
3.0511
3.0686
3.0963
4
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
1.0121
5
3.4106
3.4112
3.4137
3.4171
3.4239
3.4353
3.4482
3.4811
6
13.8076
13.8076
13.8078
13.8088
13.8096
13.9858
14.4394
∞
7
87.7380
87.7391
87.7582
87.8001
116.1792
99.9348
∞
∞
8
0.4247
0.4247
0.4247
0.4248
0.4250
0.4285
0.4698
∞
9
3.2032
3.2032
3.2032
3.2032
3.1798
3.2034
3.2036
3.2093
10
0.4268
0.4268
0.4270
0.4270
0.4274
0.4542
∞
∞
11
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
0.0024
12
4.7256
4.7256
4.7256
4.7257
4.7258
4.7258
4.7258
4.7385
13
2.1740
2.1740
2.1740
2.1740
2.1741
2.1741
2.1741
2.1749
14
2.2985
2.2986
2.2986
2.2990
2.3642
2.3642
2.3642
∞
15
2.4463
2.4465
2.4466
2.4478
∞
∞
∞
∞
4
Table 7. Appropriate matrix order, m (at 99% and 90% of the process finite impulse response steady state) versus the delay of the system and disturbance steady state. Case
Delay
99%
90%
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
5 12 28 6 6 6 6 5 3 6 3 7 2 2 2
50 88 282 52 52 52 52 50 46 52 46 154 58 16 20
30 56 168 32 32 32 32 30 26 32 26 94 30 10 8
Time for the disturbance to achieve steady state 6 30 100 8 13 8 16 10 0* 15 4 30 0* 9 10
Note: Zero here is due to it being a step disturbance.
5