2.1 Analytical computation of inertia matrix and CoG coordinates

To simplify the analytical modeling of the inertia matrix and CoG in terms of joint angles, the following assumptions are considered:

• • The mechanical components of the RSQ are assumed to be rigid.

• • The four joints are modeled as ideal pivot joints.

• • All designed mechanical components are assumed to have a uniformly distributed mass.

• • Some components such as cables, screws and nuts are considered as additional mass, uniformly distributed across the four arms.

2.1.1 Calculation of the CoG

The CoG of the RSQ is calculated from the individual $ Co{G}_{s}$ of its main components: the end arms, the middle arms and the central base. Since these components are assumed to be rigid bodies, their individual CoG $_{s}$ , denoted as $ {G}_{ij}$ for the arms and $ {G}_{\beta }$ for the central base, remain constant with respect to the local frame. Their values are obtained from SolidWorks based on the CAD model in Fig. 2, where the precise geometry and weight of each component is considered. The position of each arm’s frame relative to the base frame is computed as:

(1) \begin{align} \left\{\begin{array}{l}\overrightarrow{{O}_{\beta }{O}_{\beta }}=\left[\begin{array}{l}0\\ 0\\ 0\end{array}\right]\!,\,\,\, \overrightarrow{{O}_{\beta }{O}_{11}}=\left[\begin{array}{l}-\frac{{L}_{s}}{2}-\frac{L}{2}\mathrm{c}\mathrm{o}\mathrm{s}\left({\xi }_{11}\right)\\[5pt] -\frac{L}{2}\mathrm{s}\mathrm{i}\mathrm{n}\left({\xi }_{11}\right)\\ 0\end{array}\right]\!,\,\,\, \overrightarrow{{O}_{\beta }{O}_{12}}=\left[\begin{array}{l}-\frac{{L}_{s}}{2}-L\mathrm{c}\mathrm{o}\mathrm{s}\left({\xi }_{11}\right)-\frac{L}{2}\mathrm{c}\mathrm{o}\mathrm{s}\big(\sum _{j=1}^{2} {\xi }_{1j}\big)\\[6pt] -L\mathrm{s}\mathrm{i}\mathrm{n}\left({\xi }_{11}\right)-\frac{L}{2}\mathrm{s}\mathrm{i}\mathrm{n}\big(\sum _{j=1}^{2} {\xi }_{1j}\big)\\ 0\end{array}\right]\\[20pt] \overrightarrow{{O}_{\beta }{O}_{21}}=\left[\begin{array}{l}\frac{{L}_{s}}{2}+\frac{L}{2}\mathrm{c}\mathrm{o}\mathrm{s}\left({\xi }_{21}\right)\\[5pt] -\frac{L}{2}\mathrm{s}\mathrm{i}\mathrm{n}\left({\xi }_{21}\right)\\ 0\end{array}\right]\!,\quad \overrightarrow{{O}_{\beta }{O}_{22}}=\left[\begin{array}{l}+\frac{{L}_{s}}{2}+L\mathrm{c}\mathrm{o}\mathrm{s}\left({\xi }_{21}\right)+\frac{L}{2}\mathrm{c}\mathrm{o}\mathrm{s}\big(\sum _{j=1}^{2} {\xi }_{2j}\big)\\[6pt] -L\mathrm{s}\mathrm{i}\mathrm{n}\left({\xi }_{21}\right)-\frac{L}{2}\mathrm{s}\mathrm{i}\mathrm{n}\big(\sum _{j=1}^{2} {\xi }_{2j}\big)\\ 0\end{array}\right]\end{array}\right.\end{align}

where $\overrightarrow{{O}_{\beta }{O}_{ij}}\in {\mathbb{R}}^{3\times 1}$ denotes the position vector of the arm frame $ {O}_{ij}$ relative to the base frame $ {O}_{\beta }$ , while $ L$ and $ {L}_{s}$ represent the arm length and the central base length, respectively.

The CoG of the RSQ is computed as the centroid of the CoG of the four arms and the central base, expressed in the base frame. It is given by Equation (2).

(2) \begin{align} \begin{array}{ll}& \overrightarrow{{O}_{\beta }{O}_{G}}=\dfrac{\sum _{i,j=\mathrm{1,1}}^{i,j=\mathrm{2,2}} {M}_{ij}\overrightarrow{{O}_{\beta }{G}_{ij}}+{M}_{\beta }\overrightarrow{{O}_{\beta }{G}_{\beta }}}{2{M}_{11}+2{M}_{12}+{M}_{\beta }}\end{array}\end{align}

where $ {M}_{ij}$ and $ {M}_{\beta }$ denote the masses of the arms $ \left(ij\right)$ and the central base, respectively. The vector $\overrightarrow{{O}_{\beta }{G}_{ij}}\in {\mathbb{R}}^{3\times 1}$ represents the position of the centre of mass $ {G}_{ij}$ of each arm relative to the base frame $ {O}_{\beta }$ , while $\overrightarrow{{O}_{\beta }{O}_{G}}\in {\mathbb{R}}^{3\times 1}$ denotes the offset between the base frame and the body frame, corresponding to the global CoG of the RSQ.

The positions of $ {G}_{ij}$ are expressed in the base frame $ {O}_{\beta }$ using Huygens’ theorem, as mentioned in Equation (3).

(3) \begin{align} \begin{array}{ll}& \overrightarrow{{O}_{\beta }{G}_{ij}}=\overrightarrow{{O}_{\beta }{O}_{ij}}+{R}_{{O}_{ij}\to {O}_{\beta }}\overrightarrow{{O}_{ij}{G}_{ij}}\end{array}\end{align}

with

(4) \begin{align}\left\{ \begin{matrix} {{R_{{O_{22}}\mathop \to {O_\beta }}} = R\left( {{\xi _{21}} + {\xi _{22}}} \right){\rm{\;\;\;\;}}{R_{{O_{21}}\mathop \to {O_\beta }}} = R\left( {{\xi _{21}}} \right)} \\[4pt] {{R_{{O_{12}}\mathop \to {O_\beta }}} = {R^T}\left( {{\xi _{11}} + {\xi _{12}}} \right){\rm{\;\;\;\;}}{R_{{O_{11}}\mathop \to {O_\beta }}} = {R^T}\left( {{\xi _{11}}} \right)} \end{matrix}\right.\end{align}

and

(5) \begin{align} R\left({\xi }_{ij}\right)=\left[\begin{array}{l@{\quad}l@{\quad}l}Cos\left({\xi }_{ij}\right)& Sin\left({\xi }_{ij}\right)& 0\\[4pt] -Sin\left({\xi }_{ij}\right)& Cos\left({\xi }_{ij}\right)& 0\\[4pt] 0& 0& 1\end{array}\right]\end{align}

where $ R\in \mathrm{S}\mathrm{O}\left(3\right)$ denotes the rotation matrix, and $ {R}^{T}$ represents its transpose.

2.1.2 Calculation of the inertia matrix

To derive the inertial matrix of the RSQ, the moment of inertia of each of its main components, i.e. main base and arms, is estimated using SolidWorks software. These parameters are then combined using the parallel axis theorem to account for the overall structure.

(6) \begin{align} I\left({\xi }_{ij}\right(t))=\sum _{i,j=11}^{i,j=22} \big({R}_{{O}_{ij}\to {O}_{\beta }}{I}_{Ar{m}_{ij}/{O}_{ij}}{R}_{{O}_{i\mathrm{j}}\to {O}_{\beta }}^{-1}+{M}_{Ar{m}_{ij}}{I}_{{O}_{ij}}\big)+{I}_{base/{O}_{\beta }}-M{J}_{{O}_{G}}\end{align}

where $ I\left({\xi }_{ij}\right(t))\in {\mathbb{R}}^{3\times 3}$ denotes the inertia tensor of the RSQ expressed in the body frame $ {O}_{G}$ , which depends on the configuration variables $ {\xi }_{ij}\left(t\right)$ , $ M$ denotes the total mass of the system, $ {I}_{Ar{m}_{ij}/{O}_{ij}}$ denotes the inertia tensor of each arm expressed in its local frame $ {O}_{ij}$ and $ {I}_{base/{O}_{\beta }}$ denotes the inertia tensor of the central base expressed in the base frame $ {O}_{\beta }$ . Furthermore, $ {J}_{({O}_{ij},G)}$ denotes the Huygens (parallel-axis) matrix associated with the displacement between $ {O}_{ij}$ and $ G$ , with

(7) \begin{align} {J}_{({O}_{ij},G)}=\left[\begin{array}{c@{\quad}c@{\quad}c}{y}_{({O}_{ij},G)}^{2}& -{x}_{({O}_{ij},G)}{y}_{({O}_{ij},G)}& 0\\ -{y}_{({O}_{ij},G)}{x}_{({O}_{ij},G)}& {{x}_{({O}_{ij},G)}}^{2}& 0\\ 0& 0& {x}_{({O}_{ij},G)}^{2}+{y}_{({O}_{ij},G)}^{2}\end{array}\right]\end{align}

where ( $ {x}_{{O}_{ij}}$ , $ {y}_{{O}_{ij}}$ ), and ( $ {x}_{G}$ , $ {y}_{G}$ ) represent the coordinates of the arm frame and the CoG frame, respectively, with respect to the base frame.

2.2 Validation of the analytical formulation via SolidWorks CAD modeling

To validate the equations established for the online computation of inertial parameters, four configurations were considered: ‘O’, ‘U’, ‘V’ and ‘L’, as depicted in Tables 2 and 3. The analytical results were compared with those obtained from SolidWorks, which provides these parameters in a discrete manner.

Table 2.

CoG coordinates obtained from the analytical formulation (blue) and SolidWorks CAD modeling (red)

Table 2 Long description

A table comparing coordinates obtained from analytical formulation and SolidWorks CAD modeling for four configurations. The table has four rows and eleven columns. The columns are labeled with different parameters and configurations. The rows are labeled with specific values for each parameter. The configurations are O, V, L, and U. The parameters include ξ22, ξ11, ξ12 in degrees, xG, yG, and zG in meters. The values for each parameter and configuration are listed in the table. For example, for configuration O, the values are: ξ22, ξ11, ξ12: 80, 55, 55, 80; xG: -0.00002801, -0.00002914, -0.32798483, -0.32798424, -0.00182995, -0.00182896; yG: -0.00002801, -0.00002914, -0.28407248, -0.28407195, -0.00182995, -0.00182896; zG: -0.06420157, -0.06420168, -0.06538202, -0.06538302, -0.00182995, -0.00182896. Similar values are provided for configurations V, L, and U.

Table 3.

Inertial parameters computed in blue and those obtained from SolidWorks in red

Table 3 Long description

A table comparing inertial parameters computed in blue and those obtained from SolidWorks in red. The table has 4 rows and 5 columns. The columns are labeled I subscript xx, I subscript yy, I subscript zz, I subscript xy, and I subscript xz. The rows are labeled O, V, L, and U. The values in the table are as follows: Row O: I subscript xx, 0.07790815; I subscript yy, 0.07790808; I subscript zz, 0.05895824; I subscript xy, 0.00000207; I subscript xz, -0.00000731. Row V: I subscript xx, 0.08396465; I subscript yy, 0.08396422; I subscript zz, 0.27123847; I subscript xy, 0.00000135; I subscript xz, -0.0000792. Row L: I subscript xx, 0.02905546; I subscript yy, 0.02905514; I subscript zz, 0.49772362; I subscript xy, 0.00001434; I subscript xz, -0.00129063. Row U: I subscript xx, 0.04407592; I subscript yy, 0.04407509; I subscript zz, 0.34358318; I subscript xy, 0.00000952; I subscript xz, -0.00300911.

Tables 2 and 3 summarise the obtained results for the four configurations. The findings show strong consistency with those obtained from SolidWorks, with an average positional error for the CoG below $ 1{0}^{-6}$ m and errors in moments and products of inertia less than $ 1{0}^{-7}$ kg $ \cdot $ m $ ^{2}$ .

These results also demonstrate a considerable variation of the inertial parameters among the four studied configurations. The elements of the inertia matrix vary by up to 844% for $ {I}_{yy}$ and up to 400% for $ {I}_{\mathrm{x}x}$ and $ {I}_{zz}$ . Additionally, the CoG on the other hand moves by 0.2617 m along the Y-axis and 0.064 m along the X-axis when transitioning from the ‘O’ to the ‘L’ configuration. These substantial variations must be considered in the control design to prevent potential stability issues.

To analyze the variation of the inertia matrix parameters and the CoG across all intermediate configurations. This study excludes singular configurations presented in Fig. 6, as well as configurations resulting in rotor overlap.

Figure 4 plots the position of the CoG for the studied configurations. As the figure shows, the CoG coordinates vary between −30 cm and 30 cm along the X-axis and between −42.33 cm and 42.33 cm along the Y-axis, highlighting the importance of the online calculation of the CoG position. Indeed, this position is crucial for computing the control allocation matrix and, consequently, for ensuring effective control of the torque inputs. The range of variation of the inertia matrix parameters across the assessed configurations is illustrated in Fig. 5.

Figure 4.

Range of CoG variation with respect to $ {O}_{\beta }$ for all possible configurations of the RSQ.

Figure 5.

Range of variation of inertia matrix elements across all RSQ configurations possible.

These results highlight a significant variation in these parameters. Notably, the off-diagonal terms $ {I}_{xy}$ reach values comparable to the diagonal elements, necessitating its inclusion in both simulation and control design models. In contrast, the terms $ {I}_{xz}$ and $ {I}_{yz}$ are approximately 100 times smaller, justifying their omission to simplify the control design.

To highlight the impact of morphological reconfiguration on inertial properties and flight parameters, Table 4 compares the proposed RSQ with representative rigid multirotor platforms, namely a quadrotor, a hexarotor and an octocopter. The inertia values reported for the RSQ correspond to the principal moments of inertia, computed from the full inertia matrix for different morphological configurations. From Table 4, it can be observed that the RSQ exhibits a significantly wider range of inertia values compared to conventional rigid multirotor platforms. These inertia variations are explicitly incorporated into the proposed LPV model, whose parameters form the matrices presented in Subsection (2). As a result, the eigenvalues of the linearised system vary across configurations, leading to stronger dynamic coupling effects and configuration-dependent stability margins. This behaviour fundamentally differs from conventional multirotor models [Reference Daadi, Boulebtinai, Derrouaoui and Boudjema43–Reference Kose and Oktay45], where the inertia matrix and the corresponding state space parameters are assumed to be constant. Despite its reconfigurable architecture, the RSQ maintains comparable parameters in terms of translational speed, maximum altitude and total mass. This demonstrates that the proposed platform achieves enhanced adaptability without compromising the core flight capabilities expected from multirotor UAVs.

Table 4.

Comparison of inertial and performance parameters for linearised models of multirotor platforms

Table 4 Long description

The table compares inertial and performance parameters for linearized models of multirotor platforms, including RSQ, Quadrotor, Octorotor, and Hexarotor. It has 6 rows and 5 columns. The columns are labeled Eigen value, RSQ, Quadrotor, Octorotor, and Hexarotor. The rows are labeled I_xx [kg.m^2], I_yy [kg.m^2], I_zz [kg.m^2], Velocity [m/s], Altitude_max [m], and Poids [Kg]. Row 1: I_xx [kg.m^2], [0.028 to 0.134], 0.03, 0.009, 0.035. Row 2: I_yy [kg.m^2], [0.01 to 0.667], 0.03, 0.013, 0.035. Row 3: I_zz [kg.m^2], [0.012 to 0.667], 0.032, 0.017, 0.06. Row 4: Velocity [m/s], 0.6, 0.6, 0.8, 0.8. Row 5: Altitude_max [m], 350, 300, 350, 350. Row 6: Poids [Kg], 1.91, 1.79, 1.2, 1.3.

3.1 Dynamic simulation model

The system’s equations of motion are derived to relate the inputs (forces and torques) to the outputs (position and orientation). The model integrates all relevant physical effects, including aerodynamic drag, gyroscopic moments and configuration-dependent inertial variations. It also accounts for the drone’s ability to change its geometry during flight. This subsection is divided into two parts: translational dynamics (forces) and rotational dynamics (moments).

3.1.1 Translation dynamics

The translation dynamics are derived using Newton’s first law,

(8) \begin{align} \begin{array}{lll}\sum {F\!}_{{R}_{I}}=\dfrac{d}{dt}(M{V\!}_{{R}_{I}})& & \\ & & \end{array}\end{align}

where $ {R}_{I}$ denotes the inertial frame, $ {F\!}_{{R}_{I}}$ denotes the external forces applied to the RSQ expressed in $ {R}_{I}$ and $ {V\!}_{{R}_{I}}$ denotes the velocity of the RSQ expressed in the inertial frame.

Gravitational force

The gravitational force $ {F}_{{g\!}_{{R}_{I}}}$ in the inertial frame is given by Equation (9),

(9) \begin{align} {F}_{{g\!}_{{R}_{I}}}={\left[\begin{array}{lll}0& 0& -Mg \end{array}\right]}^{T}\end{align}

where $ g$ denotes the acceleration due to gravity.

Thrust force

For the used rotors, the thrust force relationship is defined as specified in Ref. [Reference Daadi, Boulebtinai, Derrouaoui and Boudjema43],

(10) \begin{align} {F}_{ij}=d\frac{\varpi }{S+\varpi }{u}_{ij}\end{align}

where $ {u}_{ij}$ denotes the pulse width modulation (PWM) input of rotor $ ij$ , $ d$ is a positive gain, and $ \varpi $ represents the actuator bandwidth; both parameters are specified in Table 5.

Table 5.

RSQ parameters

Table 5 Long description

A table with two columns and twelve rows. The first column is labeled 'Parameters' and the second column is labeled 'Value'. The table lists various parameters related to a rotor system and their corresponding values. The parameters include force, angular velocity, gain, revolutions per minute, mass, moments of inertia, voltage, pulse width modulation, and thrust. Each row provides a specific value for each parameter. For example, the first row lists the parameter 'Force' with a value of '12 N', and the second row lists the parameter 'Angular velocity' with a value of '15 rad/sec'. The table provides a comprehensive overview of the rotor system's specifications.

The thrust forces generated by all propellers in the CoG frame $ {R}_{G}$ are given as:

(11) \begin{align} {F}_{{P}_{{R}_{G}}}={[\begin{array}{l@{\hspace*{8pt}}l@{\hspace*{8pt}}l}0& 0& {F}_{T}={F}_{11}+{F}_{12}+{F}_{21}+{F}_{22} \end{array}]}^{T}\end{align}

Drag force

The drag force acting on the RSQ depends not only on its translational velocity but also on its configuration. In this work, the drag is modeled as the sum of the individual drag forces acting on each arm, which vary with the rotation angle of the arms in the horizontal plane.

Assumption: Since reconfiguration occurs in the horizontal plane and the roll/pitch angles remain small, the projected surface area of the arms and base along the vertical axis is approximately constant. Therefore, the drag forces acting on the four arms along the Z-axis are assumed identical and depend only on the translational velocity. The drag force is thus expressed as:

(12) \begin{align} {D\!}_{{R}_{I}}=Diag\left[\begin{array}{lll}{K}_{{\beta }_{x}}+\sum \limits_{ij=11}^{22} {K}_{ar{m}_{ij,x}}\left(\xi \right(t))& & \\[3pt] {K}_{{\beta }_{y}}+\sum \limits_{ij=11}^{22} {K}_{ar{m}_{ij,y}}\left(\xi \right(t))& & \\[3pt] {K}_{{\beta }_{z}}+4{K}_{ar{m}_{ij,\mathrm{z}}}& & \end{array}\right]\dot{\chi }\end{align}

with

\begin{align*} \dot{\chi }=[\dot{x}\,\,\,\, \dot{y}\,\,\,\, \dot{z}]^{T}\end{align*}

where $ {K}_{{\beta }_{x,y,z}}$ denotes the drag coefficients of the central base along the $ {X}_{\beta }$ , $ {Y}_{\beta }$ and $ {Z}_{\beta }$ axes, respectively, and $ {K}_{ar{m}_{ij,x,y,z}}\left(\xi \right(t))$ represents the drag coefficients of the arm $ \mathrm{a}\mathrm{r}{\mathrm{m}}_{ij}$ along the corresponding axes, which depend on the arm’s orientation $ \xi \left(t\right)$ .

According to the previous discussion, all forces applied to the RSQ are expressed in the inertial frame, as follows:

(13) \begin{align} {F\!}_{_{{R}_{I}}}={F}_{{g\!}_{{R}_{I}}}+{F}_{{P\!}_{{R}_{I}}}-{D\!}_{_{{R}_{I}}}\end{align}

From Equation (8), we can write

\begin{align*} {F\!}_{_{{R}_{I}}}={F}_{_{{R}_{G}}}=M\frac{d}{dt}({V\!}_{{R}_{I}})\end{align*}

Therefore, Equation (8) becomes

(14) \begin{align} \left\{\begin{array}{lll}\ddot{x}=\dfrac{1}{M}{F}_{{P}_{{R}_{G}}}\left(C\right(\psi \left)S\right(\theta \left)C\right(\phi )+S(\psi \left)S\right(\phi ))-\dfrac{{K}_{{\beta }_{x}}+\sum _{ij=11}^{22} {K}_{ar{m}_{ij,x}}\left(\xi \right(t))}{M}\dot{x}& & \\[8pt] \ddot{y}=\dfrac{1}{M}{F}_{{P}_{{R}_{G}}}\left(S\right(\psi \left)S\right(\theta \left)C\right(\phi )-C(\psi \left)S\right(\phi ))-\dfrac{{K}_{{\beta }_{y}}+\sum _{ij=11}^{22} {K}_{ar{m}_{ij,y}}\left(\xi \right(t))}{M}\dot{y}& & \\[8pt] \ddot{z}=\dfrac{1}{M}{F}_{{P}_{{R}_{G}}}\left(C\right(\theta \left)C\right(\phi ))-\dfrac{{K}_{{\beta }_{z}}+4{K}_{ar{m}_{ij,z}}\left(\xi \right(t))}{M}\dot{z}-g \end{array}\right.\end{align}

where $ [x \quad y \quad z]^{T}$ denotes the position of the drone expressed in the inertial frame $ {R}_{I}$ .

3.1.2 Rotation dynamics

The reconfiguration of the drone affects its inertia matrix. Consequently, considering these variations is crucial for precise modeling. By applying Newton’s second law, the equations of motion can be expressed as follows:

(15) \begin{align} \begin{array}{lll}{M}_{{R}_{G}}=I\left({\xi }_{ij}\right(t))\dot{\Omega }+\Omega \times I\left({\xi }_{ij}\right(t))\Omega +\dfrac{dI({\xi }_{ij}(t))}{dt}\Omega & & \end{array}\end{align}

where $ {M}_{{R}_{G}}$ includes all torques acting on the drone, expressed in the body frame $ {R}_{G}$ and $ \Omega $ denotes the angular velocity, also expressed in $ {R}_{G}$ .

External torques

This gives after development:

(16) \begin{align} \left\{\begin{array}{lll}{M}_{\phi }=D{Y}_{_{11}}{F}_{11}+D{Y}_{_{12}}{F}_{12}+D{Y}_{_{21}}{F}_{21}+D{Y}_{_{22}}{F}_{22}& & \\ {M}_{\theta }=-(D{X}_{_{11}}{F}_{11}+D{X}_{_{12}}{F}_{12}+D{X}_{_{21}}{F}_{21}+D{X}_{_{22}}{F}_{22})& & \\ {M}_{\psi }={\tau }_{11}-{\tau }_{12}+{\tau }_{22}-{\tau }_{21}& & \end{array}\right.\end{align}

where $ {M}_{\theta }$ , $ {M}_{\phi }$ and $ {M}_{\psi }$ denote the moments around the $ y$ , $ x$ and $ z$ axes of the body frame $ {R}_{G}$ , respectively; $ D{X}_{ij}$ and $ D{Y}_{ij}$ represent the $ x$ and $ y$ coordinates of the point of application of the force $ {F}_{ij}$ in $ {R}_{G}$ , as illustrated in Fig. 3; and $ {\tau }_{ij}$ denotes the drag torque generated by propeller $ ij$ , as given in Ref. [Reference Daadi, Boulebtinai, Derrouaoui and Boudjema43],

(17) \begin{align} {\tau }_{ij}={K}_{z}{u}_{ij}\end{align}

where $ {K}_{z}$ represents a positive gain, defined in Table 5.

Gyroscopic moments

The gyroscopic moments induced by the rotors resist changes in the drone’s attitude [Reference Daadi, Boulebtinai, Derrouaoui and Boudjema43]. These moments are expressed by Equation (18):

(18) \begin{align} \begin{array}{ll}{M}_{gyr}& =\Omega \times {J}_{rz}{\left[\begin{array}{lll}0& 0& ({u}_{11}-{u}_{12}-{u}_{21}+{u}_{22}){K}_{v}V\end{array}\right]}^{T}={\left[\begin{array}{lll}{U}_{g}q& -{U}_{g}p& 0\end{array}\right]}^{T}\end{array}\end{align}

where $ {U}_{g}=\left({u}_{11}-{u}_{12}-{u}_{21}+{u}_{22}\right){J}_{rz}{K}_{v}V$ , $ {M}_{gyr}$ denotes the vector of gyroscopic moments, $ {J}_{rz}$ represents the moment of inertia of a rotor about its rotation axis ( $ Z$ -axis), $ {\mathrm{K}}_{v}$ is a positive constant, $ V$ denotes the input voltage and $ [p \, q \, r]^{T}$ represents the angular velocity of the drone about the axes of the body frame $ {R}_{G}$ .

The RSQ parameters are presented in Table 5.

After development, the rotation dynamics of the RSQ becomes

(19) \begin{align} \left\{\begin{array}{ll}\dot{p}=& {A}_{\left(p\right)}{p}^{2}+{B}_{\left(p\right)}pq-{C}_{\left(p\right)}pr-{D}_{\left(p\right)}p-{E}_{\left(p\right)}{q}^{2}+{F}_{\left(p\right)}qr\\[4pt] & -{G}_{\left(p\right)}q+{H}_{\left(p\right)}{r}^{2}-{I}_{\left(p\right)}r+{Z}_{\left(p\right)}\\[4pt] \dot{q}=& {A}_{\left(q\right)}{p}^{2}-{B}_{\left(q\right)}pq-{C}_{\left(q\right)}pr-{D}_{\left(q\right)}p-{E}_{\left(q\right)}{q}^{2}+{F}_{\left(q\right)}qr\\[4pt] & -{G}_{\left(q\right)}q-{H}_{\left(q\right)}{r}^{2}-{I}_{\left(q\right)}r+{Z}_{\left(q\right)}\\[4pt] \dot{r}=& -{A}_{\left(r\right)}{p}^{2}+{B}_{\left(r\right)}pq+{C}_{\left(r\right)}pr-{D}_{\left(r\right)}p+{E}_{\left(r\right)}{q}^{2}-{F}_{\left(r\right)}\mathrm{q}r\\[4pt] & -{G}_{\left(r\right)}q+{H}_{\left(r\right)}{r}^{2}-{I}_{\left(r\right)}r+{Z}_{\left(r\right)}\end{array}\right.\end{align}

The time derivatives of the Euler angles are expressed in terms of the angular rates in the body frame as follows:

(20) \begin{align}\left[ \begin{array}{l} {\dot \phi } \\[3pt] {\dot \theta } \\[3pt] {\dot \psi } \end{array} \right] = \overbrace {\left[ \begin{array}{c@{\quad}c@{\quad}c} 1 & {Sin\left( \phi \right)Tan\left( \theta \right)} & {Cos\left( \phi \right)Tan\left( \theta \right)} \\[3pt] 0 & {Cos\left( \phi \right)} & { - Sin\left( \phi \right)} \\[3pt] 0 & {{{Sin\left( \phi \right)} \over {Cos\left( \theta \right)}}} & {{{Cos\left( \phi \right)} \over {Cos\left( \theta \right)}}} \\[3pt] \end{array} \right]}^\Im \cdot \left[ \begin{array}{l} p \\[3pt] q \\[3pt] r \end{array} \right]\end{align}

All these parameters are detailed in Appendix 1.

The complexity of this model arises from the reconfiguration and asymmetry of the structure, which prevent the principal axes of the RSQ from coinciding with its rotation axes. As a result, the inertia matrix becomes non-diagonal, leading to coupling terms between the system’s angular accelerations. Consequently, a rotation about the $ X$ -axis can induce unintended rotations around the $ Y$ or $ Z$ -axes, and vice versa.

3.2 Dynamic control model

The RSQ is primarily designed for object grasping and transportation, as well as for navigating through narrow spaces by exploiting its horizontal reconfiguration capability. The platform is intended to operate at low speeds in cluttered environments to minimise collision risks with the surrounding obstacles.

To facilitate the synthesis of control laws, the dynamic model in Equation (19) is simplified by adopting the following assumptions:

• • The inertia cross terms $ {I}_{xz}$ and $ {I}_{yz}$ are neglected, as they are negligible compared to the main inertia components. This simplification is supported by the inertial variation analysis shown in Fig. 5 (and confirmed in Fig. 15).

• • The roll ( $ \phi $ ) and pitch ( $ \theta $ ) angles remain small during flight.

Under these assumptions, the transformation matrix reduces to:

(21) \begin{align} \mathfrak{I}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(\mathrm{1,1},1)\end{align}

Consequently, the rotational dynamics of the RSQ in the inertial frame can be expressed as:

(22) \begin{align} \left\{\begin{array}{ll}\ddot{\phi }=& \dfrac{{I}_{xy}{M}_{\theta }-{I}_{yy}{M}_{\phi }+({I}_{yy}{\skew3\dot{I}}_{xx}-{I}_{xy}{\skew3\dot{I}}_{xy}-{I}_{xy}{U}_{g})\dot{\phi }}{{I}_{xy}^{2}-{I}_{xx}{I}_{yy}}+\dfrac{({I}_{yy}{I}_{zz}-{I}_{xy}^{2}-{I}_{yy}^{2})\dot{\theta }\dot{\psi }+({I}_{xy}{I}_{zz}-{I}_{xy}{I}_{yy}-{I}_{xx}{I}_{xy})\dot{\phi }\dot{\psi }}{{I}_{xy}^{2}-{I}_{xx}{I}_{yy}}\\[5pt] & +\dfrac{({I}_{yy}{\skew3\dot{I}}_{xy}-{I}_{xy}{\skew3\dot{I}}_{yy}-{I}_{yy}{U}_{g})\dot{\theta }}{{I}_{xy}^{2}-{I}_{xx}{I}_{yy}}\\[5pt] \ddot{\theta }=& \dfrac{{I}_{xy}{M}_{\phi }-{I}_{xx}{M}_{\theta }+({I}_{xx}{\skew3\dot{I}}_{xy}-{I}_{xy}{\skew3\dot{I}}_{xx}+{I}_{xx}{U}_{g})\dot{\phi }}{{I}_{xy}^{2}-{I}_{xx}{I}_{yy}}+\dfrac{({I}_{xx}^{2}+{I}_{xy}^{2}-{I}_{xx}{I}_{zz})\dot{\phi }\dot{\psi }+({I}_{xx}{I}_{xy}+{I}_{xy}{I}_{yy}-{I}_{xy}{I}_{zz})\dot{\theta }\dot{\psi }}{{I}_{xy}^{2}-{I}_{xx}{I}_{yy}}\\[5pt] & +\dfrac{({I}_{xx}{\skew3\dot{I}}_{yy}-{I}_{xy}{\skew3\dot{I}}_{xy}+{I}_{xy}{U}_{g})\dot{\theta }}{{I}_{x \mathrm{y}}^{2}-{I}_{xx}{I}_{yy}}\\[5pt] \ddot{\psi }=& \dfrac{{M}_{\psi }-{\skew3\dot{I}}_{zz}\dot{\psi }-{I}_{xy}({\dot{\phi }}^{2}-{\dot{\theta }}^{2})+({I}_{xx}-{I}_{yy})\dot{\phi }\dot{\theta }}{{I}_{zz}}\end{array}\right.\end{align}

The rotational dynamic model highlights the strong coupling between rotational motions through the necessity of two control inputs ${M}_{\phi }$ and $ {M}_{\theta }$ for controlling the roll or pitch. As a result, regulating the roll angle inherently requires the simultaneous adjustment of both roll and pitch angles. This interdependence complicates the synthesis of control laws, making it challenging to address each rotational axis independently. Moreover, this phenomenon can significantly impact system stability, making certain manoeuvres more difficult to control. This model is distinct from conventional drone models, which typically require only a single control input to achieve a desired motion.

The nonlinear attitude model in Equation (22), characterised by time-varying parameters, can be rewritten in the form of an LPV system as follows:

(23) \begin{align} \left\{\begin{array}{ll} \dot{x}=A(I(t),\skew3\dot{I}(t))x+B(I(t))u+{B}_{w}(I(t),\skew3\dot{I}(t))w\\ y=Cx+{B}_{\mathrm{y}u}(I(t))u+{B}_{wu}(I(t),\skew3\dot{I}(t))w\end{array}\right.\end{align}

The attitude state vector is defined as follows:

\begin{align*} x={[\begin{array}{l}\phi \quad \theta \quad \psi \quad \dot{\phi } \quad \dot{\theta } \quad \dot{\psi }\end{array}]}^{T}\in {\mathbb{R}}^{6}\end{align*}

with the control input vector defined as:

\begin{align*} u={[\begin{array}{l}{V}_{\phi} \quad {V}_{\theta } \quad {M}_{\psi }\end{array}]}^{T}\in {\mathbb{R}}^{3}\end{align*}

with $ {V}_{\phi }$ and $ {V}_{\theta }$ expressed as functions of the control moments $ {M}_{\phi }$ et $ {M}_{\theta }$ , and the inertia parameters I(t), as given in Equation (29).

In the proposed LPV formulation, the scheduling parameters are defined as the configuration-dependent inertial properties and their time derivatives, i.e.

\begin{align*} \left[\begin{array}{l}I\left(t\right)\\ \skew3\dot{I}\left(t\right)\end{array}\right]=\left[\begin{array}{l}{I}_{xx}\left(t\right) \quad {I}_{yy}\left(t\right) \quad {I}_{zz}\left(t\right) \quad {I}_{xy}\left(t\right)\\ {\skew3\dot{I}}_{xx}\left(t\right) \quad {\skew3\dot{I}}_{yy}\left(t\right) \quad {\skew3\dot{I}}_{xy}\left(t\right) \quad {\skew3\dot{I}}_{zz}\left(t\right)\end{array}\right]\end{align*}

where the inertia matrix and its time derivatives are considered within bounded sets $ I\left(t\right)\in \Pi$ and $ \skew3\dot{I}\left(t\right)\in \Pi^{\prime}$

where

(24) \begin{align} \left\{\begin{array}{ll}A(I(t),\skew3\dot{I}(t))& =\left[\begin{array}{l@{\quad}l}{0}_{3\times 3}& {\mathbf{I}}_{3}\\ {0}_{3\times 3}& {A}_{22}(I(t),\skew3\dot{I}(t))\end{array}\right]\!,\\[15pt] B(I(t))& =\left[\begin{array}{l}{0}_{3\times 3}\\ {B}_{2}(I(t))\end{array}\right]\!,\end{array}\right.\end{align}

with

(25) \begin{align} {A}_{22}(I(t),\skew3\dot{I}(t))=&\left[\begin{array}{l@{\quad}l@{\quad}l}{A}_{11}& {A}_{12}& {A}_{13}\\ {A}_{21}& {A}_{22}& {A}_{23}\\ {A}_{31}& {A}_{32}& {A}_{33}\end{array}\right]\!,{B}_{2}(I(t)) = \left[\begin{array}{l@{\quad}l@{\quad}l}\dfrac{1}{\Delta }& 0& 0\\ 0& \dfrac{1}{\Delta }& 0\\ 0& 0& \dfrac{1}{{I}_{zz}}\end{array}\right]\!,\nonumber\\ &\left[\begin{array}{l}{A}_{11}=\dfrac{{I}_{yy}{\skew3\dot{I}}_{xx}-{I}_{xy}{\skew3\dot{I}}_{xy}-{I}_{xy}{U}_{g}}{\Delta }+\dfrac{{I}_{xy}{I}_{zz}-{I}_{xy}{I}_{yy}-{I}_{xx}{I}_{xy}}{2 \Delta }\\[8pt] {A}_{12}=\dfrac{{I}_{yy}{\skew3\dot{I}}_{xy}-{I}_{xy}{\skew3\dot{I}}_{yy}-{I}_{yy}{U}_{g}}{{\Delta}}+\dfrac{{I}_{yy}{I}_{zz}-{I}_{xy}^{2}-{I}_{yy}^{2}}{2 \Delta }\\[8pt] {A}_{13}=\dfrac{{I}_{xy}{I}_{zz}-{I}_{xy}{I}_{yy}-{I}_{xx}{I}_{xy}}{2 \Delta }+\dfrac{{I}_{yy}{I}_{zz}-{I}_{xy}^{2}-{I}_{yy}^{2}}{2 \Delta }\\[8pt] {A}_{21}=\dfrac{{I}_{xx}{\skew3\dot{I}}_{xy}-{I}_{x \mathrm{y}}{\skew3\dot{I}}_{xx}+{I}_{xx}{U}_{g}}{\Delta }+\dfrac{{I}_{xx}^{2}+{I}_{xy}^{2}-{I}_{xx}{I}_{xz}}{2 \Delta }\\[8pt] {A}_{22}=\dfrac{{I}_{xx}{\skew3\dot{I}}_{yy}-{I}_{xy}{\skew3\dot{I}}_{xy}+{I}_{xy}{U}_{g}}{\Delta }+\dfrac{{I}_{xx}{I}_{xy}+{I}_{xy}{I}_{yy}-{I}_{xy}{I}_{zz}}{2 \Delta }\\[8pt] {A}_{23}=\dfrac{{I}_{xx}^{2}+{I}_{xy}^{2}-{I}_{xx}{I}_{xz}}{2 \Delta }+\dfrac{{I}_{xx}{I}_{xy}+{I}_{xy}{I}_{yy}-{I}_{xy}{I}_{zz}}{2 \Delta }\\[8pt] {A}_{31}=\dfrac{{I}_{xx}-{I}_{yy}}{2{I}_{zz}}\\[8pt] {A}_{32}=\dfrac{{I}_{xx}-{I}_{yy}}{2{I}_{zz}}\\[8pt] {A}_{33}=-\dfrac{{\skew3\dot{I}}_{zz}}{{I}_{zz}}\end{array}\right]\end{align}

and

\begin{align*} \Delta (I(t))={I}_{xy}^{2}-{I}_{xx}{I}_{yy}.\end{align*}

The performance output is selected as:

(26) \begin{align} y=Cx,\quad C=\left[\begin{array}{l@{\quad}l}{I}_{3}& {0}_{3\times 3}\end{array}\right]\!,\,\,{B}_{yu}={B}_{wu}={0}_{3\times 1}.\end{align}

The thrust forces and torques exerted by the drone are related to the motor commands (PWM signals) through the control allocation matrix, defined for the RSQ as follows:

(27) \begin{align} \left[\begin{array}{l}{F}_{P}\\[5pt] {M}_{\theta }\\[5pt] {M}_{\phi }\\[5pt] {M}_{\psi }\end{array}\right]={\left[\begin{array}{l@{\quad}l@{\quad}l@{\quad}l}\unicode{x1E11}& \unicode{x1E11}.D{Y}_{22}({\xi }_{ij}(t))& -\unicode{x1E11}.D{X}_{22}({\xi }_{ij}(t))& {K}_{z}\\[5pt] \unicode{x1E11}& \unicode{x1E11}.D{Y}_{21}({\xi }_{ij}(t))& -\unicode{x1E11}.D{X}_{21}({\xi }_{ij}(t))& -{K}_{z}\\[5pt] \unicode{x1E11}& \unicode{x1E11}.D{Y}_{11}({\xi }_{ij}(t))& -\unicode{x1E11}.D{X}_{11}({\xi }_{ij}(t))& {K}_{z}\\[5pt] \unicode{x1E11}& \unicode{x1E11}.D{Y}_{12}({\xi }_{ij}(t))& -\unicode{x1E11}.D{X}_{12}({\xi }_{ij}(t))& -{K}_{z}\end{array}\right]}^{T}\left[\begin{array}{l}{u}_{22}\\[5pt] {u}_{21}\\[5pt] {u}_{11}\\[5pt] {u}_{12}\end{array}\right]\end{align}

The elements of the control allocation matrix vary according to the positions of the rotors relative to $ {O}_{G}$ . Therefore, a real-time estimation of these parameters is required to ensure an accurate mapping between the rotor PWM and the desired moments and thrust.

Figure 6.

Singular configurations for the RSQ.

Figure 6 Long description

A diagram illustrating different configurations of a multirotor unmanned aerial vehicle. Panel A: A straight, rigid frame multirotor UAV with four propellers. Panel B: A curved configuration of the same UAV with four propellers. Panel C: A diagonal configuration of the UAV with four propellers. The diagram includes mathematical relationships between the configurations.

After a thorough analysis of the control allocation matrix it was concluded that the RSQ presents singular configurations, shown in Fig. 6. For instance, a configuration where all four motors are aligned, i.e $ {\xi }_{12}={\xi }_{22}=-2\times {\xi }_{11}=-2\times {\xi }_{12}$ , renders the drone uncontrollable, reducing the available DoF for attitude control from four to three. Another critical singular configuration, characterised by symmetry with respect to the $ {O}_{\beta }$ , impairs yaw control while leaving pitch and roll unaffected. This occurs when $ {\xi }_{12}=-{\xi }_{22}$ and $ {\xi }_{11}=-{\xi }_{21}$ . These configurations must be carefully avoided when planning transitions between different configurations to ensure the controllability of the drone.

4.1 Stability analysis

In this study, the LPV framework (Equation 23) is employed to analyse and control the nonlinear dynamics of the RSQ. The LPV approach allows representing the nonlinear or time-varying system as a family of linear systems whose dynamics depend on measurable scheduling parameters, such as the time-varying inertia matrix and its derivative. This formulation facilitates both stability analysis and controller synthesis under parameter variations.

Consider the following LPV representation of the closed-loop system:

(30) \begin{align} \left\{\!\!\begin{array}{ll}& \dot{x}\left(t\right)=A(I(t),\skew3\dot{I}(t))x\left(t\right)+B(I(t))u\left(t\right)+{B}_{w}(I(t),\skew3\dot{I}(t))W\left(t\right)\\[3pt] & y\left(t\right)={C}_{y}(I(t),\skew3\dot{I}(t))x\left(t\right)+{D}_{\mathrm{y}u}(I(t),\skew3\dot{I}(t))u\left(t\right)+{D}_{yw}(I(t),\skew3\dot{I}(t))W\left(t\right)\\[3pt] & {z}_{\infty }\left(t\right)={C}_{z}(I(t),\skew3\dot{I}(t))x\left(t\right)+{D}_{zu}(I(t),\skew3\dot{I}(t))u\left(t\right)+{D}_{zw}(I(t),\skew3\dot{I}(t))W\left(t\right)\end{array}\right.\end{align}

where $ x\left(t\right)\in {\mathbb{R}}^{12}$ denotes the state vector, $ u\left(t\right)\in {\mathbb{R}}^{4}$ is the control input and $ w\left(t\right)\in {\mathbb{R}}^{4}$ represents the disturbance input. The measured output and the controlled output are denoted by $ y\left(t\right)\in {\mathbb{R}}^{6}$ and $ z\left(t\right)\in {\mathbb{R}}^{6}$ , respectively.

The matrices $ A(\cdot )$ , $ B(\cdot )$ , $ {B}_{w}(\cdot )$ , $ {C}_{y}(\cdot )$ , $ {C}_{z}(\cdot )$ , $ {D}_{yu}(\cdot )$ , $ {D}_{uw}(\cdot )$ , $ {D}_{zu}(\cdot )$ and $ {D}_{zw}(\cdot )$ are continuous matrix-valued functions of the scheduling parameter $ I\left(t\right)$ and $ \skew3\dot{I}\left(t\right)$ . The parameter vector $ I\left(t\right)\in \Pi$ and its derivative $ \skew3\dot{I}\left(t\right)\in \Pi{^{\prime}}$ belong to compact hyperrectangles $ \Pi$ and $ \Pi{^{\prime}}$ , respectively.

Theorem 1. [Reference Mattei46, Reference Song and Yang47]: Consider system (30) with $D_{yu}(I(t),\skew3\dot{I}(t))$ a null matrix for all $ I\left(t\right)\in \Pi$ and $ \skew3\dot{I}\left(t\right)\in {\Pi}^{{^{\prime}}}$ . The system is quadratically stabilisable (QS) via static output feedback, with a prescribed convergence rate $ \alpha \gt 0$ and an $ {H}_{\infty }$ norm bound less than a positive scalar $ \gamma $ , if there exist a symmetric positive definite matrix $ P={P}^{{\rm T}}\gt 0$ and a static output feedback gain $ {K}_{s}$ such that the following linear matrix inequality (LMI) holds for all $ I\left(t\right)\in \Pi$ and $ \skew3\dot{I}\left(t\right)\in {\Pi}^{{^{\prime}}}$ :

(31) \begin{align} \left[\begin{array}{c@{\quad}c@{\quad}c}({A}_{K}{)}^{{\rm T}}P+P{A}_{K}& P({B}_{w}+B{K}_{s}{D}_{yw})& {\gamma }^{-1}({C}_{z}+{D}_{zu}{K}_{s}{C}_{y}{)}^{{\rm T}}\\ {\star}& -I& {\gamma }^{-1}{D}_{zw}^{{\rm T}}\\ {\star}& {\star}& -I\end{array}\right]\lt 0,\end{align}

where

\begin{align*}{A_K} = \tilde A + B{K_s}{C_y},\,\,\,\,{\text {and}}\,\,\,\,\tilde A = A + \beta I\end{align*}

Throughout the paper, we adopt the following standard state-space representation of a PID controller with state augmentation:

(32) \begin{align} \left\{\begin{array}{lll}{\dot{x}}_{I}={G}_{I}y& & \\ u={K}_{I}{x}_{I}+{K}_{P}y+{K}_{D}\dot{y}& & \end{array}\right.\end{align}

where $ {x}_{I}\in {\mathbb{R}}^{6}$ denotes the integral state vector and $ G\in {\mathbb{R}}^{6\times 6}$ is a constant matrix

Theorem 2. [Reference Mattei46, Reference Song and Yang47]: Consider system (30) with $ {D}_{yu}(I(t),\skew3\dot{I}(t))$ and $ {D}_{yw}(I(t),\skew3\dot{I}(t))$ a null matrix for all $ I\left(t\right)\in \Pi$ , and assume that the parameter variation rate $ \skew3\dot{I}\left(t\right)$ is bounded and belongs to a compact hyperrectangle $ \Pi{^{\prime}}$ . Then, there exists a PID controller of the form (32) that guarantees QS with a prescribed convergence rate $ \beta $ and an $ {\mathcal{H}}_{\infty }$ norm bound less than a positive scalar $ {\gamma }_{b}$ , via static output feedback if there exist a symmetric positive definite matrix $ P$ and a gain matrix

\begin{align*} {K}_{s}=\left[\begin{array}{lll}{K}_{1}& {K}_{2}& {K}_{3}\end{array}\right],\,\,\,{K}_{1},\,\,\,{K}_{3}\in {\mathbb{R}}^{{n}_{u}\times {n}_{y}},\,\,\,{K}_{2}\in {\mathbb{R}}^{{n}_{u}\times {n}_{x}}\end{align*}

such that the following LMI holds:

(33) \begin{align}\left[ \begin{array}{c@{\quad}c@{\quad}c} {{{({A_k})}^{\rm T}}P + P\left( {{A_k}} \right)} & {P\left( {{{\tilde B}_w} + \tilde B{K_y}{{\tilde D}_{yw}}} \right)} & {\gamma _b^{ - 1}{{({C_z} + {D_{zu}}{K_y}{C_y})}^{\rm{{\rm T}}}}} \\ {\rm{{\star}}} & { - I} & {\gamma _b^{ - 1}\tilde D_{zw}^{\rm T}} \\ {\rm{{\star}}} & {\rm{{\star}}} & { - I} \end{array} \right] \lt 0,\;\;\;\;\forall \;\;\;\;I\left( t \right) \in {\rm{\Pi \;\;\;\;\mathrm{and}\;\;\;\;}}\dot I\left( t \right) \in {\rm{\Pi '}}\end{align}

with

\begin{align*}{A_k} = \hat A + \tilde B{K_s}{C_y},{\rm{\;\;\;\;\;\;\;\;\mathrm{and}\;\;\;\;\hat A}} = \hat A + \beta I\end{align*}

where

(34) \begin{align}\left\{\begin{array}{l}{\hat A\left( {I\left( t \right),\dot I\left( t \right)} \right) = \tilde A\left( {I\left( t \right),\dot I\left( t \right)} \right) + \beta I,{\rm{\;\;\;\;}}\hat A\left( p \right) = \left[ \begin{array}{c} {A\left( {I\left( t \right),\dot I\left( t \right)} \right)0} \\[6pt] {{G_I}{C_y}\left( {I\left( t \right),\dot I\left( t \right)} \right)0} \\[6pt] \end{array} \right],{\rm{\;\;\;\;}}\tilde B\left( {I\left( t \right)} \right) = \left[ \begin{array}{c} {B\left( {I\left( t \right)} \right)} \\[6pt] 0 \\ \end{array} \right]}, \\[20pt] {{{\tilde B}_w}\left( p \right) = \left[ \begin{array} {c}{{B_w}\left( p \right)} \\[6pt] 0 \\ \end{array} \right],{\rm{\;\;\;\;}}{{{\rm{\tilde C}}}_y}\left( {p,\dot p} \right) = \left[ \begin{array} {c} {{C_y}\left( p \right)0} \\[6pt] {{C_y}\left( p \right) + \dot p{\dfrac{\partial {C_y}} {\partial p}}A\left( p \right)} I \\[6pt] \end{array} \right],{\rm{\;\;\;\;}}{{\tilde C}_z}\left( p \right) = \left[ \begin{array}{c} {{C_z}\left( p \right)} \ \ 0 \end{array} \right]} \\[20pt] {{\tilde D}_{yw}}\left( p \right) = \left[ \begin{array} {c} 0 \\[6pt] {{C_y}\left( p \right){B_w}\left( p \right)} \\[6pt] \end{array} \right] \end{array} \right.\end{align}

Under this assumption, a feasible set of PID gains possibly scheduled with respect to the parameter vector $ p$ is given by:

(35) \begin{align} \left\{\begin{array}{ll}{K}_{p}\left(p\right)& =\left(I-{K}_{3}(I+{C}_{y}(p)B(p){K}_{3}{)}^{-1}{C}_{y}(p)B(p)\right){K}_{1}\\[5pt] {K}_{i}\left(p\right)& =\left(I-{K}_{3}(I+{C}_{y}(p)B(p){K}_{3}{)}^{-1}{C}_{y}(p)B(p)\right){K}_{2}\\[5pt] {K}_{d}\left(p\right)& ={K}_{3}(I+{C}_{y}(p\left)B\right(p){K}_{3}{)}^{-1}\end{array}\right.\end{align}

The invertibility condition $ I-{C}_{y}\left(p\right)B\left(p\right){K}_{3}$ must hold for all $ p\in \Pi$ to ensure a valid mapping. Based on Theorems 1 and 2 [Reference Mattei46, Reference Song and Yang47], the LPV system Equation (23) is quadratically stabilisable and satisfies an $ {\mathcal{H}}_{\infty }$ performance bound $ \gamma $ if there exists a symmetric positive definite matrix P such that the LMI Equation (33) holds for all admissible $ p\in \Pi\mathrm{a}\mathrm{n}\mathrm{d}\dot{p}\in {\Pi}^{{^{\prime}}}$ . Moreover, the resulting PID controller guarantees the prescribed convergence rate and ensures robust performance across the entire set of morphological configurations of the RSQ.

By computing $ \dot{y}$ , and noting that the matrix $ C$ does not depend on the scheduling parameters $ I\left(t\right)\mathrm{a}\mathrm{n}\mathrm{d}\skew3\dot{I}\left(t\right)$ , we obtain

(36) \begin{align} \dot{y}=C\dot{x}\end{align}

Substituting the LPV system dynamics Equation (23) into Equation (36) yields

(37) \begin{align} \dot{y}=CA(I(t),\skew3\dot{I}(t))x+CB(I(t))u+C{B}_{w}(I(t),\skew3\dot{I}(t))w\end{align}

Injecting this expression into the PID control law in Equation (32) leads to:

(38) \begin{align} u={K}_{I}{x}_{I}+{K}_{P}Cx+{K}_{D}(CAx+CBu+C{B}_{w}w)\end{align}

Assuming that the matrix $ (I-{K}_{D}CB(I\left(t\right)))$ is invertible for all admissible parameters, the control input can be written as follows:

(39) \begin{align} u=(I-{K}_{D}CB{)}^{-1}[{K}_{I}{x}_{I}+{K}_{P}Cx+{K}_{D}CAx+{K}_{D}C{B}_{w}w].\end{align}

We have the following closed-loop system equations (see Equation (34)) for the definition of the new matrices:

(40) \begin{align}\left\{ {\begin{matrix} {{{\dot x}_a}} = \big( {\tilde A + \tilde B{K_s}\tilde C} \big){x_a} + \left( {{{\tilde B}_w} + \tilde B{K_s}{{\tilde D}_{yw}}} \right)w \\[5pt] z = \big( {{{\tilde C}_z} + {D_{zu}}{K_s}C} \big)x + \left( {{D_{zw}}{ + _{zu}}{K_s}{{\tilde D}_{yw}}} \right)w \end{matrix} } \right.\end{align}

where $ {x}_{a}$ denotes the augmented state vector defined as

(41) \begin{align} {x}_{a}={\left[\begin{array}{l}x\quad {x}_{I}\end{array}\right]}^{T},\quad {\dot{x}}_{I}={G}_{I}Cx\end{align}

It turns out that the synthesis of the PID can be converted into a static output feedback problem on the fictitious system.

The correspondence between the constant matrix gains $ {K}_{1},{K}_{2},{K}_{3}$ and the parameter-dependent matrix-valued functions $ {K}_{p}\left(p\right),{K}_{I}\left(p\right),{K}_{d}\left(p\right)$ assured by the invertibility of matrix $ I-CB(I(t)){K}_{3},\forall I\left(t\right)\in \Pi$ [Reference Mattei46] where

\begin{align*} {K}_{s}=\left[\begin{array}{l}{K}_{1}\\ {K}_{2}\\ {K}_{3}\end{array}\right]=\left[\begin{array}{l}(I-{K}_{3}{C}_{I}B{)}^{-1}{K}_{p}\\ (I-{K}_{3}{C}_{I}B{)}^{-1}{K}_{I}\\ (I-{K}_{3}{C}_{I}B{)}^{-1}{K}_{d}\end{array}\right].\end{align*}

The resulting augmented system with the enlarged parameter vector $ ({I}^{{\rm T}},{\skew3\dot{I}}^{{\rm T}}{)}^{{\rm T}}\in \Pi\times \Pi{^{\prime}}$ guarantees exponential stability of the closed-loop system. The proof follows directly from the equivalence between the PID control law and a static output feedback acting on the augmented LPV representation.

To guarantee the well posedness of the PID control law, the invertibility of the matrix $ \left(I-CB(I(t)){K}_{3}\right)$ is imposed for all admissible inertia matrices $ I\left(t\right)\in \Pi$ . A sufficient condition ensuring this property is given by

(42) \begin{align} \left(I-CB(I(t)){K}_{3}\right){\left(I-CB(I(t)){K}_{3}\right)}^{{\rm T}}\succ 0,\quad \forall I\left(t\right)\in \Pi.\end{align}

Under this condition, the parameter-dependent PID gains remain well defined and continuous over the entire admissible parameter set.

In this work, the feedback gains $ {K}_{s}$ are tuned using an optimisation algorithm (WOA) to explore different combinations that strictly satisfy the closed-loop LPV system stability while enforcing an $ {\mathcal{H}}_{\infty }$ performance bound specified by $ \gamma =1$ [Reference Mattei46]. The optimisation procedure ensures that the LMIs in Equation (33), as stated in Theorem 2, are strictly satisfied. Once these conditions are met, the corresponding PID controller gains are computed from the optimised feedback gains using Equation (35), while guaranteeing the invertibility of the matrix $ (I-{K}_{D}CB(I(t)))$ through verification of condition Equation (42). This ensures that the parameter-dependent PID gains remain well-defined across all morphological configurations. For each feasible combination of parameters, the closed-loop trajectory is simulated to evaluate tracking performance. The integral of squared error (ISE) is computed for each case, and the optimisation iterates until the set of gains that simultaneously guarantees closed-loop stability and achieves the best possible trajectory tracking performance is identified.

The resulting optimal PID gains are summarised in Table 6, and the proposed architecture for PID parameter tuning via feedback gains $ {K}_{s}$ using WOA is illustrated in Fig. 8. This approach ensures uniform stability margins and optimal tracking performance across all morphological configurations of the reconfigurable quadrotor.

(43) \begin{align} \left\{\begin{array}{lll}{F}_{ob\mathrm{j}}=& ISE={\int }_{0}^{+\infty } {e}^{2}\left(t\right)dt=\sum {e}_{\psi }{e}_{\psi }^{T}+{e}_{\theta }{e}_{\theta }^{T}+{e}_{\psi }{e}_{\psi }^{T}+{e}_{x}{e}_{x}^{T}+{e}_{y}{e}_{y}^{T}+{e}_{z}{e}_{z}^{T}& \end{array}\right.\end{align}

where

\begin{align*} {e}_{i}\left(t\right)={i}_{d}\left(t\right)-i\left(t\right)\end{align*}

with $ {e}_{i}(t)$ representing the error of the RMQ, and ( $ i\in \{\psi ,\phi ,\theta ,x,y,\mathrm{a}\mathrm{n}\mathrm{d}\,z\}$ ).

Table 6.

The optimised PID controller gains

Table 6 Long description

A table comparing PID controller gains across different parameters. The table has three rows and six columns. The columns are labeled Gain, F sub p, U sub x, U sub y, M sub psi, V sub psi, and V sub theta. The rows are labeled Kp, Ki, and Kd. Row 1: Kp, 20, 0.26480, 0.27010, 1.94553, 2.46799, 10. Row 2: Ki, 0.017, 0, 0, 5.07227, 1.64702, 2.51823. Row 3: Kd, 20, 0.23808, 0.24411, 10, 0.71995, 1.98547.

Figure 8.

Proposed architecture for PID parameter tuning via WOA.

The search space is defined as:

\begin{align*} {k}_{{s}_{i}}=[{K}_{{1}_{i}},{K}_{{2}_{i}},{K}_{{3}_{i}}]\in [-\mathrm{20,\,20}]\end{align*}

where $ {K}_{{1}_{i}}$ , $ {K}_{{2}_{i}}$ and $ {K}_{{3}_{i}}$ denote the output feedback gain matrix $ {k}_{{s}_{i}}$ that satisfies the LMIs (33).

The WOA parameters were selected as:

• • Maximum number of iterations: $ \mathrm{M}\mathrm{a}{\mathrm{x}}_{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}=10$ ;

• • Number of search agents: $ \mathrm{S}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{h}\mathrm{A}\mathrm{g}\mathrm{e}\mathrm{n}\mathrm{t}{\mathrm{s}}_{\mathrm{n}\mathrm{o}}=700$ ;

• • Objective function: the ISE, which accounts for tracking errors across all controlled states ( $ \phi ,\theta ,\psi ,x,y,z$ ).

After determining the gain matrix $ {k}_{{s}_{i}}$ that satisfies the LMIs, the PID controller parameters can be obtained using Equation (35).

These parameter settings were found to be sufficient for determining the appropriate controller gains, as increasing the number of iterations or search agents significantly increases computation time without improving the obtained results.

To validate the adopted control scheme, extensive simulations were conducted before moving to real-world experiments. The validation scenario involves a sequence of configuration changes, simulating the navigation through a narrow gap with a height of 20 cm and a width of 30 cm.

The drone takes off in the ‘O’ configuration and then transitions to the ‘V’ configuration while moving along the Y-axis. As it approaches the narrow passage, it progressively reconfigured into the ‘L’ configuration, which offers a narrower profile compared to the other shapes. While passing through the narrow passage, the drone adopts a ‘U’ configuration, and subsequently aligns its rear arms to achieve the inverted ‘L’ configuration. This transition helps avoid potential singularities resulting from rotor alignment. After clearing the narrow passage, the drone returns to the ‘V’ configuration, and finally transitions back to the initial ‘O’ configuration. The variation of the desired and realised joint angles during the simulated scenario is shown in Fig. 9, where the blue curve represents the displacement phases and the yellow curve indicates the hovering phases allocated for the transition between configurations. The corresponding angles for these configurations are reported in Table 2.

Figure 9.

Evolution of joint angles for the RSQ while executing the mission of penetrating a narrow space.

Figure 10 illustrates the 3D trajectory, along with the corresponding configurations adopted throughout the flight. As observed, the RSQ successfully completed the assigned trajectory with minimal tracking error.

Figure 10.

Evolution of 3D trajectory.

The time response of the drone altitude $ Z$ is depicted in Fig. 11(c), while the trajectories along the $ X$ and $ Y$ axes are illustrated in Fig. 11(a) and (b), respectively.

Figure 11.

Evolution of position and error for the RSQ.

Figure 11 Long description

The image contains four line graphs depicting the tracking performance of position and error for a system. Panel A: A line graph shows the tracking performance of the x position over time. The x-axis represents time in seconds, and the y-axis represents the x position in meters. The graph includes two lines, one for the desired position (x_d) and one for the actual position (x). Panel B: A line graph shows the tracking performance of the y position over time. The x-axis represents time in seconds, and the y-axis represents the y position in meters. The graph includes two lines, one for the desired position (y_d) and one for the actual position (y). Panel C: A line graph shows the tracking performance of the z position over time. The x-axis represents time in seconds, and the y-axis represents the z position in meters. The graph includes two lines, one for the desired position (z_d) and one for the actual position (z). Panel D: A line graph shows the tracking performance of the error positions over time. The x-axis represents time in seconds, and the y-axis represents the error in meters. The graph includes three lines, one for the error in the x position (e_x), one for the error in the y position (e_y), and one for the error in the z position (e_z).

As shown in Fig. 11(d), the tracking error along the $ X$ -axis remains bounded within the range of $ -0.14$ m to $ 0.13$ m. For the $ Y$ -axis, the error reaches 0.15 m but rapidly decreases, often canceling out with each manoeuvre. Coupling between the $ X$ and $ Y$ motions is observed, particularly between [29.5–50] s and [59–69.5] s, where fluctuations under 2 cm appear on the $ Y$ -axis due to $ X$ -axis movements. The controller does not fully mitigate these coupling effects, as the narrow profile complicates control along the $ Y$ -axis. On the other hand, the movement along the $ Y$ -axis has less influence on the $ X$ -axis, due to the important rotor to CoG distance along the $ X$ -axis, as shown in Fig. 11. The $ Z$ -axis error starts at 0.6 meters at $ t=0s$ and decreases to zero as the desired altitude $ {Z}_{d}$ is reached (see Fig. 11(d)).

For attitude control, Fig. 12 illustrates the overall angular responses. The yaw angle, shown in Fig. 12(c), accurately follows the desired trajectory, with a tracking error below $ 0.05$ deg. Similarly, the roll ( $ \phi $ ) and pitch ( $ \theta $ ) angles, depicted in Fig. 12(a) and (b), respectively, also align well with the reference commands. As reported in Fig. 12(d), the maximum tracking error for roll and pitch remains below $ 1.1{5}^{\circ }$ .

Figure 12.

Evolution of Euler angles and error for the RSQ.

Figure 13 illustrates the time evolution of the control inputs during the simulated mission. As shown in Fig. 13(a), the thrust force $ {F}_{p}$ rapidly increases to approximately 25 N before stabilising at a constant value to counteract the gravitational force. The yaw control moment $ {M}_{\psi }$ ensures proper alignment with the desired trajectory, with a peak value of 0.02 N $ \cdot $ m (see Fig. 13(d)). The roll moment $ {M}_{\phi }$ reaches a maximum of 0.13 N $ \cdot $ m, while the pitch moment $ {M}_{\theta }$ peaks at 0.3 N $ \cdot $ m, as depicted in Fig. 13(b) and (c).

Figure 13.

Control signals input evolution.

The evolution of the CoG position relative to $ {O}_{\beta }$ is shown in Fig. 14, revealing significant variations throughout the flight scenario, ranging from $ -34$ to $ -7$ cm along the $ Y$ -axis, and from $ -7$ to $ 7$ cm along the $ X$ -axis.

Figure 14.

Evolution of the CoG position relative to the base frame $ {O}_{\beta }$ .

Figure 15 illustrates the time variation of the inertia matrix elements across different configurations, highlighting notable changes in $ {I}_{xx}$ , $ {I}_{yy}$ and $ {I}_{zz}$ , as well as a non-negligible $ {I}_{xy}$ component compared to $ {I}_{xz}$ and $ {I}_{yz}$ . These variations are significant and reflect the system’s dynamic structural reconfiguration.

Figure 15.

Evolution of the inertia parameters of the RSQ.

The variations of the derivatives of the inertia matrix elements are shown in Figs 16 and 17. The results indicate that, during reconfiguration, these derivatives reach up to $ 0.4\mathrm{K}\mathrm{g}\cdot {\mathrm{m}}^{2}/\mathrm{s}$ for $ {I}_{zz}$ , and $ 0.1\mathrm{K}\mathrm{g}\cdot {\mathrm{m}}^{2}/\mathrm{s}$ for $ {I}_{xx}$ and $ {I}_{yy}$ . For the inertia product $ {I}_{xy}$ , the derivative reaches $ 0.04\mathrm{K}\mathrm{g}\cdot {\mathrm{m}}^{2}/\mathrm{s}$ . It can be observed that these variations are not negligible in the system dynamics.

Figure 16.

Derivatives of the diagonal inertia elements.

Figure 17.

Derivatives of the off-diagonal inertia elements.

Evolution of the distances between the rotor positions $ {O}_{ij}$ and the CoG over time is presented in Fig. 18. These distances define the elements of the control allocation matrix. Unlike a conventional quadrotor, where such parameters remain constant, they vary substantially in the case of the RSQ due to its transformable architecture.

Figure 18.

Evolution of the distance between the CoG frame and rotor $ _{ij}$ frame.