Parametric and Nonparametric PID Controller Tuning Method for Integrating Processes Based on Magnitude Optimum

: Integrating systems are frequently encountered in power plants, paper-production plants, storage tanks, distillation columns, chemical reactors, and the oil industry. Due to the open-loop instability that leads to an unbounded output from a bounded input, the efficient control of integrating systems remains a challenging task. Many researchers have addressed the control of integrating processes: Some solutions are based on a single closed-loop controller, while others employ more complex control structures. However, it is difficult to find one solution requiring only a simple tuning procedure for the process. This is the advantage of the magnitude optimum multiple integration (MOMI) tuning method. In this paper, we propose an extension of the MOMI tuning method for integrating processes, controlled with a two-degrees-of-freedom (2-DOF) proportional– integral–derivative (PID) controller. This extension allows for calculations of the controller parameters from either time domain measurements or from a process transfer function of an arbitrary order with a time-delay, when both approaches are exactly equivalent. The user has the option to emphasise disturbance-rejection or tracking with the reference weighting factor b or apply two different reference filters for the best overall response. The proposed extension was also compared to other tuning methods for the control of integrating processes and tested on a charge-amplifier drift-compensation system. All closed-loop responses were relatively fast and stable, all in accordance with the magnitude optimum criteria.


Introduction
This paper is an extension of a previous study [1] that proposed tuning of the proportional-integral (PI) controller for integrating processes. A proportional-integral-derivative (PID) controller, in single or cascade loops and with or without derivative/integral terms, is the most commonly used control algorithm in industry [2][3][4], although several more advanced control structures have been developed. According to the authors in Reference [5], "Based on a survey of over eleven thousand controllers in the refining, chemicals and pulp and paper industries, 97% of regulatory controllers utilise a PID feedback control algorithm". The disturbance-rejection and tracking performance depend on tuning the PID controller parameters. Despite the presence of auto-tuning algorithms [6,7], PID controllers in practice are inadequately tuned due to a lack of skill and time. According to Reference [8], "a visit rejection. A PID controller with an additional zero for the control of second-order IPs was proposed in Reference [113]. The controller tuning rules also included the controller's sampling time, Ts. In Reference [109], a new classification for a large class of stable, oscillatory, integrating, and unstable processes with or without time-delay was proposed. This was done by using the proposed ρ-ϕ parameter plane, defined by the normalised gain, ρ, and the angle, ϕ, of the Nyquist curve slope at the process' ultimate frequency. Next, PID controller tuning formulae for obtaining the desired performance/robustness trade-off in the chosen region of the ρ-ϕ parameter plane are proposed.
Nonparametric tuning methods are data-based, i.e., they directly exploit the closed-or the open-loop process datasets, without the need for a process model. This is an attractive approach for alleviating the drawbacks of process-model mismatches. A PI/PID controller tuning method that enables specification of the desired closed-loop transfer function for disturbance-rejection, while tracking, using PID controllers, can be independently improved with set-point weighting, was proposed in Reference [115]. A tuning method for PI and PID controllers, based on an extended relay feedback procedure, was proposed in Reference [117].
Most of the proposed tuning methods offer limited range of IP models, require an accurate IP model, or optimisation of the controller parameters. Additionally, some methods do not give satisfactory closed-loop responses, and some others may not work with parameter uncertainty (stability/robustness).
Magnitude optimum (MO) is a controller tuning rule based on the optimal closed-loop frequency response [118][119][120]. The closed-loop performance for tracking and disturbance-rejection is usually fast and without oscillations [121]. By using a nonparametric time domain approach with multiple integrations of the process signals, the applicability of the MO method has been extended. This extension is called the magnitude optimum multiple integration (MOMI) method [122] and retains all the advantages of the MO method. Additionally, this extension offers an extraction of the necessary data from a simple time domain experiment [121,122].
When using a one-degree-of-freedom (1-DOF) controller structure, the MO criteria result in the integrating gain of the PID controller equal to zero. Since the resulting PD controller cannot reject input disturbances, the original MOMI tuning rule cannot be used for IPs. However, in Reference [1], Kos et al. proposed an extension of the MOMI tuning method for IPs controlled by a 2-DOF PI controller.
In this paper, an extension of the MOMI tuning method for integrating processes, controlled with the 2-DOF PID controller, is proposed. Additionally, to achieve the best overall closed-loop responses (the optimal tracking and disturbance-rejection), two reference filter structures are proposed. The proposed method is nonparametric, meaning that it does not require an explicit process model, which is an attractive approach for alleviating the possible problems related to the process-model mismatch. The controller parameters can be calculated from the process' closed-or open-loop time response or from the arbitrary process transfer function with a time-delay. Both approaches are equivalent. Furthermore, the user has the additional option to emphasise disturbance-rejection or tracking with the reference weighting factor b or to simply use the reference filter for the best overall responses.
The content of this paper is organised as follows. First, the MOMI tuning method and the PID extension for IPs are presented. Section 3 presents two reference filter structures for achieving the best overall closed-loop responses. Section 4 provides examples of diverse IP models and a comparison between the proposed MOMI PID and MOMI PI extension for the IPs and between the proposed extension with and without the reference filter. Section 5 assesses the stability and robustness of the proposed tuning method. In Section 6, the proposed method is compared with other IPs tuning methods. A real-time experiment in the time domain on a charge-amplifier drift-compensation system is outlined in Section 7.

