随着汽车保有量越来越大，石油资源短缺、环境污染问题日益严重. 2016年，我国开始实施第四阶段乘用车燃油消耗标准，目标至2020年实现国产乘用车平均油耗降至5 L/100 km. 2019年，我国各大城市开始逐步实施国六排放法规^[1]. 在节能减排的大背景下，各汽车厂商采用发动机优化技术和电气化技术使车辆满足法规要求. 在全面进入电动车时代之前，混合动力技术得到了广泛的关注及应用.

能量管理策略是决定混合动力车辆燃油经济性的关键因素，目前混合动力能量管理策略主要分成3大类：基于规则的控制策略^[2-3]、局部优化的控制策略^[4-5]和全局优化控制策略^[6-8]. 其中局部优化控制策略中的模型预测控制算法^[9-10]（model predictive control, MPC）的运算速度较快，具备实时应用的潜力，因此受到越来越多的关注.

目前已有不少关于模型预测控制算法的研究，如张洁丽^[11]构建模型预测控制算法的框架，采用动态规划算法进行求解，该求解方法的速度较慢，所需的存储空间较大；Borhan等^[12]将非线性模型线性化处理成线性模型，采用线性时变模型预测控制算法进行求解，但是没有对驾驶员需求转矩/车速进行预测，直接将标准工况的车速输入到预测模型中进行求解，而实际工况不能完全准确地获得；Sun等^[13-14]采用指数预测方法、神经网络预测方法对驾驶员需求转矩进行预测，其中指数预测方法计算简单、预测精度较低，神经网络预测方法计算量较大、精度较高.

本文采用自回归模型（autoregressive exogenous model, ARX）对驾驶员需求转矩进行预测，根据预测得到的驾驶员需求转矩序列，通过纵向动力学公式求得未来的车速序列；将非线性的整车模型线性化处理，将非线性模型预测控制转化为线性时变模型预测控制（linear time varying model predictive control，LTV-MPC）问题进行求解. 根据当前时刻的状态量，在保证代价函数最小的前提下，获得预测时域内最优的发动机、电机转矩分配序列，将序列中的第1组发动机、电机转矩作用于被控系统，利用整车模型可以获得下一时刻的状态量. 不断迭代重复上述过程，可以获得整个工况的最优转矩分配序列.

本文主要的创新点如下：对比分析指数预测和自回归模型预测这2种方法；结合自回归模型预测和模型预测控制算法，分析不同预测时域下车辆燃油经济性的变化趋势；假设工况已知，分析不同预测时域下车辆燃油经济性的变化趋势；将自回归模型预测结果和工况已知时的车辆燃油经济性进行对比，分析预测信息准确时的节油潜力.

本文的研究对象为并联式混合动力车辆，如图1所示，包括发动机、电机、动力电池和其他附件等，电机处于离合器和变速器之间，可以利用变速器档位来优化电机的工作点，提升运行效率. 该结构对动力总成的改动相对较小，整车成本的增加相对较少，是当前各大汽车厂家普遍采用的混合动力结构形式.

对象车辆的主要参数如表1所示.

车辆模型包括驾驶员、转矩预测/车速预测、模型预测控制、发动机、电机、电池、变速器和纵向动力学模块，各模块之间的相互关系如图2所示. 根据驾驶员需求转矩预测模块可得预测时域内的驾驶员需求转矩，根据纵向动力学模型可得预测时域内的车速. 根据预测车速，利用模型预测控制模块计算得到的最优发动机、电机转矩，将其输入到发动机、电机控制模块计算得到发动机油耗和电池充放电量. 在实车上，模型预测控制模块被编译并导入到整车控制器中，计算得到的动力部件的需求转矩信号发送给发动机、电机控制模块控制器中，由各动力部件控制器进行动态控制. 由于在能量管理中主要考虑动力部件的稳态特性，且动力部件的瞬态响应数据不易获得，在对动力部件进行建模时主要考虑稳态特性，忽略频率高于1 Hz的瞬态特性.

驾驶员模块采用PID模型，数学模型为

(1) $ \alpha = {K_{\rm{P}}}\left( {e + {K_{\rm{I}}}\int e {\rm{d}}t + {K_{\rm{D}}}({\rm{d}}e/{\rm{d}}t)} \right). $

式中：α为油门/制动踏板开度，K_P、K_I和K_D分别为比例、积分和微分常数. 根据目标车速和实际车速之间的差值e，通过比例、积分和微分3个部分相加，可得油门/制动踏板开度. 在获得了油门/制动踏板开度之后，根据需求转矩计算模块可得驾驶员需求转矩.

电池模型用于求解电池荷电状态（state of charge, SOC），本文采用RC电池模型，计算公式为

