磁悬浮轴承技术在电力行业领域的应用主要体现在新能源发电，其中风力发电和飞轮储能对磁悬浮轴承的研究应用已十分广泛。由于磁悬浮轴承具有无摩擦的优点，在风力发电中将其励磁转子与风电机组的叶轮刚性连接，可以有效降低风电机组的启动风速，减小阻力矩，提高风能利用率和发电效率[1-2]，还可以避免润滑油的使用给风电机组带来的一系列问题，维护成本小。磁悬浮技术也被应用于偏航系统，采用磁悬浮驱动替代传统的多齿轮机械驱动，提高偏航对风精度[3-4]。以磁悬浮轴承为核心的飞轮储能作为新能源电力系统的一种新型储能技术，具有无污染、效率高、寿命长的优点[5-6]。
磁悬浮轴承-转子系统的控制技术是磁悬浮轴承高速发展的一项重要依据，主流控制技术包括PID控制、滑模控制、协同控制等。然而由于磁悬浮轴承-转子系统的复杂性，某种单一的控制算法很难保持系统的控制性能，通常需要融合不同的算法，例如模糊PID控制[7-8]、粒子群算法优化PID[9]和高阶滑模控制[10-12]等。
模型预测控制（Model Predictive Control，MPC）是一种利用模型对系统未来行为进行预测，根据预测结果对指定指标求解最优化问题的一种控制方法。文献[13]基于磁悬浮轴承的线性时不变模型，在不违反输入和状态约束的情况下构造离散时间模型预测控制，在每个采样时刻求解一个有界MPC问题，计算最优控制输入。文献[14]在磁悬浮轴承MBC500试验平台的基础上，采用MPC算法设计磁悬浮轴承径向四自由度系统控制器，通过与试验、其他控制方法的对比验证了MPC下磁悬浮轴承更优越的系统性能。文献[15]以未建模的动力学形式作为磁悬浮轴承-转子系统的外界扰动，结合扰动观测器和H∞控制，有效提高了系统的抗干扰能力。考虑MPC抗外界干扰能力较弱的问题，本文引入鲁棒H∞控制方法，将Min-Max鲁棒控制的思想代入MPC的核心环节即滚动优化算法中，通过Simulink仿真分析鲁棒模型预测控制（Robust Model Predictive Control，RMPC）下磁悬浮轴承-转子系统的鲁棒性。
1 磁悬浮轴承-转子系统模型
磁悬浮轴承-转子系统如图1所示，lₐ、lᵦ分别为磁悬浮轴承a、b到转子质心的距离，l = lₐ + lᵦ；xₐ、yₐ分别为转子在轴承a处沿x、y轴的位移，xᵦ、yᵦ分别为转子在轴承b处沿x、y轴的位移；θₓ、θᵧ分别为转子绕x、y轴旋转的角度。
根据文献[16]，四自由度径向磁悬浮轴承的状态空间表达式为：
Ẋ = AX + BU （1）
其中：
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] ]
式中：X为状态向量；U为控制向量；Y为输出向量；iₓₐ、iₓᵦ、iᵧₐ、iᵧᵦ为通入的控制电流；Kₓ、Kᵧ为位移刚度系数；Kᵢₓ、Kᵢᵧ为电流刚度系数；m为转子质量；Jₓ、Jᵧ分别为转子绕x、y轴的转动惯量；ω为转子的角速度。
2 鲁棒模型预测控制
本文基于Min-Max鲁棒模型预测控制，将不确定性因素引起“最坏”情况下的系统性能指标作为最优化目标。将系统的输入、输出约束条件转化为矩阵不等式（Linear Matrix Inequality，LMI），构造LMI最优化问题，求解状态反馈控制律，最小化LMI中的最大化问题。
利用欧拉法对（1）式的四自由度径向磁悬浮轴承的状态空间表达式离散化，采样时间Tₛ = 0.001 s。将外部干扰作为输入变量加入到系统模型中，表达式为：
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）
其中：A₁ = I + TₛA，B₁ = TₛB，C₁ = C，D₁₁ = D，A、B、C、D均为其离散化矩阵；x为系统状态；u为控制电流输入；w为所有形式的外部干扰输入；y为系统约束输出；z为控制性能输出；A₁、B₁、B₂、C₁、C₂、D₁₁、D₁₂、D₂₁、D₂₂均为参数矩阵（根据实际情况确定）。
在k时刻的预测时域N内，系统的控制输入约束、输出约束和状态约束为：
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）
由于一般情况下外部干扰输入无法直接影响到系统约束输出，所以令D₁₂ = 0。鲁棒H∞控制的设计思想是在不违反系统输出约束的前提下满足系统内部的闭环反馈稳定，优化目标是在不违反系统输出约束的前提下满足w到z的H∞范数最小，即小于正数γ：
||T_zw||∞ < γ; γ > 0 （4）
设状态反馈控制律为u = Kx，代入（2）式整理得：
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）
其中：A_cl = A₁ + B₁K，C₁_cl = C₁ + D₁₁K，C₂_cl = C₂ + D₂₁K，则存在正定对称矩阵P > 0，满足（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）
根据Schur补公式，（6）式可等价于（7）式（此处省略矩阵展开形式，详见原文推导）。对（7）式的矩阵分别左乘[x(k)ᵀ w(k)ᵀ]和右乘[x(k) w(k)]ᵀ，整理得：
x(k+1)ᵀP x(k+1) + z(k)ᵀz(k) ≤ x(k)ᵀP x(k) + γ² w(k)ᵀw(k) （8）
利用P构造李雅普诺夫函数V(x) = xᵀPx，得到耗散不等式：
V(x(k+1)) + Σₖⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σₖⁿ⁻¹ w(i)ᵀw(i) + V(x(k)) （9）
假设系统在整个预测时域内都能够满足（9）式，则有：
V(x(n)) + Σ₀ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σ₀ⁿ⁻¹ w(i)ᵀw(i) + V(x(0)) （10）
可得到（11）式的耗散不等式成立：
V(x(k+1)) + Σₖᵏ⁺ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σₖᵏ⁺ⁿ⁻¹ w(i)ᵀw(i) + V(x(k)) （11）
由于x(0) = 0，且V(x) ≥ 0，根据（11）式得：
Σ₀ⁿ⁻¹ z(i)ᵀz(i) ≤ γ² Σ₀ⁿ⁻¹ w(i)ᵀw(i) （12）
（12）式证明从w到z的H∞范数小于γ。
令Q = P⁻¹，用块对角矩阵diag{Q, I, Q, I}对（6）式的矩阵分别左乘和右乘，得（13）式；再令Y = KQ，代入（13）式可得（14）式（矩阵展开形式详见原文）。
假设系统外部干扰输入有界，则有：
Σₖᵏ⁺ⁿ⁻¹ w(i)ᵀw(i) ≤ n wₘₐₓ² （15）
式中：wₘₐₓ为外部干扰输入的最大值。给定α > 0，若系统状态属于一个椭圆域Ω₁，即：
Ω₁: {x ∈ Rⁿ | γ² n wₘₐₓ² + V(x) ≤ α} （16）
V(x) = xᵀPx，Ω₁由P、α和wₘₐₓ定义。系统闭环状态轨迹属于另一个椭圆域Ω₂，Ω₂由P和α定义：
Ω₂: {x ∈ Rⁿ | V(x) ≤ α} （17）
根据（3）式及Cauchy-Schwarz不等式，控制输入约束和输出约束可写为：
||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）
式中：K为状态反馈矩阵；α为正数。
系统约束条件可表示为：
uₘₐₓ² I - α Q⁻¹ KᵀK Q⁻¹ ≥ 0 （20）
假设存在对称矩阵X，满足X对角线元素Xᵢᵢ ≤ uₘₐₓ²（i=1,2,⋯,n）；存在对称矩阵Z，满足Z对角线元素Zᵢᵢ ≤ yₘₐₓ²（i=1,2,⋯,n）。根据Schur补公式和（13）式得出系统控制输入约束和输出约束的LMI形式：
[ [X, α Yᵀ], [α Y, Q] ] ≥ 0, i=1,2,⋯,n （21）
[ [Z, (C₁Q + D₁₁Y)ᵀ], [C₁Q + D₁₁Y, Q] ] ≥ 0 （22）
在k时刻，系统状态为x(k) ∈ Ω₁，则：
γ² n wₘₐₓ² + x(k)ᵀQ⁻¹ x(k) ≤ α （23）
结合（8）、（14）、（16）、（17）式得到LMI优化问题：
min γ²
s.t. （8）、（14）、（16）、（17）式 （24）
在每个采样时刻k（k=0,1,2,⋯）求解LMI问题，若存在解(γₖ², Qₖ, Yₖ, Xₖ, Zₖ)，则k时刻的状态反馈矩阵可表示为：
K(k) = Yₖ Qₖ⁻¹ （25）
MPC的流程如图2所示。RMPC与MPC的控制结构类似，优化目的旨在将Min-Max鲁棒控制的思想代入到MPC的核心环节即滚动优化算法中。
RMPC的滚动优化算法如下：
初始化：在k=0时刻，初始化系统状态。
迭代计算优化：在k时刻，求解（24）式LMI最优化问题。若存在解，则得到状态反馈矩阵K(k)，计算系统控制输入量uₖ = K(k)x(k)并作用于系统；若LMI最优化问题不存在解，则需增加α的值后重新求解；若持续增加α的值后仍无解，则由上一个状态反馈矩阵K(k-1)取代，即K(k) = K(k-1)。
k = k+1，循环步骤2。
磁悬浮轴承-转子系统参数见表1。
| 参数 | 数值 | 单位 | 说明 |
|---|---|---|---|
| m | 2.8 | kg | 转子质量 |
| Jₓ、Jᵧ | 0.022 | kg·m² | 转子绕x、y轴的转动惯量 |
| J | 0.0017 | kg·m² | 转子绕其他轴的转动惯量（原文未明确轴方向） |
| lₐ | 144 | mm | 磁悬浮轴承a到转子质心的距离 |
| lᵦ | 144 | mm | 磁悬浮轴承b到转子质心的距离 |
| Kᵢ（Kᵢₓ、Kᵢᵧ） | 64.5107 | N·A⁻¹ | 电流刚度系数 |
| K（Kₓ、Kᵧ） | 847855.1796 | N·m⁻¹ | 位移刚度系数 |
表1 磁悬浮轴承-转子系统参数（Tab.1 Parameters of magnetic bearing-rotor system）
3 仿真分析
根据文献[16]的模型预测控制算法，对磁悬浮轴承-转子系统进行抗干扰仿真。设置参考值为[0.1; 0.1; 0.1; 0.1]（单位：mm），系统输入约束为[2.5; 2.5; 2.5; 2.5]（单位：A），系统输出约束为[0.25; 0.25; 0.25; 0.25]（单位：mm）。0.5 s时，在磁悬浮轴承a处对转子施加一个单位脉冲信号作为外部干扰，干扰时间为0.1 s，使转子在安全气隙范围内沿x、y轴分别产生一段位移。2种控制方法下采用的参数相同，MPC的仿真
MPC的系统调节时间约为0.2 s，响应曲线有10%左右的超调量。由于设置了输入约束，系统的控制输入电流不会超过2.5 A。在0.5 s施加干扰后，各方向转子偏移，最大偏移约为0.12 mm。0.15 s后转子位移趋于稳定，但系统并未完全恢复至干扰前的稳定状态，转子出现了低频振动，即转子处于微弱的锥动状态，转子在2个旋转自由度上未达到稳定，表明MPC下该系统的抗干扰能力较弱。
由于建立的磁悬浮轴承-转子系统模型并非精确的数学模型，而是适当简化后的线性模型（标称模型），MPC是基于模型的一种反馈校正控制系统，所以无法保证快速准确地消除外部干扰的影响。本文通过增加鲁棒控制，提高系统抗干扰能力，根据（2）式的系统离散线性模型，取：
B₂ = D₁₂ = D₂₂ = I （26）
系统控制输入约束为||u||∞ ≤ 2.5 A，磁悬浮轴承-转子系统正常工作时，转子与定子之间的气隙为0.35 mm，考虑安全气隙裕度，设置系统输出约束为||y||∞ ≤ 0.25 mm；α初始值设置为α₀ = 0.1。基于RMPC，再次对系统进行干扰仿真，分别在0.5、0.7 s时加入大小相同、方向相反的扰动，仿真结果如图5所示。
在干扰输入下，RMPC的系统无超调，表明该系统具有良好的动态稳定性能，在受到外部干扰后能够快速恢复至稳定状态。转子位移变化不到0.11 mm，系统能够在0.1 s内完全恢复至受干扰前的状态，且未出现低频振动现象，有效抑制了干扰对系统的影响，提高了系统的鲁棒性。比较2种控制模式下系统的输入电流，RMPC的控制电流曲线几乎完全重叠，表明RMPC的磁悬浮轴承-转子系统在受到外界扰动后具有更强的自平衡调节能力。
4 结束语
本文建立了磁悬浮轴承-转子系统的数学模型，基于鲁棒控制的思想，在LMI优化问题的框架下，将系统外部干扰输入到控制性能输出的H∞范数作为最优化指标，求解状态反馈控制律，得出能够满足H∞性能的矩阵不等式条件。结合模型预测控制，对滚动优化算法进行改进，将原滚动优化中重复进行的二次优化问题替换为线性矩阵不等式优化问题，以系统控制输入约束和输出约束作为手段提高系统抗干扰性能。