Optimal Identification and Metaheuristic PID Control of a Two-Tank System

Abstract

In the modern optimization context, this paper introduces an optimal PID-based control strategy for a two-tank installation, namely ASTANK2. The process model was identified by using raw and spline smoothed measured data, respectively. Two PID controller configurations, a standard (regular) one (PID-R) and a non-standard one (PID-N), were considered for each type of model, resulting in four regulators. The optimal tuning parameters of each regulator were obtained by a searching approach relying on a combination of two metaheuristics. Firstly, an improved version of the Hill Climbing algorithm was employed to comprehensively explore the searching space, aiming to find fairly accurate tuning parameters. Secondly, an improved version of the Firefly Algorithm was proposed to intensively refine the search around the previously found optimal parameters. A comparative analysis between the four controllers was achieved in terms of performance and robustness. The simulation results showed that all optimal controllers yielded good performance in the presence of exogenous stochastic noise (bounded error tracking, setpoint tracking, reduced overshoot, short settling time). Robustness analysis is extensive and illustrates that the PID-R controllers are more robust to model uncertainties, whilst PID-N controllers are more robust to tracking staircase type references.

1. Introduction

Although there is a variety of complex control strategies (such as model predictive, sliding mode, adaptive control, neural and fuzzy control [1]), in control engineering practice, the Proportional–Integral–Derivative (PID) control structures are the most widely used (over 95% from industrial applications [2]), owing to their simple structure, good stability [3], high reliability, reasonable closed loop performance and robustness [4,5,6].

The key problem in designing a PID controller ensuring these features is to tune the parameters accurately and efficiently. This is a challenging task when dealing with industrial processes that usually are modelled as multi-variable systems, with nonlinear dynamics and a large number of tunable variables [5]. Adjusting controller parameters in Multi-Input–Multi-Output (MIMO) process can become a difficult problem, especially in case of strong coupling of the output channels, for which Single-Input–Single-Output (SISO) decoupled loops might affect the performance [4].

Over the years, the PID tuning methods evolved from empirical (Ziegler–Nichols [7,8], Cohen–Coon) and analytical (internal model, pole placement, gain-phase margin [9]) to adaptive and optimization-based [5]. Unlike the empirical and analytical methods that provide weak or only satisfying performance, the optimization techniques applied to PID tuning yield good performance (disturbance rejection, set-point tracking) and robustness [4,10,11]. Therefore, in recent decades, the problem of designing an optimal PID controller that drives the systems into a desired behavior has increasingly captured the interest of the automatic control community.

To find the optimal PID controller, a multi-objective optimization problem with constraints has to be solved [4,12,13,14]. Although there is no clear methodology to define the cost function to optimize [12] (usually, to minimize), most of them rely on the tracking error variation (both on transient and steady-state zones in signal dynamics). For example, root mean square error (RMSE) and integral RMSE are basic minimization criteria, from which different variants have been derived, such as integral absolute error, integral square error, integral time square error [1], integral time weighted square error and integral time absolute error [3,6]. By default, minimizing these criteria guarantees the trade-off between of a set of desired but conflicting control performance indexes (e.g., overshoot, settling time, rise time, decay ratio, steady-state error) [3,15,16]. Moreover, in real applications, since the process outputs are affected by environmental and/or measure disturbances (exogenous noises) that are stochastic in nature, the performance indexes cannot be computed by means of closed form expressions. In this case, building the multi-objective cost function is not an easy task. A compromise between opposite desired performances has to be found. Usually, a weighted sum criterion, with some experimentally set weights, is defined for this aim [12].

Two main classes of optimization techniques are mostly referred to in the literature for optimal PID controller design. The first one is represented by traditional (exact) optimization methods, being applied in context of linear/nonlinear control theory. The second one is concerned with metaheuristics.

Interestingly, some methods usually employing state space representations can be adapted to find optimal PID controllers. For example, one can mention the methods to synthesize Linear Quadratic Regulators of H-infinity controllers (with the help of the Linear Matrix Inequality optimization procedure) [11,17]. This approach has been successfully employed to design multivariable and multiloop PID controllers implemented either in continuous time [10,11] or in discrete time [18]. Nevertheless, this approach involves an elaborated mathematical foundation and complex algorithms [12]. Additionally, the proposed solutions seem to be insufficiently tested for noisy output signals.

The second class deals with non-traditional optimization techniques, also referred to as metaheuristics, mainly based on Artificial Intelligence [14], which led to evolutionary computing. In the last two decades, metaheuristics has gained increasing attention in various research domains, including control engineering, due to its capacity to find optimal solutions (often close to global optima) in multidimensional search spaces [12], with reasonable computational burden and a good trade-off between multiple performance objectives [5,14]. Many metaheuristics use the swarm intelligence revealed by the natural (individual and social) behavior of some living entities (insects, animals) [9]. A wide range of metaheuristics has been used in optimal PID design, as reported in the literature [5]. For example, the following metaheuristics were preferred: the Genetic Algorithms [19] in [6,8,16,17,20]; Particle Swarm Optimization [19] in [8,13,21]; Ant Colony Optimization [19] in [6,9]; Artificial Bee Colony [19] in [3,13,21,22]; Differential Evolution [23] in [8]; Firefly Algorithm [19] in [4,15,24]; Cuckoo Search [25] in [2,24]; Bat Algorithm Optimization [19] in [13,22]; Symbiotic Organisms Search [26] in [9]; and Grey-Wolf Optimization [27] in [3]. In some of these references, two or more metaheuristics were combined.

Intensive research has been reported in the literature for tuning the parameters of PID or Fractional Order PID [1] controllers using metaheuristics in the genuine or improved versions. Additionally, numerous comparative studies highlight the superiority of a metaheuristic over the others [6,8,21,22,24] or over empirical tuning methods [2,9,16,21].

Most contributions deal with continuous time SISO PID controllers [6,9,12,13,15,21,22,24]. To the best of our knowledge, few contributions are concerned with designing optimal PID controllers for the real-world MIMO processes [4,9,20]. In this case, the plant constructive characteristics, the operating conditions and, moreover, the noises corrupting the useful signals impose significant restrictions. Moreover, real-world processes often fail to verify some sufficient hypotheses usually employed to prove robustness of designed controllers. This involves the necessity to test robustness for each controller. Generally, the control strategy either is of the multi-loop PID type (decentralized) [4] or is based on decoupled SISO loops [3].

In this paper, two multivariable PID-based control strategies for a laboratory installation (ASTANK2 [28]) were designed and implemented starting from the identification models proposed in [29]. This resulted in two control strategies for each type of identification model. Based on the modern optimization context, to find the optimal PID parameters, a combination of two metaheuristics was proposed. The first is the Advanced Hill Climbing Algorithm (AHCA) [29], aiming to widely explore the search space and to find a sufficiently accurate optimal point. The second is an improved version of the Firefly Algorithm (FA) [19,30], employed to efficiently refine the exploitation around the previously found optimum point and to improve accuracy.

The article is structured as follows: Within the next section, the context of research development is described at length. Additionally, the optimization problem is formulated, and the methods to find optimal solutions, based on AHCA and improved FA metaheuristics, are proposed. Section 3 is entirely devoted to simulation results analysis. The investigations are directed towards three main directions: metaheuristics behavior, PID controllers performance for step references (such as bounded tracking error of step references, commands smoothness, short transient response, reduced overshoot, disturbance rejection) and robustness. The performance and robustness are comparatively analyzed. Four types of robustness are considered: to sudden change of references, to stochastic and constant disturbances, as well as to model uncertainties. To this aim, automatic procedures to define and detect fairly accurate robustness margins are introduced. The supplementary explanations are given in three appendices. Some concluding remarks, together with future developments are included into the last section of article.

2. Materials and Methods

2.1. Main Framework

The present study aims to design and implement a PID-based control strategy of a two-tank system, starting from the identification models determined in [29]. Actually, this research was already mentioned as possible future development in the end of that article.

The two-tank system, namely ASTANK2, is a laboratory installation, as illustrated in the functional diagram in Figure 1a.

