CN103246200B - 一种基于分布式模型的动车组同步跟踪控制方法 - Google Patents

Abstract

一种基于分布式模型的动车组同步跟踪控制方法，针对动力分散式动车组由若干牵引/制动控制单元相互耦合组成的结构特点，所述方法基于动车组运行过程数据和牵引/制动特性曲线，利用数据驱动建模方法，建立动车组运行过程分布式模型，采用子空间方法辨识得到模型参数。本发明提出基于分布式状态空间模型的同步跟踪控制方法，得到各单元的准确控制量，使高速动车组各单元运行速度精确跟踪目标曲线，实现高速动车组运行过程中各单元高精度同步跟踪控制,？确保动车组正点、安全、有效运行。本发明简单实用，适用于高速动车组运行过程建模和同步跟踪控制。

Description

一种基于分布式模型的动车组同步跟踪控制方法

技术领域

本发明涉及动车组运行过程各单元分布式建模与高精度同步跟踪控制方法，属动车组运行过程建模与速度同步跟踪控制领域。

背景技术

相对普通中低速内燃列车和其它交通工具，动车组具有运营速度高、燃料省、安全可靠、乘客舒适性得以改善等明显优势，在一定程度上缓解了我国交通资源紧缺的现状。列车运行控制系统是动车组安全运行的保障，研究动车组各动力单元的动态特性及其相应的控制方法具有重要意义。国内普遍采用的是动力分散式动车组，相对集中式动车组，动力分散式动车组具有加速度快、轴重轻、易于操控和安全性高等优点，已逐渐成为现代轨道交通工具的主流和未来的发展方向。随着动车组运行速度的不断提高，各个动力单元间的非线性和强耦合等特性显著加剧，如何准确辨识每个动力单元的动态特性并设计有效的控制器对实现高速动车组安全、舒适、平稳运行具有重要意义。

针对动车组的建模方法，大致可以分为集中式和分布式建模两种。集中式建模方法是把高速动车组当作单质点进行研究，较少考虑动力单元间的耦合作用。然而，随着动车组速度的提高，这种耦合效应成为影响系统控制性能改善的主要因素，采用分布式建模更能表征这种动态特性。目前运用到动车组的分布式模型通常采用基于牵引计算和运行阻力经验模型的描述方法，还无法完整刻画动车组复杂多变的动态行为。

在列车运行控制方法发展过程中，PID控制以及带参数自适应的改进PID算法，都在一定程度上取得令人满意的效果。但高速动车组的运行过程是一个复杂的多变量、非线性动态系统，PID算法难以获得理想的效果。基于模糊控制、专家系统和鲁棒自适应控制等列车运行过程智能控制算法，取得了比PID更好的控制效果，但是模糊控制方法和专家系统控制方法主要是模拟专家实际经验，当动车组运行环境发生变化或受到随机干扰时，难以在线调整；鲁棒自适应控制器是从全局稳定性的角度而设计的且忽略了系统动态性能。此外，上述控制方法主要应用在城市轨道交通等普通速度列车，较少涉及动车组运行控制方法。

发明内容

本发明的目的是，对高速动车组复杂的非线性运行过程建立有效的分布式状态空间模型描述，以各单元子空间模型输出和期望输出之间的偏差平方和为目标函数，设计动车组各单元同步跟踪控制器，实现动车组运行速度的高精度跟踪。

本发明的技术方案是：针对动力分散式动车组由若干牵引/制动控制单元相互耦合组成的结构特点，提出动车组运行过程的分布式描述方法；借助动车组牵引/制动特性曲线和实际运行数据，采用子空间模型辨识方法建立了动车组分布式状态空间模型；提出基于高速动车组分布式模型的同步跟踪控制方法，实现高速动车组各单元运行速度的高精度同步跟踪。所述技术方案具体描述为：

1、建立基于子空间辨识的高速动车组分布式模型：

高速动车组是由多节动力车和拖车固定编组在一起组成的动力分散式系统，把高速动车组分成若干个牵引/制动控制单元来进行分布式建模更能准确的分析高速动车组的受力情况和描述其运行特性。图1描述了由r个动力单元串联组成的高速动车组。其各单元运动过程动力学模型可表示为：

式中，t＝2,...,r-1,r是高速动车组动力单元数，是动力单元受到的空气阻力；是机械阻力；u_t代表各动力单元的牵引力/制动力，且u_t＝u+Δu_t；x_t代表各动力单元的位移，代表各动力单元的速度，代表各动力单元的加速度；m_t代表各动力单元的质量；代表动力单元之间的弹簧系数；b代表相邻动力单元间阻尼器的阻尼系数。c₀、c_v、c_a是阻力系数，具有高度不确定性。代表空气阻力，是速度的非线性函数，并且随着列车运行速度的增加所占的比例越来越大，其非线性特性就越明显。由式(1)可知，高速动车组第t个动力单元的非线性状态空间模型描述如下：

$\{\begin{array}{c}{x}_t\left(k+1\right)=\mathrm{\φ}\left\{x_t\left(k\right),u_t\left(k\right)\right\}+q_t\left(k\right)\\ y_t\left(k\right)=\mathrm{\θ}\left\{x_t\left(k\right),u_t\left(k\right)\right\}+d_t\left(k\right)\end{array}---\left(2\right)$

其中，为k时刻系统输入序列，为k时刻系统输出序列，为k时刻系统状态序列；和是非线性向量函数；和分别是k时刻系统不可观测的噪声序列。

然后，利用局部建模思想，在高速动车组平衡点(x_e,u_e)附近将其近似成线性系统，式(2)等效成线性离散状态空间模型，可描述为：

$\left\{\begin{array}{c}{x}_t\left(k+1\right)\≈A_t{x}_t\left(k\right)+B_t{u}_t\left(k\right)+w_t\left(k\right)\\ y_t\left(k\right)\≈C_t{x}_t\left(k\right)+D_t{u}_t\left(k\right)+v_t\left(k\right)\end{array}\right.---\left(3\right)$

其中 $A_t\∈R^{n_t\×n_t},B_t\∈R^{n_t\×l_t},C_t\∈R^{m_t\×n_t},D_t\∈R^{m_t\×l_t}$ 分别为第t个动力单元的系统矩阵，(当t＝1,r时l_t＝2；当t＝2,...，(r-1)时l_t＝3)；为k时刻第t个动力单元的过程噪声，分别为k时刻第t个动力单元的输出噪声。

子空间辨识是一种基于数据驱动的建模方法，由输入输出数据块构成的Hankel矩阵是子空间辨识算法的基础。基于第t个动力单元的输入输出数据，由(3)得到适用于子空间辨识的广义数据输入输出矩阵方程：

$Y_p={\mathrm{\Γ}}_i{X}_p+H_i^b{U}_p+H_i^e{R}_p+V_p---\left(4\right)$

$Y_f={\mathrm{\Γ}}_i{X}_f+H_i^b{U}_f+H_i^e{R}_f+V_f---\left(5\right)$

其中，下标p和f分别表示过去和未来时刻；U_p,U_f分别为过去和未来的输入Hankel矩阵。

$U_p=\left(\begin{array}{cccc}{u}_0& u_1& \mathrm{...}& u_{j-1}\\ u_1& u_2& \mathrm{...}& u_j\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ u_{i-1}& u_i& \mathrm{...}& u_{i+j-2}\end{array}\right),U_f=\left(\begin{array}{cccc}{u}_i& u_{i+1}& \mathrm{...}& u_{i+j-1}\\ u_{i+1}& u_{i+2}& \mathrm{...}& u_{i+j}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ u_{2i-1}& u_{2i}& \mathrm{...}& u_{2i+j-2}\end{array}\right)$

其中u_i＝[u_i1u_i2...u_il]^T，同理定义过去和未来的输出Hankel矩阵Y_p,Y_f，状态矩阵X_p和X_f也类似定义如下：

X_p＝[x₀x₁...x_j-1],X_f＝[x_jx_j+1...x_j+i]

广义能观性矩阵Γ_i和低维下三角Toeplitz矩阵H_i分别定义为：

根据以上的定义方法，可以定义下三角Toeplitz矩阵以及噪声矩阵R_p,R_f,V_p,V_f。

由式(4)和(5)，经过合适的变换，可得子空间预测输出

${\hat{Y}}_f=L_w{W}_p+L_u{U}_f+L_e{R}_f---\left(6\right)$

通过使用斜向投影消除噪声的影响，则预测值可以由下述最小二乘问题求解 $\underset{L_w,L_u}{\min }\left|\right|Y_f-(L_w,L_u)\left(\begin{array}{c}{W}_p\\ U_f\end{array}\right)|{|}_F^2---\left(7\right)$

其中，W_p＝[Y_pU_p]^T,符号||·||_F表示Frobenius范数。

式(7)可通过采用QR分解来实现得到：

$\left(\begin{array}{c}{W}_p\\ U_f\\ Y_f\end{array}\right)=\left(\begin{array}{ccc}{L}_{11}& 0& 0\\ L_{21}& L_{22}& 0\\ L_{31}& L_{32}& L_{33}\end{array}\right)\left(\begin{array}{c}{Q}_1^T\\ Q_2^T\\ Q_3^T\end{array}\right)---\left(8\right)$

记

式中，L＝[L_wL_u]。

将求得L_w,L_u代入式(6)，就可以记算出同理基于每个动力单元的输入输出数据，辨识出各单元对应的子空间模型，继而组合各单元子空间模型得到高速动车组分布式状态空间模型，即适用于同步跟踪控制算法的预测模型：

$\left\{\begin{array}{c}x\left(k+1\right)\≈Ax\left(k\right)+Bu\left(k\right)+w\left(k\right)\\ y\left(k\right)\≈Cx\left(k\right)+Du\left(k\right)+v\left(k\right)\end{array}\right.---\left(9\right)$

其中，u(k)∈R^l为k时刻预测模型的输入序列，y(k)∈R^m为k时刻预测模型的输出序列，x(k)∈Rⁿ为k时刻预测模型的状态序列，x(k+1)∈Rⁿ为k+1时刻预测模型的状态序列；A∈R^n×n，B∈R^n×l，C∈R^m×n，D∈R^m×l分别为预测模型的参数矩阵；w(k)∈Rⁿ为k时刻预测模型的过程噪声，v(k)∈R^m为k时刻预测模型的输出噪声。

2、高速动车组各单元同步跟踪控制方法：

本发明提出基于分布式模型的同步跟踪控制方法来实现高速动车组各单元的高精度同步跟踪控制；以模型输出和期望输出的速度偏差设计目标函数。

高速动车组状态空间模型的同步跟踪控制框图如图2所示。采用子空间模型架构，r个动力单元分别根据自身输入，相邻单元的输入以及自身输出数据(u_t-1,u_t,u_t+1,y_t)辨识出其各自的模型。在每个采样时刻，将分布式状态空间模型切入闭环系统来调整控制器的参数，从而实现各单元同步跟踪，使各单元间的速度误差满足目标要求，各单元的控制力变化平缓；动车组安全平稳运行和舒适性指标也得以改善。

基于式(9)得到的高速动车组运行过程分布式状态空间模型，设计同步跟踪控制器；以模型输出和期望输出的速度偏差设计目标函数，即：

$J=E\left\{\underset{j=N_0}{\overset{P}{\Σ}}{\[\hat{y}\left(k+j\right)-y_r\left(k+j\right)\]}^2+\underset{j=1}{\overset{N}{\Σ}}{\[r_{j}u\left(k+j\right)\]}^2\right\}---\left(10\right)$

式(10)中，为模型超前j步的最优预测值，y_r(k+j)为未来(k+j)时刻的期望输出，N₀是最小输出长度、P是预测长度、N是控制长度，一般情况下取N₀＝1。r_j为控制加权系数，约束控制量；u_f(k+j)为控制量序列。

将式(10)用矩阵形式表示为：

$J={({\hat{y}}_f-Y_r)}^{T}Q({\hat{y}}_f-Y_r)+{u_f}^T{\mathrm{Ru}}_f---\left(11\right)$

式(11)中，矩阵矩阵

Y_r＝[y_r(k+1)y_r(k+2)...y_r(k+P)]^T，矩阵

u_f＝[u_f(k+1)u_f(k+2)...u_f(k+N)]^T,(·)^T表示矩阵的转置，Q为输出误差加权矩阵，R为控制加权矩阵，Q、R均为正定对称矩阵，即Q＝Q^T＞0，R＝R^T＞0。

根据动车组预报模型(9)实现对其输出预测：

${\hat{y}}_f=l_w{w}_p+l_u{u}_f---\left(12\right)$

式(12)中， $w_p={\[y_p^T,u_p^T\]}^T\left(y_p=\[y_t^T,...,y_{t-P+1}^T\],u_p=\[u_t^T,...,u_{t-P+1}^T\]\right);$ l_w＝L_w(1:lP,:),l_u＝L_u(1:lP,1:n_u)(L_u、L_w为子空间矩阵；n_u为u_f的行数，与控制时域N相关)。

将式(12)代入式(11)得到：

J＝((l_ww_p+l_uu_f)-Y_r)^TQ((l_ww_p+l_uu_f)-Y_r)+u_f ^TRu_f(13)

由式(13)可知求目标函数的最小化问题等价于求解关于向量u_f的极小值问题。令dJ/du_f＝0可得到最优控制律u_f的表达式：

$u_f={(l_u^T{\mathrm{Ql}}_u+R)}^{-1}{l}_u^{T}Q(Y_r-l_w{w}_p)---\left(14\right)$

式(14)中，L_u为子空间矩阵，l_u＝L_u(1:lP,1:n_u)；L_w为子空间矩阵，l_w＝L_w(1:lP,:)；Q为输出误差加权矩阵；R为控制加权矩阵；表示预测时域内的输出、输入数据矩阵；表示预测时域内输出数据矩阵的转置；表示预测时域内输入数据矩阵的转置；Y_r为期望输出的矩阵形式；(·)^-1表示矩阵求逆运算。

取u_f的第一、二、…、r行分量施加到高速动车组上，即为当前时刻各个动力单元所需控制量u＝u_f(1:r,:)。

综上所述，针对动车组编组结构和运行过程特点,建立了动车组分布式状态空间模型,提出基于高速动车组分布式模型的同步跟踪控制方法,实现高速动车组各单元高精度同步跟踪控制，满足高速动车组安全、平稳、高精度跟踪控制需求。

本发明与现有技术比较的有益效果是，高速动车组各动力单元的动态特性具有非线性、强耦合的特点，一般采用集中式模型和分布式静态模型，难以建立精确的控制模型，从而影响各动力单元跟踪控制的效果。本技术方案首先根据动车组运行数据构造输入输出Hankel矩阵，然后对Hankel矩阵进行投影计算得到状态向量，通过求解一个最小二乘问题求出状态空间模型的参数矩阵[A,B,C,D]。提出基于分布式状态空间模型的同步跟踪控制方法，得到各单元的准确控制量，使高速动车组各单元运行速度精确跟踪目标曲线，实现动车组正点、安全、有效运行，保证了乘客安全。本技术方案简单实用，可为动车组乘务员操纵提供可靠的信息。

本发明适用于高速动车组运行过程建模和同步跟踪控制。

附图说明

图1是高速动车组分布式描述示意图；

图2是高速动车组状态空间模型的同步跟踪控制框图；

图3是CHR380AL型动车组七个分布式控制单元；

图4为各控制单元速度跟踪曲线(图4－1至图4－7)；

图4-1是单元1与给定速度跟踪曲线；

图4-2是单元2与给定速度跟踪曲线；

图4-3是单元3与给定速度跟踪曲线；

图4-4是单元4与给定速度跟踪曲线；

图4-5是单元5与给定速度跟踪曲线；

图4-6是单元6与给定速度跟踪曲线；

图4-7是单元7与给定速度跟踪曲线

图5为各控制单元速度跟踪误差曲线(图5－1至图5－7)；

图5-1是单元1与给定速度跟踪误差曲线；

图5-2是单元2与给定速度跟踪误差曲线；

图5-3是单元3与给定速度跟踪误差曲线；

图5-4是单元4与给定速度跟踪误差曲线；

图5-5是单元5与给定速度跟踪误差曲线；

图5-6是单元6与给定速度跟踪误差曲线；

图5-7是单元7与给定速度跟踪误差曲线；

图6为相邻控制单元速度偏差(图6－1至图6－6)；

图6-1是单元1与单元2速度偏差；

图6-2是单元2与单元3速度偏差；

图6-3是单元3与单元4速度偏差；

图6-4是单元4与单元5速度偏差；

图6-5是单元5与单元3速度偏差；

图6-6是单元6与单元3速度偏差；

图7为各单元控制力变化图(图7－1至图7－2)；

图7-1是单元1、2、3、4控制力变化曲线；

图7-2是单元5、6、7控制力变化曲线；

图8为相邻单元控制力偏差(图8－1至图8－6)；

图8-1是单元1与单元2的相邻单元控制力偏差；

图8-2是单元2与单元3的相邻单元控制力偏差；

图8-3是单元3与单元4的相邻单元控制力偏差；

图8-4是单元4与单元5的相邻单元控制力偏差；

图8-5是单元5与单元6的相邻单元控制力偏差；

图8-6是单元6与单元7的相邻单元控制力偏差。

具体实施方式

本发明实施选用十四动两拖(第一节和第十六节为拖车，其余均为动车)CHR380AL型动车组为研究对象，分布图如图3所示。采集该型号动车组运行于徐州东济南区间段的实际数据，结合高速动车组牵引/制动特性曲线，开展分布式建模及运行控制仿真验证。采用子空间模型辨识方法对每个动力单元的1600组输入输出数据样本进行辨识，得到高速动车组分布式状态空间预测模型，各动力单元的状态矩阵及初始状态分别如下所示：

单元1的模型：

A₁＝[0.99886]，B₁＝[4.2142e-005-4.1864e-005]

C₁＝[6730.3]，D₁＝[00]，K₁＝[0.00015309]，x₁(0)＝[9.0846e-010]

单元2的模型：

A₂＝[0.99862]，B₂＝[-4.5886e-0052.1897e-0064.398e-005]

C₂＝[6625.7]，D₂＝[000]，K₂＝[0.00015496]，x₂(0)＝[-2.1127e-006]

单元3的模型：

A₃＝[0.9986]，B₃＝[-0.000154235.2682e-0050.00010183]

C₃＝[6598.3]，D₃＝[000]，K₃＝[0.0001566]，x₃(0)＝[2.6075e-006]

单元4的模型：

A₄＝[0.99862]，B₄＝[-9.5415e-0055.1916e-0069.0508e-005]

C₄＝[6628.9]，D₄＝[000]，K₄＝[0.00015458]，x₄(0)＝[-3.1365e-006]

单元5的模型：

A₅＝[0.99885]，B₅＝[-3.0691e-0055.2487e-005-2.1568e-005]

C₅＝[6611.8]，D₅＝[000]，K₅＝[0.00015622]，x₅(0)＝[4.1115e-006]

单元6的模型：

A₆＝[0.99885]，B₆＝[0.00018681-3.8988e-005-0.00014753]

C₆＝[6626.1]，D₆＝[000]，K₆＝[0.00015413]，x₆(0)＝[-4.0784e-006]

单元7的模型：

A₇＝[0.99859]，B₇＝[-3.7637e-0053.7864e-005]

C₇＝[6650]，D₇＝[00]，K₇＝[0.00010622]，x₇(0)＝[5.8271e-006]

根据以上七个单元的模型，将其集成建立动车组状态空间模型，模型状态矩阵和初始状态如下：

A＝diag[A₁A₂...A₇]，B＝[B₁B₂...B₇]^T，C＝diag[C₁C₂...C₇]

D＝[D₁D₂...D₇]^T，x(0)＝[x₁(0)x₂(0)...x₇(0)]^T

采用本发明方法对CRH380AL动车组运行过程进行同步跟踪控制(初始里程为693.31km/h，终点里程为393.38km/h)。图4描述了各个控制单元的速度跟踪曲线，如图4－1至图4－7所示。相应的误差曲线如图5所示，包括图5－1至图5－7。表1为各控制单元速度跟踪误差范围。

表1各控制单元速度跟踪误差范围

控制单元号	误差(km/h)
		单元一	(-0.1150,0.0944)
单元二	(-0.0754，0.0620)
		单元三	(-0.0352,0.0289)
单元四	(-0.0704，0.0578)
		单元五	(-0.0872，0.0717)
单元六	(-0.0662，0.0544)
		单元七	(-0.0485,0.0399)

根据图4、图5和表2可知，各控制单元的输出在各个运行工况下均有良好的跟踪能力，速度跟踪误差均在±0.1km/h以内，能以较高精度满足运行控制要求。

为了更好地描述各个单元的相对运动性能，相邻单元间的速度偏差如图6，包括图6－1至图6－6所示。表2为相邻单元速度偏差变化范围；图7为各个控制单元的控制力，如图7－1和图7－2所示；图8为相邻控制单元的控制力偏差，如图8－1至图8－6所示；表3为相邻单元控制力偏差变化范围。

表2相邻单元间速度偏差变化范围

控制单元号	误差(km/h)
		单元一与二	(-0.0325,0.0396)
单元二与三	(-0.0331,0.0402)
		单元三与四	(-0.0352,0.0289)
单元四与五	(-0.0169,0.0138)
		单元五与六	(-0.0172,0.0210)
单元六与七	(-0.0146,0.0177)

由图6和表2可以看出，在牵引、惰行、恒速和制动工况下，相邻单元的速度误差在±0.05km/h以内，车辆间的相对位移较小，有效降低了车辆间的耦合效应，高速动车组平稳运行性能得以提高。

表3相邻单元控制力偏差变化范围

控制单元号	误差(kN)
		单元一与二	(-4.4168,5.1083)
单元二与三	(-4.8839,6.0403)
		单元三与四	(-4.7788,6.2111)
单元四与五	(-6.9792,5.8581)
		单元五与六	(-4.2498,3.5242)
单元六与七	(-3.7858,3.2765)

根据图7、图8和表3可知，各个控制单元的控制力在整个运行过程缓和变化，过渡比较平滑；相邻单元间的控制力误差范围在±7kN以内，乘客舒适性指标得到一定程度的提高。

Claims (1)

1.一种基于分布式模型的动车组同步跟踪控制方法，其特征是，所述方法针对动力分散式动车组由若干牵引/制动控制单元相互耦合组成的结构特点，提出动车组运行过程的分布式描述方法；借助动车组牵引/制动特性曲线和实际运行数据，采用子空间模型辨识方法建立了动车组分布式状态空间模型；提出基于高速动车组分布式模型的同步跟踪控制方法，实现高速动车组各单元运行速度的高精度同步跟踪；

所述分布式模型的同步跟踪控制方法包括以下步骤：

(1)基于每个动力单元的输入输出数据，辨识出各单元对应的子空间模型，继而组合各单元子空间模型得到高速动车组分布式状态空间模型，即适用于同步跟踪控制算法的预测模型：

$\left\{\begin{array}{c}x(k+1)\≈Ax\left(k\right)+Bu\left(k\right)+w\left(k\right)\\ y\left(k\right)\≈Cx\left(k\right)+Du\left(k\right)+v\left(k\right)\end{array}\right.---\left(1\right)$

其中，u(k)∈R^l为k时刻预测模型的输入序列；y(k)∈R^m为k时刻预测模型的输出序列；x(k)∈Rⁿ为k时刻预测模型的状态序列；A∈R^n×n，B∈R^n×l，C∈R^m×n，D∈R^m×l分别为预测模型的参数矩阵；w(k)∈Rⁿ为k时刻预测模型的过程噪声；v(k)∈R^m为k时刻预测模型的输出噪声；

(2)提出基于分布式模型的同步跟踪控制方法来实现高速动车组各单元的高精度同步跟踪控制；以模型输出和期望输出的速度偏差设计目标函数：

$J=E\{\underset{j=N_0}{\overset{P}{\mathrm{\Σ}}}{\[\widehat{y_f}(k+j)-y_r(k+j)\]}^2+\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{\[r_j{u}_f(k+j)\]}^2\}---\left(2\right)$

式(2)中，为模型超前j步的最优预测值；y_r(k+j)为未来(k+j)时刻的期望输出；N₀是最小输出长度；P是预测长度，N是控制长度，取N₀＝1；r_j为控制加权系数，约束控制量；u_f(k+j)为控制量序列；

将式(2)用矩阵形式表示为：

$J={({\hat{y}}_f-Y_r)}^{T}Q({\hat{y}}_f-Y_r)+{u_f}^T{\mathrm{Ru}}_f---\left(3\right)$

式(3)中，矩阵矩阵 $Y_r={\left[\begin{array}{cccc}{y}_r(k+1)& y_r(k+2)& ...& y_r(k+P)\end{array}\right]}^T;$ 矩阵 $u_f={\left[\begin{array}{cccc}{u}_f(k+1)& u_f(k+2)& ...& u_f(k+N)\end{array}\right]}^T;$ (·)^T表示矩阵的转置；Q为输出误差加权矩阵；R为控制加权矩阵；Q、R均为正定对称矩阵，即Q＝Q^T＞0，R＝R^T＞0；

根据动车组预报模型(1)实现对其输出预测：

${\hat{y}}_f=l_w{w}_p+l_u{u}_f---\left(4\right)$

式(4)中， $w_p={\[y_p^T,u_p^T\]}^T;$ $y_p=\[y_t^T,...,y_{t-P+1}^T\],u_p=\[u_t^T,...,u_{t-P+1}^T\];$ l_w＝L_w(1:lP,:)；l_u＝L_u(1:lP,1:n_u)，L_u、L_w为子空间矩阵；n_u为u_f的行数，与控制时域N相关；

将式(4)代入式(3)得到：

J＝((l_ww_p+l_uu_f)-Y_r)^TQ((l_ww_p+l_uu_f)-Y_r)+u_f ^TRu_f(5)

由式(5)可知求目标函数的最小化问题等价于求解关于向量u_f的极小值问题；令dJ/du_f＝0可得到最优控制律u_f的表达式：

$u_f={(l_u^T{\mathrm{Ql}}_u+R)}^{-1}{l}_u^{T}Q(Y_r-l_w{w}_p)---\left(6\right)$

取u_f的第一、二、…、r行分量施加到高速动车组上，即为当前时刻各个动力单元所需控制量u＝u_f(1:r,:)；

式(6)中，L_u为子空间矩阵，l_u＝L_u(1:lP,1:n_u)；L_w为子空间矩阵，l_w＝L_w(1:lP,:)；Q为输出误差加权矩阵；R为控制加权矩阵；表示预测时域内的输出、输入数据矩阵；表示预测时域内输出数据矩阵；表示预测时域内输入数据矩阵；Y_r为期望输出的矩阵形式；(·)^-1表示矩阵求逆运算。