Closed-Loop Configuration and Magnitude Optimum Multiple Integration
In this section, we present a method for the non-parametric identification of the process using so-called characteristic areas (process moments) [18]. These areas can be identified from a closed-or open-loop time response (in the time domain) or directly from an arbitrary rational transfer function with a pure time-delay (in the frequency domain). It is worth mentioning that both approaches (time and frequency domain) are equivalent. Later, the characteristic areas will be used for the calculation of the controller and the reference filters' parameters. Figure 1 shows a 2-DOF controller in a closed-loop configuration with a process. Signals u, d, y, and r represent the process output, input disturbance, controller output, and controller reference, respectively. The process is defined by the stable transfer function: G P (s) = K PR s 1 + b 1 s + b 2 s 2 + · · · + b m s m 1 + a 1 s + a 2 s 2 + · · · + a n s n e −sT delay , where the process time-delay is represented by T delay . As stated in Reference [119], "one possible design aim is to maintain the closed-loop magnitude response curve as flat and as close to unity for as large a bandwidth as possible" (see Figure 2). Thus, the goal is to find a controller that has a low-frequency closed-loop response as close to unity as possible [121].  This approach is called modulus optimum, Betragsoptimum, or magnitude optimum (MO), and results in a relatively optimal closed-loop time response for a majority of process models, when considering the speed of its response and robustness [121].
With a 2-DOF controller (in Figure 1), the closed-loop transfer function becomes The controller is configured so that is fulfilled for as many k as possible [119,122]. Equation (3) is thus easily fulfilled. If the controller structure contains integral term and the closed-loop response is stable, the steady-state control error is zero; hence, Equation (3) holds. The number of satisfied conditions in Equation (4) is correlated to the number of controller parameters (controller order). For the PID controller, k max = 3. Let us define the closed-loop transfer function by the following equation: Then, the following conditions must be satisfied to fulfil Expression (4) [121]: For the calculation of the PI controller parameters, the first two equations in Condition (6) should be fulfilled, and for the calculation of the PID controller parameters, the first three equations in Condition (6) should be fulfilled.
To simplify further derivations, Process (1) is developed into an infinite Taylor series [122]: where A i represents the "characteristic areas" of the process [121], which can be expressed by Process Parameters (1) [119]: This means that the characteristic areas can be obtained directly from the process transfer function. Additionally, the characteristic areas can also be calculated in the time domain, i.e., from the process time response. To this end, the steady state of the process should be changed first. Then, the multiple integrations of the process input (u(t)) and output (y(t)) signals can be calculated [119,122]: where u 0 represents the normalised process input signal, and The characteristic areas are expressed as follows: Note that ..
.. y (0) = · · · = 0 and that A 0 equals the process steady-state gain, K PR . In practice, the integration in Expression (9) can end when the process signals settle. Thus, Expression (9) can be recursively and simply calculated, therefore the process Model (1) is not needed. Therefore, using the characteristic areas, the process can be easily parameterised from any change in the process steady state [121].
In practice, certain rules need to be followed to accurately calculate the areas from the process response corrupted by noise. A few practical guidelines for calculating the characteristic areas from the process time responses in practice, including the high-frequency process noise, are provided in Reference [123]. As will be shown in Section 7, the method was also tested on a real process with present noise. Moreover, the proposed method for the PI controller [1] has also been tested on several processes with measurement noise. In all the mentioned tests, it was shown that the calculated areas and, therefore, the controller parameters, are not extensively sensitive on high-frequency process noise. However, the method still requires deeper noise analysis to be carried out in the future.
For ease of clarification, Figure 3 shows a graphical representation of the characteristic area A 1 . Signals u 0 and y 0 represent the process input and process output responses, respectively [121]. Figure 3. Characteristic area A 1 (graphical representation). This area was measured after the change of the process steady state. Figure 4 presents the closed-loop system with a 2-DOF PID controller. Signals y, d, u, and r represent the process output, disturbance, controller output, and controller reference, respectively. Parameters b and c are the proportional and derivative weighting factors, while K D , K P , K I , and T F are the derivative gain, proportional gain, integral gain, and filter time constant, respectively. The controller presented in Figure 4 is defined by the following transfer functions:

