Nonlinear Optimal Control of Magnetically Geared Induction Motors

G. Rigatos; P. Siano; M. Abbaszadeh; G. Cuccurullo; K. Ouahada

doi:10.2478/pead-2025-0009

Full Article

1.

Introduction

Magnetic gears exhibit several advantages over mechanical gears. Due to contactless operation, the constituent parts of magnetic gears are not subject to wear, which significantly reduces the need for their maintenance and increases their reliability (Habibi et al., 2024; Montegue et al., 2012). Magnetic gears are also characterised by inherent load protection under faults, torque transfer with reduced friction, low mechanical stresses and low acoustic noise (Song et al., 2022; Sun et al., 2017). In magnetic gears, energy losses are minimised and electric machines connected to them can operate efficiently across a wide range of speeds (Pop et al., 2018; Wang and Gerber, 2014). Magnetic gears are used in various applications, for instance, the traction system of electric vehicles (EVs), wind and tidal turbines and marine power generation units (Liao et al., 2023; McGilton et al., 2018). The coaxial magnetic gear configuration comprises coaxial rotating parts, and because of its very good motion transmission capability over a wide range of speeds and torques, it has become widely used in Hybrid Electric Vehicles (HEVs) and EVs (Tong et al., 2023; Xie et al., 2024). The dynamical system which is formed by connecting magnetic gears to the rotor of an electric machine (three-phase synchronous or induction machines or multi-phase synchronous and induction machines) exhibits complex non-linear dynamics (Long et al., 2023; Yang et al., 2024). The solution of the associated non-linear control problem is a non-trivial task, and so-far several non-linear control techniques have been proposed for it (e.g. sliding-mode control, non-linear model predictive control or global linearisation-based control schemes) (Druant et al., 2016; Xi et al., 2023). The application of non-linear control to magnetically geared electric machines comes also against several estimation issues due to harsh operating conditions which prevent the use of dedicated sensors for the state variables of these machines or for the mechanical load (Kumashira et al., 2004; Zhu, 2018). Magnetically-geared electric machines combine the merits of contactless motion transmission with the ability of the above-noted machines to vary their torque and speed over a wide range (Bidouche et al., 2020; Qu et al., 2011). Magnetically geared electric motors can be used efficiently in vehicles traction, while magnetically-geared electric power generators can be used for producing electric power from renewable energy sources. The related non-linear control problems are challenging and several noteworthy results have appeared in this area (Dobzhanskyi et al., 2019; Hazra et al., 2020).

The present article proposes a new non-linear optimal control method for the dynamic model of the magnetically-geared induction motor (MGIM) (Rigatos, 2016; Rigatos et al., 2024a). First, it is proven that the dynamic model of the induction motor with magnetic gears is differentially flat (Rigatos, 2015; Rigatos et al., 2022). Next, to apply this control scheme, the dynamic model of the MGIM undergoes approximate linearization around the temporary operating point (x^*, u^*) which is recomputed at each time step of the control algorithm, where x^* is the present value of the system’s state vector and u^* is the last sampled value of the control inputs vector (Basseville and Nikiforov, 1993; Rigatos and Tzafestas, 2007; Rigatos and Zhang, 2009). The linearization process is based on first-order Taylor series expansion and on the computation of Jacobian matrices (Rigatos et al., 2024b, 2025). The modelling error which is due to the truncation of higher-order terms from the Taylor series is considered to be a perturbation which is asymptotically compensated by the robustness of the control algorithm. For the approximately linearized model of the system an H-infinity feedback controller is designed (Rigatos et al., 2024b, 2025). To compute the stabilizing feedback gains of the H-infinity controller an algebraic Riccati equation has to be solved repetitively at each iteration of the control algorithm. The global stability properties of the control method are proven through Lyapunov analysis (Rigatos et al., 2024c; Toussaint et al., 2000).

2.

Dynamic Model of the Magnetically Geared Induction Motor

The diagram of a magnetically geared induction motor (MGIM) being used in the traction system of an EV is shown in Figure 1. The dynamic model of the induction motor is written in the dq asynchronously rotating frame using the assumption of field orientation (Rigatos, 2015; Rigatos et al., 2024a). This results in the following set of state equations: (1) ${\dot{ω}}_{m} = \frac{μ}{J} ψ_{r_{d}} i_{s_{q}} - \frac{T_{L}}{J}$ {\dot \omega _m} = {\mu \over J}{\psi _{{r_d}}}{i_{{s_q}}} - {{{T_L}} \over J} (2) ${\dot{ψ}}_{r_{d}} = - a ψ_{r_{d}} + {aMis}_{d}$ {\dot \psi _{{r_d}}} = - a{\psi _{{r_d}}} + aMi{s_d} (3) ${\dot{i}}_{s_{d}} = - γ i_{s_{d}} + a β ψ_{r_{q}} + n_{p} ω_{m} i_{s_{q}} + \frac{ω_{m} M i_{s_{q}}^{2}}{ψ_{r_{d}}} + \frac{1}{σ L_{s}} v_{s_{d}}$ {\dot i_{{s_d}}} = - \gamma {i_{{s_d}}} + a\beta {\psi _{{r_q}}} + {n_p}{\omega _m}{i_{{s_q}}} + {{{\omega _m}Mi_{{s_q}}^2} \over {{\psi _{{r_d}}}}} + {1 \over {\sigma {L_s}}}{v_{{s_d}}} (4) ${\dot{i}}_{s_{q}} = - γ i_{s_{q}} - β n_{p} ω_{m} ψ_{r_{d}} - n_{p} ω_{m} i_{s_{d}} - \frac{ω_{m} {Mi}_{s_{q}} i_{s_{q}}}{ψ_{r_{q}}} + \frac{1}{σ L_{s}} v_{s_{q}}$ {\dot i_{{s_q}}} = - \gamma {i_{{s_q}}} - \beta {n_p}{\omega _m}{\psi _{{r_d}}} - {n_p}{\omega _m}{i_{{s_d}}} - {{{\omega _m}M{i_{{s_q}}}{i_{{s_q}}}} \over {{\psi _{{r_q}}}}} + {1 \over {\sigma {L_s}}}{v_{{s_q}}} (5) $\dot{ρ} = n_{p} ω_{m} + \frac{{aMi}_{s_{q}}}{ψ_{r_{d}}}$ \dot \rho = {n_p}{\omega _m} + {{aM{i_{{s_q}}}} \over {{\psi _{{r_d}}}}} where, ω_m is the angular speed of the rotor, ψ_rd is the d-axis component of the magnetic flux of the rotor (due to field orientation, about the q-axis component of the magnetic flux it holds that ψ_rq = 0), i_sd is the d-axis component of the stator’s currents vector, i_sq is the q-axis component of the stator’s currents vector and ρ is the orientation angle of the magnetic field of the motor, which changes at an asynchronous rate. In the dynamic model of the induction motor parameters, α, β, γ, μ, L_s. M and σ are associated with the coefficients of the electric circuit of the stator and rotor (resistance, inductance, mutual inductance, magnetic permeability and number of poles).

The parameters of the dynamic model of the MGIM are outlined in Table 1:

Table 1.

Parameters of the MGIM dynamic model.

Parameter	Definition
ω_m, ω_L	Angular speed of the motor, load
J_m, J_L, J_g	Moment of inertia of the rotor, load, gear
T_e, T_L, T_g	Torque of the rotor, load, gear
ϕ	Angle denoting the speed difference between rotor and load
G_r	Transmission ratio of the magnetic gear
B_m, B_L, B_g	Friction coefficient at rotor, load and gear
p_o, p_m	Number of ferromagnetic pole pieces and air gaps
n_L	Sum of ferromagnetic pole pieces and air gaps
i_{s_d}, i_{s_q}	d,q axis components of the IM stator currents
R_s, R_r	Resistance of the IM’s stator, rotor
ψ_{r_d}, ψ_{r_q}	d,q axis components of the IM rotor flux
L_s, L_r	Inductance of the IM’s stator, rotor
M	Mutual inductance between IM’s stator and rotor
n_p	Number of poles of the IM’s stator
ρ	Orientation of the IM’s magnetic field
α, β, γ	$α = \frac{R_{r}}{L_{r}}, β = \frac{M}{σ L_{s} L_{s}}, γ = \frac{M^{2} R_{r}}{σ L_{s} {L_{r}}^{2}} + \frac{R_{s}}{σ L_{s}}$ \alpha = {{{R_r}} \over {{L_r}}},\beta = {M \over {\sigma {L_s}{L_s}}},\gamma = {{{M^2}{R_r}} \over {\sigma {L_s}{L_r}^2}} + {{{R_s}} \over {\sigma {L_s}}}
μ, σ	$μ = \frac{n_{p} M}{J}, σ = 1 - \frac{M^{2}}{L_{s} L_{r}}$ \mu = {{{n_p}M} \over J},\sigma = 1 - {{{M^2}} \over {{L_s}{L_r}}}

IM, Induction Motor; MGIM, magnetically-geared induction motor.

It holds that: $α = \frac{R_{r}}{L_{r}}$ \alpha = {{{R_r}} \over {{L_r}}} , $β = \frac{M}{σ L_{s} L_{r}}$ \beta = {M \over {\sigma {L_s}{L_r}}} , $γ = (\frac{M^{2} R_{r}}{σ L_{s} L_{r}^{2}} + \frac{R_{s}}{σ L_{s}})$ \gamma = \left( {{{{M^2}{R_r}} \over {\sigma {L_s}L_r^2}} + {{{R_s}} \over {\sigma {L_s}}}} \right) , $μ = \frac{n_{p} M}{J}$ \mu = {{{n_p}M} \over J} , where R_r is the rotor resistance, R_s is the stator resistance, L_r is the rotor inductance, L_s is the stator inductance, M i the mutual inductance, n_p is the number of poles and J is the moment of inertia of the rotor (Rigatos, 2015; Rigatos et al., 2024a).

The magnetic gear is composed of an inner, middle and an outer rotor. The inner rotor has a p_m number of pole pairs while the outer rotor has a p_o number of pole pairs. The middle rotor, called ‘cage rotor’, has n_L = p_o + p_m number of ferromagnetic pole pieces and air gaps. The inner rotor of the magnetic gear usually serves as the high-speed rotor connected to the prime mover. On the other hand, either the middle rotor is kept stationary and the outer rotor serves as the low-speed rotor or the middle rotor serves as the low-speed rotor and the outer rotor is kept stationary (Habibi et al., 2024). In the considered MGIM, the inner rotor is connected with the induction motor and is the high-speed part of the motion transmission system, while the middle rotor is connected with the traction system of the EV and is the low-speed part of the motion transmission system. The outer rotor is kept stationary. The electromagnetic torque of the induction motor is T_e and the IM’s turn speed is ω_m. The torque at the side of the traction system of the EV is T_L, and the turn speed of the load is ω_L. It holds the magnetic gear’s ratio is $G_{r} = \frac{ω_{m}}{ω_{L}} or G_{r} = \frac{n_{L}}{p_{m}}$ {G_r} = {{{\omega _m}} \over {{\omega _L}}}\;{\rm{or}}\;{G_r} = {{{n_L}} \over {{p_m}}} .

The dynamics of the magnetic gear is given by the following set of equations (Habibi et al., 2024): (6) ${\dot{ω}}_{m} = \frac{T_{e}}{J_{m}} - \frac{1}{J_{m} G_{r}} T_{g} \sin (ϕ) - \frac{B_{m} ω_{m}}{J_{m}}$ {\dot \omega _m} = {{{T_e}} \over {{J_m}}} - {1 \over {{J_m}{G_r}}}{T_g}sin\left( \phi \right) - {{{B_m}{\omega _m}} \over {{J_m}}} (7) ${\dot{ω}}_{L} = \frac{T_{g} \sin (ϕ)}{J_{L} + J_{g}} - \frac{1}{J_{L} + J_{g}} T_{L} - \frac{(B_{L} + B_{g}) ω_{L}}{J_{L} + J_{g}}$ {\dot \omega _L} = {{{T_g}sin\left( \phi \right)} \over {{J_L} + {J_g}}} - {1 \over {{J_L} + {J_g}}}{T_L} - {{\left( {{B_L} + {B_g}} \right){\omega _L}} \over {{J_L} + {J_g}}} (8) $\dot{φ} = p_{m} ω_{m} - n_{L} ω_{L}$ \dot \varphi = {p_m}{\omega _m} - {n_L}{\omega _L} with T_L = 0, which signifies that the load’s torque is constant or piece-wise constant. In Eq. (6), T_e = μψ_rdi_sq is the electromagnetic torque of the motor, while the torque of the magnetic gear at the rotor’s side is T_gsin(φ), which opposes to the motion of the inner rotor and the torque due to friction forces applied on the motor is B_mω_m. In Eq. (7), T_gsin(φ) is the torque of the magnetic gear, which causes the load’s turn motion, T_L is the unknown torque of the load which is taken to be piece-wise constant and (B_L + B_g)ω_L is the torque of friction forces applied to the gear and the load. In Eq. (8), φ is the ‘angle’ variable which determines the changes of the gear’s torque.

In Eqs (6)–(8), ω_m is the turn speed of the motor (high-speed rotation or HSR), ω_L is the turn speed of the load (low-speed rotation or LSR), φ is the ‘angle’ variable denoting the difference in the turn speed of the rotor from the turn speed of the load, G_r is the transmission ratio of the magnetic gear, J_m, J_g, J_L are the moments of inertia of the motor, gear and load, respectively, while B_m, B_g, B_L are the viscous damping coefficients of the motor, gear and load, respectively.

Thus, the aggregate dynamic model of the MGIM is: (9) ${\dot{ω}}_{m} = \frac{T_{e}}{J_{m}} - \frac{1}{J_{m} G_{r}} T_{g} \sin (ϕ) - \frac{B_{m} ω_{m}}{J_{m}}$ {\dot \omega _m} = {{{T_e}} \over {{J_m}}} - {1 \over {{J_m}{G_r}}}{T_g}sin\left( \phi \right) - {{{B_m}{\omega _m}} \over {{J_m}}} (10) ${\dot{ω}}_{L} = \frac{T_{g} \sin (ϕ)}{J_{L} + J_{g}} - \frac{1}{J_{L} + J_{g}} T_{L} - \frac{(B_{L} + B_{g}) ω_{L}}{J_{L} + J_{g}}$ {\dot \omega _L} = {{{T_g}sin\left( \phi \right)} \over {{J_L} + {J_g}}} - {1 \over {{J_L} + {J_g}}}{T_L} - {{\left( {{B_L} + {B_g}} \right){\omega _L}} \over {{J_L} + {J_g}}} (11) $\dot{ϕ} = p_{m} ω_{m} - n_{L} ω_{L}$ \dot \phi = {p_m}{\omega _m} - {n_L}{\omega _L} (12) ${\dot{ψ}}_{r_{d}} = - a ψ_{r_{d}} + {aMi}_{s_{d}}$ {\dot \psi _{{r_d}}} = - a{\psi _{{r_d}}} + aM{i_{{s_d}}} (13) ${\dot{i}}_{s_{d}} = - γ i_{s_{d}} + a β ψ_{r_{q}} + n_{p} ω_{m} i_{s_{q}} + \frac{ω_{m} M i_{s_{q}}^{2}}{ψ_{r_{d}}} + \frac{1}{σ L_{s}} v_{s_{d}}$ {\dot i_{{s_d}}} = - \gamma {i_{{s_d}}} + a\beta {\psi _{{r_q}}} + {n_p}{\omega _m}{i_{{s_q}}} + {{{\omega _m}Mi_{{s_q}}^2} \over {{\psi _{{r_d}}}}} + {1 \over {\sigma {L_s}}}{v_{{s_d}}} (14) ${\dot{i}}_{s_{q}} = - γ i_{s_{q}} - β n_{p} ω_{m} ψ_{r_{d}} - n_{p} ω_{m} i_{s_{d}} - \frac{ω_{m} {Mi}_{s_{q}} i_{s_{q}}}{ψ_{r_{d}}} + \frac{1}{σ L_{s}} v_{s_{q}}$ {\dot i_{{s_q}}} = - \gamma {i_{{s_q}}} - \beta {n_p}{\omega _m}{\psi _{{r_d}}} - {n_p}{\omega _m}{i_{{s_d}}} - {{{\omega _m}M{i_{{s_q}}}{i_{{s_q}}}} \over {{\psi _{{r_d}}}}} + {1 \over {\sigma {L_s}}}{v_{{s_q}}}

Eq. (5) about the orientation angle of the magnetic flux can be omitted from the above given dynamic model. Actually, the orientation of the magnetic field ρ is affected by state variables i_sq and ψ_rd, but has no impact on the rest of the state variables of the model.

Next, the state vector of the MGIM is defined as x = [x₁,x₂,x₃,x₄,x₅,x₆]^T ⇒ or x = [ω_m,ω_L,φ,ψ_rd,i_sd,i_sq]^T and the control inputs vector of the MGIM is defined as: u = [u₁, u₂]^T ⇒ u = [υ_{s_d}, υ_{s_q}]^T. This results in the state-space description: (15) $(\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ {\dot{x}}_{3} \\ {\dot{x}}_{4} \\ {\dot{x}}_{5} \\ {\dot{x}}_{6} \end{matrix}) = (\begin{matrix} \frac{μ x_{4} x_{6}}{J_{m}} - \frac{1}{J_{m} G_{r}} T_{g} \sin (x_{3}) - \frac{B_{m} x_{1}}{J_{m}} \\ \frac{T_{g} \sin (x_{3})}{J_{L} + J_{g}} - \frac{1}{J_{L} + J_{g}} T_{L} - \frac{(B_{L} + B_{g}) x_{2}}{J_{L} + J_{g}} \\ p_{m} x_{1} - n_{L} x_{2} \\ - {ax}_{4} + {aMx}_{5} \\ - γ x_{5} + a β x_{4} + n_{p} x_{1} x_{6} + \frac{x_{1} M x_{6}^{2}}{x_{4}} \\ - γ x_{6} - β n_{p} x_{1} x_{4} - n_{p} x_{1} x_{5} - \frac{x_{1} {Mx}_{5} x_{6}}{x_{4}} \end{matrix}) + (\begin{matrix} 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \\ \frac{1}{σ L_{s}} & 0 \\ 0 & \frac{1}{σ L_{s}} \end{matrix}) (\begin{matrix} u_{1} \\ u_{2} \end{matrix})$ \left( {\matrix{ {{{\dot x}_1}} \cr {{{\dot x}_2}} \cr {{{\dot x}_3}} \cr {{{\dot x}_4}} \cr {{{\dot x}_5}} \cr {{{\dot x}_6}} \cr } } \right) = \left( {\matrix{ {{{\mu {x_4}{x_6}} \over {{J_m}}} - {1 \over {{J_m}{G_r}}}{T_g}sin({x_3}) - {{{B_m}{x_1}} \over {{J_m}}}} \cr {{{{T_g}sin\left( {{x_3}} \right)} \over {{J_L} + {J_g}}} - {1 \over {{J_L} + {J_g}}}{T_L} - {{\left( {{B_L} + {B_g}} \right){x_2}} \over {{J_L} + {J_g}}}} \cr {{p_m}{x_1} - {n_L}{x_2}} \cr { - a{x_4} + aM{x_5}} \cr { - \gamma {x_5} + a\beta {x_4} + {n_p}{x_1}{x_6} + {{{x_1}Mx_6^2} \over {{x_4}}}} \cr { - \gamma {x_6} - \beta {n_p}{x_1}{x_4} - {n_p}{x_1}{x_5} - {{{x_1}M{x_5}{x_6}} \over {{x_4}}}} \cr } } \right) + \left( {\matrix{ 0 & 0 \cr 0 & 0 \cr 0 & 0 \cr 0 & 0 \cr {{1 \over {\sigma {L_s}}}} & 0 \cr 0 & {{1 \over {\sigma {L_s}}}} \cr } } \right)\left( {\matrix{ {{u_1}} \cr {{u_2}} \cr } } \right)

Using the x∈R^6×1, f(x)∈R^6×1, g(x)∈R^6×2 and u∈R^2×1, the dynamic model of the MGIM is finally written in the following concise non-linear affine-in-the-input state-space form (16) $\dot{x} = f (x) + g (x) u$ \dot x = f\left( x \right) + g\left( x \right)u

3.

Differential Flatness of the MGIM

It can be proven that the dynamic model of the MGIM, which was previously described in the state-space model of Eq. (15), is differentially flat, with flat outputs vector Y = [x₂,x₄]^T = [ω_L, ψ_rd]^T. A system is differentially flat if (i) all its state variables and its control inputs can be written as differential relations of a subset of the state vector elements, which constitute the flat outputs vector of the system, (ii) the flat outputs vector and its time-derivatives are differentially independent which means that they are not connected through relations in the form of an homogenous differential equation. From the second row of the state-space model and considering that T_L is a piece-wise constant, one solves for x₃. This gives: (17) $\begin{matrix} \sin (x_{3}) = \frac{J_{L} + J_{g}}{J_{g}} {\dot{x}}_{2} + \frac{B_{g} + B_{L}}{J_{g} + J_{L}} x_{2} + \frac{T_{L}}{J_{g} + J_{L}}] \Rightarrow \\ x_{3} = \sin^{- 1} (\frac{J_{L} + J_{g}}{J_{g}} {\dot{x}}_{2} + \frac{B_{g} + B_{L}}{J_{g} + J_{L}} x_{2} + \frac{T_{L}}{J_{g} + J_{L}}]\}) \Rightarrow \\ x_{3} = h_{3} (Y, \dot{Y}) \end{matrix}$ \matrix{ {sin\left( {{x_3}} \right) = {{{J_L} + {J_g}} \over {{J_g}}}\left[ {{{\dot x}_2} + {{{B_g} + {B_L}} \over {{J_g} + {J_L}}}{x_2} + {{{T_L}} \over {{J_g} + {J_L}}}} \right] \Rightarrow } \cr {{x_3} = si{n^{ - 1}}\left( {\left\{ {{{{J_L} + {J_g}} \over {{J_g}}}\left[ {{{\dot x}_2} + {{{B_g} + {B_L}} \over {{J_g} + {J_L}}}{x_2} + {{{T_L}} \over {{J_g} + {J_L}}}} \right]} \right\}} \right) \Rightarrow } \cr {{x_3} = {h_3}\left( {Y,\dot Y} \right)} \cr } which signifies that x₃ is a differential function of the flat outputs vector. From the third row of the state-space model, one solves for the state variable, x₁. This gives (18) $\begin{matrix} x_{1} = p_{m} x_{1} - n_{L} x_{2} \Rightarrow x_{1} = \frac{1}{p_{m}} (x_{3} + {\dot{n}}_{L} x_{2}) \\ \Rightarrow x_{1} = h_{1} (Y, \dot{Y}) \end{matrix}$ \matrix{ {{x_1} = {p_m}{x_1} - {n_L}{x_2} \Rightarrow {x_1} = {1 \over {{p_m}}}\left( {{x_3} + {{\dot n}_L}{x_2}} \right)} \cr { \Rightarrow {x_1} = {h_1}\left( {Y,\;\dot Y} \right)} \cr } which signifies that x₁ is a differential function of the flat outputs vector. From the first row of the state-space model, one solves for the state variable, x₆. This gives (19) $\begin{matrix} x_{6} = \frac{J_{m}}{μ x_{4}} {\dot{x}}_{1} + \frac{1}{J_{m} G_{r}} T_{g} \sin (x_{3}) + \frac{B_{m}}{J_{m}} x_{1}] \Rightarrow \\ x_{6} = h_{7} (Y, \dot{Y}) \end{matrix}$ \matrix{ {{x_6} = {{{J_m}} \over {\mu {x_4}}}\left[ {{{\dot x}_1} + {1 \over {{J_m}{G_r}}}{T_g}sin\left( {{x_3}} \right) + {{{B_m}} \over {{J_m}}}{x_1}} \right] \Rightarrow } \cr {{x_6} = {h_7}\left( {Y,\dot Y} \right)} \cr } which signifies that x₆ is a differential function of the flat output vector. From the fifth row of the state-space model, one solves for the state variable, x₅. This gives (20) $x_{5} = \frac{1}{aM} {\dot{x}}_{4} + {ax}_{4}] \Rightarrow x_{5} = h_{6} (Y, \dot{Y})$ {x_5} = {1 \over {aM}}\left[ {{{\dot x}_4} + a{x_4}} \right] \Rightarrow {x_5} = {h_6}\left( {Y,\dot Y} \right) which signifies that x₅ is a differential function of the flat outputs vector. Next, from the sixth row of the state-space model, one solves for the control input, u₁. This gives (21) $\begin{matrix} u_{1} = σ L_{s} {\dot{x}}_{5} + γ x_{5} - a β x_{4} - n_{p} x_{1} x_{6} - \frac{x_{1} M x_{6}^{2}}{x_{4}}] \\ \Rightarrow u_{1} = h_{u_{1}} (Y, \dot{Y}) \end{matrix}$ \matrix{ {{u_1} = \sigma {L_s}\left[ {{{\dot x}_5} + \gamma {x_5} - a\beta {x_4} - {n_p}{x_1}{x_6} - {{{x_1}Mx_6^2} \over {{x_4}}}} \right]} \cr { \Rightarrow {u_1} = {h_{{u_1}}}\left( {Y,\dot Y} \right)} \cr } which signifies that control input u₁ is a differential function of the flat outputs vector, Y. Finally, from the seventh row of the state-space model, one solves for the control input, u₂. This gives (22) $\begin{matrix} u_{2} = σ L_{s} {\dot{x}}_{6} + γ x_{7} + β n_{p} x_{1} x_{4} + n_{p} x_{1} x_{5} + \frac{x_{1} {Mx}_{5} x_{6}}{x_{4}}] \\ \Rightarrow u_{2} = h_{u_{2}} (Y, \dot{Y}) \end{matrix}$ \matrix{ {{u_2} = \sigma {L_s}\left[ {{{\dot x}_6} + \gamma {x_7} + \beta {n_p}{x_1}{x_4} + {n_p}{x_1}{x_5} + {{{x_1}M{x_5}{x_6}} \over {{x_4}}}\;} \right]} \cr { \Rightarrow {u_2} = {h_{{u_2}}}\left( {Y,\dot Y} \right)} \cr } which signifies that control input u₂ is a differential function of the flat outputs vector, Y. As a result of the above, all state variables and the control inputs of the dynamic model of the MGIM are differential functions of the flat outputs vector, Y. Consequently, this system is differentially flat.

4.

Approximate Linearisation of the MGIM

The dynamic model of the MGIM, being initially in the non-linear form $\dot{x} = f (x) + g (x) u$ \dot x = f\left( x \right) + g\left( x \right)u undergoes approximate linearisation around the temporary operating point (x^*,u^*), where x^* is the present value of the system’s state vector and u^* is the last sampled value of the control inputs vector. The linearisation is based on the use of first-order Taylor series expansion and on the Jacobian matrices of the system and takes place at each sampling instance. This gives: (23) $\dot{x} = Ax + Bu + \tilde{d}$ \dot x = Ax + Bu + \tilde d where A, B are the Jacobian matrices of the system and $\tilde{d}$ \tilde d is the cumulative disturbance term, which may comprise (i) the modelling error due to the truncation of higher order terms from the Taylor series, (ii) exogenous perturbations and (iii) sensor measurement noise of any distribution. Using that the control inputs gain matrix g(x) is time-invariant, the Jacobian matrices A, B are computed as follows: (24) ${A = \nabla_{x} f (x) + g (x) u]|}_{(x^{*}, u^{*})} \Rightarrow {A = \nabla_{x} f (x)|}_{(x^{*}, u^{*})}$ {\left. {A = {\nabla _x}\left[ {f\left( x \right) + g\left( x \right)u} \right]} \right|_{\left( {{x^*},{u^*}} \right)}} \Rightarrow {\left. {A = {\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}} (25) ${B = \nabla_{u} f (x) + g (x) u]|}_{(x^{*}, u^{*})} \Rightarrow {B = g (x)|}_{(x^{*}, u^{*})}$ {\left. {B = {\nabla _u}\left[ {f\left( x \right) + g\left( x \right)u} \right]} \right|_{\left( {{x^*},{u^*}} \right)}} \Rightarrow {\left. {B = g\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}

Computation of the Jacobian matrix ∇_xf(x) |(x^*,u^*):

First row of the Jacobian matrix ${\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{1}}{\partial x_{1}} = - \frac{B_{m}}{J_{m}}, \frac{\partial f_{1}}{\partial x_{2}} = 0, \frac{\partial f_{1}}{\partial x_{3}} = - \frac{T_{g}}{J_{m} G_{r}} \cos (x_{3}), \frac{\partial f_{1}}{\partial x_{4}} = \frac{μ x_{6}}{J_{m}}, \frac{\partial f_{1}}{\partial x_{5}} = 0, \frac{\partial f_{1}}{\partial x_{6}} = \frac{μ x_{4}}{J_{m}} .$ {\left. {{\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_1}} \over {\partial {x_1}}} = - {{{B_m}} \over {{J_m}}},{{\partial {f_1}} \over {\partial {x_2}}} = 0,{{\partial {f_1}} \over {\partial {x_3}}} = - {{{T_g}} \over {{J_m}{G_r}}}cos\left( {{x_3}} \right),{{\partial {f_1}} \over {\partial {x_4}}} = {{\mu {x_6}} \over {{J_m}}},{{\partial {f_1}} \over {\partial {x_5}}} = 0,{{\partial {f_1}} \over {\partial {x_6}}} = {{\mu {x_4}} \over {{J_m}}}.
Second row of the Jacobian matrix ${\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{2}}{\partial x_{1}} = 0, \frac{\partial f_{2}}{\partial x_{2}} = - \frac{B_{g} + B_{L}}{J_{g} + J_{L}}, \frac{\partial f_{2}}{\partial x_{3}} = \frac{T_{g}}{J_{L} + J_{g}} \cos (x_{3}), \frac{\partial f_{2}}{\partial x_{4}} = 0, \frac{\partial f_{2}}{\partial x_{5}} = 0, \frac{\partial f_{2}}{\partial x_{6}} = 0 .$ {\left. {{\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_2}} \over {\partial {x_1}}} = 0,{{\partial {f_2}} \over {\partial {x_2}}} = - {{{B_g} + {B_L}} \over {{J_g} + {J_L}}},{{\partial {f_2}} \over {\partial {x_3}}} = {{{T_g}} \over {{J_L} + {J_g}}}cos\left( {{x_3}} \right),{{\partial {f_2}} \over {\partial {x_4}}} = 0,{{\partial {f_2}} \over {\partial {x_5}}} = 0,{{\partial {f_2}} \over {\partial {x_6}}} = 0.
Third row of the Jacobian matrix ${\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{3}}{\partial x_{1}} = p_{m}, \frac{\partial f_{3}}{\partial x_{2}} = - n_{L}, \frac{\partial f_{3}}{\partial x_{3}} = 0, \frac{\partial f_{3}}{\partial x_{4}} = 0, \frac{\partial f_{3}}{\partial x_{5}} = 0, \frac{\partial f_{3}}{\partial x_{6}} = 0 .$ {\left. {{\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_3}} \over {\partial {x_1}}} = {p_m},{{\partial {f_3}} \over {\partial {x_2}}} = - {n_L},{{\partial {f_3}} \over {\partial {x_3}}} = 0,{{\partial {f_3}} \over {\partial {x_4}}} = 0,{{\partial {f_3}} \over {\partial {x_5}}} = 0,{{\partial {f_3}} \over {\partial {x_6}}} = 0.
Fourth row of the Jacobian matrix ${\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{4}}{\partial x_{1}} = 0, \frac{\partial f_{4}}{\partial x_{2}} = 0, \frac{\partial f_{4}}{\partial x_{3}} = 0, \frac{\partial f_{4}}{\partial x_{4}} = - a, \frac{\partial f_{4}}{\partial x_{5}} = aM, \frac{\partial f_{4}}{\partial x_{6}} = 0 .$ {\left. {{\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_4}} \over {\partial {x_1}}} = 0,{{\partial {f_4}} \over {\partial {x_2}}} = 0,{{\partial {f_4}} \over {\partial {x_3}}} = 0,{{\partial {f_4}} \over {\partial {x_4}}} = - a,{{\partial {f_4}} \over {\partial {x_5}}} = aM,{{\partial {f_4}} \over {\partial {x_6}}} = 0.
Fifth row of the Jacobian matrix ${\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{5}}{\partial x_{1}} = n_{p} x_{6} + \frac{M x_{6}^{2}}{x_{4}}, \frac{\partial f_{5}}{\partial x_{2}} = 0, \frac{\partial f_{5}}{\partial x_{3}} = 0, \frac{\partial f_{5}}{\partial x_{4}} = a β - \frac{x_{1} M x_{6}^{2}}{x_{4}^{2}}, \frac{\partial f_{5}}{\partial x_{5}} = - γ, \frac{\partial f_{5}}{\partial x_{6}} = n_{p} x_{1} + \frac{2 x_{1} {Mx}_{6}}{x_{4}} .$ {\left. {{\nabla _x}f\left( x \right)} \right|_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_5}} \over {\partial {x_1}}} = {n_p}{x_6} + {{Mx_6^2} \over {{x_4}}},{{\partial {f_5}} \over {\partial {x_2}}} = 0,{{\partial {f_5}} \over {\partial {x_3}}} = 0,{{\partial {f_5}} \over {\partial {x_4}}} = a\beta - {{{x_1}Mx_6^2} \over {x_4^2}},{{\partial {f_5}} \over {\partial {x_5}}} = - \gamma ,{{\partial {f_5}} \over {\partial {x_6}}} = {n_p}{x_1} + {{2{x_1}M{x_6}} \over {{x_4}}}.
Sixth row of the Jacobian matrix $\begin{array}{l} {\nabla_{x} f (x)|}_{(x^{*}, u^{*})} : \frac{\partial f_{6}}{\partial x_{1}} = - β n_{p} x_{4} - n_{p} x_{5} - \frac{{Mx}_{5} x_{6}}{x_{4}}, \frac{\partial f_{6}}{\partial x_{2}} = 0, \frac{\partial f_{6}}{\partial x_{3}} = 0, \\ \frac{\partial f_{6}}{\partial x_{4}} = - β n_{p} x_{1} + \frac{{Mx}_{1} x_{5} x_{6}}{x_{4}^{2}}, \frac{\partial f_{6}}{\partial x_{5}} = - n_{p} x_{1} - \frac{{Mx}_{1} x_{6}}{x_{4}}, \frac{\partial f_{6}}{\partial x_{6}} = - γ - \frac{{Mx}_{1} x_{5}}{x_{4}} . \end{array}$ \matrix{ {{{\left. {{\nabla _x}f\left( x \right)} \right|}_{\left( {{x^*},{u^*}} \right)}}:{{\partial {f_6}} \over {\partial {x_1}}} = - \beta {n_p}{x_4} - {n_p}{x_5} - {{M{x_5}{x_6}} \over {{x_4}}},{{\partial {f_6}} \over {\partial {x_2}}} = 0,{{\partial {f_6}} \over {\partial {x_3}}} = 0,} \hfill \cr {{{\partial {f_6}} \over {\partial {x_4}}} = - \beta {n_p}{x_1} + {{M{x_1}{x_5}{x_6}} \over {x_4^2}},{{\partial {f_6}} \over {\partial {x_5}}} = - {n_p}{x_1} - {{M{x_1}{x_6}} \over {{x_4}}},{{\partial {f_6}} \over {\partial {x_6}}} = - \gamma - {{M{x_1}{x_5}} \over {{x_4}}}.} \hfill \cr }

5.

Design of an H-Infinity Feedback Controller

5.1.

Equivalent linearised dynamics of the MGIM

After linearisation around its current operating point, the dynamic model for the MGIM is written as: (26) $\dot{x} = Ax + Bu + d_{1}$ \dot x = Ax + Bu + {d_1}

Parameter d₁ stands for the linearisation error in the MGIM’s model that was given previously in Eq. (26). The reference setpoints for the state vector of the aforementioned dynamic model are denoted by $x_{d} = x_{1}^{d}, \dots, x_{6}^{d}]$ {x_d} = \left[ {x_1^d, \cdots ,x_6^d} \right] Tracking of this trajectory is achieved after applying the control input u^*. At every time-instant, the control input u^* is assumed to differ from the control input u appearing in Eq. (26) by an amount equal to Δu, i.e. u^* = u + Δu. (27) ${\dot{x}}_{d} = {Ax}_{d} + {Bu}^{*} + d_{2}$ {\dot x_d} = A{x_d} + B{u^*} + {d_2}

The dynamics of the controlled system described in Eq. (26) can also be written as (28) $\dot{x} = Ax + Bu + {Bu}^{*} - {Bu}^{*} + d_{1}$ \dot x = Ax + Bu + B{u^*} - B{u^*} + {d_1} and by denoting d₃ = −Bu^* + d₁ as an aggregate disturbance term one obtains (29) $\dot{x} = Ax + Bu + {Bu}^{*} + d_{3}$ \dot x = Ax + Bu + B{u^*} + {d_3}

By subtracting Eq. (27) from Eq. (29), one has (30) $\dot{x} - {\dot{x}}_{d} = A (x - x_{d}) + Bu + d_{3} - d_{2}$ \dot x - {\dot x_d} = A\left( {x - {x_d}} \right) + Bu + {d_3} - {d_2}

By denoting the tracking error as e = x − x_d and the aggregate disturbance term as $L \tilde{d} = d_{3} - d_{2}$ L\tilde d = {d_3} - {d_2} , the tracking error dynamics becomes (31) $\dot{e} = Ae + Bu + L \tilde{d}$ \dot e = Ae + Bu + L\tilde d where L is the disturbance inputs gain matrix. The above linearised form of the MGIM can be efficiently controlled after applying an H-infinity feedback control scheme.

5.2.

The non-linear H-infinity control

The initial non-linear model of the MGIM is in the form (32) $\dot{x} = f (x, u) x \in R^{n}, u \in R^{m}$ \dot x = f\left( {x,\;u} \right)\;\;\;x \in {R^n},\;\;\;u \in {R^m}

Linearisation of the MGIM is performed at each iteration of the control algorithm around its present operating point (x^*,u^*) = (x(t),u(t − T_s)). The linearised equivalent of the system is described by (33) $\begin{matrix} \dot{x} = Ax + Bu + L \tilde{d} x \in R^{n}, u \in R^{m}, \tilde{d} \in R^{q} \\ y = Cx \end{matrix}$ \matrix{ {\dot x = Ax + Bu + L\tilde d\;\;\;x \in {R^n},\;\;\;u \in {R^m},\;\;\;\tilde d \in {R^q}} \cr {y = Cx} \cr } where matrices A and B are obtained from the computation of the previously defined Jacobians and vector $\tilde{d}$ \tilde d denotes disturbance terms due to linearisation errors, while L is a disturbance input gain matrix.

The problem of disturbance rejection for the linearised model cannot be handled efficiently if the classical Linear Quadratic Regulator (LQR) control scheme is applied. This is because of the existence of the perturbation term $\tilde{d}$ \tilde d . The disturbance term $\tilde{d}$ \tilde d , apart from modelling (parametric) uncertainty and external perturbation terms, can also represent noise terms of any distribution. On the contrary, in the H∞ control approach, a feedback control scheme is designed for trajectory tracking by the system’s state vector and simultaneous disturbance rejection, considering that the disturbance affects the system in the worst possible manner. The disturbances’ effects are incorporated in the following quadratic cost function: (34) $J (t) = \frac{1}{2} \int_{0}^{T} y^{T} (t) y (t) + {ru}^{T} (t) u (t) - ρ^{2} {\tilde{d}}^{T} (t) \tilde{d} (t)] dt, r, ρ > 0$ J\left( t \right) = {1 \over 2}\int_0^T {\left[ {{y^T}\left( t \right)y\left( t \right) + r{u^T}\left( t \right)u\left( t \right) - {\rho ^2}{{\tilde d}^T}\left( t \right)\tilde d\left( t \right)} \right]dt} ,\;\;\;r,\rho > 0

5.3.

Computation of the feedback control gains

For the linearised system given by Eq. (33), the cost function of Eq. (34) is defined, where coefficient r determines the penalisation of the control input and the weight coefficient ρ determines the reward of the disturbances’ effects. It is assumed that (i) the energy that is transferred from the disturbance signal $\tilde{d} (t)$ \tilde d\left( t \right) is bounded, i.e. $\int_{0}^{\infty} {\tilde{d}}^{T} (t) \tilde{d} (t) dt < \infty$ \int_0^\infty {{{\tilde d}^T}\left( t \right)\tilde d} \left( t \right)dt < \infty , (ii) matrices [A,B] and [A,L] are stabilisable, (iii) matrix [A,C] is detectable. In the case of a tracking problem, the optimal feedback control law is given by (35) $u (t) = - Ke (t)$ u\left( t \right) = - Ke\left( t \right) with e = x−x_d to be the tracking error, and $K = \frac{1}{r} B^{T} P$ K = {1 \over r}{B^T}P , where P is a positive definite symmetric matrix. As it will be proven in Section 6, matrix P is obtained from the solution of the Riccati equation (36) $A^{T} P + PA + Q - P (\frac{2}{r} {BB}^{T} - \frac{1}{ρ^{2}} {LL}^{T}) P = 0$ {A^T}P + PA + Q - P\left( {{2 \over r}B{B^T} - {1 \over {{\rho ^2}}}L{L^T}} \right)P = 0 where Q is a positive semi-definite symmetric matrix. The worst-case disturbance is given by (37) $\tilde{d} (t) = \frac{1}{ρ^{2}} L^{T} Pe (t) .$ \tilde d\left( t \right) = {1 \over {{\rho ^2}}}{L^T}Pe\left( t \right).

The solution of the H-infinity feedback control problem for the MGIM and the computation of the worst-case disturbance that the related controller can sustain, comes from superposition of Bellman’s optimality principle when considering that the induction motor with magnetic gears is affected by two separate inputs (i) the control input u and (ii) the cumulative disturbance input $\tilde{d} (t)$ \tilde d\left( t \right) . Solving the optimal control problem for u, (Gorecki 2018; Grimble and Majeski 2020; Molley et al., 2022; Tao et al., 2021) i.e. for the minimum variation (optimal) control input that achieves elimination of the state vector’s tracking error, gives $u = - \frac{1}{r} B^{T} Pe$ u = - {1 \over r}{B^T}Pe . Equivalently, solving the optimal control problem for $\tilde{d}$ \tilde d , i.e. for the worst-case disturbance that the control loop can sustain gives $\tilde{d} = \frac{1}{ρ^{2}} L^{T} Pe$ \tilde d = {1 \over {{\rho ^2}}}{L^T}Pe . The diagram of the considered control loop for the MGIM is given in Figure 2.

6.

Lyapunov Stability Analysis

