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

Abstract

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.

1. 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 to a process plant will usually show that a large number of the PID controllers are poorly tuned”. Furthermore, as stated in Reference [9], “25% of all PID controller loops use default factory settings, implying that they have not been tuned at all”.

Various tuning methods have been proposed. The first tuning rules proposed by Ziegler and Nichols [10] were followed by many tuning techniques, such as the Cohen–Coon, Chien–Hrones–Reswick, and refined Ziegler–Nichols rules. In addition to these typical tuning rules, more advanced methods have also been proposed. The latter usually use complex algorithms for the process identification [9,11,12,13,14,15,16,17,18,19]. According to the literature, the number of tuning rules for stable overdamped processes is much larger than that for integrating processes. This suggests that PID tuning for integrating processes is not an extensively studied area.

Integrating processes (IPs) involve at least one pole at their origin and are frequently encountered in industry. The main characteristics of IPs are their open-loop instability that leads to an unbounded output for a bounded input, i.e., in case of disturbances, the output drifts to extreme levels. Accordingly, the efficient IPs closed-loop control is a demanding task [20]. The IPs can be classified according to the process order. For example, typical representatives of the processes that employ pure integrator plus time-delay (ITD) transfer function models are the oil–water–gas separators in the oil industry [20], the bottom-level control in distillation columns [21], the composition control loops of high-purity distillation columns [22], heat-integrated distillation columns [23], isothermal continuous copolymerisation reactors [24], storage tanks with outlet pumps [25], vertical take-off airplanes [26], industrial injection moulding machines [27], pulp and paper plants [20], and the high-pressure steam flow to steam-turbine generators in power plants [28]. First-order integrating systems with/without delay (IFOTD) can be encountered in liquid-storage tanks [29], jacketed continuous stirred tank reactors (CSTRs) carrying out exothermic reactions [30], and paper-drum-dryer cans [31]. The most frequently encountered example of an IP with inverse response with a time-delay (INPTD) is a boiler steam drum [32].

Additionally, the open-loop responses of various industrial processes, such as liquid-storage tanks, heating boilers, liquid level systems with pumps attached to the outflow, and batch chemical reactors, are slow, with a large dominant time constant. An IP model could be used when the required settling time is shorter than the dominant time constant of a process, even if the process dynamics are inherently of a first-order lag type [18,33]. As observed in Reference [21], controller design based on IPs could provide better closed-loop performance than designs based on a first-order plus time delay process.

Many researchers have addressed controller design methods for integrating and unstable processes. Some proposed solutions are based on using a single closed-loop controller, while others employ complex control structures. Numerous IP tuning methods have been proposed so far [9,34]. These techniques can be classified according to their IP model structure, i.e., tuning methods for ITD [20,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,61,62,63,64,65,66], for IFOTD [8,32,34,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107], for higher-order IPs with time-delay [108,109,110,111,112,113,114], and nonparametric tuning methods for IPs [1,115,116,117].

First, we present the tuning methods for ITD. Controller designs based on the IMC tuning rules for PID controllers were proposed in References [62,63,64]. In References [36,37], integral–proportional derivative (I-PD) control strategies were proposed. A two-degrees-of-freedom (2-DOF) PD-PI controller, that includes a second-order reference filter for improving tracking performance, was proposed in Reference [47]. A method for the auto-tuning of PID controllers based on closed-loop identification was proposed in Reference [38]. This method consists of the closed-loop identification of the model (after the reference or the disturbance step change), a performance assessment, and retuning of the PID controller. In Reference [39], a similar model-free PI/PID tuning method that uses a simple closed-loop reference step change (proportional control gain only) for obtaining the process model was proposed. In References [20,57], a PID controller and three-state structure for achieving a time-optimal tracking response was proposed. Since this method uses an identification technique (least-squares-based) for the calculation of the process model during steady-state changes, it is not necessary to know the process model in advance. A modified smith predictor for the control of stable, integrating, and unstable processes that consists of a 2-DOF PID controller in series with a second-order filter, defined by the dead-time and an adjustable parameter, was proposed in Reference [42]. In Reference [54], a novel approach for tuning discrete-time PID controllers with higher-order derivatives with a special focus on noise attenuation filters was proposed. These controllers offer a third degree of freedom in control design devoted to the attenuation of measurement noise. In Reference [55], two integrated tuning procedures for the joint design of a PID controller in series with an nth-order binomial low-pass filter were proposed. The proposed higher-order controller represents an intermediate step between integer and fractional-order PID controllers. Then, in Reference [56], a derivation of higher-order PID controllers by generalising the Skogestad Internal Model Control (SIMC) method [8] for analytical model-based controller tuning was proposed. In Reference [59], robust and simple analytical 2-DOF PI/PID controller tuning rules based on the multiple dominant pole method (MDPM) were proposed. The method enables calculation of controller parameters while maintaining non-oscillatory closed-loop responses without overshoots.

Note that the tuning methods for ITD can also be used for other types of IPs if the process can be modelled by the ITD model. This, for example, can be achieved by using a relay feedback identification method [37,38,41,42,43,45,46,65,101].

Next, the IFOTD tuning methods will be presented. In Reference [105], an analytical PID tuning method, where the tuning of parameters is based on maximum sensitivity (M_S), was proposed. To reduce overshoots, a second-order reference filter is employed. Similar tuning rules were proposed by the same authors in Reference [106], where overshoots are reduced by cancelling some controller zeros and using fourth-order reference filtering. In Reference [97], analytical tuning rules for a 2-DOF PID controller, where overall optimal closed-loop response is obtained via two-step tuning of the closed-loop transfer function, were proposed. Controller designs based on IMC PID tuning rules were proposed in References [8,34,67,68,69,70,107]. In Reference [71], a modified IMC design based on a 2-DOF control structure that allows the separate optimisation of tracking and disturbance-rejection performance was proposed. In Reference [72], an IMC PD controller with a fractional-order filter was proposed. In Reference [73], particle swarm optimisation (PSO) and artificial bee colony (ABC) algorithms were used for PID controller tuning. Tuning rules, obtained by optimisation, were also proposed in References [74] for PI/PID, [75] for PI/PID, and [76] for I-PD. The method proposed in Reference [75] applied an additional reference filter to reduce overshoots in the tracking response. In Reference [32], analytical PI/PID controller tuning rules for IPs with non-minimum phase zeros were proposed. Calculation of the controller parameters, based on the model and the minimum Integral of the Absolute Error (IAE) criterion, was accomplished via the golden section searching algorithm. To create a 2-DOF control system, a PID controller with a unique reference filter was proposed in Reference [83]. A PID controller with set-point weighting in series with a lead/lag compensator was proposed in Reference [84]. PI/PID controller tuning rules and a method for the identification of IPs with a time-delay and inverse response from an open-loop step response were proposed in Reference [100]. In Reference [86], a methodology for the closed-loop performance assessment of a PID controller is proposed. First, the process is approximated as IFOTD from the closed-loop step response, and then the performance is evaluated by comparing it with the performance provided by the SIMC [8] tuning rule. In case the performance is not satisfactory, the PID controller is retuned using the same tuning rule. Analytical PI/PID controller tuning rules for ITD and for IFOTD based on the MDPM were proposed in Reference [102]. If the closed-loop tracking responses have higher overshoots, the method employs a 2-DOF PID controller.

Note that most tunings methods for IFOTD can also be applied to other types of IPs if the process can be modelled by the IFOTD model [8,68,69,74,75,77,78,79,80,84,87,88,91,92,93,94,95,96,98,100,106,107].

Next, we review the tuning methods for general (higher-order) IPs with time-delay. A symmetrical–optimum method for tuning PID controllers for general IPs was proposed in Reference [108]. Although this method results in higher overshoots for tracking, it is effective for disturbance 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.

2. MOMI Tuning Method for Integrating Processes

2.1. 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\left(s\right)=\frac{K_{PR}}{s}\frac{1+b_{1}s+b_2{s}^2+\cdots +b_m{s}^m}{1+a_{1}s+a_2{s}^2+\cdots +a_n{s}^n}{e}^{-s{T}_{delay}},$$

(1)

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

$$G_{CL}\left(s\right)=\frac{Y\left(s\right)}{R\left(s\right)}=\frac{G_P\left(s\right)G_{CR}\left(s\right)}{1+G_P\left(s\right)G_{CY}\left(s\right)}.$$

(2)

The controller is configured so that

$$G_{CL}\left(0\right)=1$$

(3)

$$\underset{\omega \to 0}{\lim }\left[\frac{d^{2k}{\left|G_{CL}\left(j\omega \right)\right|}^2}{d{\omega }^{2k}}\right]=0;k=1,2,\cdots ,k_{max}$$

(4)

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:

$$G_{CL}\left(s\right)=\frac{f_0+f_{1}s+f_2{s}^2+\cdots }{e_0+e_{1}s+e_2{s}^2+\cdots }.$$

(5)

Then, the following conditions must be satisfied to fulfil Expression (4) [121]:

$${{\displaystyle \sum }}_{i=0}^{2n}{\left(-1\right)}^{i+n}\left(f_i{f}_{2n-i}{e}_0{}^2-e_i{e}_{2n-i}{f}_0{}^2\right)=0;n=1,2,\dots$$

(6)

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]:

$$G_P\left(s\right)=\frac{G_{P1}}{s}=\frac{A_0-A_{1}s+A_2{s}^2-A_3{s}^3+\cdots }{s},$$

(7)

where A_i represents the “characteristic areas” of the process [121], which can be expressed by Process Parameters (1) [119]:

$$\begin{gathered} A_0=K_{PR} \\ {A}_1=K_{PR}\left(a_1-b_1+T_{delay}\right) \\ {A}_2=K_{PR}\left(b_2-a_2-T_{delay}{b}_1+\frac{T_{delay}{}^2}{2!}\right)+A_1{a}_1 \\ {A}_3=K_{PR}\left(a_3-b_3+T_{delay}{b}_2-\frac{T_{delay}{}^2{b}_1}{2!}+\frac{T_{delay}{}^3}{3!}\right)+A_2{a}_1-A_1{a}_2. \\ {A}_4=K_{PR}\left(-a_4+b_4-T_{delay}{b}_3+\frac{T_{delay}{}^2{b}_2}{2!}-\frac{T_{delay}{}^3{b}_1}{3!}+\frac{T_{delay}{}^4}{4!}\right)+A_3{a}_1-A_2{a}_2+A_1{a}_3 \\ ⋮ \\ {A}_k=K_{PR}\left({\left(-1\right)}^{k+1}\left(a_k-b_k\right)+{\displaystyle \sum }_{i=1}^k{\left(-1\right)}^{k+i}\frac{T_{delay}{}^i{b}_{k-i}}{i!}\right)+{\displaystyle \sum }_{i=1}^{k-1}{\left(-1\right)}^{k+i-1}{A}_i{a}_{k-i}. \end{gathered}$$

(8)

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]:

$$\begin{array}{c}{u}_0=\frac{u\left(t\right)-u\left(0\right)}{\Delta U}\hfill \\ I_{U1}\left(t\right)={{\displaystyle \int }}_0^t{u}_0\left(\tau \right)d\tau I_{Y1}\left(t\right)=\frac{y\left(t\right)-y\left(0\right)-\dot{y}\left(0\right)\cdot t}{\Delta U},\hfill \\ I_{U2}\left(t\right)={{\displaystyle \int }}_0^t{I}_{U1}\left(\tau \right)d\tau I_{Y2}\left(t\right)={{\displaystyle \int }}_0^t{I}_{Y1}\left(\tau \right)d\tau \hfill \\ ⋮⋮\hfill \end{array}$$

(9)

where u₀ represents the normalised process input signal, and

$$\Delta U=u(\infty )-u\left(0\right).$$

(10)

The characteristic areas are expressed as follows:

$$\begin{array}{c}{A}_0=y_0(\infty );y_1=A_0{I}_{U1}\left(t\right)-I_{Y1}\left(t\right)\hfill \\ A_1=y_1(\infty );y_2=A_1{I}_{U1}\left(t\right)-A_0{I}_{U2}\left(t\right)+I_{Y2}\left(t\right)\hfill \\ A_2=y_2(\infty );y_3=A_2{I}_{U1}\left(t\right)-A_1{I}_{U2}\left(t\right)+A_0{I}_{U3}\left(t\right)-I_{Y3}\left(t\right).\hfill \\ ⋮\hfill \end{array}$$

(11)

Note that $\ddot{y}\left(0\right)=\stackrel{⃛}{y}\left(0\right)=\cdots =0$ and that A₀ 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₁. Signals u₀ and y₀ represent the process input and process output responses, respectively [121].

2.2. Extension to Integrating Processes

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:

$$\begin{gathered} G_{CR}\left(s\right)=\frac{b{K}_D{s}^2+b{K}_{P}s+K_I}{s\left(1+s{T}_F\right)} \\ {G}_{CY}\left(s\right)=\frac{K_D{s}^2+K_{P}s+K_I}{s\left(1+s{T}_F\right)}. \end{gathered}$$

(12)

To simplify the derivations, c = b is chosen. Thereafter, this controller will be denoted PID_b. While the controller filter

$$G_F=\frac{1}{1+s{T}_F}$$

