Design of non-parametric process-specific optimal tuning rules for PID control of flow loops ?
Igor Boiko a a The
Petroleum Institute, P.O. Box 2533, Abu Dhabi, U.A.E.
.
Abstract The problem of designing optimal process-specific rules for non-parametric tuning is undertaken in the paper. It is shown that producing non-parametric processspecific optimal tuning rules for PID controllers leads to the problem that can be characterized as optimization under uncertainty. This happens due to the fact that tuning rules, unlike tuning constants, are produced not for a particular process or plant model but for a set of models from a certain domain. The novelty of the proposed approach is that the problem of obtaining optimal tuning rules for a flow process is formulated and solved as a problem of optimization of an integral performance criterion parametrized through values that define the domain of available process models. The considered non-parametric tuning assumes the use of the modified relay feedback (MRFT) recently proposed in the literature. It allows one to tune PID controller satisfying the requirements to gain or phase margins that is achieved through coordinated selection of tuning rules and test parameters. This approach constitutes a holistic approach to tuning. In the present paper, optimal tuning rules coupled with the MRFT test, for flow loops, are proposed. Final results are presented in the form of tables containing coefficients of optimal tuning rules for the PI controller, obtained for a number of specified gain margins. The produced non-parametric tuning rules well agree with the practice of loop tuning. Key words: Parametric optimization; PID control; flow process; tuning.
? Corresponding author I. Boiko. Tel. +971-2607-5505. Fax +971-2607-5194. Email address:
[email protected] (Igor Boiko).
Preprint submitted to Elsevier Science
28 October 2013
1
Introduction
Despite the availability of variety of control methods the proportional-integralderivative (PID) control remains the main control used in the process industry. Some research shows that the share of PID controllers in a typical overall plant control system is about 97%. Tuning of PID controllers remains a practically important and theoretically interesting problem. Two approaches to PID controller tuning are parametric, which is based on identification of underlying model of the process, and non-parametric that is not based on any specific model of the process. In non-parametric tuning, rules are based on measurement some important characteristics of a process such as ultimate gain and frequency proposed by Ziegler and Nichols [1], or amplitude and frequency of self-excited oscillations [2], [4]. There are a number of methods that use the same tests as [1] or [2] but different tuning rules (e.g. [3], [5] and others). Various tuning rules are designed for for integrating and non-integrating )self-regulating) processes, processes with long time delay, processes with inverse action, etc. A thorough overview of methods of tuning can be found in [6], [7], [8]. Non-parametric methods require some specially designed tuning rules that relate characteristics measured in a test over the process with tuning parameters of the controller. The very popular tuning methods and rules [1], [2] are of semi-empirical nature because they can guarantee specified stability only in the case of proportional-only controller ([6], [4]). In the case of PI or PID control, strictly speaking, there is no theoretical justification of provided stability margins. In method [4], this problem is resolved and the designed test and tuning rules can guarantee specified stability gain or phase margins for PID, PI and PD controllers, too. This is done through coordinated selection of a specific tuning rule and specific parameters of the proposed test named modified relay feedback test (MRFT), which is in fact a certain discontinuous control algorithm that can provide either lead or lag in the switching in comparison with conventional relay feedback test (RFT) [2]. Satisfaction to the specified gain of phase margin is achieved through introduction of an equality-type constraint that relates parameters of the tuning rules. Yet availability of 2 (for PI or PD) or 3 (for PID) parameters that define the tuning rule leaves some freedom in selection of these parameters because many different combinations can ensure the same stability margins. Therefore, some optimization of the tuning rules would be possible. The problem of finding optimal parameters that define tuning rules for a specific process, which will be referred to as design of process-specific optimal tuning rules, is undertaken in the present research. The goal of this approach is to design not universal tuning rules but tuning rules suitable for a specific process, for example separately for flow, pressure, level and temperature loops. The assumption is that there is more similarity between different flow loops 2
than between a flow loop and a level loop, for example, and therefore, tuning rules aimed at serving a narrower class would be more precise. The advantage of this approach is in having more suitable tuning rules for each class of processes and through this to ensure more accurate tuning. The problem of finding optimal tuning rules for a class of processes is not a conventional parametric optimization problem. What makes it unconventional is the fact that the parameters of the process model, for which the tuning rule is going to be used, are not known in advance. Yet for some classes of processes respective implied process models can be obtained along with ranges of parameters for these models. In the described situation, we have to deal with parametric optimization under uncertainty in which a cost function contains some parameters with unknown values. In other words the cost function is parametrized through variables, that describe the situation of the application of a tuning rule, having unknown values but with known ranges of possible variations. We shall call these parameters situational because they describe a particular situation of the application of a tuning rule.
In the parametric optimization theory, parametrized cost functions are considered in a number of publications in sensitivity analysis of nonlinear programming problems [9], [10], [11]. Continuity analysis of optima with parametrized cost functions was given by Berge [12]. Close to the considered problem robust optimization is usually associated with linear and nonlinear programming problem in which uncertainty is a part of the constraint [13], [14], [15] and not of the cost function. The problem of parametric optimization with a cost function having an uncertain parameter was considered in [16], [17] and other publications. It is usually referred to as optimization under uncertainty. An overview of methods of solution of this problem is presented in [18]. However, usually uncertainty in the referenced and other publications is related to several scenarios that may develop in the future but not to a continuous design parameter. To the best knowledge of the author, the specific approach developed in the paper is not presented in the literature.
The paper is organized as follows. First, an example that illustrates the considered formulation of optimization is given. Then the optimization problem is formulated rigorously. Third, the model of the flow process and the MRFT are presented. And finally the optimization problem is solved for the coefficients that define tuning rules for a PID controller and, thus, process-specific optimal tuning rules are produced. 3
2
Example illustrating the problem of engineering design under uncertainty
In publications on optimization under uncertainty [16], [17], [18] and others, uncertainty is related to possible development several different scenarios in the future. Both deterministic and stochastic approaches are used. To illustrate the problem of engineering design under uncertainty, in which uncertainty is associated with continuous parameter, and develop a deterministic robust solution to this problem we consider the following example of engineering design. Consider the problem of design of a bucket of the cylindric shape made of 1 square meter 1 of sheet metal, which can hold a maximum amount (volume) of content, subject to the condition that it can be used for water, sand, gravel, cement, etc. Of these substances, only water (or other liquid) can be poured up to the hight of the bucket. All other materials may be put higher than the height of the bucket (see Fig. 1). In other words, capacity of the bucket depends on the material that is put in it.
Fig. 1. Bucket holding dry material
We denote the height of the bucket as H > 0 and the radius of its round bottom as R > 0. Assume that all granular materials can fill the bucket making a cone-like shape at the top up to a certain angle α without running over. This angle will give the maximum amount of a granular material, which 1
We disregard the shape of the piece of sheet metal but are only concerned with the amount of material used (area of the bucket surface)
4
the bucket can hold (Fig. 1), that we refer to as the capacity of the bucket. We also assume that angle α is different for each material and can lie within the range of [0, αmax ]. Obviously, the capacity of the bucket is found as a sum of the cylindric part and the cone part as follows: 1 V = πR3 tan α + πR2 H 3
(1)
The constraint on the amount of the sheet metal used to make the bucket is given by the equation that shows the sum of the areas of the round bottom and tubular wall: 2πRH + πR2 = 1
(2)
We find the hight from (2) and substitute into (1) thus obtaining a criterion for unconstrained optimization problem: V = πR
3
1 1 1 + R tan α − 3 2 2
(3)
function V (R) is We can note that if the coefficient 13 tan α − 12 > 0 then √ monotone increasing and its maximum is reached at R = 1/ π, when H = 0. We further assume that 13 tan α − 12 < 0 or, equivalently, tan α < 32 , which results in α < 56.31◦ . We find the optimal solution by equating the derivative dV /dR to zero and solving the respective equation, which gives: R∗ = q
1
(4)
π(3 − 2 tan α)
√ If, for example, the √ angle α = 0 (water) then R∗ = 1/ 3π and the optimal hight H ∗ = R∗ = 1/ 3π. For some other values of α, the optimal radius and the volume of material held are given in Table 1. We should note that despite presenting a few discrete values we treat parameter α as a continuous variable. One can see that maximum volume of the material held by the bucket depends not only on the bucket design but also on the situation in which it is used. This situation is characterized by parameter α. We shall call it a situational parameter. Therefore, we have a case when the cost function is parametrized by a situational parameter and this parameter may take values from a certain (known) range. In other words α ∈ D, where D is the domain of possible variation of a situational parameter(s). The obtained results (Table 1) do not, however, allow us to make a decision about optimal bucket design. 5
Table 1 Optimal radius and maximum volume for different α α
4
8
12
16
20
R∗
0.334
0.342
0.352
0.362
0.374
V (R∗ )
0.111
0.114
0.117
0.121
0.125
f∗ = 1/V (R∗ )
8.99
8.77
8.53
8.28
8.02
We come to the situation of decision making, which can be done using minmax of probabilistic approach. Yet we should note that in the classic formulation of robust optimization the uncertainty is contained in the constraint allowing for application of Wald’s minmax approach. In the example being considered the uncertainty is contained in the parameter by which the cost function is parametrized. However, the minmax approach can be used in the considered problem too, which requires reformulation of the cost function aimed at compensation for the effects of the situational parameter. The following approach can be used. First, instead of maximization of the volume, we shall consider minimization of the cost function f (R) = 1/V (R) optimal values of which are given in Table 1. This will match the problem of finding optimal tuning rules considered below. Second, we need to cancel the effect of the situational parameter on the cost function by including a weight inversely proportional to the value of f ∗ . And third, another cost function must be used. Maximum deterioration of the cost function f (R) from the optimal value f ∗ due to the use of R that does not coincide with R∗ for this α can be used as a new cost function: (
f (R, α) g(R) = max α∈D f ∗ (α)
)
(5)
For every α, the following inequality holds: f (R, α) ≥ f ∗ (α), which follows from the meaning of f ∗ (α), and (5) provides the worst deterioration of optimality on the domain of situational parameters. The results of calculation of f (R, α) ≥ f ∗ (α) for different R and α are presented in Table 1. The largest value in each row is shown in bold font. For the considered example (5) would mean by how much the capacity of the bucket designed for a particular α would be smaller if it is used to hold a different material. In other words, if our bucket is optimally designed for sand, how much we loose in its capacity (in %) if we are going to use it for water, gravel, etc. Minimization of this criterion would constitute the solution of the optimization problem. The optimal solution is, therefore, the one that gives the smallest value of those that are given in bold font in Table 1. Obviously, the optimal radius in the considered example is R = 0.352 because in the worst 6
Table 2 Values of
f (R,α) f ∗ (α)
α=0
α=4
α=8
α = 12
α = 16
α = 20
R= 0.326
1.0000
1.0008
1.0033
1.0078
1.0147
1.0245
R= 0.334
1.0010
1.0000
1.0008
1.0037
1.0089
1.0170
R= 0.342
1.0038
1.0010
1.0000
1.0011
1.0046
1.0110
R= 0.352
1.0101
1.0047
1.0013
1.0000
1.0012
1.0052
R= 0.362
1.0197
1.0113
1.0052
1.0013
1.0000
1.0016
R= 0.374
1.0358
1.0234
1.0136
1.0063
1.0016
1.0000
scenario (the use of this bucket, that is optimally designed for α = 12◦ , for water for which α = 0) it gives only 1.01% of performance (bucket capacity) deterioration, while all other options give higher performance deterioration.
3
Formulation of robust optimization problem with parametrized cost function
As we saw in the example provided, a general formulation of the optimization problem involves finding the minimum of function f (x, θ). However, because of the dependence on parameter θ there is uncertainty in such a formulation. The requirement of finding the minimum of f (x, θ), which would suit the “best” to all θ from a certain range, can be formulated as the following problem of optimization under uncertainty (robust optimization). Minimum of a cost function f (x, θ) parametrized through a vector of situational parameters θ ∈ D can be found through the solution of the following problem: (
)
f (x, θ) , minimize g(x) = max θ∈D f ∗ (θ)
(6)
subject to P (x, θ) ≤ 0
(7) 7
and Q(x, θ) = 0,
(8)
where x ∈ X ⊂ Rn is a vector of decision variables, θ ∈ D ⊂ Rm is a vector of situational parameters, f (x, θ) ∈ R is a parametrized cost function, g(x) is a cost function on domain D of situational parameters, f ∗ (θ) is the solution of the minimization of f (x, θ) subject to constraints (7) and (8), for a given θ, which will be referred to as a minimizer cost for a given θ. Either of (7) or (8) or even both of these constrains may not be present. In the latter case we get a problem of unconstrained optimization. If the situational parameter θ were fixed we would formulate the problem of optimization as finding minimum of function f subject to (7) and (8). As one can see from the example given above, the situational parameter has certain effect on the optimal solution. Therefore, if we ignore this and treat situational parameters as additional decision variables we would arrive at the solution, which will favorable to specific situational parameters values. To cancel the effect of the situational parameter we use a transformation that allows us make the cost function invariant to the situational parameter in the point of minimum: ft (x, θ) =
f (x, θ) , f ∗ (θ)
(9)
However, ft (x, θ) cannot be used as a cost function because it has the whole set of minima points, with values equal to 1. Instead, function g(x) (6) has to be used. It is worth noting that, unlike function f , function g does not contain argument θ and applies to the whole domain D of the situational parameters. Minimizer costs f ∗ (θ) are found a-priori through minimization of f (x, θ) subject to (7) and (8) for every θ from D. The presented criterion of optimization is, therefore, a minimum of deviation from optimality in terms of cost f in the worst case scenario. Minimization of g(x) ensures robustness of solution because the most suitable solution for the whole domain of θ is sought for.
4
Modified relay feedback test and holistic approach to test and tuning
We stated above that optimal tuning rules can be found through the solution of parametric optimization problem with uncertainty concerning parameters 8
that participate in the cost function. These uncertain parameters are parameters of the model of the process. Uncertainty comes from the fact that flow process may be the process of liquid flow or gaseous flow, may be controlled by modulating a control valve position or changing the speed of a pump or compressor, etc. The variety of possibilities result in the fact that parameters of the process model may vary significantly. Fortunately, the model structure (type of equations) can normally be used for the whole variety, so that the differences are only parametric. 4.1
Modified relay feedback test
The objective of the present study is to find optimal tuning rules that can be used together with the modified relay feedback test (MRFT). The MRFT was proposed in publication [19]. Later [4] it was extended to the option of having a phase margin constraint and was used in the industrial software [20]. The approach of [4] involves a test, which allows one to excite oscillations in the third quadrant of the process Nyquist plot, and a coordinated selection of test and tuning parameters. It is briefly summarized below as follows. Consider the following discontinuous control: u(t) =
h if e(t) ≥ b1 or (e(t) > −b2 and u(t−) = h)
−h
if e(t) ≤ b2 or (e(t) < −b1 and u(t−) = −h)
(10)
where b1 = βemin , b2 = −βemax , emax > 0, emin > 0 are last “singular” points of the error signal corresponding to the last maximum and minimum values of e(t) after crossing the zero level, u(t−) = lim→0,>0 u(t − ) is the control value at the time immediately preceding current time t, h is the amplitude of the relay, β is a positive constant. The error signal is e(t) = r(t) − y(t), where r is the reference and y is the system output (process variable). Initial value of emax or emin can be assigned as e(0) (the choice of either emax or emin depends on the sign of e(0) ). The algorithm (10) is similar to the so-called “generalized sub-optimal” algorithm used for generating a second-order sliding mode in systems of relative degree two [21]. It was shown in [4] that if the reference signal r (t) is zero then motions in the with control (10) and process given by transfer function Wp (s) are periodic. Because the motion in the system with MRFT is periodic, the describing function (DF) method [22] can be used for finding the amplitude and the frequency of these oscillations. The DF of the the MRFT is given in [4]: 4h q 2 1 − β − jβ N (a) = πa
(11)
9
Parameters of the oscillations can be found from the harmonic balance equation: Wp (jΩ0 ) = −
1 N (a0 )
(12)
where a 0 and Ω0 are the amplitude and the frequency of the periodic motions. The negative reciprocal of the DF is given as follows: 1 πa q 2 − 1 − β + jβ =− N (a) 4h
(13)
Finding a periodic solution in a system with control (10) has a simple graphic interpretation (Fig. 4.1) as finding the point of intersection of the Nyquist plot of the process and of the negative reciprocal of the DF, which is a straight line that begins in the origin and makes a counterclockwise angle ψ = arcsin β with the negative part of the real axis. The condition of the existence of a periodic solution is the location of the Nyquist plot of the process in the third quadrant of the complex plane . This condition is always satisfied in practice because of the inevitable existence of delays in process dynamics.
Fig. 2. Finding periodic solution
In the problem of analysis, frequency Ω0 and amplitude a0 are unknown variables and are found from the complex equation (12). In non-parametric tuning, Ω0 and a0 are measured from the MRFT, and on the basis of the measurements obtained tuning parameters are computed. 10
4.2
Non-parametric tuning rules for specification on gain margin
We assume that the PID controller to be tuned is given by the so-called parallel equation:
Wc (s) = Kc 1 +
1 + Td s , Ti s
(14)
with Kc being the proportional gain , Ti the integral time constant, and Td the derivative time constant. The task of tuning is, therefore, to determine for a given process the values of Kc , Ti , and Td , which provide optimal in a certain sense performance of the control system (loop). If the tuning rules are formulated in the format of the following homogeneous tuning rules, because the tuning parameters are homogeneous functions of 4h critical gain (given by πa ) and critical period (reciprocal of critical frequency), 0 Kc = c1
4h , πa0
Tic = c2
2π , Ω0
Tdc = c3
2π , Ω0
(15)
then the rules are defined just by the coefficients c 1 , c 2 , and c 3 . The idea behind the tuning rules (15) is the possibility of scaling of tuning parameters for processes having different time scales. It can be noted that if the tuning rules are given by (15) then the closed-loop system characteristics become invariant to the time constants of the process, so that if all time constants of the process were increased by the factor ρ then the critical frequency would decrease by the same factor ρ, and the product of every time constant by the critical frequency would remain unchanged. The use of the homogeneous tuning rules along with the MRFT will allow us to fully utilize the features of non-parametric tuning methods. If homogeneous tuning rules (15) are used then the frequency response of the PID controller at the frequency Ω 0 becomes Wc (jΩ0 ) = c1
4h 1 1−j + j2πc3 πa0 2πc2
(16)
Therefore, if the tuning rules are established through the choice of parameters c1 , c2 , c3 , and the test provides self-excited oscillations of the frequency Ω 0 , which will be equal to the phase cross-over frequency ωπ of the open-loop system (including the controller), then the controller phase lag at the frequency ωπ = Ω0 will depend only on the values of c2 and c3 :
ϕc (ωπ ) = arctan 2πc3 −
1 , 2πc2 11
(17)
which directly follows from formula (16) if ωπ = Ω0 . In publication [4], the relationship was derived that would allow us to tune PID controllers with specification on gain margin for the open-loop system. If the specified gain margin is given by γm > 1 (in absolute values), then the absolute values of both sides of (16) and consideration of the definition of the gain margin in connection with |Wp (jΩ0 )| =
πa0 , 4h
(18)
leads to the following equation: γ m c1
s
1 1 + 2πc3 − 2πc2
2
= 1,
(19)
which is a constraint that is complementary to the tuning rules (15). To provide the specified gain margin, the MRFT must be carried out with parameter
β = − sin ϕc (Ω0 ) = − sin arctan 2πc3 − 1 2πc3 − 2πc
= −r
2
1 1− 2πc3 − 2πc
2
2
1 2πc2
(20)
The MRFT with parameter β calculated as (20) and tuning rules (15) satisfying the constraint (19) can ensure the desired gain margin. However, (19) is an equation containing three unknown variables, which gives us a freedom to select different combinations of parameters c1 , c2 , c3 . This feature provides and opportunity for optimization.
5
Optimization criteria and method
In the previous section, we established the tuning rules in terms of providing the specified gain margin to the system. However, if the designed controller is of PI or PID type (not purely proportional) then there is some freedom in assigning coefficients c1 and c2 or c1 , c2 and c3 , for PI or PID controller, respectively. As far as the combination of the coefficients, which determines the tuning rule, satisfies equation (19) the controller that is designed per respective tuning rules, and the test that is carried out with parameter beta satisfying conditions (20), will provide the system with the specified gain margin. Yet it is quite obvious that some combinations of the coefficient satisfying (19) would provide better results in terms of step response, for example, than others. For 12
that reason the problem of finding optimal in a certain sense tuning rules, or combination of the coefficients which define the tuning rule, is important. Its solution would be very beneficial to the practice. However, it is also obvious that the solution of the optimization problem makes sense only for certain particular type of process or plant. For example, what is optimal for flow loops may provide a poor result for level loops, and so on. Despite the fact that all flow processes differ from each other, they differ from each other less than each of them differs from a level process. This happens because all flow processes are self-regulating and have similar values of the time constants, while all level processes are integrating and usually have quite a different order of the time constants compared to the typical flow process. To mathematically define the optimization problem we need to establish the criterion of optimization (cost function), determine the constraints, and select the optimization method. The constraint is given by the equation (19) from the specification on the gain margin, with respective determination of MRFT parameter β through equation (20). The criterion of optimization should preferably be chosen as a function from the domain other than the frequency domain because usually imposing constrains on one frequency-domain characteristic leaves very little room for varying the another frequency-domain characteristic. It would be more expedient and productive to select the optimization criterion from the time domain, for example. Time-domain criteria of optimality are very popular in the practice of process control and other industrial controls. The most popular are so-called integral criteria measured through the process reaction to the step input signals, such as integral absolute error (IAE) fIAE = fISE =
R∞ 2 e (t)dt 0
R∞ 0
|e(t)| dt , integral square error (ISE)
, integral time absolute error (ITAE) fIT AE =
integral time square error (ITSE) fIT SE =
R∞
R∞ 0
t |e(t)| dt ,
te2 (t)dt , and some others. It is
0
also worth noting that the step signal can be applied either as the set point or as a disturbance. And strictly speaking optimal tuning rules generated from the consideration of those two different step responses would be different. Another aspect of the considered optimization problem is that the optimization has to be done not for a certain particular process but for a set of different processes of the same kind representing possible variations of gains and time constants in this kind of the process. This circumstance involves uncertainty in the problem of optimization and requires application of the optimization approach considered above. The particular cost function, which is the most suitable for producing optimal tuning rules, will be selected below – once the model of the process is defined. 13
Given that the tuning rules are defined only through three parameters (coefficients), we have to deal with low-order optimization problem. Respective optimization algorithms are available in commercial mathematical packages. For example, function “fminsearch” is available in Matlab and can be used for finding optimal tuning rules. Normally algorithms used for unconstrained optimization are more straightforward than those used for the constrained optimization. While the former can be relatively simple by-coordinate descent or gradient descent, the latter have to include some modifications through inclusion of penalties to account for the presence of constraints. In our case, there is an equality constrain that comes from the requirement of the gain margin. Because of the simplicity of those formulas and the possibility of expressing one coefficient through other two the presence of the constraint makes the optimization problem simpler. In fact it allows us to reduce the order of the optimization problem by one and look for optimal values of only two coefficients. Moreover, in designing optimal tuning rules for the PI controller the problem can be reduced to onedimensional, which can be solved with the use of very simple algorithms, such as dichotomy, and even overcome the problem of identifying the global minimum.
6 6.1
Optimal tuning of flow loop Implied model of flow process
We shall call the process model, which is used for finding optimal tuning rules through the solution of the parametric optimization problem (nonlinear programming problem), the implied process model. The meaning of implied process model is different from that of underlying process model which is used in methods of identification and parametric tuning. The underlying process model is usually chosen from the considerations of sufficient simplicity subject to retaining some important properties of the process. The requirement of simplicity comes from limitations of the number of tests on the process, their accuracy, and so on. Therefore, there are reasons to select the underlying process model possibly simple. The implied process model is not used for identification of parameters, and can be relatively complex. The only purpose of this model is to be used for obtaining optimal tuning rules. In fact the more precise and adequate this model the more accurate tuning rules will be obtained. Therefore, there are good reasons to use more accurate models of the processes. Then, when the derived tuning rules are applied to the actual process, tuning will be more 14
accurate. A typical flow loop is given in Fig. 6.1. The following notations are used in Fig. 6.1: FT is a flow transmitter, FIC is a flow controller (of PI or PID type), I/P is a current-to-air pressure transducer, with air supplied to the pneumatic actuator, and the actuator driving the valve.
Fig. 3. Flow loop process diagram
Flow rate q through the valve is given by the following formula: q = Cv fv (l)
s
Δpv , gs
(21)
where Cv is the valve coefficient (it is very often referred to as CV of the valve), fv (l) ∈ [0, 1] is a function characterizing valve orifice pass area at a particular valve opening, l is a valve opening (lift), Δpv is a pressure drop across the valve, and gv is a specific gravity of the fluid. We consider only the linear valve characteristic with fv = l in the subsequent analysis. Despite the fact that valve characteristic provides only the theoretical dependence of the flow rate on valve opening and the actual dependence (installed characteristic) is complex we shall assume that tuning is applied in a certain operating point, with insignificant valve travel about this point, and , therefore, possible nonlinearities will not be revealed. It is worth noting that the model (21) provides only static dependence of the flow on the valve and process parameters, and the valve opening. Some sources recommend the use of the first-order model for the flow process with the transfer function K W (s) = T s+1 . It must be noted that the first-order model is unsuitable for the derivation of tuning rules because its use results in the infinite proportional gain as an optimal solution. According to [23], the model of the flow must contain a delay. The presence of delay in the model is not only a convenient approximation but also a reflection of the fact that there is a real delay in signal propagation through the transducer-actuator-valve dynamics. Because in our selection of the implied process model we are not limited by the order and, quite the contrary, are motivated to select possibly a precise 15
model, for the flow process, instead of commonly used first-order plus dead −τ s time (FOPDT) model W (s) = Te s+1 , we are going to use the second-order plus dead time (SOPDT) model. Wp (s) =
Kp e−τ s T 2 s2 + 2ξT s + 1
(22)
This order increase serves the purpose of bringing the process response in the MRFT in compliance with the used method of analysis, which is the describing function method requiring validity of the filtering hypothesis. According to this hypothesis, the output of the process in the modified RFT must be sinusoidal (in a steady mode). Obviously, it is impossible to achieve with the FOPDT model, while easy when using the SOPDT model. However, the SOPDT model has not three but four parameters: static gain, time constant, dead time (delay), and damping coefficient, which creates the situation of an additional degree of freedom. We will approach this situation considering the purpose of the use of the higher-order model, which is bringing the response of the plant in correspondence with the limitations of the analysis method, and the fact that most of processes in the process control industry are overdamped. Therefore, we select the value of the damping coefficient corresponding to a highly damped response of the SOPDT process model. The exact selection of the damping coefficient is not very important – as far as the model provides a damped and smooth (close to sinusoidal) response in the MRFT – because the model is used in the non-parametric tuning and must reproduce just certain qualitative properties of an actual process. Strictly speaking, the value of the damping coefficient that we are going to use is smaller than 1, which provides under-damped response. The choice in favor of smaller values of the damping coefficient is motivated by better filtering properties of this model. However, the damping with this value is still high enough, so that we shall refer to that model as to a damped response. Because of the high damping properties of the SOPDT model we shall refer to it as the damped second-order plus dead time (DSOPDT) model. We shall also establish some correspondence between the DSOPDT model and the FOPDT model in terms of providing the same values of the phase characteristic at the phase cross-over frequency (for RFT) or the same frequency of the oscillations at MRFT (subject to the filtering hypothesis). This would help us to quantify the ratio between the delay present in the process dynamics and the dominant time constant of these dynamics. It would also correspond to the practice of process control, where the step test remains the most wide-spread test over the process allowing one to determine the dead time and the time constant from the “s-curve” [1] of the step response. And even if the detailed process model is more complex than the FOPDT model the knowledge of the time constant, which would be an equivalent time constant of the FOPDT model, is helpful and convenient because the dynamics of many processes are quantified in terms of the dead time-to-time constant 16
ratio. Following this approach, we also find some equivalence of the FOPDT and DSOPDT models. In the conventional RFT, oscillations are generated at the phase cross-over frequency, i.o. the frequency at which the phase characteristic of the openloop system is equal to −180◦ . Therefore, for the FOPDT process given by the transfer function Wp (s) =
Kp e−τ s , Te s + 1
(23)
where Te is the time constant of the FOPDT model equivalent to the DSOPDT model, the phase cross-over frequency ωπ is found from the following equation: −τ ωπ − arctan(Te ωπ ) = −π.
(24) −τ s
And for the DSOPDT model given by the transfer function W (s) = T 2 s2e+2ξT s+1 , the equation that determines the phase cross-over frequency can be written as 2ξT ωπ −τ ωπ − arctan 1 − T 2 ωπ 2
!
= −π.
(25)
Equating left-hand sides of (24) and (24) yields Te ω π =
2ξT ωπ 1 − T 2 ωπ 2
(26)
Now, with this equivalence established, computation of T and Te can be done as follows. If what is known is the FOPDT model and it is required to find the equivalent DSOPDT model, equation (24) must be solved for ωπ first, and after that the time constant of the DSOPDT model is found through the solution of equation (26), which is a quadratic equation for T , as follows: T =
−ξ +
q
ξ 2 + Te2 ωπ 2
Te ωπ 2
.
(27)
If, vice versa, the DSOPDT model is known and the parameters of the FOPDT model have to be found then for the given τ and T , equation (25) must be solved for ωπ first, and after that the time constant Te is found from (26) as follows: Te =
2ξT . 1 − T 2 ωπ 2
(28)
17
Table 3 Equivalence between FOPDT and DSOPDT models for RFT τ /T
1.72 1.86 2.00 2.14 2.28 2.42 2.56 2.70
τ /Te
0.1
τ /T
2.84 2.98 3.20 3.54 3.68 3.97 4.39
τ /Te
0.9
0.2
1
0.3
0.4
1.2
1.4
0.5
1.5
0.6
1.7
0.7
0.8
2
Therefore, the relationship between T of the DSOPDT model and Te of the FOPDT model for the conventional RFT can be found a-priori, which is presented in the Table 6.1. It was also mentioned earlier that the use of the equivalent FOPDT model still remains convenient for assessment of the contribution of the delay and time constants. However, in MRFT the relationship between Te of the FOPDT model and T of the DSOPDT model are more complex than those established by equations (27) and (28) because the oscillations are excited not in the frequency point ωπ . The value of the MRFT parameter β is involved in the relationship between Te and T , and this relationship is, therefore, unknown a-priori but should be found in the process of optimization. For that reason it is easier to start from the values of T and later, after finding an optimal point, compute corresponding value of Te . The FOPDT model is going to be used only in terms of “equivalence” to the DSOPDT model to determine applicable τ /Te ranges. The DSOPDT model (22) still has three parameters which can be varied (damping ξ is fixed). And, therefore, the problem of finding optimal tuning rules can hardly be solved with this model due to enormously large number of possible combinations of these three situational parameters that needs to be investigated. We apply a normalization to model (22) to reduce the number of situational parameters. Introduce the normalized Laplace variable s0 = T s, which corresponds to the following time transformation: t0 = t/T . Noting also that the process gain Kp can be assumed unity and accounted for later because the loop gain is determined by the product of the controller proportional gain Kc and the process gain Kp , we can use the following normalized DSOPDT model of the flow process: Wp (s0 ) =
0
e−τ /T s s02 + 2ξs0 + 1
(29)
Model (29) has only one situational parameter τ /T , which makes the situational parameter space one-dimensional instead of three-dimensional and significantly simplifies the solution of the optimization problem. Therefore, it is necessary to compute the value of τ /T and find an optimal solution for the model (29). It is easy to show that the same optimal solution (in terms of 18
coefficients c1 , c2 and c3 ) that is found for the model (29) also applies to the model (22). This happens because coefficients c1 , c2 and c3 establish only proportionality between the frequency of the oscillations and the time constants, and between the ultimate gain (amplitude of the oscillations) and the proportional gain: the use of the normalized frequency will produce the normalized time constants but the use of the actual frequency will produce the actual time constants. The same rule also applies to the proportional gain. Considering that the flow loop can be used in two modes of operation: as a servo system when the flow should follow the varying set point, and as a regulator in which the set point is constant and the aim of the controller is compensation for the disturbances, we will produce two different sets of tuning rules. If the operating mode is servo mode then usually the task of disturbance compensation remains valid too. However, the magnitude of these disturbances may not be significant enough to adjust the tuning rules in favor of the regulator problem-based.
6.2
Tuning for set point changes
The use of model (29) for finding optimal tuning rules for the set point change (servo problem) is straightforward. This mode of operation usually occurs when the flow loop is cascade-connected to some primary loop. The system in this case is linear, and for a fixed value of τ /T optimal values of coefficients c1 , c2 and c3 are found through application of standard optimization algorithms. We further investigate the following range of the ratios of dead time to equivalent time constant for FOPDT model: τ /Te ∈ [0.1; 1.0], which covers most typical values of this parameter for the flow process [23]. The ratio τ /T for the DSOPDT model can approximately be estimated using Table 6.1. The exact relationship, however, is found in the process of optimization after the frequency of the self-excited oscillations in the test (ultimate frequency) is determined. The time constant Te can be found form the phase balance equation −τ Ω0 − arctan(Te Ω0 ) = −π + Ψ,
(30)
which is similar to the equation (24) that was valid for the conventional RFT. The angle Ψ in (30) is the one that is created by the MRFT algorithm. It is related with the parameter β of the algorithm as follows: Ψ = arctan √ β 2 . 1−β
We express Te from (30) as follows: 1 β tan π − arctan √ − τ Ω0 Te = Ω0 1 − β2 19
!
(31)
Table 4 Optimal solutions for DSOPDT model with τ /T = 1.5 and gain margin constraint γm = 2 Criterion
c1
c2
fmin
β
τ /Te
IAE
0.444
0.308
84.8
0.643
0.320
ITAE
0.438
0.288
307.1
0.663
0.349
ISE
0.452
0.336
59.2
0.534
0.293
ITSE
0.447
0.317
108.2
0.581
0.318
The parameter τ /Te can be computed after that simply by the division of τ used in the DSOPDT model by the value of Te calculated per (31). To select the optimization criterion we find optimal solutions using the four presented earlier integral performance criteria and select the most suitable one. For this purpose we find optimal solutions for the PI controller (finding optimal values of c1 and c2 ) for a certain typical value of τ /T . This value does not have to be found precisely because the results are going to be used only for the comparison of the criteria of optimization and the selection of one criterion which will further be used. However, this “typical” value still needs to be found. What is known is the range of typical values of τ /Te in terms of the FOPDT model, which allows use to estimate corresponding range of τ /T of the DSOPDT model as belonging to τ /T ∈ [1.2; 3.8]. The value of τ /T = 1.5 (that approximately corresponds to the value of τ /Te = 0.3 of the FOPDT model) is selected as representative for the comparison of different criteria and selection of the most suitable one. Optimal solutions for this point with the use of the IAE, ISE, ITAE and ITSE cost functions are presented in Table 6.2. Unit step responses of the closed-loop systems with PI control, corresponding to the optimal solutions presented in Table 6.2 are given in figures 6.2, 6.2, 6.2, and 6.2. Analysis of these step responses shows that the type of response does not depend much on the selection of the criterion as far as the constraint (19) is included in the optimization. Therefore, any criterion of the considered four would result in approximately the same performance. Therefore, the use of the gain margin constraint imbedded into the MRFT algorithm provides a significant equalizing effect with respect to the criterion selection effect. However, we select the ISE criterion as providing slightly higher values of c2 than other criteria, which better corresponds to the current practice of flow loops tuning, and further proceed with the use of this criterion. Having selected the integral criterion we can now formulate the problem of finding optimal coefficients that define optimal tuning rules in non-parametric 20
Fig. 4. Step response of IAE-optimal with constraint γm = 2 PI controller; DSOPDT model parameters: T = 1, ξ = 0.8, τ = 1.5, equivalent τ /Te = 0.32 of FOPDT model
Fig. 5. Step response of ITAE-optimal with constraint γm = 2 PI controller; DSOPDT model parameters: T = 1, ξ = 0.8, τ = 1.5, equivalent τ /Te = 0.35 of FOPDT model
tuning with MRFT as follows: minimize g(c1 , c2 , c3 ) = max
τ /T ∈D
subject to γm c1
s
(
)
f (c1 , c2 , c3 , τ /T ) , f ∗ (τ /T )
1 1 + 2πc3 − 2πc2
2
= 1,
(32)
(33)
where c1 , c2 , c3 are decision variables, τ /T ∈ D is a situational parameter, D := {τ /T : τ /T ∈ [1.3; 2.3]} is the domain of variation of the situational parameter, f (c1 , c2 , c3 , τ /T ) =
Z∞
e2 (t, c1 , c2 , c3 , τ /T )dt
0
21
(34)
Fig. 6. Step response of ISE-optimal with constraint γm = 2 PI controller; DSOPDT model parameters: T = 1, ξ = 0.8, τ = 1.5, equivalent τ /Te = 0.29 of FOPDT model
Fig. 7. Step response of ITSE-optimal with constraint γm = 2 PI controller; DSOPDT model parameters: T = 1, ξ = 0.8, τ = 1.5, equivalent τ /Te = 0.32 of FOPDT model
is a parametrized ISE cost function, g(c1 , c2 , c3 ) is the ISE cost function on domain D of situational parameters, f ∗ (τ /T ) is the solution of the minimization problem of f (c1 , c2 , c3 , τ /T ) subject to constraints (33), for a given τ /T , γm is a specified gain margin of the loop. The MRFT parameter β is selected in accordance with γm as per (20). Error e(t) is measured in the loop reaction to the set point step change. In practical solution of this problem we limit it only to the design of optimal tuning rules for a PI controller, which agreed the practice of process control (the derivative component is not used because of high noise component in flow measurements). At first the conventional optimization problem is solved separately for each point of τ /T . Functions f ∗ (τ /T ) are presented in Fig. 6.2 for several different γm and used as weights in the criterion (32). The use of weight f ∗ (τ1/T ) allows one to provide 22
the possibility of comparison of optimal values in different points of τ /T . The whole criterion (32) thus can be interpreted as providing minimal possible deterioration of optimality on the domain D due to the use of non-optimal solution (tuning rules) in most points (except for possibly one) of the domain D. It is worth noting that the criterion (32) can also be used if the domain D is two- or higher-dimensional. In that case certain “meshing” of the domain of situational parameters needs to be done to investigate function f in various points of D. The results of the solution of the optimization problem in terms of the cost function value (34) are presented in Fig. 6.2. These values of the cost functions can now be used as reciprocal weights for the optimization on the domain D in accordance with the criterion (32). The results of solution of the optimization problem (32), (33) for each gain margin value, in terms of coefficients c1 and c2 are presented in Tables 6.2, 6.2, 6.2, and 6.2. It is worth noting that the performance deterioration on the domain compared to the performance in a particular point is insignificant. For all for considered cases it is smaller than 1% in terms of the cost function increase in any point of the domain D. The produced tuning rules well agree with the practice of flow loop tuning [20].
Fig. 8. ISE-optimal cost function values for DSOPDT model and γm ∈ [2; 5]
With results of the optimization available, the controller tuning can be described as the following step-by-step algorithm: A) Desired gain margin γm is selected. 23
Table 5 Optimal tuning rules for set point response and gain margin γm = 2 Controller
c1
c2
c3
β
P
0.500
0
0
0
PI
0.451
0.331
0
0.548
Table 6 Optimal tuning rules for set point response and gain margin γm = 3 Controller
c1
c2
c3
β
P
0.333
0
0
0
PI
0.296
0.308
0
0.603
Table 7 Optimal tuning rules for set point response and gain margin γm = 4 Controller
c1
c2
c3
β
P
0.250
0
0
0
PI
0.222
0.307
0
0.607
Table 8 Optimal tuning rules for set point response and gain margin γm = 5 Controller
c1
c2
c3
β
P
0.200
0
0
0
PI
0.177
0.299
0
0.628
B) The modified RFT with parameter β corresponding to the selection made at step A is carried out. C) The values of frequency Ω0 and amplitude a0 of the self-excited oscillations in the system are measured. D) Tuning parameters of the controller are calculated per (15).
6.3
Tuning for disturbance rejection
For the optimization based on the consideration of a disturbance effect we have to consider the nature of the disturbances in the flow loop. As shown below, a mere application of an external signal to any point of the flow loop does not correspond to the actual situation. In analysis of equation (21) we see that, with valve opening and specific gravity constant, the disturbance may only come as a change of the pressure drop across the valve due to the change of the upstream (source) pressure or downstream pressure. Both these 24
changes act in the same way: increase of the upstream pressure or decrease of the downstream pressure would cause the increase of the pressure drop across the valve, which in turn would result in the flow increase, so that the controller would have to make adjustments to decrease the flow and bring it to the set point. And vice versa, decrease of the upstream or increase of the downstream pressures result in the decrease of the pressure drop and require the controller action aimed at the flow increase. To consider this problem in more details, we use the following diagram Fig. 6.3. In this setup the flow loop is supposed to control the flow in the line, which is done through the use of the flow transmitter, flow controller and the valve. At the same time operation of the second valve creates disturbances affecting the considered flow loop. Assume that at some point a steady state is maintained and the first valve is open at 30%, and the second valve is fully closed. We also assume that the valves are identical, and that the downstream pressures for both valves are constant (due to opening to the atmosphere, for example). The flow and pressure in this steady state can be found as shown in Fig. 6.3. The operation point providing steady state flow and pressure is found as the point of intersection of the pump performance curve 1 and the system curve 2, for 30% opening of the first valve and fully closed second valve. Therefore, the flow through the valve at this time is q1 . Assume that the second valve opens to exactly the same 30% at some point. The new steady state can be found through finding the point of intersection of the pump performance curve 1 and the new system curve, which now corresponds to both valves open at 30%. Because both valves are identical and have the same opening the total flow is calculated as double of that for one valve and the system curve 3 is made of the curve 2 by expanding it in the horizontal direction by two times. The total flow q3 is found as the point of intersection of curves 1 and 3. Yet the flow through the first valve is only half of the total flow and denoted as q2 . Therefore, we can see from 6.3 that opening of the second valve creates a disturbance to the flow control loop being considered in the form of pressure change (pressure drop change), and that if the opening of the second valve happens fast (so that the flow loop does not have enough time to adjust the position of the first valve) then immediately after coming of the disturbance the flow through the valve will be lower. The flow controller will have to make adjustments in the direction of valve opening to bring the flow to the set point. One can see that the disturbance to the flow loop has not additive but multiplicative character (see formula (21)). In fact, the increase of the upstream pressure increases the loop gain and the decrease of the pressure decreases the gain. This important circumstance leads to the formulation of a criterion of optimization for finding coefficients c1 , c2 and c3 . Therefore, tuning should be done in the process conditions providing maximal pressure drop across the flow valve. Alternatively, if such conditions cannot be produced at the time of loop tuning then lower gain margins should be selected thus compensating for possible loop gain increase due to the increase of pressure drop across the 25
Fig. 9. Flow loop disturbed by valve connected to the same source
Fig. 10. Pressure and flow in the system with two valves and a pump
valve. Despite the nonlinear character of the disturbance involvement a linear still can be used if we assume that the pressure drop across the valve changes stepwise. Assume also that this change occurs at time t = 0. We shall refer to the time immediately after this change as to t = 0+ and to the time immediately before the change as to t = 0−. Obviously, because the pressure drop change is step-wise, in the dynamic model used for optimization of the tuning rules we can consider a linear model of the process with parameters corresponding to the new operating point (new steady state) Assume that the initial condition of the system at time t = 0+ is described by variables u0 , l0 , q0 (we disregard for a moment other variables implicitly containing in the process transfer function). The initial conditions do not match the ones of the new steady state (at t = 0+). Therefore, the initial 26
conditions of the deviations are Δq(0+) 6= 0, Δl(0+) 6= 0, and Δu(0+) 6= 0. There is, however, a relationships between the initial values of Δl(0+) and Δu(0+) through the static gain (due to the assumption that the system was in a steady state before application of the disturbance), and between Δu(0+) and Δq(0+) through the coefficient. We should also note that the value of the set point r did not change.qNow using the substitution u0 = u − u(0+), q l0 = l − l(0+),q 0 = q − l(0+)Cv δpgs0v , r0 = r − l(0+)Cv δpgs0v , we obtain a new equivalent system in deviations Δq 0 (0+) = 0, Δl0 (0+) = 0, Δu0 (0+) = 0, and q Δr0 (0+) = −l(0+)Cv δpgs0v . From this analysis, one can see that application of a step-wise multiplicative disturbance in the form of the pressure drop change can be equivalently analyzed as a step change in the set point using the linear model, subject to the pressure drop across the valve being considered the pressure drop in the new operating point. Therefore, the optimal tuning rules obtained above for the step point change (Tables 6.2, 6.2, 6.2, and 6.2 are fully applicable to the case of tuning for the best disturbance rejection. The approach to the system transformation considered above is valid, of course, if we assume that the dependence of flow on the valve position and on the pressure drop is purely static, which is an accurate assumption when the fluid is incompressible. For compressible fluids, there is dynamic dependence of the flow on the valve opening and pressure drop, so that optimal tuning rules aimed at the best disturbance rejection should differ from those aimed at the set point following. However, considering the equalizing effect of the MRFT revealed as very insignificant deterioration of the optimality over the domain compared to the optimality in a point, one could legitimately expect that optimal tuning rules for the disturbance tuning should not be much different from those for the set point following - even if a compressibility of the fluid in the process model is accounted for.
7
Conclusions
An approach to producing optimal process-specific tuning rules for a PID controller in a flow loop, based on the modified relay feedback test, is presented. The approach is based on the developed solution of the optimization problem involving uncertainty due to variety of possible models. It is also based on the development of a suitable model of the flow process, modes of operation and the mechanism of disturbance generation in the loop. As a result of this approach, a set of tuning rules for different specifications for the gain margin is produced. The produced tuning rules well agree with the practice of flow loop tuning. 27
References [1]
Ziegler J.G. and N.B. Nichols (1942) Optimum settings for automatic controllers. Trans. Amer. Soc. Mech. Eng., 64: 759–768.
[2]
Astrom, K. J. and T. Hagglund (1984) Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatica, 20: 645– 651.
[3]
Cohen, G.H. and G.A. Coon (1953) Theoretical considerations of retarded control, Trans. ASME, 75, 827–835.
[4]
Boiko I. (2012) Loop tuning with specification on gain and phase margins via modified second-order sliding mode control algorithm, Internat. J. System Science, Vol. 43 (1), 97–104.
[5]
Luyben W.L., and Luyben M.L. (1997) Essentials of Process Control, McGrawHill, New York.
[6]
PID Control: New Identification and Design Methods (2005). M. Johnson and M. Moradi (eds.). Springer-Verlag, London, UK.
[7]
O’Dwyer A. (2006) Handbook of PI and PID Controller Tuning Rules, London, Imperial College Press.
[8]
Visioli, A. (2006). Practical PID Control. Springer-Verlag, London, UK.
[9]
Kyparisis J. (1985) Sensitivity analysis in nonlinear programming under second order assumptions, Lecture Notes in Control and Information Science, no. 66 (A.V. Balakrishnan and M. Thoma, eds.), Springer-Verlag, New York, 74–97.
[10] Levy, A.B. and R.T. Rockafellar (1994). Sensitivity analysis of solutions to generalized equations. Trans. AMS, Vol. 345 (2), 661–671. [11] Shapiro, J.F. (1979) Mathematical programming: Structures and algorithms. New York: Wiley-Interscience [John Wiley & Sons]. 388 p. [12] Berge, C. (1963) Topological Spaces. Oliver and Boyd, Edinburgh and London (reprinted by Dover Publications, Inc., Mineola, New York 1997). [13] Falk, J.E. (1976) Exact solutions of inexact linear programs. Operations Research, 24, 783–787. [14] Ben-Tal, A. and A. Nemirovski (1998). Robust convex optimization. Math. Oper. Res., 23, 769–805. [15] El Ghaoui, L. and H. Lebret (1997) Robust solutions to least-squares problems with uncertain data. SIAM J. Matrix Analysis and Applications, 18 (4), 1035– 1064. [16] Rockafellar, R. T. & R. J.-B.Wets (1991). Scenarios and policy aggregation in optimization under uncertainty. Mathematics of Operations Research, 16 (1), 119147.
28
[17] Diwekar, U. (2008). Optimization under uncertainty. Springer Optimization and its Applications, Vol. 22, 1–54. [18] Sahinidis, N. V. (2004). Optimization under uncertainty: State-of-the-art and opportunities. Computers and Chemical Engineering, 28, 971–983. [19] Boiko I. (2008) Modified relay feedback test and its Use for Non-parametric Loop Tuning. 2008 American Control Conference, Seattle, USA, 4755–4760. [20] Boiko I., Ernyes, A., Oli, W. and E. Tamayo (2009) Modified relay feedback test: industrial loop tuner implementation and experiments, 2009 American Control Conference, Saint Louis, USA, 4693–4698. [21] Bartolini G., Ferrara A. and E. Usai (1997) Output tracking control of uncertain nonlinear second-order systems. Automatica, 33 (12), 2203–2212. [22] Atherton D.P. (1975). Nonlinear Control Engineering - Describing Function Analysis and Design, Workingham, Berks, UK: Van Nostrand Company Limited. [23] Boiko I. and S. Sayedain (2010) Analysis of dynamic nonlinearity of flow control loop through modified relay feedback test probing, Internat. J. Control, Vol. 83 (12), 2580–2587.
29