(2) $ {\rm{S\dot OC}} = \frac{{{V_{{\rm{oc}}}} - \sqrt {V_{{\rm{oc}}}^2 - 4{R_{{\mathop{\rm int}} }}{P_{{\rm{batt}}}}} }}{{2{R_{{\mathop{\rm int}} }}{Q_{{\rm{batt}}}}}}. $

式中：V_oc为电池开路电压，R_int为电池充放电内阻，Q_batt为电池容量. 其中电池功率P_batt的计算公式为

(3) $ {P_{{\rm{batt}}}} = {T_{{\rm{mot}}}}{\omega _{{\rm{mot}}}}{\eta ^{-{\rm{sgn}}\;({T_{{\rm{mot}}}})}_{\rm{mot}}} . $

式中：T_mot为电机转矩，ω_mot为电机转速，η_mot为电机效率. 为了简化计算，假设电池开路电压、电池充放电内阻和电池容量都为定值，不随电池SOC的变化而变化.

发动机模型基于台架测试数据构建，发动机万有特性如图3所示，包括比油耗等高线、转矩外特性曲线. 图中，n为转速.

电机模型基于台架测试数据构建，电机万有特性如图4所示，包括效率等高线、转矩外特性曲线.

各控制策略中的变速器换挡策略根据最小发动机比油耗曲线来确定，如图5所示. 图中，V为车速，θ为油门开度.

车辆纵向动力学用于求取车辆的纵向速度，表达式为

(4) $ M\dot V = ({T_{{\rm{eng}}}} + {T_{{\rm{mot}}}}) {i_{\rm{g}}}(\gamma ) {i_0} {\eta _{\rm{t}}}/{R_{\rm{w}}} - {F_{\rm{b}}} - {F_{\rm{r}}}. $

式中：M为整备质量，T_eng为发动机转矩， $\gamma $为变速器档位， ${i_{\rm{g}}}(\gamma )$为档位 $\gamma $对应的变速器传动比，η_t为机械传动效率，i₀为主减速器传动比，R_w为车轮半径，F_b为制动力，F_r为行驶阻力之和.

发动机和电机数据由台架测试获得，都是散点数据. 为了方便后续进行线性化处理，在Matlab中对发动机和电机数据进行多项式拟合. 将发动机油耗曲面拟合为关于发动机转矩3次方、转速4次方的多项式曲面，具体如式（5）、图6所示. 图中，ϕ为发动机瞬时油耗. 拟合数据和原数据之间的均方根误差为0.098 g/s；将发动机外特性转矩拟合为关于发动机转速的6次多项式曲线，具体如式（6）、图7所示，拟合数据和原数据之间的均方根误差为2.143 N·m；将电机效率曲面拟合为关于电机转矩2次方、电机转速5次方的多项式曲面，具体如式（7）、图8所示，拟合数据和原数据之间的均方根误差为0.012；电机外特性转矩可以拟合为关于转速的6次多项式曲线，具体如式（8）、图7所示，拟合数据和原数据之间的均方根误差为3.182 N·m. 通过该多项式拟合的结果误差都非常小，可以满足后续计算精度的要求.

(5) $ \begin{split} {{\dot m}_{{\rm{fuel}}}} = & {f_{00}} + {f_{10}}{\omega _{{\rm{eng}}}} + {f_{01}}{T_{{\rm{eng}}}} + {f_{20}}\omega _{{\rm{eng}}}^2 + {f_{11}}{\omega _{{\rm{eng}}}}{T_{{\rm{eng}}}} + \\ & {f_{02}}T_{{\rm{eng}}}^2 + {f_{30}}\omega _{{\rm{eng}}}^3 + {f_{21}}\omega _{{\rm{eng}}}^2{T_{{\rm{eng}}}} + {f_{12}}{\omega _{{\rm{eng}}}}T_{{\rm{eng}}}^2 + \\ & {f_{03}}T_{{\rm{eng}}}^3 + {f_{40}}\omega _{{\rm{eng}}}^4 + {\rm{ }}{f_{31}}\omega _{{\rm{eng}}}^3{T_{{\rm{eng}}}} + {f_{22}}\omega _{{\rm{eng}}}^2T_{{\rm{eng}}}^2 + \\ & {f_{13}}{\omega _{{\rm{eng}}}}T_{{\rm{eng}}}^3. \\[-10pt] \end{split} $

(6) $ \begin{split} T_{{\rm{eng}}}^{\max } = & {a_1}w_{{\rm{eng}}}^6 + {a_2}w_{{\rm{eng}}}^5 + {a_3}w_{{\rm{eng}}}^4 + {a_4}w_{{\rm{eng}}}^3 + \\ & {a_5}w_{{\rm{eng}}}^2 + {a_6}{w_{{\rm{eng}}}} + {a_7}.\\[-10pt] \end{split}\qquad\qquad\quad\; $

