基于Youla参数化的高速磁浮列车悬浮系统控制与优化设计

摘要

为了解决高速磁浮列车悬浮系统长期运行过程中面临的性能退化问题，通过悬浮系统的复杂动态特性分析，基于控制器Youla参数化形式，提出了一种即插即用的悬浮系统控制与优化模块化架构，并设计了基于残差驱动的在线优化算法。仿真结果表明，设计的控制与优化架构以及在线优化算法有效地提高了悬浮系统对未知扰动的鲁棒性和适应能力。

关键词

高速磁浮列车; 性能退化; Youla参数化; 残差; 在线优化

高速磁浮列车被称为零高度飞行器，彻底颠覆了传统高速轮轨列车的轮轨关系与弓网关系，从而摆脱了黏着限制、噪音振动以及车轨磨损，因此高速磁浮技术是引发轨道交通产业变革的前瞻性、颠覆性技术，也是世界科技争夺制高点^［1-3］。德国、日本、美国等国家均开展了相关技术研究^［4-6］，我国也将发展高速磁浮技术提升到国家战略层面。悬浮系统是时速600 km高速磁浮列车的核心关键系统，车辆能否在复杂工况下克服重力实现贴地飞行，直接取决于悬浮系统。然而，悬浮系统结构复杂，呈现出非线性及强动态耦合特性，是一个多输入多输出、时变的不稳定系统，因此悬浮系统的控制与优化设计面临巨大的挑战。针对磁浮列车悬浮系统，比例-积分-微分（PID）算法不依赖系统模型、实现形式简单，但是依赖调试经验^［6-7］；滑模控制等先进控制方法能够提升系统的适应能力、提高电流的响应速度，但是需要精确的模型参数^［8-10］，对控制器的计算能力和存储能力也要求较高。针对高速磁浮列车悬浮系统的特性和复杂工况，很多学者开展了大量研究工作，取得了很好的控制效果^［11-13］。然而，悬浮系统长期运行过程中会面临运行环境的不确定性以及外界扰动的强随机性、不可预测性，导致已经设计并调试良好的控制回路也会出现性能退化问题。

为了提高悬浮控制系统对未知扰动的鲁棒性和适应能力，解决长期运行过程中面临的性能退化问题，重点对悬浮系统的控制与优化问题展开研究。通过悬浮系统的复杂动态特性分析，基于控制器Youla参数化形式，提出一种即插即用的悬浮系统控制与优化模块化架构，并设计了基于残差驱动的在线优化算法。通过对悬浮控制系统的动态补偿与在线优化，实现高速磁浮列车悬浮系统的高效韧性运行。

1 高速磁浮列车悬浮系统结构

高速磁浮列车车辆总体结构分为上下2层，以车厢为主的上层部分与高铁等传统轮轨列车基本相同，而以悬浮架以及悬浮单元为主的下层结构是高速磁浮列车的特有部分。一节高速磁浮列车由4个悬浮架构成，每个悬浮架包含4个悬浮搭接结构，每个悬浮搭接结构由2个悬浮单元构成。悬浮单元是高速磁浮列车的“轮子”，支撑起整个列车车体以实现贴地飞行。每个悬浮单元为一套完整的悬浮控制系统，由悬浮电磁铁、悬浮控制器以及悬浮传感器构成。高速磁浮列车的车辆总体结构如图1所示。

图1  高速磁浮列车结构

Fig.1  Structure of high-speed maglev train

高速磁浮列车悬浮系统（见图2）采用一种专门设计的搭接结构，悬浮搭接结构将左右两侧的悬浮电磁铁通过2个托臂以及托臂连接件耦合在一起。托臂与悬浮电磁铁之间安装有叠片弹簧。当列车收到起浮指令向上运动时，悬浮电磁铁压缩叠片弹簧从而产生向上的支撑力，支撑起整个列车车体。

图2  高速磁浮列车悬浮系统结构

Fig.2  Structure of magnetic levitation system in high-speed maglev train

2 悬浮系统建模

为了分析高速磁浮列车悬浮系统的动力学特征与运行规律，需要建立相应的数学模型。然而，悬浮系统结构复杂，呈现出非线性及强动态耦合特性，是一个多输入多输出、时变的不稳定系统，因此描述悬浮系统变化规律的微分方程阶次非常高，巨量的计算不利于解析求解。为了降低求解难度，分析和提取悬浮系统蕴含的主要特性，一些不重要的因素在实际建模过程中被忽略。其中，漏磁和边缘效应被忽略，磁势被认为全部均匀降落在悬浮间隙上。同时，忽略轨道的弹性形变和悬浮电磁铁沿水平方向的滚动，只考虑悬浮电磁铁在垂直方向上的俯仰和升降运动。高速磁浮列车悬浮系统模型示意图如图3所示。

图3  悬浮系统模型示意图

Fig.3  Schematic diagram of magnetic levitation system model

根据电磁学理论，悬浮电磁铁和轨道之间产生电磁力，电磁力的解析表达式为

$$F_{\mathrm{e}\mathrm{l}}(i_{\mathrm{l}},c_{\mathrm{l}})=\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{l}}}{c_{\mathrm{l}}}\right)}^{\mathrm{2}}$$

（1）

$$F_{\mathrm{e}\mathrm{r}}(i_{\mathrm{r}},c_{\mathrm{r}})=\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{r}}}{c_{\mathrm{r}}}\right)}^{\mathrm{2}}$$

（2）

式中：$F_{\mathrm{e}\mathrm{l}}(i_{\mathrm{l}},c_{\mathrm{l}})$为左侧悬浮电磁铁产生的电磁力；$F_{\mathrm{e}\mathrm{r}}(i_{\mathrm{r}},c_{\mathrm{r}})$为右侧悬浮电磁铁产生的电磁力；$i_{\mathrm{l}}(t)$为左侧悬浮电磁铁中的电流，$i_{\mathrm{r}}(t)$为右侧悬浮电磁铁中的电流；$c_{\mathrm{l}}(t)$为左侧悬浮电磁铁在垂直方向上的位移，$c_{\mathrm{r}}(t)$为右侧悬浮电磁铁在垂直方向上的位移；μ₀为真空磁导率，N为线圈匝数，A为悬浮电磁铁的有效磁极面积。由于悬浮电磁铁是一个电感性元件，因此其电压与电流的关系式为

$$u_{\mathrm{l}}(t)=R{i}_{\mathrm{l}}(t)+\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{N_{}^{\mathrm{2}}{i}_{\mathrm{l}}(t)}{\mathrm{2}{c}_{\mathrm{l}}(t)/{\mu }_{\mathrm{0}}A}\right)$$

（3）

$$u_{\mathrm{r}}(t)=R{i}_{\mathrm{r}}(t)+\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{N_{}^{\mathrm{2}}{i}_{\mathrm{r}}(t)}{\mathrm{2}{c}_{\mathrm{r}}(t)/{\mu }_{\mathrm{0}}A}\right)$$

（4）

式中：$u_{\mathrm{l}}(t)$为施加到左侧悬浮电磁铁两端的控制电压；$u_{\mathrm{r}}(t)$为施加到右侧悬浮电磁铁两端的控制电压；$R$为悬浮电磁铁电阻。

