Proposal of a General Identification Method for Fractional-Order Processes Based on the Process Reaction Curve

Abstract

This paper aims to present a general identification procedure for fractional first-order plus dead-time (FFOPDT) models. This identification method is general for processes having S-shaped step responses, where process information is collected from an open-loop step-test experiment, and has been conducted by fitting three arbitrary points on the process reaction curve. In order to validate this procedure and check its effectiveness for the identification of fractional-order models from the process reaction curve, analytical expressions of the FFOPDT model parameters have been obtained for both situations: as a function of any three points and three points symmetrically located on the reaction curve, respectively. Some numerical examples are provided to show the simplicity and effectiveness of the proposed procedure. Good results have been obtained in comparison with other well-recognized identification methods, especially when simplicity is emphasized. This identification procedure has also been applied to a thermal-based experimental setup in order to test its applicability and to obtain insight into the practical issues related to its implementation in a microprocessor-based control hardware. Finally, some comments and reflections about practical issues relating to industrial practice are offered in this context.

1. Introduction

For the design and tuning tasks of a control system, information about the dynamic behavior of the controlled process is required, capturing it into a process mathematical model [1]. This model must provide reliable information to predict the effects that a control system will have on the behavior of the controlled process at the operating point, in particular on its output variable in servo and regulatory control.

The academic and industrial community recognizes that proportional integral derivative (PID) controller is the most widely used option in the process industry, having become an industry standard for process control [2]. Although the most commonly used models for tuning PID controllers are the first-, dual-pole, and second-order plus dead-time (FOPDT, DPPDT, SOPDT) ones, they are not able to represent the process dynamics with the required accuracy in some cases, as suggested, e.g., in [3]. Due to this, increasing the robustness of the control system can introduce issues in the well-known robustness/performance trade-off, constraining the desired performance [4]. Therefore, obtaining more accurate models of the controlled process is expected to improve performance.

In the technical literature there are a wide variety of identification procedures for integer-order models that are based on an open-loop step-test experiment; see, e.g., the identification methods detailed in [5,6,7,8] and the references cited therein. In general, these types of identification procedures are characterized by being performed with very little information about the plant, which makes them ideal for their use in industry. More specifically, some references also described identification algorithms that are based on fitting different points on the reaction curve [9,10,11]. Generally, identification methods for integer-order models that are based on the location of two or three points on the process reaction curve use FOPDT, DPPDT, and SOPDT models.

Considering both two-point and three-point methods separately, in the first group, the following methods can be considered [11,12,13,14], while in the second group, methods [15,16,17] can be examined. The identification method proposed in [11], which is used as a two-point method for FOPDT and DPPDT models, can also be used as a three-point method for SOPDT models. In a recent paper, the method proposed in [11] has been also extended for identifying a multiple-pole with dead-time model [18].

Note that in the case of two-point methods considered above, the sets of points are asymmetrical with respect to the central point on the reaction curve (50%). The only method with a symmetrical set of points is the 123c identification method proposed by Alfaro in [11]. For the considered three-point methods, refs. [15,17] present asymmetrical sets of points, while [11] for SOPDT models and ref. [16] provide sets of points that are symmetrical, although for the latter the central point does not match the central value on the reaction curve.

In recent decades, the advent of fractional calculus and new computational techniques have made possible a major academic and industrial effort focused on the transition from classical models and controllers to those described by non-integer order differential equations. Thus, fractional-order dynamic models and controllers were introduced; see, e.g., [19,20].

The apparent benefit of fractional calculus in the field of modelling has been justified from an industrial point of view. However, it has been more difficult to convey the advantages of fractional calculus on the controller side because of implementation issues [21]. That is why the adoption of fractional-order PID controllers in the industry is currently low, even though fractional-order PID controllers offer clear advantages in comparison with integer-order ones.

Recent studies have pointed to fractional-order PID controllers as an emerging trend in process control (see [21,22,23,24]), the main reason for their success being the intrinsic robustness they offer with a higher degree of freedom to operate and tune the controller parameters. It is important to note that in some of the existing fractional-order controller design methods, a simple process model has been used in order to tune the intended controller settings [25,26,27,28,29,30].

In the technical literature there is a wide range of methods for identifying fractional-order models based on the process reaction curve; however, there are not many that are analytical techniques and whose main feature is simplicity of implementation. Some strategies to estimate the FFOPDT model parameters by using step-response data have been proposed in [31]. These combine numerical computation and graphical estimation. Integral-based estimation methods, whose main feature is their robustness to the presence of measurement noise, are proposed in [32,33].

The most common approach in industrial practice is based on nonlinear optimization; see, e.g., [34,35,36,37]. These methods are generally applied by minimizing the error between the fractional-order model step response and the process reaction curve. These techniques are characterized by the fact that they require more computational effort compared to other existing analytical methods.

Despite the fact that the fractional-order model has been demonstrated to be technologically superior on multiple occasions, the industrial adoption of the fractional-order approach requires further analysis [38].

Considering all the above, the existence of identification methods for simple-structure fractional-order models is of major relevance and can be very helpful in the practical design of integer- and fractional-order control systems. There are multiple reaction curve-based methods for identifying integer-order models, as discussed previously. This is mainly due to the fact that the step response of a process has a straightforward physical interpretation and that identification methods for integer-order models based on fitting several points of the process reaction curve are very easy to implement and apply. Therefore, one may consider it appropriate to extend such methods for fractional models. To that end, an FFOPDT model identification method, which is based on fitting three points on the reaction curve, has been conducted in [39]. However, that study was restricted to considering that the central point is located in the center of the reaction curve (50%), and that the extreme points are located symmetrically with respect to this central point.

It is the authors’ opinion that it is of significant interest to extend the identification method proposed in [39] to any set of points on the process reaction curve.

Therefore, the objective of this work is, on the one hand, to validate this identification procedure for three asymmetrical points on the reaction curve and to obtain insight into the selection of such points and their influence on the accuracy of the identified fractional-order model. On the other hand, another additional objective is to test its applicability to a laboratory prototype and to obtain some insights on the practical issues related to its implementation on a microprocessor-based control hardware.

This paper is organized as follows. Section 2 is devoted to presenting some preliminaries and theoretical background. In Section 3, a generalization of the three-point step-response method for FFOPDT models is proposed. This general method is particularized for several symmetrical and asymmetrical sets of points on the process reaction curve. The results of some numerical simulations are presented in Section 4, illustrating the effectiveness and simplicity of the proposed method for both symmetrical and asymmetrical sets of points in comparison to well-known identification methods. In Section 5, the applicability of the proposed identification procedure on a laboratory prototype is verified, showing several practical issues related to its implementation on an industrial control hardware. Finally, Section 6 presents conclusions of this paper.

2. Preliminaries and Theoretical Background

Some elementary definitions and basic concepts in fractional calculus are provided in this section. Elementary ideas from fractional calculus can be found in many books, such as [40,41].

The fractional integral of order α for function f(t) is defined as:

$${}_0{}^{ }\mathrm{I}{}_{\mathrm{t}}^{\mathsf{\alpha }}\mathrm{f}(\mathrm{t})\equiv {}_0{}^{ }\mathrm{D}{}_{\mathrm{t}}^{-\mathsf{\alpha }}\mathrm{f}(\mathrm{t})=\frac{1}{\mathsf{\Gamma }(\mathsf{\alpha })}{{\displaystyle \int }}_0^{\mathrm{t}}{(\mathrm{t}-\mathsf{\tau })}^{\mathsf{\alpha }-1}\mathrm{f}(\mathsf{\tau })\mathrm{d}\mathsf{\tau }$$

(1)

where Γ(·) is the Gamma function [40], t ≥ 0, and α $\in {ℝ}^{+}$.

The α-th order Riemann-Liouville definition of fractional derivative of the given function f(t) is defined as:

$${}_0{}^{ }\mathrm{D}{}_{\mathrm{t}}^{\mathsf{\alpha }}\mathrm{f}(\mathrm{t})=\frac{1}{\mathsf{\Gamma }(\mathrm{m}-\mathsf{\alpha })}\frac{{\mathrm{d}}^{\mathrm{m}}}{{\mathrm{dt}}^{\mathrm{m}}}{{\displaystyle \int }}_0^{\mathrm{t}}{(\mathrm{t}-\mathsf{\tau })}^{\mathrm{m}-\mathsf{\alpha }-1}\mathrm{f}(\mathsf{\tau })\mathrm{d}\mathsf{\tau }$$

(2)

where m − 1 < α < m, m $\in {\mathbb{Z}}^{+}$. The subscripts 0 and t in Definitions (1) and (2) can be considered as the limits of operation and referred to as the terminals of fractional-order integration and differentiation, respectively. For simplicity, ${}_0{}^{ }\mathrm{D}{}_{\mathrm{t}}^{-\mathsf{\alpha }}$ will be denoted by ${\mathrm{D}}^{-\mathsf{\alpha }}$ and ${}_0{}^{ }\mathrm{D}{}_{\mathrm{t}}^{\mathsf{\alpha }}$ by ${\mathrm{D}}^{\mathsf{\alpha }}$.

The Laplace transform of the fractional derivative based on Riemann-Liouville is:

$$\mathrm{L}\{{\mathrm{D}}^{\mathsf{\alpha }}\mathrm{f}(\mathrm{t})\}={\mathrm{s}}^{\mathsf{\alpha }}\mathrm{L}\{\mathrm{f}(\mathrm{t})\}-{\displaystyle \sum }_{\mathrm{k}=0}^{\mathrm{m}-1}{\mathrm{s}}^{\mathrm{k}}{\mathrm{D}}^{\mathsf{\alpha }-\mathrm{k}-1}\mathrm{f}(0)$$

(3)

where m − 1 < α < m. For zero initial conditions, (3) is reduced to:

$$\mathrm{L}\{{\mathrm{D}}^{\mathsf{\alpha }}\mathrm{f}(\mathrm{t})\}={\mathrm{s}}^{\mathsf{\alpha }}\mathrm{L}\{\mathrm{f}(\mathrm{t})\}$$

(4)

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

$${\mathrm{a}}_{\mathrm{n}}{\mathrm{D}}^{{\mathsf{\alpha }}_{\mathrm{n}}}\mathrm{y}(\mathrm{t})+{\mathrm{a}}_{\mathrm{n}-1}{\mathrm{D}}^{{\mathsf{\alpha }}_{\mathrm{n}-1}}\mathrm{y}(\mathrm{t})+\cdots {\mathrm{a}}_0{\mathrm{D}}^{{\mathsf{\alpha }}_0}\mathrm{y}(\mathrm{t})={\mathrm{b}}_{\mathrm{m}}{\mathrm{D}}^{{\mathsf{\beta }}_{\mathrm{m}}}\mathrm{u}(\mathrm{t})+{\mathrm{b}}_{\mathrm{m}-1}{\mathrm{D}}^{{\mathsf{\beta }}_{\mathrm{m}-1}}\mathrm{u}(\mathrm{t})+\cdots +{\mathrm{b}}_0{\mathrm{D}}^{{\mathsf{\beta }}_0}\mathrm{u}(\mathrm{t})$$

(5)

The transfer function of incommensurate real orders corresponding to the differential Equation (5) has the following expression [40]:

$$\mathrm{G}(\mathrm{s})=\frac{\mathrm{Q}({\mathrm{s}}^{{\mathsf{\beta }}_{\mathrm{k}}})}{\mathrm{P}({\mathrm{s}}^{{\mathsf{\alpha }}_{\mathrm{k}}})}=\frac{{\mathrm{b}}_{\mathrm{m}}{\mathrm{s}}^{{\mathsf{\beta }}_{\mathrm{m}}}+{\mathrm{b}}_{\mathrm{m}-1}{\mathrm{s}}^{{\mathsf{\beta }}_{\mathrm{m}-1}}+\cdots +{\mathrm{b}}_0{\mathrm{s}}^{{\mathsf{\beta }}_0}}{{\mathrm{s}}^{{\mathsf{\alpha }}_{\mathrm{n}}}+{\mathrm{a}}_{\mathrm{n}-1}{\mathrm{s}}^{{\mathsf{\alpha }}_{\mathrm{n}-1}}+\cdots +{\mathrm{a}}_1{\mathrm{s}}^{{\mathsf{\alpha }}_1}+{\mathrm{a}}_0}$$

(6)

where $\mathrm{P}({\mathrm{s}}^{{\mathsf{\alpha }}_{\mathrm{k}}})$ and $\mathrm{Q}({\mathrm{s}}^{{\mathsf{\beta }}_{\mathrm{k}}})$ have no common zeros, a_k (k = 0, …, n), b_k (k = 0, …, m) are constants, and α_k (k = 0, …, n), β_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 > ⋯ > β₀.

Considering that the transfer function G(s) given by (6) is strictly proper, then it is BIBO stable if and only if P(s) has no root in {Re(s) ≥ 0} [19].

A particular case occurs when a real number α exists as the greatest common divisor of α_i, i = 1, …, n and β_i, i = 0, …, m. That value is referred to as the commensurate order. It holds that α_k = kα, β_k = kα, 0 < α < 1, ∀k $\in \mathbb{Z}$, and the incommensurate order system (6) can also be rewritten in commensurate form as follows:

$$\mathrm{G}(\mathrm{s})=\frac{\mathrm{Q}({\mathrm{s}}^{\mathsf{\alpha }})}{\mathrm{P}({\mathrm{s}}^{\mathsf{\alpha }})}=\frac{{{\displaystyle \sum }}_{\mathrm{k}=0}^{\mathrm{m}}{\mathrm{b}}_{\mathrm{k}}{\mathrm{s}}^{\mathrm{k}\mathsf{\alpha }}}{{{\displaystyle \sum }}_{\mathrm{k}=0}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{k}}{\mathrm{s}}^{\mathrm{k}\mathsf{\alpha }}}$$

(7)

It has been proven that the commensurate system G(s) brought in (7) 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.

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 (7) can be written in the following general form:

$$\mathrm{G}({\mathrm{s}}^{\mathsf{\alpha }})={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{r}}_{\mathrm{i}}}{{\mathrm{s}}^{\mathsf{\alpha }}+{\mathsf{\lambda }}_{\mathrm{i}}}$$

(8)

where r_i, i = 1, …, n are the corresponding residues and λ_i, i = 1, …, n are the roots of P(x) = 0.

Taking the inverse Laplace transform from (8), the impulse response of G(s^α) is obtained, which is also given in [40].

$$\mathrm{h}(\mathrm{t})={\mathrm{L}}^{-1}\left\{{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\frac{{\mathrm{r}}_{\mathrm{i}}}{{\mathrm{s}}^{\mathsf{\alpha }}+{\mathsf{\lambda }}_{\mathrm{i}}}\right\}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{r}}_{\mathrm{i}}{\mathrm{t}}^{\mathsf{\alpha }-1}{\mathrm{E}}_{\mathsf{\alpha },\mathsf{\alpha }}(-{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{t}}^{\mathsf{\alpha }})$$

(9)

where E_α,α(z) is the Mittag–Leffler function. This function is defined for an arbitrary value z as:

$${\mathrm{E}}_{\mathsf{\alpha },\mathsf{\beta }}(\mathrm{z})={\displaystyle \sum }_{\mathrm{r}=0}^{\infty }\frac{{\mathrm{z}}^{\mathrm{r}}}{\mathsf{\Gamma }(\mathsf{\alpha }\mathrm{r}+\mathsf{\beta })}$$

(10)

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

$$\mathrm{g}(\mathrm{t})={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{r}}_{\mathrm{i}}\frac{{\mathrm{E}}_{\mathsf{\alpha },1}(-{\mathsf{\lambda }}_{\mathrm{i}}{\mathrm{t}}^{\mathsf{\alpha }})-1}{{\mathsf{\lambda }}_{\mathrm{i}}}$$

(11)

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

3. Fractional First-Order Plus Dead-Time Model Identification

The general identification procedure to be presented in this section uses process information obtained from an open-loop test and is applied only for identifying a fractional-order model for processes having an S-shaped step response.

The differential equation for an FFOPDT model can be expressed as follows:

$$\mathrm{T}·{\mathrm{D}}^{\mathsf{\alpha }}\mathrm{y}(\mathrm{t})+\mathrm{y}(\mathrm{t})=\mathrm{K}·\mathrm{u}(\mathrm{t}-\mathrm{L})$$

(12)

where the initial condition y(0⁺) is generally taken as zero to obtain a transfer function model. The controlled processes under consideration can be well characterized by this model. The standard FFOPDT model can be derived from Equation (12) by taking the Laplace transform, obtaining the following transfer function:

$$\mathrm{P}(\mathrm{s})=\frac{{\mathrm{Ke}}^{-\mathrm{Ls}}}{1+{\mathrm{Ts}}^{\mathsf{\alpha }}}$$

(13)

where K is the process gain, T > 0 is the time constant, L ≥ 0 is the apparent dead-time, and α is the fractional order of the model.

FFOPDT model parameters θ_P = {K, T, L, α} can be identified using the process reaction curve, as a result of an open-loop test.

For the particular case α = 1, the FFOPDT model (13) becomes the standard FOPDT model, which has been broadly used in practice to capture the essential dynamic response of industrial processes for the purpose of control design [2].

The FFOPDT model (13) can be viewed as a generalization of the conventional FOPDT model, as discussed in [42], where the relevance of this generalization and the implications it has for both the identification of dynamic processes and the design of feedback control loops are indicated.

The step response of the considered fractional model can conveniently describe both monotonic or non-monotonic behaviors depending on the fractional order α. The step responses of FFOPDT models for increasing values of α, with 0.5 ≤ α ≤ 1.1 which is the range considered in this paper, are shown in Figure 1. As a particular case, the step response of the considered system for α = 1, which represents an FOPDT model, is depicted in red line.

In this paper, two parameters will be used to characterize process dynamics, namely T_ar and τ, which are two classical parameters well-defined for integer-order processes and can also be defined for the fractional case.

T_ar is referred to as the average residence time and has been defined in the fractional context as follows [39]:

$${\mathrm{T}}_{\mathrm{ar}}\dot{=}\frac{{{\displaystyle \int }}_0^{\infty }\mathrm{tg}(\mathrm{t})\mathrm{dt}}{{{\displaystyle \int }}_0^{\infty }\mathrm{g}(\mathrm{t})\mathrm{dt}}=\mathrm{L}+{\mathrm{T}}^{1/\mathsf{\alpha }}$$

(14)

where g(t) represents the impulse response of the fractional-order process. The average residence time is a classical index that characterizes process dynamics by indicating the time it takes the input to have a significant influence on the output [2].

τ is referred to as the normalized dead-time, the typical range of which is 0 ≤ τ ≤ 1, and has been defined in the fractional context as follows [39]:

$$\mathsf{\tau }\dot{=}\frac{\mathrm{L}}{{\mathrm{T}}_{\mathrm{ar}}}=\frac{\mathrm{L}}{\mathrm{L}+{\mathrm{T}}^{1/\mathsf{\alpha }}}$$

(15)

This parameter has been defined to characterize the degree of difficulty in controlling a process. Broadly speaking, processes with small τ can be considered to be easy to control, while processes with a larger value of τ are difficult to control. The standard definition of the classical parameters T_ar and τ for the case α = 1 corresponds to the one for an FOPDT model [2].

Although with an FFOPDT model the monotonic and non-monotonic behavior of the process can be characterized, in this work we will only focus on processes with monotonic response with 0.50 ≤ α ≤ 1.00, characterized by an S-shaped response. The importance of processes with essentially monotone step responses lies in the fact that they are very common in process control [2].

The FFOPDT model normalized step responses for different values of α and different values of τ are shown in Figure 2. It can be noticed from the figure that all responses have a common point at t = T_ar because time is normalized with respect to the average residence time.

Considering that a step signal u(t) with an amplitude Δu is applied to an FFOPDT model (13), a signal y_α(t) with an amplitude variation of Δy is obtained, as shown in Figure 3.

The corresponding FFOPDT model response to a ∆u step input change is:

$${\mathrm{y}}_{\mathsf{\alpha }}(\mathrm{t})=\{\begin{array}{cc}0,& 0 \le \mathrm{t}<\mathrm{L}\\ \mathrm{K}\left\{1-{\mathrm{E}}_{\mathsf{\alpha },1}\left[-\frac{1}{\mathrm{T}}{(\mathrm{t}-\mathrm{L})}^{\mathsf{\alpha }}\right]\right\}\mathsf{\Delta }\mathrm{u},& \mathrm{t}\ge \mathrm{L}\end{array}$$

(16)

where E_α,β is the two-parameter Mittag–Leffler function, which was previously defined in (10).

Normalizing the process output y_α(t) with respect to its final value ∆y = K·∆u and using the shifted and normalized time $\mathsf{\tau }=\frac{1}{\mathrm{T}}{(\mathrm{t}-\mathrm{L})}^{\mathsf{\alpha }}$, Equation (16) is reduced to:

$$\begin{array}{cc}{\tilde{\mathrm{y}}}_{\mathsf{\alpha }}(\mathsf{\tau })=1-{\mathrm{E}}_{\mathsf{\alpha },1}(-\mathsf{\tau }),& \mathsf{\tau }\ge 0\end{array}$$

(17)

The FFOPDT model parameters, θ_P = {K, T, L, α}, will be determined considering the step response of an FFOPDT model (16) and the normalized process output (17), respectively.

From Equation (16), the gain is given by:

$$\mathrm{K}=\frac{\mathsf{\Delta }\mathrm{y}}{\mathsf{\Delta }\mathrm{u}}$$

(18)

where ∆y is the total process output change when a step input with amplitude ∆u is applied, as indicated in Figure 3.

Since ${\tilde{\mathrm{y}}}_{\mathsf{\alpha }}(\mathsf{\tau })$ is a specific value of the normalized output, which has the following range 0 ≤ ${\tilde{\mathrm{y}}}_{\mathsf{\alpha }}(\mathsf{\tau })$ ≤ 1, its corresponding normalized time τ_x can be found using (17). Therefore, the time t_x required for the process output (16) to reach such x-point is:

$${\mathrm{t}}_{\mathrm{x}}=\mathrm{L}+{({\mathsf{\tau }}_{\mathrm{x}}\mathrm{T})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}$$

(19)

In order to obtain the rest of the FFOPDT model parameters (T, L, and α), the set of times {t_x1, t_x2, t_x3} to reach points {y_α(t_x1), y_α(t_x2), y_α(t_x3)}, respectively, on the process reaction curve are required.

Note that, without loss of generality, t_x1 < t_x2 < t_x3 and τ_x1 < τ_x2 < τ_x3 will be considered in this paper.

Considering expressions like (19), the following equations set is defined:

$$\begin{array}{c}{\mathrm{t}}_{\mathrm{x}1}=\mathrm{L}+{({\mathsf{\tau }}_{\mathrm{x}1}\mathrm{T})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}\\ {\mathrm{t}}_{\mathrm{x}2}=\mathrm{L}+{({\mathsf{\tau }}_{\mathrm{x}2}\mathrm{T})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}\\ {\mathrm{t}}_{\mathrm{x}3}=\mathrm{L}+{({\mathsf{\tau }}_{\mathrm{x}3}\mathrm{T})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}\end{array}$$

(20)

In order for the model fractional order α to be estimated, the following ratio index Δ is found by considering the previously obtained set of Equation (20):

$$\mathsf{\Delta } \dot{=}\frac{{\mathrm{t}}_{\mathrm{x}3}-{\mathrm{t}}_{\mathrm{x}1}}{{\mathrm{t}}_{\mathrm{x}2}-{\mathrm{t}}_{\mathrm{x}1}}=\frac{{\mathsf{\tau }}_{\mathrm{x}3}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}-{\mathsf{\tau }}_{\mathrm{x}1}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}}{{\mathsf{\tau }}_{\mathrm{x}2}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}-{\mathsf{\tau }}_{\mathrm{x}1}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}}$$

(21)

where τ_x1, τ_x2, and τ_x3 are normalized times defined by (19) and which can be determined using Equation (17). Note that Δ index can be interpreted as the ratio between the time difference in reaching from x₁ to x₃% and from x₁ to x₂% of the total variation of the process output, as indicated in (21). From Equation (21) it is clear that there is a dependence between fractional order α and ratio index ∆, which can be expressed as α = f₁(Δ).

For obtaining the time-based parameters (T, L) from (20), it is required to consider the following two points, {t_x1, y_α(t_x1)} and {t_x3, y_α(t_x3)}, on the process reaction curve. Then, the expressions for these parameters are:

$$\mathrm{T}={\mathrm{a}}^{\mathsf{\alpha }}{({\mathrm{t}}_{\mathrm{x}3}-{\mathrm{t}}_{\mathrm{x}1})}^{\mathsf{\alpha }}$$

(22)

and

$$\mathrm{L}={\mathrm{t}}_{\mathrm{x}3}-{\mathsf{\tau }}_{\mathrm{x}3}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}{\mathrm{T}}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}$$

(23)

where

$$\mathrm{a}=\frac{1}{{\mathsf{\tau }}_{\mathrm{x}3}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}-{\mathsf{\tau }}_{\mathrm{x}1}{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathsf{\alpha }$}\right.}}$$

(24)

The two equivalent normalized points for obtaining parameters T and L are {ỹ_α(τ_x1), τ_x1} and {ỹ_α(τ_x3), τ_x3}.

The set of Equation (25) are the expressions required to obtain the FFOPDT model parameters, θ_P = {K, T, L, α}, by using the times needed for the response to reach any three points on the process reaction curve {x₁-x₂-x₃%}.

