PID Tuning Method Based on IMC for Inverse-Response Second-Order Plus Dead Time Processes

: This work addresses a set of tuning rules for PID controllers based on Internal Model Control (IMC) for inverse-response second-order systems with dead time. The transfer function, and some time-response characteristics for such systems are ﬁrst described. Then, the IMC-based methodology is developed by using an optimization objective function that mixes performance and robustness. A correlation that minimizes the objective function and that allows the user to compute the controller’s tuning parameter is found. The obtained expressions are mathematically simple, which facilitate their application in a ten-step systematic methodology. Finally, the proposed tuning method is compared to other well-known tuning rules that have been reported in literature, for a wide range of parameters of the process. The performance achieved with the proposed method is very good not only for disturbance rejection but for set-point tracking, when considering a wide-range of parameters of the process’ transfer function.


Introduction
Inverse response or non-minimum-phase systems have at least one zero located at the right-hand side of the complex plane [1]. Their time response to a step input goes, at the beginning, towards the opposite direction of the steady-state value [2]. Such behavior appears due to the combination of two opposite-dynamics phenomena: a small-gain fast-dynamics responsible for the inverse response, and a slow-dynamics with higher gain, which dominates the transient response [3][4][5].
Since the inverse response phenomenon presents some similar characteristics to those in dead-time systems, it is usual to find tuning techniques with Smith predictor-like structures [6,7]. Other tuning techniques present the model of the non-minimum-phase zero as dead time, which allows one to use traditional tuning techniques developed for time-delayed systems [8]; conservative controllers are obtained with such an approximation [9][10][11]. Another traditional way to deal with inverse response dynamics is to reduce the order of the system and get a first-order model, which allows one to obtain PI controllers; however, these controllers can cause slower transient responses than the aforementioned techniques [12].
The development of Proportional-Integral (PI) and Proportional-Integral-Derivative (PID) controllers for inverse response processes started around the 1970s. Waller and Nygardas [13] compared lead-lag compensators to the conventional PID tuning proposed by Ziegler and Nichols [14]. Then, Scali and Rachid [15], Luyben [16], Chien et al. [17], and Sree and Chidambaram [18], among others, reported different tuning techniques for PI/PID controllers used for inverse-response systems by considering, in general, set-point changes. One can find some references, such as Chen and Seborg [19], Shamsuzzoha and Lee [20], and Pai et al. [9] that developed tuning equations for disturbance rejection. More recently, several PID tuning or controller design methods for processes whose dynamics include inverse-response, time-delay, and integrating characteristics have been developed [21][22][23][24][25][26]. Other authors gone beyond the PID controller and proposed the use of fractional control for non-minimum phase plus dead time systems [27], and Sliding Mode Controllers applied to high-order long dead-time inverse-response processes [28].
Other PID design methods based on frequency response techniques have been developed. For instance, Luyben [29] presented a PI controller tuning procedure for an inverse-response integrating process; the author states that the controller tuning process for this kind of complex dynamics is not a trivial task. Chen and Seborg [19] developed a design method for PID controllers based on direct synthesis; they demonstrated that the design method works very good for several processes, including those with non-minimum phase characteristics that appear when Pade's approximation is used for systems with dead time. Such approximation induces an additional modeling error, which implies a decrease in the system's bandwidth [30]. Alfaro et al. [31] used a nondimensional version of equations proposed by Chen and Seborg [19] and found expressions in order to compute the minimum value of the tuning parameter, τ c , which bounds the maximum sensitivity.
Internal Model Control (IMC) has been used as an alternative for designing and tuning PID controllers since 1980s [30]. It has resulted to be of particular interest in industry together with the PID algorithm [32], since the equations for the controller's parameters can be obtained from the transfer function of the process and the desired behavior of the closed-loop response; in most cases, only the closed-loop time constant is required as the user-defined tuning parameter, considering an appropriate trade-off between performance and robustness [20,[33][34][35][36][37]. Additional works, regarding IMC, that have been developed more recently can be found in [38][39][40][41][42][43][44].
This work addresses the design of a tuning methodology for PID controllers, for inverse-response second-order plus dead time processes, extending the work of Chien and Fruehauf [45], which has been recently studied in [46], by using dimensional analysis and experimental design. The development of the tuning technique considers the optimization of an objective function that combines performance and robustness indexes: the Integral of the Absolute Value of the Error (IAE) and the Integral of the change in the Manipulated Variable (IMV). The obtained expressions are simple and can be applied by using the proposed ten-step methodology. The performance of a controller tuned with this method excels the performance offered by other well-known techniques for a wide range of parameters of the process.
The first section contains a description of process' time response. Then, the design process of the tuning rules is addressed, by using dimensional analysis and experimental design. The third section contains the application methodology and the performance comparison for a controller tuned with the new rules and other tuning methods. Then, the robustness analysis for the variability of parameters is performed. Finally, the conclusions are presented.