(7) $\begin{split} {\eta _{{\rm{mot}}}} = & {p_{00}} + {p_{10}}{\omega _{{\rm{mot}}}} + {p_{01}}{T_{{\rm{mot}}}} + {p_{20}}\omega _{{\rm{mot}}}^2 + {p_{11}}{\omega _{{\rm{mot}}}}{T_{{\rm{mot}}}} + \\ & {p_{02}}T_{{\rm{mot}}}^2 + {p_{30}}\omega _{{\rm{mot}}}^3 + {p_{21}}\omega _{{\rm{mot}}}^2{T_{{\rm{mot}}}} + {p_{12}}{\omega _{{\rm{mot}}}}T_{{\rm{mot}}}^2 + \\ & {p_{40}}\omega _{{\rm{mot}}}^4 + {p_{31}}\omega _{{\rm{mot}}}^3{T_{{\rm{mot}}}} + {p_{22}}\omega _{{\rm{mot}}}^2T_{{\rm{mot}}}^2 + {p_{50}}\omega _{{\rm{mot}}}^5 + \\ & {p_{41}}\omega _{{\rm{mot}}}^4{T_{{\rm{mot}}}} + {p_{32}}\omega _{{\rm{mot}}}^3T_{{\rm{mot}}}^2. \\[-10pt] \end{split} $

(8) $ \begin{split} T_{{\rm{mot}}}^{\max } = & {b_1}w_{{\rm{mot}}}^6 + {b_2}w_{{\rm{mot}}}^5 + {b_3}w_{{\rm{mot}}}^4 + {b_4}w_{{\rm{mot}}}^3 + \\ & {b_5}w_{{\rm{mot}}}^2 + {b_6}{w_{{\rm{mot}}}} + {b_7}. \\[-10pt] \end{split} \qquad\qquad\quad\; $

式（5）中的系数如下：f₀₀=0.39，f₁₀=−2.7×10⁻³，f₀₁=−2.53×10⁻²，f₂₀=6.85×10⁻⁶，f₁₁=2.64×10⁻⁴，f₀₂= 3.24×10⁻⁴，f₃₀=1.59×10⁻⁹，f₂₁=−3.80×10⁻⁷，f₁₂=−2.23×10⁻⁶，f₀₃=−8.53×10⁻⁷，f₄₀=−6.11×10⁻¹²，f₃₁=1.73×10⁻¹⁰，f₂₂=2.39×10⁻⁹，f₁₃=5.20×10⁻⁹. 式（6）中的系数如下：a₁=−5.76×10⁻¹³，a₂=1.33×10⁻⁹，a₃=1.24×10⁻⁶，a₄=5.88×10⁻⁴，a₅=−0.15，a₆=19.67，a₇=−842.5. 式（7）中的系数为：P₀₀=0.68，P₁₀=2.23×10⁻³，P₀₁=3.78×10⁻⁴，P₂₀=−7.21×10⁻⁶，P₁₁=−3.76×10⁻⁶，P₀₂=−6.54×10⁻⁷，P₃₀=1.00×10⁻⁸，P₂₁=1.36×10⁻⁸，P₁₂=1.74×10⁻⁹，P₄₀=−4.82×10⁻¹²，P₃₁=−1.98×10⁻¹¹，P₂₂=2.23×10⁻¹³，P₅₀=−3.55×10⁻¹⁶，P₄₁=1.04×10⁻¹⁴，P₃₂=−1.82×10⁻¹⁵. 式（8）中的系数为b₁=−6.15×10⁻¹⁴，b₂=1.30×10⁻¹⁰，b₃=−9.82×10⁻⁸，b₄=3.11×10⁻⁵，b₅=−4.09×10⁻³，b₆=0.17，b₇=276.1.

联立动力电池SOC方程（式（2））、车辆动力学方程（式（4））和发动机油耗方程（式（5）），可得如式（（9））所示的整车模型表达式，需要满足发动机最大转矩（式（6））、电机最大最小转矩（式（8））的约束.

(9) $ \left. \begin{array}{l} {{\dot x }}= f{{(x,u,v),}} \\ {{y}} = g({x,u,v}). \end{array} \right\} $

式中：x为状态变量，u为控制变量，v为控制变量，y为输出量.

${x} = \left[ \begin{array}{l} V \\ {\rm{SOC}} \\ {m_{{\rm{fuel}}}} \\ \end{array} \right]{\rm{ }},{\rm{ }}{{u}} = \left[ \begin{array}{l} {T_{{\rm{eng}}}} \\ {T_{{\rm{mot}}}} \\ \end{array} \right]{\rm{ }},\;{{v}} = \left[ {{F_{\rm{r}}}} \right] ,\; {{y}} = \left[ \begin{array}{l} V \\ {\rm{SOC}} \\ {m_{{\rm{fuel}}}} \\ \end{array} \right].$

在实际车辆中，驾驶员通过油门、刹车踏板来控制车辆车速，能量管理控制器是在满足驾驶员需求转矩的前提下进行能量管理，车速根据动力部件的输入/输出转矩不断变化. 此外，预测驾驶员需求转矩时还可以结合加速踏板开度/开度变化率等因素^[15]. 基于上述优势，虽然驾驶员需求转矩和车速之间可以通过如式（4）所示的车辆纵向动力学公式进行转化，最终选择对驾驶员需求转矩进行预测.

指数函数预测公式如下所示：

10) $ {T_{{\rm{req}}}}(k + i) = {T_{{\rm{req}}}}(k)\times{\exp\;{( - i/{T_{\rm{d}}})}};\;i = 1,2, \cdots ,{T_{\rm{p}}}. $

式中：T_req（k+i）为k时刻预测得到的k+i时刻的需求转矩，T_d为指数函数的衰退因子.

为了判断转矩预测的准确度，选用误差均方根作为评判指标，表达式为

(11) $R(k) = \sqrt {\sum\limits_{i = 1}^{{T_{\rm{p}}}} {{{\left( {T_{{\rm{req}}}^{{\rm{act}}}(k + i) - T_{{\rm{req}}}^{{\rm{ref}}}(k + i)} \right)}^2}} /{T_{\rm{p}}}} ,$

(12) $ {R_{\rm{e}}} = \sqrt {\sum\limits_{k = 1}^L {{R^2}} (k)/L} . $

式中：R（k）为k时刻预测时域内的转矩均方根误差；R_e为整个循环工况的转矩均方根误差，R_e越小则预测误差越小，预测准确度越高； $ {T}_{\mathrm{r}\mathrm{e}\mathrm{q}}^{\mathrm{r}\mathrm{e}\mathrm{f}}(k+i) $为 $k$时刻预测得到的k+i时刻的需求转矩； $ {T}_{\mathrm{r}\mathrm{e}\mathrm{q}}^{\mathrm{a}\mathrm{c}\mathrm{t}}(k+i) $为k+i时刻的实际需求转矩；L为循环工况的总步长；T_p为预测步长.

基于指数函数预测的转矩误差均方根如图9所示. 可知，随着T_p的增加，R_e会不断增加，这是因为随着预测步长的增加，驾驶员需求转矩的变化更没有规律可循，预测的精度降低. 在相同的预测时域下，增大T_d可以使得R_e逐步减小，但随着T_d的增加, R_e减小的幅度也会减少，慢慢趋于一个定值.

由于指数函数预测的公式是固定的，预测得到的整体趋势都是一致的，不能很好地反映实际变化趋势，因此预测精度相对较低. 目前除了指数预测，采用较多的预测方法还有马尔科夫预测、神经网络预测等. 除了上述预测方法外，自回归模型^[16-17]被应用于车速预测，在获得较好的预测效果时计算量相对较小，因此采用该算法对驾驶员需求转矩进行预测. 该算法的原理如下：根据过去一段时间内的输入输出量采用ARX模型进行拟合，通过求取拟合误差最小，获得最优的ARX模型系数，根据ARX模型预测未来一段时间的输出.

ARX模型的表达式为

(13) $ {y_{\rm{p}}}(k + m) = \sum\limits_{i = 0}^{n - 1} {\left( {{a_i} y(k - i) + {b_i} u(k - i)} \right)} . $

式中：a_i、b_i分别为输出量和控制量的常量系数，基于当前时刻k至过去n−1时刻的控制量u（k−i）及输出量y（k−i），预测将来m时刻的输出量y_p（k+m）；控制量u为驾驶员油门开度，输出量y为需求转矩.

式（13）可以整理成最小二乘形式:

(14) $ {y_{\rm{p}}}(k+m) ={{ \varphi }}(k){{\theta}} + \xi (k). $

式中：

$\left. \begin{array}{l} {{\varphi}} (k) = [ - y(k), \cdots , - y(k - n + 1),u(k), \cdots ,u(k - n + 1)], \\ {{\theta}} = {\left[ {{a_1}, \cdots ,{a_n},{b_1}, \cdots ,{b_n}} \right]^{\rm{T}}} . \end{array} \right\}$

根据上述推导过程可知，随着时刻k的推移，可得y_p（k+1）至y_p（k+m）的预测值，整理成矩阵的形式：

15) $ {{{Y}}_{\rm{p}}}(k + m) = {{\varPhi}} (k){{\theta}} . $

式中：

$ {{{Y}}_{\rm{p}}}(k + m) = {[y(k + 1), \cdots ,y(k + m)]^{\rm{T}}}, $

$ {{\varPhi}} (k) = \left[ {\begin{array}{*{20}{c}} {{{\varphi}} (k - N + 1)} \\ {{{\varphi}} (k - N + 2)} \\ \vdots \\ {{{\varphi}} (k)} \end{array}} \right]. $

根据实际输出值和预测输出值之间的误差，构建代价函数：

(16) $ J(k) = {\left\| {{{Y}}_{\rm{p}}(k) - {{\varPhi}} (k - m){{ \theta} }} \right\|^2}. $

对代价函数进行求导后，令求导后的表达式等于零，可得最优系数向量的表达式：

(17) ${{\theta}} = {\left[ {{{{\varPhi}} ^{\rm{T}}}(k - m){{\varPhi}} (k - m)} \right]^{ - 1}}{{{\varPhi}} ^{\rm{T}}}(k - m){{Y}}_{\rm{p}}(k).$

将最优系数向量代入式（13），可得k+1～k+m的预测输出值.

对不同个数的输出量n_a和控制量n_b时的预测误差进行仿真分析，预测误差如表2所示. 表中，加粗的数据为不同预测步长下的最小预测误差，可得当预测步长N=n_a+n_b=24，ARX模型的总体预测误差最小，为了简化计算，直接取n_a=n_b=12. 将N=24时的预测结果和指数函数预测的结果（T_d=40时）进行对比，以预测步长20为例，对比结果如图10所示. 图中，T_dr为驾驶员需求转矩. 由图11的预测误差比对可得，当预测步长较短时，采用ARX模型的预测误差小于指数函数，随着预测步长的增加，2种预测方法之间的预测误差逐渐接近. 当N<200时采用ARX模型对转矩进行预测在精度上有优势，当N>200时，采用2种算法的预测精度基本一致. 后续采用ARX模型对驾驶员需求转矩进行预测，获得预测值后采用模型预测控制算法进行求解，分析燃油经济性.

由式（9）可得整车模型具有非线性特性，直接对整车模型进行求解的速度较慢，且容易陷入局部最优，不一定能够获得全局最优解. 将非线性模型进行线性化处理^[18-19]，转换成线性时变模型，以提高求解速度. 将式（9）中的状态方程f（x,u）在参考状态轨迹附近展开成泰勒级数：

(18) $\begin{split} f({{x}},{{u}},{{v}}) = & f({{{x}}_0},{{{u}}_0},{{{v}}_0}) + {\left. {\frac{{\partial f}}{{\partial {{{x}}^{\rm{T}}}}}} \right|_{{{{x}}_0}\atop {{{{u}}_0}\atop {{{v}}_0}}}}\Delta {{x}} + {\left. {\frac{{\partial f}}{{\partial {{{u}}^{\rm{T}}}}}} \right|_{{{{x}}_0}\atop {{{{u}}_0}\atop {{{v}}_0}}}}\Delta {{u}} + \\ & {\left. {\frac{{\partial f}}{{\partial {{{v}}^{\rm{T}}}}}} \right|_{{{{x}}_0}\atop {{{{u}}_0}\atop {{{v}}_0}}}}\Delta {{v}} + \alpha (\Delta {{x}},\Delta {{u}},\Delta {{v}}). \end{split} $

式中：

${\left. {\frac{{\partial f}}{{\partial {{{X}}^{\rm{T}}}}}} \right|_{{{{x}}_0}\atop {{{{u}}_0}\atop {{{v}}_0}}}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {f_1}}}{{\partial {{{X}}_1}}}}&{\dfrac{{\partial {f_1}}}{{\partial {{{X}}_2}}}}& \cdots &{\dfrac{{\partial {f_1}}}{{\partial {{{X}}_n}}}}\\ {\dfrac{{\partial {f_2}}}{{\partial {{{X}}_1}}}}&{\dfrac{{\partial {f_2}}}{{\partial {{{X}}_2}}}}& \cdots &{\dfrac{{\partial {f_2}}}{{\partial {{{X}}_n}}}}\\ \vdots & \vdots &{}& \vdots \\ {\dfrac{{\partial {f_n}}}{{\partial {{{X}}_1}}}}&{\dfrac{{\partial {f_n}}}{{\partial {{{X}}_2}}}}& \cdots &{\dfrac{{\partial {f_n}}}{{\partial {{{X}}_n}}}} \end{array}} \right],\; {{X}} \in \{{{x}},{{u}},{{v}}\}{\rm{.}}$

$ \partial f/\partial {{X}} $为非线性状态方程在参考轨迹处的雅可比 (Jacobi) 矩阵，其中X∈{x, u, v}.

采用上述线性化方法，可以将如式（9）处理成如下所示的线性时变模型预测控制的标准形式：

