Robust Model Predictive Control of Magnetic Levitation Bearings

Robust Model Predictive Control of Magnetic Bearings

The application of magnetic bearing technology in the power industry is mainly reflected in new energy power generation, among which wind power generation and flywheel energy storage have extensively researched and applied magnetic bearings. Due to the advantage of frictionlessness, when the excitation rotor of the magnetic bearing is rigidly connected to the impeller of the wind turbine, it can effectively reduce the starting wind speed of the wind turbine, decrease the resistance torque, and improve wind energy utilization and power generation efficiency [1-2]. It also avoids a series of problems caused by the use of lubricating oil in wind turbines, resulting in low maintenance costs. Magnetic suspension technology is also applied in the yaw system, where magnetic suspension drive replaces the traditional multi-gear mechanical drive to improve the yaw wind alignment accuracy [3-4]. Flywheel energy storage with magnetic bearings as the core is a new type of energy storage technology for new energy power systems, featuring no pollution, high efficiency, and long service life [5-6].

The control technology of the magnetic bearing-rotor system is an important basis for the rapid development of magnetic bearings. The mainstream control technologies include PID control, sliding mode control, coordinated control, etc. However, due to the complexity of the magnetic bearing-rotor system, a single control algorithm can hardly maintain the system's control performance, and it is usually necessary to integrate different algorithms, such as fuzzy PID control [7-8], particle swarm optimization-based PID [9], and high-order sliding mode control [10-12].

Model Predictive Control (MPC) is a control method that uses a model to predict the future behavior of the system and solves the optimization problem for specified indicators based on the prediction results. Reference [13] constructed a discrete-time model predictive control based on the linear time-invariant model of magnetic bearings without violating input and state constraints, solving a bounded MPC problem at each sampling moment to calculate the optimal control input. Reference [14] designed a controller for the radial four-degree-of-freedom system of magnetic bearings using the MPC algorithm based on the magnetic bearing MBC500 test platform. Through comparison with experiments and other control methods, it verified the superior system performance of magnetic bearings under MPC. Reference [15] took the unmodeled dynamics as the external disturbance of the magnetic bearing-rotor system, and combined the disturbance observer with H∞ control to effectively improve the anti-disturbance ability of the system. Considering the weak anti-external disturbance ability of MPC, this paper introduces the robust H∞ control method, incorporates the idea of Min-Max robust control into the core link of MPC, i.e., the rolling optimization algorithm, and analyzes the robustness of the magnetic bearing-rotor system under Robust Model Predictive Control (RMPC) through Simulink simulation.

1 Magnetic Bearing-Rotor System Model

The magnetic bearing-rotor system is shown in Figure 1. Wherein, lₐ and lᵦ are the distances from magnetic bearings a and b to the rotor centroid, respectively, and l = lₐ + lᵦ; xₐ and yₐ are the displacements of the rotor at bearing a along the x and y axes, respectively; xᵦ and yᵦ are the displacements of the rotor at bearing b along the x and y axes, respectively; θₓ and θᵧ are the rotation angles of the rotor around the x and y axes, respectively.

Figure 1 Schematic diagram of the magnetic bearing-rotor system

According to Reference [16], the state-space expression of the four-degree-of-freedom radial magnetic bearing is:

Ẋ = AX + BU  (1)

Where:

X = [xₐ, ẋₐ, xᵦ, ẋᵦ, yₐ, ẏₐ, yᵦ, ẏᵦ]ᵀ,

U = [iₓₐ, iₓᵦ, iᵧₐ, iᵧᵦ]ᵀ, Y = [xₐ, xᵦ, yₐ, yᵦ]ᵀ,

A = [0₂ₓ₄  I₂ₓ₄;  A₂₁  0₂ₓ₄], B = [0₄ₓ₄;  B₂₁],

C = [I₄ₓ₄  0₄ₓ₄], D = 0₄ₓ₄,

A₂₁ = [        [0, 0, 0, 0, Kₓ/m, 0, 0, 0],        [0, 0, 0, 0, 0, Kₓ/m, 0, 0],        [lₐKₓ/Jₓ, 0, -lᵦKₓ/Jₓ, 0, 0, 0, 0, 0],        [0, lₐKₓ/Jᵧ, 0, -lᵦKₓ/Jᵧ, 0, 0, 0, 0],        [0, 0, 0, 0, 0, 0, Kᵧ/m, 0],        [0, 0, 0, 0, 0, 0, 0, Kᵧ/m],        [0, 0, 0, 0, lₐKᵧ/Jₓ, 0, -lᵦKᵧ/Jₓ, 0],        [0, 0, 0, 0, 0, lₐKᵧ/Jᵧ, 0, -lᵦKᵧ/Jᵧ]        ]