The plant is composed of two main tanks (#1 and #2) located above an auxiliary storage reservoir, a main inverted-driven pump P (acting in the [0, 10] V range), two auxiliary on-off pumps (P1, P2) and a flexible pipeline network handled by seven manual taps (R11, R12, R13, R21, R22, R23, R). By acting the taps, different configurations of the installation can be set, such as serial/cascaded tank connection or independent tanks. The water is continuously recycled from the storage reservoir by means of the main pump P through the main vertical pipe that forks at the upper side and feeds the upper tanks. The supply flows of both tanks can be controlled by means of the two electro-valves (e-v) FC1 and FC2 (the voltage of which can also vary in the [0, 10] V range). The water can be evacuated from the main tanks either by a communication horizontal pipeline (in cases of serial configuration) or directly in the storage tank through the emptying pipes (operated by the manual taps R11, R21). The water level data are acquired by means of transducers LT1 and LT2. Additionally, the feeding water flows can be monitored thanks to FT1 and FT2 transducers. A detailed description of ASTANK2 can be found in [28].

As explained in [29], the two-tank system can be represented as a nonlinear MIMO system with $2\times 2$ I/O channels. The inputs are the voltages $u_1$ and $u_2$ of the two e-v. The outputs are the water levels in the upper tanks, $h_1$ and $h_2$, respectively. Note that the main pump voltage $U$ is kept constant to protect the pump against sudden variations or shocks.

Despite the nonlinear dynamics, in [29], a linear MIMO Box-Jenkins (BJ) identification model (idmo) of ASTANK2 were determined and implemented for simulations. In order to identify the MIMO system, two Multi-Input Single-Output (MISO) systems corresponding to each output channel $y_1$ and $y_2$, namely MISO₁ and MISO₂, respectively, were identified, starting from the experimental data. A comprehensive block diagram of the global model identification is illustrated in Figure 1b.

Each MISO component model includes three additive filters. Each input ($u_1$, $u_2$) stimulates two useful linear filters of $\mathrm{B}/\mathrm{A}$ type. A noise linear filter of the $\mathrm{C}/\mathrm{D}$ type modifies each output after processing two uncorrelated Gaussian white noises ($e_1$, $e_2$). For each output channel $y_j$ ($j\in \{1,2\}$), the MISO model relies on polynomials ${\mathrm{A}}_j$, ${\mathrm{B}}_{j1}$, ${\mathrm{B}}_{j2}$, ${\mathrm{C}}_j$ and ${\mathrm{D}}_j$, of degrees $n{a}_j$, $n{b}_{j1}$, $n{b}_{j2}$, $n{c}_j$ and $n{d}_j$, respectively, with intrinsic delays $n{k}_{j1}$, $n{k}_{j2}$, defined as:

$$\left\{\begin{array}{l}{\mathrm{A}}_j\left({\mathrm{q}}^{-1}\right)=1+a_{1,j}{\mathrm{q}}^{-1}+\cdots +a_{n{a}_j,j}{\mathrm{q}}^{-n{a}_j}\\ {\mathrm{B}}_{j1}\left({\mathrm{q}}^{-1}\right)={\mathrm{q}}^{1-n{k}_{j1}}\left(b_{1,j1}{\mathrm{q}}^{-1}+b_{2,j1}{\mathrm{q}}^{-2}\cdots +b_{n{b}_{j1},j1}{\mathrm{q}}^{-n{b}_{j1}}\right)\\ {\mathrm{B}}_{j2}\left({\mathrm{q}}^{-1}\right)={\mathrm{q}}^{1-n{k}_{j2}}\left(b_{1,j2}{\mathrm{q}}^{-1}+b_{2,j2}{\mathrm{q}}^{-2}\cdots +b_{n{b}_{j2},j2}{\mathrm{q}}^{-n{b}_{j2}}\right)\\ {\mathrm{C}}_j\left({\mathrm{q}}^{-1}\right)=1+c_{1,j}{\mathrm{q}}^{-1}+\cdots +c_{n{c}_j,j}{\mathrm{q}}^{-n{c}_j}\\ {\mathrm{D}}_j\left({\mathrm{q}}^{-1}\right)=1+d_{1,j}{\mathrm{q}}^{-1}+\cdots +d_{n{d}_j,j}{\mathrm{q}}^{-n{d}_j}\end{array}\right..$$

(1)

Two types of experimental data were employed to identify the BJ models: raw (directly acquired from the sensors) and smoothed. The smoothed data were built from the raw data by using a least squares optimal stochastic cubic spline model [29]. In both cases, the optimal structure of the BJ models was obtained by means of a metaheuristic that maximizes the Signal-to-Noise Ratio (SNR).

Table 1. Structural indices of optimal MISO BJ models and estimated noise variances.

Channel	Output Data Type	Optimal $\left(\mathit{n}\mathit{a},\mathit{n}{\mathit{b}}_1,\mathit{n}{\mathit{b}}_2,\mathit{n}\mathit{c},\mathit{n}\mathit{d}\right)$	Noise Variance $\left({\mathit{\lambda }}^2\right)$
1	Raw	(30,22,16,65,0)	0.0067
1	Spline Smoothed	(3,21,47,89,0)	0.006
2	Raw	(27,7,26,95,0)	0.006
2	Spline Smoothed	(2,2,46,95,0)	0.0063

illustrates the structural indices of each optimal MISO subsystem for both raw and smoothed data (whilst $n{k}_{j1}=n{k}_{j1}=1$), as well as the variances of both Gaussian white noises.

The simulation results showed that when using smoothed data, the useful models have significantly fewer poles than when using raw data (see the first structural index, $na$). Such models are expected to yield better performance in closed loop than the models identified from raw data.

In addition to the idmo of ASTANK2, a static input/output (I/O) map was built experimentally. This map includes quadruples $\left(u_{1,\mathrm{s}},u_{2,\mathrm{s}},y_{1,\mathrm{s}},y_{2,\mathrm{s}}\right)$ with the following meaning: if the tanks #1 and #2 need to be filled with water at constant levels $y_{1,\mathrm{s}}$ and $y_{2,\mathrm{s}}$, respectively, then, in the absence of perturbations, the corresponding constant voltages of e-v FC1 and FC2 are $u_{1,\mathrm{s}}$ and $u_{2,\mathrm{s}}$, respectively. This map is extremely useful in setting nominal functioning points of installation, around of which the dynamics of nonlinear installation are almost linear. Hence, if some reference water levels $y_{1,\mathrm{ref}}$ and $y_{2,\mathrm{ref}}$ are prescribed, then, from the static map, one can obtain the necessary voltages $u_{1,\mathrm{ref}}$ and $u_{2,\mathrm{ref}}$ to be applied on each e-v. If $y_{1,\mathrm{ref}}$ and/or $y_{2,\mathrm{ref}}$ do not correspond to the static points in the map, then the voltages $u_{1,\mathrm{ref}}$ and $u_{2,\mathrm{ref}}$ are obtained by linear interpolation. The static map includes 49 combinations of static water levels (7 for each tank), between 7 cm and 11 cm. This range was preferred to avoid water turbulences in tanks, which may provide wrong acquired data. Moreover, the linear models above were identified for the same level range.

The static map also revealed some constructive limitations of the ASTANK2 plant, to be considered as constraints in the design of any controller. The most important of them are as follows: (a) the e-v voltages are saturated at 2 V and 10 V (below 2 V, the e-v does not react any more); (b) the water levels in any tank can vary between 0 cm and 15 cm only; (c) to avoid hard saturation of voltages on e-v, the absolute difference between level references, $\left|y_{1,\mathrm{ref}}-y_{2,\mathrm{ref}}\right|$, is bounded by 2 cm; (d) to protect e-v against premature failure, any voltage cannot jump from the current value to the next one with more than $\pm 3$ V.

2.2. Selected PID Controllers and Closed Loop Configurations

The proposed control strategy relies on PID type controllers. Since the ASTANK2 plant is a multivariable system, two PID-based control configurations were designed and investigated (among many possible others) in order to meet certain performances, such as bounded tracking error of step references, smooth commands, short transient response, reduced overshoot, disturbance rejection, robustness against reference shocks, external perturbations and model uncertainties. An extended analysis of performance indicators is made in Section 3, devoted to simulation results. Note, however, that the optimal PID controllers were synthesized by means of a cost function, referred to as fitness, defined in the next subsection. The fitness does not consider all possible performance indicators, especially because many of them are defined for systems with smooth I/O signals. Alternatively, ASTANK2 installation works with noisy signals, as shown later.

The two selected PID are described by the schemes in Figure 2 and Figure 3. They are depicted for SISO continuous time systems to allow for easy understanding. In both configurations, $\mathrm{H}$ is the transfer function of system to regulate, $u$ is the input (the command), $y$ is the output, $y_r$ is the reference signal, $\mathsf{\epsilon }\equiv y_r-y$ is the tracking error, whereas $K_p>0$ (the proportionality constant), $T_i>0$ (the integral constant), $T_d\ge 0$ (the derivative constant) and $N_d>0$ (the derivative divisor) are parameters of the PID controller.

The derivative component is filtered by a first order system to yield physical realizability. The derivative divisor $N_d$ usually is an integer at least equal to 10. Thus, the time constant of filter ($T_d/N_d$) is small enough (i.e., its bandwidth is large enough) to minimally distort the derivative component of the command. Here and hereafter, $N_d$ was set to 10.

The fist PID structure, as seen in Figure 2, is a regular (standard) one, being abbreviated by PID-R. The characteristic transfer function of the PID-R controller, from tracking error $\mathsf{\epsilon }$ to command $u$, is:

$${\mathrm{H}}_{\mathrm{PID}-\mathrm{R}}(s)=K_p\left[1+\frac{1}{T_{i}s}+\frac{T_{d}s}{\frac{T_d}{N_d}s+1}\right]=K_p\frac{T_i{T}_d\left(N_d+1\right)s^2+\left(T_d+N_d{T}_i\right)s+N_d}{T_d{T}_i{s}^2+N_d{T}_{i}s}.$$

(2)

The second structure, illustrated in Figure 3, is non-standard and was considered for its main feature; the derivative component of the command is furthermore filtered by the other two components (P, I), which results in quite smooth commands (with as few spikes as possible). The acronym for this structure is PID-N, and the corresponding transfer function is:

$${\mathrm{H}}_{\mathrm{PID}-\mathrm{N}}(s)=K_p\frac{T_{i}s+1}{T_{i}s}\left[\begin{array}{c}1-\frac{T_d\left(N_d+1\right)s+N_d}{T_{d}s+N_d}\end{array}\right].$$

(3)

Unlike PID-R, the transfer function of PID-N is of the MISO type, with two inputs: the reference signal $y_r$ and the system output $y$. This allows for working with two regulation blocks of the SISO type: one on feedback, ${\mathrm{H}}_R$ (the ”pure” regulator), and another one in front of the system, ${\mathrm{H}}_{\mathbf{S}}$ (the ”sensitive” controller), providing the command $u$:

$${\mathrm{H}}_{\mathrm{PID}-\mathrm{N}}(s)={\mathrm{H}}_{\mathbf{S}}(s)\left[\begin{array}{c}1-{\mathrm{H}}_R(s)\end{array}\right],\mathrm{with}:\left\{\begin{array}{l}{\mathrm{H}}_{\mathbf{S}}(s)=K_p\frac{T_{i}s+1}{T_{i}s}\\ {\mathrm{H}}_R(s)=\frac{T_d\left(N_d+1\right)s+N_d}{T_{d}s+N_d}\end{array}\right..$$

(4)

Since the idmo of the ASTANK2 plant is discrete, it is suitable to discretize the PID regulators. A fairly accurate discretizing method is the bilinear Tustin method, in which $s\simeq \frac{2}{T_s}\frac{1-z^{-1}}{1+z^{-1}}$, for some sampling period $T_s>0$. After discretization, the transfer functions (2) and (4) become:

$${\mathrm{H}}_{\mathrm{PID}-\mathrm{R}}(s)\rightleftarrows {\mathrm{H}}_{\mathrm{PID}-\mathrm{R}}(z)=\frac{\mathrm{Q}(z)}{\mathrm{P}(z)}\mathrm{and}{\mathrm{H}}_{\mathrm{PID}-\mathrm{N}}(s)\rightleftarrows {\mathrm{H}}_{\mathrm{PID}-\mathrm{N}}(z)=\frac{{\mathrm{Q}}_{\mathbf{S}}(z)}{{\mathrm{P}}_{\mathbf{S}}(z)}\left[\begin{array}{c}1-\frac{{\mathrm{Q}}_R(z)}{{\mathrm{P}}_R(z)}\end{array}\right],$$

(5)

where $\mathrm{Q}$, $\mathrm{P}$, ${\mathrm{Q}}_{\mathbf{S}}$, ${\mathrm{P}}_{\mathbf{S}}$, ${\mathrm{Q}}_R$, ${\mathrm{P}}_R$ are polynomials of the second degree at most. Their coefficients are expressed as follows, depending on the controller type:

• P controllers ($T_i=\infty$, $T_d=0$):

  $$\left\{\begin{array}{l}{\mathrm{Q}}_{\mathbf{S}}(z)=\mathrm{Q}(z)=K_p\\ {\mathrm{P}}_{\mathbf{S}}(z)=\mathrm{P}(z)=1\end{array}\right.;\left\{\begin{array}{l}{\mathrm{Q}}_R(z)=1\\ {\mathrm{P}}_R(z)=1\end{array}\right..$$

  (6)
• PI controllers ($T_d=0$):

  $$\left\{\begin{array}{l}\mathrm{Q}(z)={\mathrm{Q}}_{\mathbf{S}}(z)=K_p\left[\left(T_s+2T_i\right)z+\left(T_s-2T_i\right)\right]\\ \mathrm{P}(z)={\mathrm{P}}_{\mathbf{S}}(z)=2T_i(z-1)\end{array}\right.;\left\{\begin{array}{l}{\mathrm{Q}}_R(z)=1\\ {\mathrm{P}}_R(z)=1\end{array}\right..$$

  (7)
• PD controllers ($T_i=\infty$):

  $$\left\{\begin{array}{l}\mathrm{Q}(z)=K_p{\mathrm{Q}}_R(z)=K_p\left\{\left[N_d{T}_s+2\left(N_d+1\right)T_d\right]z+\left[N_d{T}_s-2\left(N_d+1\right)T_d\right]\right\}\\ \mathrm{P}(z)={\mathrm{P}}_R(z)=\left(N_d{T}_s+2T_d\right)z+\left(N_d{T}_s-2T_d\right)\end{array}\right.;\left\{\begin{array}{l}{\mathrm{Q}}_{\mathbf{S}}(z)=K_p\\ {\mathrm{P}}_{\mathbf{S}}(z)=1\end{array}\right..$$

  (8)
• PID controllers (general case): $\left\{\begin{array}{l}\mathrm{Q}(z)=q_2{z}^2+q_{1}z+q_0\\ \mathrm{P}(z)=p_2{z}^2+p_{1}z+p_0\end{array}\right.$, with:

  $$\left\{\begin{array}{l}{q}_2=K_p\left[2\frac{\left(N_d+1\right)T_i{T}_d}{N_d{T}_s}+\frac{T_s}{2}+\frac{T_d}{N_d}+T_i\right]\\ q_1=K_p\left[T_s-4\frac{\left(N_d+1\right)T_i{T}_d}{N_d{T}_s}\right]\\ q_0=K_p\left[2\frac{\left(N_d+1\right)T_i{T}_d}{N_d{T}_s}+\frac{T_s}{2}-\frac{T_d}{N_d}-T_i\right]\end{array}\right.;\left\{\begin{array}{l}{p}_2=2\frac{T_i{T}_d}{N_d{T}_s}+T_i\\ p_1=-4\frac{T_i{T}_d}{N_d{T}_s}\\ p_0=2\frac{T_i{T}_d}{N_d{T}_s}-T_i\end{array}\right.;$$

  (9)

  $$\left\{\begin{array}{l}{\mathrm{Q}}_{\mathbf{S}}(z)=K_p\left[\left(T_s+2T_i\right)z+\left(T_s-2T_i\right)\right]\\ {\mathrm{P}}_{\mathbf{S}}(z)=2T_i(z-1)\\ {\mathrm{Q}}_R(z)=\left[N_d{T}_s+2\left(N_d+1\right)T_d\right]z+\left[N_d{T}_s-2\left(N_d+1\right)T_d\right]\\ {\mathrm{P}}_R(z)=\left(N_d{T}_s+2T_d\right)z+\left(N_d{T}_s-2T_d\right)\end{array}\right..$$

  (10)

For ASTANK2 installations, the selected closed loop configurations are displayed in Figure 4, for the two PID structures. There are 4 characteristic polynomials for the configuration in Figure 4a: ${\mathrm{Q}}_1$, ${\mathrm{P}}_1$, ${\mathrm{Q}}_2$ and ${\mathrm{P}}_2$. Additionally, in Figure 4b, 8 polynomials have to be set: ${\mathrm{Q}}_{\mathbf{S},1}$, ${\mathrm{P}}_{\mathbf{S},1}$, ${\mathrm{Q}}_{R,1}$, ${\mathrm{P}}_{R,1}$, ${\mathrm{Q}}_{\mathbf{S},2}$, ${\mathrm{P}}_{\mathbf{S},2}$, ${\mathrm{Q}}_{R,2}$ and ${\mathrm{P}}_{R,2}$. Their coefficients are determined by the 6 parameters of continuous time PID controllers, namely $K_{p,1}$, $T_{i,1}$, $T_{d,1}$ (for the first channel), and $K_{p,2}$, $T_{i,2}$, $T_{d,2}$ (for the second channel). To compute them, Equations (6)–(10) can be employed. Each PID controller is fed with the output and/or the tracking error of one channel, but provides a command that applies to both MISO blocks.

2.3. Optimization Problem Statement and Cost Function Definition

One seeks to find the optimal parameters of the two PID controllers by optimizing a cost function. There is an extensive discussion in the literature concerning the choice of such a cost function. Most approaches rely on the RMSE applied on the tracking error and the norm of command. Both have to be minimized. Usually, they are combined in a cost function with some weights in their original or modified versions. Basically, the cost function defined in this article, referred to as fitness, also employs the RMSE, but with some refinements adapted to ASTANK2 dynamics.

Take, for example, the output signal $y$ of Figure 5, tracking some reference $y_{\mathrm{ref}}$.

The global RMSE is then:

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{n=1}^N{\left|y[n]-y_{\mathrm{ref}}[n]\right|}^2}},$$

(11)

where $N\in {\mathbb{N}}^{\ast }$ is the number of acquired output samples. Evidently, RMSE is a cost function to minimize. Using definition (11) in the context of the ASTANK2 plant revealed several disadvantages. The most important comes from the fact the tracking error can have similar values on transient and steady-state zones, which can lead to output signals that oscillate too much around references, while RMSE takes relatively small values. Thus, separating the two zones is necessary. Since the output signal is corrupted by stochastic noises, no analytical formula can be employed to compute the transient time, $N_{\mathrm{t}}$. Therefore, $N_{\mathrm{t}}$ has to be estimated by using statistical tools.

Let $\mathsf{\epsilon }\equiv y-y_{\mathrm{ref}}$ be the tracking error. A key parameter in estimating $N_{\mathrm{t}}$ is the backward moving standard deviation (std) of tracking error:

$${\mathsf{\sigma }}_{\mathsf{\epsilon }}^{\mathrm{b}}[k]=\sqrt{\frac{1}{k+1}{\displaystyle \sum _{n=N-k}^N{\left|\mathsf{\epsilon }[n]-{\mathsf{\mu }}_{\mathsf{\epsilon }}^{\mathrm{b}}[k]\right|}^2}},\forall \hspace{0.17em}k\in \overline{1,N-1},$$

(12)

where ${\mathsf{\mu }}_{\mathsf{\epsilon }}^{\mathrm{b}}$ is the backward moving average:

$${\mathsf{\mu }}_{\mathsf{\epsilon }}^{\mathrm{b}}[k]=\frac{1}{k+1}{\displaystyle \sum _{n=N-k}^N\mathsf{\epsilon }[n]},\forall \hspace{0.17em}k\in \overline{1,N-1}.$$

(13)

For the signal in Figure 5, ${\mathsf{\sigma }}_{\mathsf{\epsilon }}^{\mathrm{b}}$ is displayed in Figure 6a. At the right side end, the std of final tracking error is computed, whereas the larger values obtained at the left side end correspond to the first values of tracking error.

Actually, the variation above was horizontally flipped, as, according to the computing formula, the first values of std are the ones that lie at the right side end. Normally, the transient time is given at the left entry of the flat variation in the graphic. However, finding this entry is rather difficult, because, as one can easily notice, any of the lags between 180 and 300 can play this role. Therefore, a penalty can be applied to the std, aiming to better isolate the transient time. Among several penalties known in the literature, the Akaike-Rissanen (A-R) criterion seems to be suitable to the variation above. The A-R penalized std is then:

$${\mathsf{\sigma }}_{\mathsf{\epsilon }}^{\mathrm{b},\mathrm{A}-\mathrm{C}}[k]=\log \left({\mathsf{\sigma }}_{\mathsf{\epsilon }}^{\mathrm{b}}[k]\right)+\frac{k}{2N}\log (N),\forall \hspace{0.17em}k\in \overline{1,N-1}.$$

(14)

In Figure 6b, the penalized A–R std is overlapped on normalized signal and reference variations (in range [0, 2]). One can notice how easy it now is to find the transient time (the minimum point of the penalized std) and how well this value corresponds to the subjective selection performed by visual inspection of signal variation.

In the framework of the ASTANK2 plant, the two output signals $y_1$ and $y_2$ are tracking the references $y_{1,\mathrm{ref}}$ and $y_{2,\mathrm{ref}}$, with errors ${\mathsf{\epsilon }}_1^y$ and ${\mathsf{\epsilon }}_2^y$, respectively. To define the fitness, the two zones of each signal dynamics are processed separately. Let $N_{\mathrm{t},1}$ and $N_{\mathrm{t},2}$ be the two transient times detected on output signals by using the procedure above. For the transient zone, the fitness can be defined as follows:

$${\mathsf{\phi }}_{\mathrm{t},j}^y\equiv \frac{100}{1+\mathsf{\sigma }\sqrt{\frac{1}{N_{\mathrm{t},j}}{\displaystyle \sum _{n=1}^{N_{\mathrm{t},j}}{\left|{\mathsf{\epsilon }}_j^y[n]\right|}^2}}}[\%],\forall \hspace{0.17em}j\in \{1,2\},$$

(15)

where $\mathsf{\sigma }>0$ is a factor that allows for scaling of the fitness. For example, the function $100/(1+x)$ (on which the definition (15) relies) takes the value of 58.82% in $x=0.7$; by using the scaling factor $\mathsf{\sigma }=0.5$, i.e., $\mathsf{\sigma }\hspace{0.17em}x$ instead of $x$, the value can be increased to 74.07%. Assuming both output signals have the same length, $N$, for the steady-state zone, the fitness is:

$${\mathsf{\phi }}_{\mathrm{s},j}^y\equiv \frac{100}{1+\mathsf{\sigma }\sqrt{\frac{1}{N-N_{\mathrm{t},j}}{\displaystyle \sum _{n=N_{\mathrm{t},j}+1}^N{\left|{\mathsf{\epsilon }}_j^y[n]\right|}^2}}}[\%],\forall \hspace{0.17em}j\in \{1,2\}.$$

(16)

Now, it is natural to consider that each fitness above has a weight proportional to the duration of its zone, which leads to the fitness of each output channel:

$${\mathsf{\varphi }}_j^y\equiv \frac{N_{\mathrm{t},j}{\mathsf{\phi }}_{\mathrm{t},j}^y+\left(N-N_{\mathrm{t},j}\right){\mathsf{\phi }}_{s,j}^y}{N},\forall \hspace{0.17em}j\in \{1,2\}.$$

(17)

In the case of commands $u_1$ and $u_2$, the references $u_{1,\mathrm{ref}}$ and $u_{2,\mathrm{ref}}$ are computed from the static map, depending on references $y_{1,\mathrm{ref}}$ and $y_{2,\mathrm{ref}}$, as explained in the end of Section 2.2. Two zones in command variation are considered as well, with transient times $M_{\mathrm{t},1}$ and $M_{\mathrm{t},2}$. Unlike for outputs, the two commands should be analyzed as a whole, since they are simultaneously applied to both MISO models that produce the outputs; therefore, on the transient zones, a unique fitness is defined by means of tracking errors ${\mathsf{\epsilon }}_1^u$ and ${\mathsf{\epsilon }}_2^u$:

$${\mathsf{\phi }}_{\mathrm{t},12}^u\equiv \frac{100}{1+\mathsf{\sigma }\sqrt{\frac{1}{M_{\mathrm{t},1}+M_{\mathrm{t},2}}\left({\displaystyle \sum _{n=1}^{M_{\mathrm{t},1}}{\left|{\mathsf{\epsilon }}_1^u[n]\right|}^2}+{\displaystyle \sum _{n=1}^{M_{\mathrm{t},2}}{\left|{\mathsf{\epsilon }}_2^u[n]\right|}^2}\right)}}[\%].$$

(18)

Similarly, the fitness on the steady-state zone is:

$${\mathsf{\phi }}_{\mathrm{s},12}^u\equiv \frac{100}{1+\mathsf{\sigma }\sqrt{\frac{1}{2N-M_{\mathrm{t},1}-M_{\mathrm{t},2}}\left({\displaystyle \sum _{n=M_{\mathrm{t},1}+1}^N{\left|{\mathsf{\epsilon }}_1^u[n]\right|}^2}+{\displaystyle \sum _{n=M_{\mathrm{t},2}+1}^N{\left|{\mathsf{\epsilon }}_2^u[n]\right|}^2}\right)}}[\%].$$

(19)

Finally, the fitness of both input channels is a weighted average of (18) and (19):

$${\mathsf{\varphi }}_{12}^u\equiv \frac{\left(M_{\mathrm{t},1}+M_{\mathrm{t},2}\right){\mathsf{\phi }}_{\mathrm{t},12}^u+\left(2N-M_{\mathrm{t},1}-M_{\mathrm{t},2}\right){\mathsf{\phi }}_{\mathrm{s},12}^u}{2N}.$$

(20)

The fitness definitions (17) and (20) can be combined into a weighted average, to obtain the overall fitness of PID-ASTANK2 scheme:

$$\Phi \equiv \mathsf{\alpha }\left({\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y\right)+\mathsf{\beta }\hspace{0.17em}{\mathsf{\varphi }}_{1,2}^u[\%],$$

(21)

where $2\mathsf{\alpha }+\mathsf{\beta }=1$. For example, $\mathsf{\alpha }=0.35$ and $\mathsf{\beta }=0.3$. The fitness of output signals has the same weight, as no channel is more important than the other. Obviously, the fitness is a cost function to maximize.

Definition (21) has an important drawback. Theoretically, a certain value of fitness, ${\Phi }_0$, can be obtained by an infinity of triplets $\left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y,{\mathsf{\varphi }}_{1,2}^u\right\}$. Thus, one of the fitness functions from the triplet can be small, whilst, in compensation, the other two are large. Usually, this happens when ${\mathsf{\varphi }}_{1,2}^u$ is small (the commands having big oscillations) and ${\mathsf{\varphi }}_1^y$ and ${\mathsf{\varphi }}_2^y$ are large (tracking of references being very fast, with small overshoot or even without overshoot). This means the e-v are worked hard in order to fasten the reference tracking. Or, such a situation is to be avoided. It is more suitable to have smoother commands, even though the references are tracked slower and/or with some overshoot. To deal with this ambiguity, the characteristics of a closed loop system can be taken into consideration. For an ASTANK2 installation, simulations have shown that good fitness values are above 75% (for definition (21) and the scaling factor $\mathsf{\sigma }=0.5$). In this case, the difference between ${\mathsf{\varphi }}_1^y$ and ${\mathsf{\varphi }}_2^y$ is less than 10% ($\left|{\mathsf{\varphi }}_1^y-{\mathsf{\varphi }}_2^y\right|<0.1$). Moreover, commands fitness exhibit two behaviors: (a) if ${\mathsf{\varphi }}_{12}^u<\frac{{\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y}{2}$, then the overall fitness can take big values, although the commands are oscillating too much (increasing the risk of e-v premature wearing out or failure); (b) if, on the contrary, ${\mathsf{\varphi }}_{12}^u\ge \frac{{\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y}{2}$ (or, even better, if ${\mathsf{\varphi }}_{12}^u\ge \max \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}$), then the commands are quite smooth, and the overall fitness can take large values. Obviously, the last case is suitable for finding optimal controllers. In order to drive the search towards smooth commands and large overall fitness values, definition (21) can be refined:

$$\Phi \equiv \left\{\begin{array}{ll}\mathsf{\beta }\frac{{\mathsf{\varphi }}_{12}^u+\max \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}}{2}+2\mathsf{\alpha }\min \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}\hfill & ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\min \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}<75\hfill \\ \mathsf{\beta }\frac{{\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y}{2}+2\mathsf{\alpha }\hspace{0.17em}{\mathsf{\varphi }}_{12}^u\hfill & ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\min \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}\ge 75\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\&\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathsf{\varphi }}_{12}^u<\frac{{\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y}{2}\hfill \\ \mathsf{\alpha }\left({\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y\right)+\mathsf{\beta }\hspace{0.17em}{\mathsf{\varphi }}_{12}^u\hfill & ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\min \left\{{\mathsf{\varphi }}_1^y,{\mathsf{\varphi }}_2^y\right\}\ge 75\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\&\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathsf{\varphi }}_{12}^u\ge \frac{{\mathsf{\varphi }}_1^y+{\mathsf{\varphi }}_2^y}{2}\hfill \end{array}\right.,$$

(22)

where $\mathsf{\alpha }\in (0.25,0.5)$. Practically, in definition (22), the former definition (21) is only employed to find optimal regulators with smooth commands—see the third case. For the other two cases, the term that produce fitness ambiguity is now multiplied by large weight (of at least 0.5), which results in smaller overall fitness values.

The fitness depends on the 6 parameters of continuous time PID controllers. Obviously, it is not possible to visually observe the characteristics of fitness in this case. Nevertheless, 6 fitness surfaces can be drawn, by varying couples of parameters, while keeping constant the other four. Such surfaces are depicted in Figure 7.

The stochastic nature of the I/O measured data induces irregular variations, with many local extremes of fitness. A mutual characteristic of all surfaces above is the alternance between smooth slopes and fractal heights or ridges, where the maximum point lies. Solving the maximization problem becomes a difficult task that cannot be accomplished by using exact methods [31]. Instead, a suitable alternative to this problem is adopting evolutionary computing techniques, referred to as metaheuristics [19]. In this article, two metaheuristics are employed: one to widely explore the search space and another one to exploit some narrow zones where the maximum seems to exist. Figure 8 suggests the search strategy. For exploration, the AHCA from [29] was employed. Its main task is to find an optimum with fair accuracy, while exploring large zones of search space. Since in case of the 6 optimal PID parameters the upper bound is unknown and the fitness exhibits smooth slopes where the gradient is not oscillating, AHCA fits very well. Once the climbers have found some optima, a vicinity around them can be isolated. Now the search can be performed within a bounded zone. To exploit this zone, the FA was selected (and improved), due to its simplicity and effectiveness. Whilst the climbers are strongly dependent on the local gradient of fitness, the fireflies can easily jump from one peak to another, regardless of the gradient (which, anyway, may oscillate too much—see the magnification on the left side of Figure 8).

It has to be outlined that both the exploration and the exploitation phases are performed by varying all 6 parameters (and not only couples of them, as Figure 7 might suggest). The improved FA (IFA) is described next.

2.4. Improved Firefly Algorithm

The FA [19,30] was inspired by the principle of mutual attraction between insects and the phenomenon of natural attenuation of light with distance. The fireflies are small flying beetles, usually luminescent at twilight, that produce a cold light to attract mates or prey. They are grouped into swarms and, apparently, fly at random. In its flight, a firefly is tempted to fly towards other shinier fireflies in the swarm, being attracted by the light signals the latter are emitting. In this way, the shinier fireflies become targets for the current one. The strongest target is the shiniest firefly, which attracts the entire swarm. Thus, the swarm is more and more grouped around the firefly emitting the most intense light. Nevertheless, the shiniest firefly flies as well and their light can fade. In this case, another firefly becomes the shiniest. Even though the swarm was grouped around a firefly, it can move away towards another shinier firefly. In this way, there is a good balance between the convergence and the diversity of firefly swarms.

Figure 9 illustrates a swarm of fireflies and its evolution from current to next positions. The current firefly is colored in blue, and the shiniest firefly is in red, whereas the other fireflies are in orange (the shinier ones) or in yellow (the less shining ones).

The perceived light attractiveness depends on the brilliance of the casting beetle. The brighter a firefly, the more fireflies it attracts. The current firefly is motivated to fly towards the shiniest one along the arrow in red (see the left side of Figure 9). Nevertheless, while tracking the target, our firefly can perceive lights from other fireflies, especially from the shinier ones, which constitute attraction points as well. Hence, the other fireflies in the swarm, especially the ones in orange, could disturb the flight along the direction colored in orange as well. By composing the two directions (in red and in orange), the final direction, in blue, is taken by the current firefly. Meanwhile, the swarm moves to new positions and changes light intensities, as suggested in the right side of Figure 9. It is possible that the shiniest firefly decreases its light intensity, whilst another firefly becomes the shiniest. All fireflies can change not only in position, but in light intensity as well. The swarm permanently tries to group around the shiniest firefly, but this becomes possible only if no other firefly becomes shinier. Otherwise, the swarm can oscillate between several positions, until a best position (of maximum light intensity) is found, eventually.

This behavior can be simulated by means of a numerical algorithm based on the following recursive equation expressing how the position of an artificial firefly changes:

$$x_p^{k+1}=x_p^k+{\displaystyle \sum _{q\in {\mathbf{S}}_p}{\mathsf{\beta }}_{p,q}^k\left(x_q^k-x_p^k\right)}+{\mathsf{\eta }}^k{\mathsf{\xi }}^k\left(\mathbf{S}\right),\forall \hspace{0.17em}p\in \overline{1,P},\forall \hspace{0.17em}k\in \mathbb{N}.$$

(23)

The notations in Equation (23) are explained next. Assume the swarm includes $P\in {\mathbb{N}}^{\ast }$ fireflies. Each firefly $p\in \overline{1,P}$ takes a position $x_p^k$ at instant $k\in \mathbb{N}$, into the search space $\mathbf{S}\subseteq {ℝ}^{nx}$, where $nx\in {\mathbb{N}}^{\ast }$ is the length of position vector. (For example, in the case of PID-ASTANK2 configurations, $nx=6$.) The search space can be specified as a cartesian product of compact intervals:

$$\mathbf{S}=\left[x_{\min ,1},x_{\max ,1}\right]\times \left[x_{\min ,2},x_{\max ,2}\right]\times \cdots \times \left[x_{\min ,nx},x_{\max ,nx}\right]={\displaystyle \prod _{i=1}^{nx}\left[x_{\min ,i},x_{\max ,i}\right]}.$$

(24)

To compute the next position, the firefly will take ($x_p^{k+1}$), and the current position is affected by two disturbances: one produced by the shinier fireflies than the current one (including the shiniest one), referred to as diversity disturbance, and another one generated by the entire swarm—the convergence disturbance. The sum term in Equation (23) stands for the diversity disturbance, where all the shinier fireflies than the current one are grouped into the set ${\mathbf{S}}_p$. This disturbance actually contributes to swarm diversity and is expressed as a linear combination of attraction directions from the current firefly to its shinier mates, ${\left\{x_q^k-x_p^k\right\}}_{q\in {\mathbf{S}}_p}$. Thus, the current firefly is not exclusively attracted by the shiniest firefly, but by the other shinier fireflies as well. This phenomenon keeps the diversity of fireflies’ swarms to some extent. The corresponding coefficients, ${\left\{{\mathsf{\beta }}_{p,q}^k\right\}}_{q\in {\mathbf{S}}_p}$, are light intensities perceived by the current firefly from the shinier ones. According to Biology, each perceived intensity can be computed as follows:

$${\mathsf{\beta }}_{p,q}^k=I_0\exp \left[-\frac{‖x_p^k-x_q^k‖}{\mathbf{D}(\mathbf{S})}\right],\forall \hspace{0.17em}p\in \overline{1,P},\forall \hspace{0.17em}q\in {\mathbf{S}}_p,\forall \hspace{0.17em}k\in \mathbb{N},$$

(25)

where $I_0$ is a biological constant standing for the maximum light intensity a firefly can produce (usually set to unit) and $\mathbf{D}(\mathbf{S})$ is the diameter of search space, defined by:

$$\mathbf{D}(\mathbf{S})=\underset{x,y\in \mathbf{S}}{\sup }\left\{‖x-y‖\right\}=\sqrt{{\displaystyle \sum _{i=1}^{nx}{\left(x_{\min ,i}-x_{\max ,i}\right)}^2}}.$$

(26)

Definition (25) shows that the perceived light intensity is higher if the fireflies are closer to each other. However, the current firefly also can be strongly attracted by a distant firefly radiating an intense light. This case is not taken into consideration in definition (25). Therefore, a first improvement consists of including the radiated light intensity into the perceived intensity formula. Let ${\Phi }_p^k$ be the light intensity radiated by firefly $p\in \overline{1,P}$ at instant $k\in \mathbb{N}$. Then, the intensity of light the other fireflies perceive should increase as ${\Phi }_p^k$ increases and the distance from the current firefly decreases. A second improvement concerns the diameter. Normally, the perceived intensity should adapt itself to the space the entire swarm of fireflies is taking at current instant, denoted by ${\mathbf{S}}^{\hspace{0.17em}k}$. Obviously, ${\mathbf{S}}^{\hspace{0.17em}k}\subseteq \mathbf{S}$ and then:

$${\mathbf{S}}^{\hspace{0.17em}k}={\displaystyle \prod _{i=1}^{nx}\left[\underset{p\in \overline{1,P}}{\min }\left\{x_{p,i}^k\right\},\underset{p\in \overline{1,P}}{\max }\left\{x_{p,i}^k\right\}\right]}={\displaystyle \prod _{i=1}^{nx}\left[x_{\min ,i}^k,x_{\max ,i}^k\right]};\mathbf{D}\left({\mathbf{S}}^{\hspace{0.17em}k}\right)=\sqrt{{\displaystyle \sum _{i=1}^{nx}{\left(x_{\max ,i}^k-x_{\min ,i}^k\right)}^2}}.$$

(27)

With the two improvements above, definition (25) can be replaced by:

$${\mathsf{\beta }}_{p,q}^k=I_0\exp \left[-\frac{‖x_p^k-x_q^k‖}{\mathbf{D}\left({\mathbf{S}}^{\hspace{0.17em}k}\right)}\frac{{\Phi }_p^k}{{\Phi }_q^k}\right],\forall \hspace{0.17em}p\in \overline{1,P},\forall \hspace{0.17em}q\in {\mathbf{S}}_p,\forall \hspace{0.17em}k\in \mathbb{N},$$

(28)

where, evidently, ${\Phi }_q^k\ge {\Phi }_p^k$, $\forall \hspace{0.17em}q\in {\mathbf{S}}_p$.

The last term of Equation (23) expresses the second disturbance, generated by the entire swarm at each instant $k\in \mathbb{N}$. The disturbance intensity is given by the scalar factor ${\mathsf{\eta }}^k$, uniformly selected at random in the interval $[-1,+1]$. Its sign encodes whether the disturbance is direct (+1) or opposite (–1). The other factor, ${\mathsf{\xi }}^k\left(\mathbf{S}\right)\in {ℝ}^{nx}$, stands for the absolute disturbance vector generated by the swarm. According to biologists, each component of this vector, ${\mathsf{\xi }}_i^k\left(\mathbf{S}\right)$ ($i\in \overline{1,nx}$), is generated at random as well, but not uniformly. This time, the random numbers generator uses Lévy’s (probability) distribution:

$$\mathrm{p}\left(u\right)=u^{-\mathsf{\lambda }},\forall \hspace{0.17em}u\in {ℝ}_{+}^{\ast },$$

(29)

where $\mathsf{\lambda }\in [0,3]$ is preset. Usually, $\mathsf{\lambda }=1.5$. To generate numbers with prescribed probability distribution, Baker’s algorithms can be implemented and run [19,32]. The Lévy’s distribution is assigned to each component of search space, i.e., to each interval $\left[x_{\min ,i},x_{\max ,i}\right]$ ($i\in \overline{1,nx}$). Thus, ${\mathsf{\xi }}_i^k\left(\mathbf{S}\right)\in \left[x_{\min ,i},x_{\max ,i}\right]$ ($i\in \overline{1,nx}$), which explains why the vector ${\mathsf{\xi }}^k$ depends on the search space $\mathbf{S}$. However, during the run of IFA, this definition revealed two caveats: (a) Lévy’s distribution over $\left[x_{\min ,i},x_{\max ,i}\right]$ yields generating the most values in the vicinity of $x_{\min ,i}$, although there is no reason the fireflies should be directed to the lower bound of the search space; and (b) as the swarm concentrates around a shiniest firefly, the disturbing direction becomes unrealistic, because the vector ${\mathsf{\xi }}^k$ does not belong anymore to the space the fireflies are moving into, i.e., ${\mathbf{S}}^{\hspace{0.17em}k}$. Therefore, two more improvements of basic FA were made: (a) Lévy’s distribution is centered around the light center of the swarm; and (b) instead of using $\mathbf{S}$ to generate the components of vector ${\mathsf{\xi }}^k$, the current search space ${\mathbf{S}}^{\hspace{0.17em}k}$ is employed. The light center can be determined by using the radiated light of each firefly:

$$x_c^k=\frac{{\displaystyle \sum _{p=1}^P{\Phi }_p^k{x}_p^k}}{{\displaystyle \sum _{p=1}^P{\Phi }_p^k}},\forall \hspace{0.17em}k\in \mathbb{N}.$$

(30)

Now Lévy’s distribution is defined as follows, for each component $i\in \overline{1,nx}$:

$$\mathrm{p}\left(u\right)={\left|u-x_{c,i}^k\right|}^{-\mathsf{\lambda }},\forall \hspace{0.17em}u\in \left[x_{\min ,i}^k,x_{\max ,i}^k\right],$$

(31)

Hence, in definition (23), the factor ${\mathsf{\xi }}^k\left(\mathbf{S}\right)$ is replaced by ${\mathsf{\xi }}^k\left({\mathbf{S}}^{\hspace{0.17em}k}\right)$. The product ${\mathsf{\eta }}^k{\mathsf{\xi }}^k\left({\mathbf{S}}^{\hspace{0.17em}k}\right)$ is a disturbance that ensures the swarm convergence to some extent.

After being computed, the next position $x_p^{k+1}$ may need to be made viable. This means bringing back all outsider components of $x_p^{k+1}$ inside the corresponding min–max interval, according to the initial search space definition (24). Assume $x_{p,i}^{k+1}\notin \left[x_{\min ,i},x_{\max ,i}\right]$ for some $i\in \overline{1,nx}$. Then, the integer below can be computed:

$$n_{p,i}^{k+1}=⌊\frac{x_{p,i}^{k+1}-x_{\min ,i}}{x_{\max ,i}-x_{\min ,i}}⌋,\mathrm{if}{x}_{p,i}^{k+1}>x_{\max ,i}\mathrm{or}{n}_{p,i}^{k+1}=⌊\frac{x_{\max ,i}-x_{p,i}^{k+1}}{x_{\max ,i}-x_{\min ,i}}⌋,\mathrm{if}{x}_{p,i}^{k+1}<x_{\min ,i}.$$

(32)

Consequently, $x_{p,i}^{k+1}$ can be replaced by the new component below:

$$x_{p,i}^{k+1}\leftarrow \left\{\begin{array}{ll}{x}_{p,i}^{k+1}-n_{p,i}^{k+1}\left(x_{\max ,i}-x_{\min ,i}\right)\hfill & ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_{p,i}^{k+1}>x_{\max ,i}\hfill \\ x_{p,i}^{k+1}+n_{p,i}^{k+1}\left(x_{\max ,i}-x_{\min ,i}\right)\hfill & ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{x}_{p,i}^{k+1}<x_{\min ,i}\hfill \end{array}\right..$$

(33)

It is easy to prove that the new component belongs now to the corresponding min-max interval.

The notation for radiated light intensity is not selected at random. In fact, $\Phi$ is the cost function to maximize among the fireflies’ swarm during the search (and will stand for fitness, in case of PID controllers optimization).

The numerical procedure of IFA (the Algorithm 1) is summarized next.

Algorithm 1. Main steps of IFA procedure
1	Input data: cost function to maximize, $\Phi$ (associated with light intensity eradiated by fireflies); search space boundaries (see Equation (24)); configuring parameters: swarm size $P$, specific light intensity $I_0$, Lévy’s constant $\mathsf{\lambda }$, maximum number of iterations $K$, surviving (aging) factor $\mathbf{S}$.
2	Initialization. Spread the fireflies approximately uniform over the search space $\mathbf{S}$ and evaluate their performance (by using $\Phi$). Save the shiniest firefly.
3	Main loop. While the number of iterations is less than $K$ and the shiniest firefly survived over less than $\mathbf{S}$ iterations, do the following: 3.1. Evaluate the current search space ${\mathbf{S}}^{\hspace{0.17em}k}$ and its diameter using Equation (27). 3.2. Evaluate the mass left of swarm, by using Equation (30). 3.3. For each firefly in the swarm: 3.3.1. For each shinier firefly: 3.3.1.1. Evaluate the perceived light intensity by using Equation (28). 3.3.1.2. Compute the corresponding disturbance, according to the second term of Equation (23). 3.3.2. Update the current position by adding all disturbances of shinier fireflies. 3.3.3. Set the swarm disturbance intensity at random, uniformly, in $[-1,+1]$. 3.3.4. Use Baker’s algorithm to generate the absolute swarm disturbing direction with Lévy’s distribution (31). 3.3.5. Update the current position by adding the convergence disturbance, according to the last term of Equation (23). 3.3.6. Make viable the new position by using Equations (32) and (33). 3.4. Update the swarm performance and determine the new shiniest firefly. 3.5. If the current shiniest firefly still is shinier than the current one, increment the surviving counter. 3.6. Otherwise, replace the current shiniest firefly by the new one and reset the surviving factor. 3.7. Increment the number of iterations.
4	Output data: optimal value of $\Phi$ and position of the shiniest firefly; performance parameters of IFA: total number of iterations, surviving counter.

As in cases of many metaheuristics, there are no theoretical results proving the IFA is convergent. Therefore, to terminate the numerical procedure, the stop tests of the main loop in Algorithm 1 are necessary. Additionally, the position of the current shiniest firefly changes only because of the convergence disturbance, and its light can fade.

3. Analysis of Simulation Results

3.1. Simulations Setting

The simulations were performed within the MATLAB-SIMULINK™ environment. The two closed loop configurations of Figure 4 were implemented as SIMULINK models, as depicted in Figure 10 and Figure 11. Unlike the schemes in Figure 4, the SIMULINK models include the characteristic limitations of ASTANK2 plant (enumerated in the end of Section 2.1). The output references are of two types: a step for optimal PID design, and a staircase for robustness analysis. In Figure 10 and Figure 11, the step references are set as follows: $y_{1,\mathrm{ref}}=8.5$ V, $y_{2,\mathrm{ref}}=9.5$ V. They can be switched, as displayed in the left side of both schemes (see the block sw). Additionally, the sampling period was set to unit($T_s=1$ s) in all simulations. Since ASTANK2 has a slow dynamic, the simulation duration was set to $T=1000T_s=1000$ s for step references. Longer durations are employed in cases of other references types.

In the case of staircase references, the input of each scheme is changed to the one in Figure 12. The staircases can be designed by the user within an M-script with two variables: one giving the height of each step (stair), and another one setting the time instants of each step (time_steps). The staircases have to comply with ASTANK2 limitations of outputs (no more than 2 cm difference between steps on different channels).

Additionally, the staircases can be scaled in a suitable interval $\left[{\mathrm{Y}}_{\min },{\mathrm{Y}}_{\max }\right]$ (specified by means of block Y_minmax). By default, the scaling interval is $\left[7,11\right]$ cm.

To the right side of schemes in Figure 10 and Figure 11, the characteristic static map of the ASTANK2 plant with linear interpolation is implemented (see the block u_ref). Thus, the command references are available during the simulation, depending on output references $y_{1,\mathrm{ref}}$ and $y_{2,\mathrm{ref}}$.

The Gaussian white noises are generated by means of corresponding blocks on each channel. Their variances are set according to the last column of Table 1, where the estimations resulted after linear identification of the ASTANK 2 plant.

Many signals are displayed in both schemes. Nevertheless, when simulation is employed to evaluate the regulators fitness, signals plotting is suppressed in order to speed up the search.

For the two metaheuristics (AHCA and IFA), the configuring parameters are listed in

Table 2. Configuring parameters of metaheuristics.

AHCA—Exploration							IFA—Exploitation
$P$	$K$	$S$	$M$	$\mathsf{\delta }$	$P_e$		$P$	$K$	$\mathbf{S}$	$I_0$	$\mathsf{\lambda }$	${\mathsf{\rho }}_{\min }$	${\mathsf{\rho }}_{\max }$
40	$30P$	30	30	1%	3			25	100	30	1	1.5	0.5	2

, where, beside $P$, $K$, $\mathbf{S}$, $I_0$, $\mathsf{\lambda }$, the following notations are used: $M$ is the maximum number of mutations per climber, $\delta$ is the climber’s altitude relative accuracy, $P_e$ is the number of elite climbers (see [29] for details of AHCA), ${\mathsf{\rho }}_{\min }$ is the relative lower bound of the search space, ${\mathsf{\rho }}_{\max }$ is the relative upper bound of the search space.

The last two notations are useful to confine the search space for IFA. Assume AHCA provided a sufficiently accurate optimal solution, $x_{\mathrm{o}}^{\mathrm{AHCA}}$. Then the boundaries of search space for fireflies are: $x_{\min }={\mathsf{\rho }}_{\min }{x}_{\mathrm{o}}^{\mathrm{AHCA}}$ and $x_{\max }={\mathsf{\rho }}_{\max }{x}_{\mathrm{o}}^{\mathrm{AHCA}}$.

3.2. Metaheuristic Optimization Results

Before running AHCA and IFA, the PID-R regulators were tuned by means of the well-known Ziegler–Nichols empirical method [7]. The resulting parameters of PID-R controllers are listed in

Table 3. Parameters of PID-R controllers tuned by the Ziegler–Nichols method.

Identification Model	Channel	${\mathit{K}}_{\mathit{u}}$	$$\begin{gathered} {\mathit{T}}_{\mathit{u}} \\ \left[\mathbf{s}\right] \end{gathered}$$	${\mathit{K}}_{\mathit{p}}$	$$\begin{gathered} {\mathit{T}}_{\mathit{i}} \\ \left[\mathbf{s}\right] \end{gathered}$$	$$\begin{gathered} {\mathit{T}}_{\mathit{d}} \\ \left[\mathbf{s}\right] \end{gathered}$$	$$\begin{gathered} \mathsf{\Phi } \\ [\%] \end{gathered}$$
<raw>	1	1.3851	355.31	0.83106	177.66	44.41	70.12
	2	1.2519	40.83	0.75114	20.41	5.10
<spline>	1	1.233625	57.09	0.740175	28.54	7.14	73.12
	2	2.479675	28.15	1.487805	14.07	3.52

(where $K_u$ is the critical gain that pushes the system to stability limit, with oscillations of period $T_u$).

After simulation with step references, the I/O signals of PID-R controllers are displayed in Figure 13 (<raw> idmo) Figure 14 (<spline> idmo). Channels 1 and 2 are represented on top and at bottom of figures, respectively. Additionally, variation of commands and outputs are depicted to the left side and right side of figures, respectively. The commands are represented in blue (left side) and the outputs in red (right side). All references are tracked in both cases, but the I/O signals are smoother for <spline> idmo than for <raw> idmo. This is proven by the corresponding fitness values in Table 3 as well (73.12% versus 70.12%). It is well known that Ziegler–Nichols method usually yields oscillating commands, but the variations for <raw> idmo are unacceptable (Figure 13). In the case of <spline> idmo (Figure 14), while the first command can be marginally accepted, the second command is still too variable. The results after running the two metaheuristics are now analyzed. In

Table 4. Running performances of metaheuristics.

AHCA—Exploration							IFA—Exploitation
PID	Identification Model	$\mathit{k}$	$\mathit{n}\mathit{s}$	$\mathit{m}$	${\mathit{n}}_{\mathsf{\Phi }}$	$$\begin{gathered} {\mathsf{\Phi }}_{\mathbf{o}} \\ [\%] \end{gathered}$$	PID	Identification Model	$\mathit{k}$	$\mathit{n}\mathit{s}$	${\mathit{n}}_{\mathsf{\Phi }}$	$$\begin{gathered} {\mathsf{\Phi }}_{\mathbf{o}} \\ [\%] \end{gathered}$$
R	<raw>	64	30	4565	5421	80.03	R	<raw>	96	30	28,595	83.90
	<spline>	61	30	5143	6127	83.08		<spline>	100	29	29,799	85.01
N	<raw>	62	30	4192	4982	81.68	N	<raw>	67	30	19,866	83.84
	<spline>	53	30	4484	5322	83.50		<spline>	66	30	19,565	85.07

, some running performances of AHCA and IFA are given, such as the final number of iterations, $k$; the surviving counter of best climber/firefly, $ns$; the total number of mutations applied to climbers, $m$; the total number of fitness evaluations $n_{\Phi }$; and the optimal fitness, ${\Phi }_{\mathrm{o}}$.

The algorithms were run on a regular computer at current technological levels. Usually, the exploration by AHCA takes about 20 min, whereas the exploitation by IFA is about 90 min long. (The optimal PID-N controllers were found to be a few minutes faster than the optimal PID-R controllers.) The difference between the two durations is due to the fact that the climbers are only employed to find the most promising zone where the optimum could exist. Inside that zone, as the fitness rather exhibits fractal variations, the climbers cannot use their compasses based on local gradient evaluation and, therefore, AHCA usually oscillates. Fireflies do not use the gradient to fly, but IFA is not as efficient as AHCA. Consequently, the configuring parameters were set to allow longer search durations and, thus, to increase the optimal point accuracy. The number of fitness evaluations, $n_{\Phi }$, is much higher for IFA than for AHCA, which means the search for optimum is more intense in the case of fireflies than in the case of climbers. One can see that both algorithms complete the search after the best optimal climber or firefly has survived at maximum (see parameters $S$ and counters $ns$ in Table 2 and Table 4, respectively) and not after maximally iterating the positions updating (see parameters $K$ and counters $k$, respectively). There is one exception; the optimal PID-R for <spline> idmo was found with a maximum number of iterations (100), but the surviving factor was almost maximum too (29 instead of 30). For AHCA, an interesting parameter is $m$, with values starting from 4192. Mutations are applied to blocked climbers on local peaks, for which the fitness cannot increase by using the gradient. The purpose of mutation is to throw the climber on another slope of fitness hyper-surface in the hope that, from there, it can continue to climb. Good exploration of search space is thus ensured. This short analysis proves that, very likely, the best local optimum was found by combining the two metaheuristics.

Looking again at Table 4, one can see that the optimal fitness of climbers is smaller than the optimal fitness of fireflies. Nevertheless, the difference is not so big. Furthermore, this difference can barely be observed when comparing the corresponding I/O signals in Figure 15, Figure 16, Figure 17 and Figure 18. (The signals are arranged exactly like in the case of Figure 13 and Figure 14.) Additionally, the optimal PID controllers for <spline> idmo are slightly fitter than the optimal ones for <raw> idmo, as expected. However, one can notice that both PID-N and PID-R controllers exhibit almost the same fitness for the same idmo, which does not meet the initial expectations. Normally, the commands the PID-N controller provides should have been smoother than the commands of the PID-R controller. Or, the command smoothness is almost the same for both PID controllers. One can notice, though, a slight improvement for PID-N commands in the case of <spline> idmo.

The parameters of optimal PID controllers found by IFA are listed in

Table 5. Parameters of optimal PID controllers found by means of IFA.

PID	Identification Model	${\mathit{K}}_{\mathit{p},1}$	${\mathit{T}}_{\mathit{i},1}\left[\mathbf{s}\right]$	${\mathit{T}}_{\mathit{d},1}\left[\mathbf{s}\right]$		${\mathit{K}}_{\mathit{p},2}$	${\mathit{T}}_{\mathit{i},2}\left[\mathbf{s}\right]$	${\mathit{T}}_{\mathit{d},2}\left[\mathbf{s}\right]$
R	<raw>	0.2413502	45.6150814	0.0001611		1.4740651	42.7184059	0.0005492
	<spline>	0.397183	76.6958348	0.00002		0.6197605	28.9826891	0.0012607
N	<raw>	0.0406009	3.1435647	45.1171123		1.1140813	37.4705655	0.0013905
	<spline>	0.2778442	55.2662226	0		0.6518653	28.9044758	0.1017565

. The optimal PID-R parameters are very different from the ones in Table 3, especially in the case of derivative time constant ($T_d$). Not only the parameters are different, but the I/O signals are as well, as revealed by Figure 15, Figure 16, Figure 17 and Figure 18. Interestingly, the best PID controller found (PID-N for <spline> idmo) works with a PI rule on the first channel (see the null $T_{d,1}$ in Table 5), whilst the derivative time constant on the second channel is at least 70 times larger than the other derivative time constants on that channel. Simulations proved that the command on channel 2 mastered very well both output channels, which, in turn, were less sensitive to the variations of command on channel 1.

Moreover, as seen in Figure 15 and Figure 17, corresponding to <raw> idmo, the command $u_2$ exhibits saturation of e-v voltage at 10 V on the first 100 s of transient zone. Though this may have speeded up the references tracking on output channels, saturating the voltages at their upper limit is to be avoided in order to protect the e-v. Fortunately, saturation does not occur in the case of <spline> idmo (Figure 16 and Figure 18), which constitute an important advantage of this model.

Another interesting characteristic concerns the outputs. In all variations above, the output signals are like step responses of a second order system, corrupted by some noises.

This allows some supplementary performance parameters (beside the fitness) to be estimated, provided that the noise can be filtered. Such a filter can be designed by identification of each noisy output with the ideal step response of order II.

In Appendix A, the identification method is briefly described. Let $y_{j,\mathrm{II}}^{(\mathsf{\zeta },\mathsf{\omega },\mathsf{\tau })}$) be the identified step response of order II corresponding to output $y_j$ ($j\in \{1,2\}$. Then, after performing identification on outputs in Figure 15, Figure 16, Figure 17 and Figure 18, the matched step responses of order II are depicted in Figure 19.

Starting from $y_{j,\mathrm{II}}^{(\mathsf{\zeta },\mathsf{\omega },\mathsf{\tau })}$ ($j\in \{1,2\}$), several specific performances can be evaluated. In Appendix B, the most important performance parameters are defined.

Table 6. Optimal PID controllers performance parameters related to second order step response.

PID	Identification Model	Channel	$\mathsf{\zeta }$	$\mathsf{\omega }$ [rad/s]	$\mathsf{\tau }$ [s]	$\mathit{T}$	$$\begin{gathered} \mathbf{RO} \\ [\%] \end{gathered}$$	$$\begin{gathered} {\mathsf{\tau }}_{\mathit{p}} \\ [\%] \end{gathered}$$	$$\begin{gathered} {\mathsf{\tau }}_{\mathit{d}} \\ \left[\mathbf{s}\right] \end{gathered}$$	$$\begin{gathered} {\mathsf{\tau }}_{\mathit{r}} \\ \left[\mathbf{s}\right] \end{gathered}$$	$$\begin{gathered} {\mathsf{\tau }}_{\mathit{s}} \\ \left[\mathbf{s}\right] \end{gathered}$$	${\mathit{N}}_{\mathbf{t}}$
R	<raw>	1	0.8243590	0.0243468	11.03	49.82	1.030554	238.98	75.8	195.32	175.57	165
		2	0.6028431	0.0309242	0	53.64	9.31346	127.33	45.98	89.89	172.81	195
	<spline>	1	0.8506410	0.0267793	9.03	43.90	0.620133	232.16	68.60	192.84	168.76	165
		2	0.6946152	0.0304711	0	47.25	4.814988	143.32	48.78	106.69	157.10	108
N	<raw>	1	0.8287348	0.0250022	11.03	48.26	0.954089	235.55	74.23	193.10	183.62	165
		2	0.6724633	0.0338145	0	43.98	5.759211	125.53	43.49	92.23	144.98	192
	<spline>	1	0.8506410	0.0262748	9.03	44.74	0.620133	236.45	69.75	196.37	171.83	165
		2	0.6787924	0.0305118	0	48.28	5.480398	140.21	48.35	103.41	159.55	194

encompasses the values of performance parameters corresponding to signals in Figure 19.

The damping factor $\mathsf{\zeta }$ in the table shows that all outputs are under damped, although the relative overshoot $\mathrm{RO}$ is noticeably bigger on the second channel compared to the first channel (where this parameter is very small). Nevertheless, the maximum RO is less than 10% of reference for all controllers, which constitutes a good performance. The natural pulsation $\mathsf{\omega }$ is slightly bigger on the second channel. This property is confirmed by the parameter $T_u$ in Table 3 (smaller for the second channel) and the small, damped oscillations of $y_2$ in Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. On channel 1, an intrinsic delay $\mathsf{\tau }$ between 9 and 11 s was detected, which corresponds to variations of $y_1$ in Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. The controllers working with <spline> idmo seem to be slightly faster than the ones synthesized for <raw> idmo, not only in terms of intrinsic delay, but also when looking at the other time parameters: shorter time delay ${\mathsf{\tau }}_d$, rise time ${\mathsf{\tau }}_r$ and settling time ${\mathsf{\tau }}_s$. However, PID-N controllers are not necessarily better than PID-R controllers, according to the parameters in Table 6. In the last column of the table, the estimated transient time $N_{\mathrm{t}}$ is shown for each controller. This parameter is compared with the settling time, ${\mathsf{\tau }}_s$, to validate the estimation method described in Section 2.3 One can notice that the approximation is fair. Nevertheless, recall that ${\mathsf{\tau }}_s$ is computed for a clearly specified tolerance tube around the reference in the steady-state zone (set to $\pm 5$%), whereas this specification is not employed in estimation of $N_{\mathrm{t}}$.

3.3. Robustness Analysis

Controllers can or cannot be robust in several respects. Four types of robustness are investigated within this subsection: (a) against sudden change of references; (b) against constant perturbations on outputs; (c) against stochastic perturbations on outputs; and (d) against idmo uncertainties (stochastic variations of model parameters). The reason such an analysis is developed hereafter is twofold. On one hand, very often, installations like ASTANK2 do not comply with the theoretical hypotheses founding the very sound robustness theory devised in the literature. The hypotheses regarding the nature and power of exogenous noises usually are violated. Since such hypotheses only are sufficient (but not necessary) conditions of robustness, the fact they are not verified does not mean the controller cannot be robust to some extent. So, robustness has to be tested in order to estimate robustness margins without employing theoretical results in this aim. On the other hand, robustness analysis of such a real-world plant (especially the robustness against model uncertainties) are very useful in practice. Knowing the robustness margins, even with satisfactory accuracy, can help in predicting the mean time before failure or the right instant to act in correcting defects (not too soon, not too late). This may result in a better maintenance of installation, with optimally small costs.

3.3.1. Robustness against Sudden Change of References

The Simulink model of Figure 10 and Figure 11 with inputs of Figure 12 were employed to test this type of robustness. More specifically, the step references were replaced by staircases. The performance of each PID controller was evaluated by means of fitness. Unlike fitness in (22), in this case, the fitness was computed according to the procedure described in Appendix C. Several staircase references were tested. Overall, if the references steps comply with ASTANK2 limitations, in all tested scenarios, the outputs succeed to track them, which means this type of robustness is accomplished (although with different fitness values, for different references).

To illustrate the tracking performance, the matrices of variations in Figure 20 and Figure 21 were built. In Figure 21, the results correspond to switched staircase references of Figure 20. Additionally, each row of matrix exhibits I/O signals of a controller, as specified on the graphical windows. (Many of the following figures have the same structure.)

The best fitness values are obtained by the PID-N controller with <raw> idmo (77.84% and 76.96%), for which the commands better follow their static references and the outputs better track the specified references. The second best is the PID-R controller with <spline> idmo (75.18% and 74.99%). Additionally, one can easily see that channel 1 is less robust than channel 2, where the references are sensibly better tracked. This phenomenon probably is due to the fact the first tank of installation, with a slanted face, induces stronger dynamic non linearities than the second tank (a regular parallelepiped). Hence, the first linear MISO idmo is not as accurate as the second MISO idmo when approximating the nonlinearities of the ASTANK2 plant.

3.3.2. Robustness against Constant Perturbations on Outputs

Two constant perturbations were artificially added to each output, while keeping constant the reference: the first one above the reference and the second one, later, beneath the reference. Additionally, the perturbations were not synchronized between channels. The simulation results are displayed by the matrix of variations in Figure 22. All four controllers succeeded in rejecting the constant perturbations, but no one had sensibly better fitness. The commands were slightly smoother in the case of <spline> idmo.

3.3.3. Robustness against Stochastic Perturbations on Outputs

The previous simulations were performed by using two sources of Gaussian white noise with variances estimated after ASTANK2 identification (see the last column of Table 1).

To investigate controllers robustness against stochastic perturbations, the variances of noises generators were varied in range $\left[0,{\mathsf{\lambda }}_{\max }^2\right]$, where ${\mathsf{\lambda }}_{\max }^2$ was experimentally set to 1 (after few simulations). The problem is how to define and estimate the robustness margin. The fitness can be used for this aim. Obviously, the fitness decreases when noises increase in power, as the variation in Figure 23a illustrates. Here, the horizontal index points to the noise variance. In this example, $N=31$ noise variance levels are employed within the interval $[0,1]$. Thus, if the index is set to $n\in \overline{1,N}$, then the two noises have the same variance ${\mathsf{\lambda }}_{1,n}^2={\mathsf{\lambda }}_{2,n}^2$.

The problem with fitness is that, apparently, the robustness margin can only be set by subjective means while visually inspecting the I/O signals of the plant. Clearly, the marginal behavior of a closed loop system is expressed by outputs that oscillate around references with approximately constant amplitudes. However, there is a fine line between various oscillating behaviors. Hence, setting the right robustness margin becomes a difficult task. To remove the subjectiveness, a supplementary parameter can be employed. Consider the tracking errors on the steady-state zone for all I/O signals (commands and outputs). Evidently, when the noise variance increases, the std of each tracking error increases as well.

Therefore, a global std on the steady-state zone can be computed as follows (with natural notations related to I/O signals):

$${\mathsf{\sigma }}^{\mathrm{ssz}}=\mathsf{\alpha }\left({\mathsf{\sigma }}_{{\mathsf{\epsilon }}_{y,1}}^{\mathrm{ssz}}+{\mathsf{\sigma }}_{{\mathsf{\epsilon }}_{y,2}}^{\mathrm{ssz}}\right)+\frac{\mathsf{\beta }}{2}\left({\mathsf{\sigma }}_{{\mathsf{\epsilon }}_{u,1}}^{\mathrm{ssz}}+{\mathsf{\sigma }}_{{\mathsf{\epsilon }}_{u,2}}^{\mathrm{ssz}}\right),$$

(34)

where the weights $\mathsf{\alpha }$ and $\mathsf{\beta }$ are the ones employed to compute fitness (see definition (21)), e.g. $\mathsf{\alpha }=0.35$ and $\mathsf{\beta }=0.3$. Such a global std varies as shown in Figure 23b. The variation is increasing with the noise index, as expected.

Since the fitness decreases and the variance increases, a natural way to quantify the robustness margin is to normalize both variations and to define a cost function by taking their product. More specifically, one can define $\mathbf{J}\equiv {\Phi }_{\mathrm{norm}}{\mathsf{\sigma }}_{\mathrm{norm}}^{\mathrm{ssz}}$, where:

$${\Phi }_{\mathrm{norm}}\equiv \frac{\Phi -\underset{n\in \overline{1,N}}{\min }\left\{\Phi [n]\right\}}{\underset{n\in \overline{1,N}}{\max }\left\{\Phi [n]\right\}-\underset{n\in \overline{1,N}}{\min }\left\{\Phi [n]\right\}}\mathrm{and}{\mathsf{\sigma }}_{\mathrm{norm}}^{\mathrm{ssz}}\equiv \frac{{\mathsf{\sigma }}^{\mathrm{ssz}}[n]-\underset{n\in \overline{1,N}}{\min }\left\{{\mathsf{\sigma }}^{\mathrm{ssz}}[n]\right\}}{\underset{n\in \overline{1,N}}{\max }\left\{{\mathsf{\sigma }}^{\mathrm{ssz}}[n]\right\}-\underset{n\in \overline{1,N}}{\min }\left\{{\mathsf{\sigma }}^{\mathrm{ssz}}[n]\right\}}.$$

(35)

Figure 23c shows the variation of cost function $\mathbf{J}$. The maximum point, $n_{\max }$, can be set to the index of robustness margin. The robustness margin is then defined as $\left\{n_{\max },\Phi \left[n_{\max }\right]\right\}$. If the exogenous noises lead to fitness values at least equal to $\Phi \left[n_{\max }\right]$, then one considers the controller is robust. Otherwise, the controller is not performing well. Nevertheless, it is possible that the closed loop system remains at stability limit even outside the robustness zone.

After simulations with step references and variable noise variances, the results in Figure 24 were obtained. In fact, in the figure, the noise variance was replaced by the estimated SNR, on each channel, computed on the steady-state zone of outputs. The fitness surfaces in Figure 24 are split in two regions by a horizontal plane determined after maximization of cost function $\mathbf{J}$. One considers that all points on the surface above the plane are due to the robust behavior of controller. The height of each separation plan is listed in

Table 7. Robustness margins of PID controllers depending on output SNR.

PID	Identification Model	Marginal Fitness [%]	Marginal SNR [dB]	Optimal Fitness [%]	SNR1 [dB]	SNR2 [dB]
R	<raw>	65.37	25.99	83.90	40.37	41.34
	<spline>	64.12	25.95	85.01	40.80	41.23
N	<raw>	66.45	25.98	83.84	40.34	41.46
	<spline>	64.01	25.96	85.07	40.80	41.25

, together with the corresponding marginal SNR. The lower the separation plan height, the larger the robustness zone of the controller. One can notice that the PID-N controller for <spline> idmo is slightly more robust than the other three controllers. Additionally, every controller can be robust if one noise is strong while the other one is weak.

Another remark concerns the robustness surface: large fitness values cover a larger area for <raw> idmo than for <spline> idmo. This means that, in general, the stochastic noises are better rejected by PID-R and PID-N controllers determined for <raw> idmo than by the other two controllers.

The I/O signals at the robustness margin are shown by the matrix of variations in Figure 25. Such behavior is considered as robust. Many other situations with similar I/O variations exist, but the robustness is decided by fitness margins in Table 7.

3.3.4. Robustness against Identification Model Uncertainties

The uncertainties are induced by changing the coefficients of idmo polynomials. Assume $\mathrm{P}$ is any of the model polynomials, expressed as:

$$\mathrm{P}\left({\mathrm{q}}^{-1}\right)=p_0+p_1{\mathrm{q}}^{-1}+\cdots +p_{np}{\mathrm{q}}^{-np}.$$

(36)

Then, the coefficient $p_i$ ($i\in \overline{1,np}$) is replaced by a (uniform) random value in the interval $\left[1-\mathsf{\eta },1+\mathsf{\eta }\right]p_i$, if $p_i\ne 0$ or in the interval $\left[-\mathsf{\eta },+\mathsf{\eta }\right]$, if $p_i=0$. The parameter $\mathsf{\eta }\in [0,1]$ is referred to as uncertainty power.

The problem is to define and determine the robustness margin of a PID regulator, while working with uncertain model polynomials and variable uncertainty powers. Obviously, the PID performance is expressed in terms of fitness. Nevertheless, to solve the problem of uncertainty margin detection, the fitness needs to be put into correspondence with the reference tracking property. One can consider that the uncertainty margin is reached as soon as the closed loop system is unable to track at least one output reference signal. In this case, the polynomials are affected by uncertainties of some power, whilst the closed loop system exhibits a marginal fitness. Thus, the uncertainty margin can be defined by the smallest fitness value for which the closed loop system still manages to track all given references. This lower limit also reveals the marginal (maximal) values of uncertainty power. For example, the robustness margin of 75% (in terms of fitness) can be obtained by affecting the plant model with uncertainties of estimated marginal (maximal) power of 10% ($\mathsf{\eta }\in [0,1]$). Increasing the uncertainty power leads to the incapacity of the closed loop system to track references.

In the following, an example is given. The output $y_1$ of ASTANK2 installation is acquired after applying uncertainties on the numerator and denominator of the useful transfer function in MISO₁ idmo (see Figure 1b again). The uncertainties were generated at random, by using the uniform pseudo-random generator of MATLAB (rand function).

The uncertainties can affect the 3 polynomials of useful transfer function, namely: ${\mathrm{A}}_1$, ${\mathrm{B}}_{1,1}$ and ${\mathrm{B}}_{1,2}$. As Equation (1) shows, the first coefficient of polynomials is normalized: 1 for ${\mathrm{A}}_1$ and 0 for ${\mathrm{B}}_{1,1}$, ${\mathrm{B}}_{1,2}$. Uncertainties cannot affect this coefficient. By inducing uncertainties on the remaining coefficients, a fitness surface can be generated, depending on two uncertainty powers: one affecting the denominator ${\mathrm{A}}_1$ and another one modifying the numerators ${\mathrm{B}}_{1,1}$, ${\mathrm{B}}_{1,2}$. The ASTANK2 idmo usually have polynomials with many coefficients (up to 50, in some cases). Trying to test all combinations of uncertainty applied on model coefficients is unrealistic. Therefore, the Monte-Carlo method can be employed to estimate the uncertainty margin. Hence, 100 uncertainty cases are uniformly generated at random on each axis. The resulted fitness surface looks like that in the left side of Figure 26. The estimated fitness margin is represented by the horizontal plane intersecting the fitness surface. All points above the plane correspond to the robust behavior of the closed loop system. The other fitness values, beneath the plane, rather reveal the incapacity of the closed loop system to track at least one reference. As the figure clearly shows, there are cases in which the robustness is preserved even though one uncertainty power is high and the other one is low. Consequently, it is quite difficult to determine the marginal powers if they are different. It is more suitable to determine such a marginal power along the main diagonal of the separation plane, where the two powers are equal. In the right side of Figure 26, the fitness variation along the diagonal (with equal uncertainty powers) is displayed. The horizontal line refers to the separation plane.

One can see that the most fitness values are beneath this line, as expected. Increasing the uncertainty power naturally leads to the loss of robustness. Therefore, to estimate the robustness margin, one should focus on the zone with small uncertainty powers. In the figure above, one can easily notice that the marginal uncertainty power is less than 10%.

The problem with fitness is the same as in the previous subsection; apparently, the robustness margin can only be set by subjective means, while visually inspecting the I/O signals of the plant. For example, the first two output signals in Figure 27 reveal that the system is unable to track the reference, but is very close to this goal. The fitness equals 73.98% with 5% uncertainty power and 74.03% with 4.9% uncertainty power, respectively. By decreasing the uncertainty power, the last two output signals in Figure 27 were obtained. At 74.37% fitness with 4.75% uncertainty power, the reference started to be tracked in the end of variation (see the third signal of figure). Furthermore, decreasing the uncertainty power to 4.6% leads to a fitness of 74.5% and to an output that succeeds in reaching the reference sooner than previously (see the last signal). In this case, the robustness margin can be set to 74.4% fitness with 4.7% uncertainty power. The estimation of uncertainty margin is fair but subjectively set.

The robustness margins can thus be estimated by varying the separation plane and analyzing the cases leading to fitness values close to the plane height. If the uncertain idmo that leads to fitness values above and beneath the plane height still allow tracking of the references, then the plane should be lowered. On the contrary, if the idmo with fitness values around the plane height yield loss of robustness, then the plane should be raised. The separation plane is placed at the right height when the idmo with fitness below that height cannot allow reference tracking, whereas for the other idmo, with fitness above the height, the references are eventually tracked. Because of the noises corrupting the data, such a balanced behavior can be discovered experimentally, as a theoretical approach on this topic is difficult (if not impossible) in this case.

For ASTANK2 installation, three scenarios can be investigated: (a) affect MISO₁ model by uncertainties and the keep the genuine MISO₂ model; (b) switch MISO models (keep the genuine MISO₁ and work with uncertain MISO₂); (c) induce uncertainties on both MISO models. After several simulations, it has been noticed that output $y_1$ is more sensitive to idmo uncertainties than output $y_2$. This means the controller loses robustness because of $y_1$ in the first. Therefore, in the following figures, only $y_1$ is displayed.

Figure 28, Figure 29 and Figure 30 display the simulation results corresponding to the three scenarios mentioned above. Each figure includes a matrix of variations, as usual, but with different representations from the previous figures.

In all figures, the fitness surface obtained by Monte-Carlo method takes the first column of matrix. The second column displays the diagonal of fitness surface. The last two columns are devoted to output $y_1$ in two cases: non robust (column 3) and robust (column 4). The corresponding fitness values of closed loop system for columns 3 and 4 are below and above the estimated marginal fitness of robustness zone, respectively.

The simulation scenarios providing the results in Figure 28, Figure 29 and Figure 30 offer further parameters allowing analyzing robustness. Hence, the robustness zone can be confined by: the marginal fitness ${\Phi }_{\min }$; the marginal uncertainty power ${\mathsf{\eta }}^{\max }$; the relative density $\mathsf{\rho }$ of robust behavior; the relative robustness margin $\mathsf{\mu }$. The first two parameters (${\Phi }_{\min }$ and ${\mathsf{\eta }}^{\max }$) are determined experimentally. The height of separation plane is set to ${\Phi }_{\min }$. The fitness values ${\Phi }_{\mathrm{down}}$ and ${\Phi }^{\mathrm{up}}$ corresponding to variations of output $y_1$ in Figure 28, Figure 29 and Figure 30 are close to each other.

Moreover, ${\Phi }_{\mathrm{down}}<{\Phi }_{\min }<{\Phi }^{\mathrm{up}}$, which proves that ${\Phi }_{\min }$ is fairly accurate. The marginal uncertainty power ${\mathsf{\eta }}^{\max }$ was determined by using the fitness diagonal drawn in column 2 of graphical matrices. The principle illustrated in Figure 26 was used in this aim.

The other two parameters ($\mathsf{\rho }$ and $\mathsf{\mu }$), although conventionally defined, allow performing a top of controllers with respect to robustness against uncertainties. Since the Monte-Carlo method is employed to draw fitness surfaces, one can compute the ratio between the number of robust behaviors and the number of all behaviors (10.000). If this ratio is multiplied by 100, the relative density $\mathsf{\rho }$ of robust behavior is estimated. Now, one controller is more robust than another if ${\Phi }_{\min }$ is as small as possible, whilst ${\mathsf{\eta }}^{\max }$ and $\mathsf{\rho }$ are as big as possible. Then, the relative robustness margin $\mathsf{\mu }$ can be defined as:

$$\mathsf{\mu }=\frac{100}{1+\mathsf{\sigma }\frac{{\Phi }_{\min }}{\mathsf{\rho }\hspace{0.17em}{\mathsf{\eta }}^{\max }}}[\%],$$

(37)

where $\mathsf{\sigma }>0$ is a constant scaling factor set to 0.04, such that $\mathsf{\mu }$ can vary in a suitable range (recall the explanation between Equations (15) and (16)).

The robustness parameters are listed in

Table 8. Parameters of robustness against model uncertainties.

PID	Identification Model	Uncertainties Affected	${\mathsf{\Phi }}_{\mathbf{min}}$ [%]	$$\begin{gathered} {\mathsf{\eta }}^{\mathbf{max}} \\ [\%] \end{gathered}$$	$$\begin{gathered} {\mathsf{\Phi }}_{\mathbf{down}} \\ [\%] \end{gathered}$$	$$\begin{gathered} {\mathsf{\Phi }}^{\mathbf{up}} \\ [\%] \end{gathered}$$	$$\begin{gathered} \mathsf{\rho } \\ [\%] \end{gathered}$$	$$\begin{gathered} \mathsf{\mu } \\ [\%] \end{gathered}$$
R	<raw>	MISO1	74.5	9	74.48	75	2.81	89.47
		MISO2	70.5	15.5	70	70.49	5.84	96.98
		MISO1 and MISO2	70	2.25	69.76	70	1.06	45.97
	<spline>	MISO1	60	2.5	58.11	60	0.7	42.03
		MISO2	65.5	1	64.74	65.52	0.45	14.68
		MISO1 and MISO2	60	1	57.54	60.22	0.17	6.49
N	<raw>	MISO1	74.5	8	74.04	74.51	2.54	87.21
		MISO2	70.5	9	70.43	70.65	4.51	93.5
		MISO1 and MISO2	68	3.5	67.71	68.47	1.3	62.65
	<spline>	MISO1	60	1	59.85	60.7	1.08	31
		MISO2	66	1	65.64	66.26	0.12	4.22
		MISO1 and MISO2	65	1	64.54	65.21	0.16	5.69

, for all four controllers, in all scenarios. Without the last column, Table 8 is rather confusing and says almost nothing. The $\mathsf{\mu }$ parameter can enlighten the closed loop behavior against model uncertainty.

The first clear insight is that the PID controllers for <raw> idmo are noticeably more robust to uncertainties than the PID controllers for <spline> idmo. This can be explained by the fact the <spline> idmo has a reduced number of poles, which are more exposed to the loss of stability in a way that the PID controllers cannot compensate. Secondly, among controllers for <raw> idmo, PID-R is slightly more robust than PID-N when coping with uncertainties that affect one MISO model at a time, whereas when both models are corrupted, PID-N seems to be the best (62.65% versus 45.97%). Thirdly, for <spline> idmo, the PID-R controller is better than the PID-N controller, in terms of robustness, in all scenarios. To conclude this analysis, one can notice that, regardless of the idmo type, robustness of MISO₂ is superior to robustness of MISO₁, as expected (since, in all simulations, output $y_1$ loses its tracking ability before output $y_2$). Moreover, another expectation is met: the relative robustness margin of the entire ASTANK2 idmo (in the last scenario) is the smallest one (compared to the other two scenarios).

4. Concluding Remarks and Future Developments

The main goal of this article was to perform an extensive analysis of metaheuristic PID control applied to a nonlinear installation such as ASTANK2, after identification with linear multi-MISO models. Two PID configurations were investigated, PID-R and PID-N, for each identification model type, <raw> and <spline>. Two metaheuristics were employed in the design of four optimal controllers: AHCA to extensively explore the search space, and IFA to intensively exploit a restricted zone around some promising optimal point found by climbers. Four improvements of classical FA were introduced in order to fasten the search and increase the final results accuracy. Defining a fair cost function (referred to as fitness in the article) is not an easy task. Many other authors coped with the same difficult problem, as reported in their publications. Unfortunately, to the best of our knowledge, no general cost function can be defined to be employed in any automatic control optimization problem. Therefore, adapting the cost function to the installation characteristics (and limitations) seemingly is a wise way to proceed.

The simulation results revealed the performance as well as the limitations of designed controllers. Overall, one can say that any of the four PID controllers can be employed in a closed loop with an ASTANK2 installation, as they exhibit almost the same performance. Basically, PID-R controllers are recommended to reject model uncertainties, whilst PID-N controllers are better in tracking staircase type references.

Various future research directions can be approached. A multivariable (non-decentralized) PID control structure likely yields superior performance. Another development is possible in the direction of metaheuristic multivariable RST control. Of course, the metaheuristic control based on state space representation of ASTANK2 installation is another interesting direction to investigate.