Enhanced Event-Based Identification Procedure for Process Control

Apr 25, 2018 - Department of Computer Sciences and Automatic Control, UNED ..... by Bajarangbali and Mahji(23,33,34) report online and off-line method...
0 downloads 3 Views 3MB Size
Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Process Systems Engineering

Enhanced Event-based Identification Procedure for Process Control José Sánchez, María Guinaldo, Antonio Visioli, and Sebastián Dormido Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b05239 • Publication Date (Web): 25 Apr 2018 Downloaded from http://pubs.acs.org on April 26, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Enhanced Event-based Identification Procedure for Process Control José Sánchez†*, María Guinaldo†, Antonio Visioli‡, Sebastián Dormido† †

Department of Computer Sciences and Automatic Control, UNED, C/ Juan del Rosal 16, 28040 Madrid, Spain ‡

Department of Mechanical and Industrial Engineering, University of Brescia, via Branze 38, 25123 Brescia, Italy

KEYWORDS: Describing function, event-based, send-on-delta, estimation, limit cycle, oscillations.

ABSTRACT: An enhanced event-based identification procedure for process control has been developed based on the information obtained from the oscillations that an event-based sampler introduces in the feedback loop. The describing function analysis is used to explain the basis of the method as the event-based sampler behaves as a static nonlinearity. Features of the method are: (a) the event-based procedure does not require a priori process information, (b) non-iterative algorithms are sufficient to derive the process parameters, (c) only one test is needed, and (d) it allows identifying the process at a user-specified phase angle in the third quadrant. The method is presented for estimation of most common transfer functions found in chemical and process industry: integrators, first and second order, as well as processes with non-minimum-phase

ACS Paragon1Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dynamics. Further, the procedure can be extended to estimate models of any structure with five or more parameters.

1. INTRODUCTION One key issue in the chemical and process industry is to get some knowledge of the dynamics of the processes to be controlled. One way of doing so is by selecting a template of a transfer function that represents a specific dynamics (e.g., second order with time delay) and to estimate the parameters of this template using different approaches. The most common approaches1 are based on the application of a step to the process and the analysis of the output, or on the replacement of the controller by a relay to introduce oscillations and analyse the limit cycle. All these approaches are specifically designed to work in time-based control loops. To promote the inclusion of event-based control engineering techniques in the process industry as currently happens with their time-based counterparts, a key issue is to produce new methods by taking into account the intrinsic characteristics provided by the event-based paradigms. For example, the limit cycles can be considered an implicit product of the event-based feedback loop context. The motivation of this paper is to provide an identification method specifically designed to take advantage of the features provided by the event-based control. Thus, the described eventbased identification procedure has not been explicitly designed for standard time-based process control architectures, as it happens with the approaches based on relays presented in the literature. The main contribution of this work is a parameter estimation method based on an event-based sampler that demonstrates that the characteristics of the sampler (describing function, hysteresis) can be exploited in the design of a procedure to estimate a model. As said, the approach presented in this paper is event-based (see Miskowitcz2 for an extended treatment of the event-based paradigm and Blevins et al.3 as an example of the relevancy that

ACS Paragon2Plus Environment

Page 2 of 71

Page 3 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

this paradigm has in process industry). This means that the control loop has an event-based sampler receiving data from the process in a synchronized way but sending them to the controller only when a certain triggering rule in the sampler is satisfied. In most event-based samplers, the transmission from the sensor to the controller occurs when the process output varies of a certain amount, δ, with respect to the last sample sent to the controller. Due to its nonlinear nature, the sampler can introduce oscillations in the loop, and this oscillatory behaviour can be explained in the frequency domain thanks to the describing function (DF) of the sampler. However, the data needed to make the estimations are provided in the temporal domain, that is, from experimental measurements taken during a test (oscillation frequency, harmonics, oscillatory gain). For this reason, the event-based procedure uses tools from both domains: from the time domain, by measuring the oscillatory signals, and from the frequency domain, by using the describing function as the mechanism to pose and solve the estimation problem. The event-based feedback control loop considered in this work is divided into three units (see Figure 1): • The event-generator unit is composed of the sensor and the event-based sampler. The purpose of this unit is to measure the process output y(t), to calculate the error e(t) between the measured signal y(t) and the set-point value yr, and to send the error to the controller. The event-based sampler is in charge of taking the error e(t) and producing an event-based sampled signal e∗(t). • The control unit C(s) receives the sampled signal e∗(t) from the event-based sampler. In this work, and without loss of generality, the control law considered is a PID controller. As the estimation approach only requires applying a proportional action to produce oscillations

ACS Paragon3Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(known as oscillatory gain K osc ), some parts of the controller could be by-passed or deactivated during the estimation phase. • The actuator unit receives u(t) from C(s) and applies it to the process G(s). The autotuning logic block in Figure 1 contains the programming code to apply the different steps of the approach explained in the present paper. During the identification phase, the logic would be in charge of deactivating the integral and derivative actions of the PID controller in Figure 1. In addition, it would modify the proportional gain and would introduce changes in reference yr to disturb the process (for example, a step or a pulse) to produce the oscillations needed to take measurements and estimate the parameters. It is important to remark that the event-based control system working in normal conditions will not be oscillating and events should be generated just in the following situations: presence of disturbances, set-point changes, time-out for safety conditions, and identification phase. The control actions u(t) will be always event triggered, regardless the working conditions. However, during the identification phase, stable oscillations will be generated and it will be necessary to sample the system at regular times to calculate the integrals required to apply the approach here described; the control actions u(t) sent to G(s) will continue being triggered by events at regular time instants always at a lower frequency than the one required to sample the system. Once the model is identified and the controller tuned (for example, by AMIGO rules or other approaches) by the autotuning logic, the system will continue working in normal conditions, that is, without oscillating.

ACS Paragon4Plus Environment

Page 4 of 71

Page 5 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 1. Control scheme of the event-based control loop. The dashed line indicates that the data sent from the sampler to the controller are not synchronous.

For the majority of the architectures of event-based feedback control loops, the sampling is done by applying a Send-On-Delta (SOD) sampling strategy to the process output4,5. When the output is quantified by a quantity multiple integer of δ, the relationship between the input and the output of the event sampling block is symmetric with respect to the origin. This strategy is known as Symmetric-Send-On-Delta (SSOD) and has received much attention from the control community, and results on different control issues have been obtained during the last years6,7,8,9,10,11,12. The key point of the relationship between e(t) and e∗(t) created by the SSOD block is that it is similar to the behaviour of a relay with hysteresis where there is an infinite number of thresholds13. This implies that the describing function of the SSOD block can be applied to explain the event-based oscillations and used to solve the problem (at hand). The paper provides an event-based parameter estimation approach using a SSOD block located at the process output that can be applied to estimate transfer functions of any order and structure. A first approach on using the describing function of the SSOD block to the parameter estimation in an event-based control loop was reported in Beschi et al.14. In that work, the process parameters are estimated considering a limit cycle generated by the SSOD block and a pre-tuned PI controller, being necessary to know at least one of the three FOPTD parameters. In other works on the parameter estimation in an event-based control loop8,15, the estimation approach is

ACS Paragon5Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

based on curving fitting and state-space using a relay in series with the SSOD block and the pretuned PI controller. From measurements obtained from the oscillatory signals and using the controller parameter values, expressions to obtain the process parameters are given. In the approach presented here, a previous tuning of the controller is not necessary and thus the method can be applied at the beginning of the design phase. The main result presented in this paper is an enhanced estimation procedure based on the oscillations introduced by an eventbased sampler such as the SSOD. The proposed approach allows the identification of the model parameters from measurements derived from the oscillations generated by the event-based sampling. These measures are the oscillation frequency, the main odd harmonics, and the DC gain (if the process is without integration). The adjective enhanced is applied to the method for four reasons: 1. Unlike previous works based on DF, the model parameters do not need to be calculated using a priori information of the process (i.e., static gain, velocity gain, or dead time); 2. Iterative methods to solve non-linear equations are not necessary; 3. The number of tests performed to obtain the measures is always one, regardless of the transfer function structure or order; 4. The procedure can be extended to estimate any type of transfer function. Likewise, no matter the real process structure is, the event-based approach allows always to obtain a first-order-plus-time-delay (FOPTD) or a second-order-plus-time-delay (SOPTD) transfer function that fit to the real process behaviour at ω = 0 and at a user-specified phase lag in the third quadrant of the Nyquist plot. Other DF-based methods do not produce FOPTD solutions if the real process structure is not a first order one or

ACS Paragon6Plus Environment

Page 6 of 71

Page 7 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

produce a solution whose behaviour at the steady-state is different if the true process is of a higher order than the transfer function to fit. A study of a basic version of this approach has been presented in Sánchez et al.16. In such event-based approach, the transfer functions to fit are limited to three parameters (static gain, one lag and dead time), reducing the range of models to identify. Another key difference is that in Sánchez et al.16 the amplitude of the limit cycle generated by the SSOD sampler must be high to produce good solutions as the method depends on data provided by the SSOD DF (the highest the number of steps of the SSOD output, the better accuracy from the DF). In this improved version, the limit cycle amplitude can be any as the approach does not use the SSOD DF expression to produce results. It must be underlined the relationship of the multi-step signal that can be produced by a SSOD sampler with the approach described in Jeon et al.17 to produce oscillations. The authors improve the estimation of the critical frequency and gain by using a multi-step input signal generated by adding sub-relay signals of different amplitudes and frequencies; the rationale is to reduce the harmonics and obtain an input signal much more similar to a sinusoidal function. In the approach here described, it is not necessary to force the SSOD sampler to produce multi-steps signals. The structure of the paper is as follows. First, the basis of the estimation approach and the problems that have been found in the literature are described in Section 2. In Section 3, the procedure is explained and solutions to the problems found when the DF is used for estimation are given; moreover, expressions to estimate the parameters of common transfer functions are derived and simulation results are presented. The generalization of the method to transfer function models with five or more parameters is explained in Section 4. Section 5 discusses some

ACS Paragon7Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ways to estimate with user-defined phase lag specifications. The paper finishes with conclusions and recommendations about the use of the approach.

2. ESTIMATION WITH THE DESCRIBING FUNCTION: BASIS AND PROBLEMS The most similar methods for estimation of transfer function parameters to the event-based approach are the relay identification methods based on its describing function. The pioneering works on the use of relay feedback for identification purposes are from the 80’s18 and the work of the process control community on this subject continues until now19,20,21,22,23. The basis of the method is that a linear system under an ideal relay control will oscillate, approximately, at its ultimate frequency, that is, ω osc ≈ ω u and the critical gain Ku is derived from the describing function of the ideal relay. It is well known that the use of the describing function information can introduce errors in K180º

of 5%-20%24. For this reason, there has been a lot of research work to improve the

information that the describing function provides and to obtain additional critical points in the Nyquist plot at different phase lags to estimate transfer functions of order higher than one. The first attempt to use the describing function of a simple relay for estimating different types of transfer functions is the AutoTune Variation method (ATV)25. In that method, it is necessary to get a priori the steady gain and the dead time from inspecting the temporal response. Since the method uses only the ultimate gain and the ultimate frequency from the relay test, that is, just one critical point, achieving solutions for some models is difficult and, according to the ATV’s authors, “there is no guarantee that any of the models will fit the data”. The ATV was improved by the ATV2 method26 being only necessary to know the dead time. In the ATV2, analytical expressions were obtained for the steady-state gain as well as time constants of transfer function

ACS Paragon8Plus Environment

Page 8 of 71

Page 9 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

models from first to third order. The bi-ATV method has been introduced by Shen et al.27; the main difference with respect to the previous ones is the use of a biased-relay to obtain the steadystate gain from the experiment. The approach to solve the equations is similar to the original ATV. Another improvement, known as ATV+28, avoids the prior knowledge of the dead time as in previous ATV versions, the ATV+ proposes to find an estimate of the delay through the determination of minimum and maximum bounds. Another modification of the ATV2 method29 provides better results but just in limiting cases: for small and large values of ωˆ u T . The DF-based approaches presented by Srinivasan and Chidambaram30 and Vivek and Chidambaram31 provide a simple procedure to estimate a FOPTD system with a single symmetrical relay but it has not been extended to higher order processes and an iterative method has to be applied to obtain the system parameters. Another method, known as phase deviation32, obtains the parameters with only a relay test plus experimental measurements of the first-, thirdand fifth-order harmonics of the process. The approaches presented by Bajarangbali and Mahji23,33,34 report on-line and off-line methods to fit stable and unstable first and second order systems using the DF analysis of a relay with hysteresis to reduce the effect of measurement noise. In other DF based contributions35,36, the information provided by the DF formulation of a relay without hysteresis is used to know the critical point and the dead time is obtained by measuring the initial difference in the response of the process output with respect to the process input. One of the main disadvantages of the relay DF analysis is that it estimates one critical point in the Nyquist map, that is, G ( jω −180 º ) , limiting the number of transfer function parameters to identify. To obtain additional critical points with phase angles different to zero, the approaches found in the literature can be classified into two groups: (a) to displace the negative reciprocal of

ACS Paragon9Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the describing function of the relay (by modifying the hysteresis width) to intersect with the process at lower frequencies in the Nyquist plot, and (b) the opposite action, that is, to displace the lower frequencies of the process toward the negative reciprocal of the describing function. The modification of the hysteresis width is the first solution found in the literature to get additional points although it was originally proposed for tuning of simple regulators18. The negative reciprocal of the describing function of a relay with hysteresis can be described as a straight line parallel to the real axis in the Nyquist map. Thus, by manipulating the hysteresis width, this line can be moved up and down in the third quadrant of the Nyquist plot, allowing the intersection of the process at a desired phase angle36. There are three different solutions to displace the lower frequencies of the process toward the reciprocal of the relay DF. The first one is by inserting a known time delay26,28,37,38,39; as the delay as the effect of rolling the process Nyquist plot around the Nyquist zero, it is possible to get new critical points at lower frequencies than the original ω u . The second solution is to insert a nonlinear block composed of an integrator plus a relay in parallel with the original relay38,40,41,42; the phase lag of the integrator produces the orthogonal displacement of the lower frequencies. The third solution is to introduce filters in the loop43,44; by designing adequately the filter it is possible to introduce lag and arbitrarily displace a point of the loop transfer function to a desired point. Additionally, there is an extensive research on relay identification methods not based on the describing function that try to improve the accuracy. These methods are based on Typskin’s locus, curve fitting, Laplace transform, state-space analysis, and hybrid approaches. Estimation procedures based on the A-locus method, a derivation of Typskin’s locus, are described by Atherton and Kaya45,46; those works present analysis to estimate stable and unstable first and second order systems with delay and integrating processes using a relay with hysteresis. An

ACS Paragon10 Plus Environment

Page 10 of 71

Page 11 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

approach based on curve fitting for parameter estimation of integrating and time delay processes is reported in Panda et al.47. With the same rationale, the Liu and Gao’s book48 describes different algorithms to identify parameters of stable, unstable, and integrating processes of first and second order. Approaches based on Laplace transform are presented by Srinivasan and Chidambaram30 and Vivek and Chidambaram31; both contributions describe analogous procedures to estimate stable first order systems using a simple relay. Majhi49 reports expressions to identify stable and unstable first and second order systems with delay by a state-space approach using a simple relay. More contributions based on the statespace analysis to identify many process dynamics are reported by Bajarangbali et al.50, Bajarangbali and Majhi51, and Pandey and Majhi52. Regardless the rationale of the relay-based identification method, one common requirement is the convergence to a stable limit cycle oscillation, producing long experiment duration during which the process must not be subject to disturbances. Recent relay-based approaches to avoid such requirement are described in Soltesz et al.53 and Berner and Soltesz54. These approaches are posed as an optimization problem aiming to minimize the output error using all the experimentation data generated during a few asymmetric relay switches. A hybrid approach is presented in Berner et al.22 where a mixture of different techniques is applied. In such procedure, the type of model to estimate is determined by analysing the waveform of the asymmetric relay used to excite the process; if a first order dynamics with delay and w/o integration is needed, a curve fitting approach is applied; if a second order or higher, then an optimization procedure is run. It must be underlined the difference of the relay-based methods with other approaches where the system is excited by oscillations as, for instance, the one presented in Ionescu et al.55. The

ACS Paragon11 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

rationale of the method is to force the system to oscillate by using as input a multisine signal and multiply the process output by a sine and cosine, respectively, of the input signal frequency. The signals are generated in open loop by a sine/consine generator (a relay is not used) and the result is the response of the system for the entire range of frequencies using one excitation signal. In the DF relay-based methods, the system is identified using the cross-over frequency. In the event-based approach presented here, the frequencies used depend of the model order. After reviewing the literature on DF-based identification, the main problems found are summarized as follows: - The describing function provides just an approximation of the Nyquist point at the frequency where the process oscillates. This approximation reduces the accuracy of the results or even produces solutions in the complex plane. It is fundamental to obtain an exact value of G( jωosc ) during a test. - In processes without integration, it is necessary to obtain G(0) as a first parameter. Indeed, the identification will provide exact results at any frequency only if the transfer function template to fit is exactly equal to the real process to identify. On the contrary, if the true process has higher order or different structure from the template and G(0) is not known, this will produce the result to be accurate around the critical frequency