悬浮电磁铁与托臂之间的叠片弹簧等效为弹簧阻尼系统，叠片弹簧作用在托臂上的力为

$$F_{\mathrm{s}\mathrm{l}}=k_{\mathrm{l}}(l_{\mathrm{l}\mathrm{0}}-(c_{\mathrm{l}}+h_{\mathrm{l}}-H))-{\eta }_{\mathrm{l}}({\dot{c}}_{\mathrm{l}}(t)-\dot{H}(t))$$

（5）

$$F_{\mathrm{s}\mathrm{r}}=k_{\mathrm{r}}(l_{\mathrm{r}\mathrm{0}}-(c_{\mathrm{r}}+h_{\mathrm{r}}-H))-{\eta }_{\mathrm{r}}({\dot{c}}_{\mathrm{r}}(t)-\dot{H}(t))$$

（6）

式中：$F_{\mathrm{s}\mathrm{l}}$表示左侧悬浮电磁铁在托臂上的作用力，$F_{\mathrm{r}\mathrm{l}}$表示右侧悬浮电磁铁在托臂上的作用力；$h_{\mathrm{l}}$表示左侧悬浮电磁铁高度，$h_{\mathrm{r}}$表示右侧悬浮电磁铁高度，$H$表示托臂在垂直方向上的位移；$k_{\mathrm{l}}$表示左侧叠片弹簧的刚度，${\eta }_{\mathrm{l}}$表示左侧叠片弹簧的阻尼，$l_{\mathrm{l}\mathrm{0}}$表示左侧叠片弹簧的固有长度；$k_{\mathrm{r}}$表示右侧叠片弹簧的刚度，${\eta }_{\mathrm{r}}$表示右侧叠片弹簧的阻尼，$l_{\mathrm{r}\mathrm{0}}$表示右侧叠片弹簧的固有长度。

根据牛顿第二定律，描述高速磁浮列车悬浮系统运动的微分方程式为

$$m_{\mathrm{l}}{\ddot{c}}_{\mathrm{l}}=m_{\mathrm{l}}g+F_{\mathrm{s}\mathrm{l}}+F_{\mathrm{d}\mathrm{l}}-F_{\mathrm{e}\mathrm{l}}$$

（7）

$$m_{\mathrm{r}}{\ddot{c}}_{\mathrm{r}}=m_{\mathrm{r}}g+F_{\mathrm{s}\mathrm{r}}+F_{\mathrm{d}\mathrm{r}}-F_{\mathrm{e}\mathrm{r}}$$

（8）

$$m_{\mathrm{b}}\ddot{H}=m_{\mathrm{b}}g+F_{\mathrm{d}\mathrm{b}}-F_{\mathrm{s}\mathrm{l}}-F_{\mathrm{s}\mathrm{r}}$$

（9）

式中：$m_{\mathrm{l}}$为左侧悬浮电磁铁等效质量，$m_{\mathrm{r}}$为右侧悬浮电磁铁等效质量，$m_{\mathrm{b}}$为托臂质量；$F_{\mathrm{d}\mathrm{l}}$为左侧悬浮电磁铁受到的外部扰动力，$F_{\mathrm{d}\mathrm{r}}$为右侧悬浮电磁铁受到的外部扰动力，$F_{\mathrm{d}\mathrm{b}}$为托臂受到的外部扰动力。

综合式（1）―（9），消去中间变量，高速磁浮列车悬浮系统的开环模型为

$$\left\{\begin{array}{l}{u}_{\mathrm{l}}(t)=R{i}_{\mathrm{l}}(t)+\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{N_{}^{\mathrm{2}}{i}_{\mathrm{l}}(t)}{\mathrm{2}{c}_{\mathrm{l}}(t)/{\mu }_{\mathrm{0}}A}\right)\\ u_{\mathrm{r}}(t)=R{i}_{\mathrm{r}}(t)+\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{N_{}^{\mathrm{2}}{i}_{\mathrm{r}}(t)}{\mathrm{2}{c}_{\mathrm{r}}(t)/{\mu }_{\mathrm{0}}A}\right)\\ m_{\mathrm{l}}{\ddot{c}}_{\mathrm{l}}=m_{\mathrm{l}}g+k_{\mathrm{l}}(l_{\mathrm{l}\mathrm{0}}-(c_{\mathrm{l}}+h_{\mathrm{l}}-H))-{\eta }_{\mathrm{l}}({\dot{c}}_{\mathrm{l}}(t)-\dot{H}(t))+F_{\mathrm{d}\mathrm{l}}-\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{l}}}{c_{\mathrm{l}}}\right)}^{\mathrm{2}}\\ m_{\mathrm{r}}{\ddot{c}}_{\mathrm{r}}=m_{\mathrm{r}}g+k_{\mathrm{r}}(l_{\mathrm{r}\mathrm{0}}-(c_{\mathrm{r}}+h_{\mathrm{r}}-H))-{\eta }_{\mathrm{r}}({\dot{c}}_{\mathrm{r}}(t)-\dot{H}(t))+F_{\mathrm{d}\mathrm{r}}-\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{r}}}{c_{\mathrm{r}}}\right)}^{\mathrm{2}}\\ \begin{array}{l}{m}_{\mathrm{b}}\ddot{H}=m_{\mathrm{b}}g+F_{\mathrm{d}\mathrm{b}}-k_{\mathrm{l}}(l_{\mathrm{l}\mathrm{0}}-(c_{\mathrm{l}}+h_{\mathrm{l}}-H))+{\eta }_{\mathrm{l}}({\dot{c}}_{\mathrm{l}}(t)-\dot{H}(t))-\\ k_{\mathrm{r}}(l_{\mathrm{r}\mathrm{0}}-(c_{\mathrm{r}}+h_{\mathrm{r}}-H))+{\eta }_{\mathrm{r}}({\dot{c}}_{\mathrm{r}}(t)-\dot{H}(t))\end{array}\end{array}\right.$$

（10）

高速磁浮列车悬浮系统是典型的非线性系统，悬浮系统模型中存在平方项等非线性环节，通过求解非线性微分方程组来解析系统性能将非常困难。高速磁浮列车悬浮系统工作点稳定，工作区间集中在工作点附近，符合非线性系统的局部特性，可以通过线性化悬浮系统模型，求解和分析系统的动态特性。

高速磁浮列车悬浮系统在平衡点处的边界约束条件为

