Influence of the Selection of Reaction Curve’s Representative Points on the Accuracy of the Identified Fractional-Order Model

Abstract

In this paper, a general procedure for identifying a fractional first-order plus dead-time (FFOPDT) model is presented. This procedure is based on fitting three arbitrary points on the process reaction curve, where process information is obtained from a simple open-loop test. A simplification of the general identification procedure is also considered, where only points symmetrically located on the reaction curve are selected. The proposed symmetrical procedure has been applied to the following sets of representative points: (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%). Analytical expressions of the corresponding FFOPDT model parameters for these sets of symmetrical points have been obtained. In order to show the effectiveness and applicability of this procedure for the identification of fractional-order models and to get insight into the influence of selection of the set of symmetrical points on the accuracy of the identified model, some numerical examples are proposed. This identification procedure gives good results in comparison with other integer- and fractional-order identification methods. Finally, some conclusions and final remarks are offered in this context.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Although it might be possible to obtain a model of the process analytically—a physical model—, it is more frequent that the dynamic information of the process is obtained through experimental tests—empirical model—from the process reaction curve, the critical gain and the period of oscillations of the control system at the limit of stability, or by feedback with a relay [1].

From the controller point of view, the controlled process includes the process itself, the final controlling element, and the measuring instrument. In control studies and for the design and tuning of the controller in a feedback control loop, it is necessary to have information about the dynamic behaviour of the controlled process, usually in the form of a reduced-order mathematical model [2]. The controlled process model provides the dynamic information between the controller and the measuring instrument output signals. This model should be simple, but, at the same time, it should provide reliable information about the controlled process behaviour at the operating point. This information usually includes gain, time constant (s), and apparent dead-time of the controlled process.

In spite of all the considerable advances in process control over the past several decades, it is common knowledge that proportional integral derivative (PID) controllers are widely and successfully used in industrial control applications, having become an industrial standard for process control [3]. It is acknowledged by the academic and industrial community that the most widely used controlled process models for this control algorithm are first-order, dual-pole, and second-order plus dead-time (FOPDT, DPPDT, SOPDT) ones; however, in some cases, these models cannot represent the process dynamics with the required accuracy. This was already suggested in, e.g., [4], where it is indicated that it is necessary to have better information about the process to obtain optimal PID tuning rules for lag-dominated processes.

Practical experience shows that controller tuning can be generally accomplished with very little information about the plant from the point of view of standard design techniques. Experimental process model identification using an open-loop step test is widely covered in the technical literature; see, e.g., the identification procedures described in [5–8].

Some papers also reported identification algorithms based on fitting several representative points in the process transient response to a step change [9–11]. In the technical literature, there are several two- or three-point identification methods for FOPDT, DPPDT, and SOPDT models based on information taken from the reaction curve.

In the case of two-point procedures, the following methods and times can be considered: [12] (35–85%), [13] (28.3–63.2%), or [14] (33–70%).

In the case of three-point methods, the following references and corresponding times can be considered: [15] (2–70–90%) and (5–70–90%), [16] (15–45–75%) by Stark, and [17] (14–55–91%).

The 123c identification method [11] uses the times sets (25–75%) for FOPDT and DPPDT and (25–50–75%) for SOPDT models, respectively. This method has been recently extended in [18] for the identification of repeated-pole plus dead-time (RPPDT) model.

In the previously mentioned references, the location of such points is diverse: in some cases, no explanation is given for the selection; in others, significant points closely related to model parameters are selected. In some methods, see, e.g., [11], these points are not arbitrarily fixed but have been selected in order to optimize the identified model parameters. It is obvious that the accuracy of the identified model depends on the selection of these points on the process reaction curve, as will be discussed later.

Over the last decades, the emergence of fractional calculus has made possible a great deal of academic and industrial effort focused on obtaining methods for more accurate modelling and identification of real-world phenomena; see, for example, [19–21]. There is a wide variety of fractional-order modelling techniques based on the reaction curve in the technical literature, the most common approach being the one based on nonlinear optimization; see, for example, [22–25]. In these methods, the fractional-order model parameters are generally obtained by minimizing the error between the process reaction curve and the fractional-order model step response. These techniques require a higher computational effort compared to existing analytical methods, whose main characteristic is the simplicity of their application. However, there are not many analytical techniques for modelling fractional-order processes based on the process reaction curve. In [26], some strategies have been proposed in order to determine the parameters of a FFOPDT model by making use of the step response data. It combines numerical computation and graphical estimation. This reference can be considered as a pioneering work in the fractional-order case. Integral-based estimating methods are proposed in [27, 28], which are robust against the presence of measurement noise. These last two methods can be considered as an extension of the area methods existing in the classical case [3].

Even though the fractional model has been proven technologically superior, industrial adoption for fractional approach requires more analysis [29].

In addition, recent surveys have outlined fractional-order PID controllers as an emerging tool in the field of process control, with the major reason for its success being the intrinsic robustness they offer at a higher degree of freedom to operate and tune the parameters [30–33]. In this context, several recent applications of this control algorithm can be considered in [34–37]. In some of the existing methodologies of designing fractional-order PID controllers, a simple model of the process is utilized to tune the parameters of the intended controller [38–43].

According to all the above, obtaining a simple-structure fractional-order model for a process is of significant importance and would be very useful in practically designing integer- and fractional-order control systems. Since the physical interpretation of a process step response is straightforward and identification algorithms for integer-order models based on fitting several representative points on the process reaction curve are easy to apply, it can be considered natural to extend such methods for the fractional-order case. For that purpose, a general identification procedure for FFOPDT models has been conducted in this paper based on fitting three points in the transient process response to a step input. In the same framework, a simplification of the general identification procedure is proposed, considering that the points are selected symmetrically on the process reaction curve. Analytical expressions of the corresponding FFOPDT model parameters for several sets of symmetrical points are obtained in order to show the applicability of the proposed procedure and to get insight into the influence of the set of symmetrical points on the accuracy of the identified model.

This paper is organized as follows: Section 2 is devoted to present mathematical fundamentals and background. In Section 3, a general procedure for identifying a FFOPDT model from any three points on the process reaction curve is presented. In Section 4, a simplification of the general identification procedure is considered, where the points are selected symmetrically on the process reaction curve. Analytical expressions of the corresponding FFOPDT model parameters for several sets of symmetrical points are also obtained in this section. Some examples are provided in Section 5 to show the effectiveness and applicability of the proposed procedure in obtaining identification methods for determining FFOPDT model parameters in comparison with other integer- and fractional-order identification methods. Furthermore, results of several numerical simulations are also illustrated in this section in order to get insight into the influence of selection of the set of symmetrical points on the accuracy of the identified model. Conclusions and final remarks are presented in Section 6.

2. Theoretical Background

In this section, some basic concepts and definitions in fractional calculus are briefly discussed.

Fractional calculus is a generalization of integration and differentiation to noninteger fundamental operator ${}_{a}D{}_t^{\alpha }$. The continuous integrodifferential operator is defined as$$\begin{array}{}\text{(1)}& {}_{a}D{}_t^{\alpha }=\begin{array}{cc}\frac{{\text{d}}^{\alpha }}{\text{d}{t}^{\alpha }},& \alpha >0,\\ 1,& \alpha =0,\\ {\displaystyle {\int }_a^t}{d\tau }^{\alpha },& \alpha <0,\end{array}\end{array}$$where a and t are bounds of the operation and α ∈$\mathrm{ℝ}$.

There are several different definitions of fractional operators [44]. One of the most used definitions of the fractional integration is the Riemann–Liouville definition:$$\begin{array}{}\text{(2)}& {}_{0}I{}_t^{\alpha }ft={}_{0}D{}_t^{-\alpha }ft=\frac{1}{\Gamma \alpha }{\displaystyle {\int }_0^t{t-\tau }^{\alpha -1}f\tau \text{d}\tau },\end{array}$$where t ≥ 0, α ∈${\mathrm{ℝ}}^{+},$ and Γ(·) is the Gamma function [44].

From relation (2), the Riemann–Liouville definition of fractional derivative of the order α can be written in the following form:$$\begin{array}{}\text{(3)}& {}_{0}D{}_t^{\alpha }ft=\frac{1}{\Gamma n-\alpha }\frac{d^n}{{\text{dt}}^n}{\displaystyle {\int }_0^t{t-\tau }^{n-\alpha -1}f\tau d\tau },\end{array}$$where n−1<α<n, n∈${\mathbb{Z}}^{+}$. In definitions (2) and (3), the subscripts 0 and t are the limits of operation and known as the terminals of fractional integration and differentiation, respectively. For the sake of simplicity, ${}_{0}D{}_t^{-\alpha }$ is denoted by $D^{-\alpha }$ and ${}_{0}D{}_t^{\alpha }$ by $D^{\alpha }$.

The Laplace transform of the Riemann–Liouville-based fractional derivative is$$\begin{array}{}\text{(4)}& L{D}^{\alpha }ft=s^{\alpha }Lft-{\displaystyle \sum _{k=0}^{n-1}{s}^k{D}^{\alpha -k-1}}f0,\end{array}$$where n − 1 ≤ α < n, which for zero initial conditions is reduced to$$\begin{array}{}\text{(5)}& L{D}^{\alpha }ft=s^{\alpha }Lft.\end{array}$$

A general fractional-order system can be described by a fractional differential equation of the form$$\begin{array}{}\text{(6)}& {\displaystyle \sum _{k=0}^n{a}_k{D}^{{\alpha }_k}yt}={\displaystyle \sum _{k=0}^m{b}_k{D}^{{\beta }_k}ut}.\end{array}$$

The corresponding fractional-order transfer function of incommensurate real orders has the following form [44]:$$\begin{array}{}\text{(7)}& Gs=\frac{Q{s}^{{\beta }_k}}{P{s}^{{\alpha }_k}}=\frac{b_m{s}^{{\beta }_m}+b_{m-1}{s}^{{\beta }_{m-1}}+\cdots +b_0{s}^{{\beta }_0}}{a_n{s}^{{\alpha }_n}+a_{n-1}{s}^{{\alpha }_{n-1}}+\cdots +a_1{s}^{{\alpha }_1}+a_0{s}^{{\alpha }_0}},\end{array}$$where P($s^{{\alpha }_k}$) and Q($s^{{\beta }_k}$) have no common zeros, $a$_k (k = 0,…,n) and b_k(k = 0, …, m) are constants, and $\alpha$_k(k = 0, …, n) and ß_k(k = 0, …, m) are arbitrary real or rational numbers and without loss of generality, they can be arranged as α_n > α_n − 1 > · · · > α₀, and β_m > ß_m − 1 > · · · > ß₀.

