Relay Identiﬁcation Using Shifting Method for PID Controller Tuning

: The aim of this study was to present a relay shifting method for relay feedback identiﬁcation of dynamical systems suitable for PID controller tuning. The proposed technique uses a biased relay to determine frequency response points from a single experiment without any assumptions about a model transfer function. The method is applicable for open-loop stable, unstable, and integration processes, even with a delay, and regardless of whether they are oscillating or non-oscillating. The core of this technique was formed by the so-called relay shifting ﬁlter. In this study, the method was applied to a parameter estimation of a second-order time-delayed (SOTD) model that can describe, with acceptable accuracy, the dynamics of most processes (even with a transport delay) near the operating point. Simultaneously, a parameter setting for the PID controller was derived based on the model parameters. The applicability of the proposed method was demonstrated on various simulated processes and tested on real laboratory apparatuses.


Introduction
In the process industry, relay feedback identification is popularly used for estimating model parameters. The relay control was originally used for process identification by Rotač [1]. This approach has been increasingly explored in recent years due to the relay feedback method suggested by Åström and Hägglund [2] for PID controller tuning. In the relay feedback test, the period of oscillation and the amplitude of the process output were measured. The critical gain and the critical frequency were determined by the describing function. For more details about the describing function, see [3]. The same parameters can also be estimated using the continuous cycling method by Ziegler and Nichols [4]. However, with the relay feedback method proposed by Åström and Hägglund, one can obtain these parameters without a priori information about the process in a controlled manner and in a shorter time. The critical frequency and the critical gain can be utilized for the estimation of one point on the process frequency response. However, the knowledge of one frequency point only allows for the estimation of two model parameters. Therefore, models with a relatively small number of parameters are predominantly used for modeling. Most often, these are first-or second-order time-delayed linear models. To obtain the model parameters of these simple models from the relay feedback test, Luyben [5] assumed that the steady-state gain was already known and estimated the dead time from the initial response. Shen et al. [6] suggested using a biased relay for the steady-state gain estimation. Li et al. [7] proposed using two relay feedback tests, wherein the second test runs with an additional transport delay. Be et al. showed that multiple frequency response points can be estimated from a single relay feedback test using a parasitic relay and the FFT algorithm [8].
In [9], Chidambaram and Sathe describe another possibility, in which the exponential attenuation factor of the expression is used in the Laplace transform. At present, more publications exist that are devoted to relay-based parametric estimation techniques for control purposes. A summary of the obtained results can be found in survey publications, for example, [9][10][11][12][13][14].
The existing relay identification methods often assume that the processes are at a steady state, use several relay feedback tests, do not consider disturbances, and pre-suppose knowledge about the transport delay or the static gain. Most methods are further associated with the type of process (underdamped/overdamped, stable/integrating/unstable). In the work by Sánchéz et al. [15], it is shown that there are only two methods for measuring several frequency response points in only one relay experiment. The first one is based on the Laplace transform of input and output, for example, [9,16], and the second method is the shifting method described in [17,18]. The extension of the basic shifting method and its use for overdamped systems is described in [15].
The aim of this paper is to present the relay shifting method for identifying various types of dynamical systems and for tuning PID controllers. First, the relay shifting method was derived, which helps to determine multiple frequency response points from a single experiment without any assumptions about the model transfer function. Then the proposed method was applied to fit the SOTD model to a wide range of processes with various properties (underdamped/overdamped, stable/integrating/unstable, and time-delayed systems). A derivation of PID controller parameters for the control of the stable process described by the SOTD model is given in Section 6. Sections 5 and 7 describe the experimental results. Finally, the essential properties of the shifting method are summarized.