$$\{\begin{array}{c}\begin{array}{c}\mathrm{K}=\frac{\mathsf{\Delta }\mathrm{y}}{\mathsf{\Delta }\mathrm{u}}\\ \mathsf{\alpha }={\mathrm{f}}_1(\mathsf{\Delta })\end{array}\\ \mathrm{T}={\mathrm{f}}_2(\mathsf{\alpha }){({\mathrm{t}}_{\mathrm{x}3}-{\mathrm{t}}_{\mathrm{x}1})}^{\mathsf{\alpha }}\\ \mathrm{L}=\max \left[{\mathrm{t}}_{\mathrm{x}3}-{\mathrm{f}}_3(\mathsf{\alpha }){\mathrm{T}}^{1/\mathsf{\alpha }}, 0\right]\end{array}$$

(25)

where ∆y is the total process output change to the step input ∆u, as shown in Figure 3; function f₁ depends on ∆ and is defined in (21) as a function of the times {t_x1, t_x2, t_x3}; functions f₂(α) = a^α and f₃(α) = τ_x3^1/α depend on τ_x1 and τ_x3, and τ_x3, respectively, and finally a-parameter is defined in (24). It can be observed from (25) that the expressions of T and L have a high dependence on the value of α. This makes it of significant importance to determine the value of α parameter accurately.

In (25), α > 0 and T > 0 are fulfilled in a natural way, since τ_x1 < τ_x2 < τ_x3 and t_x1 < t_x2 < t_x3. The following condition must be satisfied to ensure L ≥ 0:

t_x3 ≥ τ_x3^1/α

(26)

For a more detailed development of Equation (25) we refer the reader to [39], where these equations are obtained and subsequently particularized for the case in which the three points on the process reaction curve are symmetrical with respect to the central point.

Figure 4 shows the general scheme of the complete procedure for obtaining the expressions for the identification of the FFOPDT model parameters, θ_P = {K, T, L, α}, from three arbitrary points of the process reaction curve.

This procedure is summarized in the following steps, as depicted in Figure 4.

• From the normalized process output (17), ${\tilde{\mathrm{y}}}_{\mathsf{\alpha }}(\mathsf{\tau })$, the values of the normalized times {τ_x1, τ_x2, τ_x3} of the three considered points on the process reaction curve are obtained for the different values of α, 0.50 ≤ α ≤ 1.10.
• Data sets {Δ, α}, {α, aα}, and {α, (τx3)1/α} are obtained for the considered set of points (x₁-x₂-x₃%) by using the values of the corresponding normalized times.

Note that this procedure is general and admits the points x₁, x₂, and x₃ to be arbitrary. In this paper, both a symmetrical and an asymmetrical location of the points on the process reaction curve will be considered.

3.

By means of a curve-fitting procedure, the values of the parameters {pi, qi} for the rational functions α = f1(Δ), f2(α), and f3(α), respectively, are obtained.

4.

From the rational functions α = f1(Δ), f2(α), and f3(α) obtained in the previous step, the expressions for the FFOPDT model parameters (25) are completed.

5.

Once the numerical values of α, f2, and f3 are determined, the values of the FFOPDT model parameters, θ_P = {K, T, L, α}, are calculated using expressions (25) and experimental data collected from the process reaction curve, {Δy, Δu, t_x1, t_x2, t_x3}.

A detailed algorithm for estimating the FFOPDT model parameters, θ_P = {K, T, L, α}, from the information collected from the process reaction curve is provided in Section 5.4. This algorithm, which simplifies its software implementation, can serve as a complement to the scheme represented in Figure 4.

Table 1. Numeric values of normalized times τ₅, τ₁₀, τ₂₀, τ₂₅, τ₅₀, τ₅₅, τ₆₀, τ₇₅, τ₉₀, and τ₉₅, respectively, for different values of α, with 0.5 ≤ α ≤ 1.1.

α	τ5 [s]	τ10 [s]	τ20 [s]	τ25 [s]	τ50 [s]	τ55 [s]	τ60 [s]	τ75 [s]	τ90 [s]	τ95 [s]
0.5	0.0462	0.0963	0.2113	0.2781	0.7691	0.9216	1.1072	2.0516	5.5556	11.3640
0.6	0.0464	0.0965	0.2101	0.2752	0.7402	0.8810	1.0481	1.8734	4.7612	9.3440
0.7	0.0471	0.0975	0.2107	0.2749	0.7181	0.8470	0.9988	1.7127	3.9916	7.4160
0.8	0.0481	0.0993	0.2131	0.2768	0.7028	0.8220	0.9601	1.5757	3.2873	5.5890
0.9	0.0495	0.1019	0.2172	0.2811	0.6945	0.8060	0.9326	1.4665	2.7112	4.0190
1.0	0.0513	0.1054	0.2232	0.2877	0.6932	0.7990	0.9163	1.3863	2.3026	2.9960
1.1	0.0536	0.1097	0.2309	0.2967	0.6988	0.8000	0.9109	1.3334	2.0419	2.4560

includes the numerical values of the normalized times τ₅, τ₁₀, τ₂₀, τ₂₅, τ₅₀, τ₅₅, τ₆₀, τ₇₅, τ₉₀, and τ₉₅, respectively. These data will be used indistinctly as τ_x, τ₅₀, or τ_100−x in Section 3.1 for a symmetrical set of points (x-50-(100 − x)%), and τ_x1, τ_x2, or τ_x3 in Section 3.2 for an asymmetrical set of points (x₁-x₂-x₃%), respectively, and constitute the main source of data to determine the corresponding data sets {Δ, α}, {α, a^α}, and {α, (τ_x3)^1/α} or {α, (τ_100−x)^1/α}, needed to determine expressions for the functions α = f₁(Δ), f₂(α), and f₃(α).

3.1. Symmetrical Set of Points (x-50-(100 − x)%)

This case can be considered as a simplification of the general identification procedure, where only points that are symmetrically located on the process reaction curve are selected. Note that the central point or centroid will be located in the middle of the range, x₂ = 50%, (t_x2 = t₅₀, y_α(t₅₀)), as depicted in Figure 3. The remaining two points could be located arbitrarily on the reaction curve, but symmetrically located with respect to the central point (Δx₃₂ = Δx₂₁). One of the extreme points will be denoted x₁ = x, the other being x₃ = 100 − x. In this case, the times to be determined will be t_x1 = t_x and t_x3 = t_100−x, where t_x and t_100−x denote the time required to reach x% (y_α(t_x)) and (100 − x)% (y_α(t_100−x)) of the total process output change, respectively, considering the following range 0 < x < 50.

In this section, the sets of points indicated in

Table 2. Sets of symmetrical points that have been considered.

Set #	Symmetrical Points	Centroid	Distance from Centroid
1	(5-50-95%)	x2 = 50%	Δx21 = Δx32 = 45%
2	(10-50-90%)	x2 = 50%	Δx21 = Δx32 = 40%
3	(25-50-75%)	x2 = 50%	Δx21 = Δx32 = 25%

will be considered. For the case of symmetrical points, sets #1 and #2 have been chosen because the accuracy of the fractional-order model is improved for low values of x. In particular, the influence of the location of the symmetrical representative points of the reaction curve on the accuracy of the identified model is explained in detail in [39]. In addition, set #3 has been chosen because it gives very good results for integer- or close to integer-order models, see also [39]. This behavior has already been observed in the technical literature in some identification methods for integer-order models, e.g., the aforementioned 123c method [11] uses time sets (25–75%) to identify FOPDT and DPPDT models and (25-50-75%) for SOPDT models; and in [18], points (25-50-75%) are proposed to identify multiple-pole with dead-time models.

Then, the procedure summarized in Figure 4 is followed for the symmetrical case. Data sets {∆, α}, {α, a^α}, and {α, (τ_100−x)^1/α} for 0.5 ≤ α ≤ 1.1, and the functions f₁(Δ), f₂(α), and f₃(α) obtained by curve fitting for the different sets of symmetrical points, respectively, are shown in Figure 5. Note that Δ depends on normalized times {τ_x, τ₅₀, and τ_100−x} and has a significant dependence on α parameter, that f₂(α) = a^α depends on α and normalized times τ_x and τ_100−x, and that f₃(α) = (τ_100−x)^1/α depends on α and normalized time τ_100−x.

The following rational functions have been used for curve fitting of functions f₁(Δ), f₂(α), and f₃(α), respectively:

$${\mathrm{f}}_1(\mathsf{\Delta })=\frac{{\mathrm{p}}_1{\mathsf{\Delta }}^2+{\mathrm{p}}_2\mathsf{\Delta }+{\mathrm{p}}_3}{{\mathsf{\Delta }}^2+{\mathrm{q}}_1\mathsf{\Delta }+{\mathrm{q}}_2}$$

(27)

$${\mathrm{f}}_2(\mathsf{\alpha })=\frac{{\mathrm{p}}_1\mathsf{\alpha }+{\mathrm{p}}_2}{{\mathsf{\alpha }}^2+{\mathrm{q}}_1\mathsf{\alpha }+{\mathrm{q}}_2}$$

(28)

$${\mathrm{f}}_3(\mathsf{\alpha })=\frac{{\mathrm{p}}_1{\mathsf{\alpha }}^2+{\mathrm{p}}_2\mathsf{\alpha }+{\mathrm{p}}_3}{{\mathsf{\alpha }}^2+{\mathrm{q}}_1\mathsf{\alpha }+{\mathrm{q}}_2}$$

(29)

The Levenberg–Marquardt least-squares curve-fitting algorithm has been used for fitting data in all graphs in Figure 5. The values of the corresponding parameters {p_i, q_i} for functions f₁(Δ), f₂(α), and f₃(α), and each one of the selected sets of symmetrical points are shown in

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

(5-50-95%)	(10-50-90%)	(25-50-75%)
p1 = 0.4259	p1 = 0.3808	p1 = 0.2676
p2 = 38.78	p2 = 13.57	p2 = 1.756
p3 = 14.34	p3 = −3.067	p3 = −2.578
q1 = 45.33	q1 = 14.69	q1 = −0.7042
q2 = −27.8	q2 = −15.9	q2 = −1.289

,

Table 4. Parameters {p_i, q_i} of the rational function f2(α) for the different sets of symmetrical points.

(5-50-95%)	(10-50-90%)	(25-50-75%)
p1 = −0.0337	p1 = −0.05698	p1 = −0.3443
p2 = 0.0595	p2 = 0.1596	p2 = 0.7806
q1 = −2.328	q1 = −2.528	q1 = −3.017
q2 = 1.404	q2 = 1.753	q2 = 2.496

and

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

(5-50-95%)	(10-50-90%)	(25-50-75%)
p1 = 10.43	p1 = 4.518	p1 = 3.901
p2 = −22.09	p2 = −8.675	p2 = −4.833
p3 = 13.14	p3 = 5.666	p3 = 4.546
q1 = −0.586	q1 = −0.3471	q1 = 2.239
q2 = 0.07943	q2 = 0.003197	q2 = −0.6319

, respectively.

3.2. Asymmetrical Set of Points (x₁-x₂-x₃%)

In this section, the general identification procedure is applied to three arbitrary points asymmetrically located on the reaction curve (x₁, x₂, and x₃), as shown in Figure 3. This means that the times to be determined will be t_x1, t_x2, and t_x3, where they denote the time needed to reach x₁% (y_α(t_x1)), x₂% (y_α(t_x2)), and x₃% (y_α(t_x3)) of the total process output change, respectively.

The origin of the asymmetry comes not only from the fact that, in general, the centroid x₂ will not be located in the middle of the range (0–100%) of the process output, but also from the fact that the distance from the centroid x₂ to the extreme points x₁ and x₃, respectively, is different, i.e., Δx₂₁ ≠ Δx₃₂.

As discussed previously, the objective is to validate this identification procedure for the asymmetrical case and to obtain insight into the selection of representative points of the reaction curve and their influence on the accuracy of the identified model.

In this section, the sets of points indicated in

Table 6. Sets of asymmetrical points that have been considered.

Set #	Asymmetrical Points	Centroid	Distance from Centroid
4	(10-55-90%)	x2 = 55%	Δx21 = 45%, Δx32 = 35%
5	(20-60-95%)	x2 = 60%	Δx21 = 40%, Δx32 = 35%
6	(20-75-95%)	x2 = 75%	Δx21 = 55%, Δx32 = 25%

will be considered.

The selection of these three sets of points (x₁-x₂-x₃%) is based on experimentation. A large number of experiments have been performed in order to draw the following observations:

• The three sets of points have been chosen with a high x₃ value, where x₃ = 90 or 95%, because the obtained model fits better the reaction curve, especially in the final part. In this regard, it has been shown in [39] that the step response of the identified models gives a good fit with the process reaction curve for the symmetrical case, particularly in the interval [x-(100 − x)]. Due to the symmetry exhibited by this method, the interval [x-(100 − x)] is larger for lower values of x and, therefore, the step response of these models fits better the process reaction curve, which translates into a lower value in the performance index S for this fractional-order model.
• In general, the selection of x₁ affects the accuracy of the model in the initial part and, together with x₃, allows better fitting of T parameter.
• With respect to the centroid x₂, set #4 has been chosen in order to test the effect of moving the centroid x₂, increasing Δx₂₁ and decreasing Δx₃₂, in comparison with the symmetrical set #2. Sets #5 and #6 have been chosen because the effect of moving the centroid x₂, while keeping the extreme values x₁ and x₃, can be observed. In particular, set #5 allows us to analyze the effect of increasing x₁, with asymmetric distances Δx₂₁ = 40% and Δx₃₂ = 35%, while set #6 shows the effect of increasing x₂, with asymmetric distances Δx₂₁ = 55% and Δx₃₂ = 25%.

Although the sets of points could have been chosen considering other criteria, in the authors’ opinion, these are the ones that best reflect the effect of their position on the model’s accuracy.

In the following section, several examples will be used to verify that these observations are met.

Then, the procedure summarized in Figure 4 is followed for the asymmetrical case. Data sets {∆, α}, {α, a^α}, and {α, (τ_100−x)^1/α} for 0.5 ≤ α ≤ 1.1, and the functions f₁(Δ), f₂(α), and f₃(α) obtained by curve fitting for the different sets of asymmetrical points, i.e., (10-55-90%), (20-60-95%), and (20-75-95%), respectively, are shown in Figure 6. Note that Δ depends on normalized times {τ_x1, τ_x2, and τ_x3} and has a significant dependence on α, that f₂(α) = a^α depends on α and normalized times τ_x1 and τ_x3, and that f₃(α) = (τ_x3)^1/α depends on α and normalized time τ_100−x.

The Levenberg–Marquardt least-squares curve-fitting algorithm has been used for fitting data in all graphs in Figure 6. The same rational functions (27)–(29) have been used for obtaining functions f₁(Δ), f₂(α), and f₃(α), respectively, by using curve fitting.

The values of the corresponding parameters {p_i, q_i} for functions f₁(Δ), f₂(α), and f₃(α), and each one of the selected sets of asymmetrical points are shown in

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

(10-55-90%)	(20-60-95%)	(20-75-95%)
p1 = 0.3693	p1 = 0.4165	p1 = 0.3665
p2 = 9.55	p2 = 18.5	p2 = 7.191
p3 = −5.112	p3 = −15.5	p3 = −8.393
q1 = 9.53	q1 = 18.69	q1 = 5.838
q2 = −11.38	q2 = −25.6	q2 = −8.767

,

Table 8. Parameters {p_i, q_i} of the rational function f₂(α) for the different sets of asymmetrical points. Note that the parameters for (20-60-95%) and (20-75-95%) are the same since their x₁ and x₃-values are the same (x₁ = 20% and x₃ = 95%).

(10-55-90%)	(20-60-95%)	(20-75-95%)
p1 = −0.05698	p1 = −0.03498	p1 = −0.03498
p2 = 0.1596	p2 = 0.05957	p2 = 0.05957
q1 = −2.528	q1 = −2.33	q1 = −2.33
q2 = 1.753	q2 = 1.398	q2 = 1.398

and

Table 9. Parameters {p_i, q_i} of the rational function f₃(α) for the different sets of asymmetrical points. Note that the parameters for (20-60-95%) and (20-75-95%) are the same since their x₃-values are the same (x₃ = 95%).

(10-55-90%)	(20-60-95%)	(20-75-95%)
p1 = 4.518	p1 = 10.43	p1 = 10.43
p2 = −8.675	p2 = −22.09	p2 = −22.09
p3 = 5.666	p3 = 13.14	p3 = 13.14
q1 = −0.3471	q1 = −0.586	q1 = −0.586
q2 = 0.003197	q2 = 0.07943	q2 = 0.07943

, respectively.

Note that the obtained values of T, L, and α, which are determined by using Equations (25), depend on functions f₁(Δ), f₂(α), and f₃(α). These functions thus play a relevant role in the identification method as the features of normalized step responses (17) can be well characterized due to their respective contribution. It is important to emphasize that, for any different choice of times set {t_x1 t_x2, t_x3} to reach three points {y_α(t_x1), y_α(t_x2), 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 α.

4. Simulation Results

In this section, the identification method proposed in this work has been proved for several models that exhibit fractional behavior.

Process models (30), (31), and (32) have been selected to test the effectiveness of the proposed method in obtaining an FFOPDT model (13) in comparison with several identification methods for integer- and fractional-order models.

$${\mathrm{P}}_1(\mathrm{s})=\frac{1}{{(1+{\mathrm{s}}^{0.75})}^2}{\mathrm{e}}^{-0.1\mathrm{s}}$$

(30)

$${\mathrm{P}}_2(\mathrm{s})=\frac{3}{(1+3{\mathrm{s}}^{0.88})(1+2{\mathrm{s}}^{0.88})(1+{\mathrm{s}}^{0.88})}$$

(31)

$${\mathrm{P}}_3(\mathrm{s})={\mathrm{e}}^{-\sqrt{\mathrm{s}}}$$

(32)

On the one hand, process P₁ has been used to evaluate the proposed identification method for several symmetrical and asymmetrical sets of points in comparison with several well-known identification methods for integer-order models and to obtain insight into the influence of the location of asymmetrical points on the accuracy of the identified model.

On the other hand, processes P₂ and P₃ have been used to evaluate the model performance of the proposed procedure for symmetrical and asymmetrical points with other fractional-order methods.

These examples selected in this section have been also used to validate the proposed identification method for both symmetrical and asymmetrical points on the process reaction curve.

The experimental procedure followed is as follows: A step signal has been applied to the input of these processes and the reaction curve has been registered. Then the process output responses have been used to obtain the parameters of the models using the different identification methods. The sampling period used in all the experiments is Ts = 10 ms.

Finally, it is necessary to evaluate the accuracy of the model parameters that have been identified and the effectiveness of the model structure that has been adopted. A wide variety of model validation methods and fitting objective functions for system identification are available in the technical literature [5].

Although another objective function could have been used, the Mean Squared Error (MSE) has been used as a time-domain fitting criterion as a measure of performance for the identified model:

$$\mathrm{S}(\overline{\mathsf{\theta }})=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}{\displaystyle \sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{s}}}{[\mathrm{e}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})]}^2=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}{\displaystyle \sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{s}}}{[\mathrm{y}({\mathrm{kT}}_{\mathrm{s}})-{\mathrm{y}}_{\mathrm{m}}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})]}^2$$