Inverse Response Second Order Plus Dead Time System Model
The transfer function of an inverse-response second-order system, with total time delay θ, can be written as where are the transfer functions of two opposite-dynamics phenomena: a small-gain fast-dynamics part (3), and a slow-dynamics part with higher gain (2). K 1 and K 2 are the gains, and τ 1 and τ 2 are the time constants. Substituting (2) and (3) into (1) yields As it can be seen in (4), the process' transfer function only considers the case with real different poles because the dynamics comes from parallel balances [47]; this is a common case in industrial processes containing boilers, heat exchangers, tanks, distillation columns, chemical reactors, waste incinerators, among others [2,47].
Operating expressions in (4) yields the transfer function of the process which is given by where For this process, at the beginning, the time response slope's sign is opposite to the sign of the gain, and as η increases, so does the inverse-response behaviour; then, control issues arise due to the fact that the non-minimum phase zero is closer to the imaginary axis [4,15,48,49]. Adding a left-hand side pole into (5) decreases this inverse-response effect, but the settling time increases. When there are multiple non-minimum phase zeros, multiple inversions appear; however, such a system is not commonly found.

Internal Model Control
Internal Model Control (IMC) is a technique that uses a simulation of the process, running in parallel [30], as depicted in Figure 1; where: R(s) is the input reference, D(s) is the disturbance signal, Y(s) is the output of the system, G P (s) is the process transfer function, G D (s) is the transfer function that relates the output and the disturbance, G C (s) is the controller transfer function, and G M (s) is the model of the process that is simulated in parallel. The IMC methodology requires the transfer function of the model of the process be written as: where G M + (s) is a transfer function containing all the non-minimum phase zeros and dead time, and G M − (s) contains the minimum-phase elements of the system. Then, the IMC controller can be defined in terms of the invertible part of the transfer function and a user-specified low-pass filter f (s), as follows The IMC controller's structure can be reduced to a classical PI/PID control structure. For this purpose, Chien and Fruehauf [45] defined the filter's function, for first and second-order systems (without integrator), as follows They used a parallel PID controller structure which is given by These authors obtained equations for the tuning parameters, which are given for an inverse response second order plus dead time system by [45] where τ c is the tuning parameter.

Optimization Objective Function
In this work, we use the optimization of a performance index that combines the Integral of the Absolute Value of the Error (IAE) and the Integral of the change in the Manipulated Variable (IMV). The controller's main operation in this case favors disturbance rejection (load regulation); however, reference tracking performance is also evaluated. The objective function is then given by The first term in (13) is the IAE, and the second term is the IMV multiplied by a suppression factor γ. Increasing or decreasing γ allows one to modify the weight on the manipulated variable; for instance, when γ = 1, both, the IAE and IMV are of the same importance, and when γ = 0 only the IAE is considered. Several authors have used this type of objective function in order to measure the performance of different control strategies [50,51]; typical values of the suppression factor are 0 ≤ γ ≤ 8 [52,53].

Scaling and Nondimensionalization
In order to reduce the amount of parameters to be optimized, Equations (5), (10)- (12) are scaled using dimensional analysis techniques, allowing us to go from the 4 parameters contained in (5) to the 3 parameters found in Equations (18)- (20); this reduces the number of simulation runs and computation time. Nondimensional equations are written in the form proposed in [49,54] and given bŷ where:η This can be interpreted as scaling in time or frequency domains. The typical values for the ratios between the parameters can be found in literature and are given by [9,49] A central composite experimental design (uniform and circumscribed for rotatability) was designed in order to provide appropriated values of τ 2 for each simulation run. This type of experimental design is selected since it can take account of curvature in a response [55], i.e., a quadratic surface can be adjusted to the data provided. Gutierrez and de la Vara [56] stated that in order to get accurate predictions, values of R 2 ≥ 70 % are expected; since the selection of parameters out of ranges provided in expressions (22)-(24) cannot guarantee such R 2 value, we developed the tuning method for the given ranges, which were selected as in [9,49].
Castellanos and Castrillón [57] presented an expression that allows one to compute the limit value that can be assigned to the controller's parameter in order to get a fast response, without compromising stability. Additionally, they reported a typical polynomial-like correlation that allows one to obtain the value of the tuning parameter that minimizes only the IAE. Such expressions are given by [57] whereτ cult is the minimum allowable value of the tuning parameter; i.e., lower values ofτ cult make the system exhibit unstable behavior; andτ cI AE is the value of the tuning parameter that minimizes the IAE.
In the present work we develop the equations for γ > 0 which allow one to use the combined index.