$$\left\{\begin{array}{l}{F}_{\mathrm{e}\mathrm{l}}(i_{\mathrm{l}\mathrm{0}},c_{\mathrm{l}\mathrm{0}})=m_{\mathrm{l}}g+k_{\mathrm{l}}(l_{\mathrm{l}\mathrm{0}}-(c_{\mathrm{l}\mathrm{0}}+h_{\mathrm{l}}-H_{\mathrm{0}}))+F_{\mathrm{d}\mathrm{l}\mathrm{0}}=\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{l}\mathrm{0}}}{c_{\mathrm{l}\mathrm{0}}}\right)}^{\mathrm{2}}\\ F_{\mathrm{e}\mathrm{r}}(i_{\mathrm{r}\mathrm{0}},c_{\mathrm{r}\mathrm{0}})=m_{\mathrm{r}}g+k_{\mathrm{r}}(l_{\mathrm{r}\mathrm{0}}-(c_{\mathrm{r}\mathrm{0}}+h_{\mathrm{r}}-H_{\mathrm{0}}))+F_{\mathrm{d}\mathrm{r}\mathrm{0}}=\frac{{\mu }_{\mathrm{0}}{N}_{}^{\mathrm{2}}A}{\mathrm{4}}{\left(\frac{i_{\mathrm{r}\mathrm{0}}}{c_{\mathrm{r}\mathrm{0}}}\right)}^{\mathrm{2}}\\ m_{\mathrm{b}}g=k_{\mathrm{l}}(l_{\mathrm{l}\mathrm{0}}-(c_{\mathrm{l}\mathrm{0}}+h_{\mathrm{l}}-H_{\mathrm{0}}))+k_{\mathrm{r}}(l_{\mathrm{r}\mathrm{0}}-(c_{\mathrm{r}\mathrm{0}}+h_{\mathrm{r}}-H_{\mathrm{0}}))-F_{\mathrm{d}\mathrm{b}\mathrm{0}}\\ F_{\mathrm{d}\mathrm{b}\mathrm{0}}=F_{\mathrm{d}\mathrm{l}\mathrm{0}}+F_{\mathrm{d}\mathrm{r}\mathrm{0}}\end{array}\right.$$

（11）

式中：$c_{\mathrm{l}\mathrm{0}}$为左侧悬浮电磁铁在平衡点处的稳态间隙，$c_{\mathrm{r}\mathrm{0}}$为右侧悬浮电磁铁在平衡点处的稳态间隙；$i_{\mathrm{l}\mathrm{0}}$为左侧悬浮电磁铁在平衡点处的稳态电流，$i_{\mathrm{r}\mathrm{0}}$为右侧悬浮电磁铁在平衡点处的稳态电流；$F_{\mathrm{d}\mathrm{l}\mathrm{0}}$为左侧悬浮电磁铁在平衡点处的静态扰动力，$F_{\mathrm{d}\mathrm{r}\mathrm{0}}$为右侧悬浮电磁铁在平衡点处的静态扰动力，$F_{\mathrm{d}\mathrm{b}\mathrm{0}}$为托臂在平衡点处的静态扰动力；$H_{\mathrm{0}}$为悬浮系统工作在平衡点处时托臂与轨道间垂向距离。

高速悬浮系统模型在平衡工作点处进行泰勒级数展开，如下所示：

$$\left\{\begin{array}{l}\mathrm{\Delta }{u}_{\mathrm{l}}(t)=R\mathrm{\Delta }{i}_{\mathrm{l}}(t)+\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}}\right)\mathrm{\Delta }{\dot{i}}_{\mathrm{l}}(t)-\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{l}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}^{\mathrm{2}}}\right)\mathrm{\Delta }{\dot{c}}_{\mathrm{l}}(t)\\ \mathrm{\Delta }{u}_{\mathrm{r}}(t)=R\mathrm{\Delta }{i}_{\mathrm{r}}(t)+\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}}\right)\mathrm{\Delta }{\dot{i}}_{\mathrm{r}}(t)-\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{r}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}^{\mathrm{2}}}\right)\mathrm{\Delta }{\dot{c}}_{\mathrm{r}}(t)\\ m_{\mathrm{l}}\mathrm{\Delta }{\ddot{c}}_{\mathrm{l}}=k_{\mathrm{l}}(\mathrm{\Delta }H-\mathrm{\Delta }{c}_{\mathrm{l}}(t))-{\eta }_{\mathrm{l}}(\mathrm{\Delta }{\dot{c}}_{\mathrm{l}}(t)-\mathrm{\Delta }\dot{H}(t))-\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{l}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}^{\mathrm{2}}}\right)\mathrm{\Delta }{i}_{\mathrm{l}}(t)+\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{l}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}^{\mathrm{3}}}\right)\mathrm{\Delta }{c}_{\mathrm{l}}(t)\\ m_{\mathrm{r}}\mathrm{\Delta }{\ddot{c}}_{\mathrm{r}}=k_{\mathrm{r}}(\mathrm{\Delta }H-\mathrm{\Delta }{c}_{\mathrm{r}}(t))-{\eta }_{\mathrm{r}}(\mathrm{\Delta }{\dot{c}}_{\mathrm{r}}(t)-\mathrm{\Delta }\dot{H}(t))-\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{r}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}^{\mathrm{2}}}\right)\mathrm{\Delta }{i}_{\mathrm{r}}(t)+\left(\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{r}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}^{\mathrm{3}}}\right)\mathrm{\Delta }{c}_{\mathrm{r}}(t)\\ m_{\mathrm{b}}\mathrm{\Delta }\ddot{H}={\eta }_{\mathrm{l}}(\mathrm{\Delta }{\dot{c}}_{\mathrm{l}}(t)-\mathrm{\Delta }\dot{H}(t))-k_{\mathrm{l}}(\mathrm{\Delta }H-\mathrm{\Delta }{c}_{\mathrm{l}}(t))+{\eta }_{\mathrm{r}}(\mathrm{\Delta }{\dot{c}}_{\mathrm{r}}(t)-\mathrm{\Delta }\dot{H}(t))-k_{\mathrm{r}}(\mathrm{\Delta }H-\mathrm{\Delta }{c}_{\mathrm{r}}(t))\end{array}\right.$$

（12）

由上述分析和讨论可知，线性化后高速磁浮列车悬浮系统模型的阶次为8阶，状态变量如下所示：

$$\mathit{x}=(x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}},x_{\mathrm{5}},x_{\mathrm{6}},x_{\mathrm{7}},x_{\mathrm{8}}{)}^{\mathrm{T}}=(\mathrm{\Delta }{c}_{\mathrm{l}},\mathrm{\Delta }{\dot{c}}_{\mathrm{l}},\mathrm{\Delta }{i}_{\mathrm{l}},\mathrm{\Delta }{c}_{\mathrm{r}},\mathrm{\Delta }{\dot{c}}_{\mathrm{r}},\mathrm{\Delta }{i}_{\mathrm{r}}{,\mathrm{\Delta }H,\mathrm{\Delta }\dot{H})}^{\mathrm{T}}$$

悬浮系统的输入为

$$\mathit{u}=[\begin{array}{cc}{u}_{\mathrm{1}}& u_{\mathrm{2}}\end{array}]=[\begin{array}{cc}\mathrm{\Delta }{u}_{\mathrm{l}}& \mathrm{\Delta }{u}_{\mathrm{r}}\end{array}]$$

悬浮系统的状态方程为