19) $ \left. \begin{array}{l} {{\dot x}} = \tilde {{Ax}} + {{\tilde {{B}_{u}}{{u}}}} + {{\tilde {{B}_{v}}{{v}}}} + \tilde {{F}} , \\ {{y}} = \tilde {{Cx}} + {{\tilde {{D}}}_{u}}{{u}} + {{\tilde {{D}}}_{v}}{{v}} + \tilde {{G}} . \end{array} \right\} $

各时变系数矩阵的具体表达式为

$\begin{split} & \tilde {{A}} = {\left(\frac{{\partial f}}{{\partial {{x}}}}\right)_{({{{x}}_0},{{{u}}_0},{{{v}}_0})}},\;{{\tilde {{B}}}_u} = {\left(\frac{{\partial f}}{{\partial {{u}}}}\right)_{\left({{{x}}_0},{{{u}}_0},{{{v}}_0}\right)}}, \\ & {{\tilde {{B}}}_v} = {\left(\frac{{\partial f}}{{\partial {{v}}}}\right)_{({{{x}}_0},{{{u}}_0},{{{v}}_0})}},\;\tilde {{C}} = {\left(\frac{{\partial g}}{{\partial {{x}}}}\right)_{\left({{{x}}_0},{{{u}}_0},{{{v}}_0}\right)}}, \\ & {{\tilde {{D}}}_u} = {\left(\frac{{\partial g}}{{\partial {{u}}}}\right)_{({{{x}}_0},{{{u}}_0},{{{v}}_0})}},\;{{\tilde {{D}}}_v} = {\left(\frac{{\partial g}}{{\partial {{v}}}}\right)_{\left({{{x}}_0},{{{u}}_0},{{{v}}_0}\right)}}, \end{split} $

$\begin{split} & \tilde {{F}} = f({{{x}}_0},{{{u}}_0},{{{v}}_0}) - \tilde {{A}}{{{x}}_0} - {{\tilde {{B}}}_u}{{{u}}_0} - {{\tilde {{B}}}_v}{{{v}}_0}, \\& \tilde {{G}} = g({{{x}}_0},{{{u}}_0},{{{v}}_0}) - \tilde {{C}}{{{x}}_0} - {{\tilde {{D}}}_u}{{{u}}_0} - {{\tilde {{D}}}_v}{{{v}}_0}. \end{split} $

式（19）中的各个变量和式（9）中的一致，其中状态量x=[V, SOC, m_fuel]，控制量u=[T_eng, T_mot]，可测干扰量v=[F_r]，输出量y和状态量x一致. 将式（19）中的 $ \tilde {{{F}}} $、 $ \tilde {{{G}}} $和v合并为一项，作为可测干扰项. 将线性化之后的公式进行离散化处理，可得

(20) $ \left. \begin{array}{l} {{x}}(k + 1) = {{Ax}}(k) + {{{B}}_u}{{u}}(k) + {{{B}}_v}{{v}}(k), \\ {{y}}(k) = {{Cx}}(k) + {{{D}}_u}{{u}}(k) + {{{D}}_v}{{v}}(k). \end{array} \right\} $

线性模型预测控制的目标函数包括输出量和参考值之间的误差，控制量约束为

(21) $\begin{split} & \mathop {\min\; J}\limits_{{{u}}(t)} = \\ & \int_t^{t + \Delta t} [{({{y}}(\tau ) - {{r}}(\tau ))^{\rm{T}}}{{Q}}({{y}}(\tau ) - {{r}}(\tau )) + {{{u}}^{\rm{T}}}(\tau ){{Ru}}(\tau )]{\rm{d}}\tau . \end{split}$

在使得目标函数最小的同时，还需满足控制量、控制量增量和输出量的约束：

(22) $\left. \begin{aligned} & {{u}}_i^{\min } \leqslant {{u}}(k + i|k) \leqslant {{u}}_i^{\max }; \\ & \Delta {{u}}_i^{\min } \leqslant \Delta {{u}}(k + i|k) \leqslant \Delta {{u}}_i^{\max }; \\ & - {{\varepsilon}} + {{y}}_i^{\min } \leqslant {{y}}(k + i + 1|k) \leqslant {{y}}_i^{\max } + {{\varepsilon}},\; {{\varepsilon}} \geqslant {\bf{0}}. \end{aligned} \right\}$

当基于ARX模型对驾驶员需求转矩进行预测时，发动机转矩和电机转矩在每个时刻须等于驾驶员需求转矩，因此需要添加如下所示的等式约束：

(23) ${T_{{\rm{eng}}}}(k + i|k) + {T_{{\rm{mot}}}}(k + i|k) = {T_{{\rm{req}}}}(k + i|k).$

将式（23）的等式约束转化成如下所示的不等式约束：

(24) $\left. \begin{array}{l} {T_{{\rm{eng}}}}(k + i|k) + {T_{{\rm{mot}}}}(k + i|k) \leqslant {T_{{\rm{req}}}}(k + i|k), \\ {T_{{\rm{eng}}}}(k + i|k) + {T_{{\rm{mot}}}}(k + i|k) \geqslant - {T_{{\rm{req}}}}(k + i|k). \end{array} \right\}$

将式（20）进行迭代，可得预测区域内的各个输出量：

(25) $\begin{split} {{y}} &(k + i + 1|k) = {{C}}[{{{A}}^{i + 1}}{{x}}(k) + \\ & \sum\limits_{l = 0}^i {{{{A}}^i}} {{{B}}_u}\left( {{{u}}(k - 1) + \sum\limits_{j = 0}^l \Delta {{u}}(k + j|k)} \right) + \\ & {{{B}}_v}{{v}}(k + l|k)] + {{{D}}_v}{{v}}(k). \end{split} $

将利用式（22）求得的输出量代入目标函数中，通过整理获得标准二次规划的形式：

(26) $\left[ {\Delta {{{U}}^{{\rm{opt}}}},{{\varepsilon}} } \right] = \mathop {\arg \min }\limits_{\Delta {{U}},{{\varepsilon}} } \frac{1}{2}\Delta {{{U}}^{\rm{T}}}{{H}}\Delta {{U}} + {{{F}}^{\rm{T}}}\Delta {{U}}.$

同时满足下列不等式约束：

(27) $ {{{G}}_u}\Delta {{U}} + {{{G}}_\varepsilon }{{\varepsilon}} \leqslant {{W}} + {{Sx}}(k). $

通过求解式（27）的标准二次规划方程式，可以得到 $\Delta {{{U}}}^{{\rm{opt}}}$，由此可得k时刻的最优控制量^[20]：

(28) $ {{u}}(k) = {{u}}(k - 1) + \Delta {{{U}}^{{\rm{opt}}}}. $

针对NEDC循环工况，采用ARX模型对驾驶员需求转矩进行预测，采用线性时变模型预测控制（LTVMPC）进行求解. 将仿真结果和基于规则（rule based，RB）的控制策略进行对比分析，如图12、13所示. 其中，RB参考Advisor中的电辅助控制策略.

从图12可知，采用2种算法在各时间点t时，车速曲线表现基本一致，都可以很好地跟随目标车速V. 采用RB时SOC马上下降到目标值（SOC=0.6）附近，并一直在目标值上下浮动，当SOC低于目标值时会让发动机带动电机给电池充电；采用MPC时SOC缓慢下降，最终在制动时SOC回到目标值.

从图13可知，采用RB时，在部分时刻由于SOC低于设定值，发动机需要带动电机发电，因此发动机转矩维持在较高水平；采用MPC时，在代价函数（式（21））中考虑输出量和参考量之间误差、发动机转矩变化幅值，通过调节发动机转矩变化值的权重，可得更平缓的发动机转矩曲线，电机在部分时刻给发动机提供助力，调整发动机运行工况点.

针对NEDC、UDDS和WLTC 3个标准循环工况，将基于ARX模型预测的LTV-MPC和RB这2种控制策略进行仿真对比分析，其中对不同预测步长的MPC控制结果进行分析，分析预测步长对控制结果的影响，为了保证跟随效果，选取控制步长为预测步长的一半.

为了修正各种控制策略在仿真始末因SOC的不同而引起燃油消耗量变化，方便对比控制策略间的真实油耗，采用线性矫正法对燃油消耗量结果进行修正. 具体方法如下：通过调整起始SOC，经过仿真可以获得仿真结束时的SOC，将仿真结束时的SOC减去起始时的SOC定义为δSOC，获得δSOC>0和δSOC<0时的油耗数据点，用直线对各数据点进行拟合，通过直线表达式可得δSOC=0时的油耗，以此消除SOC出现偏差对油耗的影响.

如图14所示为针对NEDC工况采用基于规则和模型预测控制2种控制策略的仿真结果对比图. 可知，随着预测步长的增加，采用MPC控制策略时，车辆百公里油耗F呈现先小幅上升后下降再上升的趋势. 当N＜100时，采用MPC控制策略时的车辆百公里燃油消耗均低于RB控制策略. 当N=2时，采用MPC控制策略比采用RB控制策略时的车辆百公里油耗减少了4.34%；当N=10~50时，百公里油耗小幅上升，采用MPC控制策略比采用RB控制策略时的车辆百公里油耗提高了3%左右；当N=100时，采用MPC控制策略时的车辆百公里油耗上升明显，相对于采用RB控制策略只减少了0.37%.