Through Lyapunov stability analysis, it will be shown that the proposed non-linear control scheme assures H∞ tracking performance for the MGIM, and that, in case of bounded disturbance terms, asymptotic convergence to the reference setpoints is achieved. The tracking error dynamics for the MGIM are written in the form of Eq. (31) $\dot{e} = Ae + Bu + L \tilde{d}$ \dot e = Ae + Bu + L\tilde d , where in the MGIM’s case, L = ∈R^6×6 is the disturbance inputs gain matrix. Variable $\tilde{d}$ \tilde d denotes model uncertainties and external disturbances of the model of the MGIM. The following Lyapunov equation is considered (38) $V = \frac{1}{2} e^{T} Pe$ V = {1 \over 2}{e^T}Pe where e = x − x_d is the tracking error. By differentiating with respect to time, one obtains (39) $\begin{matrix} \dot{V} = \frac{1}{2} {\dot{e}}^{T} Pe + \frac{1}{2} e^{T} P \dot{e} \Rightarrow \dot{V} = \frac{1}{2} {Ae + Bu + L \tilde{d}]}^{T} Pe + \frac{1}{2} e^{T} P Ae + Bu + L \tilde{d}] \Rightarrow \\ \dot{V} = \frac{1}{2} e^{T} A^{T} + u^{T} B^{T} + {\tilde{d}}^{T} L^{T}] Pe + \frac{1}{2} e^{T} P Ae + Bu + L \tilde{d}] \Rightarrow \\ \dot{V} = \frac{1}{2} e^{T} A^{T} Pe + \frac{1}{2} u^{T} B^{T} Pe + \frac{1}{2} {\tilde{d}}^{T} L^{T} Pe + \frac{1}{2} e^{T} PAe + \frac{1}{2} e^{T} PBu + \frac{1}{2} e^{T} PL \tilde{d} \end{matrix}$ \matrix{ {\dot V = {1 \over 2}{{\dot e}^T}Pe + {1 \over 2}{e^T}P\dot e \Rightarrow \dot V = {1 \over 2}{{\left[ {Ae + Bu + L\tilde d} \right]}^T}Pe + {1 \over 2}{e^T}P\left[ {Ae + Bu + L\tilde d} \right] \Rightarrow } \cr {\dot V = {1 \over 2}\left[ {{e^T}{A^T} + {u^T}{B^T} + {{\tilde d}^T}{L^T}} \right]Pe + {1 \over 2}{e^T}P\left[ {Ae + Bu + L\tilde d} \right] \Rightarrow } \cr {\dot V = {1 \over 2}{e^T}{A^T}Pe + {1 \over 2}{u^T}{B^T}Pe + {1 \over 2}{{\tilde d}^T}{L^T}Pe + {1 \over 2}{e^T}PAe + {1 \over 2}{e^T}PBu + {1 \over 2}{e^T}PL\tilde d} \cr }

The previous equation is rewritten as (40) $\dot{V} = \frac{1}{2} e^{T} (A^{T} P + PA) e + (\frac{1}{2} u^{T} B^{T} Pe + \frac{1}{2} e^{T} PBu) + (\frac{1}{2} {\tilde{d}}^{T} L^{T} Pe + \frac{1}{2} e^{T} PL \tilde{d})$ \dot V = {1 \over 2}{e^T}\left( {{A^T}P + PA} \right)e + \left( {{1 \over 2}{u^T}{B^T}Pe + {1 \over 2}{e^T}PBu} \right) + \left( {{1 \over 2}{{\tilde d}^T}{L^T}Pe + {1 \over 2}{e^T}PL\tilde d} \right)

Assumption

For a given positive definite matrix Q and coefficients r and ρ, there exists a positive definite matrix P, which is the solution of the following matrix equation (41) $A^{T} P + PA = - Q + P (\frac{2}{r} {BB}^{T} - \frac{1}{ρ^{2}} {LL}^{T}) P$ {A^T}P + PA = - Q + P\left( {{2 \over r}B{B^T} - {1 \over {{\rho ^2}}}L{L^T}} \right)P

Moreover, the following feedback control law is applied to the system (42) $u = - \frac{1}{r} B^{T} Pe$ u = - {1 \over r}{B^T}Pe

By substituting Eqs (41) and (42), one obtains (43) $\begin{matrix} \dot{V} = \frac{1}{2} e^{T} - Q + P (\frac{2}{r} {BB}^{T} - \frac{1}{ρ^{2}} {LL}^{T}) P] e + e^{T} PB (- \frac{1}{r} B^{T} Pe) + e^{T} PL \tilde{d} \Rightarrow \\ \dot{V} = - \frac{1}{2} e^{T} Qe + \frac{1}{r} e^{T} {PBB}^{T} Pe - \frac{1}{2 ρ^{2}} e^{T} {PLL}^{T} Pe - \frac{1}{r} + e^{T} {PBB}^{T} Pe + e^{T} PL \tilde{d} \end{matrix}$ \matrix{ {\dot V = {1 \over 2}{e^T}\left[ { - Q + P\left( {{2 \over r}B{B^T} - {1 \over {{\rho ^2}}}L{L^T}} \right)P} \right]e + {e^T}PB\left( { - {1 \over r}{B^T}Pe} \right) + {e^T}PL\tilde d \Rightarrow } \cr {\dot V = - {1 \over 2}{e^T}Qe + {1 \over r}{e^T}PB{B^T}Pe - {1 \over {2{\rho ^2}}}{e^T}PL{L^T}Pe - {1 \over r} + {e^T}PB{B^T}Pe + {e^T}PL\tilde d} \cr } which, after intermediate operations, gives (44) $\begin{matrix} \dot{V} = - \frac{1}{2} e^{T} Qe - \frac{1}{2 ρ^{2}} e^{T} {PLL}^{T} Pe + e^{T} PL \tilde{d} \Rightarrow \\ \dot{V} = - \frac{1}{2} e^{T} Qe - \frac{1}{2 ρ^{2}} e^{T} {PLL}^{T} Pe + \frac{1}{2} e^{T} PL \tilde{d} + \frac{1}{2} {\tilde{d}}^{T} L^{T} Pe \end{matrix}$ \matrix{ {\dot V = - {1 \over 2}{e^T}Qe - {1 \over {2{\rho ^2}}}{e^T}PL{L^T}Pe + {e^T}PL\tilde d \Rightarrow } \cr {\dot V = - {1 \over 2}{e^T}Qe - {1 \over {2{\rho ^2}}}{e^T}PL{L^T}Pe + {1 \over 2}{e^T}PL\tilde d + {1 \over 2}{\tilde d^T}{L^T}Pe} \cr }

Lemma

The following inequality holds (45) $\frac{1}{2} e^{T} L \tilde{d} + \frac{1}{2} \tilde{d} L^{T} Pe - \frac{1}{2 ρ^{2}} e^{T} {PLL}^{T} Pe \leq \frac{1}{2} ρ^{2} {\tilde{d}}^{T} \tilde{d}$ {1 \over 2}{e^T}L\tilde d + {1 \over 2}\tilde d{L^T}Pe - {1 \over {2{\rho ^2}}}{e^T}PL{L^T}Pe \le {1 \over 2}{\rho ^2}{\tilde d^T}\tilde d

Proof

The binomial ${(ρ α - \frac{1}{ρ} b)}^{2}$ {\left( {\rho \alpha - {1 \over \rho }b} \right)^2} is considered. Expanding the left part of the above inequality, one gets (46) $\begin{matrix} ρ^{2} a^{2} + \frac{1}{ρ^{2}} b^{2} - 2 ab \geq 0 \Rightarrow \frac{1}{2} ρ^{2} a^{2} + \frac{1}{2 ρ^{2}} b^{2} - ab \geq 0 \Rightarrow \\ ab - \frac{1}{2 ρ^{2}} b^{2} \leq \frac{1}{2} ρ^{2} a^{2} \Rightarrow \frac{1}{2} ab + \frac{1}{2} ab - \frac{1}{2 ρ^{2}} b^{2} \leq \frac{1}{2} ρ^{2} a^{2} \end{matrix}$ \matrix{ {{\rho ^2}{a^2} + {1 \over {{\rho ^2}}}{b^2} - 2ab \ge 0 \Rightarrow {1 \over 2}{\rho ^2}{a^2} + {1 \over {2{\rho ^2}}}{b^2} - ab \ge 0 \Rightarrow } \cr {ab - {1 \over {2{\rho ^2}}}{b^2} \le {1 \over 2}{\rho ^2}{a^2} \Rightarrow {1 \over 2}ab + {1 \over 2}ab - {1 \over {2{\rho ^2}}}{b^2} \le {1 \over 2}{\rho ^2}{a^2}} \cr }