$$\left\{\begin{array}{l}{\dot{x}}_{\mathrm{1}}=x_{\mathrm{2}}\\ {\dot{x}}_{\mathrm{2}}=\frac{p_{\mathrm{5}}-k_{\mathrm{l}}}{m_{\mathrm{l}}}{x}_{\mathrm{1}}+\frac{-{\eta }_{\mathrm{l}}}{m_{\mathrm{l}}}{x}_{\mathrm{2}}+\frac{-p_{\mathrm{2}}}{m_{\mathrm{l}}}{x}_{\mathrm{3}}+\frac{k_{\mathrm{l}}}{m_{\mathrm{l}}}{x}_{\mathrm{7}}+\frac{{\eta }_{\mathrm{l}}}{m_{\mathrm{l}}}{x}_{\mathrm{8}}\\ {\dot{x}}_{\mathrm{3}}=\frac{p_{\mathrm{2}}}{p_{\mathrm{1}}}{x}_{\mathrm{2}}+\frac{-R}{p_{\mathrm{1}}}{x}_{\mathrm{3}}+\frac{\mathrm{1}}{p_{\mathrm{1}}}{u}_{\mathrm{1}}\\ {\dot{x}}_{\mathrm{4}}=x_{\mathrm{5}}\\ {\dot{x}}_{\mathrm{5}}=\frac{p_{\mathrm{6}}-k_{\mathrm{r}}}{m_{\mathrm{r}}}{x}_{\mathrm{4}}+\frac{-{\eta }_{\mathrm{r}}}{m_{\mathrm{r}}}{x}_{\mathrm{5}}+\frac{-p_{\mathrm{5}}}{m_{\mathrm{r}}}{x}_{\mathrm{6}}+\frac{k_{\mathrm{r}}}{m_{\mathrm{r}}}{x}_{\mathrm{7}}+\frac{{\eta }_{\mathrm{r}}}{m_{\mathrm{r}}}{x}_{\mathrm{8}}\\ {\dot{x}}_{\mathrm{6}}=\frac{p_{\mathrm{4}}}{p_{\mathrm{3}}}{x}_{\mathrm{5}}+\frac{-R}{p_{\mathrm{3}}}{x}_{\mathrm{6}}+\frac{\mathrm{1}}{p_{\mathrm{3}}}{u}_{\mathrm{2}}\\ {\dot{x}}_{\mathrm{7}}=x_{\mathrm{8}}\\ {\dot{x}}_{\mathrm{8}}=\frac{k_{\mathrm{l}}}{m_{\mathrm{b}}}{x}_{\mathrm{1}}+\frac{{\eta }_{\mathrm{l}}}{m_{\mathrm{b}}}{x}_{\mathrm{2}}+\frac{k_{\mathrm{r}}}{m_{\mathrm{b}}}{x}_{\mathrm{4}}+\frac{{\eta }_{\mathrm{r}}}{m_{\mathrm{b}}}{x}_{\mathrm{5}}-\frac{k_{\mathrm{l}}+k_{\mathrm{r}}}{m_{\mathrm{b}}}{x}_{\mathrm{7}}-\frac{{\eta }_{\mathrm{l}}+{\eta }_{\mathrm{r}}}{m_{\mathrm{b}}}{x}_{\mathrm{8}}\end{array}\right.$$

（13）

悬浮系统的输出方程为

$$\mathit{y}=(y_{\mathrm{1}},y_{\mathrm{2}},y_{\mathrm{3}}{)}^{\mathrm{T}}=(\mathrm{\Delta }{c}_{\mathrm{l}},\mathrm{\Delta }{c}_{\mathrm{r}}{,\mathrm{\Delta }H)}^{\mathrm{T}}$$

高速磁浮列车悬浮系统的状态空间表达式为

$$\left\{\begin{array}{l}\dot{\mathit{x}}=\mathit{A}\mathit{x}+\mathit{B}\mathit{u}\\ \mathit{y}=\mathit{C}\mathit{x}\end{array}\right.$$

（14）

$$\dot{\mathit{x}}=\left(\begin{array}{cccccccc}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \frac{p_{\mathrm{5}}-k_{\mathrm{l}}}{m_{\mathrm{l}}}& \frac{-{\eta }_{\mathrm{l}}}{m_{\mathrm{l}}}& \frac{-p_{\mathrm{2}}}{m_{\mathrm{l}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{k_{\mathrm{l}}}{m_{\mathrm{l}}}& \frac{{\eta }_{\mathrm{l}}}{m_{\mathrm{l}}}\\ \mathrm{0}& \frac{p_{\mathrm{2}}}{p_{\mathrm{1}}}& \frac{-R}{p_{\mathrm{1}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{p_{\mathrm{6}}-k_{\mathrm{r}}}{m_{\mathrm{r}}}& \frac{-{\eta }_{\mathrm{r}}}{m_{\mathrm{r}}}& \frac{-p_{\mathrm{5}}}{m_{\mathrm{r}}}& \frac{k_{\mathrm{r}}}{m_{\mathrm{r}}}& \frac{{\eta }_{\mathrm{r}}}{m_{\mathrm{r}}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{p_{\mathrm{4}}}{p_{\mathrm{3}}}& \frac{-R}{p_{\mathrm{3}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \frac{k_{\mathrm{l}}}{m_{\mathrm{b}}}& \frac{{\eta }_{\mathrm{l}}}{m_{\mathrm{b}}}& \mathrm{0}& \frac{k_{\mathrm{r}}}{m_{\mathrm{b}}}& \frac{{\eta }_{\mathrm{r}}}{m_{\mathrm{b}}}& \mathrm{0}& -\frac{k_{\mathrm{l}}+k_{\mathrm{r}}}{m_{\mathrm{b}}}& -\frac{{\eta }_{\mathrm{l}}+{\eta }_{\mathrm{r}}}{m_{\mathrm{b}}}\end{array}\right)\mathit{x}+\left(\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \frac{\mathrm{1}}{p_{\mathrm{1}}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \frac{\mathrm{1}}{p_{\mathrm{3}}}\\ \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{u}$$

（15）

$$\mathit{y}=\left[\begin{array}{llllllll}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]\mathit{x}$$

（16）

$$\begin{array}{c}{p}_{\mathrm{1}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}},p_{\mathrm{2}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{l}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}^{\mathrm{2}}}\\ p_{\mathrm{3}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}},p_{\mathrm{4}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{r}\mathrm{0}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}^{\mathrm{2}}}\\ p_{\mathrm{5}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{l}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{l}\mathrm{0}}^{\mathrm{3}}},p_{\mathrm{6}}=\frac{{\mu }_{\mathrm{0}}A{N}_{}^{\mathrm{2}}{i}_{\mathrm{r}\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{c}_{\mathrm{r}\mathrm{0}}^{\mathrm{3}}}\end{array}$$

式（14）中：$\mathit{A}$为悬浮系统状态矩阵；$\mathit{B}$为悬浮系统控制矩阵；$\mathit{C}$为悬浮系统输出矩阵。

高速磁浮列车悬浮系统建模过程中涉及的参数如表1所示。

表1  高速磁浮列车悬浮系统参数

Tab.1  Parameters of magnetic levitation system in high-speed maglev train

符号	物理意义	参考取值	单位
μ	真空磁导率	4π×10-7	H•m-1
N	电磁铁线圈匝数	270
A	电磁铁磁极面积	0.115	m2
R	电磁铁线圈电阻	3.84	Ω
$m_{\mathrm{l}}$	左侧电磁铁质量	300	kg
$m_{\mathrm{r}}$	右侧电磁铁质量	300	kg
$m_{\mathrm{b}}$	托臂等效质量	450	kg
$F_{\mathrm{d}\mathrm{b}0}$	车厢作用于托臂上的静态扰动力	24.5	kN
$c_{\mathrm{l}0}$	左侧设定间隙	12	mm
$c_{\mathrm{r}0}$	右侧设定间隙	12	mm
$i_{\mathrm{l}0}$	左侧稳态电流	31	A
$i_{\mathrm{r}0}$	右侧稳态电流	31	A

3 基于Youla参数化的控制与优化方案设计

3.1　镇定控制器及其Youla参数化形式

对于可控可观的线性时不变离散系统$\mathit{G}(z)$，其离散状态空间表达式的最小实现形式^［14］为

$$\left\{\begin{array}{l}\mathit{x}(k+\mathrm{1})=\mathit{A}\mathit{x}(k)+\mathit{B}\mathit{u}(k)\\ \mathit{y}(k)=\mathit{C}\mathit{x}(k)+\mathit{D}\mathit{u}(k)\end{array}\right.$$

（17）

式中：$\mathit{D}$为直接传递矩阵。对于任意反馈回路，控制输入$\mathit{u}(k)\in {\mathbf{R}}^l$，系统输出$\mathit{y}(k)\in {\mathbf{R}}^n$，系统状态变量$\mathit{x}(k)\in {\mathbf{R}}^n$。系统的传递函数为

$$\mathit{G}(z)=\frac{\mathit{y}(z)}{\mathit{u}(z)}=\mathit{C}{(z\mathit{I}-\mathit{A})}^{-\mathrm{1}}\mathit{B}+\mathit{D}$$

（18）

式中：I为单位矩阵。如果反馈控制系统是适定的，则在$\mathrm{R}{\mathrm{H}}_{\mathrm{\infty }}$域内$\mathit{G}(z)$均可以实现左右互质因式分解^［15］。

取${\mathit{A}}_F=\mathit{A}+\mathit{B}\mathit{F}$，其中${\mathit{A}}_F$的特征根全部限定在以原点为中心的单位圆内。重新选择新的系统控制输入变量$\mathit{\omega }(k)=\mathit{u}-\mathit{F}\mathit{x}(k)$，系统输出变量为$\mathit{u}(k)$，则

$$\left\{\begin{array}{l}\mathit{x}(k+\mathrm{1})={\mathit{A}}_F\mathit{x}(k)+\mathit{B}\mathit{\omega }(k)\\ \mathit{u}(k)=\mathit{F}\mathit{x}(k)+\mathit{\omega }(k)\end{array}\right.$$

（19）

输入$\mathit{\omega }(z)$到系统输出$\mathit{u}(z)$的传递函数为

$$\mathit{M}(z)=\left[\begin{array}{cc}\mathit{A}+\mathit{B}\mathit{F}& \mathit{B}\\ \mathit{F}& \mathit{I}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}+\mathit{B}\mathit{F}& \mathit{B}& \mathit{F}& \mathit{I}\end{array}\right]$$

（20）

取${\mathit{C}}_F=\mathit{C}+\mathit{D}\mathit{F}$，其中${\mathit{A}}_F$的特征根全部限定在以原点为中心的单位圆内。重新选择新的系统控制输入变量$\mathit{\omega }(k)=\mathit{u}-\mathit{F}\mathit{x}(k)$，系统输出变量为$\mathit{y}(k)$，则

$$\left\{\begin{array}{l}\mathit{x}(k+\mathrm{1})={\mathit{A}}_F\mathit{x}(k)+\mathit{B}\mathit{\omega }(k)\\ \mathit{y}(k)={\mathit{C}}_F\mathit{x}(k)+\mathit{D}\mathit{\omega }(k)\end{array}\right.$$

（21）

输入$\mathit{\omega }(z)$到系统输出$\mathit{y}(z)$的传递函数为

$$\mathit{N}(z)=\left[\begin{array}{cc}\mathit{A}+\mathit{B}\mathit{F}& \mathit{B}\\ \mathit{C}+\mathit{D}\mathit{F}& \mathit{D}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}+\mathit{B}\mathit{F}& \mathit{B}& \mathit{C}+\mathit{D}\mathit{F}& \mathit{D}\end{array}\right]$$

（22）

$\mathit{M}(z)$可逆的条件下，则满足

$$\mathit{G}(z)=\mathit{N}(z){\mathit{M}}^{-\mathrm{1}}(z)$$

（23）

此时，被控对象$\mathit{G}(z)$成功实现了右互质分解。依据对偶原理，$\mathit{G}(z)$的对偶系统可以表示为${\mathit{G}}^{\mathrm{T}}(z)=\left[\begin{array}{cccc}{\mathit{A}}^{\mathrm{T}}& {\mathit{C}}^{\mathrm{T}}& {\mathit{B}}^{\mathrm{T}}& {\mathit{D}}^{\mathrm{T}}\end{array}\right]$，定义${\mathit{A}}_L=\mathit{A}-\mathit{L}\mathit{C}$，${\mathit{B}}_L=\mathit{B}-\mathit{L}\mathit{D}$，其中${\mathit{A}}_L$的特征根全部限定在以原点为中心的单位圆内。同理，定义矩阵$\widehat{\mathit{M}}(z)$和$\widehat{\mathit{N}}(z)$满足下列关系：

$$\widehat{\mathit{M}}\left(z\right)=\left[\begin{array}{cc}\mathit{A}-\mathit{L}\mathit{C}& -\mathit{L}\\ \mathit{C}& \mathit{I}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}-\mathit{L}\mathit{C}& -\mathit{L}& \mathit{C}& \mathit{I}\end{array}\right]$$

（24）

$$\widehat{\mathit{N}}\left(z\right)=\left[\begin{array}{cc}\mathit{A}-\mathit{L}\mathit{C}& \mathit{B}-\mathit{L}\mathit{D}\\ \mathit{C}& \mathit{D}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}-\mathit{L}\mathit{C}& \mathit{B}-\mathit{L}\mathit{D}& \mathit{C}& \mathit{D}\end{array}\right]$$

（25）

$\widehat{\mathit{M}}\left(z\right)$可逆的条件下，则满足

$$\mathit{G}(z)={\widehat{\mathit{M}}}^{-\mathrm{1}}(z)\widehat{\mathit{N}}(z)$$

（26）

此时，被控对象$\mathit{G}(z)$成功实现了左互质分解。对于被控对象传递函数$\mathit{G}(z)$，存在8个稳定的传递函数矩阵$\mathit{M}(z)$、$\mathit{N}(z)$、$\mathit{X}(z)$、$\mathit{Y}(z)$、$\widehat{\mathit{M}}(z)$、$\widehat{\mathit{N}}(z)$、$\widehat{\mathit{X}}(z)$以及$\widehat{\mathit{Y}}(z)$满足Bezout等式^［16］：

$$\left[\begin{array}{cc}\mathit{X}(z)& \mathit{Y}(z)\\ -\widehat{\mathit{N}}(z)& \widehat{\mathit{M}}(z)\end{array}\right]\left[\begin{array}{cc}\mathit{M}(z)& -\widehat{\mathit{Y}}(z)\\ \mathit{N}(z)& \widehat{\mathit{X}}(z)\end{array}\right]=\left[\begin{array}{cc}\mathit{I}& \mathrm{0}\\ \mathrm{0}& \mathit{I}\end{array}\right]$$

（27）

所有镇定控制器$\mathit{K}(z)$都可以通过左右互质分解的形式实现Youla参数化，分别为

$$\mathit{K}(z)=-(\mathit{X}(z)-\mathit{Q}(z)\widehat{\mathit{N}}{(z))}^{-\mathrm{1}}(\mathit{Y}(z)+\mathit{Q}(z)\widehat{\mathit{M}}(z))$$

（28）

$$\mathit{K}(z)=-(\widehat{\mathit{Y}}(z)+\mathit{M}(z)\mathit{Q}(z))(\widehat{\mathit{X}}(z)-\mathit{N}(z)\mathit{Q}{(z))}^{-\mathrm{1}}$$

（29）

式中：$\mathit{Q}(z)\in \mathrm{R}{\mathrm{H}}_{\mathrm{\infty }}$是任意稳定的传递函数矩阵，称之为Youla参数矩阵。

3.2　基于Youla参数化的控制与优化模块化架构

高速磁浮列车悬浮控制系统采用标准的反馈控制系统结构，其动态结构如图4所示。图4中：$\mathit{K}(z)$为设计的悬浮控制器；$\mathit{G}(z)$为被控对象悬浮电磁铁；$\mathit{\omega }$为参考输入，即悬浮设定间隙；$\mathit{e}$为悬浮控制系统的跟踪误差；$\mathit{u}$为控制输入，即悬浮电磁铁上的电压；$\mathit{y}$为系统输出，即测量得到的悬浮间隙。

图4  悬浮系统标准反馈控制结构

Fig.4  Standard feedback control structure of magnetic levitation system

悬浮控制器经过Youla参数化之后，控制输入$\mathit{u}$可以表示为

$$\mathit{u}=(\mathit{X}(z)-\mathit{Q}(z)\widehat{\mathit{N}}{(z))}^{-\mathrm{1}}(\mathit{Y}(z)+\mathit{Q}(z)\widehat{\mathit{M}}(z))\mathit{e}$$

（30）

经过Youla参数化后的控制器阶次较高，各个环节并没有实际的物理意义，因此在实际应用过程中缺乏可行性，并不能直接用于悬浮控制器设计。

对于可控可观的线性时不变离散系统$\mathit{G}(z)$，其状态观测器的空间表达式为

$$\left\{\begin{array}{l}\widehat{\mathit{x}}(k+\mathrm{1})=\mathit{A}\widehat{\mathit{x}}(k)+\mathit{B}\mathit{u}(k)+\mathit{L}(\mathit{y}(k)-\widehat{\mathit{y}}(k))\\ \widehat{\mathit{y}}(k)=\mathit{C}\widehat{\mathit{x}}(k)+\mathit{D}\mathit{u}(k)\end{array}\right.$$

（31）

式中：$\widehat{\mathit{x}}(k)$为系统状态量的观测量。定义残差$\mathit{r}(k)=\mathit{y}(k)-\widehat{\mathit{y}}(k)$，$\mathit{r}(k)$是实际系统的输出值与估计值之间的差值，则可以得到

$$\begin{array}{l}\widehat{\mathit{x}}(k+\mathrm{1})=(\mathit{A}-\mathit{L}\mathit{C})\widehat{\mathit{x}}(k)+\\ (\mathit{B}-\mathit{L}\mathit{D})\mathit{u}(k)+\mathit{L}\mathit{y}(k)\end{array}$$

（32）

分别对式（31）以及式（32）进行$z$变换，则可以得到

$$\left\{\begin{array}{l}\widehat{\mathit{x}}(z)={(z\mathit{I}-\mathit{A}+\mathit{L}\mathit{C})}^{-\mathrm{1}}((\mathit{B}-\mathit{L}\mathit{D})\mathit{u}(z)+\mathit{L}\mathit{y}(z))\\ \widehat{\mathit{y}}(z)=\mathit{C}\widehat{\mathit{x}}(z)+\mathit{D}\mathit{u}(z)\end{array}\right.$$

（33）

$$\left\{\begin{array}{l}\widehat{\mathit{M}}\left(z\right)=\left[\begin{array}{cc}\mathit{A}-\mathit{L}\mathit{C}& -\mathit{L}\\ \mathit{C}& \mathit{I}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}-\mathit{L}\mathit{C}& -\mathit{L}& \mathit{C}& \mathit{I}\end{array}\right]=\mathit{I}-\mathit{C}{(z\mathit{I}-\mathit{A}+\mathit{L}\mathit{C})}^{-\mathrm{1}}\mathit{L}\\ \widehat{\mathit{N}}\left(z\right)=\left[\begin{array}{cc}\mathit{A}-\mathit{L}\mathit{C}& \mathit{B}-\mathit{L}\mathit{D}\\ \mathit{C}& \mathit{D}\end{array}\right]=\left[\begin{array}{cccc}\mathit{A}-\mathit{L}\mathit{C}& \mathit{B}-\mathit{L}\mathit{D}& \mathit{C}& \mathit{D}\end{array}\right]=\mathit{D}+\mathit{C}{(z\mathit{I}-\mathit{A}+\mathit{L}\mathit{C})}^{-\mathrm{1}}(\mathit{B}-\mathit{L}\mathit{D})\end{array}\right.$$

（34）

系统残差可以表示为

$$\mathit{r}(z)=\widehat{\mathit{M}}(z)\mathit{y}(z)-\widehat{\mathit{N}}(z)\mathit{u}(z)$$

（35）

将式（33）―（35）代入式（30），最终得到以下表示式：

$$\mathit{u}(z)=\mathit{F}{(z\mathit{I}-\mathit{A}+\mathit{L}\mathit{C})}^{-\mathrm{1}}((\mathit{B}-\mathit{L}\mathit{D})\mathit{u}(z)+\mathit{L}\mathit{y}(z))-\mathit{Q}(z)(\widehat{\mathit{M}}(z)\mathit{y}(z)-\widehat{\mathit{N}}(z)\mathit{u}(z))=$$

$$\mathit{F}\widehat{\mathit{x}}(z)-\mathit{Q}(z)\mathit{r}(z)$$

（36）

此时，控制器的结构形式分解成2个部分，即基于观测器估计量的反馈控制部分$\mathit{F}\widehat{\mathit{x}}(z)$以及基于残差驱动的动态补偿部分$-\mathit{Q}(z)\mathit{r}(z)$。观测器估计得到的相关状态信息往往比较敏感，而且包含各种噪声，不利于实际应用。

存在反馈控制器${\mathit{K}}_{\mathrm{0}}(z)$使闭环控制系统稳定，按照式（36）所示的结构形式，反馈控制器${\mathit{K}}_{\mathrm{0}}(z)$可以表示为

$${\mathit{u}}_{\mathrm{0}}(z)=\mathit{F}\widehat{\mathit{x}}(z)-{\mathit{Q}}_{\mathrm{0}}(z)\mathit{r}(z)$$

（37）

任意镇定控制器$\mathit{K}(z)$均可以采用如式（36）所示的结构形式进行表示：

$$\mathit{u}(z)=\mathit{F}\widehat{\mathit{x}}(z)-{\mathit{Q}}_{\mathrm{c}}(z)\mathit{r}(z)$$

（38）

式中：${\mathit{Q}}_{\mathrm{c}}(z)\in \mathrm{R}{\mathrm{H}}_{\mathrm{\infty }}$为对应的Youla参数矩阵。如果令${\mathit{Q}}_{\mathrm{c}}(z)={\mathit{Q}}_{\mathrm{0}}(z)-\mathit{Q}(z)$，任意镇定控制器$\mathit{K}(z)$可以表示为

$$\mathit{u}(z)=\mathit{F}\widehat{\mathit{x}}(z)-{\mathit{Q}}_{\mathrm{0}}(z)\mathit{r}(z)+\mathit{Q}(z)\mathit{r}(z)=$$

$${\mathit{u}}_{\mathrm{0}}(z)+\mathit{Q}(z)\mathit{r}(z)={\mathit{u}}_{\mathrm{0}}(z)+{\mathit{u}}_{\mathrm{r}}(z)$$

（39）

此时，控制器的结构形式分解成2个模块：既有控制器${\mathit{K}}_{\mathrm{0}}(z)$以及基于残差驱动的动态补偿器。控制量$\mathit{u}$由两部分组成：既有控制器产生的${\mathit{u}}_{\mathrm{0}}$以及动态补偿器产生的${\mathit{u}}_{\mathrm{r}}$。因此，悬浮控制器可以进行模块化设计，既有控制器${\mathit{K}}_{\mathrm{0}}(z)$能够完全继承和保留现有的控制器设计与调试经验，可以直接利用现有控制器的主体结构与控制算法。基于残差驱动的动态补偿器能够根据悬浮系统残差的变化，动态调整控制量的输出。当悬浮控制系统工作在期望状态时，悬浮系统的残差为零，此时动态补偿器控制量的输出为零。当悬浮控制系统工作在非期望状态时，动态补偿器会对悬浮控制系统进行动态补偿，最终实现较为理想的控制状态。高速磁浮列车悬浮控制系统的动态结构如图5所示。

图5  基于Youla参数化的悬浮系统控制与优化模块化架构

Fig.5  Integrated control and optimization architecture of magnetic levitation system based on Youla parameterization

如图5所示，提出的基于Youla参数化的悬浮系统控制与优化模块化架构物理意义明确。该架构不改变既有控制器，只通过设计动态补偿器就可以实现对悬浮控制系统的动态补偿与在线优化，即插即用，实现形式简单灵活。

4 基于残差驱动的在线优化算法设计

动态补偿器中Youla参数矩阵$\mathit{Q}(z)$是悬浮系统实现在线优化的关键，如图5所示，Youla参数矩阵$\mathit{Q}(z)$直接决定动态补偿器控制量的输出。Youla参数矩阵$\mathit{Q}(z)$的输入规范型可以表示为

$$\left\{\begin{array}{l}{\mathit{x}}_{\mathrm{r}}(k+\mathrm{1})={\mathit{A}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}}){\mathit{x}}_{\mathrm{r}}(k)+{\mathit{B}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}})\mathit{r}(k)\\ {\mathit{u}}_{\mathrm{r}}(k)={\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}}}){\mathit{x}}_{\mathrm{r}}(k)\end{array}\right.$$

（40）

采用最优二次型性能指标$J$作为Youla参数矩阵$\mathit{Q}(z)$动态调整的依据，在能量消耗最小的情况下使得系统的跟踪误差最小。J的计算式如下所示：

$$J^{(j)}=$$

$$\frac{\mathrm{1}}{\mathrm{2}n}{\displaystyle \sum _{k=k_j}^{n+k_j-\mathrm{1}}}(\mathit{e}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{e}}(k)\mathit{e}(k)+\mathit{u}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{u}}(k)\mathit{u}(k))$$