B₂₁ = [        [Kᵢₓ/m, 0, 0, 0],        [0, Kᵢₓ/m, 0, 0],        [lₐKᵢₓ/Jₓ, -lᵦKᵢₓ/Jₓ, 0, 0],        [0, 0, lₐKᵢₓ/Jᵧ, -lᵦKᵢₓ/Jᵧ],        [0, 0, Kᵢᵧ/m, 0],        [0, 0, 0, Kᵢᵧ/m],        [0, 0, lₐKᵢᵧ/Jₓ, -lᵦKᵢᵧ/Jₓ],        [0, 0, 0, 0]        ]

Where: X is the state vector; U is the control vector; Y is the output vector; iₓₐ, iₓᵦ, iᵧₐ, iᵧᵦ are the input control currents; Kₓ, Kᵧ are the displacement stiffness coefficients; Kᵢₓ, Kᵢᵧ are the current stiffness coefficients; m is the rotor mass; Jₓ, Jᵧ are the moments of inertia of the rotor around the x and y axes, respectively; ω is the angular velocity of the rotor.

2 Robust Model Predictive Control

This paper is based on Min-Max robust model predictive control, taking the system performance index under the "worst-case" caused by uncertain factors as the optimization objective. The input and output constraint conditions of the system are transformed into Linear Matrix Inequalities (LMI), and an LMI optimization problem is constructed to solve the state feedback control law and minimize the maximization problem in LMI.

The Euler method is used to discretize the state-space expression of the four-degree-of-freedom radial magnetic bearing in Equation (1), with a sampling time Tₛ = 0.001 s. External disturbance is added to the system model as an input variable, and the expression is:

x(k+1) = A₁x(k) + B₁u(k) + B₂w(k)
       y(k) = C₁x(k) + D₁₂w(k)
       z(k) = C₂x(k) + D₂₁u(k) + D₂₂w(k)  (2)

Where: A₁ = I + TₛA, B₁ = TₛB, C₁ = C, D₁₁ = D, and A, B, C, D are all discrete matrices; x is the system state; u is the control current input; w is all forms of external disturbance input; y is the system constraint output; z is the control performance output; A₁, B₁, B₂, C₁, C₂, D₁₁, D₁₂, D₂₁, D₂₂ are all parameter matrices (determined according to actual conditions).

Within the prediction horizon N at time k, the control input constraints, output constraints, and state constraints of the system are:

uₘᵢₙ ≤ u(k+i) ≤ uₘₐₓ,  i = 0,1,...,N-1
       yₘᵢₙ ≤ y(k+i) ≤ yₘₐₓ,  i = 1,2,...,N
       xₘᵢₙ ≤ x(k+i) ≤ xₘₐₓ,  i = 1,2,...,N  (3)

Since external disturbance input generally cannot directly affect the system constraint output, D₁₂ = 0 is set. The design idea of robust H∞ control is to satisfy the internal closed-loop feedback stability of the system without violating the system output constraints. The optimization objective is to satisfy that the H∞ norm from w to z is minimized (i.e., less than the positive number γ) without violating the system output constraints:

||T_zw||∞ < γ; γ > 0  (4)

Let the state feedback control law be u = Kx, substitute it into Equation (2) and rearrange to get:

x(k+1) = A_cl x(k) + B₂ w(k)
       y(k) = C₁_cl x(k) + D₁₂ w(k)
       z(k) = C₂_cl x(k) + D₂₁_cl w(k)  (5)

Where: A_cl = A₁ + B₁K, C₁_cl = C₁ + D₁₁K, C₂_cl = C₂ + D₂₁K. Then there exists a positive definite symmetric matrix P > 0 that satisfies Equation (6):

[        [P⁻¹, 0, A_clᵀ, C₂_clᵀ],        [0, I, B₂ᵀ, D₂₁_clᵀ],        [A_cl, B₂, P, 0],        [C₂_cl, D₂₁_cl, 0, I]        ] > 0  (6)

According to the Schur complement formula, Equation (6) is equivalent to Equation (7) (the matrix expansion form is omitted here, see the original text for detailed derivation). Multiply the matrix of Equation (7) by [x(k)ᵀ w(k)ᵀ] on the left and [x(k) w(k)]ᵀ on the right, respectively, and rearrange to get:

x(k+1)ᵀP x(k+1) + z(k)ᵀz(k) ≤ x(k)ᵀP x(k) + γ² w(k)ᵀw(k)  (8)

Construct the Lyapunov function V(x) = xᵀPx using P, and obtain the dissipative inequality:

V(x(k+1)) + Σₖⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σₖⁿ⁻¹ w(i)ᵀw(i) + V(x(k))  (9)

Assuming that the system can satisfy Equation (9) within the entire prediction horizon, we have:

V(x(n)) + Σ₀ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σ₀ⁿ⁻¹ w(i)ᵀw(i) + V(x(0))  (10)

Then the dissipative inequality of Equation (11) holds:

V(x(k+1)) + Σₖᵏ⁺ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σₖᵏ⁺ⁿ⁻¹ w(i)ᵀw(i) + V(x(k))  (11)

Since x(0) = 0 and V(x) ≥ 0, according to Equation (11):

Σ₀ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σ₀ⁿ⁻¹ w(i)ᵀw(i)  (12)

Equation (12) proves that the H∞ norm from w to z is less than γ.

Let Q = P⁻¹, multiply the matrix of Equation (6) by the block diagonal matrix diag{Q, I, Q, I} on the left and right, respectively, to get Equation (13); then let Y = KQ, substitute it into Equation (13) to get Equation (14) (the matrix expansion form is omitted here, see the original text for details).

Assuming that the system's external disturbance input is bounded, we have:

Σₖᵏ⁺ⁿ⁻¹ w(i)ᵀw(i) ≤ n wₘₐₓ²  (15)

Where: wₘₐₓ is the maximum value of the external disturbance input. Given α > 0, if the system state belongs to an elliptical domain Ω₁, i.e.:

Ω₁: {x ∈ Rⁿ | γ² n wₘₐₓ² + V(x) ≤ α}  (16)

V(x) = xᵀPx, Ω₁ is defined by P, α, and wₘₐₓ. The system's closed-loop state trajectory belongs to another elliptical domain Ω₂, which is defined by P and α:

Ω₂: {x ∈ Rⁿ | V(x) ≤ α}  (17)

According to Equation (3) and the Cauchy-Schwarz inequality, the control input constraints and output constraints can be written as:

||u(k)||∞ = ||Kx(k)||∞ ≤ √(λₘₐₓ(KᵀK) λₘₐₓ(x(k)ᵀx(k))) ≤ √(λₘₐₓ(KᵀK) α / λₘᵢₙ(P)) ≤ uₘₐₓ  (18)

||y(k)||∞ = ||(C₁ + D₁₁K)x(k)||∞ ≤ √(λₘₐₓ((C₁ + D₁₁K)ᵀ(C₁ + D₁₁K)) α / λₘᵢₙ(P)) ≤ yₘₐₓ  (19)

Where: K is the state feedback matrix; α is a positive number.

The system constraint conditions can be expressed as:

uₘₐₓ² I - α Q⁻¹ KᵀK Q⁻¹ ≥ 0  (20)

Assume that there exists a symmetric matrix X whose diagonal elements satisfy Xᵢᵢ ≤ uₘₐₓ² (i=1,2,⋯,n); and a symmetric matrix Z whose diagonal elements satisfy Zᵢᵢ ≤ yₘₐₓ² (i=1,2,⋯,n). According to the Schur complement formula and Equation (13), the LMI forms of the system's control input constraints and output constraints are derived:

[        [X, α Yᵀ],        [α Y, Q]        ] ≥ 0,  i=1,2,⋯,n  (21)

[        [Z, (C₁Q + D₁₁Y)ᵀ],        [C₁Q + D₁₁Y, Q]        ] ≥ 0  (22)

At time k, the system state is x(k) ∈ Ω₁, then:

γ² n wₘₐₓ² + x(k)ᵀQ⁻¹ x(k) ≤ α  (23)

Combining Equations (8), (14), (16), and (17), the LMI optimization problem is obtained:

min γ²
       s.t. Equations (8), (14), (16), (17)  (24)

At each sampling moment k (k=0,1,2,⋯), solve the LMI problem. If there exists a solution (γₖ², Qₖ, Yₖ, Xₖ, Zₖ), the state feedback matrix at time k can be expressed as:

K(k) = Yₖ Qₖ⁻¹  (25)