Extension to Integrating Processes
To simplify the derivations, c = b is chosen. Thereafter, this controller will be denoted PID b . While the controller filter is actually part of the PID b controller transfer function. To simplify the derivation of the PID b controller parameters, the controller filter is considered part of the Process (1). Therefore, the aforementioned filter must be added to Expression (1) every time the PID b controller parameters are calculated with the final formula. The controller filter is used for reduction of the controller's output high-frequency noise. The controller's high-frequency gain equals to K HF = K D T F [18]. Besides defining the controller filter time constant in advance, it can also be calculated based on the desired high-frequency gain as follows: Initially, the controller parameters (including derivative gain K D ) should be calculated with T F = 0. Next, the controller parameters can be calculated again by applying Equation (14). The procedure can be repeated until the T F value settles. In practice, up to three iterations are usually sufficient.
In the following derivation, the PID b controller parameters will be calculated. The IP has been modelled by the transfer Function (1). Next, Process (7) and the Controller (12) transfer functions are inserted into a closed-loop transfer Function (2). In this way, the following closed-loop transfer Function (5) parameters are obtained: To simplify the derivation of the PID b controller parameters, in the next step, all closed-loop transfer function parameters, represented by Expressions (15) and (16), are divided by A 0 K I . Consequently, the following closed-loop transfer function Parameters (5) are obtained: where For calculation of the three PID b controller parameters, the first three equations in Condition (6) should be solved. The first condition of Expression (6) gives the following result: and the second condition in Expression (6) provides the following result: The third condition of Expression (6) gives the following fourth-order polynomial: where To calculate the PID b parameters, the integral time constant T I , represented by Expressions (22) and (23), should be solved first. The correct result is the highest positive real number. In the next step, integral gain, K I , should be calculated using Equation (21). In the last step, the proportional, K P , and derivative, K D , gains can be calculated using the following expression: Since calculating the PID b controller parameters requires solving the fourth-order polynomial, the final derivation of the tuning formula is not presented. Regardless, the tuning formula remains analytical. To facilitate the implementation of the proposed tuning method and save the reader's time and effort, all MATLAB files used for the calculation of the PID b parameters are available online in Reference [124].

Remark 1.
Weighting factor b influences the disturbance-rejection and tracking performance. When b approaches 1, the integral gain K I (21) approaches 0 [1]. For b = 1, we obtained the following PID controller: Thus, we obtained a PD controller that cannot reject the process input disturbances. For this reason, increasing the value of weighting factor b degrades the disturbance-rejection performance. In practice, the value of b should be limited within the interval 0 < b < 0.9 [1]. An illustrative example showing the tracking and disturbance rejection performance and its relationship to factor b is presented in Section 2.3. (21) and (24) result in stable and fast closed-loop responses for most IP models. However, the closed-loop stability of an arbitrary process model is still not guaranteed. Section 5 provides additional information on closed-loop stability and robustness.

Remark 2. The PID b controller parameters calculated by Expressions
As explained in Section 2.1, the characteristic areas A 0 to A 4 can be obtained directly from the process transfer function or calculated from the process time response when changing the process steady state. The following steps are required for PID b controller tuning:

1.
Choose the appropriate filter time constant, T F , of the PID b controller.

2.
Determine the characteristic areas A 0 -A 4 . If the process is given in a Laplace form, use Expression (8) (frequency domain). Otherwise, change the steady-state of the process, capture the controller and process output signals, and calculate the characteristic areas using Expressions (9)-(11) (time domain). The initial values of the controller and the process output signals (before the first change of the process input signal) can be approximated from the measurements. Do not forget to add the filter time constant, T F , to the process transfer Function (1) before calculating the areas. If the time domain approach is used, use the first-order filter with a time constant that equals T F to filter the process output signal.

3.
Select the reference weighting factor b. As seen in Remark 1, the recommended values are 0 ≤ b < 0.9.

4.
Use Expressions (21) and (24) to calculate the PID b parameters K P , K I , and K D .

Remark 3.
For obtaining the a-priori chosen high-frequency controller gain K HF , the exact value of T F can be calculated according to the procedure described below Equation (14).

Illustrative Example
The following second-order delayed IP is chosen to illustrate the proper design: First, the PID b controller filter time constant T F = 0.01 s is added to Process (26). Using Expression (8), the following characteristic areas are calculated from Process (26)'s transfer-function with a filter: In the next step, using Expressions (21) and (24), the PID b controller parameters are calculated for different values of reference weighting factor b. The calculated parameters presented in Table 1 are applied to the PID b controller. The closed-loop responses (reference step change r = 1 and input disturbance step change d = 1) are shown in Figure 5.  (26)). Note that, for clarity, the controller output response is zoomed in for the reference step change. As seen in Figure 5, the reference weighting factor b has a significant impact on the closed-loop responses. Tracking (reference-following) performance is improved at higher values of b. Thus, in cases where good tracking is desired, choose a value of b ≥ 0.5. However, high values of b degrade the disturbance-rejection performance, i.e., the disturbance-rejection performance seriously deteriorates at b > 0.7. The most effective disturbance-rejection is obtained with b ≤ 0.5. An acceptable compromise between both tracking and disturbance-rejection seems to be b = 0.5.

Optimal Tracking and Disturbance-Rejection
Besides using the reference-weighting parameter b (note that c = b), the closed-loop tracking response can be improved by applying a reference filter. The main advantage of using a reference filter is that such a filter can substantially improve the tracking response without deterioration of the disturbance-rejection response (the optimal one, as calculated with b = 0). In this way, the best overall response can be achieved. A similar approach was proposed in Reference [125] for non-integrating processes. The filter (G F ) in the closed-loop configuration with the PID controller (G C ) and the process (G P ) is shown in Figure 6. For the proposed tuning method with a reference filter, the controller parameters are calculated for the weighting factor b = 0 to retain the good disturbance rejection properties, while the controller structure is 1-DOF PID (b = 1).
controller process Figure 6. The reference filter with the PID controller and the process.
Here, two controller structures with reference filters are proposed according to the following filter order: The PID controller with the second-order filter (PID f2 ) and the higher-order filter (PID fh ).