（41）

式中：n表示窗口宽度；$k_j$表示窗口初值；括号中$j$表示迭代次数。权重矩阵${\mathit{W}}_{\mathit{e}}(k)$与参数向量$\mathit{\theta }$并不相关，其物理意义为误差$\mathit{e}(k)$的权重值。权重矩阵${\mathit{W}}_{\mathit{u}}(k)$与参数向量$\mathit{\theta }$也不相关，其物理意义为控制量$\mathit{u}(k)$的权重值。

Youla参数矩阵$\mathit{Q}(z)$与参数向量$\mathit{\theta }$存在一一对应的映射关系，如果参数向量$\mathit{\theta }$沿着二次型性能指标$J$对$\mathit{\theta }$的负梯度方向更新，二次型性能指标$J$就将以最快的速度达到最小。动态补偿器中Youla参数矩阵$\mathit{Q}(z)$的参数向量$\mathit{\theta }$更新律为

$${\mathit{\theta }}_{\mathit{i}}^{(\mathit{j})}={\mathit{\theta }}_{\mathit{i}}^{(\mathit{j}-\mathrm{1})}-{\lambda }_i\nabla J^{(j)}$$

（42）

式中：$i$为参数向量$\mathit{\theta }$的分量标号；${\lambda }_i$为参数更新的步长；$\nabla J^{(j)}$为性能指标$J^{(j)}$的梯度。$\nabla J^{(j)}$的计算式如下所示：