ω ≈−180 º but with discrepancies at ω = 0 (see Example 1). This is due to the fact that the template is fitted with lesser degrees of freedom than the true process. So, the fitting will be exact around ω osc but will present discrepancies in frequencies close to zero. - If the process has an integrator, as the identification is based on the behavior in the third quadrant, the results around ω ≈ −180 º will be accurate even when the transfer function has a different structure to the real process. However, if the structure is different, the results

ACS Paragon12 Plus Environment

Page 12 of 71

Page 13 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

will differ at low frequencies, improving the results if the identification is done at frequencies around ω−135º . - The classical DF-analysis only allows the derivation of two process parameters with the information of one experiment. If more parameters are needed, they must be measured from the process reaction curve or a second experiment must be performed. A solution is to devise a procedure to get in just one test as many values of G ( jω ) (for instance,

G ( jω osc ) , G ( j 3ω osc ) , G (5ω osc ) , etc.) as needed to obtain the model parameters. At least one point must be situated in the third quadrant of the Nyquist map, that is

G( jωosc ) , and another point must be G(0) when needed, depending on the transfer function template to adjust. - It is not possible to select a user-defined phase lag for the identification. For example, recommendations on the margin phase of the Nyquist point to use in the identification depending on the process features are given by Friman and Waller40. Åström and Hägglund56 recommend estimating the behavior of the process at ω −135º or ω −180º depending on whether PI or PID control are to be applied, respectively.

ACS Paragon13 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 71

Figure 2. Differences between true process and estimated transfer function when G(0) is not considered in the fitting.

Example 1: An example of not taking into account G(0) can be found in Wang et al.32. As the identification

is

done

using

a

Nyquist

point

close

to

ω −180º

and

the

process G ( s) = [(− s + 1) ( s + 1) 4 ] e −4 s is estimated using a second order model, the fitting is good enough around the Nyquist point but with a 24.6% error in the steady-state gain. The estimated transfer function is Gˆ ( s ) = [1 (2.779 s 2 + 2.282 s + 1.246)] e −6.706 s (see Figure 2).

ACS Paragon14 Plus Environment

Page 15 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 3. Relationship between e(t) and e∗(t) in a SSOD block.

3. THE EVENT-BASED IDENTIFICATION APPROACH The event-based sampling used in this work, known as SSOD, is a variant of the original SOD. The main difference is that the SSOD is centered in zero, so that the relationship between the input and output is symmetric with respect to the origin. The sampled output signal e∗(t) can only assume values multiple of δ, namely e∗(t) = iδ with i ∈ Z. The sampled signal e∗(t) changes its value to the upper quantization level when the input signal e(t) increases more than δ, or to the lower quantization level when e(t) decreases more than δ. The behaviour of a SSOD block can be mathematically described as

 (i + 1)δ if e(t ) ≥ (i + 1)δ and e * (t − ) = iδ  e * (t ) =  iδ if e(t ) ∈ [(i − 1)δ , (i + 1)δ ] and e * (t − ) = iδ (i − 1)δ if e(t ) ≤ (i − 1)δ and e * (t − ) = iδ 

ACS Paragon15 Plus Environment

(1)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 71

The event-based sampling of a signal using SSOD produces an oscillatory symmetric signal similar to that generated by a relay with hysteresis where there is an infinite number of δ thresholds (i.e., ±δ, ±2δ, ±3δ, …, ±nδ), as shown in Figure 3.

Figure 4. Nyquist plot of −1 N( A, δ ) . The describing function of a SSOD block is given by the following equation57,

m −1 δ δ 2δ  1 + 1 −  m  + 2∑ 1 −  k  N ( A, δ ) = πA  A  A  k =1  2

2

 2δ 2 − j 2 m πA  

(2)

where A is the amplitude of a sinusoidal input signal, and m =  A δ  . In Figure 4, the portrait of − 1 N ( A, δ ) is represented for A ∈ [δ , ∞ ) . Each intersection of the process G ( jω ) with an arc in Figure 4 will produce an oscillation of a different amplitude: intersections with the lower arc produce oscillations with A ∈ [δ ,2δ ) , the next A ∈ [2δ ,3δ ) , and so on. As the condition for the existence of limit cycles is given by

ACS Paragon16 Plus Environment

Page 17 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

G ( jω osc ) = −

1 N ( A, δ )

(3)

the application of a specific gain C ( s ) = K osc to the process that produces an oscillatory behaviour will mean that (3) is being satisfied. That is,

Gol ( jω osc ) = K osc G ( jω osc ) = −

1 N ( A, δ )

(4)

It must be noticed in Figure 4 that the first point of − 1 N ( A, δ ) in the Nyquist plot, the C1, corresponds to an oscillation of amplitude δ and phase angle ϕosc = 45º . So, by varying K osc it will be possible to increase the amplitude and also to reduce the phase angle from 45º to near 0º. This point will be discussed in more detail in Section 5. As said before, the reciprocal of the describing function is just an approximation and its information can introduce errors in the estimation of the parameters. The solution adopted in this work to get accurately Gol ( jω osc ) during a test is presented by Vivek and Chidambaram31. As y(t) and u(t) are periodic and piecewise, using the Laplace transform of both, it can be written

2π ωosc

Y ( jωosc ) = Gol ( jωosc ) = U ( jωosc )

∫ ∫

y(t )e − jωosct dt

0 2π ωosc

(5) u(t )e

− jωosct

dt

0

and following (5), and as indicated in Wang et al.32, it is possible to obtain the harmonics

ACS Paragon17 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2π ωosc

Gol ( jnωosc ) =

∫ ∫

Page 18 of 71

y (t )e − jnωosct dt , n = 1,3,5,...

0 2π ωosc

u(t )e

− jnωosc t

(6)

dt

0

Expression (6) allows the solution of two problems: (i) getting the value of Gol ( jω osc ) that represents the exact point in the Nyquist plot where the intersection with − 1 N ( A, δ ) takes place, and (ii) obtaining additional points to work out the equations in just one test by calculating as many harmonics Gol ( jnω osc ) as needed. The last problem to solve is the estimation of the dc gain. Equation (5) cannot be applied when

ωosc = 0 due to the fact that the oscillations produced by the SSOD block are symmetric and the integration of the semi periods will be zero. However, when a biased relay test is used, the process static gain can be derived from (5)22,27,33,45 as

2π ωosc

Gol (0) =

∫ y(t) dt ∫ u(t) dt 0 2π ωosc

(7)

0

so

K=

Gol (0) K osc

(8)

Producing an asymmetric oscillation in the event-based approach during the test is enough to add a bias to the output of the SSOD block. It is important to notice that without bias, the signals

ACS Paragon18 Plus Environment

Page 19 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

are symmetric and even harmonics become zero. To avoid the dependence of the bias value, the harmonics used for

the estimation in the paper are the odd, i.e.,

Gol ( jω osc ) ,

Gol ( j 3ω osc ) , Gol ( j 5ω osc ) , ... Hence, to integrate the identification approach it will be necessary to include the previous expressions in the programming of the SSOD block. An example of code to make the eventbased sampling, the calculation of K with (8) and the three first odd harmonics with (6) is given in Listing I. The first if-structure corresponds to the event-based sampling, and the second to the detection of the period and computation of harmonics and oscillation frequency. It is supposed that the code in Listing I is run at a regular sampling period h. This sampling period h in the event-based approach is relevant because it is used to approximate the integrals of the control action and process output (see Listing I). Sampling period values of 0.01 and 0.001 have been used in the simulations. In general, the results are similar in both cases, but the estimation improves as h is reduced because the approximations of the integrals (6) and (7) are more accurate.

3.1 MODELS The event-based identification procedure can be adapted to most of the transfer function models that can be found in process industry by just obtaining the magnitude and argument expressions of the transfer function to fit. Expressions for a first order plus time delay (FOPTD), an overdamped second order plus time delay (SOPTD-1), an underdamped second order plus time delay (SOPTD-2), and an integrating process with inverse response and time delay (IPIRTD) transfer function are given in this paper.

ACS Paragon19 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 71

Model 1: FOPTD Ke − Ls G1 ( s ) = Ts + 1

(9)

Model 2: SOPTD-1 G2 (s) =

Ke − Ls (Ts + 1) 2

(10)

Model 3: SOPTD-2 Ke − Ls G3 ( s ) = 2 as + bs + 1

(11)

K (1 − T1 s )e − Ls (1 − T1 s )e − Ls G4 ( s) = = s (T2 s + 1) as 2 + bs

(12)

Model 4: IPIRTD

When K is needed for estimation of Models 1, 2, and 3, it is obtained directly from the test using (8) by adding a bias. The dead time L is derived from the argument expression