The Second-Order Filter
The closed-loop transfer function G CL from the reference to the process output is the following: The purpose of the reference filter G F is to make the reference following transfer Function (28) as closely as possible to the desired function with the PD controller (b = 1): where K P1 and K D1 represent the controller Parameters (25) of the PD controller. The structure of reference filter G F is chosen as follows: The denominator of the Filter (30) is chosen to cancel the controller zeros in the closed-loop transfer Function (28), similar to the approach proposed by Medarametla and Manimozhi [106]. The numerator of the filter (parameters b F1 and b F2 ) is then calculated by equating the first two characteristic areas (A 0 and A 1 ) of the closed-loop transfer Function (28) (using Expression (8) with G CL instead of G P1 in the Expression (7)) with the equivalent areas of the desired closed-loop transfer function G CL1 when using the PD Controller (29). The following filter parameters are thus obtained:

The Higher-Order Filter
To achieve an even faster closed-loop tracking response than that with the PD Controller (25), a higher-order filter can be applied. Here, the calculation of the filter parameters is based on the process model. However, this does not mean that the main advantage of the proposed tuning method is lost (note that the process model is not required for the calculation of the controller parameters). Instead, the process model will be estimated from the characteristic areas, where the only required parameter from the user is the estimated process time-delay, which can be easily acquired in practice from the process time-response. The chosen process model is a second-order process with zero and the estimated pure time-delay where T dm and K M represent the process time-delay and steady-state gain, respectively. Note that the model parameters (except for the estimated time-delay, T dm ) can be calculated directly from the equations for the characteristic areas A 0 to A 3 (8): For minimum-phase processes, the model transfer function can be simplified by setting b 1m = 0, which in turn simplifies the identification procedure. Namely, all the process model parameters, including the process time-delay T dm , can be calculated analytically [126] from Expression (8): where T dm is calculated from the above third-order equation. The result (among three solutions) is the lowest positive real number.
The closed-loop transfer function from the reference to the process output, when process G P is replaced by the process model G M , is the following: We can now define the desired closed-loop transfer function as follows: where T CL is the desired closed-loop time constant. Note that parameter T CL is related to the speed of the response of the PD controller, so it will be automatically calculated, and no user input is required. Let us define the closed-loop time constant T CL0 equal to the integral of the control error at the unit stepwise reference change: where K P1 is the proportional gain of the PD controller. The desired closed-loop time constant T CL is then defined as follows: where parameter K CL defines the speed-up coefficient of the closed-loop time constant in relation to that obtained with the PD controller. According to the results of the experiments using different process models, the parameter K CL = 3 is found to be a good compromise between tracking speed and robustness.
To equate Expressions (35) and (36), the time-delay in the denominator of (35) is expressed by the second Pade's approximant: By Equations (35) and (36) and taking into account (40) in the denominator of (35), the following reference filter transfer function is obtained: where n is calculated from Expression (37). Parameters b F1 to b F5 in the numerator are calculated from the denominator of the closed-loop transfer Function (35) without the reference filter G F : by equating the terms with the same power of the Laplace variable s. The resulting parameters are the following: where T F is the controller filter time constant.

Remark 4.
The reference filter is not in the feedback part of the closed-loop system transfer function. Therefore, as long as the reference filter is stable, it does not change the stability of the entire closed-loop system. By observing the higher-order reference filter Denominator (41), it can be seen that majority of the filter poles are on the left-hand side of the complex plane. The only part of the denominator filter that might become unstable is the following term: By using simple derivation, it can be shown that the mentioned Term (44) is stable if all of the controller gains are of the same sign. The same applies for the second-order reference Filter (30), which also contains the same Term (44).
To facilitate the implementation of the proposed reference filters and to save the reader's time and effort, all MATLAB files for the calculation of the second-and higher-order filter parameters are available online in Reference [124]. Note that, as with the MOMI tuning method, the reference filter parameters can be calculated from the measured characteristic areas, which can be acquired from the process steady-state time response.

Illustrative Examples
The proposed method (PID b , PID f2 , and PID fh ) was tested on some process models. Additionally, it was compared to the MOMI PI controller tuning method for integrating processes [1].

MOMI PID b Controller
In this section, the proposed method without a reference filter has been tested on some common process models. The chosen IPs were a delayed second-order process, a pure integrator plus time-delay process, a fourth-order process, and a non-minimum phase process. These processes are represented by the following transfer functions: First, the a-priori chosen PID controller filter time constant T F = 0.01 s was added to Processes (45). Then, using Expression (8), the characteristic areas A 0 to A 4 were calculated. The values are presented in Table 2. The PID b controller parameters for different values of b (b = 0, b = 0.5, and b = 0.7) were calculated from Equations (21) and (24) using the calculated characteristic areas from Table 2. The calculated parameters for all selected IP models are presented in Table 3. The calculated controller parameters were applied to a PID b controller. Figures 7-10 show the closed-loop responses for the input disturbance step change (d = 1) and for the reference step change (r = 1). From the experimental results, it can be seen that the proposed tuning method provided relatively fast and stable closed-loop responses. Additionally, all responses were highly damped with only minor overshoots. The closed-loop responses are in accordance with the magnitude optimum criterion. Exceptions include the sluggish disturbance-rejection responses at b = 0.7. As expected, the best tracking performance was achieved with b = 0.7, and the best disturbance-rejection performance was achieved with b = 0 (see also the illustrative example in Section 2.3). The Bode plots of the closed-loop amplitude gains (between the reference and the process output) are shown in Figure 11. This figure shows that all the processes have a flat amplitude low-frequency response, which is in accordance with the MO criterion (see Figure 2). Moreover, the reference weighting factor b has an influence on the closed-loop bandwidth, i.e., the closed-loop bandwidth decreases by decreasing factor b. Therefore, the tracking performance improves by increasing factor b, which is in accordance with the findings in Section 2.3 (illustrative example).

Comparison with the MOMI Tuning Method for PI Controllers
The performance of the proposed MOMI PID b method was compared with the MOMI PI method for IPs [1] on four process models. Note that this PI controller also includes the reference weighting factor b on the proportional term. The following IPs were chosen: First, the parameters for the PI controller were acquired. To that end, the characteristic areas A 0 -A 2 , presented in Table 2, were calculated. The PI controller parameters for two different values of b (b = 0 and b = 0.5) were then calculated via the following Equations [1]: using the calculated characteristic areas from Table 4. The calculated PI controller parameters for all four IP models are presented in Table 5. Next, the parameters for the PID b controller were acquired. Note that the a priori chosen PID b controller filter time constant T F = 0.01 s (13) was added to Processes (46). Then, by using Expression (8), the characteristic areas A 0 to A 4 presented in Table 6 were calculated. The PID b controller parameters for different values of b (b = 0 and b = 0.5) were calculated from Equations (21) and (24) using the calculated characteristic areas from Table 6. The calculated PID b controller parameters for all four IP models are presented in Table 7. The calculated controller parameters presented in Tables 5 and 7 were applied to the PI and PID b controllers. Figures 12-15 show the compared closed-loop responses for the input disturbance step change (d = 1) and for the reference step change (r = 1). This comparison shows the superior tracking and disturbance-rejection performance of the PID b controller.  Figure 13. The PI and the PID b controller closed-loop responses for the process G P2 . Note that, for clarity, the controller output response is zoomed in for the reference step change.  Figure 15. The PI and the PID b controller closed-loop responses for the process G P4 . Note that, for clarity, the controller output response is zoomed in for the reference step change.

MOMI PID with Reference Filters
Besides using the reference weighting factor b, the 2-DOF controller can also be realised via a PID with an additional reference filter (see Section 3).
Here, the proposed PID b controller will be compared to the proposed PID f2 and PID fh controllers. For this purpose, the same process models used in Section 4.1 (45) are used here.
The chosen controller filter time constant was the same as before (T F = 0.01 s). The characteristic areas A 0 to A 4 are given in Table 2. The PID b controller parameters calculated using various values of b are given in Table 3. Note that, for the PID f2 and PID fh controllers, the proportional, integral, and derivative gains for PID b with b = 0 were used.
The reference filter parameters for the PID f2 and PID fh controllers were calculated from Equation (31) (PID f2 controller) and Equation (43) (PID fh controller) using the calculated characteristic areas from Table 2 and the calculated PID b controller parameters (b = 0) from Table 3. The values of the PD controller parameters, the estimated second-order process model parameters, and the chosen speed-up coefficient of the closed-loop time constant are shown in Table 8. The calculated reference filter parameters for the PID f2 and PID fh controllers are given in Tables 9 and 10, respectively.  Figures 16-19 show a comparison of the closed-loop responses for input disturbance step change (d = 1) and the reference step change (r = 1). This comparison shows the superior tracking abilities of the PID f2 and PID fh controllers. The tracking speed of the PID fh controller is considerably greater than that of the other methods. Note that the disturbance-rejection performance of the PID f2 and PID fh controllers is equal to that of the PID b controller with b = 0. On this basis, we can conclude that the proposed tuning method for the PID f2 and PID fh controller gives the best overall response (optimal tracking and disturbance-rejection).  Figure 19. The comparison of process G P4 and the controller output responses using the PID b , PID f2 , and PID fh controllers.

Stability and Robustness
In general, the stability of the closed-loop system can be calculated from the roots of the closed-loop transfer function denominator (den CL ): where Clearly, the denominator consists of a pure time-delay, which complicates the stability analysis. The time-delay in den CL (49) can be developed into an infinite-order Taylor series. However, the analysis of infinite-order series using, e.g., the Routh-Hurwitz stability criterion, is far from straightforward and outside the scope of our research.
However, the stability and robustness of closed-loop systems can be tested on process models that cover most of the integrating processes in various industries. One such process model is the second-order integrating process with a time-delay: The process G P (s) (50) also covers pure time-delay and the first-order integrating processes. The stability and robustness of the closed-loop system will be tested with different ratios of T delay /a 1 and factor α. Note that this process has real poles only when α ≤ 0.25.
When choosing the PID b controller reference weighting factor b = 0, the closed-loop system is stable for any positive values of α and T delay . The system is also stable for all the processes with real process poles for b = 0, b = 0.5, and b = 0.9. As shown by Figure 20, increasing factor b has a negative impact on system stability. This system is stable for any time-delay with b = 0.5 if α ≤ 1.5, and α ≤ 1.05 with b = 0.9. The overall stability improves by increasing the time-delay or decreasing b and the factor α.  The closed-loop system's robustness was tested by measuring the maximum sensitivity function, M S [1,11,127]. A higher value of M S indicates lower system robustness since the open-loop transfer function G P G CY is closer to the critical point (−1 + 0i). The most common values of M S for integrating processes are usually higher than those for non-integrating processes [11] and are typically between 1.7 and 3.
The maximum sensitivity functions are calculated for different values of b (0, 0.5, and 0.9) and are shown in Figures 21-23. Note that the M S values for processes with real poles are below 3.15 for b = 0, below 2.45 for b = 0.5, and below 1.85 for b = 0.9. The expected robustness of the closed-loop system for the processes with real poles is, therefore, satisfactory and relatively high. Processes with complex poles result in lower robustness, which is expected. The robustness of the proposed tuning method (PID b and PID fh controller) is next studied for the following process model: First, the PID b controller filter time constant T F = 0.01 s was added to Process (51). Using Expression (8), the following characteristic areas were calculated from Process (51)'s transfer function: In the next step, the PID b controller parameters for reference weighting factors b = 0.5 and b = 0 were calculated using Expressions (21) and (24): b = 0.5, K P = 0.46336, K I = 0.062992, K D = 0.5563, b = 0, K P = 0.51394, K I = 0.083446, K D = 0.58265. (53) Additionally, the parameters of the PID fh controller reference filter were calculated using Expression (43) The calculated parameters were applied to the PID b and PID fh controllers. The robustness was studied by assuming a ±10% change in the Process (51) steady-state gain and time-delay. Figures 24  and 25 show the corresponding responses for the input disturbance step change (d = 1) and the reference step change (r = 1) for the PID b controller with b = 0.5 and for the PID fh controller, respectively. As can be seen from the process output responses, a change in the process gain or time-delay does not deteriorate the tracking or disturbance-rejection responses. Nevertheless, the PID fh controller seems to be slightly more sensitive to process perturbations (especially the increased process gain).

Comparisons with Other Methods
As far as we know, no other tuning method for IPs offers tuning based on both the measurement of a simple process open-or closed-loop time response in the time domain and a general process transfer function with an arbitrary order with time-delay. This represents a significant advantage over alternative tuning methods for PID controllers.
Nevertheless, in this section, the results of the simulations of five different IP models that often appear in the literature are compared with those of other PID tuning methods for IPs. Notably, the selected tuning methods are not model-free controller tuning methods. Some methods cannot be used on arbitrary IP models. However, all the chosen methods provide analytical tuning formulas for the calculation of the controller parameters.
The proposed tuning method for the PID f2 and PID fh controllers was compared to five other methods [18,77,97,103,106]. The first method is Ali and Majhi's method [77] (hereafter denoted as Ali-Majhi) for 1-DOF PID controllers. This method achieves optimal load disturbance-rejection for IPs with Integral of the Squared Error (ISE) criterion minimisation with Nyquist curve constraints, e.g., the slope has a specified inclination at the gain crossover frequency. The second method is that of Taguchi and Araki [97] (hereafter denoted as Taguchi) for 2-DOF PID controllers. This method optimises the settling time while keeping the overshoot under 20%. The third one is Medarametla and Manimozhi's method [106] (hereafter denoted as Medarametla), which is based on the loop sensitivity transfer function. This controller consists of a 1-DOF PID controller with a second-order Lead/Lag filter and a fourth-order reference filter that cancels some controller zeros. The fourth method is that of Anil and Padma Sree [103] (hereafter denoted as Anil-Padma), which employs a pole-placement strategy, and the tuning parameter is the M S . The controller consists of a 2-DOF PID controller with a first-order Lead/Lag filter. The last is Åström and Hägglund's method [18]. This method is based on selecting a maximal sensitivity value M S = 1.4 (hereafter denoted as Åström 1.4) or M S = 2.0 (hereafter denoted as Åström 2). The method employs a 2-DOF PID controller. Note that these methods employ different controller structures.
The closed-loop tracking (r = 1) and disturbance-rejection (d = 1) responses were measured on all five process models. In all cases, the a priori chosen controller filter time constant was T F = 0.01, with the exception of the Ali-Majhi method, where the controller filter time constant was set to the proposed value of T F = T D × 0.1.
The calculated controller parameters (for all presented tuning methods and the proposed tuning method) for all five IP models are found in the Tables S1-S10 of the Supplementary Materials Section 1.
The tracking and disturbance-rejection performance of the tuning methods can be evaluated by the following criteria: where RVy is the relative variance of the process output signal: and y 1 and y 2 are the first and the second extrema of the process output signal. The RVy for tracking becomes 1 for a monotone process output signal. For any other response (overshoots, undershoots, or oscillations), the value will increase. The RVy for disturbance becomes 1 for the process output signal, which has only one maximum and one minimum. This is an ideal disturbance-rejection response.
Since the RVy measure favours monotonic responses, while IAE favours fast responses, the combined measure (RVy × IAE) offers the best compromise between speed and monotonicity. The RVy measures are based on the TV0y and TV2y measures introduced in References [48,128]. Note that the TV0y and TV2y measures are 0 for monotone responses (or responses with one maximum and one minimum). Since we multiply the RVy with IAE, the RVy is modified to give values of 1 instead of 0 under ideal circumstances.
Like RVy × IAE, the IAT 2 E measure favours fast and non-oscillatory responses.
Note that the measurements of the criteria were taken after the process output completely settled, meaning that the simulations were longer than those shown for the closed-loop responses below.

