Model predictive control of a magnetic levitation system using two-level state feedback

Abstract

The magnetic levitation system is a critical part to guarantee safe and reliable operations of a maglev train. In this paper, the control strategy is proposed for the magnetic levitation system based on the model predictive control incorporating two-level state feedback. Taking advantage of the measurable state variables, that is, air gap, electromagnet acceleration, and control current through high-resolution sensor measurement, the first-level nonlinear state feedback is to linearize the unstable nonlinear magnetic levitation system, and the second-level linear state feedback is to further stabilize the system and improve the dynamic performances, which together provide a stable prediction model. The simulation results demonstrate that the proposed control strategy can ensure high-precision air gap control and favorable disturbance resistance ability.

Introduction

The electromagnetic suspension system for a magnetic levitation (maglev) train utilizes the electromagnetic attractive forces between the train body-equipped levitation electromagnets and their above electromagnetic components to make the train body levitated above the guide rails for train longitudinal motions. The distance between the electromagnets and the electromagnetic components is generally maintained around 8–12 mm, which is called the air gap. The active air gap control is the operation precondition and pivotal part for a maglev train.

Some studies have been made about magnetic levitation control. Sinha et al. ¹ employed neural networks to design a levitation controller. MacLeod and Goodall ² developed a linear–quadratic–Gaussian (LQG)-based optimal control using a frequency-domain weighting approach. Tsunashima et al. ³ derived the stability criteria for magnetic levitation control and analyzed the effect of a dead zone on the control performance. Lin and Gau ⁴ applied the fuzzy control approach to levitation control. Bittar and Sales ⁵ developed the H₂ and H_∞ controller of a levitation system. Sinha and Pechev ⁶ employed the adaptive control to design a levitation controller. Al-Muthairi and Zribi ⁷ designed the slide mode controller using the continuous approximation approach of switching control. Shen et al. ⁸ established the vertical and lateral control model of a high-speed maglev train. Molero et al. ⁹ put forward the slide mode control of an electromagnetic levitation system using exact linearization. Han et al. ¹⁰ established a multi-body dynamic model of a maglev vehicle to predict the ride quality, including the effects of system parameters of levitation control. Zhou and Duan ¹¹ proposed a novel levitation control to remove the interactions among redundant and general actuators. Zhang et al. ¹² performed the nonlinear analysis of a maglev control system with time-delayed feedback. Yang et al. ¹³ employed a high-order sliding-mode differentiator to attenuate high-order matched and mismatched disturbances for maglev suspension vehicles. Xu et al. ¹⁴ proposed a robust levitation control strategy with guaranteed magnitude limitation of the air gap under uncertain conditions. Duka et al. ¹⁵ put up an internal model control-based proportional–integral–derivative (PID) controller using a linear model to represent the nonlinear dynamics of a magnetic levitation system. Klaučo et al. ¹⁶ proposed an inner proportional–summation–difference controller for a magnetic levitation system, but using model predictive control (MPC) as a reference governor. However, the direct external disturbances, for example, being originated from a crowd of passengers’ boarding and alighting a maglev train, are less involved in the current studies.

MPC has proved in practice as an applicable and popular control approach in industry.^17–20 It employs a prediction model to predict plant performances in the prediction horizon under a series of presumed control inputs in the control horizon. The optimization control objective is to make the predicted outputs approach the reference inputs in the prediction horizon. Through solving the optimization problem, the control sequences can be found out. The first control input is actually imposed on the plant at the beginning of next period. At this instant, the actual output is sampled, the new prediction and optimization processes will be performed again. Compared with other control approaches, the control mechanism of prediction, optimization, and feedback over the rolling horizon makes the control performances have preferable tracking and disturbance-resisting abilities. MPC can directly deal with disturbances and have many successful industry applications, which drives us to develop a potentially implementable model predictive controller. MPC depends on a stable prediction model; however, the magnetic levitation system is unstable and nonlinear. Under these circumstances, nonlinear state feedback technique ²¹ is incorporated into the control strategy to linearize the unstable nonlinear levitation system; furthermore, a linear state feedback is to, on one hand, stabilize the system and, on the other hand, to improve the dynamic performances. Through the two-level state feedback, a stable prediction model is provided to MPC which globally ameliorates the transient performances especially under disturbances. This paper proposes a control scheme to apply MPC to an unstable plant.

The remainder of this paper is organized as follows. At first, the magnetic levitation system is modeled. And then, the state feedback process and the control performances are elaborated. Furthermore, the MPC control strategy incorporating two-level state feedback is outlined for the magnetic levitation system, and the control performances are demonstrated. Finally, the conclusions are drawn.

Levitation system model

Levitation system structure

Figure 1 demonstrates the structure of a levitation system, which includes vehicle body, rail and electromagnetic component, electromagnet, coil, sensors, levitation controller, and chopper. The air gap is detected through the eddy current sensor with a resolution of micron level. The gap sensor is based on the principle of electromagnetic induction, that is, the alternately varied current through the coils of an electromagnet can engender eddy current through a closed loop. The vertical movement velocity and acceleration are acquired through acceleration sensors, which utilize the movements of a reference object to produce a voltage proportional to the acceleration. The velocity is an integral of the acquired acceleration. The control current is sampled by current sensors through the principle of electromagnetic induction. The gap, acceleration, and current signals are fed back to the levitation controller. The chopper is the actuator of levitation control system with great voltage and current outputs. The air gap is controlled through adjusting the output voltage or current of the chopper. When the electromagnet coil is input with current, the generated attractive force counteracts the weight of vehicle body to make it move upward. If the gap is greater than the set point, for example, 1 cm, the input voltage or current will be increased, otherwise decreased.

Figure 1.

Levitation system structure.

Plant model

Figure 2 shows the forces a magnetic levitation system is subject to. In Figure 2, $F(i,\delta )$ is the electromagnetic attractive force, m is the mass of vehicle body, g is the gravity acceleration, and f_d is the external disturbance. $\delta (t)$ is the air gap, ${\mathrm{\Phi }}_m$ is the magnetic flux through the air gap, ${\mathrm{\Phi }}_T$ is the main magnetic flux, ${\mathrm{\Phi }}_L$ is the leak magnetic flux, and $A_m$ is the polar area of electromagnet core. $u(t)$ is the control voltage, and $i(t)$ is the control current.

Figure 2.

Force analysis.

According to the Ampere circuital theorem, we have

$$2\delta \left(t\right)\frac{B_m}{{\mu }_0}=Ni\left(t\right)$$

(1)

where $B_m$ is the gap magnetic flux density, ${\mu }_0$ is the gap magnetic conductivity, and N is the number of turns of electromagnet coils. Ignoring the leak magnetic flux ${\mathrm{\Phi }}_L$ , the gap magnetic density is expressed as

$$B_m=\frac{{\mathrm{\Phi }}_m}{A_m}=\frac{{\mathrm{\Phi }}_T}{A_m}=\frac{{\mu }_{0}Ni\left(t\right)}{2\delta \left(t\right)}$$

(2)

At instant t, the electromagnetic attractive force is

$$F(i,\delta )=\frac{B_m^2{A}_m}{{\mu }_0}=\frac{{\mu }_0{N}^2{A}_m}{4}{\left[\frac{i\left(t\right)}{\delta \left(t\right)}\right]}^2$$

(3)

According to Kirchhoff’s electrical law, the voltage equilibrium equation of electromagnet coil loop is established as follows

$$\begin{array}{c}\multicolumn{1}{c}{u\left(t\right)=Ri\left(t\right)+\frac{d\left(Li\left(t\right)\right)}{dt}}\\ \multicolumn{1}{c}{=Ri\left(t\right)+\frac{d}{dt}\left[N{\mathrm{\Phi }}_T(i,\delta )\right]}\\ \multicolumn{1}{c}{=Ri\left(t\right)+\frac{{\mu }_0{N}^2{A}_m}{2\delta \left(t\right)}\frac{\mathrm{di}\left(t\right)}{dt}-\frac{{\mu }_0{N}^2{A}_{m}i\left(t\right)}{2\delta {\left(t\right)}^2}\frac{d\delta \left(t\right)}{dt}}\end{array}$$