The following substitutions are carried out: $a = \tilde{d}$ a = \tilde d and b = e^TPL, and the previous relation becomes (47) $\frac{1}{2} {\tilde{d}}^{T} L^{T} Pe + \frac{1}{2} e^{T} PL \tilde{d} - \frac{1}{2 ρ^{2}} e^{T} {PLL}^{T} Pe \leq \frac{1}{2} ρ^{2} {\tilde{d}}^{T} \tilde{d}$ {1 \over 2}{\tilde d^T}{L^T}Pe + {1 \over 2}{e^T}PL\tilde d - {1 \over {2{\rho ^2}}}{e^T}PL{L^T}Pe \le {1 \over 2}{\rho ^2}{\tilde d^T}\tilde d

Eq. (47) is substituted in Eq. (44), and the inequality is enforced, thus giving (48) $\dot{V} \leq - \frac{1}{2} e^{T} Qe + \frac{1}{2} ρ^{2} {\tilde{d}}^{T} \tilde{d}$ \dot V \le - {1 \over 2}{e^T}Qe + {1 \over 2}{\rho ^2}{\tilde d^T}\tilde d

Eq. (48) shows that the H_∞ tracking performance criterion is satisfied. The integration of $\dot{V}$ \dot V from 0 to T gives (49) $\begin{matrix} \int_{0}^{T} \dot{V} (t) dt \leq - \frac{1}{2} \int_{0}^{T} {e‖}_{Q}^{2} dt + \frac{1}{2} ρ^{2} \int_{0}^{T} {\tilde{d}‖}^{2} dt \Rightarrow \\ 2 V (T) + \int_{0}^{T} {e‖}_{Q}^{2} dt \leq 2 V (0) + ρ^{2} \int_{0}^{T} {\tilde{d}‖}^{2} dt \end{matrix}$ \matrix{ {\int_0^T {\dot V\left( t \right)dt} \le - {1 \over 2}\int_0^T {\left\| e \right\|_Q^2dt + {1 \over 2}{\rho ^2}} \int_0^T {{{\left\| {\tilde d} \right\|}^2}dt} \Rightarrow } \cr {2V\left( T \right) + \int_0^T {\left\| e \right\|_Q^2dt \le 2V\left( 0 \right)} + {\rho ^2}\int_0^T {{{\left\| {\tilde d} \right\|}^2}dt} } \cr }

Moreover, if there exists a positive constant M_d > 0 such that $\int_{0}^{\infty} {\tilde{d}‖}^{2} dt \leq M_{d}$ \int_0^\infty {{{\left\| {\tilde d} \right\|}^2}dt \le {M_d}} , then one gets $\int_{0}^{\infty} {e‖}_{Q}^{2} dt \leq 2 V (0) + ρ^{2} Md$ \int_0^\infty {\left\| e \right\|_Q^2dt} \le 2V\left( 0 \right) + {\rho ^2}Md . Thus, the integral $\int_{0}^{\infty} {e‖}_{Q}^{2} dt$ \int_0^\infty {\left\| e \right\|_Q^2dt} is bounded. Moreover, V (T) is bounded, and from the definition of the Lyapunov function V in Eq. (38), it becomes clear that e(t) will also be bounded since e(t) ∈ Ω_e = {e|e^TPe ≤ 2V (0) + ρ²Md}. According to the above and with the use of Barbalat’s Lemma, one obtains lim_t→∞ e(t) = 0.

Through the stages of the stability proof, one arrives at Eq. (48), which shows that the H-infinity tracking performance criterion holds. By selecting the attenuation coefficient ρ to be sufficiently small and in particular to satisfy $ρ^{2} < {e‖}_{Q}^{2} / {\tilde{d}‖}^{2}$ {\rho ^2} < \left\| e \right\|_Q^2/{\left\| {\tilde d} \right\|^2} , one has that the first derivative of the Lyapunov function is upper bounded by 0. This condition holds at each sampling instance, and consequently, global stability for the control loop can be concluded.

7.

Simulation Tests

The global stability properties of the control method and the elimination of the state vector’s tracking error, which were previously proven through Lyapunov analysis, are further confirmed through simulation experiments. In the implementation of the proposed non-linear optimal control method for the MGIM, the algebraic Riccati equation of Eq. (41) has to be solved in each sampling period. Indicative values about the parameters of the dynamic model of the MGIM have been as follows: (i) induction motor: L_d = 1.1mH, L_q = 1.1 mH, τ_s = 5, M = 40.3 mH, σ = 2, μ = 40, J_m = 0.5 kgm², B_m = 0.01, n_p = 4, α = 0.025, β = 0.3, γ = 0.7, J_L = 0.2 kgm², and (ii) magnetic gear: J_g = 0.1 kgm², G_r = 2, B_g = 0.01, p_m = 10, n_L = 20, B_L = 0.01. The obtained results are depicted in Figures 3–6. In the obtained diagrams, the real values of the system’s state vector are depicted in blue colour, the reference setpoints are printed in red colour, while the state estimates, which are provided by the H-infinity Kalman Filter, are plotted in green colour. It can be noticed that in all test cases, fast and accurate tracking of setpoints was achieved by the state variables of the MGIM, and this was done under moderate variations of the control inputs.

To elaborate on the tracking performance and on the robustness of the proposed non-linear optimal control method for the MGIM, Tables 2 and 3 are given which provide information about the accuracy of tracking of the reference setpoints by the state variables of the MGIM.

Table 2.

Tracking RMSE for the MGIM in the disturbance-free case

	RMSE_x1	RMSE_x2	RMSE_x3	RMSE_x4	RMSE_x5	RMSE_x6
Test₁	0.0052	0.0026	0.0064	0.0037	0.0001	0.0002
Test₂	0.0041	0.0020	0.0064	0.0063	0.0002	0.0003

MGIM, magnetically geared induction motor; RMSE, Root Mean Square Error.

Table 3.

Tracking RMSE for the MGIM in the case of disturbances

Δa%	RMSE_x1	RMSE_x2	RMSE_x3	RMSE_x4	RMSE_x5	RMSE_x6
0%	0.0052	0.0026	0.0064	0.0037	0.0001	0.0002
10%	0.0057	0.0029	0.0064	0.0014	0.0001	0.0003
20%	0.0062	0.0031	0.0064	0.0007	0.0001	0.0003
30%	0.0066	0.0033	0.0064	0.0027	0.0002	0.0001
40%	0.0069	0.0035	0.0065	0.0046	0.0002	0.0003
50%	0.0073	0.0036	0.0065	0.0064	0.0002	0.0003
60%	0.0075	0.0038	0.0065	0.0081	0.0002	0.0003

MGIM, magnetically geared induction motor; RMSE, Root Mean Square Error.

8.

Conclusions

A non-linear optimal control method has been proposed for the dynamic model of the MGIM. Using the field-orientation assumption and a description in the asynchronously rotating dq reference frame, the dynamic model of the MGIM has been formulated, and differential flatness properties have been proven about it. To apply the proposed non-linear optimal control method, the dynamic model of the MGIM has undergone approximate linearisation with the use of first-order Taylor-series expansion and through the computation of the associated Jacobian matrices. The linearisation process was taking place at each sampling instance around the temporary operating point (x^*,u^*), where x^* is the present value of the system’s state vector and u^* is the last sampled value of the control inputs vector. For the approximately linearised model of the system, an H-infinity feedback controller was designed. To compute the feedback gains of the H-infinity controller, an algebraic Riccati equation had to be solved repetitively at each time step of the control algorithm. The global stability properties of the control scheme have been proven through Lyapunov analysis.

Nonlinear Optimal Control of Magnetically Geared Induction Motors

Full Article

Paradigm

My account