(13)

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}=\frac{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:

$$T_F=\frac{K_D}{K_{HF}}$$

(14)

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:

$$\begin{gathered} e_0=A_0{K}_I \\ {e}_1=-A_1{K}_I+A_0{K}_P \\ {e}_2=A_2{K}_I-A_1{K}_P+A_0{K}_D+1 \\ {e}_3=-A_3{K}_I+A_2{K}_P-A_1{K}_D \\ {e}_4=A_4{K}_I-A_3{K}_P+A_2{K}_D \\ ⋮ \end{gathered}$$

(15)

$$\begin{gathered} f_0=A_0{K}_I \\ {f}_1=-A_1{K}_I+A_{0}b{K}_P \\ {f}_2=A_2{K}_I-A_{1}b{K}_P+A_{0}b{K}_D \\ {f}_3=-A_3{K}_I+A_{2}b{K}_P-A_{1}b{K}_D \\ {f}_4=A_4{K}_I-A_{3}b{K}_P+A_{2}b{K}_D, \\ ⋮ \end{gathered}$$

(16)

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₀K_I. Consequently, the following closed-loop transfer function Parameters (5) are obtained:

$$\begin{gathered} e_0=1 \\ {e}_1=G_{10}+T_I \\ {e}_2=G_{20}+G_{10}{T}_I+K_{DI}+G_{00}{K}_{I1} \\ {e}_3=G_{30}+G_{20}{T}_I+G_{10}{K}_{DI} \\ {e}_4=G_{40}+G_{30}{T}_I+G_{20}{K}_{DI} \\ ⋮ \end{gathered}$$

(17)

$$\begin{gathered} f_0=1 \\ {f}_1=G_{10}+b{T}_I \\ {f}_2=G_{20}+G_{10}b{T}_I+b{K}_{DI} \\ {f}_3=G_{30}+G_{20}b{T}_I+G_{10}b{K}_{DI} \\ {f}_4=G_{40}+G_{30}b{T}_I+G_{20}b{K}_{DI}, \\ ⋮ \end{gathered}$$

(18)

where

$$T_I=\frac{K_P}{K_I},K_{i1}=\frac{1}{K_I},K_{Di}=\frac{K_D}{K_I},G_{00}=\frac{1}{A_0},G_{10}=-\frac{A_1}{A_0},G_{20}=\frac{A_2}{A_0},G_{30}=-\frac{A_3}{A_0},G_{40}=\frac{A_4}{A_0}.$$

(19)

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:

$$K_{Di}=\frac{0.5\left(1-b^2\right)T_I{}^2-\frac{G_{00}}{K_I}}{\left(1-b\right)},$$

(20)

and the second condition in Expression (6) provides the following result:

$$K_I=-\frac{G_{00}{\left(b-1\right)}^2(-2G_{10}{}^2-2G_{10}{T}_I+2G_{20}+b{T}_I{}^2\left(1+b\right))}{0.25{\left(b^2-1\right)}^3{T}_I{}^4}.$$

(21)

The third condition of Expression (6) gives the following fourth-order polynomial:

$$\alpha T_I{}^4+\beta T_I{}^3+\left({\gamma }_1+{\gamma }_2\right)T_I{}^2+\delta T_I+\epsilon =0,$$

(22)

where

$$\begin{array}{c}\mathsf{\alpha }=\left(1-b\right){\left(1+b\right)}^2\left(G_{10}{}^2\left(1-b\right)+2b{G}_{20}\right),\hfill \\ \beta =8\left(\left(1-b^2\right)G_{10}{}^3+\left(2b^2+b-1\right)G_{10}{G}_{20}-b\left(1+b\right)G_{30}\right),\hfill \\ {\gamma }_1=8\left(-\left(b^2-3\right)G_{10}{}^4+\left(3b^2+b-6\right)G_{10}{}^2{G}_{20}-2\left(b^2+b-1\right)G_{10}{G}_{30}\right),\hfill \\ {\gamma }_2=8\left(\left(1-b^2\right)G_{20}{}^2+b\left(1+b\right)G_{40}\right),\hfill \\ \delta =16\left(2G_{10}{}^5-6G_{10}{}^3{G}_{20}+3\left(G_{10}{}^2{G}_{30}+G_{10}{G}_{20}{}^2\right)-G_{10}{G}_{40}-G_{20}{G}_{30}\right),\hfill \\ \epsilon =16\left(G_{10}{}^6-4G_{10}{}^4{G}_{20}+2G_{10}{}^3{G}_{30}+4G_{10}{}^2{G}_{20}{}^2-G_{10}{}^2{G}_{40}-2G_{10}{G}_{20}{G}_{30}-G_{20}{}^3+G_{20}{G}_{40}\right).\hfill \end{array}$$

(23)

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:

$$\begin{gathered} K_P=K_I{T}_I, \\ {K}_D=\frac{0.5\left(1-b^2\right)T_I{}^2{K}_I-G_{00}}{\left(1-b\right)}. \end{gathered}$$

(24)

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:

$$\begin{gathered} K_P=\frac{K_D{A}_0+0.5}{A_1}, \\ {K}_I=0, \\ {K}_D=\frac{A_3}{2\left(A_1{A}_2-A_0{A}_3\right)}. \end{gathered}$$

(25)

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.

Remark 2.

The PID_b controller parameters calculated by Expressions (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.

As explained in Section 2.1, the characteristic areas A₀ to A₄ 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:

• Choose the appropriate filter time constant, T_F, of the PID_b controller.
• Determine the characteristic areas A₀–A₄. 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.
• Select the reference weighting factor b. As seen in Remark 1, the recommended values are 0 ≤ b < 0.9.
• 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).

2.3. Illustrative Example

The following second-order delayed IP is chosen to illustrate the proper design:

$$G_P\left(s\right)=\frac{e^{-0.2s}}{s{\left(1+s\right)}^2}.$$

(26)

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:

$$A_0=1, A_1=2.21, A_2=3.4421, A_3=4.6758, A_4=5.9095.$$

(27)

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. The PID_b controller parameters.

b	KP	KI	KD
0	0.93011	0.23327	0.85427
0.2	0.93383	0.2486	0.85466
0.5	0.85191	0.19158	0.84114
0.7	0.75396	0.11619	0.8251
0.9	0.64303	0.03635	0.80683

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.

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.

3. 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).

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).

3.1. The Second-Order Filter

The closed-loop transfer function G_CL from the reference to the process output is the following:

$$G_{CL}\left(s\right)=\frac{G_F\left(K_I+K_{P}s+K_D{s}^2\right)G_P}{s+\left(K_I+K_{P}s+K_D{s}^2\right)G_P}.$$

(28)

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):

$$G_{CL1}\left(s\right)=\frac{\left(K_{P1}+K_{D1}s\right)G_P}{s+\left(K_{P1}+K_{D1}s\right)G_P}.$$

(29)

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:

$$G_F\left(s\right)=\frac{1+b_{F1}s+b_{F2}{s}^2}{\left(1+\frac{K_P}{K_I}s+\frac{K_D}{K_I}{s}^2\right)}.$$

(30)

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₀ and A₁) 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:

$$\begin{gathered} b_{F1}=\frac{K_P}{K_I}-\frac{1}{A_0{K}_{P1}}, \\ {b}_{F2}=\frac{\left(A_0{}^2{K}_D{K}_{P1}{}^2+A_0{K}_{P1}{}^2-A_0{K}_P{K}_{P1}+A_0{K}_I{K}_{D1}-A_1{K}_I{K}_{P1}+K_I\right)}{A_0{}^2{K}_I{K}_{P1}{}^2}. \end{gathered}$$

(31)

3.2. 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

$$G_M\left(s\right)=\frac{K_M\left(1+b_{1m}s\right)e^{-T_{dm}s}}{s\left(1+a_{1m}s+a_{2m}{s}^2\right)},$$

(32)

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₀ to A₃ (8):

$$\begin{gathered} K_M=A_0, \\ {a}_{1m}=\frac{2A_1\left(A_2-\frac{A_0{T}_{dm}^2}{2}\right)-2A_0\left(A_3-\frac{A_0{T}_{dm}^3}{6}\right)+\left(A_1-A_0{T}_{dm}\right)\left(A_0{T}_{dm}^2-2A_1{T}_{dm}\right)}{A_0^2{T}_{dm}^2-2A_0{A}_1{T}_{dm}-2A_0{A}_2+2A_1^2}, \\ {a}_{2m}=\frac{\left(A_2-\frac{A_0{T}_{dm}^2}{2}\right)(2A_2-A_0{T}_{dm}^2)-2\left(A_1-A_0{T}_{dm}\right)\left(A_3-\frac{A_0{T}_{dm}^3}{6}\right)+\left(A_1-A_0{T}_{dm}\right)\left(A_1{T}_{dm}^2-2A_2{T}_{dm}\right)}{A_0^2{T}_{dm}^2-2A_0{A}_1{T}_{dm}-2A_0{A}_2+2A_1^2}, \\ {b}_{1m}=\frac{2A_1\left(A_2-\frac{A_0{T}_{dm}^2}{2}\right)-2A_0\left(A_3-\frac{A_0{T}_{dm}^3}{6}\right)-2\left(\frac{A_1}{A0}-T_{dm}\right)\left(A_1^2-A_0{A}_2\right)}{A_0^2{T}_{dm}^2-2A_0{A}_1{T}_{dm}-2A_0{A}_2+2A_1^2}. \end{gathered}$$

(33)

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):

$$\begin{gathered} K_M=A_0, \\ \frac{T_{dm}^3}{6}-\frac{T_{dm}^2{A}_1}{2A_0}+T_{dm}\left(\frac{A_1^2}{A_0^2}-\frac{A_2}{A_0}\right)-\frac{A_1^3}{A_0^3}+\frac{2A_1{A}_2}{A_0^2}+\frac{A_3}{A_0}=0, \\ {a}_{1m}=\frac{A_1}{A_0}-T_{dm}, \\ {a}_{2m}=-\frac{A_2+A_1{a}_{1m}}{A_0}+\frac{T_{dm}^2}{2}, \\ {b}_{1m}=0, \end{gathered}$$

(34)

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:

$$G_{CL}\left(s\right)=\frac{G_F{K}_M\left(K_D{s}^2+K_{P}s+K_I\right)\left(1+b_{1m}s\right)e^{-T_{dm}s}}{\left(1+a_{1m}s+a_{2m}{s}^2\right)s^2+K_M\left(K_D{s}^2+K_{P}s+K_I\right)\left(1+b_{1m}s\right)e^{-T_{dm}s}}.$$

(35)

We can now define the desired closed-loop transfer function as follows:

$$G_{CLD}\left(s\right)=\frac{e^{-T_{dm}s}}{{\left(1+\frac{T_{CL}}{n}s\right)}^n},$$

(36)

where

$$\begin{array}{cc}n=2& T_{dm}>0\\ n=3& T_{dm}=0\end{array}.$$

(37)

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:

$$T_{CLO}={{\displaystyle \int }}_{t=0}^{t=\infty }e\left(t\right)dt=\frac{1}{A_0{K}_{P1}},$$

(38)

where K_P1 is the proportional gain of the PD controller. The desired closed-loop time constant T_CL is then defined as follows:

$$T_{CL}=\frac{T_{CL0}}{K_{CL}}=\frac{1}{A_0{K}_{P1}{K}_{CL}},$$

(39)

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:

$$e^{-T_{dm}s}\approx \frac{1-\frac{T_{dm}}{2}s+\frac{T_{dm}^2}{12}{s}^2}{1+\frac{T_{dm}}{2}s+\frac{T_{dm}^2}{12}{s}^2}.$$

(40)

By Equations (35) and (36) and taking into account (40) in the denominator of (35), the following reference filter transfer function is obtained:

$$G_F=\frac{1+b_{F1}s+b_{F2}{s}^2+b_{F3}{s}^3+b_{F4}{s}^4+b_{F5}{s}^5}{\left(1+\frac{K_P}{K_I}s+\frac{K_D}{K_I}{s}^2\right)\left(1+\frac{T_{dm}}{2}s+\frac{T_{dm}^2}{12}{s}^2\right){\left(1+\frac{T_{CL}}{n}s\right)}^n},$$

(41)

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:

$$\begin{gathered} de{n}_{CL}=\frac{1}{A_0{K}_I}\left(1+a_{1m}s+a_{2m}{s}^2\right)\left(1+\frac{T_{dm}}{2}s+\frac{T_{dm}^2}{12}{s}^2\right)s^2+\left(1+\frac{K_{P }}{K_I}s+\frac{K_D}{K_I}{s}^2\right)\times \\ \left(1+b_{1m}s\right)\left(1-\frac{T_{dm}}{2}s+\frac{T_{dm}^2}{12}{s}^2\right) \end{gathered}$$

(42)

by equating the terms with the same power of the Laplace variable s. The resulting parameters are the following:

$$\begin{gathered} b_{F1}=b_{1m}-\frac{T_{dm}}{2}+\frac{K_P}{K_I}, \\ {b}_{F2}=\frac{T_{dm}\left(T_{dm}-6b_{1m}\right)}{12}+\frac{ \left(A_0{K}_D-\frac{A_0{K}_P{T}_{dm}}{2}+A_0{K}_P{b}_{1m}+1\right)}{A_0{K}_I}, \\ {b}_{F3}=\frac{K_D{b}_{1m}-\frac{T_{dm}\left(K_D+K_P{b}_{1m}\right)}{2} +\frac{T_{dm}^2\left(K_P+K_I{b}_{1m}\right)}{12}}{K_I}+\frac{\frac{T_{dm}}{2}+T_F+a_{1m}}{A_0{K}_I}, \\ {b}_{F4}=\frac{\frac{T_{dm}^2\left(K_D+K_P{b}_{1m}\right)}{12}-\frac{\left(K_D{T}_{dm}{b}_{1m}\right)}{2}}{K_I}+\frac{a_{2m}+T_F{a}_{1m}+\frac{T_{dm}\left(T_F+a_{1m}\right)}{2}+\frac{T_{dm}^2}{12}}{A_0{K}_I}, \\ {b}_{F5}=\frac{\frac{T_{dm}^2\left(T_F+a_{1m}\right)}{12}+\frac{T_{dm}\left(a_{2m}+T_F{a}_{1m}\right)}{2}+T_F{a}_{2m}}{A_0{K}_I}+\frac{K_D{T}_{dm}^2{b}_{1m}}{12K_I}. \end{gathered}$$

(43)

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:

$$\left(1+\frac{K_P}{K_I}s+\frac{K_D}{K_I}{s}^2\right).$$

(44)

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.

4. 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].

4.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:

$$\begin{array}{c}{G}_{P1}\left(s\right)=\frac{e^{-0.2s}}{s\left(1+s\right)\left(1+0.2s\right)}\hfill \\ G_{P2}\left(s\right)=\frac{e^{-s}}{s} \hfill \\ G_{P3}\left(s\right)=\frac{1}{s{\left(1+s\right)}^4}\hfill \\ G_{P4}\left(s\right)=\frac{1-2s}{s{\left(1+s\right)}^2} .\hfill \end{array}$$

(45)

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₀ to A₄ were calculated. The values are presented in

Table 2. The calculated characteristic areas for the process models with the additional filter time constant T_F = 0.01 s.

	A0	A1	A2	A3	A4
GP1	1	1.41	1.5141	1.5365	1.541
GP2	1	1.01	0.5101	0.17177	0.04338
GP3	1	4.01	10.0401	20.1004	35.201
GP4	1	4.01	7.0401	10.0704	13.1007

.

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 PID_b controller parameters.

	b = 0			b = 0.5			b = 0.7
	KP	KI	KD	KP	KI	KD	KP	KI	KD
GP1	2.6269	1.3991	1.4661	2.2778	1.1385	1.418	1.8731	0.63473	1.365
GP2	1.128	0.47494	0.33942	1.0017	0.3258	0.30958	0.89095	0.18652	0.28414
GP3	0.37914	0.04554	0.57814	0.3454	0.03499	0.55729	0.30728	0.02075	0.53411
GP4	0.33133	0.03935	0.39492	0.2785	0.02479	0.34757	0.24024	0.01344	0.31634

.

The calculated controller parameters were applied to a PID_b controller. Figure 7, Figure 8, Figure 9 and Figure 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).

4.2. 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:

$$\begin{array}{c}{G}_{P1}\left(s\right)=\frac{e^{-0.5s}}{s\left(1+4s\right)}\hfill \\ G_{P2}\left(s\right)=\frac{e^{-0.5s}}{s{\left(1+2s\right)}^2}\hfill \\ G_{P3}\left(s\right)=\frac{e^{-0.5s}}{s{\left(1+s\right)}^4}\hfill \\ G_{P4}\left(s\right)=\frac{e^{-4s}}{s\left(1+0.5s\right)} .\hfill \end{array}$$