Relay Shifting Method Background
Let us consider a time-invariant process that can be described by the frequency response function G P (jω) around its operating point, where j is the imaginary unit and ω is the angular frequency. The process is kept near the operating point by a two-position asymmetric relay; see Figure 1, where y denotes the controlled variable, w is the desired variable, e is the control error and u is the manipulated variable. Figure 2 shows the static characteristics of the relay with hysteresis, where the parameters u A , u B are selected so that there is a stable oscillation with the period T p (T p = T 1 + T 2 , T 1 = T 2 ) (see Figure 3). The parameters ε A , ε B help to reduce the influence of a noisy environment. The hysteresis helps reduce the effects of measurement noise by enabling noiseless control action and increases the period of oscillation. It is recommended to adjust the hysteresis to 2-3 times the noise level [19]. The experimental time corresponds to approximately 3 to 4 periods of T p [13].   The core of the method consists of a simple calculation of the auxiliary variables u n (t) and y n (t), n = 2, 3, . . . , N from the measured courses y 1 (t) ≡ y(t) and u 1 (t) ≡ u(t) using the following relations: The calculated variables y n (t) and u n (t), n = 2, 3, . . . , N are periodic with the period T P /n if the measured variables u 1 (t) and y 1 (t) are periodic with the period T P . In our case, we assume t > t L , where t L is the time when the stable oscillation is achieved. In Figure 4, the time courses of the measured variables u 1 (t) and y 1 (t) are depicted together with the calculated variables y n (t) and u n (t), n = 2, 3 for illustration. The relationships (1) and (2) describe the shifting filter with the following frequency transfer function: F n (jω) = n for ω = n · 2π/T P 0 for ω = k · 2π/T P , k = 1, 2, · · · , n − 1 The shifting filter F n filters out lower harmonics and amplifies the n-th harmonic. The frequency response points G(jω k ), k = 1, 2, . . . , N can be determined by one of the following procedures ((a) or (b)): (a) where The integrals in the numerators are calculated by numerical integration. The great advantage of the relay shifting method is that without assuming any model structure, it is possible to estimate the frequency response points G(jω k ), k = 1, 2, . . . , N from a single relay test.
(b) The position G(jω k ) can also be determined without the utilization of Relation (5) using the approximation of the periodic signals u k (t) and y k (t) for k = 1, 2, . . . , N, t > t L by harmonic functions and by determining their amplitude ratio and phase shift (see [17] for example).
The real and imaginary parts (Re[G(jω k )] and Im[G(jω k )]) or the amplitude ratio |G(jω k )| and the phase shift arg(G(jω k )) determine two conditions for each frequency response point G(jω k ). This gives us 2N conditions for fitting a model from N frequency response points G(jω k ), k = 1, 2, . . . , N.
If we know the static gain K of a proportional system, then the frequency response point G(0) is determined too. This point can also be estimated from a single relay test by the following formula (see [6,20]) if the operating point (u 0 ,y 0 ) is known since the shifting method uses the biased relay.

Model Fitting
Let us assume that a process can be described by a model transfer function, as follows: where s is the complex variable in the Laplace transform, a 0 , a 1 , . . . , a n , b 0 , b 1 , . . . , b m ∈ R, and T d is a transport delay. Then the parameters a 0 , a 1 , . . . , a n , b 0 , b 1 , . . . , b m , and T d can be estimated using the frequency response points G(jω k ), k = 1, 2, . . . , N by iterative calculation if the selected value of N is ≥(n + m + 3)/2. If the process is at a steady state, we can obtain the time delay from the initial output response to the relay test and determine the other parameters by the least-squares method (LSM), for example, [21]. The LSM can also be used for parameter fitting if a rational (Padé) approximation of the transport delay term is substituted. The previous theoretical result is not directly usable for practical applications because we must consider the initial conditions, the precision of measurement, the position of the frequency response points, the choice of the linear model structure, the influence of disturbances, numerical aspects, and so on.
Our goal was to find a mathematical model that will be used to design a PID controller. For this purpose, we used a second-order time-delayed (SOTD) model with a transfer function.
The SOTD model is very versatile and can be used to describe most non-oscillating or oscillating proportional systems with or without a transport delay (see [22,23]). This model has only four parameters-a 2 , a 1 , K, and T d , which can be estimated from the frequency response points G(0) and G(jω k ), k = 1, 2, . . . , N. The model parameters can be estimated bŷ where the vector of model parameters Θ is the set of acceptable values of the model parameters. For a stable system, τ m is the maximum possible transport delay (see Figure 3) and "T" denotes the transpose operator. By denoting the real and imaginary part of the complex values G(jω i ), Then,θ = argmin where ∆τ is the selected accuracy of the estimation τ, the matrix Z is assumed with the maximum column rank. In this way, we can estimate the model parameters from 2N equations (generalization in comparison with [24]). A similar approach can also be used for fitting linear models with more parameters. The proposed procedure, unlike that of [15], does not require the solution of a system of nonlinear equations and is also applicable to the identification of underdamped systems. Model (15) has four parameters and, therefore, only two frequency response points are sufficient for the parameter estimation (N = 2). In case we know the values of some parameters, for example, the steady-state gain or the dead time, the number of estimated parameters will be reduced, and a similar approach may be used.

Relay Shifting Method with Additional Delay or Integrator
Using the shifting method, we obtained the frequency response points for higher frequencies. We know from classical control theory that the values of the process transfer function for PID tuning are decisive in the range where the phase delay is from 0 • to 180 • and, consequently, for lower frequencies. To get the frequency response points for lower frequencies, we used an additional integrator or a delay (see Figure 5). The additional integrator or the delay can be used only if a steady-state oscillation can be achieved with the relay feedback control. To a certain extent, the frequencies of the obtained frequency response points can also be reduced by increasing the relay hysteresis.

Examples of Using the Relay Shifting Method
In this section, several different types of processes are identified using the relay shifting method. The proportional processes are modeled by transfer function (9) and the integrating process is modelled by transfer function (16).
In all the simulated examples, we chose N = 2 for the parameter estimation and the biased relay had the following parameters (see Figure 2).

Lag-Dominated Non-Oscillatory Proportional Process
The lag-dominated non-oscillatory proportional process considered in [22] is given by the following transfer function: The process is controlled by a relay (see Figure 1). The relay output u and the process output y are depicted in Figure 6.
The Nyquist frequency responses of process P 1 (s) and model M 1 (s) are shown in Figure 7. In the same figure, model M S (s) and the points G(0), G(jω 1 ), and G(jω 2 ) are depicted as well, where The model M S (s) was derived by Skogestad [22] directly from (18) using the half rule. The unit step responses of process P 1 (s), model M 1 (s), and model M S (s) are shown in Figure 8.

Lag-Dominated Oscillatory Proportional Process
The lag-dominated oscillatory proportional process considered in [23] is given by the following transfer function: where ξ 1 = 0.5, ω n1 = 1.1180, ξ 2 = 0.6, ω n2 = 2.2361, ξ 3 = 0.65, ω n3 = 3.1623, ξ 4 = 0.5, and ω n4 = 1.4142. The process is controlled by a relay (see Figure 1). The relay output u and the process output y are depicted in Figure 9. Using the relay shifting method, we can determine from the observed stable courses u and y that T P = 6.44s, and the SOTD model, estimated according to Sections 2 and 3, is The Nyquist frequency responses of process P 2 (s) and model M 2 (s) are shown in Figure 10. The model M R (s) was derived by Novella-Rodriguez et al. [23] directly from (25). If we use an additional delay or an integrator (see Figure 5), we can receive better positions of the points G(jω 1 ) and G(jω 2 ) for the model fitting. In case we use the additional delay D = 3 s, we obtain the time courses u and y depicted in Figure 11. Using the relay shifting method, we can determine from the observed stable courses u and y that T P = 11.3s, and the SOTD model, estimated according to Sections 2 and 3, is The Nyquist frequency responses of process P 2 (s), model M 2D (s), model M R (s), and points G(jω 1 ) and G(jω 2 ) are shown in Figure 12. The unit step responses of process P 2 (s) and models M 2 (s), M 2D (s), and M R (s) are shown in Figure 13.

Delay-Dominated Non-Oscillatory Proportional Process
The delay-dominated non-oscillatory proportional process considered in [11] (p. 165) is given by the following transfer function: The process is controlled by a relay (see Figure 1). The relay output u and the process output y are depicted in Figure 14. Using the relay shifting method, we can determine from the observed stable courses u and y that T P = 10.90s, and the SOTD model, estimated according to Sections 2 and 3, is The Nyquist frequency responses of process P 3 (s) and model M 3 (s) are shown in Figure 15. In the same figure, the model M LG (s) and the points G(jω 1 ) and G(jω 2 ) are depicted as well, where M LG (s) = 1.000 · exp(−4.2771s) (0.6183s + 1) 2 (43) The model M LG (s) was derived by Liu and Gao [11] (p. 165) with algorithm RS-SU1. The unit step responses of process P 3 (s), model M 3 (s), and model M LG (s) are shown in Figure 16.