The strictly proper transfer function G(s) given by (7) is BIBO stable if and only if P(s) has no root in {Re(s) ≥ 0} [45].

In particular case, when there is a real number α as the greatest common divisor of α_k, k = 1, …, n, and ß_k, k = 1, …, m, it is called the commensurate order. It holds that α_k = kα and β_k = kα, 0 < α < 1, $\forall kϵ\mathrm{\hspace{0.17em}\mathbb{Z}}$, and incommensurate order system (7) can also be rewritten in the commensurate form as follows:$$\begin{array}{}\text{(8)}& Gs=\frac{Q{s}^{\alpha }}{P{s}^{\alpha }}=\frac{{\displaystyle {\sum }_{k=0}^m{b}_k{s}^{k\alpha }}}{{\displaystyle {\sum }_{k=0}^n{a}_k{s}^{k\alpha }}}.\end{array}$$

It has been proven that the commensurate system G(s) brought in (8) is BIBO stable if all the roots of polynomial equation P(x) = 0, in which x = s^α are positioned out of the sector |arg(x)| ≤ απ/2 [45].

Considering n > m, the function G(s) becomes a proper rational function in the complex variable s^α, and if it is supposed that roots of P(x) = 0 are distinct, the partial fraction expansion of transfer function (8) can be written in the following general form:$$\begin{array}{}\text{(9)}& G{s}^{\alpha }={\displaystyle \sum _{i=1}^n\frac{r_i}{s^{\alpha }+{\lambda }_i}},\end{array}$$where $\lambda$_i, i = 1, …, n are the roots of P(x) = 0 and r_i, i = 1, …, n are the corresponding residues. Taking inverse Laplace transform from (9) results in the impulse response of G(s^α) which is given in [44].$$\begin{array}{}\text{(10)}& ht=L^{-1}{\displaystyle \sum _{i=1}^n\frac{r_i}{s^{\alpha }+{\lambda }_i}}={\displaystyle \sum _{i=1}^n{r}_i{t}^{\alpha -1}{E}_{\alpha ,\alpha }-{\lambda }_i{t}^{\alpha }},\end{array}$$where E_α,α(z) denotes the so-called two-parameter Mittag–Leffler function, which for an arbitrary value z is defined as$$\begin{array}{}\text{(11)}& E_{\alpha ,\beta }z={\displaystyle \sum _{r=0}^{\infty }\frac{z^r}{\Gamma \alpha r+\beta }}.\end{array}$$

Integrating the right-hand side of (10), the following step response of the transfer function G(s^α) is obtained:$$\begin{array}{}\text{(12)}& gt={\displaystyle \sum _{i=1}^n{r}_i\frac{E_{\alpha ,1}-{\lambda }_i{t}^{\alpha }-1}{{\lambda }_i}}.\end{array}$$

Each component of the step response $g$(t) in (12) converges to its final value in a similar way as function t^–α does, as has been shown in [41].

3. General Identification Method

A general identification procedure for identifying a fractional-order model from the process reaction curve is presented in this section.

The controlled processes considered have an S-shaped step response, and they can be well characterized by a FFOPDT model, for which the differential equation can be expressed as$$\begin{array}{}\text{(13)}& T·D^{\alpha }yt+yt=K·ut-L,\end{array}$$where initial conditions are generally taken as zero to obtain the standard FFOPDT transfer function model:$$\begin{array}{}\text{(14)}& Ps=\frac{K{e}^{-Ls}}{1+T{s}^{\alpha }},\end{array}$$where K is the process gain, T > 0 is the apparent time constant, L ≥ 0 is the apparent dead-time, and α is the fractional order of the model.

The set of parameters$$\begin{array}{}\text{(15)}& {\theta }_P=K,T,L,\alpha ,\end{array}$$represents the FFOPDT model parameters, which will be identified in this paper using information taken from the process reaction curve.

The standard FOPDT$$\begin{array}{}\text{(16)}& Ps=\frac{K{e}^{-Ls}}{1+Ts},\end{array}$$can be considered as a particular case of the FFOPDT model (14) with α = 1. The step responses of FFOPDT models for increasing values of α, from α = 0.2 to 1.8, are shown in Figure 1. The step response of the considered system for α = 1 is represented in the dashed line.

[figure(s) omitted; refer to PDF]

FOPDT model (16) has been broadly used in practice to capture the essential dynamic response of industrial processes for the purpose of control design [3]. In this paper, the set of parameters in (15) characterize the dynamic behaviour of the considered controlled process. The FFOPDT model (14) can be considered as a generalization of the classical FOPDT, and the relevance of this generalization has great implications in both the identification of dynamic processes as well as in the controller parameter design of dynamic feedback loops [46].

It is also useful to consider a few parameters to characterize process dynamics. T_ar is the average residence time, which can be defined in the context of a fractional-order model as$$\begin{array}{}\text{(17)}& T_{ar}\dot{=}\frac{{\displaystyle {\int }_0^{\infty }tgt\text{d}t}}{{\displaystyle {\int }_0^{\infty }gt\text{d}t}}=L+T^{1/\alpha },\end{array}$$where $g$(t) is the impulse response of the system. The average residence time is essentially a rough measure of how long it takes the input to have a significant influence on the output.

In the same context, the normalized dead-time τ, which has the property 0 ≤ τ ≤ 1, can also be extended for a FFOPDT model as$$\begin{array}{}\text{(18)}& \tau \dot{=}\frac{L}{T_{ar}}=\frac{L}{L+T^{1/\alpha }}.\end{array}$$

This parameter can be used to characterize the difficulty of controlling a process. Roughly speaking, processes with small τ are easy to control, and the difficulty in controlling the system increases as τ increases. For the particular case α = 1, T_ar and τ parameters correspond to their standard definition for a FOPDT model [3].

Although the step response of the considered fractional model can conveniently describe both monotonic or nonmonotonic behaviours depending on the fractional-order α, this identification procedure is only applied to a process having an S-shaped step response, since systems with essentially monotone step responses are very common in process control [3]. Furthermore, it is suggested in [4] that the class of processes where PID is suitable can be characterized as having essentially monotone step responses. One way to characterize such processes is to introduce the monotonicity index:$$\begin{array}{}\text{(19)}& {\alpha }_m\dot{=}\frac{{\displaystyle {\int }_0^{\infty }gt\text{d}t}}{{\displaystyle {\int }_0^{\infty }gt\text{d}t}},\end{array}$$where $g$(t) is the impulse response of the system. Systems with α = 1 have monotone step responses, and systems with α > 0.8 are considered essentially monotone.

Figure 2 shows the normalized step responses of FFOPDT model (14) for different values of the fractional-order α and different values of the normalized dead‐time τ. Notice that all curves intersect at one point t = T_ar because of the normalization.

[figure(s) omitted; refer to PDF]

3.1. Model Identification Using the Controlled Process Reaction Curve

A step signal u(t) with amplitude Δu will be considered as input and a signal y_α(t) with an amplitude Δy as the response of the system, as shown in Figure 3. FFOPDT model (14) response to a ∆u step input change is$$\begin{array}{}\text{(20)}& y_{\alpha }t=\begin{array}{c}0,0\le t<L,\\ K1-E_{\alpha ,1}-\frac{1}{T}{t-L}^{\alpha }\Delta u,& t\ge L,\end{array}\end{array}$$where E_α,β is the two-parameter Mittag–Leffler function defined in (11).

[figure(s) omitted; refer to PDF]

From process reaction curve (20), the gain is given by$$\begin{array}{}\text{(21)}& K=\frac{\Delta y}{\Delta u},\end{array}$$where ∆u is the amplitude of the input signal and ∆y is the total process output change, as shown in Figure 3.

The process output y_α(t) can be normalized to its final value ∆y = K·∆u and using the shifted and normalized time $\tau =1/T{t-L}^{\alpha }$, and equation (20) is reduced to the following expression:$$\begin{array}{}\text{(22)}& {\tilde{y}}_{\alpha }\tau =1-E_{\alpha ,1}-\tau ,\hspace{1em}\tau \ge 0.\end{array}$$

If ${\tilde{y}}_{\alpha }\tau$ is the normalized process output, ${\tau }_x$ is the normalized time required to reach some specific output normalized value ${\tilde{y}}_{\alpha }{\tau }_x$ which is between 0 and 1, or 0% and 100% of the process output total change. The value of τ_x can be easily obtained using (22).

Correspondingly, the time t_x required for the process output to reach x% of the process output total change is$$\begin{array}{}\text{(23)}& t_x=L+{{\tau }_{x}T}^{1/\alpha }.\end{array}$$

As the remaining FFOPDT model parameters {T, L, α} must be obtained, it is necessary to determine the times {t_x1, t_x2, t_x3} to reach three points {y_α(t_x1), y_α(t_x2), y_α(t_x3)} on the process reaction curve. Considering equation (23), the following equations set is defined:$$\begin{array}{}\text{(24)}& \begin{array}{c}{t}_{x1}=L+{{\tau }_{x1}T}^{1/\alpha },\\ t_{x2}=L+{{\tau }_{x2}T}^{1/\alpha },\\ t_{x3}=L+{{\tau }_{x3}T}^{1/\alpha }.\end{array}.\end{array}$$

The so-called ratio index Δ in [18] can be adapted to this fractional framework accordingly and can be used to determine the fractional order α:$$\begin{array}{}\text{(25)}& \Delta \dot{=}\frac{t_{x3}-t_{x1}}{t_{x2}-t_{x1}}=\frac{{\tau }_{x3}^{1/\alpha }-{\tau }_{x1}^{1/\alpha }}{{\tau }_{x2}^{1/\alpha }-{\tau }_{x1}^{1/\alpha }},\end{array}$$where τ_x1, τ_x2, and τ_x3 are normalized times and can be obtained using equation (22). Equation (25) can be expressed as a function relating the fractional order α and ratio index ∆, where α = f₁(Δ).