Tuning Parameter Computation
It is desired to obtain a new expression in order to relate the controller's tuning parameter that minimizes the objective function and the ultimate value that can be chosen in order to obtain a stable response:τ A simulation environment was developed using MATLAB R /Simulink R . This computational tool allows one to find the value of the objective function for each data set of parameters of the process' transfer function and γ. Since we used a central composite experimental design for rotatability, the optimization tool correspond to the steepest-descent gradient method [58].
The starting value ofτ c for each simulation was obtained from (25); the minimum assigned value for this parameter matched the limit value to obtain a stable system. The simulations considered unit-step changes in the disturbance variable (disturbance rejection). At the end of the simulation process, we obtained data shown in Table A1 (Appendix A) that contains the value ofτ c that minimizes OF for each data set of the process' transfer function (from the experimental design) and values of 0 ≤ γ ≤ 8. We found that no significant changes inτ cOF were obtained for γ > 4; therefore, the range 0 ≤ γ ≤ 4 was selected and Figure 2 showsτ cOF for different values ofτ cult considering such range for γ. Using the information provided in Figure 2, data were fitted using a linear regression, and the polynomial represents a straight line passing through the origin. The correlations obtained are given by Equations (28) to (32) are of the formτ where m OF is the slope of the line that relatesτ cOF andτ cult ; m OF addtionally represents a stability margin for the system. When γ increases, the value of the slope increases, andτ cOF increases, which make the system more robust. This characteristic allows getting a safety factor for the closed-loop system with respect to an ultimate valueτ cult ; if this value is reached, the time response would acquire sustained oscillations (marginal stability).
In order to measure the goodness of fit, some statistics need to be computed: the correlation coefficient (R), that measures the linear correlation between two variables; the coefficient of determination (R 2 ), that indicates the proportion of the variance that can be explained from the obtained model; and Durbin-Watson, that is used to check if there is a linear autocorrelation among the residuals of the model. Table 1 shows the statistics results; it can be noticed that values greater than 0.98 were obtained for R, and values greater than 97 % were obtained for R 2 ; this shows a strong goodness of fit. Values between 2 and 2.5 were obtained for Durbin-Watson, which shows that there is no linear autocorrelation among the model's residuals. It is necessary to obtain an expression that relates the value of m OF with γ. This is done using the value m OFi found in the linear Equations (28)- (32) for different values of γ. So, for 0 ≤ γ ≤ 4 we get the plot shown in Figure 3.
Using a nonlinear regression for the plot in Figure 3 we get the expression m OF = f (γ), as follows The best fit was given by a sigmoid function; the goodness-of-fit is measured again by using R and R 2 , both with values greater than 0.99 and 99%, respectively, as shown in Table 2.   Figure 3).
Equation (34) shows the relation between the stability margin m OF and the suppression factor γ. Then, an expression that allows computing the nondimensional optimal value of the controller's tuning parameter in a direct way is given bŷ Finally, Equations (10)-(12) are used to compute the PID controller's parameters using τ cOF , with units of time, that is computed as

Application Methodology
The application of the Castellanos-Castrillón-Vásquez (CCV) PID tuning method for inverse-response second-order plus dead time systems can be summarized in the following ten steps; Figure 4 shows the methodology in the algorithm form.

1.
Determine the transfer function of the inverse response second order plus dead time system and write it in the form of (5).

2.
Verify that the transfer function's parameters are within the specified ranges, established in (22)- (24); this method is valid only for such ranges.

5.
Select the value of γ, depending on the user's needs. If γ = 0, OF = I AE.
Check if the controller meets the desired performance.
If the desired performance is not met, the user can go back to Step 5 and change the value of γ adjusting the index in (13), which accounts for performance of the controlled variable and behaviour of the manipulated variable; high values γ are useful to protect the actuators, providing more conservative responses. Step 1 Step 2 Step 3 Step 4 Step 5 Step 6 Step 7 Step 8 Step 9 Step