如图15所示为针对UDDS工况，采用基于规则和模型预测控制2种控制策略的仿真结果对比图. 可知，随着预测步长的增加，采用MPC控制策略时，F呈现逐步上升的趋势. 当N<20时，采用MPC控制策略的车辆百公里油耗低于RB控制策略，当N>20时，采用MPC控制策略的车辆百公里油耗高于RB控制策略. 当N=2时，采用MPC控制策略的车辆百公里油耗比RB减少了1.95%；当N=20时，采用MPC控制策略的车辆百公里油耗基本和RB持平；当预测步长增加到50时，采用MPC控制策略时的车辆百公里油耗进一步上升，比RB增加了1.14%；当N=100时，采用MPC控制策略时的车辆百公里油耗比RB增加了4.28%.

如图16所示为WLTC工况下采用基于规则和模型预测控制2种控制策略的仿真结果对比图. 可知，随着预测步长的增加，采用MPC控制策略时的F均低于采用RB控制策略，但是和上述2种工况的趋势不同，车辆百公里油耗随着预测步长的增加，呈现先下降、后逐步上升的趋势. 当N=2时，采用MPC控制策略比采用RB控制策略时的车辆百公里油耗减少了1.36%；当N=10时，采用MPC控制策略时的车辆百公里油耗明显下降，比RB减少了2.53%；当N=100时，采用MPC控制策略时的车辆百公里油耗只比采用RB控制策略时减少了0.37%.

基于ARX转矩预测的MPC能量管理控制策略，由于预测误差随着预测步长的增加而变大，同时线性化处理会产生一定的误差，造成F随着N的增加呈现逐步增加的趋势. 在NEDC、UDDS工况下取N=2时可得最优燃油经济性，在WLTC工况下取N=10可得最优燃油经济性，因此采用该方法时，取10步以下的短时域预测步长可得最优的燃油经济性.

由于对驾驶员需求转矩进行预测时会产生预测误差，为了探究采用模型预测控制算法时的车辆燃油经济性，将标准工况的车速作为已知量输入给模型预测控制算法，相当于预测准确度为100%.

将采用准确的运行工况车速信息时和采用ARX转矩预测时，利用MPC进行能量管理控制的燃油经济性进行对比. 如图17、18所示分别为NEDC、WLTC标准循环工况下采用预测值和准确值的模型预测控制仿真结果对比图. 由于采用ARX模型对驾驶员需求转矩进行预测时，需要通过驾驶员模型获得驾驶员的历史操作信息，采用准确值时可以直接获得未来车速，不需要加入驾驶员模型，此区别造成了在N=2、采用准确值时F比采用预测值时高. 当N>10时，2种预测方法的百公里油耗趋势都是先下降，后缓慢上升. 当预测步长较短时，由于预测值和准确值之间的误差较小，采用预测值时的百公里油耗和采用准确值时较接近；随着预测步长的增加，预测误差不断增加，采用预测值时的百公里油耗和采用准确值时的差距缓慢增大.

如图19所示为UDDS下采用预测值和准确值时，采用MPC进行能量管理控制时的结果对比图. 可知，2种方法的F总体趋势相反，采用预测值时逐渐上升，采用准确值时逐渐下降. 当N较短时，由于预测值和真实值之间的误差较小，采用预测值时的百公里油耗和采用准确值时较接近；随着N的增加，预测误差不断增加，采用预测值时的百公里油耗和采用准确值时的差距不断变大.

由于预测结果会出现误差，基于预测结果和准确值的车辆燃油经济性存在一定的差距. 为了能够进一步提高车辆燃油经济性，需要进一步地提升预测精度.

（1）采用驾驶员需求转矩预测的模型预测控制算法的燃油经济性比基于规则的控制算法更好.

（2）基于驾驶员需求转矩预测的模型预测控制能量管理控制算法，在NEDC工况下，随着预测步长的增加，车辆百公里油耗呈现先上升后稍下降随后上升的趋势；在WLTC工况下，随着预测步长的增加，车辆百公里油耗呈现逐步上升的趋势；在UDDS工况下，随着预测步长的增加，车辆百公里油耗呈现先下降后逐步上升的趋势.

（3）基于模型预测控制算法，针对预测的驾驶员需求转矩和真实的驾驶员需求转矩2种情况进行分析. 在NEDC工况下，随着预测步长的增加，2种情况下的车辆百公里油耗都大致呈现先下降随后上升的趋势；在WLTC工况下，随着预测步长的增加，2种情况下的车辆百公里油耗都呈现先下降后上升的趋势；在UDDS工况下，随着预测步长的增加，基于预测的驾驶员需求转矩的车辆百公里油耗呈现逐步上升的趋势，基于准确的驾驶员需求转矩的车辆百公里油耗呈现逐步下降的趋势.

（4）当采用驾驶员需求转矩预测的模型预测控制时，燃油经济性和工况已知的模型预测控制有一定的差距，可以通过提高预测精度来提升燃油经济性.