Case 1
The next pure integrator plus time-delay is chosen as The tracking and the disturbance-rejection closed-loop responses between the aforementioned tuning rules are compared in Figures 26 and 27. The criteria values of the tuning rules are presented in Figure 28.  Figure 26. The process G P1 closed-loop responses. The PID controllers were tuned by the PID f2 , PID fh , Taguchi [97], Anil-Padma [103], and Medarametla [106] methods.   According to the IAE criterion, the proposed controllers ranked first, and according to the IAT 2 E criterion, they ranked first and third in tracking performance. The proposed PID fh controller also ranked first for the RVy × IAE criterion. The Ali-Majhi method showed a large overshoot, and the Åström 2 method exhibited a somewhat oscillatory response.
The Taguchi method obtained the lowest criteria for disturbance-rejection. The proposed tuning method was ranked second. However, the Taguchi method exhibited a slightly oscillatory response.
Concerning sensitivity M S , Anil-Padma and Medarametla suggested M S = 2 for processes without zero, while the proposed method gave values below 3.3.

Case 2
The following slightly delayed first-order process is chosen:  Figure 29. The process G P2 closed-loop responses. The PID controllers were tuned by the PID f2 , PID fh , Taguchi [97], Anil-Padma [103], and Medarametla [106] methods.   The lowest criteria for tracking was obtained by the proposed PID fh controller. The Medarametla method ranked second, and the PID f2 ranked third. The Taguchi and the Ali-Majhi method exhibited a somewhat oscillatory response.
For the disturbance-rejection performance, the criteria favoured the Medarametla method. However, the proposed method was ranked second, close to Medarametla. Like tracking, the Taguchi and Ali-Majhi methods exhibited oscillatory responses.
Concerning sensitivity M S , Anil-Padma and Medarametla suggested M S = 2 for processes without zero, while the proposed method gave values below 3.

Case 3
The next second-order process with time-delay is chosen: The tracking and disturbance-rejection closed-loop responses between the tuning rules are compared in Figure 32. Some methods are missing because they are not suitable for second-order process models. The criteria values of the tested tuning methods are presented in Figure 33.   The tracking criteria gave the lowest values for the proposed controllers. The Åström methods ranked second. Again, the Taguchi method exhibited an oscillatory response.
Similarly, the lowest criterion for disturbance-rejection was given for the proposed controllers, with the Åström 2 method close behind.
Concerning sensitivity M S , the proposed method gave values below 2.5. As expected, the M S values for the Åström methods were lower since that tuning method is based on M S values.

Case 4
The following high-order non-minimum phase process is chosen: An experiment using this process was already carried out by Medarametla and Manimozhi [106] (Medarametla), and Anil and Padma Sree [103] (Anil-Padma).
The tracking and the disturbance-rejection closed-loop responses between the aforementioned tuning rules are compared in Figures 34 and 35. Some methods are missing because they are not suitable for higher-order process models. The criteria values are presented in Figure 36.   The lowest IAE and the IAT 2 E criteria for tracking were obtained with the PID fh controller. The Medarametla method was ranked second, and the proposed PID f2 controller ranked third. Similarly, the RVy × IAE criterion favoured the proposed PID fh controller, the Anil-Padma method was ranked second, and the Medarametla method ranked third.
Similar results were obtained for the disturbance-rejection performance. The proposed PID f2 and PID fh controllers were ranked first, and the Medarametla method was ranked second.
Concerning sensitivity M S , Anil-Padma and Medarametla suggested M S = 2.81 for this process, while the proposed method gave values below 3.2. As expected, the M S values under the Åström methods were lower since the tuning method is based on the M S values.

Case 5
The next high-order process is chosen: The tracking and the disturbance-rejection closed-loop responses between the chosen tuning rules are compared in Figure 37. The criteria values are presented in Figure 38.  Figure 37. The process G P5 closed-loop responses. The PID controllers were tuned by the PID f2 , PID fh , Åström   Only the proposed and Åström methods were appropriate for the high-order process. According to the criteria, the proposed controllers (PID f2 and PID fh ) ranked first in tracking and disturbance-rejection responses. As expected, the M S values under the Åström methods were lower since the tuning method is based on the M S values.
All experiment results showed that the proposed PID f2 and PID fh controllers, compared to some other tuning methods based on similar tuning procedures, result in very good tracking and disturbance-rejection performance for a wide range of process models. The proposed controllers were ranked mostly first. The closed-loop responses were stable and fast, with only slight overshoots.

Real-Time Experiment
This section presents a real-time experiment in the time domain. Note that the proposed method for the PI controller has already been tested on several plants. These experiments were presented in a previous study [1], where the tuning rule for a PI controller for integrating processes was tested. In this section, the proposed method will be tested on a charge-amplifier drift compensation system.