(33)

where $\mathrm{e}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})$ is the difference between the process reaction curve and the step response of the identified model, $\mathrm{y}({\mathrm{kT}}_{\mathrm{s}})$ and ${\mathrm{y}}_{\mathrm{m}}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})$, respectively, $\overline{\mathsf{\theta }}$ is the vector of process model parameters, N_S is the number of collected samples, T_S is the sampling period, and N_ST_S is the time length of the dynamic (transient) response.

The Mean Absolute Error (MAE), the expression of which is (34), has also been calculated in the following examples for illustrative purposes.

$$\mathrm{E}(\overline{\mathsf{\theta }})=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}{\displaystyle \sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{s}}}\left|\mathrm{e}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})\right|=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}{\displaystyle \sum }_{\mathrm{k}=1}^{{\mathrm{N}}_{\mathrm{s}}}\left|\mathrm{y}({\mathrm{kT}}_{\mathrm{s}})-{\mathrm{y}}_{\mathrm{m}}({\mathrm{kT}}_{\mathrm{s}},\overline{\mathsf{\theta }})\right|$$

(34)

In the same context, it may be interesting to evaluate the goodness of fit of the identified model at different intervals of the reaction curve. For this reason, the index ${\mathrm{S}}_{{\mathrm{x}}_{\mathrm{i}}–{\mathrm{x}}_{\mathrm{j}}}(\overline{\mathsf{\theta }})$ is introduced in this paper, where $\overline{\mathsf{\theta }}$ is the vector of identified model parameters. This performance index represents the time-domain performance index $\mathrm{S}(\overline{\mathsf{\theta }})$ restricted to the interval [x_i–x_j] on the reaction curve, where x_i and x_j are two specific points on the reaction curve as depicted in Figure 3. Note that if the step response of the identified model is divided into p intervals [x_i–x_j], the model performance index $\mathrm{S}(\overline{\mathsf{\theta }})$ and the p performance indices ${\mathrm{S}}_{\mathrm{k}}(\overline{\mathsf{\theta }})$ (k = 1, …, p) of each of the intervals satisfy the following expression:

$$\mathrm{S}(\overline{\mathsf{\theta }})=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{p}}{\mathrm{S}}_{\mathrm{k}}(\overline{\mathsf{\theta }})·{\mathrm{N}}_{\mathrm{sk}}$$

(35)

where p $\in {\mathbb{Z}}^{+}$ is the number of intervals into which the model step response is divided, and N_Sk is the number of samples of the model step response in the corresponding k interval.

The simulation results obtained in this section have been performed using the FOTF MATLAB toolbox, which is a set of built-in functions that extends the control toolbox to deal with fractional-order systems [43]. For a deeper knowledge of the FOTF toolbox, the reader is referred to the reference text in [44].

4.1. Example 1

In this example, the fractional-order process model (30) is selected. This model is a lag-dominated fractional second-order process plus dead-time (FSOPDT), providing some modelling deviation from the model structure selected in the proposed identification method, which is the FFOPDT dynamics.

This same model but for different α-values in the range 0.60 ≤ α ≤ 1.00 has been used in [39], on the one hand, as a batch of processes to validate the proposed identification procedure for different sets of symmetrical points and, on the other hand, to obtain insight into the influence of the location of symmetrical points on the accuracy of the identified model.

The objective in this case is to validate the identification procedure for asymmetrical points on the process reaction curve and to obtain insight into the influence of the location of such asymmetrical points on the accuracy of the fractional-order identified model. A comparison between the proposed method and several well-known identification methods for integer-order models is also provided in this section.

The process reaction curve for this model and the step-input signal are shown in Figure 7.

The process information summarized in

Table 10. Process information collected from the reaction curve for fractional-order model identification of process P₁.

Symmetrical Methods			Asymmetrical Methods
Method #1: (5-50-95%)	Method #2: (10-50-90%)	Method #3: (25-50-75%)	Method #4: (10-55-90%)	Method #5: (20-60-95%)	Method #6: (20-75-95%)
Δu = 1.00
Δy = 1.00
t5 = 0.3020 s	t10 = 0.4540 s	t25 = 0.9300 s	t10 = 0.4540 s	t20 = 0.7620 s	t20 = 0.7620 s
t50 = 2.0910 s	t50 = 2.0910 s	t50 = 2.0910 s	t55 = 2.4400 s	t60 = 2.8590 s	t75 = 4.9730 s
t95 = 29.1410 s	t90 = 13.2640 s	t75 = 4.9730 s	t90 = 13.2640 s	t95 = 29.1410 s	t95 = 29.1410 s

for the proposed identification method is collected from data in Figure 7.

With the information collected from Table 10, the following FFOPDT model parameters ${\overline{\mathsf{\theta }}}_{1,\mathrm{i}}=\{{\mathrm{K}}_{1,\mathrm{i}},{\mathrm{T}}_{1,\mathrm{i}},{\mathrm{L}}_{1,\mathrm{i}},{\mathsf{\alpha }}_{1,\mathrm{i}}\}$, for i = 1, …, 6, have been obtained in