The flow of MPC is shown in Figure 2. The control structure of RMPC is similar to that of MPC, and the optimization purpose is to incorporate the idea of Min-Max robust control into the core link of MPC, i.e., the rolling optimization algorithm.

Figure 2 Flow chart of Model Predictive Control

The rolling optimization algorithm of RMPC is as follows:

Initialization: At k=0, initialize the system state.

Iterative calculation and optimization: At time k, solve the LMI optimization problem in Equation (24). If a solution exists, obtain the state feedback matrix K(k), calculate the system control input uₖ = K(k)x(k) and apply it to the system; if there is no solution to the LMI optimization problem, increase the value of α and solve again; if there is still no solution after continuously increasing α, replace it with the previous state feedback matrix K(k-1), i.e., K(k) = K(k-1).

Set k = k+1, and repeat Step 2.

The parameters of the magnetic bearing-rotor system are shown in Table 1.

Table 1 Parameters of the magnetic bearing-rotor system

3 Simulation Analysis

Based on the model predictive control algorithm in Reference [16], an anti-disturbance simulation is carried out for the magnetic bearing-rotor system. The reference value is set to [0.1; 0.1; 0.1; 0.1] (unit: mm), the system input constraint is [2.5; 2.5; 2.5; 2.5] (unit: A), and the system output constraint is [0.25; 0.25; 0.25; 0.25] (unit: mm). At 0.5 s, a unit pulse signal is applied to the rotor at magnetic bearing a as an external disturbance, with a disturbance duration of 0.1 s, causing the rotor to generate displacements along the x and y axes within the safe air gap range. The same parameters are used for the two control methods, and the simulation results of MPC are as follows.

Figure 3 Rotor displacement and control current under MPC

The system adjustment time of MPC is about 0.2 s, and the response curve has an overshoot of about 10%. Due to the input constraint setting, the system's control input current will not exceed 2.5 A. After applying the disturbance at 0.5 s, the rotor shifts in all directions, with a maximum shift of about 0.12 mm. After 0.15 s, the rotor displacement tends to be stable, but the system does not fully return to the stable state before the disturbance, and the rotor exhibits low-frequency vibration, i.e., the rotor is in a weak conical motion state, and the rotor does not reach stability in the two rotational degrees of freedom, indicating that the system has weak anti-disturbance ability under MPC.

Since the established mathematical model of the magnetic bearing-rotor system is not an accurate mathematical model but a properly simplified linear model (nominal model), and MPC is a model-based feedback correction control system, it cannot guarantee to eliminate the influence of external disturbances quickly and accurately. This paper improves the system's anti-disturbance ability by adding robust control. According to the discrete linear model of the system in Equation (2), set:

B₂ = D₁₂ = D₂₂ = I  (26)

The system control input constraint is ||u||∞ ≤ 2.5 A. When the magnetic bearing-rotor system works normally, the air gap between the rotor and the stator is 0.35 mm. Considering the safe air gap margin, the system output constraint is set to ||y||∞ ≤ 0.25 mm; the initial value of α is set to α₀ = 0.1. Based on RMPC, the system is simulated again with disturbances, and disturbances of the same magnitude but opposite directions are added at 0.5 s and 0.7 s, respectively. The simulation results are shown in Figure 5.

Figure 4 Rotor displacement and control current under RMPC

Under disturbance input, the RMPC system has no overshoot, indicating that the system has good dynamic stability performance and can quickly return to a stable state after being subjected to external disturbances. The rotor displacement change is less than 0.11 mm, and the system can fully return to the state before the disturbance within 0.1 s without low-frequency vibration, which effectively suppresses the influence of disturbances on the system and improves the robustness of the system. Comparing the system input currents under the two control modes, the control current curves of RMPC almost completely overlap, indicating that the magnetic bearing-rotor system under RMPC has stronger self-balancing adjustment ability after being subjected to external disturbances.

4 Conclusion

This paper establishes the mathematical model of the magnetic bearing-rotor system. Based on the idea of robust control, under the framework of the LMI optimization problem, the H∞ norm from the system's external disturbance input to the control performance output is taken as the optimization index to solve the state feedback control law, and the matrix inequality conditions that can satisfy the H∞ performance are derived. Combined with model predictive control, the rolling optimization algorithm is improved by replacing the quadratic optimization problem repeatedly performed in the original rolling optimization with a linear matrix inequality optimization problem, and the system's anti-disturbance performance is improved by using the system's control input constraints and output constraints as means.