The time-based parameters {T, L} can be solved from (24) by considering two points, {y_α(t_x1), t_x1)}and {y_α(t_x3), t_x3)}, on the process reaction curve, for which the equivalent normalized points are {ỹ_α(τ_x1), τ_x1} and {ỹ_α(τ_x3), τ_x3}. Then, the expressions for these parameters are$$\begin{array}{}\text{(26)}& T=a^{\alpha }{t_{x3}-t_{x1}}^{\alpha },\text{(27)}& L=t_{x3}-{\tau }_{x3}^{1/\alpha }{T}^{1/\alpha },\end{array}$$where$$\begin{array}{}\text{(28)}& a=\frac{1}{{\tau }_{x3}^{1/\alpha }-{\tau }_{x1}^{1/\alpha }}.\end{array}$$

To conclude, the set of equations in (29) comprises the general expressions for determining the parameters of the FFOPDT model, θ_P = {K, T, L, α}, using the times required for the response to reach any three points on the reaction curve.$$\begin{array}{}\text{(29)}& \begin{array}{c}K=\frac{\Delta y}{\Delta u},\\ \alpha =f_1\Delta ,\\ T=f_2\alpha {t_{x3}-t_{x1}}^{\alpha },\\ L=\text{max}{t}_{x3}-f_3\alpha T^{1/\alpha },0,\end{array}\end{array}$$where f₁ depends on ∆, which is defined as a function of the generic three points in (25), f₂(α) = a^α and f₃(α) = τ_x3^1/α are functions that depend on the normalized times τ_x1 and τ_x3, and τ_x3, respectively, and a parameter is defined in (28). Notice in (29) how the values of T and L have a high dependence on the value of α, which has a significant influence on the shape of the response. This fact emphasizes the importance of determining the value of α parameter accurately, as will be discussed later.

In (29), α > 0 and T > 0 are fulfilled in a natural way, since τ_x1 < τ_x2 < τ_x3 and t_x1 < t_x2 < t_x3. Another condition that must be fulfilled in order to meet L ≥ 0 is$$\begin{array}{}\text{(30)}& t_{x3}\ge {{\tau }_{x3}T}^{1/\alpha }.\end{array}$$

4. Symmetrical Method

In the general development of the identification procedure presented in the previous section, it has been assumed that the required three points on the reaction curve could be any value. In this section, a simplification of the general method in which these points can be arbitrarily selected but located symmetrically on the response curve is presented.

Note that, as depicted in Figure 4, the central point will be located in the middle of the range and will correspond to the time needed to reach 50% of the process output total change on the reaction curve y_α(t₅₀) (t_x2 = t₅₀). The remaining points could be located arbitrarily but symmetrically placed with respect to the central point. One of the symmetrical points will be denoted as x, the other being 100‒x. This means that the times to be determined will now be t_x1 = t_x and t_x3 = t_100‒x, where t_x and t_100‒x denote the time needed to reach x% and (100‒x)% of the process output total change, respectively, where 0 < x < 50.

[figure(s) omitted; refer to PDF]

Adapting equations (29) with these symmetry considerations, the following expressions are obtained:$$\begin{array}{}\text{(31)}& \begin{array}{c}K=\frac{\Delta y}{\Delta u},\\ \alpha =f_1\Delta ,\\ T=f_2\alpha {t_{100-x}-t_x}^{\alpha },\\ L=\text{max}{t}_{100-x}-f_3\alpha T^{1/\alpha },0,\end{array}\end{array}$$where f₁ is the function that relates fractional order α with the times ratio ∆.

In addition, Δ is defined as a function of the symmetrical three-points as$$\begin{array}{}\text{(32)}& \Delta \dot{=}\frac{t_{100-x}-t_x}{t_{50}-t_x}=\frac{{{\tau }_{100-x}}^{1/\alpha }-{{\tau }_x}^{1/\alpha }}{{{\tau }_{50}}^{1/\alpha }-{{\tau }_x}^{1/\alpha }},\end{array}$$and it depends on the normalized times τ_x, τ₅₀, and τ_100‒x. f₂(α) = a^α, where$$\begin{array}{}\text{(33)}& a=\frac{1}{{{\tau }_{100-x}}^{1/\alpha }-{{\tau }_x}^{1/\alpha }},\end{array}$$and f₃ (α) = (τ_100‒x)^1/α are functions that depend on the normalized times τ_x and τ_100‒x, and τ_100‒x, respectively.

In this paper, the following sets of symmetrical points will be considered: (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%), which correspond to the following values of x = 5, 10, 15, 20, 25, and 30, respectively. The objective is to get insight into the influence of the selection of representative points on the accuracy of the obtained fractional-order model, as has been mentioned previously.

Tables 1–6 contain the values of the normalized times τ_x, τ₅₀, and τ_100‒x, for values of x = 5, 10, 15, 20, 25, and 30, respectively.

Table 1

Numeric values of normalized times τ₅, τ₅₀, and τ₉₅ for different values of α.

α	τ5 (s)	τ50 (s)	τ95 (s)
0.50	0.0021	0.5916	126.3400
0.55	0.0040	0.5990	69.3350
0.60	0.0060	0.6058	41.4000
0.65	0.0090	0.6150	26.3090
0.70	0.0130	0.6231	17.5400
0.75	0.0172	0.6340	12.1030
0.80	0.0225	0.6435	8.6100
0.85	0.0286	0.6560	6.2560
0.90	0.0355	0.6670	4.6920
0.95	0.0430	0.6810	3.6640
1.00	0.0513	0.6932	2.9980
1.05	0.0603	0.7080	2.5600
1.10	0.0699	0.7220	2.2633

Table 2

Numeric values of normalized times τ₁₀, τ₅₀, and τ₉₀ for different values of α.

α	τ10 (s)	τ50 (s)	τ90 (s)
0.50	0.0093	0.5916	30.8550
0.55	0.0142	0.5990	19.7430
0.60	0.0203	0.6058	13.4810
0.65	0.0276	0.6150	9.6760
0.70	0.0359	0.6231	7.2250
0.75	0.0453	0.6340	5.5750
0.80	0.0558	0.6435	4.4280
0.85	0.0670	0.6560	3.6140
0.90	0.0791	0.6670	3.0310
0.95	0.0919	0.6810	2.6100
1.00	0.1059	0.6932	2.3040
1.05	0.1195	0.7080	2.0810
1.10	0.1350	0.7220	1.9150

Table 3

Numeric values of normalized times τ₁₅, τ₅₀, and τ₈₅ for different values of α.

α	τ15 (s)	τ50 (s)	τ85 (s)
0.50	0.0229	0.5916	13.2100
0.55	0.0321	0.5990	9.3140
0.60	0.0428	0.6058	6.9304
0.65	0.0546	0.6150	5.3770
0.70	0.0677	0.6231	4.3117
0.75	0.0817	0.6340	3.5590
0.80	0.0966	0.6435	3.0097
0.85	0.1122	0.6560	2.6050
0.90	0.1285	0.6670	2.3004
0.95	0.1452	0.6810	2.0730
1.00	0.1626	0.6932	1.8973
1.05	0.1802	0.7080	1.7660
1.10	0.1984	0.7220	1.6626

Table 4

Numeric values of normalized times τ₂₀, τ₅₀, and τ₈₀ for different values of α.

α	τ20 (s)	τ50 (s)	τ80 (s)
0.50	0.0450	0.5916	7.0390
0.55	0.0588	0.5990	5.3510
0.60	0.0750	0.6058	4.2490
0.65	0.0907	0.6150	3.4920
0.70	0.1082	0.6231	2.9493
0.75	0.1261	0.6340	2.5520
0.80	0.1448	0.6435	2.2526
0.85	0.1639	0.6560	2.0260
0.90	0.1834	0.6670	1.8501
0.95	0.2031	0.6810	1.7160
1.00	0.2232	0.6932	1.6096
1.05	0.2434	0.7080	1.5290
1.10	0.2639	0.7220	1.4642

Table 5

Numeric values of normalized times τ₂₅, τ₅₀, and τ₇₅ for different values of α.

α	τ25 (s)	τ50 (s)	τ75 (s)
0.50	0.0780	0.5916	4.2100
0.55	0.0965	0.5990	3.4020
0.60	0.1170	0.6058	2.8480
0.65	0.1370	0.6150	2.4510
0.70	0.1590	0.6231	2.1580
0.75	0.1793	0.6340	1.9370
0.80	0.2021	0.6435	1.7662
0.85	0.2224	0.6560	1.6350
0.90	0.2450	0.6670	1.5310
0.95	0.2659	0.6810	1.4510
1.00	0.2880	0.6932	1.3870
1.05	0.3095	0.7080	1.3380
1.10	0.3314	0.7220	1.3000

Table 6

Numeric values of normalized times τ₃₀, τ₅₀, and τ₇₀ for different values of α.

α	τ30 (s)	τ50 (s)	τ70 (s)
0.50	0.1242	0.5916	2.6917
0.55	0.1478	0.5990	2.2940
0.60	0.1716	0.6058	2.0079
0.65	0.1952	0.6150	1.7980
0.70	0.2188	0.6231	1.6376
0.75	0.2422	0.6340	1.5160
0.80	0.2655	0.6435	1.4191
0.85	0.2884	0.6560	1.3450
0.90	0.3114	0.6670	1.2854
0.95	0.3341	0.6810	1.2410
1.00	0.3567	0.6932	1.2041
1.05	0.3792	0.7080	1.1780
1.10	0.4016	0.7220	1.1568

Hereafter, considering values of the corresponding normalized times, data sets {Δ, α}, {α, a^α}, and {α, (τ_100‒x)^1/α} are obtained for each one of the considered sets of symmetrical points. From these data, the analytical expressions for f₁(Δ), f₂(α), and f₃(α) are estimated using curve fitting. Finally, fractional-order model parameters {T, L, α} can be determined from expressions (31) and experimental values t_x, t₅₀, and t_100‒x collected from the reaction curve.