$$\nabla J^{(j)}=\frac{\partial J^{(j)}}{\partial {\mathit{\theta }}^{(j)}}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=k_j}^{n+k_j-\mathrm{1}}}\left(\mathit{e}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{e}}\left(k\right)\frac{\partial \mathit{e}(k)}{\partial {\mathit{\theta }}_i^{(j)}}+\mathit{u}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{u}}(k)\frac{\partial \mathit{u}(k)}{\partial {\mathit{\theta }}_i^{(j)}}\right)$$

（43）

求取关于参数向量${\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}}$的微分，则可以得到

$$\left\{\begin{array}{l}\frac{\partial \mathit{e}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}=-{\mathit{C}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}\\ \frac{\partial \mathit{u}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}=\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+{\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}}}^{(j-\mathrm{1})})\end{array}\right.$$

（44）

$$\left\{\begin{array}{l}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k+\mathrm{1})}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}={\mathit{A}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+{\mathit{B}}_{\mathrm{o}}\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+{\mathit{B}}_{\mathrm{o}}{\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}}}^{(j-\mathrm{1})})\frac{\partial {\mathit{x}}_{\mathrm{r}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}\\ \begin{array}{l}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k+\mathrm{1})}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}={\mathit{L}}_{\mathrm{c}}{\mathit{C}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+(\mathit{A}-{\mathit{L}}_{\mathrm{c}}\mathit{C}+\mathit{B}\mathit{F})\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}\\ \frac{\partial {\mathit{x}}_{\mathrm{r}}(k+\mathrm{1})}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}={\mathit{A}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}}^{(j-\mathrm{1})})\frac{\partial {\mathit{x}}_{\mathrm{r}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+\frac{\partial {\mathit{A}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}}^{(j)})}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}{\mathit{x}}_{\mathrm{r}}(k)+\frac{\partial {\mathit{B}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}}}^{(j)})}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}\mathit{r}(k)\end{array}\end{array}\right.$$

