μ-Synthesis FO-PID for Twin Rotor Aerodynamic System

Abstract

$\mu$-synthesis is a NP-hard optimization problem based on the generalized Robust Control framework which manages to find a controller which fulfills both robust stability and robust performance. In order to solve such problems, nonsmooth optimization techniques are employed to find nearly-optimal parameters values. However, the free parameters available for tuning must be involved only in classical arithmetic operations, which leads to a problem for the fractional-order operator or for its integer-order approximation, exponential operations being involved. The main goal of the current article consists of presenting a possibility to integrate a fixed-structure multiple-input-multiple-output (MIMO) fractional-order proportional-integral-derivative (FO-PID) controller in the $\mu$-synthesis optimization problem. The solution consists in a possibility to find a set of tunable parameters isomorphic with the fractional-order such that the coefficients involved in the approximation of the fractional element, along with the formulation of a fixed-structure mixed-sensitivity loop shaping $\mu$-synthesis control problem. The proposed design procedure is applied to a twin rotor aerodynamic system (TRAS) using both MATLAB numerical simulation and practical experiments on laboratory scale equipment. Moreover, a comparison with the unstructured $\mu$-synthesis is performed, highlighting the advantages of the proposed solution: simpler form and guaranteed robust stability and performance.

1. Introduction

One of the fundamental problems studied in Control Engineering concerns robustness, which characterizes the sensitivity of the closed loop system to the variation of plant parameters. One of the most used performance measures is the ${\mathcal{H}}_{\infty }$ norm. Starting from the approach of synthesizing a ${\mathcal{H}}_{\infty }$ controller by solving two Algebraic Riccati Equations (AREs) as in [1], a more numerically stable solution can be obtained using Popov triplets [2]. Alternatively, due to the limitations of this approach represented by the impossibility of solving singular problems, the AREs were replaced with Algebraic Riccati Inequalities (ARIs) and were solved using Linear Matrix Inequalities (LMIs) [3]. The last two approaches have been recently implemented in open-source manners in [4,5]. However, the classical ${\mathcal{H}}_{\infty }$ control problem manages to ensure nominal stability and nominal performance only. In order to consider dynamic and parametric uncertainties, the plant is formulated as an upper linear fractional transform with such an uncertainty block and the $\mu$-synthesis can be used for computing a robust controller based on the classical D–K iterations [6]. The major concern about these methods consists of the fact that the controller is usually of high order. However, imposing a fixed structure leads to a non-convex problem which cannot be approximated as in the case of $\mu$-synthesis. The solutions, initially proposed for ${\mathcal{H}}_{\infty }$ problem [7], and then for $\mu$-synthesis [8] as well, are based on nonsmooth optimization techniques. A CACSD toolbox that manages to offer an end-to-end solution for designing a robust controller starting from a given plant is presented in [9].

The most well-known controller structure, which is highly used in industry, is the proportional-integral-derivative (PID) regulator. Its form is generally given as an example for fixed-structure controllers and nonsmooth optimization methods were designed around it [10]. An extension with two extra degrees of freedom is represented by fractional-order PID (FO-PID), which improves the robustness of the closed loop system. As tuning methods, the well-known methods used for designing integer-order controllers were extended for fractional-order controllers as well. As such, two generalized versions of Kessler’s magnitude methods were presented in [11,12], while a fractional-order internal model controller with event-based implementation was developed in [13]. A fractional-order integrator was used as a model for the servo problem in [14], while the same structure was used as a speed controller for a DC motor in [15]. In [16] crone control methodologies were presented, along with LMI formulation for the ${\mathcal{H}}_{\infty }$ fractional-order control problem. An artificial bee colony optimization for a MIMO FO-PID controller design by solving the mixed-sensitivity $\mu$-synthesis control problem is presented in [17].

The twin rotor aerodynamic system (TRAS) is a well-known benchmark system used to illustrate the control methods designed in literature. A two degrees of freedom (2-DOF) discrete-time $\mu$-synthesis controller of order 24 was presented in [18]. A decentralized fixed-structure PID controller designed using ${\mathcal{H}}_{\infty }$ is presented in [19], along with a comparison between the full-order ${\mathcal{H}}_{\infty }$ controller. After the linearization and decoupling steps, 2-DOF continuous and discrete-time controllers were designed using ${\mathcal{H}}_{\infty }$ in [20]. A hybrid architecture using both ${\mathcal{H}}_{\infty }$ and Iterative Learning Control is described in [21]. A linear quadratic regulator (LQR) for MIMO TRAS problem was designed using particle swarm optimization in [22], while a frequency-based PID controller was combined with a lead compensator designed using root locus in [23]. An approach that further details the controller implementation with quantization aspects taken into consideration for the same family of processes is presented in [24].

In this paper, we present a design procedure that manages to optimize the controller parameters instead of tuning them. As such, we present a method for finding the parameters of a MIMO fractional-order PID (FO-PID) robust controller by solving a fixed-structure mixed-sensitivity loop shaping $\mu$-synthesis control problem. Although the resulting control problem is nonconvex in terms of the controller’s free parameters, the nonsmooth optimization techniques implemented in MATLAB’s Robust Control Toolbox can be used. However, the realp object used for these free parameters does not support exponential operations necessary in the approximation of a fractional-order element. Therefore, we present in this paper an algorithm to construct the approximation function of a fractional-order element using integer-order elements and supported arithmetic operations applied on a free parameter isomorphic to the desired fractional order. As such, we successfully manage to formulate the problem of optimizing the parameters of a MIMO FO-PID such that the available techniques can be used. Moreover, we illustrate our design method on the twin rotor aerodynamic system stand, having both MATLAB simulations and physical experiments.

The remainder of this paper is organized as follows: Section 2 summarizes the main mathematical background in terms of available results in both Robust and Fractional-Order Control, along with the description of the proposed method in terms of the algorithm for approximation of the fractional order element and of the optimization problem; Section 3 starts with the presentation of the simplified nonlinear mathematical model of the TRAS system, the linearized mathematical model around an equilibrium point, and a list of parameters with their numerical values and tolerances which manages to encompass the nonlinearities; in Section 4, the numerical results are presented, starting from the augmentation step, followed by the proposed structure of the controller and the obtain results in MATLAB and on the experimental stand; Section 5 presents the discussions of the obtained results and a comparison with another method for solving the optimization problem, while in Section 6 there are some conclusions and possible research directions.

2. Proposed Method

In this section, the mathematical background for the proposed controller design method in terms of Robust Control Framework in Section 2.1 and Fractional-Order Control Framework in Section 2.2 is firstly presented, while in Section 2.3 the method for optimizing the controller parameters using a different approach as against the procedure presented in [17] is described.

2.1. Robust Control

The generalized Robust Control Framework [25] has, besides the control input vector $\mathbf{u}\in {\mathbb{R}}^{n_u}$, two extra inputs: the exogenous input vector $\mathbf{w}\in {\mathbb{R}}^{n_w}$ and disturbance input vector $\mathbf{d}\in {\mathbb{R}}^{n_d}$. Additionally, besides the output vector $\mathbf{y}\in {\mathbb{R}}^{n_y}$, the generalized plant contains two extra outputs: the performance vector $\mathbf{z}\in {\mathbb{R}}^{n_z}$ and the disturbance output $\mathbf{v}\in {\mathbb{R}}^{n_v}$. The input and output disturbance vectors encompass both parametric and unstructured uncertainties, which are generally modeled by the following set:

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

(1)

where $I_n$ denotes the identity matrix of order n.

The uncertainty block $\mathbf{\Delta }$ is interconnected with the generalized plant $P_{\Delta }$ via an upper linear fractional transformation (ULFT), while the controller K is interconnected via a lower linear fractional transformation (LLFT) with $P_{\Delta }$, as noticed in Figure 1.

The state-space representation of the generalized plant $P_{\Delta }$ is:

$$P_{\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).$$

(2)

For robustness analysis, the singular value notion used for ${\mathcal{H}}_{\infty }$ synthesis was extended to the structural singular value, defined for the LLFT interconnection between the plant $P_{\Delta }$ and the controller K according to the uncertainty block $\Delta$ as:

$${\displaystyle {\mu }_{\Delta }\left(\mathtt{LLFT}(P_{\Delta },K)\right)=\underset{\omega \in {\mathbb{R}}_{+}}{sup}\frac{1}{{\displaystyle \underset{\Delta \in \Delta }{min}\left\{\overline{\sigma }(\Delta )\right|det(I-\mathtt{LLFT}(P_{\Delta },K)\left(j\omega \right)\Delta )=0\}}}.}$$

(3)

Given that the problem of explicitly computing such structural singular values is NP-hard, an approximation must be used. The classical ${\mathcal{H}}_{\infty }$ control problem can be extended to the following optimization problem:

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

(4)

which can be considered solved if there is a controller K such that ${\mu }_{\Delta }\left(\mathtt{LLFT}(P_{\Delta },K)\right)<1$, according to main loop theorem. As such, an upper bound is necessary for ${\mu }_{\mathbf{\Delta }}(·)$ [6]:

$${\mu }_{\mathbf{\Delta }}\left(\mathtt{LLFT}(P_{\Delta },K)\left(j\omega \right)\right)\le \underset{D\in \mathit{D}}{inf}\overline{\sigma }\left(D·\mathtt{LLFT}(P_{\Delta },K)\left(j\omega \right)·D^{-1}\right),$$

(5)

where the set $\mathit{D}$ is defined according to the uncertainty block $\mathbf{\Delta }$ as follows [6]:

$$\mathit{D}=\left\{\mathrm{diag}\left(D_1,\dots ,D_s,d_1{I}_{m_1},\dots ,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\}.$$

(6)

Summarizing, robust stability and robust performance are achieved through a controller K obtained as a solution of the optimization problem (4) which manages to obtain an objective value lower than 1. But this NP-hard problem can be approximated by 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_{\Delta },K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right).$$

(7)

As already known, the last optimization problem can be solved using the so-called D–K iteration [9,17]. This iterative procedure starts with a fixed D (usually considered the unitary system) and alternatively computes the controller K, by solving the ${\mathcal{H}}_{\infty }$ control problem with fixed D, and the D-scale factor, by solving the Parrot problem, as defined in [6], for each point from a frequency set $\Omega =\{{\omega }_l={\omega }_1<\cdots <{\omega }_N={\omega }_u\}$ followed an approximation of the obtained solutions with a minimum phase system. Therefore, after setting the initial D-scale step as $D=I$, the following steps are successively applied:

1:

The D-scale step is fixed and the controller can be computed as:

$$K=arg\underset{K\mathrm{stab}.}{inf}{\parallel \mathtt{LLFT}(P_{\Delta },K)\parallel }_{\infty }.$$

(8)

2:

The controller K is fixed and the following set of convex problems must be solved:

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

(9)

for a given frequency range $\Omega$ and, then, a stable minimum phase transfer matrix $D\left(s\right)$ is fitted.

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, the maximum number of iterations is reached, or the improvement after a prescribed number of steps is under an imposed tolerance.

2.2. Fractional-Order Control

The domain of Fractional-Order Control has recently gained more attention due to their robustness. The fractional integral operator used in Control Engineering is [26]:

$${\mathcal{I}}_{\alpha }\left\{f\left(t\right)\right\}=\frac{1}{\Gamma \left(\alpha \right)}{\int }_0^t{(t-\tau )}^{\alpha -1}f\left(\tau \right)\mathrm{d}\tau ,t>0,\alpha \in {\mathbb{R}}_{+},$$

(10)

where $\Gamma (·):{\mathbb{C}}_{+}\to \mathbb{C}$ is the Euler Gamma function. In a similar manner with the integer order integral operator, the fractional order integral operator ${\mathcal{I}}_{\alpha }$ has the following Laplace transform:

$$\mathcal{L}\left\{{\mathcal{I}}_{\alpha }\left\{f\left(t\right)\right\}\right\}\left(s\right)=s^{-\alpha }\mathcal{L}\left\{f\left(t\right)\right\}\left(s\right).$$

(11)

As previously stated, the fractional-order calculus can be used to extend the classical 3-DOF proportional-integral-derivative (PID) controller to a fractional-order PID (FO-PID) having two extra DOF $\lambda ,\mu \in {\mathbb{R}}_{+}$ – the order of the integral operator and the order of the derivative operator, respectively. As such, based on the error signal $\epsilon \left(t\right)$, the command signal $c\left(t\right)$ has the following expression:

$$c\left(t\right)=K_P·\epsilon \left(t\right)+K_I·{\mathcal{I}}_{\lambda }\left\{\epsilon \left(t\right)\right\}+K_D·{\mathcal{I}}_{-\mu }\left\{\epsilon \left(t\right)\right\},$$

(12)

where $c\left(t\right)$ would be $u\left(t\right)$ and $\epsilon \left(t\right)$ would be $r\left(t\right)-y\left(t\right)$ according to the generalized framework from Figure 1, while the differences between the transfer functions of these two controllers are:

$$H_{PID}\left(s\right)=K_P+\frac{K_I}{s}+K_{D}s\Rightarrow H_{FO-PID}\left(s\right)=K_P+\frac{K_I}{s^{\lambda }}+K_D{s}^{\mu }.$$

(13)

The main drawback of the FO-PID revolves around the implementation of the fractional order elements. One possible solution is the Oustaloup recursive approximation (ORA) introduced in crone toolbox [16]. The approximation of a fractional-order element $s^{\lambda }$ with an integer-order one is detailed for $\lambda \in (0,1)$, but it can be easily extended for $\lambda \in \mathbb{R}$. The ORA representation receives as inputs three parameters: the order N of the LTI system which approximates the fractional-order element, along with the lower bound ${\omega }_l$ and the upper bound ${\omega }_u$ of the frequency range where the approximation is valid. The LTI approximation is:

$$s^{\lambda }\approx \prod _{k=1}^N\frac{1+s/{\mathring{\omega }}_k}{1+s/{\widehat{\omega }}_k},$$

(14)

where the poles and zeros frequencies can be computed using 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}},$$

(15)

followed by the recursive relations:

$$\begin{array}{cc}\hfill {\mathring{\omega }}_1& ={\omega }_l\sqrt{\eta },\hfill \end{array}$$

(16a)

$$\begin{array}{cc}\hfill {\widehat{\omega }}_k& ={\mathring{\omega }}_k·\epsilon ,k=\overline{1,N},\hfill \end{array}$$

(16b)

$$\begin{array}{cc}\hfill {\mathring{\omega }}_{k+1}& ={\widehat{\omega }}_k·\eta ,k=\overline{1,N-1}.\hfill \end{array}$$

(16c)

The MATLAB object realp used for fixed-structure robust synthesis does not allow the use of operations other than classical arithmetic operations. Therefore, the recursive fractional-order approximation (14) cannot be used as is in order to compute the fractional-order of the integrative and derivative effects. In Section 2.3 we will give a possible implementation in order to use the realp object for optimizing the controller parameters.

2.3. Controller Design Procedure

Although the controller which results by solving the quasi-convex problem (7) manages to fulfill the robust stability and robust performance, the major drawback consists in the fact that the controller is of high-order and cannot be easily implemented. As such, the problem should be constrained to use a specific controller structure. After imposing a fixed-structure family $\mathcal{K}$, the problem (7) can be written as:

$$\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_{\Delta },K)\left(j\omega \right)·{\left(D\left(j\omega \right)\right)}^{-1}\right).$$

(17)

The above problem is non-convex in terms of the free tuning parameters of the controller $K\in \mathcal{K}$. However, the problem (17) can also be solved using the D–K iteration approach, where the K step from (8) is replaced with the following $K_{\mathcal{K}}$ step:

$$K=arg\underset{\begin{array}{c}K\in \mathcal{K}\\ K\mathrm{stab}\end{array}}{inf}{\parallel \mathtt{LLFT}(P_{\Delta },K)\parallel }_{\infty }.$$

(18)

In the MATLAB environment there exists the realp object which can be used to construct a desired family of controllers $\mathcal{K}$ and then the closed loop system contains both uncertainties and free tunable parameters alike. Using nonsmooth optimization techniques presented in [8] and implemented in [25], the fixed-structure $\mu$-synthesis control problem can be solved. For the purpose of this paper, we consider the fixed structure controller family:

$$\mathcal{K}=\left\{K_{\theta }\left(s\right)=\left(\begin{array}{cccc}{K}_{1,1}\left(s\right)& K_{1,2}\left(s\right)& \dots & K_{1,n_y}\left(s\right)\\ K_{2,1}\left(s\right)& K_{2,2}\left(s\right)& \dots & K_{2,n_y}\left(s\right)\\ ⋮& ⋮& \ddots & ⋮\\ K_{n_u,1}\left(s\right)& K_{n_u,2}\left(s\right)& \dots & K_{n_u,n_y}\left(s\right)\end{array}\right)\left|\theta \in D\right.\right\},$$

(19)

where each controller $K_{i,j}$ has the form:

$$K_{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)}},$$

(20)

having the free parameters:

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

(21)

However, the tunable parameters ${\lambda }^{(i,j)}$ and ${\mu }^{(i,j)}$ cannot be used as realp objects, due to exponential operations not supported. As a solution, ORA is used with the tunable parameter being ${\theta }_{\lambda }\equiv \sqrt{\eta }$ from (15). The transfer function (14) can be implemented using ${\theta }_{\lambda }$ as in Algorithm 1.

Algorithm 1: Construct Fractional-Order Element

Therefore, the tunable parameters for each controller $K_{i,j}\left(s\right)$ are:

$${\widehat{\theta }}_{i,j}=\left(\begin{array}{ccccc}{K}_P^{(i,j)}& K_I^{(i,j)}& K_D^{(i,j)}& {\theta }_{{\lambda }^{(i,j)}}& {\theta }_{{\mu }^{(i,j)}}\end{array}\right)\in {\mathbb{R}}^5,$$

(22)

with the special mention that the parameters ${\theta }_{{\lambda }^{(i,j)}}$ and ${\theta }_{{\mu }^{(i,j)}}$ must be in the domain $\left[1,{\left(\frac{{\omega }_u}{{\omega }_l}\right)}^{\frac{1}{N}}\right]$. If a desired fractional order $\lambda$ is out of the admissible domain, extra integrator/derivative terms can be added. Therefore, the fixed-structure $\mu$-synthesis control problem can be solved in MATLAB from the desired family $\mathcal{K}$ from (19).

Additionally, the control problem will be posed in a mixed-sensitivity loop shaping $\mu$-synthesis formulation. The main reason for this choice consists in the fact that the mixed-sensitivity loop shaping allows an adequate trade-off between robustness and performance. In the optimization process, the following functions will be used for the loop shaping procedure: the sensitivity function S, the complementary sensitivity function T, and the control effort $KS$. For each performance function, a set of performance outputs are considered, while the performance inputs are considered as the references.

On one hand, large magnitude in the open loop system implies good reference tracking, disturbance rejection, and unstable plant stabilization. On the other hand, small magnitude of the open loop system ensures robust stability and mitigation of measurement noise. Moreover, a small magnitude of the control effort is necessary to relieve actuator stress. Although all these magnitude requirements seem to lead to an impossible combination, the target frequency ranges for each component are disjunctive. Through the loop shaping mechanism, the engineer is supposed to find three weighting functions, one for each of the previously-mentioned closed loop performances and the frequency performance imposed by the weighting functions is strongly correlated to the corresponding time performance.

For the sensitivity function, the frequency performance indicators of the weighting function are the minimum bandwidth frequency ${\omega }_B$, which is inversely proportional with the rise time, the maximum magnitude $A_S$ at low frequencies, which imposes the maximum steady-state error, the peak magnitude $M_S$, which limits the overshoot of the system, along with the imposed slope $n_S$ of the sensitivity function at low and medium frequencies [9]:

$$W_S\left(s\right)={\left(\frac{\frac{1}{M_S^{1/n_S}}s+{\omega }_B}{s+{\omega }_B{A}_S^{1/n_S}}\right)}^{n_S}.$$

(23)

Similarly, the complementary sensitivity’s weighting function can be constructed using the peak amplitude $M_T$, the maximum magnitude at high frequencies $A_T$, the minimum bandwidth ${\omega }_{BT}$ and the roll-off $n_T$:

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

(24)

The control effort is generally weighted by imposing the magnitude at low and high frequencies, along with an intermediate point of interest. However, the main goal is to maintain the control effort in the range given by the saturation of the physical actuator. For MIMO systems, the weighting matrices are diagonal concatenations of the weighting functions described above. Now the optimization problem that needs to be solved for the proposed method is the mixed-sensitivity fixed-structure loop shaping $\mu$-synthesis:

$$\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}.\parallel \left(\begin{array}{ccc}{W}_{S}S& W_{T}T& W_{KS}KS\end{array}\right){\parallel }_{\infty }<1\hfill \end{array}.$$

(25)

3. Mathematical Model of a TRAS

The TRAS model is of sixth order with four inputs and two outputs. The state variables considered are the rotational speed of the tail rotor (${\omega }_h$), the rotational speed of the main rotor (${\omega }_v$), the azimuth velocity of TRAS beam (${\Omega }_h$), the pitch velocity of TRAS beam (${\Omega }_v$), the azimuth position (${\alpha }_h$), and the pitch position (${\alpha }_v$), the state vector being:

$$\mathbf{x}={\left(\begin{array}{cccccc}{\omega }_h& {\omega }_v& {\Omega }_h& {\Omega }_v& {\alpha }_h& {\alpha }_v\end{array}\right)}^{\top }\in {\mathbb{R}}^6.$$

(26)

There are two control inputs, $u_h$ and $u_v$, representing the normalized horizontal and vertical DC-motor PWM duty cycles, while the considered outputs will be the azimuth and pitch positions of the TRAS beam:

$$\mathbf{u}=\left(\begin{array}{cc}{u}_h& u_v\end{array}\right)\in {\mathbb{R}}^2,\mathbf{y}=\left(\begin{array}{cc}{\alpha }_h& {\alpha }_v\end{array}\right)\in {\mathbb{R}}^2.$$

(27)

The TRAS model is strongly nonlinear even under some simplifying assumptions, as stated in [27]. One simplification is regarding to the characteristics of the two rotors: their models are supposed to be of first order containing the moment of inertia and the velocity gain for each rotor. Moreover, the angular velocities of the TRAS beam is influenced by the aerodynamic force of each rotor, which is nonlinear in terms of its rotational speed, by the aerodynamic damping torque and by the cross momentum. Moreover, the azimuth velocity is strongly influenced by the pitch angle position, while the pitch velocity is influenced by the pitch angle as well by the return torque. The nonlinear model after some simplifying assumptions can be written as:

$${\dot{\omega }}_h=-\frac{1}{I_h}{f}_1\left({\omega }_h\right)+\frac{1}{I_h}{u}_h$$

(28)

$${\dot{\Omega }}_h=\frac{l_t·}{k_1·{cos}^2\left({\alpha }_v\right)+k_2}{f}_2\left({\omega }_h\right)·cos\left({\alpha }_v\right)-\frac{k_{fh}}{k_1·{cos}^2\left({\alpha }_v\right)+k_2}{\Omega }_h-\frac{k_{vh}}{k_1·{cos}^2\left({\alpha }_v\right)+k_2}cos\left({\alpha }_v\right)·u_v$$

(29)

$${\dot{\alpha }}_h={\Omega }_h$$

(30)

$${\dot{\omega }}_v=-\frac{1}{I_v}{f}_3\left({\omega }_v\right)+\frac{1}{I_v}{u}_v$$

(31)

$${\dot{\Omega }}_v=\frac{l_m}{J_v}{f}_4\left({\omega }_v\right)-\frac{k_{fv}}{J_v}{\Omega }_v-\frac{k_{3}cos\left({\alpha }_v\right)+k_{4}sin\left({\alpha }_v\right)+k_{5}sin\left({\alpha }_v\right)cos\left({\alpha }_v\right)}{J_v}+\frac{k_{hv}}{J_v}{u}_h$$

(32)

$${\dot{\alpha }}_v={\Omega }_v$$

(33)

All parameters of both linearized an nonlinear systems are described in

Table 1. Twin rotor aerodynamic system physical parameters, values and tolerances.

Symbols	Description	Nominal Value	Tolerance
$I_h$	moment of inertia of the tail rotor	1/37,000 (kg·m${}^2$)	-
$I_v$	moment of inertia of the main rotor	1/6100 (kg·m${}^2$)	-
$J_h$	moment of inertia with respect to the vertical axis	0.0268 (kg·m${}^2$)	$\pm 10[\%]$
$J_v$	moment of inertia with respect to the horizontal axis	0.0268 (kg·m${}^2$)	-
$k_{H_h}$	velocity gain of the tail rotor	$7.0742\times {10}^3$ (rad/s)	$\pm 10[\%]$
$k_{H_v}$	velocity gain of the main rotor	$5.1574\times {10}^3$ (rad/s)	$\pm 10[\%]$
$k_{F_h}$	thrust coefficient of the tail rotor	$1.3218\times {10}^{-4}$ (Ns/rad)	$\pm 10[\%]$
$k_{F_v}$	thrust coefficient of the main rotor	$2.0124\times {10}^{-4}$ (Ns/rad)	$\pm 10[\%]$
$k_{f_h}$	friction coefficient in the vertical axis	$5.889\times {10}^{-3}$ (Nms/rad)	$\pm 5[\%]$
$k_{f_v}$	friction coefficient in the horizontal axis	$1.271\times {10}^{-2}$ (Nms/rad)	$\pm 5[\%]$
$k_{hv}$	coefficient of the cross moment from tail rotor to pitch angle	$4.175\times {10}^{-3}$ (Nm)	$\pm 5[\%]$
$k_{vh}$	coefficient of the cross moment from main rotor to azimuth angle	$-1.782\times {10}^{-2}$	$\pm 5[\%]$
$R_v$	coefficient of the return torque	$9.360078\times {10}^{-2}$ (Nm)	$\pm 10[\%]$
$l_t$	length of the tail part of the beam	$0.2165$ (m)	-
$l_m$	length of the main part of the beam	$0.202$ (m)	-
$k_1$	coefficient of $J_h$	$2.379\times {10}^{-2}$ (kg·m${}^2$)	-
$k_2$	coefficient of $J_h$	$3.009\times {10}^{-3}$ (kg·m${}^2$)	-
$k_3$	coefficient of $R_v$	$5.006\times {10}^{-2}$ (Nm)	-
$k_4$	coefficient of $R_v$	$9.361\times {10}^{-2}$ (Nm)	-
$k_5$	auxiliary coefficient	$0.010624$ (Nm)	-

. The first step of the linearization process is to find approximations for the functions $f_1$ and $f_3$ such that the two systems from inputs to rotational speeds of the rotors are of first order. In order to obtain this scenario, these functions are estimated as $f_1\left({\omega }_h\right)=k_{Hh}·{\omega }_h$ and $f_3\left({\omega }_v\right)=k_{Hv}·{\omega }_v$, while the nonlinerity is treated using the sector bound technique, being included in the tolerance of each velocity gain. Moreover, the forces developed by each axis are also nonlinear in terms of rotational speeds of the rotors and can be approximated $f_2\left({\omega }_h\right)=k_{Fh}·{\omega }_h$ and $f_4\left({\omega }_v\right)=k_{Fv}·{\omega }_v$, where the trust coefficients encompass the nonlinearities in their tolerances. All sector bound nonlinearites described above are depicted in Figure 2.