Table 11. Fractional-order model settings for the symmetrical and asymmetrical sets of points (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively, applied to Example 1.

Symmetrical Methods			Asymmetrical Methods
Method #1: (5-50-95%)	Method #2: (10-50-90%)	Method #3: (25-50-75%)	Method #4: (10-55-90%)	Method #5: (20-60-95%)	Method #6: (20-75-95%)
K1,1 = 1.0000	K1,2 = 1.0000	K1,3 = 1.0000	K1,4 = 1.0000	K1,5 = 1.0000	K1,6 = 1.0000
T1,1 = 2.3137 s	T1,2 = 2.2492 s	T1,3 = 2.1890 s	T1,4 = 2.1963 s	T1,5 = 2.0223 s	T1,6 = 1.8807 s
L1,1 = 0.2745 s	L1,2 = 0.2782 s	L1,3 = 0.3901 s	L1,4 = 0.2791 s	L1,5 = 0.3139 s	L1,6 = 0.2921 s
α1,1 = 0.7791	α1,2 = 0.7888	α1,3 = 0.8088	α1,4 = 0.7836	α1,5 = 0.7580	α1,6 = 0.7463

for the different sets of points proposed in Table 2 and Table 6.

The FFOPDT model step responses for (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively, are compared with the process reaction curve and illustrated in Figure 8, Figure 9 and Figure 10. In each figure the corresponding representative points, symmetrical (x-50-(100−x)%) and asymmetrical (x₁-x₂-x₃%), respectively, are also displayed.

Figure 8, Figure 9 and Figure 10 show that the step responses of the fractional-order models identified with the proposed method, for both the symmetrical and asymmetrical case, give good fit with the process reaction curve, which confirms the validity of this identification method also for the asymmetrical case. Note that the identification method for the symmetrical case has been studied in detail and its validity for the identification of FFOPDT models has also been confirmed in [39].

Table 12. FOPDT, DPPDT, and SOPDT model settings obtained for the considered integer-order identification methods.

FOPDT		DPPDT		SOPDT
Alfaro [11] (25–75%)	Vitecková [14] (33–70%)	Alfaro [11] (25–75%)	Vitecková [14] (33–70%)	Stark [16] (15-45-75%)	Jahanmiri–Fallahi [15] (2-70-90%)
K = 1.00	K = 1.00	K = 1.00	K = 1.00	K = 1.0000	K1,6 = 1.0000
T = 3.68 s	T = 3.51 s	T = 2.24 s	T = 2.33 s	T1 = 3.52 s	T1 = 5.61 s
L = 0.00 s	L = 0.00 s	L = 0.00 s	L = 0.00 s	T2 = 0.34 s	T2 = 0.0072 s
-	-	-	-	L = 0.00 s	L = 0.20 s

shows the process parameters identified for FOPDT and DPPDT models obtained using the methods proposed by Alfaro in [11], and by Vitecková et al. in [14], and the ones for SOPDT obtained using methods proposed by Stark in [16] and by Jahanmiri and Fallahi in [15], using two or three points from the process reaction curve.

The corresponding integer-order model step responses for the considered classical identification methods are compared with the process reaction curve and illustrated in Figure 11 for FOPDT and DPPDT models and in Figure 12 for SOPDT models.

Figure 11 and Figure 12 illustrate that FOPDT, DPPDT, and SOPDT models, respectively, approximate process P₁ with insufficient accuracy compared to FFOPDT models obtained with the proposed method. It has been illustrated from Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 that the proposed identification method outperforms significantly the considered methods for integer-order models.

Table 13. Comparison between the performance indexes obtained for the proposed method with different sets of points, symmetrical and asymmetrical, and the ones for several well-known methods for integer-order models. The number of samples N_S for the whole range and for each interval are also displayed.

i	Method	Set of Points	S0–50$({\overline{\mathsf{\theta }}}_{1,i})$	S50–100$({\overline{\mathsf{\theta }}}_{1,i})$	S0–100$({\overline{\mathsf{\theta }}}_{1,i})=\mathrm{S}({\overline{\mathsf{\theta }}}_{1,i})$	$\mathbf{E}({\overline{\mathsf{\theta }}}_{1,i})$
1	FFOPDT Proposed #1	(5-50-95%)	3.07 × 10−4	2.24 × 10−5	2.48 × 10−5	2.4 × 10−3
2	FFOPDT Proposed #2	(10-50-90%)	4.63 × 10−4	1.45 × 10−5	1.83 × 10−5	3.6 × 10−3
3	FFOPDT Proposed #3	(25-50-75%)	4.75 × 10−4	5.36 × 10−5	5.72 × 10−5	6.7 × 10−3
4	FFOPDT Proposed #4	(10-55-90%)	6.86 × 10−4	1.09 × 10−5	1.65 × 10−5	3.2 × 10−3
5	FFOPDT Proposed #5	(20-60-95%)	1.30 × 10−3	3.46 × 10−6	1.44 × 10−5	1.2 × 10−3
6	FFOPDT Proposed #6	(20-75-95%)	3.20 × 10−3	6.86 × 10−6	3.34 × 10−5	9.7 × 10−4
7	FOPDT Alfaro	(25-75%)	-	-	7.32 × 10−4	2.19 × 10−2
8	FOPDT Vitecková	(33-70%)	-	-	7.54 × 10−4	2.20 × 10−2
9	DPPDT Alfaro	(25-75%)	-	-	1.50 × 10−3	2.51 × 10−2
10	DPPDT Vitecková	(33-70%)	-	-	1.60 × 10−3	2.53 × 10−2
11	SOPDT Stark	(15-45-75%)	-	-	8.26 × 10−4	2.26 × 10−2
12	SOPDT Jahanmiri–Fallahi	(2-70-90%)	-	-	1.40 × 10−3	2.35 × 10−2
Number of samples			NS1 = 210	NS2 = 24,791	NS = NS1 + NS2 = 25,001	NS = 25,001

shows the values of the time-domain performance indexes $\mathrm{S}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$ and $\mathrm{E}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$ for process P₁, and the ones corresponding to the intervals [0–50%], ${\mathrm{S}}_{0-50}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$, and [50–100%], ${\mathrm{S}}_{50-100}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$, of the total process output change, respectively, for the different models considered in this example. In this table, i = 1, …, 6 represents the different sets of points considered in the proposed identification method, and i = 7, …, 12 represents the different identification methods for integer-order models.

Time-domain performance indexes in both intervals, ${\mathrm{S}}_{0-50}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$ and ${\mathrm{S}}_{50-100}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})$, give information about the effect of the location of the different set of points on the accuracy of the identified models. Since the step response of the identified models can be divided into the two aforementioned intervals, the accuracy of the identified models can be quantified and the one for the first and second half of the total interval can be also determined from data in Table 13. The number of samples for the model step responses in the whole range, N_S, and in both intervals, N_S1 and N_S2, respectively, are also provided in Table 13.

Note that from expression (35) the following relation that must be fulfilled can be particularized for this case:

$$\mathrm{S}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})={\mathrm{S}}_{0-100}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}\left[{\mathrm{S}}_{0-50}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}1}+{\mathrm{S}}_{50-100}({\overline{\mathsf{\theta }}}_{1,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}2}\right]$$

(36)

The results of Table 13 for the proposed identification method in terms of S_0–50(${\overline{\mathsf{\theta }}}_{1,\mathrm{i}}$), S_50–100(${\overline{\mathsf{\theta }}}_{1,\mathrm{i}}$), and S(${\overline{\mathsf{\theta }}}_{1,\mathrm{i}}$), for i = 1, …, 6, are also shown graphically in Figure 13.

Next, the results obtained in Table 13 and illustrated in Figure 13 for the FFOPDT models identified for process P₁ will be analyzed in order to gain insights into the location of the points for the asymmetrical case and their effect on the accuracy of the identified model.

The following observations can be drawn from the above results:

• Methods with low values of x₁ for the asymmetrical case, or x in the symmetrical case, allow obtaining low values of S_0–50, as can be seen in Figure 13 for methods #1, #2, and #4, compared to #5 and #6, which give worse results in that interval.
• Methods with high values of x₃ for the asymmetrical case allow better fitting of the model to the reaction curve in the final part of the response, as verified by the low values of S_50–100 for methods #4, #5, and #6. For the symmetrical case, a high value of (100 − x) implies a low value of x, thus better fitting the model to the reaction curve in the initial and final parts of the response, as illustrated by the low values of S_0–50 and S_50–100, respectively, for methods #1 and #2.
• Comparing methods #2 and #4, it can be seen that the effect of increasing the value of x₂, while keeping x₁ and x₃ constant, is to reduce the value of S as a result of a reduction in the value of S_50–100, even though the value of S_0–50 is slightly increased.
• Comparing methods #5 and #6, the effect of increasing the value of x₂ can also be observed, while keeping x₁ and x₃ constant. Note that the value of x₁ in this case is higher than for methods #2 and #4. The value of S for method #5 is reduced as a result of a very good fit in the interval [50–100%], as shown in Figure 13. The value of S corresponding to method #6 increases substantially due to a poorer fit in the interval [0–50%] despite the good result for S_50–100.

As discussed previously, the step response of the identified model can be divided into different intervals to evaluate the effect that the selection of different sets of points (x₁-x₂-x₃%) has on the model’s accuracy in each interval.

In this regard,

Table 14. Comparison between methods #2 and #4 in terms of time-domain model performance indices for different intervals of the process reaction curve for both identified models. The time-domain performance index for the whole response is also included in this table. N_S represents the number of points for each interval.

Interval #	Interval	Method #2: (10-50-90%)	Method #4: (10-55-90%)	NS
1	[0–10%]	${\mathrm{S}}_{0-10}({\overline{\mathsf{\theta }}}_{1,2})$ = 2.61 × 10−4	${\mathrm{S}}_{0-10}({\overline{\mathsf{\theta }}}_{1,4})$ = 2.64 × 10−4	46
2	[10–55%]	${\mathrm{S}}_{10-55}({\overline{\mathsf{\theta }}}_{1,2})$ = 4.31 × 10−4	${\mathrm{S}}_{10-55}({\overline{\mathsf{\theta }}}_{1,4})$ = 6.76 × 10−4	198
3	[55–90%]	${\mathrm{S}}_{55-90}({\overline{\mathsf{\theta }}}_{1,2})$ = 8.11 × 10−5	${\mathrm{S}}_{55-90}({\overline{\mathsf{\theta }}}_{1,4})$ = 4.90 × 10−5	1117
4	[90–100%]	${\mathrm{S}}_{90-110}({\overline{\mathsf{\theta }}}_{1,2})$ = 1.14 × 10−5	${\mathrm{S}}_{90-100}({\overline{\mathsf{\theta }}}_{1,4})$ = 9.06 × 10−6	23,674
-	[0–100%]	$\mathrm{S}({\overline{\mathsf{\theta }}}_{1,2})$ = 1.83 × 10−5	$\mathrm{S}({\overline{\mathsf{\theta }}}_{1,4})$ = 1.65 × 10−5	25,001

allows the comparison of method #2 with #4 for the intervals [0–10%], [10–55%], [55–90%], and [90–100%], respectively. The comparison of methods #5 and #6 for the intervals [0–20%], [20–50%], [50–75%], and [75–100%], respectively, is shown in

Table 15. Comparison between methods #5 and #6 in terms of time-domain model performance indices for different intervals of the process reaction curve for both identified models. The time-domain performance index for the whole response is also included in this table. N_S represents the number of points for each interval.

Interval #	Interval	Method #5: (20-60-95%)	Method #6: (20-75-95%)	NS
1	[0–20%]	${\mathrm{S}}_{0-20}({\overline{\mathsf{\theta }}}_{1,5})$ = 7.65 × 10−4	${\mathrm{S}}_{0-20}$ = 1.60 × 10−3	77
2	[20–50%]	${\mathrm{S}}_{20-50}({\overline{\mathsf{\theta }}}_{1,5})$ = 1.60 × 10−3	${\mathrm{S}}_{20-50}({\overline{\mathsf{\theta }}}_{1,6})$ = 4.10 × 10−3	133
3	[50–75%]	${\mathrm{S}}_{50-75}({\overline{\mathsf{\theta }}}_{1,5})$ = 9.31 × 10−5	${\mathrm{S}}_{50-75}({\overline{\mathsf{\theta }}}_{1,6})$ = 5.69 × 10−4	288
4	[75–100%]	${\mathrm{S}}_{75-100}({\overline{\mathsf{\theta }}}_{1,5})$ = 2.40 × 10−6	${\mathrm{S}}_{75-100}({\overline{\mathsf{\theta }}}_{1,6})$ = 2.45 × 10−7	24,503
-	[0–100%]	$\mathrm{S}({\overline{\mathsf{\theta }}}_{1,5})$ = 1.44 × 10−5	$\mathrm{S}({\overline{\mathsf{\theta }}}_{1,6})$ = 3.34 × 10−5	25,001

.

4.2. Example 2

The model (31) proposed in [31] is considered in this example. This process exhibits a higher-order lag-dominated fractional-order dynamics. The proposed identification method for symmetrical and asymmetrical sets of points are compared with other well-recognized methods for FFOPDT models, which are also based on the process reaction curve, and the optimal FOPDT.

The step-input signal and the process reaction curve for this model are shown in Figure 14.

The process information summarized in

Table 16. Process information collected from the reaction curve for fractional-order model identification of process P₂.

Symmetrical Methods			Asymmetrical Methods
Method #1: (5-50-95%)	Method #2: (10-50-90%)	Method #3: (25-50-75%)	Method #4: (10-55-90%)	Method #5: (20-60-95%)	Method #6: (20-75-95%)
Δu = 1.00
Δy = 3.00
t5 = 1.4140 s	t10 = 2.0290 s	t25 = 3.5600 s	t10 = 2.0290 s	t20 = 3.0630 s	t20 = 3.0630 s
t50 = 6.3850 s	t50 = 6.3850 s	t50 = 6.3850 s	t55 = 7.1044 s	t60 = 7.9180 s	t75 = 11.4060 s
t95 = 33.9290 s	t90 = 20.8030 s	t75 = 11.4060 s	t90 = 20.8030 s	t95 = 33.9290 s	t95 = 33.9290 s

for the different identification methods is collected from data in Figure 14.

With the information collected from Table 16, the following FFOPDT model parameters ${\overline{\mathsf{\theta }}}_{2,\mathrm{i}}=\{{\mathrm{K}}_{2,\mathrm{i}},{\mathrm{T}}_{2,\mathrm{i}},{\mathrm{L}}_{2,\mathrm{i}},{\mathsf{\alpha }}_{2,\mathrm{i}}\}$, for i = 1, …, 6, have been obtained in