（45）

动态补偿器中Youla参数矩阵$\mathit{Q}(z)$关于参数向量${\mathit{\theta }}_{{\mathit{A}}_r{\mathit{B}}_r}$的更新律为

$${\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}={\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j-\mathrm{1})}-{\lambda }_i\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=k_j}^{n+k_j-\mathrm{1}}}\left(-{\mathit{C}}_{\mathrm{o}}\mathit{e}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{e}}(k)\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+\mathit{u}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{u}}(k)\left(\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{A}}_{\mathrm{r}}{\mathit{B}}_{\mathrm{r}},i}^{(j)}}+{\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}}}^{(j-\mathrm{1})})\right)\right)$$

（46）

求取关于参数向量${\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}}}$的微分，则可以得到

$$\left\{\begin{array}{l}\frac{\partial \mathit{e}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}=-{\mathit{C}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}\\ \frac{\partial \mathit{u}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}=\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+\frac{\partial {\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)})}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}{\mathit{x}}_{\mathrm{r}}(k)\end{array}\right.$$

（47）

$$\left\{\begin{array}{l}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k+\mathrm{1})}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}={\mathit{A}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+{\mathit{B}}_{\mathrm{o}}\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+{\mathit{B}}_{\mathrm{o}}\frac{\partial {\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)})}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}{\mathit{x}}_{\mathrm{r}}(k)+{\mathit{B}}_{\mathrm{o}}{\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j-\mathrm{1})})\frac{\partial {\mathit{x}}_{\mathrm{r}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}\\ \frac{\partial {\mathit{x}}_{\mathrm{c}}(k+\mathrm{1})}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}={\mathit{L}}_{\mathrm{c}}{\mathit{C}}_{\mathrm{o}}\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+(\mathit{A}-{\mathit{L}}_{\mathrm{c}}\mathit{{\rm X}}+\mathit{B}\mathit{F})\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}\end{array}\right.$$