Δ values for different values of α, 0.5 ≤ α ≤ 1.1, have been obtained in Figure 5, according to expression (32) and using the normalized times τ_x, τ₅₀, and τ_100‒x, for x = 5, 10, 15, 20, 25, and 30, obtained from Tables 1–6. Note that Δ index represents the ratio between the time difference in reaching from x to (100‒x)% and from x to 50% of the total variation of the process output, as indicated in (32).

[figure(s) omitted; refer to PDF]

The data set {Δ, α} allows establishing a relationship between α and Δ in the form of an analytic function, α = f₁ (Δ).

Figure 5 shows the data sets {∆, α} obtained from equation (32) for different values of α, where 0.5 ≤ α ≤ 1.1, and their corresponding results of least-squares curve fitting for the different sets of symmetrical points. Note that Δ depends on normalized times {τ_x, τ₅₀, and τ_100‒x}, as indicated in (32), and has a strong dependence on α parameter.

Data is fitted using the Levenberg–Marquardt least-squares curve-fitting algorithm in all graphs in Figure 5, and the following rational function is considered in all cases:$$\begin{array}{}\text{(34)}& \alpha =f_1\Delta =\frac{p_1{\Delta }^2+p_2\Delta +p_3}{{\Delta }^2+q_1\Delta +q_2},\end{array}$$where Table 7 shows the values of the parameters {p_i, q_i} for function f₁(Δ) and each one of the selected sets of symmetrical points.

Table 7

Parameters {p_i, q_i} of the rational function f₁(Δ) for the different sets of symmetrical points.

5–50–95%	10–50–90%	15–50–85%	20–50–80%	25–50–75%	30–50–70%
p1 = 0.5786	p1 = 0.4882	p1 = 0.4767	p1 = 0.4613	p1 = 0.4424	p1 = 0.4174
p2 = 93.2	p2 = 109.6	p2 = 38.15	p2 = 17.04	p2 = 8.106	p2 = 3.942
p3 = 133.9	p3 = 201.7	p3 = 24.63	p3 = −0.0419	p3 = −3.652	p3 = −3.208
q1 = 114.3	q1 = 152.6	q1 = 52.91	q1 = 23.18	q1 = 10.40	q1 = 4.384
q2 = 26.59	q2 = 30.39	q2 = −29.59	q2 = −22.93	q2 = −14.00	q2 = −8.02

It is important to note that for increasing values of x, Δ range becomes shorter. f₁(Δ) is a nonlinear function whose slope is steeper for increasing values of x, indicating that the sensitivity is greater.

Figure 6 shows the data sets {α, a^α}, where a parameter is obtained from equation (33), for different values of α, 0.5 ≤ α ≤ 1.1, and their corresponding results of least-squares curve fitting for the different sets of symmetrical points. Note that f₂(α) = a^α depends on α and normalized times τ_x and τ_100‒x.

[figure(s) omitted; refer to PDF]

In Figure 6, it can be seen that the amplitude of f₂ increases for increasing values of x. On the other hand, for decreasing values of x, function f₂ presents a more pronounced curvature for small values of α.

Data is fitted using the Levenberg–Marquardt least-squares curve-fitting algorithm in all graphs in Figure 6, and the following rational function is used in all cases:$$\begin{array}{}\text{(35)}& f_2\alpha =\frac{p_1\alpha +p_2}{{\alpha }^2+q_1\alpha +q_2},\end{array}$$where the values of the parameters {p_i, q_i} for function f₂(α) and each one of the selected sets of symmetrical points are shown in Table 8.

Table 8

Parameters {p_i, q_i} of the rational function f₂(α) for the different sets of symmetrical points.

5–50–95%	10–50–90%	15–50–85%	20–50–80%	25–50–75%	30–50–70%
p1 = 0.06589	p1 = 0.2417	p1 = 0.9406	p1 = 9746	p1 = 1.715·104	p1 = 3.691·104
p2 = −0.03054	p2 = −0.0986	p2 = −0.3488	p2 = −3199	p2 = −3675	p2 = −426.5
q1 = −2.23	q1 = −2.338	q1 = −2.722	q1 = −5405	q1 = −1.074·104	q1 = −3.276·104
q2 = 1.333	q2 = 1.651	q2 = 2.745	q2 = 1.445·104	q2 = 2.544·104	q2 = 6.349·104

Figure 7 shows data sets {α, (τ_100−x)^1/α}, where normalized times τ_100‒x can be obtained from Tables 1−6, for different values of α, 0.5 ≤ α ≤ 1.1, and their corresponding results of least-squares curve fitting for the different sets of symmetrical points. Note that f₃(α) = (τ_100−x)^1/α depends on α and normalized time τ_100‒x. The range of amplitudes in f₃ is reduced for increasing values of x, as can be shown in Figure 7. Data is fitted using the Levenberg–Marquardt least-squares curve-fitting algorithm in all graphs in Figure 7, and the following rational function is used for all cases:$$\begin{array}{}\text{(36)}& f_3\alpha =\frac{p_1{\alpha }^2+p_2\alpha +p_3}{{\alpha }^2+q_1\alpha +q_2},\end{array}$$where the values of the parameters {p_i, q_i} for function f₃(α) and each one of the selected sets of symmetrical points are shown in Table 9.

[figure(s) omitted; refer to PDF]

Table 9

Parameters {p_i, q_i} of the rational function f₃(α) for the different sets of symmetrical points.

5–50–95%	10–50–90%	15–50–85%	20–50–80%	25–50–75%	30–50–70%
p1 = 5.492	p1 = 3.351	p1 = 1.7	p1 = 1.302	p1 = 1.243	p1 = 1.315
p2 = −11.05	p2 = −6.134	p2 = −2.584	p2 = −1.629	p2 = −1.368	p2 = −1.306
p3 = 6.092	p3 = 3.36	p3 = 1.42	p3 = 0.846	p3 = 0.666	p3 = 0.6238
q1 = −1.154	q1 = −0.9959	q1 = −0.9406	q1 = −0.8693	q1 = −0.7536	q1 = −0.5318
q2 = 0.3351	q2 = 0.2491	q2 = 0.2235	q2 = 0.1919	q2 = 0.1433	q2 = 0.05728

It is important to emphasize that estimated values of α, T, and L obtained from expressions (31) depend on functions f₁(Δ), f₂(α), and f₃(α). Thus, these functions play a significant role in the process identification procedure as the features of normalized step responses (22) can be well characterized due to their contribution. Note that, for any different choice of times set {t_x1, t_x2, t_x3} to reach three points {y_α(t_x1), y_α(t_x2), and y_α(t_x3)} on the reaction curve, the accuracy of the identification results only depends upon the fitting precision. In this context, an accurate determination of α-value is of primary importance since, subsequently, functions f₂ and f₃—and therefore T and L—depend on the estimated value of α.

In the curve-fitting process, it has been found that, in general, functions f₁(Δ), f₂(α), and f₃(α) for low values of x require more parameters than those with higher x values to be able to fit the data with the rational function. Although some of the data sets could have been adjusted with fewer parameters, for each function (f₁(Δ), f₂(α), and f₃(α)), the same expressions with the same number of parameters have been used for curve fitting, respectively. In this manner, comparison of estimated values of α, T, and L, which directly depend on functions f₁, f₂, and f₃ for different sets of symmetrical points, can be made in the same conditions. In this section, it has been verified that rational expressions used fit well data sets in Figures 5–7.

Remark 1.

Expressions (34)–(36), with parameters in Tables 7–9, for functions f₁, f₂, and f₃ have been obtained for 0.5 ≤ α ≤ 1.1. This range includes the dynamics for the majority of the representative processes encountered in process control, which can be characterized as having essentially monotone step responses; see, for example, the test batch considered in [4]. For values of α outside this range, the considered curve fittings are not valid. A significant fact observed is that the step response of a FFOPDT system for α values less than 0.5 becomes extremely slow, as can be seen in Figure 1. Note that values of the normalized time τ_100‒x increase as α decreases, especially in a very significant way for small values of x.

Although the step response of a FFOPDT model is overdamped only for α values between 0 and 1, the value 1.1 has also been included since it has been observed in numerous experiments that for some processes, generally lag-dominated, with high-order integer-order transfer functions, values of α are slightly greater than 1. In this way, a better fit of the transient response of the model is obtained in relation to the reaction curve of the process. Specifically, the behavior of the step response of a FFOPDT system when α = 1.1 can be observed in Figure 2.

4.1. Algorithm for Determining FFOPDT Model Parameters

To facilitate software implementation of the proposed identification method, an algorithm is developed in Table 10. In this method, the variation in input signal Δu and process output Δy and times necessary to reach x% (t_x), 50% (t₅₀), and (100‒x)% (t_100‒x) of the total change of the process output on the reaction curve must be collected in order to determine FFOPDT model parameters θ_P = {K, T, L, α}.

Table 10

Algorithm for determining FFOPDT model parameters θ_P = {K, T, L, α} when three symmetrical points {t_x, t₅₀, t_100‒x} on the reaction curve are specified.

Algorithm: identification method for three symmetrical points

Input: {tx, t50 and t100‒x}, Δy, and Δu from the process reaction curve

Output: FFOPDT model parameters θP = {K, T, L, α}

1: Calculate the process gain K using (21).

2: Collect tx, t50, and t100‒x from the process reaction curve.

3: Calculate the value of Δ using equation (32).

4: Calculate the value of α = f1(Δ) using equation (34).

5: Calculate the value of functions f2(α) and f3(α) using equations (35) and (36), respectively.

6: Determine the value of T = f2(α)·(t100‒x – tx)α.

7: Determine the value of L = t100‒x – f3(α)·T1/α.

5. Illustrative Examples

In this section, the proposed identification method for FFOPDT models was tested for several sets of symmetrical points on the reaction curve and compared with other identification methods for integer- and fractional-order models.