Table 17. Fractional-order model settings for the symmetrical and asymmetrical sets of points (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively, applied to Example 2.

Symmetrical Methods			Asymmetrical Methods
Method #1: (5-50-95%)	Method #2: (10-50-90%)	Method #3: (25-50-75%)	Method #4: (10-55-90%)	Method #5: (20-60-95%)	Method #6: (20-75-95%)
K2,1 = 3.0000	K2,2 = 3.0000	K2,3 = 3.0000	K2,4 = 3.0000	K2,5 = 3.0000	K2,6 = 3.0000
T2,1 = 6.4066 s	T2,2 = 6.6381 s	T2,3 = 6.6817 s	T2,4 = 6.3850 s	T2,5 = 5.1733 s	T2,6 = 4.6914 s
L2,1 = 1.2638 s	L2,2 = 1.3955 s	L2,3 = 1.6203 s	L2,4 = 1.4329 s	L2,5 = 2.4042 s	L2,6 = 2.5096 s
α2,1 = 0.9189	α2,2 = 0.9470	α2,3 = 0.9802	α2,4 = 0.9391	α2,5 = 0.8901	α2,6 = 0.8759

for the different sets of points, i.e., (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively.

Process (31) is also approximated by an FFOPDT model following the method proposed by Tavakoli-Kakhki in [31], by an FFOPDT model obtained using the optimization-based method proposed by Guevara et al. in [34], and by the optimal FOPDT model, where model parameters, ${\overline{\mathsf{\theta }}}_{2,\mathrm{i}}$ for i = 7, 8, and 9, respectively, are given in

Table 18. FFOPDT model parameters obtained using methods proposed by Tavakoli-Kakhki [31] and by Guevara et al. [34], and optimal FOPDT model parameters, respectively, for process P₂.

Method #7: FFOPDT Tavakoli-Kakhki [31]	Method #8: FFOPDT Guevara et al. [34]	Method #9: FOPDT optimal
K2,7 = 3.00	K2,8 = 3.0000	K2,9 = 3.0000
T2,7 = 6.30 s	T2,8 = 5.6285 s	T2,9 = 8.7412 s
L2,7 = 1.00 s	L2,8 = 1.8833 s	L2,9 = 0.0000 s
α2,7 = 0.92	α2,8 = 0.9263	-

.

Note that the method proposed in [34] is based on optimization, where the function to be minimized is $\mathrm{E}(\overline{\mathsf{\theta }})$ (34). In this method, the approximation of the fractional term s^α is developed using the Oustaloup method; see [19]. In contrast, the parameters obtained in the optimal FOPDT model are those that minimize the function $\mathrm{S}(\overline{\mathsf{\theta }})$ (33).

The step responses of the considered approximated models for (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively, are compared with the process reaction curve and illustrated in Figure 15, Figure 16 and Figure 17. The corresponding representative points, symmetrical (x-50-(100 − x)%) and asymmetrical (x₁-x₂-x₃%), respectively, are also displayed in these figures. Moreover, the process reaction curve and the step responses of the FFOPDT model proposed by Tavakoli-Kakhki, the FFOPDT model obtained using the method proposed by Guevara et al., and the optimal FOPDT model, respectively, are shown in Figure 18.

In this example, a high-order process model has been utilized to demonstrate the effectivity and applicability of the proposed identification procedure in the task of estimating FFOPDT-type model parameters.

Similar to the previous example, Figure 15, Figure 16 and Figure 17 illustrate that the models obtained with the proposed method, for both the symmetrical and asymmetrical case, provide a good fit with the process reaction curve also for a higher-order lag-dominated fractional-order process in comparison with the results obtained using the well-recognized method proposed by Tavakoli-Kakhki, and even in comparison with optimization-based methods for FFOPDT and FOPDT models. These results are shown in Figure 18.

Figure 15a,b and Figure 16a show that the proposed method for the symmetrical case gives good results in the interval [x-(100 − x)%]. For this reason, better results are obtained for low values of x, as discussed in [39]. Figure 16b and Figure 17a,b show that by moving the points x₁, x₂, and x₃ of the reaction curve it is possible to obtain an FFOPDT model of which the step response can be better fitted to the reaction curve in certain intervals, as will be shown below.

Table 19. Comparison between the performance indexes obtained for the proposed method with different sets of points, symmetrical and asymmetrical, and the ones for several methods for FFOPDT and FOPDT models. The number of samples N_S for the whole range and for each interval are also displayed.

i	Method	Set of Points	S0–50$({\overline{\mathsf{\theta }}}_{2,i})$	S50–100$({\overline{\mathsf{\theta }}}_{2,i})$	S0–100$({\overline{\mathsf{\theta }}}_{2,i})=\mathrm{S}({\overline{\mathsf{\theta }}}_{2,i})$	$\mathbf{E}({\overline{\mathsf{\theta }}}_{2,i})$
1	FFOPDT Proposed #1	(5-50-95%)	5.80 × 10−3	8.78 × 10−4	1.10 × 10−3	2.40 × 10−2
2	FFOPDT Proposed #2	(10-50-90%)	2.80 × 10−3	1.20 × 10−3	1.30 × 10−3	3.34 × 10−2
3	FFOPDT Proposed #3	(25-50-75%)	4.40 × 10−3	3.33 × 10−3	3.30 × 10−3	5.14 × 10−2
4	FFOPDT Proposed #4	(10-55-90%)	3.80 × 10−3	9.58 × 10−4	1.10 × 10−3	3.02 × 10−2
5	FFOPDT Proposed #5	(20-60-95%)	2.78 × 10−2	4.53 × 10−4	1.60 × 10−3	1.92 × 10−2
6	FFOPDT Proposed #6	(20-75-95%)	2.88 × 10−2	1.01 × 10−4	1.30 × 10−3	1.24 × 10−2
7	FFOPDT Tavakoli-Kakhki [31]	-	2.02 × 10−2	5.25 × 10−4	1.40 × 10−3	2.47 × 10−2
8	FFOPDT Guevara et al. [34]	-	8.30 × 10−3	8.23 × 10−4	1.10 × 10−3	2.82 × 10−2
9	FOPDT optimal	-	4.55 × 10−2	9.07 × 10−4	2.80 × 10−3	2.95 × 10−2
Number of samples			NS1 = 639	NS2 = 14362	NS = NS1 + NS2 = 15,001	NS = 15,001

shows the values of the time-domain performance indexes S($\overline{\mathsf{\theta }}$_2,i) and E($\overline{\mathsf{\theta }}$_2,i) for process P₂, and the ones corresponding to the intervals [0–50%], S_0–50($\overline{\mathsf{\theta }}$_2,i), and [50–100%], S_50–100($\overline{\mathsf{\theta }}$_2,i), of the total process output change, respectively, for the different models considered in this example. In this table, i = 1, …, 6 represents the different sets of points considered in the proposed identification method, and i = 7, 8, and 9 represents methods #7 and #8 for FFOPDT models, and #9 for FOPDT, respectively.

Note that from expression (35) the following relation that must be fulfilled can be particularized for this case:

$$\mathrm{S}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})={\mathrm{S}}_{0–100}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})=\frac{1}{{\mathrm{N}}_{\mathrm{s}}}[{\mathrm{S}}_{0–50}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}1}+{\mathrm{S}}_{50–100}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}2}]$$

(37)

Table 19 shows that the accuracy of models obtained with methods #1, #5, and #6 is better in terms of E than the one obtained for methods #7 and #8. The proposed method for all the sets of points except #3 provides models with lower values of E than for the optimal FOPDT model.

For the sake of simplicity in the interpretation of results, information taken from Table 19 is depicted in Figure 19.

Figure 19 illustrates that the accuracy of the models obtained with the proposed method for the symmetrical and asymmetrical case is similar to the one obtained using the method proposed by Guevara et al., and even better than that obtained with a well-known identification method such as the one proposed by Tavakoli-Kakhki, or than the optimal FOPDT model. Note that the method proposed by Guevara et al. is based on optimization, while the method proposed in this paper is analytical.

Specifically, the value of the time-domain model performance index of methods #1, #2, #4, and #6 are lower than the one proposed by Tavakoli-Kakhki, confirming the effectiveness of the proposed method for both the symmetrical and the asymmetrical case. Furthermore, the proposed method not only gives better results than the Tavakoli-Kakhki method in terms of S, but in the authors’ opinion it is easier to apply.

The S-value of the model obtained with method #5 is slightly higher than the one by Tavakoli-Kakhki, while the one obtained with method #3 is substantially higher than the rest. This fact has already been discussed above and also in [39] for the symmetrical case and the choice of this method among the symmetrical methods allows us to determine that the accuracy of models obtained is improved for low values of x in the symmetrical case, and high values of x₃ for the asymmetrical case.

Results in Table 19 confirm the same observations taken from Example 1:

• Methods with low values of x₁ for the asymmetrical case, or x in the symmetrical case, present low values of S_0–50, as can be seen in Figure 19 for methods #1, #2, and #4.
• Methods with high values of x₃ for the asymmetric case present a lower value of S_50–100, as can be observed for method #4, and especially for #5 and #6. For the symmetric case, methods #1 and #2, with low values of x and high values of (100 − x), present good values of S_0–50 and S_50–100, as expected.
• Comparison of methods #2 and #4 yields the same conclusions as for Example 1.
• The observations drawn from Example 1 for the comparison of methods #5 and #6 are also extensible to Example 2.

The comparison of method #7 with the methods providing the best results in terms of S is then performed by dividing the reaction curve and the respective step responses of the obtained models into different intervals.

It can be extracted from data taken from Table 19 that the best results in terms of S are those obtained with models determined using method #1 for the symmetrical case and #4 for the asymmetrical case. Then, the step responses of the models obtained using methods #1, #4, and #7 are divided into four intervals, i.e., [0–10%], [10–50%], [50–90%], and [90–100%], respectively.

Table 20. Comparison between methods #1, #4, and #7 in terms of time-domain model performance indexes for different intervals of the process reaction curve for both identified models. The time-domain performance index for the whole response is also included in this table. N_S represents the number of points for each interval.

Interval #	Interval	Method #1: (5-50-95%)	Method #4: (10-55-90%)	Method #7	NS
1	[0–10%]	${\mathrm{S}}_{0–10}({\overline{\mathsf{\theta }}}_{2,1})$ = 2.80 × 10−3	${\mathrm{S}}_{0–10}({\overline{\mathsf{\theta }}}_{2,4})$ = 5.70 × 10−3	${\mathrm{S}}_{0–10}({\overline{\mathsf{\theta }}}_{2,7})$ = 4.80 × 10−3	203
2	[10–50%]	${\mathrm{S}}_{10–50}({\overline{\mathsf{\theta }}}_{2,1})$ = 7.20 × 10−3	${\mathrm{S}}_{10–50}({\overline{\mathsf{\theta }}}_{2,4})$ = 2.90 × 10−3	${\mathrm{S}}_{10–50}({\overline{\mathsf{\theta }}}_{2,7})$ = 2.73 × 10−2	436
3	[50–90%]	${\mathrm{S}}_{50–90}({\overline{\mathsf{\theta }}}_{2,1})$ = 6.10 × 10−3	${\mathrm{S}}_{50–90}({\overline{\mathsf{\theta }}}_{2,4})$ = 1.30 × 10−3	${\mathrm{S}}_{50–90}({\overline{\mathsf{\theta }}}_{2,7})$ = 2.30 × 10−3	1442
4	[90–100%]	${\mathrm{S}}_{90–100}({\overline{\mathsf{\theta }}}_{2,1})$ = 3.00 × 10−4	${\mathrm{S}}_{90–100}({\overline{\mathsf{\theta }}}_{2,4})$ = 9.25 × 10−5	${\mathrm{S}}_{90–100}({\overline{\mathsf{\theta }}}_{2,7})$ = 3.26 × 10−4	12,920
-	[0–100%]	$\mathrm{S}({\overline{\mathsf{\theta }}}_{2,1})$ = 1.10 × 10−3	$\mathrm{S}({\overline{\mathsf{\theta }}}_{2,4})$ = 1.10 × 10−3	$\mathrm{S}({\overline{\mathsf{\theta }}}_{2,7})$ = 1.40 × 10−3	15,001

provides a comparison of these methods for the aforementioned intervals.

Note that from expression (35) the following relation must be fulfilled:

$$\begin{array}{ll}\mathrm{S}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})& ={\mathrm{S}}_{0-100}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})\\ & =\frac{1}{{\mathrm{N}}_{\mathrm{s}}}[{\mathrm{S}}_{0-10}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}1}+{\mathrm{S}}_{10-50}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}2}+{\mathrm{S}}_{50-90}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}3}\\ & +{\mathrm{S}}_{90-100}({\overline{\mathsf{\theta }}}_{2,\mathrm{i}})·{\mathrm{N}}_{\mathrm{s}4}]\end{array}$$

(38)

Figure 20 represents graphically the information taken from Table 20, which simplifies the interpretation of results.

Information about the accuracy of the model obtained at the different intervals of the reaction curve can be extracted from Table 20. Specifically, it can be seen that method #1 provides a model that fits the reaction curve very well in the intervals [0–10%], [10–50%] and [90–100%], while the fit is worse in the interval [50–90%], as can also be seen in Figure 15a. This method provides the best results of the three compared methods in the intervals [0–10%] and [90–100%].

The model for method #4 fits very well in the intervals [10–50%] and [50–90%], where it presents the lowest values of the three compared methods. In contrast, the value of S increases outside these intervals, as can also be seen in Figure 16b.

The model proposed by Tavakoli-Kakhki provides a good fit in the intervals [50–90%] and [90–100%], while the results for the rest of intervals are worse, particularly for [10–50%].