arg Gol ( jω osc ) of each model. To get the remaining parameters, the expressions of magnitude Gol ( jnω osc ) = K osc G ( jnω osc ) must be determined and solved with n =1 for Models 1 and 2, n =1, 3 for Model 3, and n=1, 3, 5 for Model 4. Notice that the experimental value of Gol ( jnωosc ) is obtained from the test (see expression (6) and Listing S1). More technical details are provided in Sánchez et al.16

The following expressions are the result of solving the magnitude and argument equations for the models. For the sake of simplicity, Cn represents Gol ( jnω osc ) and argC1 corresponds to

arg Gol ( jω osc ) .

ACS Paragon20 Plus Environment

Page 21 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Model 1: FOPTD

T=

L=−

2 − C12 K 2 K osc

(13)

ω osc C1

arg C1 + arctan( ω osc T )

ω osc

(14)

Model 2: SOPTD-1 T=

L=−

C1 (KK osc − C1 )

(15)

ω osc C1

2 arctan (T ω osc ) + arg C1

ω osc

(16)

Model 3: SOPTD-2

a=

2 2 2 2 2 2 2 2 1 2C1 K K osc − 18C3 K K osc+16C1 C3 2 12 ωosc C1C3

(17)

2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 − 2C1 K K osc+162C 3 K K osc − 160C1 C 3 + 24C1C 3 2C1 K K osc − 18C 3 KK osc+16C1 C 3 b= (18) 12 ωosc C1C 3

L=

(

)

2 arctan 2 − bωosc , − aωosc + 1 − arg C1 ωosc

Model 4: IPIRTD

ACS Paragon21 Plus Environment

(19)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1 T1 = ωosc

K a = osc 2 ωosc

b=

L=−

K osc ωosc

Page 22 of 71

3C12 C32 − 25C12 C52 + 150C32 C52 − 75C12 C32+225C12 C52 − 150C32 C52

(20)

− 2C12 + 27C32 − 25C52 − 225C12 C32 + 675C12 C52 − 450C32 C52

(21)

2C12 − 243C32 + 625C52 − 225C12 C32 + 675C12 C52 − 450C32 C52

(22)

( (

)

)

2 arctan 2 − T1 aωosc − b ,−ωosc(T1b+a) + arg C1 ωosc

(23)

where

K=

a 1 , T2 = b b

(24)

3.2. OUTLINE OF THE PROCEDURE The procedure is based on applying a gain K osc to G(s ) to produce a limit cycle oscillating at

ωosc as a consequence of the triggering of events when Gol (s) crosses the thresholds of the SSOD sampler. As shown before, mathematically, the limit cycle is a consequence of the intersection of Gol (s) with − 1 N ( A, δ ) at ωosc . The intersection point corresponds to the value of Gol ( jω osc ) . First, some considerations have to be done before explaining the complete procedure: (i) the process is supposed to be at the steady state with the SSOD block centered in zero, that means,

ACS Paragon22 Plus Environment

Page 23 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

e*(t)=0, (ii) the integral and derivative parts of the controller are disabled, and (iii) K osc is set to 1. The estimation procedure can be divided into the following steps:

1. If Models 1, 2 or 3 are going to be identified, add a small bias, e.g., ≈ 0.1δ to the SSOD block output. 2. To reach a stable oscillation, move out the process from its steady state introducing a square pulse of amplitude, for example, yr = ymax_r and increasing K osc . 3. Once the oscillation is stable, annotate K, K osc , ωosc , and the harmonics Gol ( jnω osc ) needed to solve the equations depending on the model. 4. Use expressions: - (13) and (14) for FOPTD model. - (15) and (16) for SOPTD-1 model. - (17), (18), and (19) for SOPTD-2 model. - (20), (21), (22), (23), and (24) for IPIRTD model. 5. Set the bias to zero and active the PID controller to proceed with the normal operation.

Remark 1: Model 3 allows detecting if the process corresponds to a dynamics of first order instead of second order. In this case, the procedure will generate a value for a close to zero.

Remark 2: Model 2 and 3 can produce the same result if the process corresponds to an overdamped second order system. In case of being underdamped, Model 2 produces complex values for the time lag and the correct results are produced by Model 3.

Remark 3: Model 4 identifies processes without inverse response by producing a value for T1 close to zero.

ACS Paragon23 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

3.3 SIMULATION EXAMPLES In this section, the proposed event-based identification procedure is applied to processes frequently found in industry. Table S1 shows the results when the structure of the true process and the model to fit are similar. In all the cases, the identification procedure has been done pursuing an oscillation with the theoretical maximum phase angle that the SSOD block allows. This maximum phase angle is 45º and corresponds to the angle of the first point of the SSOD DF plot (point C1 in Figure 4). Obviously, as the DF analysis is an approximation, the true point where the oscillation appears is different to the theoretical C1 point and depends on the process dynamics47. For this reason, the phase angles in Table S1 differ from 45º. When possible, the proposed method is compared with a previous version of this approach described by Sánchez et al.16. The results generated with the previous event-based approach have been obtained by oscillations of similar amplitude to the used for the present version (between δ and 2δ ) in order to be fair and show the improvement of the new approach. In the approach reported in Sánchez et al.16, it was necessary to force limit cycles of amplitude higher 5δ to get a similar accuracy to the results obtained with this procedure. As the example corresponding to Case 1 is often found in the literature, the result of the proposed method is compared with the results of different references. Models 2 and 3 produce the same parameters for Cases 2 and 3 as both real processes correspond to overdamped second systems. However, in Case 4, Model 2 produces complex values as the true process is underdamped so the solution is provided by Model 3. The application of Model 4 in Case 5 identifies exactly the true process dynamics if the positive zero is deleted, as an imaginary value

T1 = 0.0101 j results; it can be appreciated that the event-based procedure produces results of the

ACS Paragon24 Plus Environment

Page 24 of 71

Page 25 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

same quality as more elaborated methods based on state-space51 and curve-fitting47 approaches. In the simulations to detect Models 1, 2 and 3, the added bias to the SSOD output was 0.1δ . Sampling periods of 0.01s and 0.001s have been used to run the simulations. In general, the results are very similar in both cases but those obtained with h=0.001s are presented as the accuracy improves when the sampling frequency is increased because the integrations are more accurate. Table S2 presents the results obtained when the order of the true process is higher or equal than the transfer function to fit. The results are compared with solutions from other procedures based on different rationales. Case 7 is an example that is used by many authors to present outcomes, and for such reason so many results are included. Case 8 corresponds to the fitting of a FOPTD and SOPTD models. Cases 9, 10, and 11 to second order processes, and cases 12 and 13 to processes with integration and inverse response. In all the situations, the results produced by the event-based procedure generated accurate solutions. Figure 5 shows the Nyquist plot of two cases when the event-based procedure improves the results of other contrasted methods particularly for the steady and velocity gains.

ACS Paragon25 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(a)

(b)

(c)

ACS Paragon26 Plus Environment

Page 26 of 71

Page 27 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(d)

Figure 5. Nyquist curves and limit cycles of cases 11 (a, b) and 13 (c, d). It can be observed in (b) the asymmetric oscillation produced by a bias of 0.1δ .

4. EXTENSION OF THE PROCEDURE TO OTHER TRANSFER FUNCTION MODELS The identification procedure described in the previous section is based on solving a set of linear equations derived by equating the magnitude and argument expressions of the model to fit to some Nyquist points obtained from the experimental test with the oscillations of the true process. The procedure can be generalized to any transfer function order but, unfortunately, the length of the parameters expressions can be too long in some cases, especially, with five parameters or more. A solution to extend the approach to fit transfer function models with five or more parameters is to solve the equations system by the least-square method (LSM). The procedure for fitting the following transfer function

G (s) =

(−b1 s + 1)e − Ls a1 s 2 + a 2 s + a 3

ACS Paragon27 Plus Environment

(25)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 71

can be explained as follows. Note that (25) corresponds to a second order system with inverse response, which is a process that is frequently found in the chemical process industry and corresponds to a sum of two first order systems with inverse sign gains. The number of unknowns in (25) is four, as L is not taken into account because it is obtained by solving the equation corresponding to the argument of G ( jω osc ) . Therefore, four experimental points are required from the test, that is, G(0) , G ( jω osc ) , G ( j 3ω osc ) , and

G ( j 5ωosc ) . From (25), the magnitude expression of the steady gain and the three first harmonics of Gol (s) is

Gol ( jnω osc ) =

2 2 K osc (b12 n 2ω osc + 1) , n = 0,1,3,5 2 4 4 2 2 2 a1 n ω osc + n ω osc (a 2 − 2a1a 3 ) + a 32

(26)

Relating the magnitudes of the four experimental points with their corresponding magnitude expressions obtained from (26), and reordering the expressions, we obtain the following set of equations