The following process models were selected in order to show the effectiveness of the proposed procedure in finding an approximated four-parameter model (14) and to get insight into the influence of the selected set of symmetrical points on the accuracy of the identified model:$$\begin{array}{}\text{(37)}& P_{1,i}s=\frac{K_1}{{1+T_1{s}^{{\lambda }_i}}^2}{e}^{-L_{1}s},\text{(38)}& P_{2}s=\frac{K_2}{{1+T_2{s}^{{\lambda }_2}}^n}.\end{array}$$

On the one hand, processes P_1,i are used to compare the proposed identification procedure for several symmetrical points with some two- and three-point identification methods for FOPDT, DPPDT, and SOPDT models and to get insight into the influence of the location of symmetrical points on the accuracy of the fractional-order identified model.

On the other hand, process P₂ is used to compare the model performance of the method proposed in this paper with other integer- and fractional-order methods.

A step-change input signal was applied to these processes in order to register the process reaction curve. The output responses of the processes were then used to calculate the parameters of the corresponding models. Note that the sample period used in all the experiments was set to T_S = 10 ms.

Given the experimental data from an identification test, it is necessary to verify the effectiveness of the model structure adopted for fitting and the accuracy of the corresponding model parameters. A number of fitting objective functions and model validation methods using a step excitation signal for system identification have been presented in the literature [8].

Without loss of generality, the following time-domain fitting criterion was used to evaluate the performance of the model identification:$$\begin{array}{}\text{(39)}& S\overline{\theta }=\frac{1}{N_s}{\displaystyle \sum _{k=1}^{N_s}{ek{T}_s,\overline{\theta }}^2}=\frac{1}{N_s}{\displaystyle \sum _{k=1}^{N_s}{yk{T}_s-y_{m}k{T}_s,\overline{\theta }}^2},\end{array}$$where $\overline{\theta }$ is the vector of process model parameters, $ek{T}_s,\overline{\theta }$ is the open-loop step response difference between the actual process and the identified process model output signal, $yk{T}_s$ and $y_{m}k{T}_s,\overline{\theta }$, respectively, T_S is the sampling period, and N_ST_S is the time length of the dynamic (transient) response. Although in the context of a step response test, N_ST_S may be taken as the settling time that is usually defined as the time to move into an error band of 2% or 5% with respect to the final steady-state output deviation in response to a step change, the upper limit of the time interval is the instant in which the response has reached its final value.

This time-domain criterion is the Mean Squared Error (MSE) and gives a measure of the average of the squares of errors. How accurately a model reproduces the process reaction curve is evaluated with this index.

In this context, a relative performance index is more important than their absolute values since the goodness of an identification method with respect to another can be quantified.

In this paper, the following relative performance index is used:$$\begin{array}{}\text{(40)}& S_{x,y}{\overline{\theta }}_x,{\overline{\theta }}_y=\frac{S_x{\overline{\theta }}_x}{S_y{\overline{\theta }}_y},\end{array}$$where $S_x{\overline{\theta }}_x$ and $S_y{\overline{\theta }}_y$ are the MSE values in (39) and ${\overline{\theta }}_x$ and ${\overline{\theta }}_y$ are the vectors of process model parameters for x and y models, respectively.

The simulation of these illustrative examples has been implemented using MATLAB and FOTF (Fractional Order Transfer Function) toolbox. FOTF is a control toolbox for fractional-order systems developed by Xue et al. [47] which extends many MATLAB built-in functions to deal with fractional-order models. The reference text for the FOTF toolbox is [48], where the author explains thoroughly all its commands and applications.

5.1. Example 1

The fractional-order process model (37), where K₁ = 1, T₁ = 1 s, L₁ = 0.1 s, and λ_i = 0.60, 0.65, …, 1.00, is considered in this example. Model (37) is a fractional second-order process with dead time (FSOPDT) to provide some modelling error or deviation from FFOPDT dynamics, which is the model structure selected in the proposed identification method.

The procedure that has been followed with this example is as follows:

(1) Nine process plants P_1,i, where i = 1, …, 9, and λ_i = 0.60, 0.65, …, 1.00, are considered.

(2) The FFOPDT model parameters ${\overline{\theta }}_{ij}$ for each process P_1,i are obtained using the proposed identification method for each of the sets of points. ${\overline{\theta }}_{ij}$ is expressed as follows:$$\begin{array}{}\text{(41)}& {\overline{\theta }}_{ij}=K_{ij},T_{ij},L_{ij},{\alpha }_{ij},\end{array}$$

  where i = 1, …, 9 represents the value of λ_i process parameter, as indicated in Table 11, and j = 1, …, 6 represents the identification method employed, as indicated in Table 12. In short, a set of 54 process model parameters is obtained.

(3) The step responses for each of the 54 models are obtained, and their respective values of the model performance index $S{\overline{\theta }}_{ij}$ are determined.

Table 11

Set of process fractional orders, sampling periods, number of samples, and time length used in FFOPDT step responses for processes P_1,i (i = 1, …, 9).

i	λi	TS (s)	NS	NSTS (s)
1	0.60	0.01	58278	582.77
2	0.65	0.01	29312	293.12
3	0.70	0.01	15816	158.15
4	0.75	0.01	8960	89.59
5	0.80	0.01	5224	52.23
6	0.85	0.01	3067	30.66
7	0.90	0.01	1760	17.60
8	0.95	0.01	965	9.64
9	1.00	0.01	594	5.93

Table 12

Set of symmetrical points used for FFOPDT model identification.

j	Method#	Set of points
1	1	5–50–95%
2	2	10–50–90%
3	3	15–50–85%
4	4	20–50–80%
5	5	25–50–75%
6	6	30–50–70%

This batch of processes is used, on the one hand, to validate the proposed identification method for the different sets of symmetrical points and, on the other hand, to get insight into the influence of the selection of representative points on the accuracy of the obtained fractional-order model.

Note that Table 11 contains sampling period, number of samples, and time length that have been considered in step responses for each of the processes. Table 12 lists the set of points considered for the identification of the FFOPDT models. Figure 8 illustrates the values of the model performance index $S{\overline{\theta }}_{ij}$ that have been obtained with the proposed procedure for each process (i = 1, …, 9) and using the different methods (j = 1, …, 6).

[figure(s) omitted; refer to PDF]

Note that the values $S{\overline{\theta }}_{ij}$ in Figure 8 decrease substantially when they move from method ^#6 to ^#1. This behavior occurs for all the processes considered, except for those where λ_i is close to 1, that is, when the process approaches one of the integer order. This behaviour has already been observed in the technical literature. The 123c method [11] uses time sets (25–75%) to identify FOPDT and DPPDT models and (25–50–75%) for SOPDT models. In [18], points (25–50–75%) are proposed to identify RPPDT models. In these two papers, these points have not been chosen arbitrarily but rather those that optimize the identified model parameters have been selected, i.e., FOPDT, DPPDT, SOPDT, and RPPDT models, respectively. In Figure 8, the lowest value of $S{\overline{\theta }}_{ij}$ for P_1,9, which is an integer-order model, corresponds to method ^#5 (25–50–75%).

Figure 8 also shows that the proposed identification method gives good fit between the corresponding identified models and the actual process reaction curves. This is confirmed by the low values in time-domain model performance indices $S{\overline{\theta }}_{ij}$ for all processes P_1,i.

From the results obtained in Figure 8, it can be observed that method ^#1 is the one that provides the best results and, therefore, the set of points (5–50–95%) is the one that is recommended to be selected when using this symmetrical identification procedure. In this regard, relative performance indices between method ^#6 (j = 6) and method ^#1 (j = 1), $S_{\mathrm{6,1}}{\overline{\theta }}_{i6},{\overline{\theta }}_{i1}$, for i = 1, …, 9, are represented in Figure 9, where it can be noticed that method ^#1 is approximately 2.0, 2.25, 2.5, 3.5, 4.0, 2.75, and 1.5 times better in terms of S than the one obtained with method ^#6, for i = 1, …, 7, respectively.

[figure(s) omitted; refer to PDF]

The case λ₅ = 0.8 for process P_1,5 is considered among processes P_1,i. The process reaction curve for this model is shown in Figure 10. From this data, the information is summarized in Table 13.

[figure(s) omitted; refer to PDF]

Table 13

Process information for fractional-order model identification of process P_1,5.

Method #1: 5–50–95%	Method #2: 10–50–90%	Method #3: 15–50–85%	Method #4: 20–50–80%	Method #5: 25–50–75%	Method #6: 30–50–70%
Δu = 1
Δy = 1
t5 = 0.3310 s	t10 = 0.4890 s	t15 = 0.6410 s	t20 = 0.7940 s	t25 = 0.9540 s	t30 = 1.1250 s
t50 = 1.9900
t95 = 19.123 s	t90 = 9.7230 s	t85 = 6.7390 s	t80 = 5.2090 s	t75 = 4.2430 s	t70 = 3.5600 s

Considering data from Table 13, the FFOPDT process parameters ${\overline{\theta }}_{i5}=K_{i5},T_{i5},L_{i5},{\alpha }_{i5}$, where i = 1, …, 6, have been obtained in Table 14 for each one of the proposed sets of symmetrical points in Table 12.

Table 14

FFOPDT model parameters for the considered sets of representative points (5–50-95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%), respectively.

Method #1: 5–50–95%	Method #2: 10–50–90%	Method #3: 15–50–85%	Method #4: 20–50–80%	Method #5: 25–50–75%	Method #6: 30–50–70%
K1 = 1.0000	K2 = 1.0000	K3 = 1.0000	K4 = 1.0000	K5 = 1.0000	K6 = 1.0000
T1 = 2.3318 s	T2 = 2.3011 s	T3 = 2.2834 s	T4 = 2.2688 s	T5 = 2.2666 s	T6 = 2.2418 s
L1 = 0.3742 s	L2 = 0.4724 s	L3 = 0.4041 s	L4 = 0.4057 s	L5 = 0.3901 s	L6 = 0.4175 s
α1 = 0.8719	α2 = 0.8882	α3 = 0.8992	α4 = 0.9067	α5 = 0.9122	α6 = 0.9143

The proposed identification method will be compared with various two- and three-point identification methods for FOPDT, DPPDT, and SOPDT models based on information also obtained from the reaction curve.

In the case of two-point methods [11,14], both FOPDT and DPPDT models are used.

Table 15 contains the parameters of the FOPDT model identified using the two-point identification methods proposed by Alfaro and by Viterková, respectively. Table 16 shows the parameters of the DPPDT model obtained using the same methods.

Table 15

FOPDT model parameters obtained using Alfaro’s and Viterková’s two-point methods, respectively.

Method #7: Alfaro 1 (25–75%)	Method #8: Viterková 1 (33–70%)
K7 = 1.0000	K8 = 1.0000
T7 = 2.9936 s	T8 = 2.8959 s
L7 = 0.0923 s	L8 = 0.0757 s

Table 16

DPPDT model parameters obtained using Alfaro’s and Viterková’s two-point methods, respectively.

Method #9: Alfaro 2 (25–75%)	Method #10: Viterková 2 (33–70%)
K9 = 1.0000	K10 = 1.0000
T9 = 1.8997 s	T10 = 1.8468 s
L9 = 0 s	L10 = 0 s

In the case of three-point methods in [15] and Stark in [16], the SOPDT model is considered.

Table 17 contains the parameters of the SOPDT model obtained using the three-point identification methods proposed by Jahanmiri and Fallahi and by Stark, respectively.

Table 17

SOPDT model parameters obtained using Stark’s and Jahanmiri and Fallahi’s three-point methods, respectively.

Method #11: Stark (15–45–75%)	Method #12: Jahanmiri and Fallahi 1 (2–70–90%)	Method #13: Jahanmiri and Fallahi 2 (5–70–90%)
K11 = 1.0000	K12 = 1.0000	K13 = 1.0000
T11a = 2.8764 s	T12a = 3.9984 s	T13a = 4.0453 s
T11b = 0.3714 s	T12b = 0.0362 s	T13b = 0.0362 s
L11 = 0 s	L12 = 0.2220 s	L13 = 0.3310 s

Figures 11–14 compare the process reaction curve with the step responses of the FFOPDT models obtained with the proposed method for different sets of symmetrical points and those of the FOPDT, DPPDT, and SOPDT models proposed by Alfaro, Viterková, Stark, and Jahanmiri and Fallahi, respectively. The step responses have been grouped in graphs based on the similarity in the points on the reaction curve that have been used in each of the identification methods.

[figure(s) omitted; refer to PDF]

Table 18 shows the values of the time-domain model performance index $S\overline{\theta }$ for the different identification methods applied to process P_1,5. The relative performance index values $S_{x,y}{\overline{\theta }}_x,{\overline{\theta }}_y$, which allow the comparison between different identification methods, are listed in Table 19.

Table 18

Time-domain model performance indices for process P_1,5.

Identification method#	Set of points	Model	S
#1: proposed	(5–50–95%)	FFOPDT	S1 = 7.78·10−5
#2: proposed	(10–50–90%)	FFOPDT	S2 = 1.51·10−4
#3: proposed	(15–50–85%)	FFOPDT	S3 = 1.79·10−4
#4: proposed	(20–50–70%)	FFOPDT	S4 = 2.18·10−4
#5: proposed	(25–50–75%)	FFOPDT	S5 = 2.47·10−4
#6: proposed	(30–50–70%)	FFOPDT	S6 = 2.69·10−4
#7: Alfaro 1	(25–75%)	FOPDT	S7 = 5.98·10−4
#8: Viterková 1	(33–70%)	FOPDT	S8 = 6.30·10−4
#9: Alfaro 2	(25–75%)	DPPDT	S9 = 1.40·10−3
#10: Viterková 2	(33–70%)	DPPDT	S10 = 1.30·10−3
#11: Stark	(15–45–70%)	SOPDT	S11 = 6.66·10−4
#12: Jahanmiri and Fallahi 1	(2–70–90%)	SOPDT	S12 = 9.25·10−4
#13: Jahanmiri and Fallahi 2	(5–70–90%)	SOPDT	S13 = 1.10·10−3

Table 19

Relative performance index-based comparison between some considered identification methods for process P_1,5.

Compared methods	Sx,y = Sx/Sy
Method #2-method #1	S2,1 = S2/S1 = 1.94
Method #3-method #1	S3,1 = S3/S1 = 2.30
Method #4-method #1	S4,1 = S4/S1 = 2.81
Method #5-method #1	S5,1 = S5/S1 = 3.17
Method #6-method #1	S6,1 = S6/S1 = 3.45
Method #7-method #1	S7,1 = S7/S1 = 7.69
Method #7-method #5	S7,5 = S7/S5 = 2.42
Method #8-method #1	S8,1 = S8/S1 = 8.10
Method #8-method #6	S8,6 = S8/S6 = 2.35
Method #9-method #1	S9,1 = S9/S1 = 17.76
Method #9-method #5	S9,5 = S9/S5 = 5.60
Method #10-method #1	S10,1 = S10/S1 = 16.75
Method #10-method #6	S10,6 = S10/S6 = 4.85
Method #11-method #1	S11,1 = S11/S1 = 8.57
Method #11-method #3	S11,3 = S11/S3 = 3.72
Method #12-method #1	S12,1 = S12/S1 = 11.89
Method #12-method #2	S12,2 = S12/S2 = 6.12
Method #13-method #1	S13,1 = S13/S1 = 14.07
Method #13-method #2	S13,2 = S13/S2 = 7.25

Figures 11–14 illustrate that the step response of the models identified with the proposed method for the different sets of representative points give a good fit with the process reaction curve, particularly in the interval [x − (100 ‒ x)].

For smaller values of x, the interval [x − (100 ‒ x)] is larger and, therefore, the step response for these models fits better with the process reaction curve, which translates into a lower value in the model performance index S, as can be seen in Table 18.

This is mainly due to the fact that the value of α identified for smaller values of x is closer to its optimal value.

Note that α value of the model has a great influence on the shape of the step response for a FFOPDT model, as can be seen in Figure 2.

Regarding two-point methods, Figures 13 and 14 show step responses of the models obtained by Alfaro’s (25–75%) and Viterková’s (33–70%) methods, both for FOPDT and DPPDT, respectively, compared to those of the proposed method for (25–50–75%) and (30–50–70%).

Figures 11 and 12 show that Alfaro’s and Viterková’s methods fit the reaction curve very well in the ranges (25–75%) and (33–70%), respectively. However, the behaviour of both models outside these ranges deviates from the reaction curve, resulting in a higher S-value.

The value of S for method ^#1 is 7.7 and 8.1 times lower than the ones proposed by Alfaro and Viterková for FOPDT models, respectively. In turn, the value of S for methods ^#5 and ^#6 is 2.4 and 2.3 times lower than methods proposed by Alfaro and Viterková for FOPDT models, respectively.

Method ^#1 is 17.8 and 16.7 times better in terms of S than Alfaro’s and Viterková’s methods for DPPDT models, respectively. On the other hand, methods ^#5 and ^#6 are 5.6 and 4.8 times better than those of Alfaro and Viterková for DPPDT models. The responses of DPPDT models for Alfaro and Viterková have higher values in the model performance index shown in Table 18 in comparison with the same methods for FOPDT models.

Regarding three-point methods, Figure 11 shows step responses of Jahanmiri and Fallahi’s methods (2–70–90%) and (5–70–90%) for a SOPDT model, in comparison with those of the proposed method for (5–50–95%) and (10–50–90%).

The value of S for method ^#1 is 11.9 and 14.1 times lower than the one for methods ^#12 and ^#13, respectively. If both three-point methods are compared with method ^#2, the latter is 6.1 and 7.2 times better than methods ^#12 and ^#13, respectively.

Figure 12 shows the step responses of Stark’s method (15–45–75%) for a SOPDT model, in comparison with those of the proposed method for (15–50–85%) and (20–50–80%). The value of S for methods ^#1 and ^#3 is 8.6 and 3.7 times lower than the one for method ^#11.

There is an improvement in terms of S of 1.9, 2.3, 2.8, 3.2, and 3.4 times for process P₁ in using method ^#1 (5–50–95%) with respect to methods ^#2, ^#3, ^#4, ^#5, and ^#6, respectively, as illustrated in Table 19.

5.2. Example 2

In this example, the proposed identification method is compared with another identification method for the FFOPDT model.

The higher-order lag-dominated fractional-order process model (38), where K₂ = 2, T₂ = 1 s, λ₂ = 0.85, and n = 5, is considered in this example as proposed in [26].

The process reaction curve for this model is shown in Figure 15. From this data, the process information needed for model identification using the proposed procedure in different sets of times is summarized in Table 20.

[figure(s) omitted; refer to PDF]

Table 20

Process information for fractional-order model identification of process P₂.

Method #1: 5–50–95%	Method #2: 10–50–90%	Method #3: 15–50–85%	Method #4: 20–50–80%	Method #5: 25–50–75%	Method #6: 30–50–70%
Δu = 1
Δy = 2
t5 = 1.7910 s	t10 = 2.3450 s	t15 = 2.8000 s	t20 = 3.2180 s	t25 = 3.6240 s	t30 = 4.0320 s
t50 = 5.9370 s
t95 = 34.790 s	t90 = 19.200 s	t85 = 14.170 s	t80 = 11.580 s	t75 = 9.9300 s	t70 = 8.7400 s

Considering data from Table 20, the FFOPDT process parameters ${\overline{\theta }}_i=K_i,T_i,L_i,{\alpha }_i$, where i = 1, …, 6, have been obtained in Table 21 for each one of the proposed sets of representative points on the reaction curve, i.e., (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%), respectively.

Table 21

FFOPDT model parameters for the considered sets of representative points (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%), respectively.

Method #1: 5–50–95%	Method #2: 10–50–90%	Method #3: 15–50–85%	Method #4: 20–50–80%	Method #5: 25–50–75%	Method #6: 30–50–70%
K1 = 2.0000	K2 = 2.0000	K3 = 2.0000	K4 = 2.0000	K5 = 2.0000	K6 = 2.0000
T1 = 5.3737 s	T2 = 5.4182 s	T3 = 5.3090 s	T4 = 5.2366 s	T5 = 5.2330 s	T6 = 5.1383 s
L1 = 1.9617 s	L2 = 2.1230 s	L3 = 2.0287 s	L4 = 2.0965 s	L5 = 2.0972 s	L6 = 2.1926 s
α1 = 0.9099	α2 = 0.9397	α3 = 0.9555	α4 = 0.9656	α5 = 0.9741	α6 = 0.9767

Process (38) is approximated by a FFOPDT model following three different strategies proposed in [26], where model process parameters ${\overline{\theta }}_7-{\overline{\theta }}_9$ are given in Table 22.

Table 22

FFOPDT model parameters for the different strategies proposed in [26].

Method #7: Tavakoli–Kakhki 1	Method #8: Tavakoli–Kakhki 2	Method #9: Tavakoli–Kakhki 3
K7 = 2.00	K8 = 2.00	K9 = 2.00
T7 = 5.00 s	T8 = 5.00 s	T9 = 4.48 s
L7 = 1.50 s	L8 = 0.69 s	L9 = 1.50 s
α7 = 0.85	α8 = 0.85	α9 = 0.85

The process reaction curve and the corresponding step responses of the considered approximated models for (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%), respectively, are compared and illustrated in Figures 16–21. Moreover, the step responses for FFOPDT models proposed in [26] for process P₂ in comparison with the process reaction curve are shown in Figure 22.

[figure(s) omitted; refer to PDF]

Table 23 shows the values of the time-domain model performance index $S\overline{\theta }$ for the different identification methods applied to process P₂. The relative performance index values $S_{x,y}{\overline{\theta }}_x,{\overline{\theta }}_y$, which show the comparison between different identification methods, are listed in Table 24.

Table 23

Time-domain model performance indices for process P₂.

Identification method#	Set of points	Model	S
#1: proposed	(5–50–95%)	FFOPDT	S1 = 5.3·10−4
#2: proposed	(10–50–90%)	FFOPDT	S2 = 7.4·10−4
#3: proposed	(15–50–85%)	FFOPDT	S3 = 9.1·10−4
#4: proposed	(20–50–70%)	FFOPDT	S4 = 1.2·10−3
#5: proposed	(25–50–75%)	FFOPDT	S5 = 1.4·10−3
#6: proposed	(30–50–70%)	FFOPDT	S6 = 1.5·10−3
#7: Tavakoli–Kakhki’s method 1	—	FFOPDT	S7 = 9.5·10−4
#8: Tavakoli–Kakhki’s method 2	—	FFOPDT	S8 = 1.6·10−3
#9: Tavakoli–Kakhki’s method 3	—	FFOPDT	S9 = 5.5·10−4

Table 24

Relative performance index-based comparison between some considered fractional-order model identification methods for process P₂.

Compared methods	Sxy = Sx/Sy
Method #2-method #1	S21 = S2/S1 = 1.40
Method #3-method #1	S31 = S3/S1 = 1.74
Method #4-method #1	S41 = S4/S1 = 2.19
Method #5-method #1	S51 = S5/S1 = 2.57
Method #6-method #1	S61 = S6/S1 = 2.80
Method #7-method #1	S71 = S7/S1 = 1.81
Method #8-method #1	S81 = S8/S1 = 2.97
Method #9-method #1	S91 = S9/S1 = 1.06

In this example, a high-order fractional-order lag-dominated process has been used in order to illustrate the applicability and effectivity of the proposed identification procedure for identifying FFOPDT models.

Figures 16–21 show that the proposed method gives good results in the interval [x − (100 − x)] for all the considered sets of representative points. We need to recall that the method tries to place three symmetric points (x − 50 − (100 − x)%) on the reaction curve.

The location of the representative points on the reaction curve is very good for higher values of x, with a slight deviation for x = 5, as can be seen in Figures 16–21.

However, the obtained model fits the reaction curve better for low values of x, since α-value obtained is closer to its optimal value, as has been commented previously. This results in lower values of S for low values of x in the considered sets of representative points, as shown in Table 23.

There is an improvement in S of 1.4, 1.7, 2.2, 2.6, and 2.8 times in using method ^#1 (5–50–95%) with respect to methods ^#2, ^#3, ^#4, ^#5, and ^#6, respectively, as illustrated in Table 24.

Figure 22 shows that Tavakoli–Kakhki’s methods give similar results in terms of S to some of the proposed methods, and this can be quantified in Table 23.

Method ^#1 is 1.8, 3.0, and 1.1 times better in terms of S than Tavakoli–Kakhki’s methods ^#1, ^#2, and ^#3, respectively.

Furthermore, the proposed method not only gives better results than Tavakoli–Kakhki’s method in terms of S, but in the authors' opinion, it is easier to apply.

Remark 2.

In the second example, the comparison is restricted to analytical methods for FFOPDT models with a simplicity and computational effort similar to the one proposed in this paper. However, there are not many analytical techniques for modelling fractional-order processes based on the process reaction curve. Therefore, the comparison of the proposed symmetrical method is done with the method in [26], which is a well-recognized reference and where three strategies have been proposed in order to determine the parameters of a FFOPDT model by making use of the step response data. It combines numerical computation and graphical estimation.

In the introduction, there is a review of fractional-order model identification methods that are based on the process reaction curve; however, all of them use optimization to obtain the model parameters. In these cases, the accuracy of the model is improved at the cost of increased computational effort. In the authors' opinion, the comparison of the proposed method with these ones using optimization would not be made under equal conditions.

Remark 3.

Since processes are generally nonlinear, the dynamic characteristics of the FFOPDT model—gain, time constant, dead time, and fractional order—will vary when the operating point of the control system changes, due to a modification of the set point or due to the effect of disturbances.

The implicit uncertainty in the nominal model must therefore be taken into consideration.

There are generally two approaches in the technical literature when considering the parametric uncertainty of the plant in an identification method.

The first approach is to consider the uncertainty explicitly in the identification process. This generally makes the identification procedure more complicated.

The second one involves considering possible model uncertainties and variations in the dynamic characteristics of the controlled process in the design of the control system; see, e.g., [3] for integer-order controllers and [45] for the fractional case. A usual way of applying this second approach is to guarantee a certain degree of robustness—relative stability—of the control system designed to ensure its stability in the presence of variations in the process characteristics.

In this paper, a procedure for the identification of a reduced-order fractional model of the FFOPDT type is presented, where the main use of the identified model is for control purposes. Accordingly, the uncertainty in the parameters of the identified model will be considered in the design of the control system.

In the industrial context, large process industries have hundreds or thousands of control loops. For this reason, simplicity is a fundamental characteristic when identifying a process model for control purposes. Considering this second approach, it is possible to maintain the simplicity of the identification procedure.

6. Conclusions

In this paper, a general procedure to identify a FFOPDT model for industrial processes having S-shaped step responses and based on fitting three arbitrary points on the process reaction curve has been proposed. A simplification of the general identification procedure has also been considered, where only points symmetrically located on the reaction curve have been selected. The symmetrical method provides an efficient way to obtain the parameters of the model, by requiring the selection of the optimal location of only one of the points (x), since the other is immediately established by the requirement of symmetry.

Some numerical examples have been used in order to show the effectiveness and applicability of this procedure for the identification of fractional-order models and to get insight into the influence of selection of the set of symmetrical points on the accuracy of the identified model.

The proposed symmetrical procedure has been applied to the following sets of representative points: (5–50–95%), (10–50–90%), (15–50–85%), (20–50–80%), (25–50–75%), and (30–50–70%). The results of this paper verify that the accuracy of the identified fractional-order model is sensitive to the selection of the set of symmetrical points on the reaction curve and it has been discussed that a more accurate model is obtained for low values of x. New insights into this selection of set of points have been offered in the context of the proposed symmetrical procedure.

This identification procedure gives good results in comparison with other integer- and fractional-order identification methods.

Another aspect to highlight is that the proposed method is analytical, which facilitates its applicability in terms of a lower computational effort compared to complicated identification algorithms generally based on optimization. In the industrial context, large process industries typically have hundreds or thousands of control loops. For this reason, simplicity is of primary interest when identifying a process model for control purposes.

It is the opinion of the authors that this type of identification methods, where simplicity is emphasized, will encourage their industrial use and will help in bridging the gap between theoretical research on fractional models and its practical application in the process industry.

This expectation is the primary motivation of the present study.

Acknowledgments

The authors would like to thank the Basque Government for its partial funding support through the TRUSTIND ELKARTEK R & D project (ref. KK-2020/00054).

References

[1] T. Liu, Q. G. Wang, H. P. Huang, "A tutorial review on process identification from step or relay feedback test," Journal of Process Control, vol. 23 no. 10, pp. 1597-1623, DOI: 10.1016/j.jprocont.2013.08.003, 2013.

[2] G. F. Franklin, J. D. Powell, A. F. Emami-Naeini, Feedback Control of Dynamic Systems, 2019.

[3] K. J. Åström, T. Hägglund, Advanced PID Control, International Society of Automation, 2006.

[4] K. J. Åström, T. Hägglund, "Revisiting the Ziegler-Nichols step response method for PID control," Journal of Process Control, vol. 14 no. 6, pp. 635-650, DOI: 10.1016/j.jprocont.2004.01.002, 2004.

[5] H. P. Huang, J. C. Jeng, "Process reaction curve and relay methods identification and PID tuning," PID Control: New Identification and Design Methods,DOI: 10.1007/1-84628-148-2, 2005.

[6] A. Visioli, "Practical PID control," Advances in Industrial Control Series,DOI: 10.1007/1-84628-586-0, 2006.

[7] K. K. Tan, Q.-G. Wang, C. C. Hang, T. Hägglund, "Advances in PID control," Advances in Industrial Control Series,DOI: 10.1007/978-1-4471-0861-0, 1999.

[8] T. Liu, F. Gao, Industrial Process Identification and Control Design. Step-Test and Relay-Experiment-Based Methods,DOI: 10.1007/978-0-85729-977-2, 2012.

[9] G. P. Rangaiah, P. R. Krishnaswamy, "Estimating second-order dead time parameters from underdamped process transients," Chemical Engineering Science, vol. 51 no. 7, pp. 1149-1155, DOI: 10.1016/S0009-2509(96)80013-3, 1996.