After this first step, the nonlinear model can be now linearized around an equilibrium point. The forced equilibrium point has been chosen such that the outputs are ${\overline{\alpha }}_h={\overline{\alpha }}_v$ = 0 [rad], i.e., plant stabilization problem. In order to obtain this point, the state vector has the rest of the components ${\overline{\omega }}_h=-1336$ [rad/s], ${\overline{\omega }}_v=1803.45$ [rad/s], ${\overline{\Omega }}_h=0$ [rad/s], ${\overline{\Omega }}_v=0$ [rad/s], while the input vector has the components ${\overline{u}}_h=-0.1492$ and ${\overline{u}}_v=0.30559$. According to [27], the moment of inertia with respect to the horizontal axis is constant, while around the vertical axis the moment of inertia is nonlinear, having the expression $J_h=k_1·{cos}^2\left({\alpha }_v\right)+k_2$. In practice, we will consider this parameter uncertain, having the nominal value ${\overline{J}}_h=k_1·cos\left({\overline{\alpha }}_v\right)+k_2$, along with a tolerance of $\pm 10[\%]$. The uncertainties from the thrust coefficients of the tail and the main rotors are necessary in order to compensate the nonlinearity of the aerodynamic forces from these rotors. The friction coefficients in the axes and the cross moments coefficients also present uncertainties in order to compensate the nonlinearities presented in the angular velocity parts and the interconnections between the two rotations. The return torque coefficient is a nonlinear function in terms of pitch position and velocity, which can be approximated by an uncertain parameter having the nominal value ${\overline{R}}_v=k_{3}sin\left({\overline{\alpha }}_v\right)-k_{4}cos\left({\overline{\alpha }}_v\right)$, and a tolerance of $\pm 10\%$. As such, the linearized state-space model can now be written as:

$$\dot{\mathbf{x}}\left(t\right)=\left(\begin{array}{cccccc}-\frac{1}{I_h·k_{H_h}}& 0& 0& 0& 0& 0\\ \frac{l_t·k_{F_h}·cos\left({\overline{\alpha }}_v\right)}{J_h}& -\frac{k_{f_h}}{J_h}& 0& 0& 0& -\frac{l_t·k_{F_h}·sin\left({\overline{\alpha }}_v\right)+k_{vh}·sin\left({\overline{\alpha }}_v\right)·{\overline{u}}_v}{J_h}\\ 0& 1& 0& 0& 0& 0\\ 0& 0& 0& -\frac{1}{k_{H_v}·I_v}& 0& 0\\ 0& 0& 0& \frac{l_m·k_{F_v}}{J_v}& -\frac{k_{f_v}}{J_v}& -\frac{R_v+2k_{5}cos\left(2{\overline{\alpha }}_v\right)}{J_v}\\ 0& 0& 0& 0& 1& 0\end{array}\right)\mathbf{x}\left(t\right)+\left(\begin{array}{cc}\frac{1}{I_h}& 0\\ 0& -\frac{k_{vh}·cos\left({\overline{\alpha }}_v\right)}{J_h}\\ 0& 0\\ 0& \frac{1}{I_v}\\ \frac{k_{hv}}{J_v}& 0\\ 0& 0\end{array}\right)\mathbf{u}\left(t\right);$$

(34a)

$$\mathbf{y}\left(t\right)=\left(\begin{array}{cccccc}0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)\mathbf{x}\left(t\right)+\left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right)\mathbf{u}\left(t\right).$$

(34b)

The singular values of the twin rotor aerodynamic system plant having the parameters presented in Table 1, before augmentation, are presented in Figure 3.

4. Numerical Results

The controller design procedure proposed in this paper will be applied on a twin rotor aerodynamic system (TRAS). The physical stand from INTECO [27] is presented in Figure 4. The numerical values of the parameters described in Section 3 are presented in Table 1, along with their nominal values and tolerances.

In order to illustrate the power of the proposed method, a comparison between the numerical simulations for the linearized system using MATLAB and the experimental results on the physical stand has been performed. For the numerical results, the block diagram is presented in Figure 5, where the reference signals ${\mathbf{w}}_1\equiv \mathbf{r}={\left(\begin{array}{cc}{\alpha }_h^{\star }& {\alpha }_v^{\star }\end{array}\right)}^{\top }$ are considered the inputs of the linearized system, while the performance output vector is:

$$\mathbf{z}=\left(\begin{array}{cccccc}{z}_{S,{\alpha }_h}& z_{S,{\alpha }_v}& z_{T,{\alpha }_h}& z_{T,{\alpha }_v}& z_{KS,{\alpha }_h}& z_{KS,{\alpha }_v}\end{array}\right).$$

(35)

For the numerical simulation part, the plant augmentation has been done with the following weighting functions parameters: ${\omega }_{B,{\alpha }_h}=0.2$ [rad/s], ${\omega }_{B,{\alpha }_v}=0.05$ [rad/s], $A_{S,{\alpha }_v}=A_{S,{\alpha }_h}=1\times {10}^{-2}$, $M_{S,{\alpha }_v}=M_{S,{\alpha }_h}=2$, $n_{S,{\alpha }_v}=n_{S,{\alpha }_h}=1$ (the reference is considered to be a unity step signal), ${\omega }_{BT,{\alpha }_h}=20$ [rad/s], ${\omega }_{B,{\alpha }_v}=5$ [rad/s], $A_{T,{\alpha }_v}=A_{T,{\alpha }_h}=1\times {10}^{-2}$, $M_{T,{\alpha }_v}=M_{T,{\alpha }_h}=2$, $n_{T,{\alpha }_v}=n_{T,{\alpha }_h}=1$, while the $DC$ component of the control effort weighting functions is 1, being the maximum value of the command signal, and the maximum value at high-frequency is of magnitude 5. The weighting functions result as follows:

$$W_S\left(s\right)=\left(\begin{array}{cc}{W}_{S,{\alpha }_h}\left(s\right)& 0\\ 0& W_{S,{\alpha }_v}\left(s\right)\end{array}\right),\mathrm{where}{W}_{S,{\alpha }_h}\left(s\right)=\frac{0.5s+0.2}{s+2\times {10}^{-3}},W_{S,{\alpha }_v}\left(s\right)=\frac{0.5s+0.05}{s+5\times {10}^{-4}},$$

(36)

$$W_T\left(s\right)=\left(\begin{array}{cc}{W}_{T,{\alpha }_h}\left(s\right)& 0\\ 0& W_{T,{\alpha }_v}\left(s\right)\end{array}\right),\mathrm{where}{W}_{T,{\alpha }_h}\left(s\right)=\frac{s+20}{0.01s+40},W_{T,{\alpha }_v}\left(s\right)=\frac{s+5}{0.01s+10},$$

(37)

$$W_{KS}\left(s\right)=\left(\begin{array}{cc}{W}_{KS,{\alpha }_h}\left(s\right)& 0\\ 0& W_{KS,{\alpha }_v}\left(s\right)\end{array}\right),\mathrm{where}{W}_{KS,{\alpha }_h}\left(s\right)=W_{KS,{\alpha }_v}\left(s\right)=\frac{0.2s+0.8532}{s+0.8532}.$$

(38)

As noted in Figure 3, the frequency range is between ${\omega }_l=1\times {10}^{-2}$ [rad/s] and ${\omega }_u=1\times {10}^3$ [rad/s], which will be also used for ORA, along with the order of approximation $N=5$. Using the augmented plant presented in Figure 5, the fixed-structure mixed-sensitivity loop shaping $\mu$-synthesis problem (25) is solved using the musyn command from MATLAB with the following specifications: the maximum number of D–K iterations is 10, the threshold for the upper bound of the ${\mu }_{\mathbf{\Delta }}\left(\mathtt{LLFT}(P_{\mathrm{aug}},K)\right)$ is 1, and the maximum number of iterations for asserting the lack of progress is 4.

The fixed-structured $\mu$-synthesis control problem was solved using three D–K iterations, having the upper bound of the structured singular value ${\mu }_{\mathbf{\Delta }}\left(\mathtt{LLFT}(P_{\mathrm{aug}},K)\right)\le 0.9902<1$, which means that the resulting FO-PID controller manages to fulfill both robust stability and robust performance. The resulting FO-PID controller is:

$$K_{FO-PID}^{{\theta }^{\star }}\left(s\right)=\left(\begin{array}{cc}0.1149+0.0603·s^{-1.267}+0.0909·\frac{s^{1.1442}}{0.1154s+1}& -7.1329+5.0864\times {10}^3\frac{s^{1.0001}}{712.97s+1}\\ 0.0315-0.0832\frac{s^{1.2251}}{27.3377s+1}& -0.0297+0.1013·s^{-1.0001}+0.0232·\frac{s^{1.2851}}{0.149s+1}\end{array}\right),$$

(39)

where the low-pass component needs to be added in order to implement the derivative element of order greater than 1 having one extra degree of freedom for each such element. The results obtained after each step are summarized 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
Peak Value of$\mu$ (FO-PID)	$2.657$	$1.066$	$1.007$	$\mathbf{0.9902}$

, where after x steps the controller design problem has been successfully solved. The upper bound of the structural singular value is presented in Figure 6.

In order to illustrate the frequency-domain performance, the sensitivity function, complementary sensitivity function and control effort are presented in Figure 7. The nominal plant has been analyzed along with 100 Monte Carlo simulations for the given uncertainty range. Also, in order to underline that the control problem has been successfully solved, the weighting functions are also depicted and it can be noticed that all the simulated functions are under the imposed thresholds. Additionally, the Bode magnitude characteristics of the resulting controller are provided in Figure 8.

The time-domain performance of the lower linear fractional transform between the linearized plant and the controller are presented using a step response in Figure 9. In a similar manner, the nominal plant is illustrated along with 100 Monte Carlo simulations. The rise time for the azimuth position varies between $0.796$ [s] and $1.05$ [s], having a settling time between $14.8$ [s] and $16.1$ [s] and an overshoot between $16.9[\%]$ and $24.2[\%]$, with no steady-state error. Similarly, the rise time for the pitch position is between $6.79$ [s] and $11.7$ [s], having a settling time between $10.7$ [s] and 19 [s], with no overshoot and no steady-state error.

Numerical results will further be compared with the experimental results. The first set of experiments, shown in Figure 10, have been made for a square reference with an amplitude of $\pm 0.1$ [rad] and a period of 100 [s] for both axes. The initial conditions were varied in practice in order to illustrate the capability of the method. It can be noticed that for the azimuth position, the practical overshoot is a bit higher than in the linear case, with a comparable settling time and near-zero steady-state error due to the quantization effects. The pitch position presents overshoot for the initial step, while for the second step the behavior is similar to that of the linear system, with no overshoot, no steady-state error and comparable settling time.

The second set of experiments are made for the stabilization problem, where the reference for both axes is ${\alpha }_h^{\star }={\alpha }_v^{\star }=0.1$ [rad]. It can be noticed that the azimuth position presents an initial overshoot comparable to that obtained for the linear system, while the second part of the oscillation is more aggressive, underling the influence of the nonlinear components. In a similar manner, the pitch position presents an overshoot along with several oscillations, while the settling time is comparable with the linear case’s. The experimental results are shown in Figure 11. Moreover, three different disturbances have been applied after 50 [s]: a perturbation on the vertical axis which leads the pitch position at the maximum value (blue), a perturbation on the horizontal axis with the same characteristics (cyan), and another small perturbation on the horizontal axis (black). It can be noticed that all disturbances have been successfully rejected.