Control System for Charge-Amplifier Drift Compensation
The proposed tuning method was applied to a control system for charge-amplifier drift compensation [129]. As stated in Reference [129], "A charge amplifier is an electronic current integrator that is frequently employed for converting electrical charges or electrical currents into voltage signals. The charge amplifier is very sensitive to DC drift, since the DC component in the input signal leads to a steady accumulation of charge in the feedback capacitor until the output voltage saturates". As explained in Reference [129], the DC drift is compensated by an additional DC component with an opposite value for the charge-amplifier input. Since a charge-amplifier is an electronic integrator, the control process (charge amplifier) has an integrating character and the process input and output are DC voltages. Note that the charge amplifier is intended for low-level electric charge measurements. Therefore, it is a very precise low-noise instrument with stable characteristics.
The studied drift compensation control system was part of a custom-made measuring system for the low-frequency and high-temperature polarisation measurements of dielectric materials ( Figure 39). This system consists of a lock-in amplifier Stanford Research Systems SR830, a commercial charge-amplifier Kistler 5018A, and a temperature-regulated furnace (for measuring heated samples). More information about this dielectric-material characterisation concept can be found in References [130,131].

Charge amplifier (DUT) Sample
PID b or PID f controller r u Figure 39. Custom-made measurement system for low-frequency and high-temperature polarisation measurements of dielectric materials with a charge-amplifier drift compensation control system. The PID b , PID f2 , and PID fh controllers were realised in the software package LabVIEW, and data acquisition was realised with an NI USB-6001 card. The process output and input signals were limited to −10-10 V. The sampling time of the closed-loop control was 5 ms.
First, a step change (input voltage signal 0 to −1 V) was applied to the process input. Figure 40 shows the process response. Then, the following characteristic areas (with and without the added PID controller filter time constant, T F ) were calculated using Expressions (9) Based on the calculated characteristic Areas (62), by using Expression (34), the following second-order model of the process was calculated [126]: G M (s) = −8.5309e −0.0179s s(1 + 0.1081s + 0.0035s 2 ) .
The comparison of the model (G M ) and the actual process time responses are given in Figure 40.
Next, using Expressions (25) and (32)- (43), the reference filter parameters for the PID fh controller b F5 = 3.8442 × 10 −6 , b F4 = 2.1293 × 10 −4 , b F3 = 4.1389 × 10 −3 , b F2 = 4.5521 × 10 −2 , b F1 = 0.30178, T dm = 1.7874 × 10 −2 , K CL = 1.2, T CL = 0.1244, n = 2, (67) and for the PID f2 controller b F2 = 1.3029 × 10 −2 , b F1 = 0.16143, were calculated. Note that for the PID f2 and PID fh controllers, the proportional, integral, and derivative gains for PID b with b = 0 (65) were used. The process closed-loop responses for the proposed PID b , PID f2 , and PID fh controllers with a reference step change (r = 0.5) and an artificially added process input step-disturbance (d = 1) are shown in Figure 41. As can be seen, the proposed method results in the efficient control of drift compensation. The comparison illustrates the superior tracking performance of the PID f2 and PID fh controllers. As expected, the disturbance-rejection performance of the PID f2 and PID fh controllers is equal to that of the PID b controller with b = 0. This is confirmed by the criteria values presented in Figure 42. The IAT 2 E measure is, in general, not recommended for noisy measured signals (real plants). However, the process output noise was negligible in our case.

Conclusions
This paper presented an adaptation of the MOMI tuning method for the PID controller for the integrating processes. Additionally, we presented a reference filter that can substantially improve the tracking response while retaining the optimal disturbance-rejection response (i.e., the best overall response was achieved). The proposed method was tested on several different IP models. The closed-loop responses showed that the proposed method offers fast responses while keeping the system stable. The proposed PID f2 and PID fh controllers were also compared to other tuning methods for IPs. The closed-loop responses were evaluated by the IAE, IAT 2 E, and RVy × IAE criteria. The comparison with existing methods for IPs revealed that the proposed method can deliver enhanced performance.
The proposed method is also nonparametric, meaning that it does not require an explicit process model, making it an attractive approach for alleviating the drawbacks of the process-model mismatch. The controller parameters can be calculated from a process closed-or open-loop time response or from the general-order transfer function with a time-delay. Both approaches are equivalent. Furthermore, the user can emphasise disturbance-rejection or tracking with the reference weighting factor b or use the reference filter for the best overall response.
The proposed method was also tested on a control system for charge-amplifier drift compensation. The closed-loop responses in all experiments were fast and sufficiently damped. In future work, we will concentrate on achieving an overall optimal closed-loop response with a special focus on attenuating the controller's output high-frequency noise. Future work will also include more vigorous testing of the proposed method under process noise.
Author Contributions: Conceptualisation, T.K. and D.V.; software, T.K. and D.V.; validation, T.K. and D.V.; formal analysis, investigation, data curation, and writing, all authors; supervision, D.V. and M.H. All authors have read and agreed to the published version of the manuscript.
Funding: This work was carried out within research programs P2-0001 and PR-07603, financed by the Slovenian Research Agency, and by grants APVV SK-IL-RD-18-0008 and VEGA 1/0745/19.