(46)

First, the parameters for the PI controller were acquired. To that end, the characteristic areas A₀–A₂, 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]:

$$\begin{gathered} K_I=0.5A_0\left(1-b^2\right)K_P{}^2, \\ {K}_P=\frac{-A_1+\sqrt{A_1{}^2+\xi }}{\xi }, \\ \xi =\left(1-b^2\right)\left(A_0{A}_2-A_1{}^2\right), \end{gathered}$$

(47)

using the calculated characteristic areas from

Table 4. The calculated characteristic areas for the process models.

	A0	A1	A2
GP1	1	4.5	18.125
GP2	1	4.5	14.125
GP3	1	4.5	12.125
GP4	1	4.5	10.25

. The calculated PI controller parameters for all four IP models are presented in

Table 5. The calculated PI controller parameters.

	b = 0		b = 0.5
	KI	KP	KI	KP
GP1	0.0065	0.1142	0.0048	0.1134
GP2	0.0073	0.1211	0.0052	0.1182
GP3	0.0078	0.1253	0.0055	0.121
GP4	0.0084	0.1298	0.0058	0.1239

.

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₀ to A₄ presented in

Table 6. The calculated characteristic areas for the process models with the additional filter time constant T_F = 0.01 s.

	A0	A1	A2	A3	A4
GP1	1	4.51	18.170	72.703	290.813
GP2	1	4.51	14.170	38.663	97.973
GP3	1	4.51	12.170	25.643	46.592
GP4	1	4.51	10.295	15.895	18.721

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 PID_b controller parameters.

	b = 0			b = 0.5
	KP	KI	KD	KP	KI	KD
GP1	4.4939	1.7606	4.7353	3.4387	1.3668	4.4883
GP2	0.44499	0.054126	0.82917	0.40688	0.044121	0.81421
GP3	0.31757	0.033077	0.52446	0.28765	0.024817	0.50058
GP4	0.25619	0.024309	0.34998	0.22777	0.01677	0.32033

.

The calculated controller parameters presented in Table 5 and Table 7 were applied to the PI and PID_b controllers. Figure 12, Figure 13, Figure 14 and Figure 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.

4.3. 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₀ to A₄ 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 intermediate calculated parameters for the reference filters calculations.

Process	KP1	KD1	KM	b1m	Tdm	a1m	a2m	KCL
GP1	1.2651	1.2838	1	0	0.2	1.2	0.2	3
GP2	0.74265	0.25007	1	0	1	0.0101	0	3
GP3	0.249005	0.498512	1	0	0.935822	3.06418	2.69459	3
GP4	0.193831	0.277263	1	−1.99879	0.011	1.99021	0.982906	3

. The calculated reference filter parameters for the PID_f2 and PID_fh controllers are given in

Table 9. The PID_f2 controller reference filter parameters.

Process	bF2	bF1
GP1	0.59091	1.0871
GP2	0.52881	1.0284
GP3	9.2833	4.3089
GP4	5.3174	3.2611

and

Table 10. The PID_fh controller reference filter parameters.

Process	bF5	bF4	bF3	bF2	bF1	Tdm	TCL	n
GP1	0.019465	0.24389	0.83779	1.5782	1.7776	0.2	0.26348	2
GP2	0.0036331	0.25639	0.93567	1.716	1.8749	1	0.44884	2
GP3	33.5171	93.952	72.4423	30.8295	7.85696	0.935822	1.33866	2
GP4	0.390273	25.8753	30.9495	18.5847	6.41596	0.011	1.71971	2

, respectively.

Figure 16, Figure 17, Figure 18 and Figure 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).

5. 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):

$$G_{CL}=\frac{G_P{G}_{CR}}{1+G_P{G}_{CY}}=\frac{nu{m}_{CL}}{de{n}_{CL}},$$

(48)

where

$$de{n}_{CL}=s\left(1+a_{1}s+a_2{s}^2+\dots \right)\left(1+s{T}_F\right)e^{s{T}_{delay}}+\left(1+b_{1}s+b_2{s}^2+\dots \right)\left(K_D{s}^2+K_{P}s+K_I\right).$$

(49)

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:

$$G_P\left(s\right)=\frac{K_{PR}{e}^{-T_{delay}s}}{s\left(1+a_{1}s+\alpha a_1^2{s}^2\right)}.$$

(50)

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₁ 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_PG_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 Figure 21, Figure 22 and Figure 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:

$$G_{P1}\left(s\right)=\frac{e^{-s}}{s{\left(1+s\right)}^2}.$$

(51)

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:

$$A_0=1, A_1=3.01, A_2=5.5301, A_3=8.22197, A_4=10.9572.$$

(52)

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):

$$\begin{gathered} 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. \end{gathered}$$

(53)

Additionally, the parameters of the PID_fh controller reference filter were calculated using Expression (43) (for chosen K_CL = 3):

$$\begin{gathered} b_{F5}=8.2389,b_{F4}=25.8477,b_{F3}=27.1014,b_{F2}=15.97,b_{F1}=5.6589, \\ {T}_{dm}=1,T_{\mathit{CL}}=1.0155,n=2. \end{gathered}$$