（48）

动态补偿器中Youla参数矩阵$\mathit{Q}\left(z\right)$关于参数向量${\mathit{\theta }}_{{\mathit{C}}_r}$的更新律为

$${\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}={\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}-{\lambda }_i\frac{\mathrm{1}}{n}{\displaystyle \sum _{k=k_j}^{n+k_j-\mathrm{1}}}\left(-{\mathit{C}}_{\mathrm{o}}\mathit{e}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{e}}(k)\frac{\partial {\mathit{x}}_{\mathrm{o}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+\mathit{u}{(k)}^{\mathrm{T}}{\mathit{W}}_{\mathit{u}}(k)\left(\mathit{F}\frac{\partial {\mathit{x}}_{\mathrm{c}}(k)}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}+\frac{\partial {\mathit{C}}_{\mathrm{r}}({\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)})}{\partial {\mathit{\theta }}_{{\mathit{C}}_{\mathrm{r}},i}^{(j)}}{\mathit{x}}_{\mathrm{r}}(k)\right)\right)$$

（49）

综上，式（46）以及式（49）给出了动态补偿器中Youla参数矩阵$\mathit{Q}(z)$关于参数向量$\mathit{\theta }$的全部更新律。基于残差驱动的在线优化算法流程如图6所示。首先进行初始化，其次设置参数向量$\mathit{\theta }$更新步长、误差$\mathit{e}(k)$的权重值、控制量$\mathit{u}(k)$的权重值以及窗口宽度，然后迭代计算性能指标关于参数向量的下降梯度，最后按照梯度下降的方向更新Youla参数矩阵$\mathit{Q}(z)$对应的参数向量。

图6  基于残差驱动的在线优化算法流程

Fig.6  Flow chart of online optimization algorithm for magnetic levitation system based on residual

5 仿真验证

高速磁浮列车悬浮系统长期运行过程中面临运行环境的不确定性和外界扰动的强随机性、不可预测性，高速磁浮轨道长定子下表面会出现随机不平顺。如图7所示，轨道台阶和梁端切向扭转角为轨道随机不平顺的典型代表，对于高速运行过程中的悬浮系统特别不利。为验证所提出的控制与优化模块化架构以及基于残差驱动的在线优化算法，对悬浮系统在轨道随机不平顺扰动下的在线优化情况进行仿真验证。

图7  高速磁浮轨道随机不平顺示意图

Fig.7  Schematic diagram of random track irregularity in high-speed maglev line

当高速磁浮列车以时速600 km通过幅值为1 mm的轨道台阶时，悬浮系统优化前后的悬浮间隙波动如图8所示。对比发现，优化后的悬浮间隙波动明显变小，悬浮系统性能得到改善。

图8  通过幅值1 mm轨道台阶时的悬浮间隙

Fig.8  Suspension gap when passing track step with amplitude of 1 mm

当高速磁浮列车以时速600 km通过梁端切向扭转角0.000 4 rad的轨道梁时，悬浮系统优化前后的悬浮间隙波动如图9所示。对比发现，优化后的悬浮间隙波动明显变小，悬浮系统性能得到改善。

图9  通过梁端切向扭转角为0.000 4 rad时的悬浮间隙

Fig.9  Suspension gap when tangential torsional angle of passing beam end is 0.000 4 rad

由图8和图9可知，悬浮系统受到未知扰动时，悬浮间隙都会出现一定波动，但是采用标称控制器的悬浮系统波动更为明显。仿真结果表明，随着参数向量沿着二次型性能指标的负梯度方向不断更新，基于Youla参数矩阵$\mathit{Q}(z)$的动态补偿器能够通过对既有悬浮控制器的动态补偿，实现悬浮系统性能的在线优化。

6 结语

以高速磁浮列车悬浮系统为研究对象，针对悬浮系统长期运行过程中面临的性能退化问题，开展了悬浮系统的控制与优化问题研究。通过悬浮系统的复杂动态特性分析，建立了悬浮系统的数学模型。采用左右互质分解的形式实现了控制器的Youla参数化，提出了一种即插即用的、由既有控制器和动态补偿器构成的悬浮系统控制与优化模块化架构。通过基于残差驱动的在线优化算法设计，利用动态补偿器完成了对悬浮控制系统的动态补偿与在线优化。仿真果表明，设计的控制与优化模块化架构以及在线优化算法有效地提高了悬浮系统对未知扰动的鲁棒性和适应能力。

作者贡献声明

翟明达：研究方案的制定，算法设计与仿真，论文撰写。

李晓龙：算法设计与仿真验证。

龙志强：参与研究的构思、设计。

窦峰山：论文修订。