µ -Synthesis for Fractional-Order Robust Controllers

: The current article presents a design procedure for obtaining robust multiple-input and multiple-output (MIMO) fractional-order controllers using a µ -synthesis design procedure with D – K iteration. µ -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 µ -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 µ -synthesis approach in MATLAB does not offer this.


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 H 2 and H ∞ norms defined in frequency domain as a performance measure. To solve H 2 /H ∞ 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 H 2 /H ∞ 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 µsynthesis approach extends the H ∞ optimization in order to obtain a robust controller for the uncertain plant which includes parametric and dynamic uncertainties [7,8]. µ-synthesis of the NP-hard µ-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 µ-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 µ-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 µ-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.

Materials and Methods
In this section we present a controller synthesis procedure which manages to find a fixed structure fractional-order controller using the µ-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.

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]: where Γ(α) : C + → C is the Euler Gamma function and the order of the integral operator is the complex parameter α ∈ C + . This extension develops a new area of research in the with t > 0 and α ∈ R + , having the Laplace transform [34]: where F(s) = L{ f (t)}(s). 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 λ ∈ R + of the integrator and the order µ ∈ R + of the differentiator, with the resulting transfer function: The time domain expression of the command signal c(t) can be expressed using the error signal ε(t) as: One of the major issues of such a fractional element, i.e., J λ c or J −µ c , 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: where the frequency values of the singularities are obtained based on the desired fractional order λ ∈ (0, 1), the integer order of the approximation N, along with the frequency range in which the approximation is valid [ω l , ω u ]. Using the following two coefficients: the above mentioned frequencies can be computed using: For the rest of the possible real values of λ, the approximation can be easily extended as: for λ ∈ (−1, 0) by inverting the relation (6), while for |λ| ≥ 1, the components could be the integer part [λ] and the fractional part {λ}, with the fractional part only approximated using (6).

Robust Control
The Robust Control framework assumes to minimize the H 2 /H ∞ norm of the lower linear fractional transformation (LLFT) of a plant P and a controller K: where the plant P can be written, in general form, as: where the signals involved in the above relation will be further detailed in a more general context. One approach to solving H 2 /H ∞ 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 δI, where δ is the maximum bound of the parameter for a physical variable, and unstructured, represented by a full block ∆ ∈ R m×m . The latter illustrates neglected or unknown dynamics uncertainties.
In the mixed-scenario case, the following set is considered: 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 ∈ C N×N the structured singular value with respect to the set ∆ is: if there exists ∆ ∈ ∆ such that the matrix I − M∆ is rank deficient, otherwise 0. For an LTI system described by the transfer matrix M(s) and an upper linear fractional transformation (ULFT) connection shown in Figure 1 (left), the structured singular value µ ∆ (M) can be defined as: Now, considering M(s) = LLFT(P, K)(s) 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: where three types of input signals are present-the command input u ∈ R n u , the performance input w ∈ R n w , and the disturbance input d ∈ R n d -and three types of outputs-the measurements vector y ∈ R n y , the performances vector z ∈ R n z , and the disturbance output v ∈ R n v . Besides the well-known H 2 /H ∞ methods, a controller that meets robust stability and robust performance alike can be computed using the µ-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 ∆, fulfills: sup Therefore, the robust control problem can be written as: 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]: where σ denotes the largest singular value, and the set D is defined as: Based on this upper bound, a good approximation of the initial non-convex problem can be employed by solving the following quasi-convex problem: If the scaling factor, represented by the system D, is fixed, then the problem (21) is nothing but a H ∞ 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 Ω = {ω 1 , . . . , ω N } using the following generalized eigenvalue problem: min γ, where from the solution X = (D(jω i )) * · D(jω i ), the matrix D(jω i ) 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(s). In summary, an iterative algorithm which solves the µ-synthesis problem starts by setting D = I, with the following steps applied successively: 1: Fix D and solve the H ∞ optimal problem to find a controller K: 2: Fix the controller K and solve the set of convex problems: for a given frequency range Ω and, then, fit a stable minimum phase transfer matrix D(s).
Steps 1 and 2 are executed in a loop sequence until the difference between two consecutive H ∞ norms is less than a prescribed tolerance or the maximum number of iterations is reached.