(54)

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. Figure 24 and Figure 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).

6. 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\times 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:

$$\begin{gathered} \mathrm{IAE}={{\displaystyle \int }}_0^{\infty }\left|\mathrm{e}\right|\left(t\right)dt, \\ \mathrm{I}\mathrm{A}{\mathrm{T}}^2\mathrm{E}={{\displaystyle \int }}_0^{\infty }{t}^2\cdot \left|\mathrm{e}\right|\left(t\right)dt, \\ \mathrm{RVy}\times \mathrm{IAE}, \end{gathered}$$

(55)

where RVy is the relative variance of the process output signal:

$$\begin{gathered} \mathrm{RVy}\left(\mathrm{tracking}\right)= \frac{{{\displaystyle \int }}_0^{\infty }\left|\frac{dy\left(t\right)}{dt}\right|dt}{\left|y(\infty )-y\left(0\right)\right|}, \\ \mathrm{RVy}\left(\mathrm{disturbance}\right)= \frac{{{\displaystyle \int }}_0^{\infty }\left|\frac{dy\left(t\right)}{dt}\right|dt}{\left|y_1-y\left(0\right)\right|+\left|y_2-y_1\right|+\left|y(\infty )-y_2)\right|}, \end{gathered}$$

(56)

and y₁ and y₂ 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 ($\mathrm{RVy}\times \mathrm{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 $\mathrm{RVy}\times \mathrm{IAE}$, the IAT²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

$$G_{P1}\left(s\right)=\frac{e^{-s}}{s}.$$

(57)

The tracking and the disturbance-rejection closed-loop responses between the aforementioned tuning rules are compared in Figure 26 and Figure 27. The criteria values of the tuning rules are presented in Figure 28.

According to the IAE criterion, the proposed controllers ranked first, and according to the IAT²E criterion, they ranked first and third in tracking performance. The proposed PID_fh controller also ranked first for the $\mathrm{RVy}\times \mathrm{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:

$$G_{P2}\left(s\right)=\frac{e^{-0.35s}}{s\left(1+1.35s\right)}.$$

(58)

Figure 29 and Figure 30 compare the closed-loop tracking and disturbance-rejection responses for all mentioned tuning rules. The criteria values of the aforementioned tuning rules are presented in Figure 31.

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:

$$G_{P3}\left(s\right)=\frac{e^{-s}}{s{\left(1+1.5s\right)}^2} .$$

(59)

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:

$$G_{P4}\left(s\right)=\frac{0.5\left(1-0.5s\right)e^{-0.7s}}{s\left(1+0.4s\right)\left(1+0.1s\right)\left(1+0.5s\right)}.$$

(60)

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 Figure 34 and Figure 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²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 $\mathrm{RVy}\times \mathrm{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:

$$G_{P4}\left(s\right)=\frac{1}{s{\left(1+0.25s\right)}^8}.$$

(61)

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.

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.

7. 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].

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)–(11):

$$A_0=-8.5309, A_1=-1.0748, A_2=-0.0880, A_3=-5.7974\times {10}^{-3},A_4=-3.1956\times {10}^{-4}.$$

(62)

$$\begin{gathered} T_F=0.02, \\ {A}_{0\mathrm{F}}=-8.5309, A_{1\mathrm{F}}=-1.2454, A_{2\mathrm{F}}=-0.11291, A_{3\mathrm{F}}=-8.0556\times {10}^{-3},A_{4\mathrm{F}}=-4.8067\times {10}^{-4}. \end{gathered}$$

(63)

Based on the calculated characteristic Areas (62), by using Expression (34), the following second-order model of the process was calculated [126]:

$$G_M\left(s\right)=\frac{-8.5309e^{-0.0179s}}{s\left(1+0.1081s+0.0035s^2\right)}.$$

(64)

The comparison of the model (G_M) and the actual process time responses are given in Figure 40.

The next PID_b controller parameters for b = 0 and b = 0.5 were calculated via Equations (21) and (24) using the calculated characteristic Areas (63):

$$b=0, K_P=-1.1655, K_I=-3.7512, K_D=-0.063854,$$

(65)

$$b=0.5, K_P-1.0688, K_I=-2.8916, K_D=-0.061832.$$

(66)

Next, using Expressions (25) and (32)–(43), the reference filter parameters for the PID_fh controller

$$\begin{gathered} b_{F5}=3.8442\times {10}^{-6}, b_{F4}=2.1293\times {10}^{-4}, b_{F3}=4.1389\times {10}^{-3}, b_{F2}=4.5521\times {10}^{-2}, \\ {b}_{F1}=0.30178,T_{dm}=1.7874\times {10}^{-2},K_{\mathit{CL}}=1.2,T_{\mathit{CL}}=0.1244,n=2, \end{gathered}$$

(67)

and for the PID_f2 controller

$$b_{F2}=1.3029\times {10}^{-2}, b_{F1}=0.16143,$$

(68)

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²E measure is, in general, not recommended for noisy measured signals (real plants). However, the process output noise was negligible in our case.

8. 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²E, and $\mathrm{RVy}\times \mathrm{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.