2 C 02 a 32 = K osc 4 2 2 2 2 C12ω osc a12 + C12ω osc (a 22 − 2a1 a 3 ) + C12 a 32 − K osc b12 = K osc ω osc 4 2 2 2 2 C 32 3 4 ω osc a12 + C 32 3 2 ω osc (a 22 − 2a1 a 3 ) + C 32 a 32 − K osc 3 2 ω osc b12 = K osc 4 2 2 2 2 C 52 5 4 ω osc a12 + C 52 5 2 ω osc (a 22 − 2a1 a 3 ) + C 52 a 32 − K osc 5 2 ω osc b12 = K osc

ACS Paragon28 Plus Environment

(27)

Page 29 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Where Cn , with n=0,1,3,5, represents the magnitudes of the points Gol ( jnωosc ) , n=0,1,3,5, measured in the experimental test (see Listing I). The set of equations (27) contains no linear combinations of the parameters. However, we can define a vector

  2 a x= 2   

 a12  − 2a1 a 3   a32  2 b1 

(28)

such as the system is transformed into a linear set of equations

z = Ux

(29)

where

 0  2 4 2 C ω a U =  12 4osc 41 C 3 3 ω osc  2 4 4 C 5 5 ω osc

0

C 02

2 C12ω osc

C12

2 C 32 3 2 ω osc

C 32

2 C 52 5 2 ω osc

C 52

  2 2 ω osc − K osc  2 2  − K osc 3 2 ω osc 2 2  5 2 ω osc − K osc 

2  K osc   2  K  z =  osc 2   K osc  2   K osc 

which can be solved. Solutions are obtained by

ACS Paragon29 Plus Environment

0

(30)

(31)

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(

x = UT U

)

−1

UT z

Page 30 of 71

(32)

Once the unknowns are obtained from (32), we use the expression corresponding to the argument of G ( jω osc ) to get L. Manipulating (25), the expression of the argument is

2 2 − a3 b1 − a 2 ) , − ((a 2 b1 + a1 )ωosc − a3 )) − ωosc L arg Gol ( jωosc ) = − arctan(−ωosc (a1b1ωosc

(33)

so

L=

2 2 − a3b1 − a2 ), − ((a2 b1 + a1 )ωosc − a3 ) ) − arg Gol ( jωosc ) arctan(ωosc(a1b1ωosc ωosc

(34)

Example 2: The proposed identification algorithm is applied to the identification of a transfer function given by

G (s) =

(− s + 1)e − s as 2 + (1 + a) s + 1

(35)

for three different values of a={0.2, 0.35, 0.5}. This plant is used by Balaguer et al.59 as a numerical example to show the quality of an algorithm specifically designed for identification of second order inverse response systems using the step reaction curve; in that algorithm, the values of L and K are directly observed from the step response so, in this case, the results are exact. The comparison with the results from the event-based procedure is presented in Table S3 and in Figure 6.

ACS Paragon30 Plus Environment

Page 31 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 6. Comparative of the Nyquist and step response of the event-based estimation (solid line) with the true process for a=0.35 (dashed line).

5. SELECTING THE IDENTIFICATION PHASE LAG As explained before, by modifying K osc it is possible to change the identification frequency but there is a limit. According to Figure 4, the theoretical range of identification that permits the SSOD block goes from 45º to 0º. That corresponds to identifying the process with an oscillation of amplitude ± δ

to an oscillation of infinite amplitude (interpreted as instability).

Unfortunately, as K osc is increased, it becomes much more sensible being difficult not to enter in instability. That means that it is not possible to identify the process at ω −180º using as event-

ACS Paragon31 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

based sampler the current SSOD. Let analyse in this section why and how this restriction must be avoided in some situations in order to identify the process with a phase angle close to zero. The event-based identification approach using the SSOD sampler produces a transfer function model that corresponds exactly to the dynamics of the true process at the oscillation frequency, that is, G ( jωosc ) = Gˆ ( jωosc ) . If the order of the process is equal to the transfer function model to fit, the parameters identification will be exact and the behaviour of Gˆ ( s ) will be equal to G (s ) in all the frequencies range. However, if the order of process and model are different, it can be necessary to modify the oscillation frequency to find a more suitable approximation to the pursued control purposes. It is known that for PI control the identification should be done at an oscillation frequency that corresponds to a phase lag of 45º due to the phase lag that the PI control introduces. For PID control the phase lag should be 0º as the controller provides a phase lead46. As seen before, one of the features of the SSOD sampling is that an increase of the oscillation amplitude by K osc produces an increase in the oscillation frequency and the reduction of the phase angle. Unfortunately, there is a limit in the reduction of the phase angle with the current SSOD sampler. It can be explained by observing the Nyquist plot of − 1 N ( A, δ ) in Figure 4. By increasing K osc , Gol (s ) intersects with − 1 N ( A, δ ) in points that represent oscillations of different amplitudes and frequencies. To stay stable in one or another oscillation will depend on the operation or initial conditions of the system. However, as the amplitude of the oscillation is larger, the phase lag of the intersection point used for the identification will be closer to zero.

ACS Paragon32 Plus Environment

Page 32 of 71

Page 33 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(a)

(b)

Figure 7. Situations of G ( s ) = e − s (s + 1) with Kosc = 2 (a) and Kosc = 2.3 (b).

A theoretical example is shown in Figure 7. With K osc = 2 , the intersection point of

Gol ( s) = K osce− s (s + 1) in the arc segment corresponds to an oscillation of amplitude 1.1δ and

ωosc = 1.61 with a phase angle of 30º (point red in Figure 7.a); if K osc is increased to 2.13, the

ACS Paragon33 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 71

intersection point corresponding to an oscillation of amplitude 3.1δ (red point in Figure 7.b) and

ωosc = 1.891 presents a phase angle of 17º. However, in Figure 7.b it can be appreciated that the

Gol (s) can become unstable as the critical point -1 is to the right of the Nyquist curve. Example 3: By applying K osc = 2.53 to G ( s) = ( s + 1) −4 , the result of the identification is

1.002e −1.097 s ˆ G (s) = (1.594 s + 1) 2

(36)

with ϕ osc ≈ 35º , δ = 1 , and bias = 0.1δ . This experimental setup produces an asymmetric oscillation of amplitude [-1, +1.18] and ω osc = 0.7297 as a consequence of the control actions triggered for crossings the thresholds ± δ of the SSOD sampler. If K osc = 3.7 , the phase angle is reduced to 21º, the frequency is increased to ω osc = 0.8315 and the new result is

−1.081s

1.0007e Gˆ ( s ) = (1.642 s + 1) 2

(37)

The differences in the phase angles of both results can be appreciated in Figure 8. The point where the true process and the result coincide corresponds to the oscillation frequency used in the identification procedure.

ACS Paragon34 Plus Environment

Page 35 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(a)

(b)

Figure 8. (a) Details of the differences in the results of the identification of ( s + 1) −4 when the oscillation frequency is increased by modifying K osc . (b) General view of the Nyquist map. Dashed line corresponds to the true process.

Figure 9. Relationship between e(t) and e*(t) in a SSOD block where the hysteresis has been set to zero. As the hysteresis in the current SSOD block behaves as a phase-shifting function, a solution is to modify the event-based sampler by reducing the hysteresis to zero during the identification

ACS Paragon35 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 71

phase. In this way, the SSOD becomes a multi-level relay with dead-zone (Figure 9) with the following describing function60,61:

4δ  m (2k − 1) δ  ∑ 1 −  N ( A, δ ) =  πA  k =1 2 A  

2

   

(38)

where m =  A δ  . As the phase shifting is removed from the sampler, the representation of

− 1 N ( A, δ ) in the Nyquist plot becomes a straight line in the real axis, and hence, any intersection of Gol (s) with − 1 N ( A, δ ) will correspond to a point with phase angle of 0º. Example 4: Using the SSOD sampler without hysteresis, G ( s) = ( s + 1) −4 is now identified as

−1.047 s

1.005e Gˆ ( s ) = (1.733s + 1) 2

(39)

with ϕ osc ≈ 1º , ω osc = 0.990 , δ = 1 , bias = 0.1δ , and K osc = 2.53 (see Figure 8). So, as shown in Example 4, the SSOD sampler without hysteresis allows the identification with a phase angle near zero. If the identification must be done with another phase angle, the hysteresis must be included into the SSOD sampler, varying K osc to modify the oscillation amplitude and the phase lag.

Remark 4. It is important to underline that the approach can be applied in the more traditional relay-feedback framework. The main difference is that a standard relay approach will generate the oscillations near the ultimate frequency with a phase angle ϕosc ≈ 0º . So, the adaptation of the procedure to the typical setup based on a simple relay with adaptable hysteresis is direct.

ACS Paragon36 Plus Environment

Page 37 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

To analyse the quality of the identification approach as part of an autotuning technique, it is necessary to study the influence on the controller tuning when the phase angle of the critical point used to identify is modified. To do that a comparison of PI and PID controllers is made for FOPTD models of the process G( s) = ( s + 1) −4 as it is considered a balanced lag and delay process22. PI and PID controllers were tuned for the models using the AMIGO-step rules46. The models, the controller parameters and the robustness measures are listed in Table S4. For comparison purposes, as in Beschi et al.8, the maximum sensitivities of the actual open loop C ( s )G ( s ) (denoted as Ms) and of the modelled open loop C ( s )Gˆ ( s ) (denoted as Msa) are presented in Table S4. It can be appreciated that the values of maximum sensitivity obtained are quite similar and any of the models is sufficiently accurate to tune the controller effectively. Regarding the tuning quality, to choose one or another Gˆ ( s ) for application of the AMIGO rules is not very relevant as it can be appreciated in Table S4. As it can be appreciated in Figure 10, the frequency responses around the crossover frequency in the four cases are very similar and the models obtained by the event-based approach present the same behaviour as the model Gˆ 0 ( s) obtained by the more advanced method described in Berner et al.22. Thus, to select a model Gˆ ( s ) will depend on how close the behaviours of the model and process at a specific frequency, that is, at

ω osc , are required.

ACS Paragon37 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 38 of 71

(b)

(a)

Figure 10. Comparison between the Nyquist plots of C ( s )G ( s ) using the tuning parameters obtained for each of the four models by AMIGO rules: (a) with PI (b) with PID. The system tuned with Gˆ 0 ( s) is almost coincident with the system tuned with Gˆ 2 ( s ) .

6. INFLUENCE OF THE NOISE Event-based control systems play a major role when wireless components are used, as they represent a practical solution to reduce the exchange of information between the different involved control agents. However, as for wired systems, the presence of noise is a critical issue. One of the advantages of the SSOD sampling is that the hysteresis of the sampler (see Figure 3) reduces the influence of the measurement noise by producing a control action free of noise in most of the cases. However, a high level of noise, for example, with n0 > 0.5δ produces false event triggerings and the calculation of, for example, the frequency of the limit cycle would

ACS Paragon38 Plus Environment

Page 39 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

become impossible. Such fact does not allow the application of the identification method without a previous filtering of the process output. To analyse the practical influence of noise on the parameter estimation approach, simulations have been conducted for the balanced process G ( s) = ( s + 1) −4 . FOPTD and SOPTD models in the presence of noise have been obtained and compared with the noise-free models presented in the previous section. To introduce noise in the simulations, band-limited white noise was connected to the process output with different levels, n0 = [0.1δ ,0.3δ ,0.5δ ] with δ = 1 and h=0.01s. As a consequence of the variability that the noise can introduce in the integrals (12) and (13), the approach was modified to use the mean value of the integrals in the last five consecutive periods of the simulations instead of filtering the input to the SSOD block. The results are presented in Table S5 and some of the simulations are illustrated in Figure 11. The noise does not deteriorate excessively the accuracy of the estimates. With respect to Gˆ 0 ( s) models, the maximum variation in the parameters is of 20% in T with the higher level of noise

n0 = 0.5δ . However, the dead time and the steady gain do not present significant variations. It must be noticed that it was necessary for n0 = 0.5δ to increase the oscillatory gain from

K osc = 2.53 to K osc = 3.7 because the limit cycle vanished, not reaching the necessary stable oscillation (Figures 11.b and 11.c).

ACS Paragon39 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 71

(a)

(b)

(c)

Figure 11. Simulations for the identification of G ( s) = ( s + 1) −4 in the presence of different levels of noise: (a) n0 = 0.3δ and K osc = 2.53 , (b) n0 = 0.5δ and K osc = 2.53 , (c) n0 = 0.5δ and

K osc = 3.7 .

It must be noticed that the use of frequencies 3ωosc and 5ωosc can introduce problems to produce estimations of SOPTD-2 and IPIRTD models as they need the points Gol(j3ωosc) and Gol (j5ωosc) to generate the parameters. To analyze that, two situations have been checked: (1) to identify

a

IPIRTD

model

G A ( s ) = 0.6 (−0.4s + 1) [ s(s + 1)] ,

of

a

and

(2)

process to

of

identify

the a

same process

order, of

that

higher

is, order

G B ( s ) = ( −5s + 1)e −0.5 s [ s (s + 1)( s 2 + s + 1)] . These two processes are the cases 6 and 13,

identified in Tables S1 and S2, respectively, without noise. ACS Paragon40 Plus Environment

Page 41 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Regarding the estimation of G A (s) , for n0 = [0.1δ ,0.3δ ] with δ = 1 , K A _ osc = 1.5 , and h=0.01s the identified models are Gˆ A _ 0.1δ ( s ) = 0.567 (−0.255s + 1)e −0.323s [ s (0.761s + 1)] and Gˆ A _ 0.3δ = 0.443 (−0.326 s + 1)e −0.317 s [ s (0.501s + 1)] ,

K B _ osc = 0.126 , the models are

respectively.

For

G B (s)

Gˆ B _ 0.1δ = 0.920 (−6.131s + 1)e −0.945 s [ s (0.995s + 1)]

with and

Gˆ B _ 0.3δ = 0.835 ( −7.091s + 1)e −0.486 s [ s (1.227 s + 1)] . In general, the estimations are correct in the bandwidth around the oscillation frequency. It must be noticed that the estimations of

G B (s) are better than of G A (s) . The reason is G B (s) owns a slower dynamics and therefore the oscillation frequency is lower than that of G A (s) , reducing the impact of the noise in the calculations of Gol ( j 3ω osc ) and Gol ( j 5ω osc ) . For example, for n0 = 0.1δ , these frequencies are

ω A _ osc = 0.7 and ω B _ osc = 0.2 . It can be observed than the when the limit cycle is of amplitude ± δ , the introduction of a bias in presence of noise can improve the estimation. The bias forces an asymmetric output that produces a reduction of the oscillation frequency. So, for example, the model of G A (s) by adding a bias of 0.1δ with n0 = 0.1δ is Gˆ A _ 0.1δ ( s ) = 0.602 (−0.3838s + 1)e −0.001 [ s (1.0191s + 1)] with ω A _ osc = 0.25 (notice that, without bias, ω A _ osc was 0.7)

CONCLUSIONS An enhanced event-based method for identification of stable transfer function models of any order and structure has been presented and explained in detail. The basis of the method is the oscillations that an event-based sampler can generate in the feedback control loop, so the describing function is the tool used to explain the oscillation phenomena and the procedure of

ACS Paragon41 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

estimation. Instead of using the approximations of the intersection points that the describing function of the event-based sampler provides, on-line measurements of the process and control signals are used to obtain the harmonics needed to solve the linear equations and generate the estimations. In this way, the number of experiments to make the estimations is always one, regardless of the order or structure of the model transfer function to fit. Moreover, the identification procedure is not iterative so that the computational cost is very low. That issue, the good results in presence of noise, and the event-based nature of the method convert the procedure into a valuable tool for working in networked control systems. It is relevant to underline that the event-based procedure described here should be considered as the generic case of an identification method based on a relay. By replacing the send-on-delta block by a simple relay, the identification procedure will work in a similar way but the identification will be done near the ultimate frequency; and by using a relay with hysteresis and modifying its value, the identification will be possible at a user-specified phase angle in the third quadrant from 90º to 0º. Further works will be focused in providing an autotuning method for event-based PID controllers based on this identification strategy and continue the study of the impact of the bias in the identification result in presence of noise.

ASSOCIATED CONTENT Supporting Information Listing S1, Tables T1, T2, T3, T4, and T5 (PDF)

AUTHOR INFORMATION Corresponding Author *E-mail: [email protected]. Phone: +34 913 987 146

ACS Paragon42 Plus Environment

Page 42 of 71

Page 43 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

ORCID José Sánchez: 0000-0002-6702-3771 María Guinaldo: 0000-0002-7043-6673 Antonio Visioli: 0000-0002-9246-5715 Sebastián Dormido: 0000-0002-2405-8771

Notes The authors declare no competing financial interest.

ACKNOWLEDGEMENTS This work has been funded by the Spanish Ministry of Economy and Competitiveness under contract DPI2014-55932-C2-2-R and DPI2017-84259-C2-2-R.

ABBREVIATIONS A = Amplitude of a limit cycle D = Relay output

δ = Threshold value of an event-based block DF = Describing Function FOPTD = First Order Plus Time Delay

G ( 0 ) = Process gain at stationary state

Gˆ ( s ) = Estimated transfer function IPIRTD = Integrating Process with Inverse Response and Time Delay

K osc = Gain applied to force a limit cycle K u = Controller critical gain K −180 º = Process gain corresponding to a phase lag of -180º

ϕ osc = Phase angle of the critical point PID = Proportional-Integral-Derivative

ACS Paragon43 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

SOPTD = Second Order Plus Time Delay SSOD = Symmetric Send-On-Delta

ω osc = Oscillation frequency of a limit cycle ω −180 º = Oscillation frequency corresponding to a phase margin of -180º ω −135 º = Oscillation frequency corresponding to a phase margin of -135º ωˆ u = Estimated ultimate frequency ω u = Ultimate frequency

REFERENCES (1) Liu, T.; Wang, Q. G.; Huang; H. P. A Tutorial Review on Process Identification from Step or Relay Feedback Test. J. Process Control 2013, 23, 1597-1623. (2) Event-Based Control and Signal Processing; Miskowicz, M, Ed.; CRC Press: Boca Raton, FL, 2015. (3) Blevins, T.; Chen, D.; Nixon, M., Wojsznis, W.K. Wireless Control Foundation Continuous and Discrete Applications; International Society of Automation (ISA): Research Triangle Park, NC, 2014. (4) Vasyutynskyy, V.; Kabitzsh, K. Implementation of PID Controller with Send-on-Delta Sampling. Proceedings of the International Control Conference: Glasgow, UK, Aug 30-Sept 1, 2006; University of Strathclyde Publishing, 2006. (5) Miskowicz, M. Send-on-Delta: An Event-Based Data Reporting Strategy. Sensors 2006, 6, 49-63.

ACS Paragon44 Plus Environment

Page 44 of 71

Page 45 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(6) Chacón, J.; Sánchez, J.; Visioli, A.; Yebra, L.; Dormido, S. Characterization of Limit Cycles for Self-Regulating and Integral Processes with PI Control and Send-on-Delta Sampling. J. Process Control 2013, 23, 826-838. (7) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A. Stability Analysis of Symmetric Send-onDelta Event-Based Control Systems. Proceedings of the American Control Conference; Washington, DC, June 17-19, 2013; IEEE, 2013; 1771-1776. (8) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A. Closed-Loop Automatic Tuning Technique for an Event-Based PI Controller. Ind. Eng. Chem. Res. 2015, 54, 6362-6370. (9) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A. Tuning of Symmetric Send-on-delta Proportional-Integral Controllers. IET Control Theory Appl. 2014, 8, 248-259. (10) Romero, J.A.; Sanchís, R.; Arrebola, E. Experimental Study of Event Based PID Controllers with Different Sampling Strategies. Application to Brushless DC Motor Networked Control System. Proceedings of the International Conference on Information, Communication and Automation Technologies; Sarajevo, Bosnia & Hezergovina, Oct 29-31, 2015; IEEE, 2015. (11) Pawlowski, A.; Beschi, M.; Guzmán, J. L.; Visioli, A.; Berenguel, M.; Dormido S. Application of SSOD-PI and PI-SSOD Event-Based Controllers to Greenhouse Climatic Control. ISA Trans. 2016, 65, 525-536. (12) Ruiz, A.; Beschi, M.; Visioli, A.; Dormido, S.; Jiménez J. E. A Unified Event-Based Control Approach for FOPTD and IPTD Processes Based on the Filtered Smith Predictor. J. Franklin Inst. 2017, 354, 1239-1264.

ACS Paragon45 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(13) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A. Characterization of Symmetric Send-onDelta PI Controllers. J. Process Control 2012, 22, 1930-1945. (14) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A.; Yebra, L. J. Event-based PI plus Feedforward Control Strategies for a Distributed Solar Collector Field. IEEE Trans. Control Syst. Technol. 2014, 22, 1615-1622. (15) Beschi, M.; Dormido, S.; Sánchez, J.; Visioli, A. An Event-based PI Controller Autotuning Technique for Integral Processes. Proceedings of the 1st IEEE International Conference on Event-Based Control, Communication, and Signal Processing; Krakow, Poland, June 17-19; IEEE, 2015. (16) Sánchez, J.; Guinaldo, M.; Visioli, A.; Dormido, S. Identification of Process Transfer Function Parameters in Event-based PI Control Loops. ISA Trans. 2018, 75, 157-171. (17) Jeon, C.H; Cheon, Y. J.; Kim, J. S.; Lee, L.; Sung, S. W. Relay Feedback Methods Combining Sub-relays to Reduce Harmonics. J. Process Control 2010, 20, 228-234. (18) Åström, K.J.; Hägglund, T. Automatic Tuning of Simple Regulators with Specifications on Phase and Amplitude Margins. Automatica 1984, 20, 645-651. (19) Chidambaram, M.; Sathe, V. Relay Autotuning for Identification and Control; Cambridge University Press: Delhi, India, 2014. (20) Mehta, U.; Ananthanarayanan, R. Optimal Autotuning of Nonminimum Phase Processes under Relay Control. Proc. Inst. Mech. Eng. I 2015, 229, 1-14.

ACS Paragon46 Plus Environment

Page 46 of 71

Page 47 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(21) Hofreiter, M. Shifting Method for Relay Feedback Identification. Proceedings of the 8th IFAC Conference on Manufacturing Modelling, Management and Control; Troyes, France, June 28-30, 2016; IFAC, 2016. (22) Berner, J.; Hägglund T.; Åström, K. Improved Relay Autotuning using Normalized Time Delay. Proceedings of the American Control Conference; Boston, MA, July 6-8, 2016; IEEE, 2016; 1869-1875. (23) Bajarangbali, R.; Majhi, S. Estimation of First and Second Order Systems. Proc. Natl. Acad. Sci., India, Sect. A Phys. Sci. 2017. DOI: 10.1007/s40010-017-0357-6. (24) Chang, R. C.; Shen, S. H.; Yu, C. C. Derivation of Transfer Function from Relay Feedback Systems. Ind. Eng. Chem. Res. 1992, 31, 855-860. (25) Luyben, W. L. Derivation of Transfer Functions for Highly Nonlinear Distillation Columns. Ind. Eng. Chem. Res. 1987, 26, 2490-2495. (26) Li, W.; Eskinat, E.; Luyben, W. L. An Improved Autotune Identification Method. Ind. Eng. Chem. Res. 1991, 3, 1530-1541. (27) Shen, S. H.; Wu, J.S.; Yu, C .C. Use of Biased-Relay Feedback for System Identification. AIChE J. 1996, 42, 1174-1180. (28) Scali, C.; Marchetti, G.; Semino, D. Relay with Additional Delay for Identification and Autotuning of Completely Unknown Processes. Ind. Eng. Chem. Res. 1999, 38, 1987-1997. (29) Srinivasan, K.; Chidambaram, M. An Improved Autotune Identification Method. Chem. Biochem. Eng. Q. 2004, 18, 249-256.

ACS Paragon47 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(30) Srinivasan, K.; Chidambaram, M. Modified Relay Feedback Method for Improved System Identification. Comput. Chem. Eng. 2003, 27, 727-732. (31) Vivek, S.; Chidambaram, M. Identification Using Single Symmetrical Relay Feedback Test. Comput. Chem. Eng. 2005, 29, 1625-1630. (32) Wang, P.; Gu, D.; Zhang, W. Modified Relay Feedback Identification Based on Describing Function Analysis. Ind. Eng. Chem. Res. 2007, 46, 1538-1546. (33) Bajarangbali, R.; Majhi, S. Relay Based Identification of Systems.
Int. J. Scientific Eng. Res. 2012, 3, 1-4. (33) Marchetti, G.; Scali, C. Use of Modified Relay Techniques for the Design of Model-Based Controllers for Chemical Processes. Ind. Eng. Chem. Res. 2000, 39, 3325-3334. (34) Bajarangbali, R.; Majhi S. Modeling of Stable and Unstable Second Order Systems with Time Delay. Proceedings of the Annual IEEE India Conference; Bombay, Munday, India, Dec. 13-15, 2013; IEEE, 2014. (35) Pandey, S.; Majhi, S. Limit Cycle Based Identification of Second Order Processes with Time Delay. Proceedings of the Indian Control Conference; Hyderabad, India, Jan 4-6, 2016; IEEE, 2016; 438-443. (36) Ghorai, P.; Majhi, S.; Pandey, S. Dynamic Model Identification of a Real-Time Simple Level Control System. J. Control Decision 2016, 3, 248-266. (37) Tan, K. K.; Lee, T. H.; Wang, Q. G. An Enhanced Automatic Tuning Procedure for PI/PID Controllers for Process Control. AIChE J. 1996, 42, 2555-2562.

ACS Paragon48 Plus Environment

Page 48 of 71

Page 49 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(39)

Leva

A.;

Bascetta,

L.;

Schiavo,

F.

Model-Based

Proportional−Integral/

Proportional−Integral−Derivative (PI/PID) Autotuning with Fast Relay Identification. Ind. Eng. Chem. Res. 2006, 45, 4052-4062. (40) Friman, M.; Waller, K. V. A Two-Channel Relay for Autotuning. Ind. Eng. Chem. Res.

1997, 36, 2662-2671. (41) Wang, Y.; Shao, H. PID Autotuner Based on Gain- and Phase-Margin Specifications. Ind. Eng. Chem. Res. 1999, 38, 3007-3012. (42) Sung, S. W.; Lee, J. Two-Channel Relay Feedback Method under Static Disturbances. Ind. Eng. Chem. Res. 2006, 45, 4071-4074. (43) Lee, J.; Kim, J. S.; Byeon, J.; Sung, S. W.; Edgar, T. F. Relay Feedback Identification for Processes under Drift and Noisy Environments. AIChE J. 2011, 57, 1809-1816. (44) Byeon, J.; Kim, J.; Sung, S.; Ryoo, W.; Lee, J. Third Quadrant Nyquist Point for the Relay Feedback Autotuning of PI Controllers. Korean J. Chem. Eng. 2011, 28, 342-347. (45) Kaya, I.; Atherton, D. P. Parameter Estimation from Relay Autotuning with Asymmetric Limit Cycle Data. J. Process Control 2001, 11, 429-439. (46) Kaya, I. Parameter Estimation for Integrating Processes Using Relay Feedback Control under Static Load Disturbances. Ind. Eng. Chem. Res. 2006, 45, 4726-4731. (47) Panda, R. C.; Vijayan, V.; Sujatha, V.; Deepa, P.; Manamali, D.; Mandal, A. B. Parameter Estimation of Integrating and Time Delay Processes Using Single Relay Feedback Test. ISA Trans, 2011, 50, 529-537.

ACS Paragon49 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(48) Liu, T.; Gao, F. Industrial Process Identification and Control: Design Step-test and Relay-experiment-based Methods. Springer-Verlag: London, UK, 2012. (49) Majhi, S. Relay Based Identification of Processes with Time Delay. J. Process Control

2007, 17, 93-101. (50) Bajarangbali, R;. Majhi, S.; Pandey, S. Identification of FOPDT and SOPDT Process Dynamics Using Closed Loop Test. ISA Trans. 2014, 53, 1223-1231. (51) Bajarangbali, R.; Majhi, S. Identification of Integrating and Critically Damped Systems with Time Delay. Control Theory Tech. 2015, 13, 29-36. (52) Pandey, S.; Majhi, S. Limit Cycle-Based Exact Estimation of FOPDT Process Parameters under Input/output Disturbances: A State-Space Approach. Int. J. Syst. Sci. 2017, 48, 118-128. (53) Soltesz, K.; Mercader, P.; Baños, A. An Automatic Tuner with Short Experiment and Probabilistic Plant Parameterization. Int. J. Robust Nonlin 2016, 27, 1857-1873. (54) Berner, J.; Soltesz, K. Short and Robust Experiments in Relay Autotuners. Proceedings of the 22nd IEEE International Conference on Emerging Technologies and Factory Automation; Limassol, Cyprus, Sept 12-15, 2017; IEEE: 2018. (55) Ionescu, C. M.; De Keyser, R.; Robayo, F.; Milica B.; Naumovic, M. B. The Transfer Function Analyzer Revisited. Proceedings of 18th Mediterranean Conference on Control & Automation; Marrakech, Morocco, June 23-25, 2010; IEEE: 2010. (56) Åström, K. J.; Hägglund, T. Advanced PID Control; International Society of Automation (ISA): Research Triangle Park, NC, 2006.

ACS Paragon50 Plus Environment

Page 50 of 71

Page 51 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(57) Romero, J.A.; Sanchís, R. A New Method for Tuning PI Controllers with Symmetric Send-on-Delta Sampling Strategy. ISA Trans. 2016, 64, 161-173. (58) Gu, D.; Ou, L.; Wang, P.; Zhang, W. Relay Feedback Autotuning Method for Integrating Processes with Inverse Response and Time Delay. Ind. Eng. Chem. Res. 2006, 45, 3119-3132. (59) Balaguer, P.; Alfaro, V.; Arrieta, O. Second Order Inverse Response Process Identification from Transient Step Response. ISA Trans. 2011, 50, 231-238. (60) Gelb, A.; Van der Velde, W. E. Multiple-Input Describing Functions and Nonlinear System Design; McGraw-Hill: New York, NY, 1968. (61) Gibson, J. E. Nonlinear Automatic Control; McGraw-Hill: New York, NY, 1963.

ACS Paragon51 Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Note: The TOC/ABS has been designed to fit in an area 3.81 cm x 8.46 cm but the size has been increased the 100% in this document.

ACS Paragon52 Plus Environment

Page 52 of 71

Page 53 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 1. Control scheme of the event-based control loop. The dashed line indicates that the data sent from the sampler to the controller are not synchronous. 86x29mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2. Differences between true process and estimated transfer function when G(0) is not considered in the fitting. 95x75mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 54 of 71

Page 55 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 3. Relationship between e(t) and e∗(t) in a SSOD block. 61x60mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4. Nyquist plot of -1/N(A,δ). 84x66mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 56 of 71

Page 57 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 5.a. Nyquist curves and limit cycles of cases 11 (a, b) and 13 (c, d). It can be observed in (b) the asymmetric oscillation produced by a bias of 0.1δ. 99x79mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5.b. Nyquist curves and limit cycles of cases 11 (a, b) and 13 (c, d). It can be observed in (b) the asymmetric oscillation produced by a bias of 0.1δ. 89x32mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 58 of 71

Page 59 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 5.c. Nyquist curves and limit cycles of cases 11 (a, b) and 13 (c, d). It can be observed in (b) the asymmetric oscillation produced by a bias of 0.1δ. 100x79mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5.d. Nyquist curves and limit cycles of cases 11 (a, b) and 13 (c, d). It can be observed in (b) the asymmetric oscillation produced by a bias of 0.1δ. 89x30mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 60 of 71

Page 61 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 6. Comparative of the Nyquist and step response of the event-based estimation (solid line) with the true process for a=0.35 (dashed line). 85x121mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7.a. Situations of G(s)=e-s/(s+1) with Kosc=2 (a) and Kosc=2.3 (b). 88x75mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 62 of 71

Page 63 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 7.b. Situations of G(s)=e-s/(s+1) with Kosc=2 (a) and Kosc=2.3 (b). 89x69mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8. (a) Details of the differences in the results of the identification of (s+1)-4 when the oscillation frequency is increased by modifying Kosc. (b) General view of the Nyquist map. Dashed line corresponds to the true process. 89x67mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 64 of 71

Page 65 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 8. (a) Details of the differences in the results of the identification of (s+1)-4 when the oscillation frequency is increased by modifying Kosc. (b) General view of the Nyquist map. Dashed line corresponds to the true process. 86x67mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9. Relationship between e(t) and e*(t) in a SSOD block where the hysteresis has been set to zero. 61x60mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 66 of 71

Page 67 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 10. Comparison between the Nyquist plots of C(s)G(s) using the tuning parameters obtained for each of the four models by AMIGO rules: (a) with PI (b) with PID. The system tuned with G0(s) is almost coincident with the system tuned with G2(s). 84x82mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10. Comparison between the Nyquist plots of C(s)G(s) using the tuning parameters obtained for each of the four models by AMIGO rules: (a) with PI (b) with PID. The system tuned with G0(s) is almost coincident with the system tuned with G2(s). 83x81mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 68 of 71

Page 69 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 11. Simulations for the identification of (s+1)-4 in the presence of different levels of noise: (a) n0=0.3δ and Kosc=2.53, (b) n0=0.5δ and Kosc=2.53, (c) n0=0.5δ and Kosc=3.7. 79x37mm (300 x 300 DPI)

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 11. Simulations for the identification of (s+1)-4 in the presence of different levels of noise: (a) n0=0.3δ and Kosc=2.53, (b) n0=0.5δ and Kosc=2.53, (c) n0=0.5δ and Kosc=3.7. 79x37mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 70 of 71

Page 71 of 71 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 11. Simulations for the identification of (s+1)-4 in the presence of different levels of noise: (a) n0=0.3δ and Kosc=2.53, (b) n0=0.5δ and Kosc=2.53, (c) n0=0.5δ and Kosc=3.7. 79x37mm (300 x 300 DPI)

ACS Paragon Plus Environment