Finally, the third set of experiments, depicted in Figure 12, illustrates the behaviour of the proposed method for an operating point far from the forced equilibrium point used in the linearization procedure and controller synthesis. As such, a step of ${\alpha }_h^{\star }={\alpha }_v^{\star }=0.7$ [rad] has been initially applied, with a different pair ${\alpha }_h^{\star }={\alpha }_v^{\star }=-0.3$ [rad] applied at the moment $t_1=50$ [s]. For the horizontal axis, the overshoot is a bit higher than in the linear case, but with similar settling time and no steady-state error. Also, for the vertical axis, the overshoot is negligible for the first step and zero for the second step, having comparable settling times and no steady-state error. Therefore, the controller can be used for other operating points.

5. Discussion

In order to compare the current iteration of FO-PID $\mu$-synthesis with previous methods, the fixed-structure part of the $\mu$-synthesis mixed-sensitivity loop shaping control problem (25) is solved using the artificial bee colony (ABC) approach presented in [17]. The hyperparameters used for this experiment are: the swarm dimension $N=50$, the maximum number of ABC cycles 50, the maximum number of cycles with no improvements 10, the limit for the abandonment counter 10, the maximum number of D–K iterations 10, the maximum window length for assessing lack of progress 4, while the parameters for the cost functions are $\alpha =1$ and $\beta ={10}^5$. Using this setup, the fixed-structure mixed-sensitivity loop shaping $\mu$-synthesis problem (25) is solved using five D–-K iterations. The resulting controller is:

$$K_{FO-PID,ABC}^{{\theta }^{\star }}\left(s\right)=\left(\begin{array}{cc}0.1642+0.0892·s^{-1.1834}+0.0913·\frac{s^{1.1209}}{0.1173s+1}& -0.0016+0.8355\frac{s^{1.0001}}{104.8s+1}\\ 0.0106-6.769\times {10}^{-4}\frac{s^{1.4938}}{100s+1}& -0.0741+0.0981·s^{-1.185}+0.001·\frac{s^{1.157}}{0.0735s+1}\end{array}\right),$$

(40)

Additionally, an experiment with unstructured $\mu$-synthesis has also been performed, leading to an upper bound of the structured singular value of $99.86$ and a controller of 34th order, which means that robust stability and robust performance are not guaranteed, with the controller additionally of high order. The optimization algorithm has been stopped after three iterations because the diverging stopping criteria has been reached. A summary of the obtained results with the proposed method, the ABC method [17] and the unstructured $\mu$-synthesis is presented in

Table 3. A comparison between the evolution of the structural singular values in the D–K iteration procedure used to solve the mixed-sensitivity fixed structure $\mu$-synthesis problem for the case study.

D–K Iteration Number	1	2	3	4	5
Peak Value of$\mu$ (FO-PID)	$2.657$	$1.066$	$1.007$	$\mathbf{0.9902}$	-
Peak Value of$\mu$ (unstructured)	100	$99.8$	$99.4748$	-	-
Peak Value of$\mu$ (ABC approach [17])	$105.7741$	$2.4095$	$1.2452$	$1.1255$	$\mathbf{0.9989}$

.

As such, the unstructured version of the $\mu$-synthesis control problem could not be solved, resulting a high-order controller which does not guarantee robust stability and performance. On the other hand, both remaining methods managed to solve the control problem described in (25). The new method introduced in this paper manages to solve the problem faster due to the advantages of the nonsmooth optimization techniques implemented in MATLAB.

As a future iteration, we propose to find a decentralized controller having a nonlinear component and a linear and time-invariant (LTI) robust component. The nonlinear component needs to be designed such that the lower linear fractional transform interconnection between the plant and such a component is asymptotically stable, as in [28], where the passivity-based control framework has been extended for quasi-linear input-affine systems. Additionally, the LTI robust controller can be designed using the proposed method, the decentralized controller managing to ensure robust stability and performance for the nonlinear system. Moreover, the fractional-order element can be approximated using the presented ORA method only for $s^{\lambda }$, with $\lambda \in (0,1)$. As such, another research direction is to find a method to integrate all positive values of $\lambda \in {\mathbb{R}}_{+}$. On the other hand, the presented methods were considered in the continuous-time domain, although, for practical implementation, the controller must be discretized and also quantized. A starting point for this aspect could be the work presented in [24].

6. Conclusions

The current paper presents an algorithm which manages to integrate the MIMO fractional-order PID (FO-PID) controller in the fixed-structure mixed-sensitivity loop shaping $\mu$-synthesis control problem by constructing an element isomorphic with the fractional order. In order to expose the method capacity and potential, a twin rotor aerodynamic system experimental stand has been utilized. After the simplified nonlinear and linearized models were presented, the linear system has been augmented with weighting functions which managed to impose the desired performance. The fixed-structure $\mu$-synthesis control problem has been successfully solved using four D–K iterations, resulting a controller which manages to ensure both robust stability and robust performance. A comparative analysis between the results obtained with the designed controller used for the linearized plant and for the practical experimental stand has also been performed.

As future work, the proposed design method will be added into a next iteration of the toolbox initially proposed in [9] in order to automatically perform the fractional-order fixed-structure $\mu$-synthesis. Also, the proposed method can be integrated into a control scheme with a decentralized controller having an extra nonlinear component which ensures that the robust stability and robust performance are also guaranteed for the nonlinear system.