Integrating Process
The integrating process considered by Liu and Gao [11] (p. 212) is given by the following transfer function: The process is controlled by a relay (see Figure 1). The relay output u and the process output y are depicted in Figure 17. Using the relay shifting method, we can determine from the observed stable courses u and y that T P = 50.1s, (45) The model M 4LGI (s) was derived by Liu and Gao [11] (p. 212) with algorithm RI-S1. The unit step responses of process P 4 (s), model M 4 (s), and model M 4LGI (s) are shown in Figure 19.

Unstable Process with Delay
The same approach can be used for unstable processes if there is sustained oscillation in an asymmetric relay feedback test. The unstable process with delay considered by Liu and Gao [11] (p. 235) is given by the following transfer function: This unstable process is controlled by the biased relay where (see Figure 2). The relay output u and the process output y are depicted in Figure 20. From these time courses, we can determine by the shifting method, The estimated model is The Nyquist frequency characteristics of process P 5 (s) and model M 5 (s) are shown in Figure 21 The model M 5LG (s) was derived in [11] (p. 235) with algorithms RU-SA and RU-SB, which result in very similar models.
The unit step responses of process P 5 (s), model M 5 (s), and model M 5LG (s) are shown in Figure 22.

PID Controller Tuning
If we know the process model, we need to determine the parameters of the controller. There are a number of recommendations for PID controller tuning that can be used for control of a process described by transfer function (9) (see, for example, Åström and Hägglund [25], Skogestad and Grimholt [26], Novella-Rodriguex et al. [23], and O'Dwyer [27]).
We can also use the following tuning rules, which are very simple, and they should work well on a wide range of stable processes (non-oscillatory/oscillatory, lag/delay dominated) described by process model (9).
Let the transfer function of the ideal PID controller be where r P ,r I , and r D are adjustable controller parameters. The open-loop transfer function is and the closed-loop transfer function is The parameter r I can be determined by the Nyquist stability criterion using the openloop transfer function (60) and by a choice of the gain margin m A or the phase margin γ. For common control systems, the recommended values are as follows (see, for example, Víteček and Vítečková [28], and Landau and Zito [29]): m A ∈ 2, 5 , γ ∈ π/6, π/3 .
With this approach, we can easily derive that and The value of the parameter r I can also be chosen according to the denominator Q(s) of the control transfer function (61). The characteristic equation of the control system has the transcendental form The control system is stable for This equation also follows from (63) and (64) for m A = 1 and γ = 0. With the parameter r I , we can change the behavior of the closed-loop system. The value of the parameter r I can be chosen, for example, from the characteristic Equation (65) to obtain the dou-ble real dominant pole s 1,2 . This pole can be determined from Equation (65) and the following equation: dQ(s) From (65) and (67), it follows and By comparing, for example, (63) and (69), it is obvious that the adjustment according to (69) corresponds to the gain margin m A = 4.27.
To summarize the presented results for PID tuning so far, at first, we determine the parameter r I (using one of the relations (63) or (64)), and then we determine using relation (59) The resulting relationships for adjusting the PID controller are computationally simple and easily understandable; they are not based on the approximation of the transport delay and are usable for stable oscillatory or non-oscillatory systems.

Examples of PID Tuning to Control Stable Processes Identified by the Relay Shifting Method
In all examples, the parameter r I was estimated according to relation (69), and the parameters r P , r D by relations (70), (71). The PID control was applied directly to the original processes, and the models, obtained by the relay shifting method, were used only for the PID tuning. The processes were controlled by the PID 2DOF controller which is described by where U(s), E(s), and Y(s) are the Laplace images of u(t), e(t), and y(t), respectively. The PID 2DOF controller eliminates "derivative kick" and is less sensitive to measurement noise. In all examples, the parameter N = 10. Let us consider the closed-loop system given in Figure 23. The simulated unit step responses of closed-loop control systems for setpoints w and change of load disturbances d are shown in the following examples.

PID Control of Process P 1
The model M 1 (23) of the process P 1 (18) was obtained by the relay shifting method. The relations (69), (70), and (71) determine the parameters r I , r P , and r D of the PID controller, respectively.
The unit step responses y of the closed-loop control system for the setpoint change w and the change of load disturbance d are shown in Figure 24. The unit step of the load disturbance d was at the time t = 20 s.

PID Control of Process P 2
The model M 2 (30) of the process P 2 (25) was obtained by the relay shifting method. The relations (69)-(71) determine the parameters r I , r P , and r D of the PID controller, respectively.
The controlled variables y for the setpoint change w and the change of load disturbance d are shown in Figure 25. The unit step of the load disturbance d was at the time t = 40 s.

PID Control of Process P 3
The model M 3 (42) of the process P 3 (37) was obtained by the relay shifting method. The relations (69)-(71) determine the parameters r I , r P , and r D of the PID controller.
The unit step responses y of the closed-loop control system for the setpoint change w and the change of load disturbance d are shown in Figure 26. The unit step of the load disturbance d was at the time t = 40 s.

Verification and Implementation
The proposed identification method for PID tuning was tested on both simulation examples and on real laboratory apparatuses, namely "Water Levitation", "Air Aggregate", and "Air Levitation" (see [30,31]), and the results of the relay shifting method were compared with the results of other identification approaches to verify the practical benefits of the use of this method in the work by Hornychová and Hofreiter [32]. The experiments were realized both in the Matlab/Simulink environment and physically, using the industrial PLC Tecomat Foxtrot produced by TECO a.s. (see TECO Advanced Automation [33]). The comprehensive development tool Mosaic was used for writing and debugging the relay identification and PID tuning programs for the PLC Tecomat Foxtrot (see Vaněk and Slabý [34]).

Conclusions
The presented relay shifting method has the following properties:

•
The frequency points G(jω i ), i = 1, 2,· · · , N are determined from one relay feedback test without any prior knowledge of the model.

•
The relay shifting method is very versatile, the same procedure can be used for the model parameter estimation of stable/unstable/integrating, non-oscillatory/oscillatory, and lag/delay dominated processes and the noisy environment if there are stable oscillations for the relay feedback test.

•
The static load disturbance has no effect on the locations of the frequency points G(jω i ), i = 1, 2,· · · , N.
• By using the SOTD model, it is possible to estimate the steady-state gain even in the presence of a constant load disturbance. • The shifting filter described by (1) and (2) filters out the lower harmonics and amplifies the n-th harmonic, and the frequency response point G(jω n ) can be determined by relation (5) or by using the approximation of the periodic signals u n (t) and y n (t) by harmonic functions.

•
The relay shifting filter helps to eliminate parasitic signals.

•
The relay shifting method is appropriate only for processes describable by linear time-invariant models.

•
The position of the points G(jω i ), i = 1, 2, . . . , N significantly affects the estimation of the model parameters. An additional integrator or a delay helps to improve the parameter estimation. • Using the SOTD model, it is possible to model processes that can be described by high-order models.

•
The relay shifting method is primarily suggested for the automatic tuning of controllers.
The combination of the shifting method and the additional delay or the additional integrator in the relay feedback loop allows for the estimation of the frequency response points G(jω k ), k = 1, 2, . . . , N, from a single relay test. To obtain more frequency response points, we can select different values of the additional delay or the integral gain of the additional integrator for each relay feedback test. This approach makes it possible to acquire other frequency response points.
In this paper, we have summarized the relay shifting method that can be used for estimating the parameters of simpler and more complex linear models with different structures describing systems where steady oscillations can be achieved by means of relay control. We presented a simple original algorithm for estimating the parameters of the SOTD model from N points of the frequency response. The relay shifting method for open-loop stable, unstable, and integration processes, even with a delay, and regardless of whether they are oscillating or non-oscillating, was demonstrated by simulation examples. The article also introduced the design and simulation testing of a simple PID controller adjustment for the control of a stable system described by an SOTD model without using the delay approximation.