(4)

where R is the coil resistance.

Specify the downward force is positive. The equilibrium condition of vertical forces is described as

$$m\frac{d^2\delta \left(t\right)}{d{t}^2}=\mathrm{mg}-F(i,\delta )+f_d\left(t\right)$$

(5)

Define the state variables $s=[s_1{s}_2{s}_3{]}^{\mathrm{T}}$ , where $s_1=\delta ,s_2={\stackrel{·}{s}}_1,\mathrm{and}{s}_3=i$ . Equations (3)–(5) jointly yield the following state-space equation

$$\{\begin{array}{c}\multicolumn{1}{c}{{\stackrel{·}{s}}_1=s_2}\\ \multicolumn{1}{c}{{\stackrel{·}{s}}_2=g-\frac{{\mathrm{ks}}_3^2}{{\mathrm{ms}}_1^2}+\frac{f_d}{m}}\\ \multicolumn{1}{c}{{\stackrel{·}{s}}_3=\frac{s_2{s}_3}{s_1}-\frac{R{s}_1{s}_3}{2k}+\frac{s_1}{2k}u}\end{array}$$

(6)

where $k={\mu }_0{N}^2{A}_m/4$ .

The output equation is defined as

$$y=s_1$$

(7)

Plant stability

Equation (6) is linearized around the equilibrium point. When the plant approaches the equilibrium state $\stackrel{·}{s}=0$ , that is, the speed and the acceleration of an electromagnet are zero, and the current variation rate is zero. The equilibrium point $(s_0,u_0)$ , where $s_0={[{\delta }_{0}0i_0]}^T$ , is solved as

$$\{\begin{array}{c}\multicolumn{1}{c}{s_1={\delta }_0=\frac{u_0}{R}\sqrt{\frac{k}{\mathrm{mg}}}}\\ \multicolumn{1}{c}{s_2=0}\\ \multicolumn{1}{c}{s_3=i_0=\frac{u_0}{R}}\\ \multicolumn{1}{c}{u=u_0}\end{array}$$

(8)

Define $s\prime =\mathrm{\Delta }s=s-s_0,u\prime =\mathrm{\Delta }u=u-u_0,\mathrm{and}y\prime =y-{\delta }_0$ . At the equilibrium point $(s_0,u_0)$ , equations (6) and (7) are linearized into

$$\{\begin{array}{c}\stackrel{·}{s}\prime =A\prime s\prime +B\prime u\prime \hfill \\ \multicolumn{1}{c}{y\prime =C\prime x\prime }\end{array}$$

(9)

where

$$A\prime =\left(\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{\frac{2g}{{\delta }_0}}& 0\hfill & -\frac{2}{{\delta }_0}\sqrt{\frac{\mathrm{kg}}{m}}\hfill \\ \multicolumn{1}{c}{0}& \sqrt{\frac{\mathrm{mg}}{k}}\hfill & -\frac{R{\delta }_0}{2k}\hfill \end{array}\right),B\prime =\left(\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{\frac{{\delta }_0}{2k}}\end{array}\right),\mathrm{and}C\prime =\left(\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right)$$

Through solving the roots of the characteristic equation $det(\lambda I-A\prime )=0$ , we can find that the linearized system has a root locating at the right-half s-plane. According to the Lyapunov stability principles, the original nonlinear plant is not stable, which implies that the plant should be stabilized at first through local feedback. In terms of the plant described by equation (6), the state variables are measurable. Therefore, state feedback is employed to stabilize the plant, which will be addressed in the following section.

State feedback

State feedback process

In order to linearize the plant, new measurable state variables are defined as x = [x₁x₂x₃]^T, where $x_1=s_1,x_2=s_2,\mathrm{and}{x}_3=g-{\mathrm{ks}}_3^2/({\mathrm{ms}}_1^2)$ . x₁, x₂, and x₃ have the physical meanings. x₁ stands for the air gap, x₂ represents the electromagnet movement speed, and x₃ has the acceleration unit because if f_d = 0, actually $x_3={\stackrel{·}{s}}_2$ , which indirectly describing the coil current s₃. For the state feedback, at first let f_d = 0, but the plant pole placement will aid in resisting the disturbance afterward. According to the characteristics of equation (6), the following state-space description is defined

$$\{\begin{array}{c}{\stackrel{·}{x}}_1=s_2\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_2=s_3}\\ \multicolumn{1}{c}{{\stackrel{·}{x}}_3=\alpha \left(s\right)+\beta \left(s\right)u}\end{array}$$

(10)

From equations (6) and (10), it follows that

$$\{\begin{array}{c}\alpha \left(s\right)=\frac{{\mathrm{Rs}}_3^2}{m{s}_1}\hfill \\ \multicolumn{1}{c}{\beta (s)=-\frac{s_3}{m{s}_1}}\end{array}$$

(11)

Applying the state feedback $u=(-\alpha (s)+v)/\beta (s)$ , where v is the system input after introducing this state feedback, the following linearization system can be obtained

$$\{\begin{array}{c}\stackrel{·}{x}=A_{1}x+B_{1}v\hfill \\ \multicolumn{1}{c}{y=C_{1}x}\end{array}$$

(12)

where

$$A_1=\left(\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill \end{array}\right),B_1=\left(\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{1}\end{array}\right),\mathrm{and}{C}_1=\left(\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right)$$

Because $\mathrm{rank}[B_1{A}_1{B}_1{A}_1^2{B}_1]=3$ , the system after the first-level state feedback is controllable. At the equilibrium point $x_0={[{\delta }_0,0,0]}^T$ , we further introduce the second-level state feedback to improve the transient performances as $v=k_1(x_1-{\delta }_0)+k_2{x}_2+k_3{x}_3+w$ , where w is the system input after the second-level state feedback. The resulted characteristic equation is calculated as

$$det(\lambda I-A_1-B_{1}K)={\lambda }^3-k_3{\lambda }^2-k_2\lambda -k_1=0$$

(13)

where $K=[k_1{k}_2{k}_3]$ . If the expected pole placement is as ${\lambda }_{1,2,3}={\lambda }_1,{\lambda }_2,{\lambda }_3$ , we have

$$\{\begin{array}{c}{k}_1={\lambda }_1{\lambda }_2{\lambda }_3\hfill \\ \multicolumn{1}{c}{k_2=-{\lambda }_1{\lambda }_2-{\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_3}\\ \multicolumn{1}{c}{k_3={\lambda }_1+{\lambda }_2+{\lambda }_3}\end{array}$$

(14)

Consequently, the ultimate state feedback is described as

$$\begin{array}{c}\multicolumn{1}{c}{u=-\frac{m{s}_1}{s_3}v+R{s}_3}\\ \multicolumn{1}{c}{=-\frac{m{s}_1}{s_3}[k_1(s_1-{\delta }_0)+k_2{s}_2+k_3(g-\frac{{\mathrm{ks}}_3^2}{{\mathrm{ms}}_1^2})+w]}\\ \multicolumn{1}{c}{+R{s}_3}\end{array}$$

(15)

Substituting equation (15) into (10), the stabilized control object is described as

$$\{\begin{array}{c}\stackrel{·}{x}=(A_1+B_{1}K)x+B_{1}w\hfill \\ \multicolumn{1}{c}{y=C_{1}x}\end{array}$$

(16)

State feedback performances

The simulation parameters are configured as: m = 750 kg, R = 0.7 Ω, N = 350, A_m = 0.0151 m², δ₀ = 1 cm, and μ₀ = 4π × 10⁻⁷ H/m (Henry per meter). Suppose the overshoot $\sigma \le 15\%$ and the settling time $t_s\le 0.25\mathrm{s}$ , we can obtain the dominant poles ${\lambda }_{1,2}=-18.6\pm i24.8$ . The third pole is empirically determined such that it has the distance to the imaginary axis 5–10 times of that a dominant pole has. Hence, it is set as ${\lambda }_3=-186$ . According to equation (14), K is calculated as $K=[-172,980-7660-220{]}^T$ . The initial gap is configured as $\delta =1.6\mathrm{cm}$ . At 0.5 s, a 1000 N downward disturbance force is imposed, and it is canceled at 1 s. The air gap is controlled around 1 cm.

Figure 3 demonstrates the state feedback performances. Because the closed-loop poles locate at the left-half s-plane, the air gap can stabilize around the set point 1 cm with a zero steady-state error during 0.215–0.5 s from the initial gap 1.6 cm. Both speed and acceleration approach zero during 0.215–0.5 s. The control current and the voltage maintain around 35.6 A and 24.9 V during 0.215–0.5 s, respectively. However, during 0.5–1.0 s when 1000 N downward disturbance force is imposed on the system, the air gap stabilizes around 1.17 cm during 0.685–1.0 s. Both speed and acceleration also approach zero during 0.685–1.0 s. The control current and the voltage keep around 44.3 A and 31.0 V during 0.685–1.0 s, respectively. During 1.0–1.5 s when the disturbance is canceled, the air gap, the control current, and the voltage can resume the steady-state values during 1.19–1.5 s. Figure 3 implies that the state feedback can make the plant stable and has a certain degree of tracking capability, but the performance of resisting disturbances requires further improving.

Figure 3.

State feedback performances of a magnetic levitation system: (a) gap, (b) speed, (c) acceleration, (d) current, and(e) voltage.

MPC

System structure

In this section, we design the MPC controller to further improve the control performance at resisting disturbances. The total control system structure is demonstrated in Figure 4. The control objective of the MPC controller is to optimize the control sequences in the control horizon such that the outputs Y_p(k) predicted by the prediction model converge the reference inputs Y_r(k) in the prediction horizon. w(k) is the controller output at instant kT_s, where T_s is a sampling period and is omitted hereinafter for description simplicity. w(k) is transferred to the plant according to equation (15). The state feedback is performed according to equation (15) to form a stabilized linear system represented by equation (16).

Figure 4.

Control system structure of a magnetic levitation system.

Prediction model

Given a sampling period T_s, the stabilized system described by equation (16) can be discretized into

$$\{\begin{array}{c}x(k+1)=\mathrm{Ax}\left(k\right)+\mathrm{Bw}\left(k\right)+\mathrm{Gd}\left(k\right)\hfill \\ \multicolumn{1}{c}{y\left(k\right)=\mathrm{Cx}\left(k\right)}\end{array}$$

(17)

where the disturbance d(k) is directly added according to equation (6). $A=e^{(A_1+B_{1}K)T_s},B={\int }_0^{T_s}{e}^{(A_1+B_{1}K)t}dt{B}_1,\mathrm{and}C=C_1$ . $x(k)=[x_1(k)x_2(k)x_3(k){]}^{\mathrm{T}}\in R^3\mathrm{and}y(k)\in R^1$ , where R stands for the real number domain, and k denotes an instant. $A\in R^{3\times 3},B\in R^3,\mathrm{and}C\in R^{1\times 3}$ . w(k) is the input of the stabilized system or the output of the MPC controller. According to equations (6) and (16), the disturbance part is modeled as G = [0 1/m 0]^T and d(k) = f_d.

The prediction horizon is P, and the control horizon is M, where $M\le P$ . The integer i represents time increment and gets a value from 1 to P, which is denoted as i = 1, 2,…, P. At instant k, the current state variables x(k), the input w(k − 1) at instant k − 1, and the disturbance d(k − 1) at instant k − 1 are known for prediction. The disturbances $\hat{d}(k+i-1)$ are predicted using the current disturbance input d(k − 1). ²¹ Provided the predicted inputs $\hat{w}(k+i-1)$ , the state variables $\hat{x}(k+i)$ are predicted as

$$\hat{x}(k+i)=A^{i}x\left(k\right)+[A^{i-1}{A}^{i-2}\cdots I]B\left[\begin{array}{c}\hat{w}\left(k\right)\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\hat{w}(k+P-1)}\end{array}\right]+[A^{i-1}{A}^{i-2}\cdots I]\mathrm{Gd}(k-1)$$

(18)

where $\hat{w}(k+i-1)=\hat{w}(k+M-1)$ when $M\le i\le P$ .

Define the input augment as $\mathrm{\Delta }\hat{w}(k+i-1)=\hat{w}(k+i-1)-\hat{w}(k+i-2)$ , where $\hat{w}(k-1)=w(k-1)$ . Consequently, the predicted inputs are transformed into

$$\begin{array}{c}\multicolumn{1}{c}{\hat{w}\left(k\right)=\mathrm{\Delta }\hat{w}\left(k\right)+w(k-1)}\\ \multicolumn{1}{c}{\hat{w}(k+1)=\mathrm{\Delta }\hat{w}(k+1)+\mathrm{\Delta }\hat{w}\left(k\right)+w(k-1)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\hat{w}(k+M-1)=\mathrm{\Delta }\hat{w}(k+M-1)+\cdots }\\ \multicolumn{1}{c}{+\mathrm{\Delta }\hat{w}\left(k\right)+w(k-1)}\end{array}$$

(19)

Substituting equation (19) into (18), we can obtain the following the matrix form

$$X_p\left(k\right)=A_{P}x\left(k\right)+B_{\mathrm{\Delta }}\mathrm{\Delta }{w}_p\left(k\right)+B_{w}w(k-1)+B_{d}d(k-1)$$

(20)

where

$$\begin{array}{c}\multicolumn{1}{c}{X_P\left(k\right)=\left[\begin{array}{c}\begin{array}{c}\hat{x}(k+1)\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\hat{x}(k+M)}\end{array}\hfill \\ \multicolumn{1}{c}{\hat{x}(k+M+1)}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\hat{x}(k+P)}\end{array}\right],\mathrm{\Delta }{w}_P\left(k\right)=\left[\begin{array}{c}\mathrm{\Delta }\hat{w}\left(k\right)\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\mathrm{\Delta }\hat{w}(k+M-1)}\end{array}\right],A_P=\left[\begin{array}{c}A\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{A^M}\\ \multicolumn{1}{c}{A^{M+1}}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{A^P}\end{array}\right]}\\ \multicolumn{1}{c}{B_{\mathrm{\Delta }}=\left[\begin{array}{ccc}B\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}{A}^{i}B}& \cdots \hfill & B\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^M{A}^{i}B}& \cdots \hfill & \mathrm{AB}+B\hfill \\ \multicolumn{1}{c}{⋮}& ⋮⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}{A}^{i}B}& \cdots \hfill & \sum _{i=0}^{P-M}{A}^{i}B\hfill \end{array}\right],B_w=\left[\begin{array}{c}B\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}{A}^{i}B}\\ \multicolumn{1}{c}{\sum _{i=0}^M{A}^{i}B}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}{A}^{i}B}\end{array}\right],\mathrm{and}{B}_d=\left[\begin{array}{c}G\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}{A}^{i}G}\\ \multicolumn{1}{c}{\sum _{i=0}^M{A}^{i}G}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}{A}^{i}G}\end{array}\right]}\end{array}$$

In equation (20), $X_p(k)\in R^{3P},\mathrm{\Delta }{w}_P(k)\in R^M,A_P\in R^{3P\times 3},B_{\mathrm{\Delta }}\in R^{3P\times M},B_w\in R^{3P},\mathrm{and}{B}_d\in R^{3P}$ .

Using equation (20), we can obtain the predicted output $\hat{y}(k+i)$ at instant k + i, which is denoted as

$$Y_P\left(k\right)=\mathrm{\Psi }x\left(k\right)+\mathrm{\Theta }\mathrm{\Delta }w\left(k\right)+{\mathrm{\Gamma }}_{w}w(k-1)+{\mathrm{\Gamma }}_{G}d(k-1)$$

(21)

where

$$\begin{array}{c}\multicolumn{1}{c}{Y_P\left(k\right)=\left[\begin{array}{c}\hat{y}(k+1)\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\hat{y}(k+p)}\end{array}\right],\mathrm{\Psi }=\left[\begin{array}{c}\mathrm{CA}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{C{A}^P}\end{array}\right]}\\ \multicolumn{1}{c}{\mathrm{\Theta }=\left[\begin{array}{ccc}\mathrm{CB}\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}C{A}^{i}B}& \cdots \hfill & \mathrm{CB}\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{M}C{A}^{i}B}& \cdots \hfill & \mathrm{CAB}+\mathrm{CB}\hfill \\ \multicolumn{1}{c}{⋮}& ⋮⋮\hfill & ⋮\hfill \\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}C{A}^{i}B}& \cdots \hfill & \sum _{i=0}^{P-M}C{A}^{i}B\hfill \end{array}\right],{\mathrm{\Gamma }}_w=\left[\begin{array}{c}\mathrm{CB}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}C{A}^{i}B}\\ \multicolumn{1}{c}{\sum _{i=0}^{M}C{A}^{i}B}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}C{A}^{i}B}\end{array}\right],\mathrm{and}{\mathrm{\Gamma }}_G=\left[\begin{array}{c}\mathrm{CG}\hfill \\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{M-1}C{A}^{i}G}\\ \multicolumn{1}{c}{\sum _{i=0}^{M}C{A}^{i}G}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{\sum _{i=0}^{P-1}C{A}^{i}G}\end{array}\right]}\end{array}$$

In equation (21), $Y_p(k)\in R^P,\mathrm{\Psi }\in R^{P\times 3},\mathrm{\Theta }\in R^{P\times M},{\mathrm{\Gamma }}_w\in R^P,\mathrm{and}{\mathrm{\Gamma }}_G\in R^P$ .

Optimization objective

The optimization objective is to make the predicted outputs Y_p(k) approach the reference inputs Y_r(k) in the prediction horizon P, which is denoted as

$$minJ\left(k\right)=\Vert Y_p\left(k\right)-Y_r\left(k\right){\Vert }_Q^2+\Vert \mathrm{\Delta }{w}_P\left(k\right){\Vert }_{R_{\mathrm{\Delta }}}^2$$

(22)

where $Y_r(k)=[y_r(k+1)\dots y_r(k+P){]}^T$ . $Q\mathrm{and}{R}_{\mathrm{\Delta }}$ are the weighting matrices, which are represented as $Q=\mathrm{diag}(Q(1),Q(2),\dots ,Q(P))\mathrm{and}{R}_{\mathrm{\Delta }}=\mathrm{diag}(R_{\mathrm{\Delta }}(0),R_{\mathrm{\Delta }}(1),\dots ,R_{\mathrm{\Delta }}(M-1))$ , respectively.

Define an error as $E_s(k)=Y_r(k)-\mathrm{\Psi }x(k)-{\mathrm{\Gamma }}_{w}w(k-1)-{\mathrm{\Gamma }}_{G}d(k-1)$ which is known and calculable at instant k, equation (22) can be transformed into

$$minJ\left(k\right)=\frac{1}{2}\mathrm{\Delta }w{\left(k\right)}^{T}H\mathrm{\Delta }w\left(k\right)+f^T\mathrm{\Delta }w\left(k\right)+F$$

(23)

where $H=2({\mathrm{\Theta }}^{T}Q\mathrm{\Theta }+R_{\mathrm{\Delta }})\mathrm{and}f=-2{\mathrm{\Theta }}^{T}Q{E}_s(k)$ . $F=E_s^T(k)Q{E}_s(k)$ , which does not depend on $\mathrm{\Delta }w(k)$ . Hence, equation (23) is equivalent to

$$minV\left(k\right)=\frac{1}{2}\mathrm{\Delta }w{\left(k\right)}^{T}H\mathrm{\Delta }w\left(k\right)+f^T\mathrm{\Delta }w\left(k\right)$$

(24)

Optimization process

The necessary condition satisfying equation (24) is

$${\nabla }_{\mathrm{\Delta }w\left(k\right)}V=H\mathrm{\Delta }w\left(k\right)+f=0$$

(25)

where ${\nabla }_{\mathrm{\Delta }w(k)}$ denotes the gradients related to $\mathrm{\Delta }w(k)$ . It follows that

$$\mathrm{\Delta }w{\left(k\right)}_{\mathrm{opt}}=-H^{-1}f$$

(26)

where $\mathrm{\Delta }w(k{)}_{\mathrm{opt}}$ denotes the optimized $\mathrm{\Delta }w(k)$ , $H\in R^{M\times M}$ is the constant reversible matrix, and $f\in R^M$ is the constant vector. The optimization procedure is represented as follows:

• Step 1. Initialize M, P, Q, R, Y_r(k), A, B, C, and G.

• Step 2. k = 0. Initialize x(k), w(k − 1), f_d(k − 1), and Δw(k).

• Step 3. Calculate $E_s(k),H,f,\mathrm{and}\mathrm{\Delta }{w}_p(k)$ .

• Step 4. $w(k)=w(k-1)+\mathrm{\Delta }\hat{w}(k)$ . Output w(k).

• Step 5. k =k + 1. If k is less than a specified instant, then go to Step 3. Otherwise, stop.

Simulation results

The MPC parameters are configured as P = 30; M = 4; Q = I_P × P, where I_P × P denotes P × P unit matrix; R = e⁻¹⁰I_P × P, which implies slightly weighting the inputs; Y_r(k) = 0.01I_P × P, that is, the control objective is to maintain the gap at 1 cm. The initial gap is configured as $\delta =1.6\mathrm{cm}$ . At 0.5 s, the 1000 N downward disturbance force is imposed on the system, and it is canceled at 1 s. The air gap is controlled toward 1 cm.

The simulation results are displayed in Figure 5. Comparing Figure 3 with Figure 5, we can learn that the MPC controller has faster dynamic performances. During 0.08–0.5 s, the air gap and other state variables have approached their steady-state values. The promptness of the MPC controller has improved 62.8% compared with that of state feedback. However, the MPC controller can promptly resist the external disturbance and maintain the air gap around 1 cm during 0.5–1.5 s. Only 0.01 cm air gap fluctuation can be observed within 0.035 s. The fluctuations of speed, acceleration, current, and voltage are much smaller during 0.5–1.5 s compared with those in Figure 3. Both speed and acceleration approach zero during 0.5–1.5 s. The control current and the voltage are kept around 37.9 A and 26.5 V during 0.61–1.0 s, respectively. During 1.09–1.5 s, the control current and the voltage are maintained around 35.6 A and 24.9 V, respectively. Figure 5 manifests that the combined MPC and state feedback control approach can, on one hand, make the plant output be stable and, on the other hand, promptly track the reference input and effectively prevent the effects from external disturbances. The promptness and disturbance-resisting ability have been improved.

Figure 5.

MPC performances of a magnetic levitation system: (a) gap, (b) speed, (c) acceleration, (d) current, and (e) voltage.

Conclusion

During the operation process of a maglev train, the vertical downward disturbances often occur because of passengers’ boarding and alighting from the train. MPC is a computer-based optimization control approach, which can directly deal with uncertain disturbances and maintain favorable set point tracking performances. However, MPC cannot be directly applied to the magnetic levitation system because it is of open-loop instability. In view of this aspect, two-level state feedback is designed. One is the nonlinear state feedback to endow the magnetic levitation system with linear dynamics, and the other is the linear state feedback to improve transient responses through the closed-loop pole placement. Based upon the constructed prediction model, MPC is utilized to further improve the transient and disturbance-resisting performances. Ultimately, the control strategy of MPC incorporating two-level state feedback is developed with prompt set point control and disturbance-resisting performances. In this study, the linearized prediction model described by equation (9) has been attempted to be employed by MPC, but stable control performances cannot be achieved. Therefore, two-level state feedback is designed and incorporated into MPC to form the control structure as shown in Figure 4. The air gap, electromagnet acceleration, and coil current are measurable using available high-resolution sensors. Consequently, the proposed control strategy provides a feasible computer-based control approach for the magnetic levitation system.