μ-Synthesis for Fractional-Order Robust Controllers

Abstract

The current article presents a design procedure for obtaining robust multiple-input and multiple-output (MIMO) fractional-order controllers using a $\mu$-synthesis design procedure with D–K iteration. $\mu$-synthesis uses the generalized Robust Control framework in order to find a controller which meets the stability and performance criteria for a family of plants. Because this control problem is NP-hard, it is usually solved using an approximation, the most common being the D–K iteration algorithm, but, this approximation leads to high-order controllers, which are not practically feasible. If a desired structure is imposed to the controller, the corresponding K step is a non-convex problem. The novelty of the paper consists in an artificial bee colony swarm optimization approach to compute the nearly optimal controller parameters. Further, a mixed-sensitivity $\mu$-synthesis control problem is solved with the proposed approach for a two-axis Computer Numerical Control (CNC) machine benchmark problem. The resulting controller using the described algorithm manages to ensure, with mathematical guarantee, both robust stability and robust performance, while the high-order controller obtained with the classical $\mu$-synthesis approach in MATLAB does not offer this.

1. Introduction

One of the active problems with major impact which have been studied for years in Control Theory refers to robustness. Robustness encompasses the sensitivity of a control system with respect to both internal and external disturbances. Several robust methods have been developed in order to achieve robust performance and stability in the presence of uncertainties. Robust control problems use ${\mathcal{H}}_2$ and ${\mathcal{H}}_{\infty }$ norms defined in frequency domain as a performance measure. To solve ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ control problems, there are several approaches. One possible solution is presented in [1] and is based on Algebraic Riccati Equations (AREs). A more numerically stable approach to solve ARE was developed using Popov triplets in [2], approach recently implemented in an open-source manner in [3], with an iterative refinement method presented in [4], but an ARE-based solution presents a limitation due to the impossibility of solving singular problems. An alternative way which manages to solve such problems was introduced in [5], where AREs were replaced by Algebraic Riccati Inequalities (ARIs). ARIs are solved through Linear Matrix Inequalities (LMIs), while regular assumptions are no longer needed due to LMI system versatility. An open-source solver for Robust Control problems using LMIs is presented in [6].

The ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ approach designs a suitable controller for the nominal plant, therefore only nominal stability and nominal performance are fulfilled. Additionally, generalized Robust Control framework allows to impose robust stability and robust performance, which cover the previous two aspects for an entire family of physical processes. As such, the $\mu$-synthesis approach extends the ${\mathcal{H}}_{\infty }$ optimization in order to obtain a robust controller for the uncertain plant which includes parametric and dynamic uncertainties [7,8]. $\mu$-synthesis is based on using structured singular values to quantify robustness margins and also on using linear fractional descriptions of the control problem containing the nominal plant model and uncertainty weighting functions. In [9], the authors present the $\mu$-synthesis problem used with the so-called D–K iteration, which, in essence, provides two steps iteratively repeated until the robust performance stops improving: designing a ${\mathcal{H}}_{\infty }$ control law and $\mu$ analysis on closed-loop system.

$\mu$-synthesis is, however, a non-convex problem and the D–K iteration represents only an approximation, without any convergence guarantees. Another significant concern of D–K iteration is that the method generates high order controllers. In order to solve this issue, various approaches based on fixed structure controller are proposed in different papers [10,11,12]. The method presented in [10] uses nonsmooth techniques for ${\mathcal{H}}_{\infty }$ synthesis. Then, using the same technique, the $\mu$-synthesis was solved using D–K iteration and the result are presented in [13].

The main issue which appears when controller structure constraints are imposed is that the optimization problem in no longer convex and also, $\mu$-synthesis is in general, considered nondeterministic polynomial time hard (NP-hard). A possible solution to that are swarm optimization algorithms. There are different approaches presented in papers [14,15,16] based on Genetic Algorithms (GA) and Particle Swarm Optimization (PSO). A solution for imposing a fixed structure controller, such as low-order or decentralized, was proposed in paper [14], which splits the problem in two parts: the convex part, solved using the classical ARE approach, and the non-convex part, solved using GA. The same authors proposed in [16] a new technique based on an evolutionary D–$K{D}_0$ iteration method, which combines the classical D-step with a $K{D}_0$ algorithm based on also a GA. Authors of the paper [15] propose an evolutionary approach to solve the $\mu$-synthesis problem without order reduction by using an improved PSO.

Artificial Bee Colony (ABC) can also be used to solve complex optimization problems with constraints and could possibly outperform the other approaches and return the best solution in shorter execution time. The initial idea was presented in [17] as an extension of another metaheuristic algorithm, namely Honey Bee Swarm (HBS). The efficiency of the algorithm for several state-of-the-art optimization problems, along with an improvement for the stopping criterion, were underlined in [18].

One possible fixed structure controller is fractional-order proportional-integral-derivative (FO-PID). FO-PID is one of the most remarkable fractional order techniques with great interest in research [19]. It is used to generalize the classical PID control by adding extra degrees of freedom [20]. Compared to the integer PID controller, FO-PID brings the advantage of improving the robustness and providing better performance. In [21], the authors propose a FO-PID controller for a fractional-order plant model, presenting an analysis in both frequency and time domains, proving that achieving better control performance is one of the advantages of this approach. The FO-PID was used on a benchmark problem, i.e., the speed control of a DC motor [22], obtaining good results in terms of performance and robustness. Two generalized versions of Kessler’s magnitude method with fractional order controllers were developed in [23,24]. A detailed comparison between classical modelling approaches and a fractional integrator approximation as a control baseline model is presented to the servo problem in [25]. A graphical method was developed in [26], while in [27] a fractional order internal model controller with event-based implementation was developed. Other applications comprise in an optimal FO-PID controller for a PMSM speed control, presented in [28], while a robust controller for a steam turbine was developed in [29].

In this paper, we present a new technique to design FO-PID robust controllers using $\mu$-synthesis. The novelty of the current approach consists in implementing an algorithm able to find the nearly optimal values of the controller parameters using an artificial bee colony optimization. The cost function to be minimized is the ${\mathcal{H}}_{\infty }$ norm of the closed-loop system, when stable, and, otherwise, a large value affinely dependent on the largest real part of the eigenvalues of the closed-loop state matrix. Therefore, the non-convex part of the NP-hard $\mu$-synthesis problem is solved using such a swarm optimization, while the D step is solved using the classical LMI technique. More than that, the realp object from MATLAB’s Robust Control Toolbox embodies a limitation because it cannot be used as an exponent, necessary in approximating the fractional element with an integer-order system. Our approach manages to deal with this limitation in order to obtain the controller parameters, because the only information necessary in our approach is the range of the controller’s parameters stored in such variables, which can be replaced with a simpler software object. Additionally, using our approach, a numerical example illustrates that the resulting controller manages to fulfill the robust stability and performance, while the controller obtained using unstructured $\mu$-synthesis does not rigorously guarantee these specifications. Therefore, our method proposes a general framework able to synthesize arbitrary fixed structure fractional order controllers by optimizing their parameters in terms of robustness and performance, surpassing the well-established approach of manually tuning them for a desired problem, harnessing the Robust Control framework’s design strongness.

We illustrate our proposed method on a benchmark problem: obtaining a controller for a two axis computer numerical control (CNC) system by solving a mixed-sensitivity $\mu$-synthesis problem. Several control methods using the state space model of the machine are presented in literature. A comparison between the classical pole-placement method and Linear Quadratic Regulator (LQR) method is presented in [30]. For the LQR problem, an energy-based minimization algorithm was proposed, which proves to be suitable in terms of stability and robustness. As presented in [31], a different approach is recommended, in order to obtain a PI controller for each axis, using the state-feedback control algorithm. In this case, the state space model is augmented with an extra state which represents the integral of the position. The PI regulator parameters are obtained from the state-feedback gains.

The paper is organized in four sections. Section 2 introduces several ideas relevant to the proposed method, such as a mathematical foundation of FO-PID, continuing with the fundamental robust control problem based on $\mu$-synthesis, and, finally, the ABC optimization algorithm. After that, the last subsection focuses on solving the non-convex problem of computing fixed structure controllers. Section 3 illustrates an application of the proposed method for position control of a two axis CNC machine, along with numerical results. In Section 4, the previously mentioned results are compared with those obtained using the well-established algorithms from MATLAB. Finally, conclusions are presented in Section 5.

2. Materials and Methods

In this section we present a controller synthesis procedure which manages to find a fixed structure fractional-order controller using the $\mu$-synthesis technique from the Robust Control framework, where the non-convex subproblem involved in the classical D–K iteration is replaced by a swarm optimization algorithm. First, the mathematical background comprised in fractional-order control, robust control and artificial bee colony optimization is underlined, while the fourth subsection presents the proposed design procedure using all the mechanisms briefly described.

2.1. Fractional-Order Controller

The classical integer-order calculus was extended by Riemann and Liouville to fractional-order calculus by introducing the fractional integral operator [32]:

$${\mathcal{I}}_a^{\alpha }\left\{f\left(t\right)\right\}=\frac{1}{\Gamma \left(\alpha \right)}{\int }_a^{t}f\left(\tau \right){(t-\tau )}^{\alpha -1}d\tau ,$$

(1)

where $\Gamma \left(\alpha \right):{\mathbb{C}}_{+}\to \mathbb{C}$ is the Euler Gamma function and the order of the integral operator is the complex parameter $\alpha \in {\mathbb{C}}_{+}$. This extension develops a new area of research in the Control System domain. A commonly used definition of this operator was introduced in [33] as:

$${\mathcal{J}}_C^{\alpha }\left\{f\left(t\right)\right\}=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)d\tau ={\mathcal{I}}_0^{\alpha }\left\{f\left(t\right)\right\},$$

(2)

with $t>0$ and $\alpha \in {\mathbb{R}}_{+}$, having the Laplace transform [34]:

$$\mathcal{L}\left\{{\mathcal{J}}_c^{\alpha }\left\{f\left(t\right)\right\}\right\}\left(s\right)=s^{-\alpha }F\left(s\right),$$

(3)

where $F\left(s\right)=\mathcal{L}\left\{f\left(t\right)\right\}\left(s\right)$. One of the most common controller structures used in practice is the proportional-integral-derivative (PID) controller, having three degrees of freedom. Using fractional-order calculus, two new degrees of freedom can be added to a PID, having the notation PI^λD^μ. As such, the fractional order PID (FO-PID) has, as extra degrees of freedom, the order $\lambda \in {\mathbb{R}}_{+}$ of the integrator and the order $\mu \in {\mathbb{R}}_{+}$ of the differentiator, with the resulting transfer function:

$$H_c\left(s\right)=K_P+\frac{K_I}{s^{\lambda }}+K_D{s}^{\mu }.$$

(4)

The time domain expression of the command signal $c\left(t\right)$ can be expressed using the error signal $\epsilon \left(t\right)$ as:

$$c\left(t\right)=K_p·\epsilon \left(t\right)+K_I·{\mathcal{J}}_c^{\lambda }\left\{\epsilon \left(t\right)\right\}+K_D·{\mathcal{J}}_c^{-\mu }\left\{\epsilon \left(t\right)\right\}.$$

(5)

One of the major issues of such a fractional element, i.e., ${\mathcal{J}}_c^{\lambda }$ or ${\mathcal{J}}_c^{-\mu }$, is its implementation. In order to solve this problem, the Oustaloup recursive approximation (ORA) was introduced [35], and allows the approximation of the fractional-order element with an LTI system of pre-specified order N:

$$s^{\lambda }=\prod _{k=1}^N\frac{1+s/{\omega }_{z,k}}{1+s/{\omega }_{p,k}},$$

(6)

where the frequency values of the singularities are obtained based on the desired fractional order $\lambda \in (0,1)$, the integer order of the approximation N, along with the frequency range in which the approximation is valid $[{\omega }_l,{\omega }_u]$. Using the following two coefficients:

$$\epsilon ={\left(\frac{{\omega }_u}{{\omega }_l}\right)}^{\frac{\lambda }{N}}\mathrm{and}\eta ={\left(\frac{{\omega }_u}{{\omega }_l}\right)}^{\frac{1-\lambda }{N}},$$

(7)

the above mentioned frequencies can be computed using:

$$\begin{array}{cc}\hfill {\omega }_{z,1}& ={\omega }_l\sqrt{\eta },\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\omega }_{p,n}& ={\omega }_{z,n}·\epsilon ,n=\overline{1,N},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {\omega }_{z,n+1}& ={\omega }_{p,n}·\eta ,n=\overline{1,N-1}.\hfill \end{array}$$

(10)

For the rest of the possible real values of $\lambda$, the approximation can be easily extended as: for $\lambda \in (-1,0)$ by inverting the relation (6), while for $\left|\lambda \right|\ge 1$, the components could be the integer part $\left[\lambda \right]$ and the fractional part $\left\{\lambda \right\}$, with the fractional part only approximated using (6).

2.2. Robust Control

The Robust Control framework assumes to minimize the ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ norm of the lower linear fractional transformation (LLFT) of a plant P and a controller K:

$$P_o=\mathtt{LLFT}(P,K)=P_{11}+P_{12}K{(I-P_{22}K)}^{-1}{P}_{21},$$

(11)

where the plant P can be written, in general form, as:

$$P:\left(\begin{array}{cc}{P}_{11}& P_{12}\\ P_{21}& P_{22}\end{array}\right)=\left(\begin{array}{ccc}A& B_w& B_u\\ C_z& D_{zw}& D_{zu}\\ C_y& D_{yw}& D_{yu}\end{array}\right),$$

(12)

where the signals involved in the above relation will be further detailed in a more general context. One approach to solving ${\mathcal{H}}_2$/${\mathcal{H}}_{\infty }$ problems is using AREs, which presents several limitations that are removed by the LMI approach. However, this framework can be used to ensure nominal stability and nominal performance only, but, the plant P must also contain an augmented model of the real process in which uncertainties are also present. There are two classic uncertainties types: parametric, represented by $\delta I$, where $\delta$ is the maximum bound of the parameter for a physical variable, and unstructured, represented by a full block $\mathsf{\Delta }\in {\mathbb{R}}^{m\times m}$. The latter illustrates neglected or unknown dynamics uncertainties. In the mixed-scenario case, the following set is considered:

$$\mathbf{\Delta }=\left\{\mathsf{\Delta }=\mathrm{diag}\left({\delta }_1{I}_{n_1},\cdots ,{\delta }_s{I}_{n_s},{\mathsf{\Delta }}_1,\cdots ,{\mathsf{\Delta }}_f\right)|{\delta }_k\in \mathbb{R},{\mathsf{\Delta }}_j\in {\mathbb{R}}^{m_j\times m_j},k=\overline{1,s},j=\overline{1,f}\right\}.$$

(13)

In the Robust Control field, one of the main tools used for robustness analysis is the structured singular value, defined as follows.

Definition 1.

For a square matrix $M\in {\mathbb{C}}^{N\times N}$ the structured singular value with respect to the set Δ is:

$${\displaystyle {\mu }_{\mathbf{\Delta }}\left(M\right)=\frac{1}{{\displaystyle \underset{\mathsf{\Delta }\in \mathbf{\Delta }}{min}\left\{\overline{\sigma }\left(\mathsf{\Delta }\right)\right|det(I-M\mathsf{\Delta })=0\}}},}$$

(14)

if there exists $\mathsf{\Delta }\in \mathbf{\Delta }$ such that the matrix $I-M\mathsf{\Delta }$ is rank deficient, otherwise 0.

For an LTI system described by the transfer matrix $M\left(s\right)$ and an upper linear fractional transformation (ULFT) connection shown in Figure 1 (left), the structured singular value ${\mu }_{\mathsf{\Delta }}\left(M\right)$ can be defined as:

$${\mu }_{\mathbf{\Delta }}\left(M\left(s\right)\right)=\underset{\omega \in {\mathbb{R}}_{+}}{sup}{\mu }_{\mathbf{\Delta }}\left(M\left(j\omega \right)\right).$$

(15)

Now, considering $M\left(s\right)=\mathtt{LLFT}(P,K)\left(s\right)$ as being the lower linear fractional transformation (LLFT) between plant P and controller K, the connection illustrated in Figure 1 (right) results. The generalized plant structure is:

$$P_{\mathsf{\Delta }}\left(s\right)=\left(\begin{array}{ccc}{P}_{vd}\left(s\right)& P_{vw}\left(s\right)& P_{vu}\left(s\right)\\ P_{zd}\left(s\right)& P_{zw}\left(s\right)& P_{zu}\left(s\right)\\ P_{yd}\left(s\right)& P_{yw}\left(s\right)& P_{yu}\left(s\right)\end{array}\right)\iff P_{\mathsf{\Delta }}:\left(\begin{array}{c}\dot{\mathbf{x}}\left(t\right)\\ \mathbf{v}\left(t\right)\\ \mathbf{z}\left(t\right)\\ \mathbf{y}\left(t\right)\end{array}\right)=\left(\begin{array}{cccc}A& B_d& B_w& B_u\\ C_v& D_{vd}& D_{vw}& D_{vu}\\ C_z& D_{zd}& D_{zw}& D_{zu}\\ C_y& D_{yd}& D_{yw}& D_{yu}\end{array}\right)\left(\begin{array}{c}\mathbf{x}\left(t\right)\\ \mathbf{d}\left(t\right)\\ \mathbf{w}\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right),$$

(16)

where three types of input signals are present—the command input $\mathbf{u}\in {\mathbb{R}}^{n_u}$, the performance input $\mathbf{w}\in {\mathbb{R}}^{n_w}$, and the disturbance input $\mathbf{d}\in {\mathbb{R}}^{n_d}$—and three types of outputs—the measurements vector $\mathbf{y}\in {\mathbb{R}}^{n_y}$, the performances vector $\mathbf{z}\in {\mathbb{R}}^{n_z}$, and the disturbance output $\mathbf{v}\in {\mathbb{R}}^{n_v}$.

Besides the well-known ${\mathcal{H}}_2/{\mathcal{H}}_{\infty }$ methods, a controller that meets robust stability and robust performance alike can be computed using the $\mu$-synthesis framework. Robust stability implies that a specific controller manages to stabilize all the processes described by the upper linear fractional transformation (ULFT) presented in Figure 1 (left), while robust performance implies that the controller is able to impose the desired closed-loop performance in the worst case scenario. In order to have a mathematical guarantee that a controller K meets the robust stability and performance, the Main Loop theorem can be used. It implies that closed-loop system meets robust stability and performance if and only if the structural singular value of the LLFT of the plant and controller, with respect to $\mathbf{\Delta }$, fulfills:

$$\underset{\omega \in {\mathbb{R}}_{+}}{sup}{\mu }_{\mathbf{\Delta }}\left(\mathtt{LLFT}(P,K)\left(j\omega \right)\right)<1.$$

(17)

Therefore, the robust control problem can be written as:

$$\underset{K\mathrm{stab}.}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}{\mu }_{\mathbf{\Delta }}\left(\mathtt{LLFT}(P,K)\left(j\omega \right)\right),$$

(18)

which is not convex. More than that, the structural singular values are hard to be explicitly computed. In practice, there are various bounds which can be used to approximate the structural singular value. One of the most used approximations of the upper bound is in [9]:

$${\mu }_{\mathbf{\Delta }}\left(M\right)\le \underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(DM{D}^{-1}\right),$$

(19)

where $\overline{\sigma }$ denotes the largest singular value, and the set $\mathit{D}$ is defined as:

$$\mathit{D}=\left\{\mathrm{diag}\left(D_1,\cdots ,D_s,d_1{I}_{m_1},\cdots ,d_f{I}_{m_f}\right)|D_k=D_k^{\top }\in {\mathbb{R}}^{n_k\times n_k},d_j>0,k=\overline{1,s},j=\overline{1,f}\right\}.$$

(20)

Based on this upper bound, a good approximation of the initial non-convex problem can be employed by solving the following quasi-convex problem:

$$\underset{K\mathrm{stab}.}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}\underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(D\left(j\omega \right)·\mathtt{LLFT}(P,K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right).$$

(21)

If the scaling factor, represented by the system D, is fixed, then the problem (21) is nothing but a ${\mathcal{H}}_{\infty }$ optimization problem. On the other hand, fixing the controller K, the D scaling step can now be obtained by solving a Parrot problem for a desired frequency set $\mathsf{\Omega }=\{{\omega }_1,\cdots ,{\omega }_N\}$ using the following generalized eigenvalue problem:

$$\begin{array}{cc}\hfill & min\gamma ,\hfill \\ \hfill \mathrm{s}.\mathrm{t}.{\left(\mathtt{LLFT}(P,K)\left(j{\omega }_i\right)\right)}^{*}·& X·\mathtt{LLFT}(P,K)\left(j{\omega }_i\right)\le {\gamma }^{2}X,\hfill \end{array}$$

(22)

where from the solution $X={\left(D\left(j{\omega }_i\right)\right)}^{*}·D\left(j{\omega }_i\right)$, the matrix $D\left(j{\omega }_i\right)$ can be extracted using a singular value decomposition. After all Parrot problems are solved, a minimum phase system is found in order to approximate the analytical solution $D\left(s\right)$. In summary, an iterative algorithm which solves the $\mu$-synthesis problem starts by setting $D=I$, with the following steps applied successively:

1: 

Fix D and solve the ${\mathcal{H}}_{\infty }$ optimal problem to find a controller K:

$$K=arg\underset{K\mathrm{stab}.}{inf}{\left|\right|\mathtt{LLFT}(P,K)\left|\right|}_{\infty }.$$

(23)

2: 

Fix the controller K and solve the set of convex problems:

$$D\left(j\omega \right)=arg\underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(D·\mathtt{LLFT}(P,K)\left(j\omega \right)·D^{-1}\right),$$

(24)

for a given frequency range $\mathsf{\Omega }$ and, then, fit a stable minimum phase transfer matrix $D\left(s\right)$.

Steps 1 and 2 are executed in a loop sequence until the difference between two consecutive ${\mathcal{H}}_{\infty }$ norms is less than a prescribed tolerance or the maximum number of iterations is reached.

2.3. Artificial Bee Colony Optimization

The artificial bee colony (ABC) optimization is a nature-inspired algorithm used to minimize a cost function:

$$f:D\to \mathbb{R},\mathrm{where}D=[l{b}_1,u{b}_1]\times [l{b}_2,u{b}_2]\times \cdots \times [l{b}_d,u{b}_d]\subset {\mathbb{R}}^d.$$

(25)

The ABC algorithm mimics the behaviour of real honeybees, where each food source represents a possible solution of the optimization problem described above. The location and the amount of nectar correspond to the design variables and the cost function, respectively. The bees are divided in two main groups: employed and unemployed bees, while the unemployed bees could be of two types as well: onlooker and scout bees. The employed bees are the ones that investigate the food source and return to the hive to inform the others by performing the waggle dance; the onlooker bees are the ones that watch the dance and decide whether or not a food source is worthy of being searched or not; the scout bees are former employed bees that have abandoned their previous food source, due to lack of nectar, and which now search for a new one. The Best Solution is represented by the food source, and the quality (or cost) of the solution is represented by the amount of nectar.

The number of employed bees coincides with number of the onlooker bees and represents the dimension of the swarm problem, denoted by N. The employed bees start the foraging process by randomly searching an initial position ${\mathbf{x}}_i^{\left(0\right)}$ in the domain D:

$${\mathbf{x}}_i^{\left(0\right)}=\left(\begin{array}{c}l{b}_1+{\varphi }_{i,1}^{\left(0\right)}·(u{b}_1-l{b}_1)\\ l{b}_2+{\varphi }_{i,2}^{\left(0\right)}·(u{b}_2-l{b}_2)\\ ⋮\\ l{b}_d+{\varphi }_{i,d}^{\left(0\right)}·(u{b}_d-l{b}_d)\end{array}\right)\in D,$$

(26)

where ${\varphi }_{i,\overline{1,d}}^{\left(0\right)}\in [-1,1]$ are random numbers. After this initialization step, the first Best Solution appears.

Each employed bee searches a new food source based on the location of the current food source ${\mathbf{x}}_i^{\left(k\right)}$ and another food source ${\mathbf{x}}_j^{\left(k\right)}$ randomly selected. The new possible position is:

$${\overline{\mathbf{x}}}_i^{\left(k\right)}=\mathtt{sat}\left({\mathbf{x}}_i^{\left(k\right)}+\varphi \odot ({\mathbf{x}}_i^{\left(k\right)}-{\mathbf{x}}_j^{\left(k\right)})\right),$$

(27)

where $\varphi \in {[-1,1]}^d$ is an array of random numbers, ⊙ is the element-wise multiplication, and sat is the saturation function that does not allow the position to be outside the searching domain D. Now, the position of the ith employed bee for the next iteration will be:

$${\mathbf{x}}_i^{(k+1)}=argmin\{f\left({\mathbf{x}}_i^{\left(k\right)}\right),f\left({\overline{\mathbf{x}}}_i^{\left(k\right)}\right)\},$$

(28)

which means that an employed bee will never choose a source with less nectar. If the position for the next iteration will not be changed, the abandonment counter of the ith employed bee increments.

The onlooker bees use the information shared by each employed bee and choose a location around the position of the employed bees. The fitness value of a solution ${\mathbf{x}}_i^{\left(k\right)}$ is given by:

$${\displaystyle logW\left(i\right)=\frac{-f\left({\mathbf{x}}_i^{\left(k\right)}\right)}{\frac{1}{N}{\sum }_{j}f\left({\mathbf{x}}_j^{\left(k\right)}\right)},}$$

(29)

and the probability that ith bee’s source will be selected by an onlooker bee is:

$${\displaystyle p_i=\frac{W\left(i\right)}{{\sum }_{j}W\left(j\right)}.}$$

(30)

Now, using a roulette wheel selection method, the onlooker bee will choose a source i and, using the same searching technique as an employed bee, a new position is computed using (27). If the outlooker bee founds a better solution than the ith employed bee, they change their roles, otherwise the abandonment counter for the ith source increments. After this step, we have N employed bees and N unemployed bees. However, if the abandonment counter for the ith source exceeds a threshold, the ith employed bee becomes a scout and tries to find a new location using the relation (26).

After every loop corresponding to employed, outlooker and scout bees, it is checked if there is a food source with a solution better than the last one. The algorithm is over when the maximum number of cycles is reached or when there is no improvement of the Best Solution after a prescribed number of cycles.

2.4. Proposed Method

The solution of the problem described in (21) using the classical D–K iteration with ${\mathcal{H}}_2/{\mathcal{H}}_{\infty }$ framework leads to a high-order controller. As a solution to this issue, a fixed structure family of controllers $\mathcal{K}$ can be considered, and the problem (21) becomes:

$$\underset{\begin{array}{c}K\in \mathcal{K}\\ K\mathrm{stab}\end{array}}{inf}\underset{\omega \in {\mathbb{R}}_{+}}{sup}\underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(D\left(j\omega \right)·\mathtt{LLFT}(P,K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right).$$

(31)

The new problem has the disadvantage of being non-convex in terms of the parameters of the proposed controller structure. However, the problem described in (31) will also be solved using a D–K iterative procedure, but the non-convex part of the problem has fewer degrees of freedom and it will be managed with an ABC optimization.

Starting with the model of the process G, a plant P is obtained after the augmentation step:

$$P:\left(\begin{array}{c}\dot{\mathbf{x}}\left(t\right)\\ \mathbf{z}\left(t\right)\\ \mathbf{y}\left(t\right)\end{array}\right)=\left(\begin{array}{ccc}A& B_w& B_u\\ C_z& D_{zw}& D_{zu}\\ C_y& D_{yw}& D_{yu}\end{array}\right)\left(\begin{array}{c}\mathbf{x}\left(t\right)\\ \mathbf{w}\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right),$$

(32)

where the meaning of the inputs and outputs remains the same as in Section 2.2. Considering the advantages of the fractional-order controllers, the proposed structure of the controller is:

$$K^{\theta }\left(s\right)=\left(\begin{array}{cccc}{C}_{FO}^{(1,1)}\left(s\right)& C_{FO}^{(1,2)}\left(s\right)& \cdots & C_{FO}^{(1,n_y)}\left(s\right)\\ C_{FO}^{(2,1)}\left(s\right)& C_{FO}^{(2,2)}\left(s\right)& \cdots & C_{FO}^{(2,n_y)}\left(s\right)\\ ⋮& ⋮& \ddots & ⋮\\ C_{FO}^{(n_u,1)}\left(s\right)& C_{FO}^{(n_u,2)}\left(s\right)& \cdots & C_{FO}^{(n_u,n_y)}\left(s\right)\end{array}\right),$$

(33)

where $C_{FO}^{(i,j)}$ is a fractional order controller from the ith input to the jth output and has the form:

$$C_{FO}^{(i,j)}\left(s\right)=K_P^{(i,j)}+\frac{K_i^{(i,j)}}{s^{{\lambda }^{(i,j)}}}+K_D^{(i,j)}{s}^{{\mu }^{(i,j)}},$$

(34)

with the tunable parameters described by the vector:

$${\theta }^{(i,j)}={\left(\begin{array}{ccccc}{K}_P^{(i,j)}& K_i^{(i,j)}& {\lambda }^{(i,j)}& K_D^{(i,j)}& {\mu }^{(i,j)}\end{array}\right)}^{\top }\in D_{ABC}^{(i,j)}\subset {\mathbb{R}}^5.$$

(35)

Using these considerations, the desired family of the fixed structure controllers can be described as follows:

$$\mathcal{K}=\left\{K^{\theta }\left(s\right)\left|\theta \in D_{ABC}\equiv D_{ABC}^{(1,1)}\times D_{ABC}^{(1,2)}\times \cdots \times D_{ABC}^{(n_u,n_y)}\right.\right\},$$

(36)

where all parameters are stored in a single vector $\theta$ describing all degrees of freedom of the tunable controller. Using the ORA mechanism, each component $K^{\theta }\left(s\right)\in \mathcal{K}$ has a state-space representation:

$$K^{\theta }:\left(\begin{array}{c}{\dot{\mathbf{x}}}_K\left(t\right)\\ \mathbf{u}\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{A}_K\left(\theta \right)& B_K\left(\theta \right)\\ C_K\left(\theta \right)& D_K\left(\theta \right)\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_K\left(t\right)\\ \mathbf{y}\left(t\right)\end{array}\right)$$

(37)

Next, we denote by $P_o^{\theta }$ the closed-loop system represented by the lower linear fractional transformation between the augmented plant P and controller $K^{\theta }$, which can be represented as:

$$P_o^{\theta }=\mathtt{LLFT}(P,K^{\theta }):\left(\begin{array}{c}{\dot{\mathbf{x}}}_o\left(t\right)\\ {\mathbf{z}}_o\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{A}_o\left(\theta \right)& B_o\left(\theta \right)\\ C_o\left(\theta \right)& D_o\left(\theta \right)\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_o\left(t\right)\\ {\mathbf{w}}_o\left(t\right)\end{array}\right),$$

(38)

where state vector, input vector and output vector of the closed-loop system are:

$${\mathbf{x}}_o=\left(\begin{array}{c}\mathbf{x}\\ {\mathbf{x}}_K\end{array}\right),{\mathbf{w}}_o\equiv \mathbf{w}\mathrm{and}{\mathbf{z}}_o\equiv \mathbf{z}.$$

(39)

In Algorithm 1 the main steps necessary to obtain the parameters of the controller having the structure (33) are presented. The inputs of the algorithm are the closed-loop plant $P_o^{\theta }$, containing the tunable controller parameters $\theta$, and the parameters $\alpha$ and $\beta$ which describe the cost function that needs to be minimized using an ABC optimization, but, $P_o^{\theta }$ must also contain the varying parameters and unmodelled dynamics. Therefore, the plant $P_{\mathsf{\Delta }}$ obtained after the augmentation step with uncertainties $\mathsf{\Delta }$ has the form presented in (16). Thus, the first step made in Algorithm 1 consists in transforming the closed-loop system $P_o^{\theta }$ in the generalized closed-loop system $P_{o,\mathsf{\Delta }}^{\theta }$ described as follows:

$$P_{o,\mathsf{\Delta }}^{\theta }=\mathtt{LLFT}(P_{\mathsf{\Delta }},K^{\theta }):\left(\begin{array}{c}{\dot{\mathbf{x}}}_o\left(t\right)\\ {\overline{\mathbf{z}}}_o\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{\overline{A}}_o\left(\theta \right)& {\overline{B}}_o\left(\theta \right)\\ {\overline{C}}_o\left(\theta \right)& {\overline{D}}_o\left(\theta \right)\end{array}\right)\left(\begin{array}{c}{\mathbf{x}}_o\left(t\right)\\ {\overline{\mathbf{w}}}_o\left(t\right)\end{array}\right),$$

(40)

where the new extended performance inputs and outputs are:

$${\overline{\mathbf{w}}}_o=\left(\begin{array}{c}\mathbf{d}\\ {\mathbf{w}}_o\end{array}\right){\overline{\mathbf{z}}}_o=\left(\begin{array}{c}\mathbf{v}\\ {\mathbf{z}}_o\end{array}\right)$$

(41)

As mentioned in Section 2.2, the structured singular value can be bounded using two D-scaling factors, one for the left and one for the right scaling, denoted $D_L$ and $D_R$, respectively. As it can be notice in Figure 2, a new closed-loop plant ${\overline{P}}_o^{\theta }$ is obtained:

$${\overline{P}}_o^{\left(\theta \right)}\left(s\right)=D_L\left(s\right)·P_{o,\mathsf{\Delta }}^{\theta }\left(s\right)·D_R^{-1}\left(s\right),$$

(42)

having the state-space representation:

$${\overline{P}}_o^{\theta }:\left(\begin{array}{c}{\dot{\widehat{\mathbf{x}}}}_o\left(t\right)\\ {\widehat{\mathbf{z}}}_o\left(t\right)\end{array}\right)=\left(\begin{array}{cc}{\widehat{A}}_o\left(\theta \right)& {\widehat{B}}_o\left(\theta \right)\\ {\widehat{C}}_o\left(\theta \right)& {\widehat{D}}_o\left(\theta \right)\end{array}\right)\left(\begin{array}{c}{\widehat{\mathbf{x}}}_o\left(t\right)\\ {\widehat{\mathbf{w}}}_o\left(t\right)\end{array}\right).$$

(43)

All the above mentioned plants are presented in Figure 2. First, the generalized plant $P_{\mathsf{\Delta }}$ has a LLFT connection with the controller $K^{\theta }$, resulting the closed-loop plant $P_{o,\mathsf{\Delta }}^{\theta }$, having the input vector ${\overline{\mathbf{w}}}_o$ and the output vector ${\overline{\mathbf{z}}}_o$. After the D-scaling step, a new plant ${\overline{P}}_o^{\theta }$ is obtained, having the input vector ${\widehat{\mathbf{w}}}_o$ and the output vector ${\widehat{\mathbf{z}}}_o$.

Before starting the while loop of Algorithm 1, an initialization of the generalized closed-loop plant with D-scale ${\overline{P}}_o^{\theta }$ is performed with the initial scale factors $D_L=I_{n_w}$ and $D_R=I_{n_z}$, as seen in line 2. As noticed in line 4, the K step is performed using this generalized plant having as degrees of freedom the tunable parameters $\theta$. In order to compute the controller parameters ${\theta }^{*}$, the ABC optimization will be used. The cost function to be minimized is:

$$f:D_{ABC}\to {\mathbb{R}}_{+},f\left(\theta \right)=\left\{\begin{array}{cc}\left|\right|{\overline{P}}_o^{\theta }|{|}_{\infty },\hfill & \mathrm{if}{\overline{P}}_o^{\theta }\mathrm{is}\mathrm{stable}\hfill \\ \alpha {\lambda }_{max}\left({\widehat{A}}_o\right)+\beta ,\hfill & \mathrm{if}{\overline{P}}_o^{\theta }\mathrm{is}\mathrm{unstable}\hfill \end{array}\right.,$$

(44)

where the operator ${\lambda }_{max}$ is defined by:

$${\lambda }_{max}:{\mathbb{R}}^{M\times M}\to \mathbb{R},{\lambda }_{max}\left(A\right)=max\left\{\mathrm{Re}\left(\lambda \right)\right|\lambda \in \mathsf{\Lambda }\left(A\right)\}.$$

(45)

Algorithm 1: Fixed Structure $\mu$-Synthesis.

The procedure computeKstep used to obtain the controller parameters is briefly presented in Algorithm 2 and will be described below. The inputs of the algorithm are the generalized closed-loop plant with D-scaling step ${\overline{P}}_o^{\theta }$, along with the parameters $\alpha$ and $\beta$ which describe the cost function (44). The first step of this routine consists in computing the domain $D_{ABC}$ of the cost function based on the constraints of the tunable parameters. The main limitation is represented by the fractional orders of the integral and derivative effects of the controller which must remain in $(0,1)$, according to ORA.

In the second step of routine Algorithm 2, an initial population is created. Let N be the dimension of the swarm problem. This parameter can be given as input, but as a good practice, this can be chosen 100 times larger than the number of tunable parameters. The initialization step consists in randomly generating the positions ${\theta }_1^{\left(0\right)},{\theta }_2^{\left(0\right)},\cdots ,{\theta }_N^{\left(0\right)}$ of the food sources for each employed bee in the domain $D_{ABC}$ using relation (26). After the initialization step, the first Best Solution (BS) is computed as:

$${\theta }^{★}=argmin\left\{f\left({\theta }_i^{\left(k\right)}\right)\left|i=\overline{1,N}\right.\right\}\Rightarrow BS=f\left({\theta }^{★}\right).$$

(46)

Using this initial population, the main while loop starts. In line 4, the employed bees step is performed. In this step, each employed bee searches a new position around using relation (27), resulting N new possible positions for the next step, denoted by ${\overline{\theta }}_1^{\left(k\right)},{\overline{\theta }}_2^{\left(k\right)},\cdots ,{\overline{\theta }}_N^{\left(k\right)}$, and the proposed positions of the employed bees at step $k+1$ will be:

$${\widehat{\theta }}_i^{(k+1)}=argmin\{f\left({\theta }_i^{\left(k\right)}\right),f\left({\overline{\theta }}_i^{\left(k\right)}\right)\},i=\overline{1,N}.$$

(47)

If, for a specific food source i, the proposed position coincides with the previous position, the abandonment counter increments.

Using the proposed solutions of the above step, each onlooker bee selects a new possible solution based on the roulette wheel selection mechanism and relation (27), resulting a new set of proposed solutions: ${\tilde{\theta }}_1^{\left(k\right)},{\tilde{\theta }}_1^{\left(k\right)},\cdots ,{\tilde{\theta }}_N^{\left(k\right)}$. If the ith onlooker bee has a better solution than ith employed bee, they exchange their roles, otherwise the abandonment counter for the ith food source increments. After this step, the new set of proposed positions for the employed bees are:

$${\widehat{\theta }}_i^{(k+1)}=argmin\{f\left({\tilde{\theta }}_i^{\left(k\right)}\right),f\left({\widehat{\theta }}_i^{(k+1)}\right)\},i=\overline{1,N}.$$

(48)

Algorithm 2: Compute K Step.

After performing the outlooker bees step from line 5, the abandonment counter for each active food source will be interrogated. If the abandonment counter of the ith position exceeds a prescribed threshold, denoted by LIMIT, the employed bee becomes a scout bee and its new position is obtained using (26). As a good practice, this paper proposed an improved mechanism of converting employed bees in scout bees. If the value of the cost function at the proposed position $f\left({\widehat{\theta }}_i^{(k+1)}\right)$ is over $\beta$, this means that, in the case of a good calibration of parameters $\alpha$ and $\beta$, the solution corresponds to an unstable closed-loop system and can be dropped. The positions of the employed bees for the next iteration become ${\theta }_1^{(k+1)},{\theta }_2^{(k+1)},\cdots ,{\theta }_N^{(k+1)}$.

The last step of the main while loop, presented in line 7, consists in computing the Best Solution after this new iteration:

$${\theta }^{★}=argmin\left\{f\left({\theta }_{old}^{★}\right),f\left({\theta }_1^{\left(k\right)}\right),f\left({\theta }_2^{\left(k\right)}\right),\cdots ,f\left({\theta }_N^{\left(k\right)}\right)\right\}\Rightarrow BS=f\left({\theta }^{★}\right),$$

(49)

and, then, checking if there exist any improvements after this step. In order to have a good trade-off between execution time and solution accuracy, it can be useful to establish a threshold for the number of steps when no improvement appears in Best Solution and mark if there exists such an improvement. Being a metaheuristic optimization algorithm, the runtime can be made deterministic and theoretically finite by imposing the variable $MNC$, without convergence guarantees. In practice, such methods have been successfully employed in various global optimization problems and, being that it is an offline optimization, it is not problematic for our approach.

Returning to Algorithm 1, after the K step is performed, a new parameter vector ${\theta }^{★}$ results, which leads to a new generalized closed-loop plant ${\overline{P}}_o^{{\theta }^{★}}$ and to a new uncertain plant $P_{o,\mathsf{\Delta }}^{{\theta }^{★}}$. The closed-loop plant $P_{o,\mathsf{\Delta }}^{{\theta }^{★}}$ is used to compute the next D-scale factors. From $P_{o,\mathsf{\Delta }}^{{\theta }^{★}}$ we extract a nominal plant $P_{o,nom}^{{\theta }^{★}}$ by fixing the tunable parameters of the controller with the values determined in step 4.

The computeDstep routine from line 7 receives as input this nominal plant and returns the left and the right D-scale factors. Based on the poles and the transmission zeros of the nominal plant $P_{o,nom}^{{\theta }^{★}}$, a set $\mathsf{\Omega }=\{{\omega }_1,{\omega }_2,\cdots ,{\omega }_F\}$ of frequencies is generated. Then, we need to get the frequency response data for each scaling factor by solving the following generalized eigenvalue problem:

$$\begin{array}{cc}\hfill & min\gamma ,\hfill \\ \hfill \mathrm{s}.\mathrm{t}.\overline{\sigma }& \left(D_L\left(j{\omega }_i\right)·P_{o,nom}^{{\theta }^{★}}\left(j{\omega }_i\right)·D_R^{-1}\left(j{\omega }_i\right)\right)<\gamma ,\hfill \end{array}$$

(50)

for each $i=\overline{1,F}$, which is nothing but a Parrot problem which can be solved point by point using LMI techniques, as  mentioned in Section 2.2. Once the frequency response data points are obtained for each value in $\mathsf{\Omega }$, we need to fit two minimum phase systems, one for each scaling factor, then perform the D-scaling step, giving a new generalized closed-loop system:

$${\overline{P}}_o^{{\theta }^{★}}\left(s\right)=D_L\left(s\right)·P_{o,\mathsf{\Delta }}^{★}\left(s\right)·D_R^{-1}\left(s\right).$$

(51)

In a similar manner with the computeKstep routine, there are possible stopping criteria which can be used. First, a threshold for the maximum number of D–K iterations can be imposed. Another important stopping condition appears if the upper bound of the structural singular value is less than 1, because this fact already guarantees that the controller ensures robust stability and robust performance. In accordance with the allowed maximum number of steps, a stopping criterion could be to check if there are any improvements after a certain number of steps.

3. Numerical Results

In this section we illustrate how the proposed method can be used on a benchmark problem. The process is represented by a Computer Numerical Control (CNC) machine with two orthogonal axis which are operated by two servo DC motors. A Trio Motion Coordinator family controller was used for the CNC motors. The programming language was Trio Basic, which provides various functions such as linear, circular and helical interpolation, variable speed and acceleration profile functions and control functions to ensure smooth and synchronized motions.

A mathematical model for each axis was determined on the basis of measured data: angular speeds ${\omega }_x$ and ${\omega }_y$, and angular positions ${\theta }_x$ and ${\theta }_y$. The state space mathematical model of the machine is described as follows:

$$G:\left(\begin{array}{c}{\dot{\omega }}_x\\ {\dot{\theta }}_x\\ {\dot{\omega }}_y\\ {\dot{\theta }}_y\\ {\theta }_x\\ {\theta }_y\end{array}\right)=\left(\begin{array}{cccccc}-\frac{1}{T_{mx}}& 0& 0& 0& \frac{K_{mx}}{T_{mx}}& K_{xy}\\ 1& 0& 0& 0& 0& 0\\ 0& 0& -\frac{1}{T_{my}}& 0& K_{yx}& \frac{K_{my}}{T_{my}}\\ 0& 0& 1& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\end{array}\right)\left(\begin{array}{c}{\omega }_x\\ {\theta }_x\\ {\omega }_y\\ {\theta }_y\\ u_x\\ u_y\end{array}\right),$$

(52)

where the state vector is $\mathbf{x}={\left(\begin{array}{cccc}{\omega }_x& {\theta }_x& {\omega }_y& {\theta }_y\end{array}\right)}^{\top }$, the input vector is $\mathbf{u}={\left(\begin{array}{cc}{u}_x& u_y\end{array}\right)}^{\top }$ and the output vector is $\mathbf{y}={\left(\begin{array}{cc}{\theta }_x& {\theta }_y\end{array}\right)}^{\top }$. The model parameters, along with their nominal values and uncertainty range are detailed in

Table 1. Nominal values and uncertainty ranges of the CNC model parameters.

Parameter	Nominal Value	Uncertainty Range	Parameter	Nominal Value	Uncertainty Range
$T_{mx}$	$0.02448$	$\pm 10\%$	$T_{my}$	$0.01139$	$\pm 10\%$
$K_{mx}$	$25.8017$	$\pm 10\%$	$K_{my}$	$25.1494$	$\pm 10\%$
$K_{xy}$	$26.65$	$\pm 10\%$	$K_{yx}$	$24.46$	$\pm 10\%$

.

For this control problem a mixed-sensitivity loop shaping technique is used, which provides a good trade-off between performance and robustness. In order to use this technique, a new plant model will be obtained after the augmentation process, as in Figure 3. The performance inputs for the resulting augmented plant P are the references for both axis $\mathbf{w}\equiv \mathbf{r}={\left(\begin{array}{cc}{r}_x& r_y\end{array}\right)}^{\top }$. For the augmentation procedure, the closed-loop transfer functions need to be weighted from the reference signals $\mathbf{r}$ to their corresponding error signals $\mathbf{e}$, output signals $\mathbf{y}$ and command signals $\mathbf{u}$, which are named: sensitivity function $S={(I+GK)}^{-1}$, complementary sensitivity function $T=I-S$ and control effort $KS$, respectively. The performance output vector is composed from the weighted outputs of these three vector-valued functions:

$$\mathbf{z}=\left(\begin{array}{c}{\mathbf{z}}_S\\ {\mathbf{z}}_T\\ {\mathbf{z}}_{KS}\end{array}\right),\mathrm{where}{\mathbf{z}}_S=\left(\begin{array}{c}{z}_{S,x}\\ z_{S,y}\end{array}\right),{\mathbf{z}}_T=\left(\begin{array}{c}{z}_{T,x}\\ z_{T,y}\end{array}\right)\mathrm{and}{\mathbf{z}}_{KS}=\left(\begin{array}{c}{z}_{KS,x}\\ z_{KS,y}\end{array}\right).$$

(53)

In order to ensure good disturbance rejection, good reference tracking and stabilization of an unstable plant, the sensitivity function must have small magnitude, which means that the magnitude of the open loop must be large. On the other hand, in order to ensure mitigation of measurement noise and robust stability, the complementary sensitivity function must have small magnitudes, which implies that the magnitude of the open loop system must be small. More than that, the command signal must have small magnitude, which implies that the control effort transfer function must also have small magnitude. However, although these requirements seem to be conflicting, the frequency range where the magnitude of the open loop must be small and high are mostly disjunctive: the magnitude should be high for low frequencies and should be low for high frequencies. In order to ensure the desired shape of these three functions, three weighting functions must be designed, respectively.

The sensitivity function is a very good indicator of the closed-loop system tracking performance and has the advantage of being sufficient to consider only the magnitude. Typical specifications for sensitivity weighting functions are: minimum bandwidth frequency ${\omega }_B^{★}$, maximum steady-state error A and maximum peak magnitude M, imposed by the following model [36]:

$$W_S\left(s\right)=\frac{\frac{1}{M}s+{\omega }_B^{★}}{s+{\omega }_B^{★}A}.$$

(54)

In a similar manner, the weighting function for the complementary sensitivity must be designed using the following specifications: the maximum peak amplitude $M_T$, the maximum value for high frequencies $A_T$, the minimum bandwidth ${\omega }_{BT}^{★}$ and the roll-off n, formulated as:

$$W_T\left(s\right)=\frac{{\left(s+{\omega }_{BT}^{★}\right)}^n}{{\left(A_T^{1/n}s+{\omega }_{BT}^{★}{M}_T^{1/n}\right)}^n}.$$

(55)

For the control effort weighting function, the main performance specifications are the maximum value of the magnitude at low and high frequencies, denoted by $DC$ and $HF$, respectively, and an intermediate point of interest. Sometimes, the main goal is simply to maintain the command signal under a prescribed value due to physical limitations of the system or other causes.

The major advantage of this approach consists in sculpting the relevant loop functions to impose performances implying good tracking and dynamic behaviour. Of great use are the rise time limitation through ${\omega }_B$, steady-state error through A, while the roll-off slope of the closed-loop system imposed using n is directly coupled with sensor noise characteristics. These performances are specified for different frequency ranges, using the adequately selected weighting functions presented above.

Being a MIMO system with two inputs and two outputs, all weighting functions must be described by $2\times 2$ transfer matrices, but, following a standard decoupling procedure, the weighting functions will be $2\times 2$ diagonal transfer matrices. For the sensitivity, we consider two nearly similar weighting transfer functions, one for each axis, having the maximum bandwidth ${\omega }_{B,x}^{★}=3\left[\mathrm{rad}/\mathrm{s}\right]$ and ${\omega }_{B,y}^{★}=5\left[\mathrm{rad}/\mathrm{s}\right]$, the maximum steady-state error $A_x=A_y={10}^{-2}$ and the desired maximum sensitivity peak $M_x=M_y=1.5$, resulting in:

$$W_S\left(s\right)=\left(\begin{array}{cc}{W}_{S,x}\left(s\right)& 0\\ 0& W_{S,y}\left(s\right)\end{array}\right),\mathrm{where}{W}_{S,x}\left(s\right)=\frac{0.6667s+3}{s+0.03}\mathrm{and}{W}_{S,y}\left(s\right)=\frac{0.6667s+5}{s+0.05}.$$

(56)

For the complementary sensitivity weighting function, the same $2\times 2$ diagonal structure approach will be used, imposing the same parameters for both axis: maximum complementary bandwidth ${\omega }_{BT,x}^{★}=50$ [rad/s], maximum peak magnitude $M_{T,x}=M_{T,y}=1.5$, maximum magnitude at high frequencies $A_{T,x}=A_{T,y}={10}^{-2}$ and roll-off $n_x=n_y=1$, resulting in:

$$W_T\left(s\right)=\left(\begin{array}{cc}{W}_{T,x}\left(s\right)& 0\\ 0& W_{T,y}\left(s\right)\end{array}\right),\mathrm{where}{W}_{T,x}\left(s\right)=W_{T,y}\left(s\right)=\frac{s+50}{0.01s+75},$$

(57)

while the control effort weighting function being designed to encompass only the physical limitation of the command signal (between $-1$ and 1), resulting the transfer matrix:

$$W_{KS}\left(s\right)=\left(\begin{array}{cc}{W}_{KS,x}\left(s\right)& 0\\ 0& W_{KS,y}\left(s\right)\end{array}\right),\mathrm{where}{W}_{KS,x}\left(s\right)=W_{KS,y}\left(s\right)=1.$$

(58)

The proposed structure of the controller is a decentralized one with two $P{I}^{\lambda }{D}^{\mu }$ controllers:

$$\mathcal{K}=\left\{K^{\theta }\left(s\right)=\left(\begin{array}{cc}{K}_{P,x}+K_{I,x}{s}^{-{\lambda }_x}+K_{D,x}{s}^{{\mu }_x}& 0\\ 0& K_{P,y}+K_{I,y}{s}^{-{\lambda }_y}+K_{D,y}{s}^{{\mu }_y}\end{array}\right)\left|\theta \in D\subset {\mathbb{R}}^{10}\right.\right\},$$

(59)

where the tunable parameters are:

$$\theta ={\left(\begin{array}{cccccccccc}{K}_{P,x}& K_{I,x}& {\lambda }_x& K_{D,x}& {\mu }_x& K_{P,y}& K_{I,y}& {\lambda }_y& K_{D,y}& {\mu }_y\end{array}\right)}^{\top }.$$

(60)

The problem to be solved with the proposed method is the mixed-sensitivity fixed structure $\mu$-synthesis one, described as:

$$\begin{array}{cc}\hfill \underset{\begin{array}{c}K\in \mathcal{K}\\ K\mathrm{stab}\end{array}}{min}& \underset{\omega \in {\mathbb{R}}_{+}}{sup}\underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(D\left(j\omega \right)·\mathtt{LLFT}(P,K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right),\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.{\left|\left|\begin{array}{c}{W}_{S}S\\ W_{T}T\\ W_{KS}KS\end{array}\right|\right|}_{\infty }<1.\hfill \end{array}$$

(61)

The settings for computeKstep used in the experiments are: the swarm dimension $N=1000$, the maximum number of cycles $MNC=50$, the maximum number of cycles with no improvements $NOIMP=10$, the limit for the abandonment counter $LIMIT=10$. The parameters necessary to describe the cost function (44) are $\alpha =1$ and $\beta ={10}^5$. The maximum number of D–K iterations is $MAX\_ITER=10$ and the maximum window length for assessing lack of progress is 4. Using this setup, the mixed-sensitivity fixed structure $\mu$-synthesis problem (61) is solved using four D–K iterations, as  noticed in

Table 2. The evolution of the structural singular value in the D–K iteration procedure used to solve the mixed-sensitivity fixed structure $\mu$-synthesis problem for the case study—FO-PID structure.

D–K Iteration Number	1	2	3	4
Number of ABC Iterations	50	39	50	45
Peak Value of$\mathit{\mu }$	$1.8956$	$1.0185$	$1.0021$	$0.9822$

. The fractional order controller has been approximated using ORA with the following parameters: the frequency range is $[{\omega }_l,{\omega }_u]$ = $[1\times {10}^{-4},1\times {10}^3]$$[\mathrm{rad}/\mathrm{s}]$ and the order of the approximation is $N=3$. Given that, the resulting controller is:

$$K^{{\theta }^{★}}\left(s\right)=\left(\begin{array}{cc}0.01+7.1116·s^{-0.7648}+0.1133·s^{0.0909}& 0\\ 0& 0.1464+10·s^{-0.9926}+0.1344·s^{0.0549}\end{array}\right).$$

(62)

The imposed shape of the sensitivity, complementary sensitivity and control effort functions for both axis are depicted in Figure 4, along with the obtained shapes of those functions with the resulting controller for 100 Monte Carlo simulations. It can be noticed that all resulting Bode diagrams are under the prescribed shapes, which is guaranteed by the fact that the upper bound of the structured singular value ${\mu }_{\mathsf{\Delta }}\left(P_{o,\mathsf{\Delta }}^{{\theta }^{★}}\right)$ is less than 1. The time-domain performances of the systems $r_x\to {\theta }_x$ and $r_y\to {\theta }_y$ are depicted in Figure 5, which are correlated with the frequency-domain performances. As such, the minimum values of the bandwidths for both x and y axis are over ${\omega }_{B,x}^{★}=3$ [rad/s] and ${\omega }_{B,y}^{★}=5$ [rad/s], which means that the rise time is less than $0.33\left[s\right]$ and $0.2\left[s\right]$, respectively. For the system $r_x\to {\theta }_x$, the rise time of the nominal system is $0.248\left[s\right]$ and varies from $0.227\left[s\right]$ to $0.281\left[s\right]$, while the settling time is $0.558\left[s\right]$ and varies from $0.496\left[s\right]$ to $0.664\left[s\right]$. On the other hand, for the system $r_y\to {\theta }_y$, the rise time of the nominal system is $0.211\left[s\right]$ and varies from $0.19\left[s\right]$ to $0.235\left[s\right]$, having the settling time for the nominal system $0.405\left[s\right]$ and varying from $0.363\left[s\right]$ to $0.454\left[s\right]$. According to the shape of the actual obtained sensitivity functions, there is no overshoot and no steady-state error for neither of the experiments presented using Monte Carlo simulations. Moreover, the reciprocal axis influence is small, as resulted from numerical simulations, where the peak amplitude from $r_y$ to ${\theta }_x$ varies from $0.01$ to $0.0143$, and the peak amplitude from $r_x$ to ${\theta }_y$ varies from $3.55\times {10}^{-3}$ to $4.97\times {10}^{-3}$, respectively.

4. Discussion

As noticed, the resulting controller contains a small fractional D term, which may be negligible. As such, the controller may be resynthesized with an imposed diagonal $P{I}^{\lambda }$ structure:

$$\mathcal{K}=\left\{K^{\theta }\left(s\right)=\left(\begin{array}{cc}{K}_{P,x}+K_{I,x}{s}^{-{\lambda }_x}& 0\\ 0& K_{P,y}+K_{I,y}{s}^{-{\lambda }_y}\end{array}\right)\left|\theta \in D\subset {\mathbb{R}}^6\right.\right\},$$

(63)

where the tunable parameters are:

$$\theta ={\left(\begin{array}{cccccc}{K}_{P,x}& K_{I,x}& {\lambda }_x& K_{P,y}& K_{I,y}& {\lambda }_y\end{array}\right)}^{\top }.$$

(64)

The control problem remains the same as in (61), maintaining the constraints as in (56)–(58). The settings for this experiment were kept the same as for the previous case: $N=1000$, $MNC=50$, $NOIMP=10$, $LIMIT=10$, $\alpha =1$, $\beta ={10}^5$, $MAX\_ITER=10$, $[{\omega }_l,{\omega }_u]=[1\mathrm{e}-4,1\mathrm{e}3][\mathrm{rad}/\mathrm{s}]$ and $N=3$. Using this setup, the mixed-sensitivity fixed structure $\mu$-synthesis problem (61) is solved using four D–K iterations, as  noticed in

Table 3. The evolution of the structural singular value in the D–K iteration procedure used to solve the mixed-sensitivity fixed structure $\mu$-synthesis problem for the case study—FO-PI structure.

D–K Iteration Number	1	2	3	4
Number of ABC Iterations	30	30	24	30
Peak Value of$\mathit{\mu }$	$3.0676$	$1.1961$	$1.0026$	$0.9959$

. The resulting controller is:

$$K_{FO-PI}^{{\theta }^{★}}\left(s\right)=\left(\begin{array}{cc}0.2229+5.9392·s^{-0.7792}& 0\\ 0& 0.3949+10·s^{-0.9733}\end{array}\right).$$

(65)

The imposed s hape of the sensitivity, complementary sensitivity and control effort functions for both axis are depicted in Figure 6, along with the obtained shapes of those functions with the resulting controller for 100 Monte Carlo simulations. It can be noticed that all resulting Bode diagrams are under the prescribed shapes, which is guaranteed by the fact that the upper bound of the structured singular value ${\mu }_{\mathsf{\Delta }}\left(P_{o,\mathsf{\Delta }}^{{\theta }^{★}}\right)$ is less than 1.

The proposed method results are compared with those obtained using the musyn routine from MATLAB [37]. The musyn routine manages to solve both fixed structure and free-structure $\mu$-synthesis control problems. The first solution of the problem (61) is obtained using the proposed method with results already presented in Section 3. Starting from the same control problem, the structured $\mu$-synthesis version of the musyn routine with the same stopping criteria was used. The resulting controller after 2 iterations achieved its best performance ${\mu }_{\mathsf{\Delta }}\left(P_{o,\mathsf{\Delta }}^{{\theta }^{★}}\right)=0.9842$, which means that there is a mathematical guarantee for robust stability and robust performance. The transfer matrix of this controller is:

$$K^{{\theta }^{★}}\left(s\right)=\left(\begin{array}{cc}0.0101+7.4689·s^{-0.7860}+0.1729·s^{0.0542}& 0\\ 0& 0.0881+5.4840·s^{-0.9749}+0.2031·s^{0.0348}\end{array}\right).$$

(66)

In opposition with the above approaches, the controller obtained with $\mu$-synthesis procedure from MATLAB without imposing a fixed structure is of order 74 and, after 10 iterations, the best achieved performance is ${\mu }_{\mathsf{\Delta }}\left(P_{o,\mathsf{\Delta }}^{{\theta }^{★}}\right)=1.003$. The frequency-domain data for structured singular values corresponding to these three numerical simulations are presented in Figure 7. As mentioned, the peak value of the $\mu$ values for the unstructured problem is over the critical value 1, while the FO-PID controller manages to fulfill all requirements. A comparison between the frequency responses of the controllers is shown in Figure 8.

The integer-order approximation of a fractional element using ORA contains two exponential terms, as seen in (7), which presently cannot be treated using the realp objects in MATLAB. The actual solution used for this paper is to approximate the exponential function using Taylor series truncation. As stated in the Introduction, our approach requires only the range of the controller parameters and, thus, another object could be used instead of the entire structure of realp necessary in nonsmooth optimization-based algorithms used in MATLAB hinfstruct and systune routines. Based on this limitation, the toolbox presented in [38] will be extended in order to address the previously mentioned mathematical issues, due to the fact that only the fixed-structure ${\mathcal{H}}_{\infty }$ synthesis algorithm uses the realp object, which can be replaced by our approach, where the NP-hard non-convex problem is solved using an ABC swarm optimization.

Compared to other approaches previously described in the introduction, where LQR methods were considered for CNC machines, the advantages of the proposed method consist in the generality of the method and the flexibility of the robustness, with the possibility to compensate measurement noise, unmodelled dynamics and input-output disturbances. Therefore, while LQR requires the complete state vector to be measured at runtime, the proposed method requires only the provided measurements, modelled by the actual outputs of the system. Although the LQR method can be augmented with a state estimator in order to obtain output feedback, the main limitation of this approach is that the model of the plant must be accurate, while in the Robust Control framework, utilized in our method, an entire family of uncertain plants can be taken into consideration at the design phase. Moreover, it is difficult to impose exact limitations on the maximum allowed command signals, using the energy-based approach, which generally is an intrinsic limitation of the execution element.

5. Conclusions and Future Work

The current paper presents a new design method for fixed structure fractional-order controllers using the Robust Control framework. The proposed method manages to return nearly optimal parameters of a MIMO FO-PID as a solution for a mixed-sensitivity fixed structure $\mu$-synthesis control problem. Although the $\mu$-synthesis control problem is NP-hard, the D–K iteration algorithm represents a good approximation which allows to convert it into a P-hard problem. However, the returned controller is of high order, which means that an order reduction must be performed in order to implement the control law. Therefore, the imposed structure is an increasingly explored approach, although such a problem presents a non-convex component for the K step. Our approach consists of an artificial bee colony swarm optimization as a solution to this non-convex fixed structure ${\mathcal{H}}_{\infty }$ subproblem. This solution requires only the range of the controller parameters, as opposed to the nonsmooth optimization-based approach from MATLAB’s Robust Control Toolbox, where the parameters can be used only in polynomial structured expressions, which is an inherent limitation when fractional-order controllers are desired. Further, the case study of mixed-sensitivity $\mu$-synthesis position control problem for both axis of a CNC machine manages to underline the strong points of the method and of the imposed structure of the controller: it provides a mathematical guarantee for robust stability and robust performance, while the unstructured version of $\mu$-synthesis from MATLAB does not manage to offer it.

As future work, we want to include our method in a toolbox, as stated in the Discussion section, which starts from an initial process model and a set of desired performances, and manages to automatically obtain the augmented plant, the controller decoupling transfer matrix, the optimal values of the controller parameters, followed by closed-loop simulations and validations. Additionally, we propose to extend this technique for certain types of nonlinear systems.