[10] H.-P. Huang, M.-W. Lee, C.-L. Chen, "A system of procedures for identification of simple models using transient step response," Industrial & Engineering Chemistry Research, vol. 40 no. 8, pp. 1903-1915, DOI: 10.1021/ie0005001, 2001.

[11] V. M. Alfaro, "Low-order model identification from the process reaction curve” (in Spanish)," Ciencia y Tecnologia, vol. 24 no. 2, pp. 197-216, 2006.

[12] W. K. Ho, C. C. Hang, L. S. Cao, "Tuning of PID controllers based on gain and phase margin specifications," Automatica, vol. 31 no. 3, pp. 497-502, DOI: 10.1016/0005-1098(94)00130-B, 1995.

[13] C. L. Smith, Digital Computer Process Control, 1972.

[14] M. Viterková, A. Vitecek, L. Smutny, "Simple PI and PID controllers tuning for monotone self-regulation plants," IFAC Proceedings Volumes, vol. 33 no. 4, pp. 259-264, DOI: 10.1016/S1474-6670(17)38254-X, 2000.

[15] A. Jahanmiri, H. R. Fallahi, "New methods for process identification and design of feedback controller," Chemical Engineering Research and Design, vol. 75 no. 5, pp. 519-522, DOI: 10.1205/026387697524038, 1997.

[16] R. A. Mollenkamp, Introduction to Automatic Process Control, 1984.

[17] G. P. Rangaiah, P. R. Krishnaswamy, "Estimating second-order plus dead time model parameters," Industrial & Engineering Chemistry Research, vol. 33 no. 7, pp. 1867-1871, DOI: 10.1021/ie00031a029, 1994.

[18] V. M. Alfaro, R. Vilanova, "Control of high-order processes: repeated-pole plus dead-time models' identification," International Journal of Control, vol. 24,DOI: 10.1080/00207179.2021.1954240, 2021.

[19] C. Ionescu, A. Lopes, D. Copot, J. A. T. Machado, J. H. T. Bates, "The role of fractional calculus in modeling biological phenomena: a review," Communications in Nonlinear Science and Numerical Simulation, vol. 51, pp. 141-159, DOI: 10.1016/j.cnsns.2017.04.001, 2017.

[20] A. Narang, S. Shah, T. Chen, "Continuous-time model identification of fractional-order models with time delays," IET Control Theory & Applications, vol. 5 no. 7, pp. 900-912, DOI: 10.1049/iet-cta.2010.0718, 2011.

[21] J.-D. Gabano, T. Poinot, H. Kanoun, "Identification of a thermal system using continuous linear parameter-varying fractional modelling," IET Control Theory & Applications, vol. 5 no. 7, pp. 889-899, DOI: 10.1049/iet-cta.2010.0222, 2011.

[22] E. Guevara, H. Meneses, O. Arrieta, R. Vilanova, A. Visioli, F. Padula, "Fractional order model identification: computational optimization," Proceedings of the IEEE 20th Conference on Emerging Technologies & Factory Automation (ETFA),DOI: 10.1109/ETFA.2015.7301630, .

[23] H. Malek, Y. Luo, Y. Chen, "Identification and tuning fractional order proportional integral controllers for time delayed systems with a fractional pole," Mechatronics, vol. 23 no. 7, pp. 746-754, DOI: 10.1016/j.mechatronics.2013.02.005, 2013.

[24] B. B. Alagoz, A. Tepljakov, A. Ates, E. Petlenkov, C. Yeroglu, "Time-domain identification of one noninteger order plus time delay models from step response measurements," International Journal of Modeling, Simulation, and Scientific Computing, vol. 10 no. 01,DOI: 10.1142/S1793962319410113, 2019.

[25] S. Ahmed, "Parameter and delay estimation of fractional order models from step response," IFAC-PapersOnLine, vol. 48 no. 8, pp. 942-947, DOI: 10.1016/j.ifacol.2015.09.091, 2015.

[26] M. Tavakoli-Kakhki, M. Haeri, M. Saleh Tavazoei, "Simple fractional order model structures and their applications in control system design," European Journal of Control, vol. 16 no. 6, pp. 680-694, DOI: 10.3166/ejc.16.680-694, 2010.

[27] M. Tavakoli-Kakhki, M. Saleh Tavazoei, "Estimation of the order and parameters of a fractional order model from a noisy step response Data," Journal of Dynamic Systems, Measurement, and Control, vol. 136 no. 3,DOI: 10.1115/1.4026345, 2014.

[28] M. Tavakoli-Kakhki, M. Saleh Tavazoei, A. Mesbahi, "Parameter and order estimation from noisy step response data," IFAC Proceedings Volumes, vol. 46 no. 1, pp. 492-497, DOI: 10.3182/20130204-3-FR-4032.00033, 2013.

[29] A. Tepljakov, B. B. Alagoz, C. Yeroglu, E. A. Gonzalez, S. H. HosseinNia, E. Petlenkov, A. Ates, M. Cech, "Towards industrialization of FOPID controllers: a survey on milestones of fractional-order control and pathways for future developments," IEEE Access, vol. 9, pp. 21016-21042, DOI: 10.1109/ACCESS.2021.3055117, 2021.

[30] A. Tepljakov, B. B. Alagoz, C. Yeroglu, E. Gonzalez, S. H. HosseinNia, E. Petlenkov, "FOPID controllers and their industrial applications: a survey of recent results 1 1This study is based upon works from cost action CA15225, a network supported by cost (European cooperation in science and technology)," IFAC-PapersOnLine, vol. 51 no. 4, pp. 25-30, DOI: 10.1016/j.ifacol.2018.06.014, 2018.

[31] I. Birs, C. Muresan, I. Nascu, C. Ionescu, "A survey of recent advances in fractional order control for time delay systems," IEEE Access, vol. 7, pp. 30951-30965, DOI: 10.1109/ACCESS.2019.2902567, 2019.

[32] A. A. Dastjerdi, B. M. Vinagre, Y. Chen, S. H. HosseinNia, "Linear fractional order controllers; A survey in the frequency domain," Annual Reviews in Control, vol. 47, pp. 51-70, DOI: 10.1016/j.arcontrol.2019.03.008, 2019.

[33] P. Shah, S. Agashe, "Review of fractional PID controller," Mechatronics, vol. 38, pp. 29-41, DOI: 10.1016/j.mechatronics.2016.06.005, 2016.

[34] S. Momani, R. El-Khazali, I. M. Batiha, "Tuning PID and PI λ D δ controllers using particle swarm optimization algorithm via El-Khazali’s approach," Proceedings of the 45th International Conference on Application of Mathematics in Engineering and Economics (Amee’19), .

[35] I. M. Batiha, R. El-Khazali, S. Momani, "Dynamic responses of PI λ , PD δ , and PI λ d δ controllers for prosthetic hand model using PSO algorithm," Proceedings of the 2019 IEEE International Symposium on Signal Processing and Information Technology (ISSPIT),DOI: 10.1109/ISSPIT47144.2019.9001835, .

[36] M. L. Frikh, F. Soltani, N. Bensiali, N. Boutasseta, N. Fergani, "Fractional order PID controller design for wind turbine systems using analytical and computational tuning approaches," Computers & Electrical Engineering, vol. 95,DOI: 10.1016/j.compeleceng.2021.107410, 2021.

[37] R. El-Khazali, I. M. Batiha, S. Momani, "The drug administration via fractional-order PI λ d δ -controller," Progress in Fractional Differentiation and Applications, vol. 8 no. 1, pp. 53-62, DOI: 10.18576/pfda/080103, 2022.

[38] C. A. Monje, B. M. Vinagre, V. Feliu, Y. Chen, "Tuning and auto-tuning of fractional order controllers for industry applications," Control Engineering Practice, vol. 16 no. 7, pp. 798-812, DOI: 10.1016/j.conengprac.2007.08.006, 2008.

[39] Y. Luo, Y. Q. Chen, C. Y. Wang, Y. G. Pi, "Tuning fractional order proportional integral controllers for fractional order systems," Journal of Process Control, vol. 20 no. 7, pp. 823-831, DOI: 10.1016/j.jprocont.2010.04.011, 2010.

[40] H. Li, Y. Luo, Y. Chen, "A fractional order proportional and derivative (FOPD) motion controller: tuning rule and experiments," IEEE Transactions on Control Systems Technology, vol. 18 no. 2, pp. 516-520, DOI: 10.1109/TCST.2009.2019120, 2010.

[41] M. Tavakoli-Kakhki, M. Haeri, "Fractional order model reduction approach based on retention of the dominant dynamics: application in IMC based tuning of FOPI and FOPID controllers," ISA Transactions, vol. 50 no. 3, pp. 432-442, DOI: 10.1016/j.isatra.2011.02.002, 2011.

[42] J. J. Gude, E. Kahoraho, "Simple tuning rules for fractional PI controllers," Proceedings of the IEEE 14th Conference on Emerging Technologies & Factory Automation,DOI: 10.1109/ETFA.2009.5347157, .

[43] J. J. Gude, E. Kahoraho, "Modified ziegler-nichols method for fractional PI controllers," Proceedings of the IEEE 15th Conference on Emerging Technologies & Factory Automation,DOI: 10.1109/ETFA.2010.5641074, .

[44] I. Podlubny, Fractional Differential Equations, 1999.

[45] C. A. Monje, Y. Chen, B. M. Vinagre, D. Xue, V. Feliu, "Fractional-order systems and controls. fundamentals and applications," Advances in Industrial Control Series,DOI: 10.1007/978-1-84996-335-0, 2010.

[46] C. I. Muresan, C. M. Ionescu, "Generalization of the FOPDT model for identification and control purposes," Processes, vol. 8 no. 6, pp. 682-699, DOI: 10.3390/pr8060682, 2020.

[47] Y. Chen, I. Petras, D. Xue, "Fractional order control—a tutorial," Proceedings of the American Control Conference,DOI: 10.1109/ACC.2009.5160719, .

[48] D. Xue, Fractional-Order Control Systems: Fundamentals and Numerical Implementations,DOI: 10.1515/9783110497977-202, 2017.