Discussion on Performance Analysis
This section contains a comparison for the performance of a controller tuned with the new proposed CCV method and three other well-known techniques. There are several PI/PID controller tuning techniques [59], which are used for processes with specific dynamics characteristics. Regarding, inverse response processes, one can find that even recent literature still uses classical methods for comparison purposes or reference them as still valid and used for adjusting controller parameters. One can find recent literature in which new tuning methods, developed by Waller and Nygardas (WN) [13] and Chien, Chung, Chen, and Chuang (CCCC) [17] for inverse response systems, are compared to traditional methods; see for instance [60,61]. Additionally, one can find several works that compare the performance with one of the most traditional tuning methods developed by Ziegler and Nichols (ZN) [14] during the 1940s, see for instance [62][63][64]. The performance indexes are given by where %TO: percentage of transmitter output, %CO: percentage of controller output, UT: units of time.
A central composite experimental design, as explained in the scaling section, was implemented in order to select a set of parameters of the process' transfer function to test the performance and robustness of the closed-loop, for each controller, over a wide range of values. These values are within the defined ranges of the ratios given by (22)- (24). A tag P i was assigned to each set, and their values are shown in Table 3; these parameters are a subset extracted from the experimental region. For each simulation, the system was excited with a unit-step input for changes in the disturbance and changes in the set-point. A value of γ = 4 was assigned for all the simulations in order to penalize the behavior of the manipulated variable, considered in the second term of the optimization objective function (13). Table 4 contains the nondimensional controller parameters computed for each simulation set. In order to assess the general performance of the control system (disturbance rejection and set-point tracking), with a PID controller tuned with the new and the traditional selected methods, a scoring system was established. For each performance index, a score between 0 and 4 was given to each tuning method; the controller with the best performance received 4 points, the next received 3 points, and so forth. When a tuning method generated a marginally stable or unstable response, it received 0 points in all indexes. Table 5 contains the results of the performance indexes for each tuning method, when considering disturbance rejection for the set P 4 of parameters of the process' transfer function. Then, scores were assigned for each controller, when analyzing each performance index; Table 6 contains the disturbance rejection scores for the set P 4 .  All units are points.
Although the controllers were adjusted for disturbance rejection, the set-point tracking capabilities were also tested. Table 7 shows the results of the performance indexes for each tuning method, when considering set-point tracking for the set P 4 of parameters of the process' transfer function. Then, scores were assigned for each controller when analyzing each performance index; Table 8 contains the set-point tracking scores for the set P 4 .  All units are points. Figures 5 and 6 show the time response for the parameter set P 4 , for both disturbance rejection and set-point tracking, respectively.
It can be noticed that the CCV controller has a high settling time; this occurs because of the high value assigned to the suppression factor (γ = 4), in order to penalize the manipulated variable and increase the robustness, smooths the actuator's behavior. Additionally, it can be noticed how the CCV controller decreases the inverse response effect. This procedure was followed for the other parameter sets. In Appendix B, Tables A2-A9 contain the performance indexes for disturbance rejection and set-point tracking for the remaining sets of the parameters P 1 , P 2 , P 3 , and P 5 . In Appendix C, Figures A1-A8 show the time response for the remaining sets of parameters: P 1 , P 2 , P 3 , and P 5 .
All the points obtained by each tuning method using the scoring system were added, by considering each performance index. Table 9 shows the total scores for disturbance rejection and Table 10 shows the total scores for set-point tracking. Table 9. Total scores for disturbance rejection for the set P 4 . Total   CCV  17  9  19  17  62  WN  14  16  12  15  57  ZN  6  6  5  6  23  CCCC  13  19  14  12  58 All units are points. As it can be noticed, the proposed CCV method shows the best over-all results for disturbance rejection, when considering a wide range of parameters of the process' transfer function. Although the controller was tuned for disturbance rejection, it still offers a very good performance for set-point tracking, when considering the same wide range of parameters of the process' transfer function. In both operations (servo and regulatory), the CCV method achieves the best performance in both the MP and I MV indexes, which indicates that this tuning method benefits the durability of the actuator.

Method ISE I AE I MV MP
Although IMC is based on a pole-zero cancellation, which can lead to poor performance in the load regulation operation [65] for lag-dominant processes [66], in Appendix D we provide another simulation for a set of parameters with small dead time and high time constant (η = 0.1, θ = 0.01, τ 2 = 0.9, τ 1 = 1) that shows the wide application of the CCV method; this set is in the limits of the experimental region, but that was not given within the sets selected by the experimental design.

Discussion on Robustness Analysis
Romagnoli and Palazoglu [6] stated that there are different techniques to evaluate uncertainties in a control system. One of them consists in the selection of a range of the parameters of the process' transfer function, as in [67], in which Arbogast et al. proposed a robust stability factor metric to examine the effect of plant-model mismatch in the process gain, dead time, and time constant for self-regulating processes. Such method uses the relation of the values of parameters that bring the system's response to marginally stable conditions and the values of the parameters used for tuning.
Following such technique, in this work the robustness analysis consists on the study of the variability of parametersη andθ from Table 3 by considering a value for the suppression factor γ = 4. The analysis was done using a Simulink scenario in which small increments for the parameter of interest were done while keeping constant all the other parameters in the plant and the controller. Simulations stopped when the output reached a sustained-oscillation response (marginal stability conditions), hence, the ultimate value was found for each parameter. Results are found in Table 11. The table indicates the value of the parameter that leads the response to reach sustained oscillations for each controller and for a particular set of parameters. Therefore, the higher the ultimate value the more robust is the technique with respect to changes in the parameter that is being analyzed.
The CCV controller is more robust with respect to variations inθ, specially in the case when the dominant time constant and dead time are almost equal, or even equal. The ratio between dead time and the dominant time constantθ = θ τ 1 , is known as the uncontrollability parameter and when it is near 1, it is more difficult to control the process [68]. Therefore, the CCV tuning method becomes a technique suitable for high values of such parameter. On the contrary, for low values ofθ, the CCV controller does not achieve the better results. With respect toη, the CCV and CCCC controllers are more robust; both with similar results. However, the latter gets better results for low values ofθ. Table 11. Robustness for sets P 1 , P 2 , P 3 , P 4 , P 5 .

Method
Set

Conclusions
In this work, we proposed an IMC-based PID tuning method for inverse-response second-order plus dead time systems. The tuning rules are based on the optimization of an objective function that combines performance and robustness. The tuning method has been presented by using an easy-to-follow ten-step methodology with equations that are mathematically simple.
A correlation that allows one to compute the value of the tuning parameter τ c that minimizes the objective function has been found. The tuning parameter τ c , affects the stability of the closed-loop control system. Small values of τ c increase the speed response of the system, but also produce an oscillatory response, to the point that the system can become unstable. Nondimensionalization reduced the number of parameters, which allows the reduction of simulations runs, saving computation time. The central composite experimental design allowed the authors to determine the appropriate number of simulations, and obtain goodness of fit for the proposed model.
The performance of a PID controller, tuned with the proposed CCV and other literature-existing tuning rules, was evaluated. The performance achieved with the proposed CCV method was excellent, not only for disturbance rejection but for set-point tracking, when considering a wide-range of parameters of the process' transfer function. For both operations (servo and regulatory), the CCV method achieves the best performance in both the MP and IMV indexes, which indicates that this tuning method benefits the durability of the actuator. This is a very important result regarding continuous plant operation in industrial processes, since it can help avoiding unexpected plant stops caused by actuator failures.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix B. Performance Results for Remaining Sets of Parameters
Tables A2-A9 contain performance results for remaining sets of parameters.

Appendix D. Additional Simulation in the Limits of the Experimental Region
We provide another simulation for a set of parameters with small dead time and high time constant that shows the wide application of the CCV method, see Table A10; this set is in the limits of the experimental region, but that was not given within the sets selected by the experimental design. Table A10. Simulation parameters set for lag-dominant simulation (with τ 1 = 1).

Parametersτ 2ηθ
P 6 0.900 0.100 0.01 The system was excited with a unit-step input for changes in the disturbance and changes in the set-point. A value of γ = 4 was assigned the simulation in order to penalize the behavior of the manipulated variable (as in all previous simulations), considered in the second term of the optimization objective function (13). Table A11 contains the nondimensional controller parameters computed for this simulation set. Tables A12 and A13 contain the performance indexes for the set of parameters P 6 . It can be noticed that the CCV controller offers the lower index for the work of the manipulated value, since we chose the suppression factor (γ = 4), to smooth the actuator's behavior and increase the robustness. Additionally, the CCV controller response do not exhibit overshoot and show the smallest inverse-response effect, see Figures A9 and A10. This proves the wide application of the CCV tuning method.