Artificial Bee Colony Optimization
The artificial bee colony (ABC) optimization is a nature-inspired algorithm used to minimize a cost function: 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 x (0) i in the domain D: . . .
where φ 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 x (k) i and another food source x (k) j randomly selected. The new possible position is: where φ ∈ [−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: 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 x (k) i is given by: and the probability that ith bee's source will be selected by an onlooker bee is: 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.

Proposed Method
The solution of the problem described in (21) using the classical D-K iteration with H 2 /H ∞ framework leads to a high-order controller. As a solution to this issue, a fixed structure family of controllers K can be considered, and the problem (21) becomes: 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: 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: where C FO is a fractional order controller from the ith input to the jth output and has the form: with the tunable parameters described by the vector: Using these considerations, the desired family of the fixed structure controllers can be described as follows: where all parameters are stored in a single vector θ describing all degrees of freedom of the tunable controller. Using the ORA mechanism, each component K θ (s) ∈ K has a state-space representation: Next, we denote by P θ o the closed-loop system represented by the lower linear fractional transformation between the augmented plant P and controller K θ , which can be represented as: where state vector, input vector and output vector of the closed-loop system are: 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 , containing the tunable controller parameters θ, and the parameters α and β which describe the cost function that needs to be minimized using an ABC optimization , but, P θ o must also contain the varying parameters and unmodelled dynamics. Therefore, the plant P ∆ obtained after the augmentation step with uncertainties ∆ has the form presented in (16). Thus, the first step made in Algorithm 1 consists in transforming the closed-loop system P θ o in the generalized closed-loop system P θ o,∆ described as follows: where the new extended performance inputs and outputs are: 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 P θ o is obtained: having the state-space representation: All the above mentioned plants are presented in Figure 2. First, the generalized plant P ∆ has a LLFT connection with the controller K θ , resulting the closed-loop plant P θ o,∆ , having the input vector w o and the output vector z o . After the D-scaling step, a new plant P  Before starting the while loop of Algorithm 1, an initialization of the generalized closedloop plant with D-scale P θ o 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 θ. In order to compute the controller parameters θ * , the ABC optimization will be used. The cost function to be minimized is: where the operator λ max is defined by: Algorithm 1: Fixed Structure µ-Synthesis.
Input: P θ o , α, β Output: K θ 1 get uncertain closed-loop plant P θ o,∆ as in (40) 2 set D L = I n w and D R = I n z and compute P check if improvement exists and increase N iter 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 P θ o , along with the parameters α and β 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 θ N 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: 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 θ N , and the proposed positions of the employed bees at step k + 1 will be: 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) 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 (θ (k+1) i ) is over β, this means that, in the case of a good calibration of parameters α and β, the solution corresponds to an unstable closed-loop system and can be dropped. The positions of the employed bees for the next iteration become θ The last step of the main while loop, presented in line 7, consists in computing the Best Solution after this new iteration: 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 θ results, which leads to a new generalized closed-loop plant P . . , ω 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: min γ, for each i = 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 Ω, 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: 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.

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 ω x and ω y , and angular positions θ x and θ y . The state space mathematical model of the machine is described as follows: where the state vector is x = ω x θ x ω y θ y , the input vector is u = u x u y and the output vector is y = θ x θ y . The model parameters, along with their nominal values and uncertainty range are detailed in Table 1. 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 w ≡ r = r x r y . For the augmentation procedure, the closed-loop transfer functions need to be weighted from the reference signals r to their corresponding error signals e, output signals y and command signals 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: 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 ω B , maximum steady-state error A and maximum peak magnitude M, imposed by the following model [36]: 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 ω BT and the roll-off n, formulated as: 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 ω 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 × 2 transfer matrices , but, following a standard decoupling procedure, the weighting functions will be 2 × 2 diagonal transfer matrices. For the sensitivity, we consider two nearly similar weighting transfer functions, one for each axis, having the maximum bandwidth ω B,x = 3 [rad/s] and ω B,y = 5 [rad/s], 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: , where W S,x (s) = 0.6667s + 3 s + 0.03 and W S,y (s) = 0.6667s + 5 s + 0.05 .
For the complementary sensitivity weighting function, the same 2 × 2 diagonal structure approach will be used, imposing the same parameters for both axis: maximum complementary bandwidth ω 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: 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: The proposed structure of the controller is a decentralized one with two PI λ D µ controllers: where the tunable parameters are: θ = K P,x K I,x λ x K D,x µ x K P,y K I,y λ y K D,y µ y .
The problem to be solved with the proposed method is the mixed-sensitivity fixed structure µ-synthesis one, described as: 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 NOI MP = 10, the limit for the abandonment counter LI MIT = 10. The parameters necessary to describe the cost function (44) are α = 1 and β = 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 µ-synthesis problem (61) is solved using four D-K iterations, as noticed in Table 2. The fractional order controller has been approximated using ORA with the following parameters: the frequency range is [ω l , ω u ] = [1 × 10 −4 , 1 × 10 3 ] [rad/s] and the order of the approximation is N = 3. Given that, the resulting controller is: 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 µ ∆ P θ o,∆ is less than 1. The time-domain performances of the systems r x → θ x and r y → θ y are depicted in Figure 5, which are correlated with the frequency-domain performances. . 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 θ x varies from 0.01 to 0.0143, and the peak amplitude from r x to θ y varies from 3.55 × 10 −3 to 4.97 × 10 −3 , respectively.

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 PI λ structure: where the tunable parameters are: θ = K P,x K I,x λ x K P,y K I,y λ y .
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, NOI MP = 10, LI MIT = 10, α = 1, β = 10 5 , MAX_ITER = 10, [ω l , ω u ] = [1e-4, 1e3] [rad/s] and N = 3. Using this setup, the mixed-sensitivity fixed structure µ-synthesis problem (61) is solved using four D-K iterations, as noticed in Table 3. The resulting controller is: Table 3. The evolution of the structural singular value in the D-K iteration procedure used to solve the mixed-sensitivity fixed structure µ-synthesis problem for the case study-FO-PI structure. 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 µ ∆ P θ o,∆ 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 µ-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 µ-synthesis version of the musyn routine with the same stopping criteria was used. The resulting controller after 2 iterations achieved its best performance µ ∆ (P θ o,∆ ) = 0.9842, which means that there is a mathematical guarantee for robust stability and robust performance. The transfer matrix of this controller is: In opposition with the above approaches, the controller obtained with µ-synthesis procedure from MATLAB without imposing a fixed structure is of order 74 and, after 10 iterations, the best achieved performance is µ ∆ (P θ o,∆ ) = 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 µ 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 H ∞ synthesis algorithm uses the realp object, which can be replaced by our approach, where the NP-hard nonconvex problem is solved using an ABC swarm optimization.

D-K Iteration Number
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.

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 µ-synthesis control problem. Although the µ-synthesis control problem is NPhard, 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 H ∞ 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 mixedsensitivity µ-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 µ-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.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: