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

#### 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 . The continuous integrodifferential operator is defined aswhere a and *t* are bounds of the operation and *α* ∈.

There are several different definitions of fractional operators [44]. One of the most used definitions of the fractional integration is the Riemann–Liouville definition:where *t* ≥ 0, *α* ∈ 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:where *n*−1<*α*<*n*, *n*∈. 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, is denoted by and by .

The Laplace transform of the Riemann–Liouville-based fractional derivative iswhere *n* − 1 ≤ *α* < *n*, which for zero initial conditions is reduced to

A general fractional-order system can be described by a fractional differential equation of the form

The corresponding fractional-order transfer function of incommensurate real orders has the following form [44]:where P() and Q() have no common zeros, _{k} (*k* = 0,…,*n*) and *b*_{k}(*k* = 0, …, *m*) are constants, and _{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} > · · · > *α*_{0}, and *β*_{m} > *ß*_{m} _{−} _{1} > · · · > *ß*_{0}.

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, , and incommensurate order system (7) can also be rewritten in the commensurate form as follows:

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:where _{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].where *E*_{α,α}(*z*) denotes the so-called two-parameter Mittag–Leffler function, which for an arbitrary value *z* is defined as

Integrating the right-hand side of (10), the following step response of the transfer function *G*(*s*^{α}) is obtained:

Each component of the step response (*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 aswhere initial conditions are generally taken as zero to obtain the standard FFOPDT transfer function model: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 parametersrepresents the FFOPDT model parameters, which will be identified in this paper using information taken from the process reaction curve.

The standard FOPDTcan 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.

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 aswhere (*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

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:where (*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.

##### 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 iswhere *E*_{α,β} is the two-parameter Mittag–Leffler function defined in (11).

From process reaction curve (20), the gain is given bywhere ∆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 , and equation (20) is reduced to the following expression:

If is the normalized process output, is the normalized time required to reach some specific output normalized value 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

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:

The so-called ratio index Δ in [18] can be adapted to this fractional framework accordingly and can be used to determine the fractional order *α*: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*_{1}(Δ).

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 arewhere

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.where *f*_{1} depends on ∆, which is defined as a function of the generic three points in (25), *f*_{2}(*α*) = a^{α} and *f*_{3}(*α*) = *τ*_{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

#### 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*_{50}) (*t*_{x2} = *t*_{50}). 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.

Adapting equations (29) with these symmetry considerations, the following expressions are obtained:where *f*_{1} is the function that relates fractional order *α* with the times ratio ∆.

In addition, Δ is defined as a function of the symmetrical three-points asand it depends on the normalized times *τ*_{x}, *τ*_{50}, and *τ*_{100‒x}. *f*_{2}(*α*) = a^{α}, whereand *f*_{3} (*α*) = (*τ*_{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}, *τ*_{50}, and *τ*_{100‒x}, for values of *x* = 5, 10, 15, 20, 25, and 30, respectively.

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*_{1}(Δ), *f*_{2}(*α*), and *f*_{3}(*α*) are estimated using curve fitting. Finally, fractional-order model parameters {*T*, *L*, *α*} can be determined from expressions (31) and experimental values *t*_{x}, *t*_{50}, 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}, *τ*_{50}, 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).

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

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}, *τ*_{50}, 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:where Table 7 shows the values of the parameters {*p*_{i}, *q*_{i}} for function *f*_{1}(Δ) and each one of the selected sets of symmetrical points.

It is important to note that for increasing values of *x*, Δ range becomes shorter. *f*_{1}(Δ) 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*_{2}(*α*) = a^{α} depends on *α* and normalized times *τ*_{x} and *τ*_{100‒x}.

In Figure 6, it can be seen that the amplitude of *f*_{2} increases for increasing values of *x*. On the other hand, for decreasing values of *x*, function *f*_{2} 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:where the values of the parameters {*p*_{i}, *q*_{i}} for function *f*_{2}(*α*) and each one of the selected sets of symmetrical points are shown in Table 8.

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*_{3}(*α*) = (*τ*_{100−x})^{1/α} depends on *α* and normalized time *τ*_{100‒x}. The range of amplitudes in *f*_{3} 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:where the values of the parameters {*p*_{i}, *q*_{i}} for function *f*_{3}(*α*) and each one of the selected sets of symmetrical points are shown in Table 9.

It is important to emphasize that estimated values of *α*, *T*, and *L* obtained from expressions (31) depend on functions *f*_{1}(Δ), *f*_{2}(*α*), and *f*_{3}(*α*). 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*_{2} and *f*_{3}—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*_{1}(Δ), *f*_{2}(*α*), and *f*_{3}(*α*) 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*_{1}(Δ), *f*_{2}(*α*), and *f*_{3}(*α*)), 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*_{1}, *f*_{2}, and *f*_{3} 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*_{1}, *f*_{2}, and *f*_{3} 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*_{50}), 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*, *α*}.

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

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*_{2} 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:where is the vector of process model parameters, is the open-loop step response difference between the actual process and the identified process model output signal, and , respectively, *T*_{S} is the sampling period, and *N*_{S}*T*_{S} is the time length of the dynamic (transient) response. Although in the context of a step response test, *N*_{S}*T*_{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:where and are the MSE values in (39) and and 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} = 1, *T*_{1} = 1 s, *L*_{1} = 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 for each process *P*_{1,i} are obtained using the proposed identification method for each of the sets of points. is expressed as follows: 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 are determined.

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 that have been obtained with the proposed procedure for each process (*i* = 1, …, 9) and using the different methods (*j* = 1, …, 6).

Note that the values 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 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 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), , 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.

The case *λ*_{5} = 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.

Considering data from Table 13, the FFOPDT process parameters , where *i* = 1, …, 6, have been obtained in Table 14 for each one of the proposed sets of symmetrical points in Table 12.

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.

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.

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.

Table 18 shows the values of the time-domain model performance index for the different identification methods applied to process *P*_{1,5}. The relative performance index values , which allow the comparison between different identification methods, are listed in Table 19.

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*_{1} 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} = 2, *T*_{2} = 1 s, *λ*_{2} = 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.

Considering data from Table 20, the FFOPDT process parameters , 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.

Process (38) is approximated by a FFOPDT model following three different strategies proposed in [26], where model process parameters are given in Table 22.

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*_{2} in comparison with the process reaction curve are shown in Figure 22.

Table 23 shows the values of the time-domain model performance index for the different identification methods applied to process *P*_{2}. The relative performance index values , which show the comparison between different identification methods, are listed in Table 24.

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.

#### Data Availability

The data sets used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare no conflicts of interest.

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