4.3. Example 3

The model (32), which has been proposed in [45] as part of a collection of systems that are suitable for testing PID controllers, is considered in this example.

Physically, this transfer function describes how temperature evolves over time in a solid medium, which assumes that the temperature is distributed in only one spatial coordinate and the heat is transferred in the direction in which the temperature decreases. This ideal transfer function contains the irrational operator $\sqrt{s}$ which may cause difficulties in a further analytical analysis.

In general, temperature control is a widespread application in the process industry, and, in particular, thermal conduction is a very common dynamic in process control.

Traditionally, the FOPDT model has been utilized to approximate this dynamic; however, it has insufficient accuracy. Since a growing body of evidence has suggested that fractional-order models are able to describe dynamic processes with higher accuracy, the proposed identification method will be used to more accurately model this process, which exhibits fractional behavior.

In order to evaluate the effectiveness of the proposed method, the estimated models for the symmetrical and asymmetrical cases will be compared with the optimal FOPDT and FFOPDT models.

Following the same procedure used in previous examples, the following FFOPDT model parameters $\overline{\mathsf{\theta }}$_3,i = {K_3,i, T_3,i, L_3,i, α_3,i} for i = 1, …, 6, have been obtained in

Table 21. Fractional-order model settings for the symmetrical and asymmetrical sets of points (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively, and optimal FFOPDT and FOPDT model parameters.

Method #1: (5-50-95%)	Method #2: (10-50-90%)	Method #3: (25-50-75%)	Method #4: (10-55-90%)	Method #5: (20-60-95%)	Method #6: (20-75-95%)	#7: FFOPDT Optimal	#8: FOPDT Optimal
K3,1 = 1.00	K3,2 = 1.00	K3,3 = 1.00	K3,4 = 1.00	K3,5 = 1.00	K3,6 = 1.00	K3,7 = 1.02	K3,8 = 0.91
T3,1 = 1.26 s	T3,2 = 1.26 s	T3,3 = 1.22 s	T3,4 = 1.22 s	T3,5 = 1.13 s	T3,6 = 1.08 s	T3,7 = 1.23 s	T3,8 = 2.24 s
L3,1 = 0.00 s	L3,2 = 0.60 s	L3,3 = 0.22 s	L3,4 = 0.67 s	L3,5 = 0.00 s	L3,6 = 0.00 s	L3,7 = 0.0001 s	L3,8 = 0.00 s
α3,1 = 0.5368	α3,2 = 0.5459	α3,3 = 0.5598	α3,4 = 0.5401	α3,5 = 0.5206	α3,1 = 0.5140	α3,7 = 0.5000	-

for the different sets of points, i.e., (5-50-95%), (10-50-90%), (25-50-75%), (10-55-90%), (20-60-95%), and (20-75-95%), respectively. Model parameters for optimal FFOPDT and FOPDT models have also been calculated. Process (32) has also been approximated by the optimal FFOPDT and FOPDT models, where model parameters $\overline{\mathsf{\theta }}$_3,i for i = 7 and 8, respectively, are given in Table 21. Note that the parameters obtained in the optimal models are those that minimize the function $\mathrm{S}(\overline{\mathsf{\theta }})$ (33). Since the Taylor series expansion for the irrational transfer function (32) give a polynomial of half-order integrators s^0.5, α_3,7 has been considered to be 0.5.

In this example, an ideal thermal conduction process model has been used to demonstrate the effectiveness of the proposed identification procedure in the task of estimating the FFOPDT model parameters, compared to the optimal FOPDT model, which is traditionally used to model this type of process. Also, the results are compared with the optimal FFOPDT model to verify that with the proposed method it is possible to obtain comparable results to the latter with much less computational effort.

The step responses of the considered approximated models in Table 21 are compared with the process reaction curve and illustrated in Figure 21, Figure 22, Figure 23 and Figure 24.

Figure 21, Figure 22 and Figure 23 illustrate that the models obtained with the proposed method, for both the symmetrical and asymmetrical case, provide a good fit with the process reaction curve compared to the results obtained using the optimal FFOPDT, which is shown in Figure 24a. It can also be observed that the proposed method outperforms the results obtained using the optimal FOPDT model, which is shown in Figure 24b.

Table 22. Comparison between the performance indexes obtained for the proposed method with different sets of points, symmetrical and asymmetrical, and the ones for optimal methods for FFOPDT and FOPDT models. The number of samples N_S is also displayed.

i	Method	Set of Points	$\mathbf{S}({\overline{\mathsf{\theta }}}_{3,\mathrm{i}})$	$\mathbf{E}({\overline{\mathsf{\theta }}}_{3,\mathrm{i}})$
1	FFOPDT Proposed #1	(5-50-95%)	1.00 × 10−4	3.00 × 10−3
2	FFOPDT Proposed #2	(10-50-90%)	3.60 × 10−4	5.70 × 10−3
3	FFOPDT Proposed #3	(25-50-75%)	8.97 × 10−5	8.90 × 10−3
4	FFOPDT Proposed #4	(10-55-90%)	4.19 × 10−4	5.40 × 10−3
5	FFOPDT Proposed #5	(20-60-95%)	1.41 × 10−4	2.00 × 10−3
6	FFOPDT Proposed #6	(20-75-95%)	1.73 × 10−4	2.50 × 10−3
7	FFOPDT optimal	-	1.47 × 10−4	6.40 × 10−3
8	FOPDT optimal	-	1.40 × 10−3	3.08 × 10−2
Number of samples			NS = 15,001	NS = 15,001

shows the values of the time-domain performance indexes S($\overline{\mathsf{\theta }}$_3,i) and E($\overline{\mathsf{\theta }}$_3,i) for process P₃, for the different models considered in this example. In this table, i = 1, …, 6 represents the different sets of points considered in the proposed identification method, and i = 7 and 8 represents optimal FFOPDT and FOPDT models, respectively.

Table 22 shows that methods #1, #5, and #6 provide E values comparable to that obtained for the optimal FFOPDT model. On the other hand, it is also observed that the FOPDT model approximates the process dynamics with insufficient accuracy. In fact, for example, methods #1, #5, and #6 have E values that are 10.3, 15.4, and 12.3 times lower than that obtained by the optimal FOPDT model.

5. Experimental Results

One of the objectives of this paper is to verify the applicability of the proposed identification procedure in a laboratory prototype and to obtain insight into the practical issues related to its implementation on industrial control hardware.

In this section, a thermal process-based hardware-in-the-loop experimental setup will be used to validate the effectiveness of the proposed identification procedure implemented on microprocessor-based control hardware.

This section is organized as follows. Firstly, the laboratory prototype to be used is described. After that, the controlled process is configured to test the identification procedures proposed in this work for both the symmetrical and the asymmetrical case. Since the prototype can be set up in three different configurations, this section describes the selected configuration and defines the different components and variables that make up the controlled process. Then, the hardware control architecture to be used for the operation of the prototype and the practical implementation of the model estimation procedures are briefly presented. Next, the identification procedures are applied to the thermal laboratory prototype using a microprocessor-based control hardware, confirming the applicability of these procedures in industrial equipment. This constitutes one of the main results of this paper. Finally, some remarks and final comments about practical issues in the microprocessor-based implementation are offered in this context.

5.1. Description of the Prototype

Recently, a thermal-based experimental setup has been designed and built in the Faculty of Engineering, University of Deusto. The considered prototype consists of two clearly different parts, as can be seen in the 3D-model layout in Figure 25.

• The upper part of the equipment, where there is a methacrylate duct, the head of a 3D-printer extruder, and an air fan in front of the hot end. User LEDs, an LCD display, a user button, and four BNC output connectors to display the main process variables on an oscilloscope, can be found outside of the enclosure.
• The inner part of the enclosure, where the power supply and all the hardware and electrical components necessary for the correct operation of the experimental setup can be found. The connection of the input and output signals to the control hardware has also been arranged through a standard 34-way IDC connector, which is placed on one side of the box.

A detailed description of this laboratory equipment, which is named Deusto Heater Experimental Setup (Deusto HES), can be found in [46].

5.2. Reconfigurable Controlled Process

The thermal process considered in this paper takes place on the upper part of the prototype. The equipment brings about the thermodynamic process of temperature control in a 3D-printer extruder head.

It should be noted that the head of the 3D-printer is not a conventional component, but the dimensions of the heat block have been conveniently modified. The head of the 3D-printer is just used as a heating element (final controlling element) and there is no type of extrusion process, as takes place in a real 3D-printer.

The control objective in this process is to achieve a temperature in the heat block T(t) that reaches the value set by the user T^SET, despite the disturbances that may take part in the process.

If it is assumed that the physical properties do not vary, only the effects of the heating power and the air fan on the controlled variable T(t) should be considered, as shown in Figure 25. Accordingly, from a thermodynamical point of view, the only source of heat in the system is represented by the heat conduction originated in the resistor hole, which exhibits a fractional behavior [47]. On the other hand, another phenomenon that takes place in the process is the forced convection on the heat block generated by the airflow due to the air fan [48,49]. As both phenomena can be controlled by the control hardware, from a control point of view the controlled process configuration can be set using software components, as will be discussed below.

In this paper, the heating element is used as the final control element while keeping constant the speed of the air fan. Figure 26 represents the block diagram of the controlled process in the selected configuration, while the different components, variables, and units are listed in

Table 23. Main process variables and components in the considered controlled process configuration.

Process Variables or Components	Controlled Process Configuration #1
Controlled variable	Temperature in the heat block T(t) [°C]
Manipulated variable	Power delivered to the heat block by the heating resistance P(t) [W]
Measured variables	Temperature measured by the thermistor Tm(t) [V] Rotational speed of fan ωF(t) [V]
Control signal	Output of the controller uH(t) [%]
Final control element	Heating resistance
Measurement devices	Temperature transmitter (TT) and Frequency transmitter (ST)
Disturbances	Ambient temperature Ta(t) [°C] and command signal to air fan uF(t) [%]

. However, two additional configurations are available with this equipment, constituting a reconfigurable controlled process, as is explained in detail in [50].

5.3. Control Hardware

In this paper, NI myRIO-1900 equipment was used as a control hardware device, although any other type of microprocessor or hardware device could be easily incorporated.

Figure 27 illustrates the scheme of the hardware architecture used for the operation of this prototype and for the practical implementation of control and model estimation algorithms. This hardware architecture has been proposed in [50] and the use of this control hardware for the implementation of integer- and fractional-order PID controllers is explained in detail. In particular, the design and experimental validation of the fractional-order controllers implemented in several control technologies were applied to this thermal-based experimental setup in order to demonstrate the effectiveness of the proposed hardware architecture.

The flexibility required in the industrial context is offered by this control hardware and the configuration of the proposed control architecture provides the following technologies for the practical implementation of model identification procedures:

• DAQ mode;
• Microprocessor-based mode;
• FPGA-based mode.

LabVIEW has been used as a programming language in this prototype. LabVIEW, which has become a de facto standard in industry, especially in the areas of industrial measurement and control, with the latest and most advanced programming techniques and transparency in the access to hardware devices, facilitates not only the application of different modelling techniques, but also the practical implementation of conventional and advanced control algorithms.

Depending on the control technology that is being used, the communication technologies between the control hardware and the local or remote PC are programmed with the following software components, as indicated in Figure 27 and detailed in [46].

• LabVIEW;
• LabVIEW RT;
• LabVIEW RT and LabVIEW FPGA

Without loss of generality, in this paper only the microprocessor-based hardware mode has been used to operate the prototype and estimate fractional-order model parameters applying the proposed identification procedure.

5.4. Model Estimation

For the purposes of this paper, it is enough to deal with the static and dynamic characteristics of the controlled process, estimating an FFOPDT model from the process reaction curve at a certain operating point, as will be considered in this section. To do this, a LabVIEW-based application has been implemented.

Consider now that an open-loop step-test experiment is applied to the aforementioned thermal laboratory setup with the controlled process configuration detailed in Figure 26.

Initially, the fan and the control signal to the heating element are at u_F = 10% and u_H = 30%, respectively.

At instant t = 0 s, a step change of amplitude Δu_H = 30% is applied, from u_H = 30% to 60%, while the value of u_F is kept constant, as shown in Figure 28a. The measured temperature in the heat block, which ranges from 60.5 to 102.5 °C (ΔT_m = 42 °C), is recorded in Figure 28b. This output is a noisy signal and that noise is reduced with a first-order low-pass filter that usually is used for filtering purposes. The sample time used in this experiment was T_S = 100 ms.

From the filtered output signal shown in Figure 28b, the process information in

Table 24. Process information for methods #2 and #4, respectively, collected from the process reaction curve for fractional-order model identification of the thermal process.

Symmetrical	Asymmetrical
Method #2: (10-50-90%)	Method #4: (10-55-90%)
Δu = ΔuH = 30%
Δy = ΔTm = 42 °C
t10 = 16.8000 s
t50 = 53.3000 s	t55 = 59.9000 s
t90 = 174.5000 s

is obtained.

Table 25. Model parameters for the thermal process at the considered operating point obtained using the identification methods #2 and #4 for FFOPDT models, using methods proposed by Alfaro and Vitecková for FOPDT models, and using methods proposed by Stark and Jahanmiri and Fallahi for SOPDT models, respectively.

Method #2 (10-50-90%)	Method #4 (10-55-90%)	FOPDT [11] (25–75%)	FOPDT [14] (33–70%)	SOPDT [16] (15-45-75%)	SOPDT [15] (2-70-90%)
K4,1 = 1.40 °C/%	K4,2 = 1.40 °C/%	K4,3 = 1.40 °C/%	K4,4 = 1.40 °C/%	K4,5 = 1.40 °C/%	K4,6 = 1.40 °C/%
T4,1 = 49.52 s	T4,2 = 48.43 s	T4,3 = 64.26 s	T3,4 = 63.37 s	T4a,5 = 59.92 s	T3,4 = 70.92 s
L4,1 = 11.51 s	L4,2 = 11.64 s	L4,3 = 10.20 s	L3,4 = 10.25 s	T4b,5 = 13.68 s	T3,4 = 0.092 s
α4,1 = 0.9462	α4,2 = 0.9430	-	-	L4,5 = 0.00 s	L4,6 = 9.10 s

contains the different model parameters for the thermal process at the considered operating point obtained using the identification methods #2 and #4 for FFOPDT models, using methods proposed by Alfaro in [11] and Vitecková et al. in [14] for FOPDT models, and using methods proposed by Stark in [16] and Jahanmiri and Fallahi in [15] for SOPDT models, respectively.

In this example, the thermal-based experimental setup has been utilized in order to illustrate the effectivity and applicability of the proposed identification method in an industrial environment.

The step responses of the estimated models for the symmetrical (10-50-90%) and asymmetrical (10-55-90%) case, and the ones for FOPDT and SOPDT models obtained using the considered classical methods, are compared with the process reaction curve and illustrated in Figure 29, Figure 30 and Figure 31, respectively. The corresponding representative points are also displayed in these figures.

Figure 29 illustrates that the FFOPDT models obtained using methods #2 and #4 give a good fit to the thermal process reaction curve in comparison to the results obtained for the FOPDT and SOPDT models, which are shown in Figure 30 and Figure 31, respectively. Note that the computational effort of all the methods considered in this experimental example is similar since all are analytical methods based on the process reaction curve.

Table 26. Comparison of the approximated models obtained with methods #2 and #4 and models for FOPDT and SOPDT in terms of time-domain model performance indices for the thermal-based process at a certain operating point. The number of samples N_S is also displayed.

i	Identification Method	Set of Points	$\mathbf{S}({\overline{\mathsf{\theta }}}_{4,\mathit{i}})$	$\mathbf{E}({\overline{\mathsf{\theta }}}_{4,\mathit{i}})$
1	FFOPDT Proposed method #2	(10-50-90%)	7.59 × 10−5	6.90 × 10−3
2	FFOPDT Proposed method #4	(10-55-90%)	6.62 × 10−5	6.00 × 10−3
3	FOPDT Alfaro	(25-75%)	7.69 × 10−4	2.50 × 10−2
4	FOPDT Vitecková	(33-70%)	8.47 × 10−4	2.59 × 10−2
5	SOPDT Stark	(15-45-75%)	1.20 × 10−3	3.20 × 10−2
6	SOPDT Jahanmiri and Fallahi	(2-70-90%)	8.35 × 10−4	2.51 × 10−2
Number of samples			NS = 4001	NS = 4001

shows the values of the time-domain performance indexes S($\overline{\mathsf{\theta }}$_4,i) and E($\overline{\mathsf{\theta }}$_4,i) for the different approximated models considered in this example.

In this table, one can observe that the result obtained by method #4 is slightly better in terms of E than the one obtained using method #2. The results obtained by fractional-order identification methods outperform those obtained using integer-order models. Specifically, the E values obtained by using methods #2 and #4 are 3.6, 3.7, 4.6, and 3.6 times lower than those obtained by using the methods of Alfaro, Vitecková, Stark, and Jahanmiri and Fallahi, respectively.

LabVIEW-Based Implementation

As discussed previously, a LabVIEW-based application has been developed to operate the prototype and implement the proposed identification procedure, demonstrating its applicability on a microprocessor-based hardware device.

Figure 32 illustrates the front panel of the LabVIEW-based application for fractional-order model estimation using the proposed identification procedure, considering both the symmetrical and the asymmetrical case.

Once the identification method has been selected, and according to the identification procedure, the variation in input signal Δu and process output Δy, and times required to reach x₁% (t_x1), x₂% (t_x2), and x₃% (t_x3) of the total process output change on the reaction curve, must be taken in order to obtain FFOPDT model parameters θ_P = {K, T, L, α}.

The algorithm illustrated in Algorithm 1 has been developed in order to simplify software implementation of the identification procedure proposed in this paper.

Algorithm 1. Three-Point Identification Method for FFOPDT Models

Input: Selected symmetrical or asymmetrical method and {Δu, Δy, tx1, tx2, tx3} collected from the process reaction curve

Output: Fractional-order model parameters θP = {K, T, L, α}

1: Select the corresponding symmetrical or asymmetrical method.

2: Collect process data {Δu, Δy, tx1, tx2, tx3} from the process reaction curve.

3: Obtain the process gain K using (18).

4: Obtain the value of the times ratio Δ by using Equation (21).

5: Obtain the value of α = f1(Δ) by using Equation (27).

6: Obtain the value of functions f2(α) and f3(α) by using Equations (28) and (29), respectively.

7: Calculate the value of T using Equation (22).

8: Calculate the value of L using Equation (23).

This application presents the following features, which can be observed in Figure 32:

• Selection of the identification procedure, which presents the following options:
  • Symmetrical case (x-50-(100−x)%): (5-50-95%), (10-50-90%), or (25-50-75%).
  • Asymmetrical case (x₁-x₂-x₃%): (10-55-90%), (20-60-95%), or (20-75-95%).
• Determination of process data {Δu, Δy, t_x1, t_x2, t_x3} from the process reaction curve.
• Estimation of the FFOPDT model parameters: θ_P = {K, T, L, α}.
• Graphs for registering control signal u_H(t) [%], command signal to air fan u_F(t) [%], process reaction curve T_m(t) [°C], representative points of the process reaction curve {(t_x1, y_α(t_x1)), (t_x2, y_α(t_x2)), (t_x3, y_α(t_x3))}, and step response of the identified model.
• Export the experimental data in Excel or text-format.

5.5. Remarks and Final Comments

Practical application of the proposed fractional-order model identification algorithm on laboratory equipment allows gaining experience about practical issues related to its implementation on industrial control hardware. Some remarks and comments on industrial practice in this context are discussed below.

Remark 1: Measurement noise. The main disadvantage provided by the use of feedback is that measurement noise is injected into the loop. Noise generally generates undesirable motion of the final controlling elements, which may cause wear and possible breakdown.

It is common practice that open-loop model identification procedures use process information based on a noise-free process reaction curve; see, e.g., [5] for integer-order systems, and [31] for fractional-order systems. However, it is usual in an industrial environment that the controlled process feedback signal includes measurement noise that must be properly filtered for model identification and control purposes.

It is important to take into account that the filter dynamics will be an integral part of the controlled process to be identified. As this will influence and add a lag to the loop dynamics, a measurement filter must be set before any controlled process model identification and/or controller tuning [18].

In this context, it is very difficult to derive general conclusions about the influence of measurement noise on the identification procedure, because this will depend on the measurement noise and filter characteristics.

The results of the experiment conducted in this section serve to demonstrate the applicability of the proposed identification method on laboratory equipment and its implementation on industrial control hardware, which presents similarities with the industrial environment.

Remark 2: Model uncertainty. In the industrial context, processes exhibit nonlinearities, i.e., the dynamic characteristics of the process and therefore those of the identified FFOPDT model—gain, time constant, apparent dead-time and fractional order—will change with the operating point. The operating point of a control system may vary due to a change in the setpoint or as a result of the effect of disturbances.

Therefore, it must be considered that there is an implicit uncertainty in the nominal model.

In general, there are two approaches in the technical literature when considering the parametric uncertainty of the plant in an identification method; see, e.g., [39].

• The first approach is to incorporate the uncertainty explicitly into the model. This typically makes the identification procedure more complicated.
• The second approach consists of taking into account the potential changes in the controlled process dynamics and model uncertainties in the design phase of the controller; see, e.g., [2] for integer-order controllers and [19] for fractional-order controllers. A common application of this second approach is to ensure a certain degree of robustness of the designed control system to guarantee its stability under variations in the process characteristics.

In the context of this work, the primary use of the identified fractional-order model is for control purposes. Consequently, the approach to be used will be the second one.

Rules of thumb in selecting sets of points: In Section 3.1 and Section 3.2 some rules of thumb, which are based on experimentation, have been provided to obtain insight about the selection of the points for the symmetric and asymmetric case, respectively, in the context of the proposed identification method.

Based on the results obtained in examples 1–3 in Section 4, a recommendation on the selection of the set of points for the proposed identification method can be provided:

• In the symmetrical case, the set of points in method #1 gives the best results in terms of the performance index E.
• In the asymmetrical case, the set of points in methods #5 and #6 give quite similar results in terms of E, although #6 is slightly better.

Final comments and discussion: This section contains a discussion about the usefulness of, and the improvement obtained by, using an FFOPDT model instead of an integer-order model, with FOPDT, DPPDT, and SOPDT models being those most commonly used in the industrial context.

In this paper, an identification method for FFOPDT models has been presented. This method is focused on processes characterized by an S-shaped response (monotonic) that exhibits a fractional behavior.

There is no exact definition to describe a fractional dynamic behavior in the technical literature. According to [51], a system has a fractional behavior if its input and output are linked by a function of the form t^υ−1, υ $\in ℝ$, υ < 1.

Although an implicit link exists in the literature between fractional behaviors and fractional-based models, they are two distinct concepts. One designates a property or a particular behavior of a physical system, while the other refers to a model class that can capture fractional behaviors.

Fractional behaviors appear in many physical-, biological-, or thermal-based processes, among others [19,20]. In spite of the slow dynamics, real processes exhibiting fractional behavior with values of α in the range proposed in this paper can be found in electrical engineering, motion controls, and process control, with thermal-based processes being the most important type encountered in the process industry [47].

Since these kinds of behaviors are ubiquitous, the existence of methods for identification of simple-structure fractional-order models is of significant interest.

The main issue for adopting fractional-order models in industry can be summarized in the form of the following question: “Is the additional effort to consider a fractional-order model worth it, in order to obtain a more accurate model and, eventually, a better control performance?”.

In the industrial context, the apparent benefit of using fractional calculus has been justified in the literature in terms of a more precise modelling; see, e.g., [20,34,36].

More specifically in the context of this work, the proposed identification method for symmetrical and asymmetrical sets of points has been compared with several two- and three-point identification methods for integer-order models, which are based on data collected from the process reaction curve. It has been illustrated that the proposed method outperforms significantly methods for integer-order models. In some cases, even better results are obtained with the estimated fractional-order model than with the optimal integer-order model.

If the process reaction curve exhibits fractional behavior, the estimated FFOPDT model will more accurately fit the reaction curve compared to the integer-order model. Note that the proposed identification method allows characterizing the existence of fractional behavior in measured data collected from the process reaction curve.

The identification procedure presented in this paper is intended for control purposes, as considered previously. In the industrial context, large process industries have hundreds or thousands of control loops. In order for such an identification procedure to have a significant impact in the industrial environment, simplicity is a fundamental feature.

In the context of this work, there is a trade-off between accuracy and computational effort. Although this identification procedure provides good results in comparison with other well-known integer- and fractional-order identification methods, it is possible to find methods that improve further the accuracy of the estimated fractional-order model at a cost of more complex algorithms or a higher computational effort. 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.

6. Conclusions

In this paper, an identification procedure for FFOPDT models, which is based on the information obtained from three points on the process reaction curve, is presented. This identification procedure has been validated for both cases, considering three symmetrical and asymmetrical points on the reaction curve.

Three simulation examples, with processes exhibiting different dynamics, have been used to verify the simplicity and effectiveness of this identification procedure for the symmetrical and asymmetrical cases. Good results have been obtained in comparison with other well-recognized integer- and fractional-order identification methods which are also based on information taken from the process reaction curve, especially when simplicity is emphasized.

A thermal process-based experimental setup has also been used, where the proposed identification procedure has been implemented in a microprocessor-based control hardware, confirming the applicability of this method in an industrial equipment.

Besides effectiveness, the main characteristic exhibited by this method is simplicity. This feature is of significant importance in the industrial context, where easy-to-implement identification methods are required.

Some comments and reflections have been offered in the context of industrial practice.

It is worth pointing out that the identification method proposed in this paper is restricted to values of the fractional order in the range 0.5 ≤ α ≤ 1.0, and as a future work, this method is being extended using the same methodology for processes with underdamped step response, extending the range of the fractional order to 1.0 ≤ α ≤ 2.0.