纯电动汽车的行驶里程受到车载动力电池组容量的限制。目前，动力电池存在能量密度低、材料价格贵等问题，单纯依靠增大电池容量来提高行驶里程，会大幅增加汽车的自重和制造成本，增程式电动汽车(extended-range electric vehicles，E-REV)应运而生。E-REV利用能量辅助单元(auxiliary power unit, APU)为车辆提供电能补充，延长行驶里程^[1]。图 1为一种增程式电动汽车的动力系统结构示意图。APU由发动机与发电机组成，通过整流装置(AC/DC)与动力电池组耦合。该动力系统具有纯动力电池组驱动、APU和动力电池组混合驱动、制动能量回收和外接充电等多种能量流动模式^[2-3]。APU和动力电池组构成的双能量源为E-REV的能量管理提供了更多的灵活性，同时也带来了挑战。

E-REV与混合动力汽车(hybrid electric vehicles，HEV)、插电式混合动力汽车(plug-in hybrid electric vehicles，PHEV)的能量管理策略相似。然而，E-REV的发动机功率较小，APU不能单独满足车辆的大功率需求，在动力电池组的电量耗尽时，车辆的动力性会受到影响。因此，E-REV的能量管理策略需要考虑管理动力电池组的SOC消耗，防止行驶过程中提前耗尽电量；同时，又要在行程结束时消耗完电量，以充分利用E-REV可以外接充电的优势。常见的能量策略包括基于规则的(rule-based)和基于优化的(optimization-based)策略两大类^[4-8]。前者通过设定逻辑门限值或者模糊规则来控制动力系统部件的输出功率^[9-12]，后者通过建立并求解能量优化问题的数学模型来优化功率分配策略。LIN等^[6]最先将确定性动态规划算法(deterministic dynamic programming，DDP)用来求解HEV的能量优化分配问题。DDP算法依赖于完整的工况信息，然而在实际行驶过程中，车辆未来的需求功率是未知的。针对这一问题，LIN等^[13]建立了需求功率的Markov预测模型，提出了基于随机动态规划的HEV能量管理策略，摆脱了对完整工况信息的依赖。近年来，基于模型预测控制(model predictive control，MPC)的能量管理策略在HEV^[14-18]和PHEV^[19]上得到了应用，相关的研究工作多是基于已知工况展开的。随机模型预测控制(stochastic model predictive control，SMPC)^[20]方法的出现，使得MPC方法在用于能量管理时，不再依赖于完整的工况信息^[21]。

E-REV能量管理的实质是在保证车辆动力性的前提下，合理分配APU与动力电池组的输出功率，从而提高整车的燃油经济性。本文以一辆12 m增程式城市客车为研究对象，针对其行驶路线相对固定、一天行驶里程可估计的特点，利用分段优化的方法，将全行程的整车能量优化问题转化为瞬时的阶段性优化问题，主要工作包括：1) 统计工况信息，建立需求功率的Markov预测模型，基于需求功率的当前值，预测其未来N-1步的值；2) 实时估计车辆在剩余里程中的能量需求，将行程结束时期望的SOC值转化为阶段性的终端约束条件，并结合需求功率的预测序列，设计SOC消耗的参考曲线；3) 以费用最小为目标，同时引入跟踪SOC消耗参考曲线的惩罚项，建立动力系统瞬时功率分配的阶段性优化模型；4) 引入模型预测控制方法，在每个控制周期内，更新工况信息，重新预测需求功率序列，利用改进的DDP算法对优化模型进行求解，滚动优化并调整功率分配策略，降低需求功率的预测误差对能量管理效果的影响。

1 动力系统建模

基于控制周期ΔT，对动力电池组、APU、驱动电机以及整车进行离散化建模。

动力电池组的SOC按式(1) 演化^[22]，

$ {\rm{SOC}}\left( {k + 1} \right) = {\rm{SOC}}\left( k \right)-\frac{{{I_{{\rm{bat}}}}\left( k \right)}}{{{Q_{{\rm{bat}}}}}}\frac{{\Delta T}}{{3{\rm{ }}600}}. $

(1)

其中：Q_bat为动力电池组的额定容量; I_bat(k)为动力电池组的充放电电流(正值表示放电),

$ {I_{{\rm{bat}}}}\left( k \right) = \frac{{{V_{{\rm{ocv}}}}\left( k \right)-\sqrt {V_{{\rm{ocv}}}^2\left( k \right)-4{R_{{\rm{bat}}}}\cdot1{\rm{ }}000{P_{{\rm{bat}}}}\left( k \right)} }}{{2{R_{{\rm{bat}}}}}}. $

(2)

其中: R_bat为动力电池组的内阻；P_bat(k)为充放电功率；V_ocv(k)为开路电压，通过SOC(k)查表确定，

$ {V_{{\rm{ocv}}}}\left( k \right) = {V_{{\rm{ocv}}}}({\rm{SOC}}\left( k \right)). $

(3)

发动机的效率如图 2所示，颜色较深区域为高效区间，对应的燃油消耗率较低。为了提高APU的发电效率，本文选择了18个工作点(停机为第1个工作点)，如图 2中白色圆点所示。各工作点对应的APU输出功率和燃油消耗率如表 1所示。将APU的燃油消耗率$ {{\dot m}_{{\rm{f, steady}}}} $定义为其输出功率P_APU(k)的函数,

$ {{\dot m}_{{\rm{f, steady}}}} = {{\dot m}_{{\rm{f, steady}}}}\left( {{P_{{\rm{APU}}}}\left( k \right)} \right). $

(4)

APU通过AC/DC与动力电池组一起为驱动电机提供电能，动力系统在任意时刻都满足式(5) 的约束，

$ {P_{{\rm{APU}}}}\left( k \right) + {P_{{\rm{bat}}}}\left( k \right) = {P_{{\rm{req}}}}\left( k \right). $

(5)

其中，P_req(k)是由车辆需求功率确定的驱动电机输入端的需求功率。

车辆的能量消耗E和位移d分别为：

$ E\left( {k + 1} \right) = E\left( k \right) + {P_{{\rm{req}}}}(k)\frac{{\Delta T}}{{3{\rm{ }}600}}, $

(6)

$ d\left( {k + 1} \right) = d\left( k \right) + {u_{\rm{a}}}(k)\frac{{\Delta T}}{{3600}}. $

(7)

其中，u_a(k)为k时刻的车速。

由式(5) 给定的约束条件可知，当给定未来需求功率的预测序列时，控制APU的输出功率可实现动力系统的功率分配。定义系统的输入为u(k)=P_APU(k)，输出为y(k)=P_bat(k)，状态量为x(k)=[SOC(k) d(k) E(k)]^T，它们分别满足如下约束：

$ \mathit{\boldsymbol{X}} \buildrel \Delta \over = {\rm{\{ }}\mathit{\boldsymbol{x}}:{\rm{SO}}{{\rm{C}}_{{\rm{min}}}} \le {x_1} \le {\rm{SO}}{{\rm{C}}_{{\rm{max}}}}\}, $

(8)

$ U \buildrel \Delta \over = \left\{ \begin{array}{l} 0, {\rm{ }}11.50, {\rm{ }}15.13, {\rm{ }}19.19, {\rm{ }}20.35, {\rm{ }}24.10, \\ 25.60, {\rm{ }}28.99, {\rm{ }}32.24, {\rm{ }}35.30, {\rm{ }}37.60, {\rm{ }}39.92, \\ 42.20, {\rm{ }}44.50, {\rm{ }}48.50, {\rm{ }}50.80, {\rm{ }}53.20, {\rm{ }}55.97 \end{array} \right\}, $

(9)

$ Y \buildrel \Delta \over = \{ y:{P_{{\rm{bat}}, {\rm{min}}}} \le y \le {P_{{\rm{bat, max}}}}\} . $

(10)

其中：P_{bat, min}=-50 kW，P_{bat, max}=150 kW，SOC_min= 0.25，SOC_max=0.95。

2 未来需求功率预测 2.1 Markov预测模型

车辆在行驶过程中，需求功率P_req受到交通拥堵情况和司机驾驶习惯等因素的影响，具有随机性。基于文[13]的工作，需求功率序列为一个齐次Markov链，即下一个时刻需求功率的状态值P_j，只与当前状态值P_i有关。

$ \begin{array}{l} P\{ {P_{{\rm{req}}}}\left( {k + 1} \right) = {P_j}|{P_{{\rm{req}}}}\left( k \right) = {P_i}\} = {\pi _{i, j}}, \\ \;\;\;\;\;\;\;\;\;\sum\limits_{j = 1}^s {{\pi _{i, j}} = 1} . \end{array} $

(11)

其中：π_{i, j}为需求功率的取值从状态i转移到状态j的概率，对应的状态转移概率矩阵为

$\mathit{\boldsymbol{ \boldsymbol{\varPi} = }}{{\rm{(}}{\pi _{i, j}})_{s \times s}}. $

(12)

其中，s为需求功率离散化后的状态数。

2.2 转移概率矩阵估计

转移概率矩阵П的估计方法如下：

1) 基于仿真工况，计算需求功率序列。

本文中，需求功率P_req定义为驱动电机输入端的功率，其与轮边需求功率P_wh的关系为

$ {P_{{\rm{req}}}} = \left\{ \begin{array}{l} \frac{{{P_{{\rm{wh}}}}}}{{{\eta _{\rm{T}}}{\eta _{\rm{M}}}}}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{P_{{\rm{wh}}}} \ge 0;\;\\ {P_{{\rm{wh}}}}{\eta _{\rm{T}}}{\eta _{\rm{M}}}\left( {1-\beta } \right){\beta _{{\rm{elc}}}}, \;\;{P_{{\rm{wh}}}} < 0. \end{array} \right. $

(13)

其中: η_T为传动系统的效率，η_M为电机的效率，β为前后轴的机械制动力分配系数(β=0.4)，β_elc表示叠加在后轴上的电制动力矩与后轴机械制动力矩的比例。P_wh是车速和加速度的函数,

$ {P_{{\rm{wh}}}} = \frac{{Gf{u_{\rm{a}}}}}{{3{\rm{ }}600}} + \frac{{{C_{\rm{D}}}Au_{\rm{a}}^{^3}}}{{76{\rm{ }}140}} + \frac{{\delta m{u_{\rm{a}}}}}{{3{\rm{ }}600}}a. $

(14)

其中: G为整车重量，m为整车质量，f为滚动阻力系数，C_D为空气阻力系数，A是车辆迎风面积，δ为整车旋转质量换算系数，a为车辆加速度。

2) 将需求功率的取值范围划分为s个状态区间，每个状态区间对应的需求功率的状态值P_j(j=1, 2, …, s)用落在该状态区间的所有需求功率的平均值表示。

3) 统计需求功率序列中从状态i转移到状态j的频数f_{i, j}(i=1, 2, …, s；j=1, 2, …, s)，并在统计结束后，计算需求功率的转移概率,

$ {\pi _{i, j}} = {f_{i, j}}/\sum\limits_{j = 1}^s {{f_{i, j}}} . $

(15)

本文的仿真工况为中国城市公交工况，车速序列如图 3所示，对应的需求功率序列如图 4所示。将需求功率的取值范围分为20个状态区间，依次标为状态1, 2，…，20，落在每个状态区间的需求功率的平均值依次为：-46.53 kW，-36.75 kW，-22.68 kW，-8.92 kW，-0.39 kW，3.53 kW，13.00 kW，20.36 kW，28.36 kW，34.67 kW，41.02 kW，48.11 kW，54.65 kW，62.49 kW，72.51 kW，88.03 kW，100.61 kW，114.57 kW，127.35 kW，140.02 kW。通过统计，得到各状态间的转移概率如图 5所示。

2.3 需求功率序列预测

基于图 5的转移概率矩阵和当前需求功率的状态，利用文[21]中的极大似然估计方法，构建一条未来需求功率的预测序列，相关定义说明如下：

1) $ \mathcal{J} $={T₁, T₂, …, T_N}：预测序列，T_i表示需求功率的状态。当驾驶员踩下油门或制动踏板时，当前的需求功率P_req(k)已知，根据需求功率的状态区间划分，获取当前需求功率的状态T₁。基于需求功率的Markov预测模型，利用最大似然估计方法，构建未来第2到第N步的需求功率状态，并得到需求功率取值的预测序列，{P_{req, 1}, P_{req, 2}, …, P_{req, N}}={P_T1, P_T2, …, P_TN}，其中：P_T1=P_req(k)；P_Ti，i∈{2, 3, …, N}是基于预测状态T_i得到的需求功率的状态值。

2) pre(T_i)：T_i在预测序列$ \mathcal{J} $中的前一个节点，例如，pre(T₂)=T₁。

3) succ(T_i, j)：Markov链中，需求功率的状态T_i能够到达的下一个状态j，j∈1, 2, …, s。

4) $ \mathcal{S} $：预测序列中属于末端的节点，其定义为$ \mathcal{S} \buildrel \Delta \over = \left\{ {{T_i}|{\rm{succ}}\left( {{T_i}, j} \right) \in \mathcal{J}, j \in \left\{ {1, 2, \cdots, s} \right\}} \right\} $，显然，$ \mathcal{S} \subset \mathcal{J} $。

5) π_Ti：预测序列中，从T₁到达T_i的概率。

3 SOC消耗参考曲线设计

本文采用MPC方法对阶段性的功率分配问题进行滚动优化。在每个控制周期内，基于需求功率的预测序列，设计一条SOC消耗参考曲线{SOC_refⁱ(k)|i=1, 2, …, N}，k为当前控制周期数。以预测时域长度N=5为例，图 6给出了3个相邻控制周期内SOC消耗参考曲线的示例。考虑到每个控制周期内优化得到的N个控制量只会执行第1个，将第k个控制周期内SOC消耗参考曲线的初值SOC_ref¹(k)设定为第k-1个控制周期内SOC消耗参考曲线上第2个点的值^[23]。

车辆平均每km能量消耗e为

$ \bar e\left( k \right) = \left\{ \begin{array}{l} {{\bar e}_0}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;d\left( k \right) = 0;\\ \max \left\{ {\frac{{E\left( k \right)}}{{d\left( k \right)}}, {{\bar e}_0}} \right\}, d\left( k \right) > 0. \end{array} \right. $

(16)

其中，e₀为12 m增程式城市客车在中国城市工况下的平均每km能量消耗，kWh/km。基于总的行程d₀，预估剩余行程下车辆的能量需求E_rest为

$ {E_{{\rm{rest}}}}\left( k \right) = \bar e\left( k \right)\left( {{d_0}-d\left( k \right)} \right). $

(17)

在第k个控制周期内，基于E_rest(k)，将SOC消耗参考曲线上剩余的可用SOC按照未来需求功率的预测值进行分配，其未来第i步的增量ΔSOC_refⁱ(k)为

$ \begin{array}{l} \Delta {\rm{SOC}}_{{\rm{ref}}}^{^i}\left( k \right) = \\ \frac{{{\rm{SOC}}_{{\rm{ref}}}^{^1}\left( k \right) - {\rm{SO}}{{\rm{C}}_{{\rm{ref,end}}}}}}{{{E_{{\rm{rest}}}}\left( k \right)}}\sum\limits_{j = 1}^{i - 1} {\left( {\frac{{\Delta T}}{{3{\rm{ }}600}}{P_{{\rm{req,}}j}}} \right)} . \end{array} $

(18)

其中：P_{req, j}为当前控制周期内预测的未来第j步的需求功率值，SOC_{ref, end}为行程结束时SOC的期望值。由此可以得到单个控制周期内第i步的SOC消耗参考值为：

$ \begin{array}{l} {\rm{SOC}}_{{\rm{ref}}}^{^i}\left( k \right) = {\rm{SOC}}_{{\rm{ref}}}^{^1}\left( k \right)-\Delta {\rm{SOC}}_{{\rm{ref}}}^{^i}\left( k \right), \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 2, {\rm{ }}3, \cdots, {\rm{ }}N. \end{array} $

(19)

其中：SOC_ref¹(k)为第k个控制周期内SOC消耗参考曲线的初值，

$ {\rm{SOC}}_{{\rm{ref}}}^{^1}\left( k \right) = \left\{ \begin{array}{l} {\rm{SO}}{{\rm{C}}(1)}, \;\;\;\;\;\;\;\;\;\;k = 1;\\ {\rm{SOC}}_{{\rm{ref}}}^2\left( {k-1} \right), k > 1. \end{array} \right. $

(20)

本文中，MPC的预测时域长度N=20。图 7给出了第101、102以及103个周期内的SOC消耗参考曲线。

4 整车能量管理的分段优化建模

以预测序列{P_T1, P_T2, …, P_TN}中每个需求功率出现的概率为加权因子，对系统的单步费用进行概率累加，从而构建以总费用最低为目标的瞬时功率分配最优化模型。为了简化符号的表述，在优化模型中，所有的T_i将简化为i来表示，即用x_i，P_{req, i}，P_{APU, i}，P_{bat, i}，π_i，pre(i)来分别表示x_Ti，P_Ti，P_{APU, Ti}，P_{bat, Ti}，π_Ti，pre(T_i)。模型如下：

$ \begin{array}{l} \min \;J = \sum\limits_{i \in \mathcal{J}\backslash \mathcal{S}} {{\pi _i}\alpha } P_{{\rm{bat}},i}{\left( {{\eta _{{\rm{bat}}}}} \right)^{-{\rm{sign}}\left( {{P_{{\rm{bat}}, i}}} \right)}}\frac{{\Delta T}}{{3{\rm{ }}600}} + \\ \sum\limits_{i \in \mathcal{J}\backslash S} {{\pi _i}\beta {{\dot m}_{\rm{f}}}} \left( {{P_{{\rm{APU, }}i}}} \right){P_{{\rm{APU, }}i}}\frac{{\Delta T}}{{3{\rm{ }}600}}\frac{1}{{{\rho _{{\rm{fuel}}}}}} + \\ \sum\limits_{i \in \mathcal{J}\backslash \{ {T_1}\} } {{\pi _i}\gamma |{\rm{SOC}}{_i}-{{\rm{SOC}}_{{\rm{ref,}}k}^i}|} . \end{array} $

(21)

系统约束：

$ \begin{array}{l} {\mathit{\boldsymbol{x}}_1} = \mathit{\boldsymbol{x}}\left( k \right), \\ {P_{{\rm{req}}, 1}} = {P_{{\rm{req}}}}\left( k \right), \\ {u_{{\rm{a}}, 1}} = {u_{\rm{a}}}\left( k \right), \\ {\mathit{\boldsymbol{x}}_i} = {\left[{{\rm{SO}}{{\rm{C}}_i}\;\;{d_i}\;\;{E_i}} \right]^{\rm{T}}}.\\ {\rm{SO}}{{\rm{C}}_i} = {\rm{SO}}{{\rm{C}}_{{\rm{pre}}\left( i \right)}} -\frac{{{I_{{\rm{bat, pre}}\left( i \right)}}}}{{{Q_{{\rm{bat}}}}}}\frac{{\Delta T}}{{3{\rm{ }}600}}, \\ {d_i} = {d_{{\rm{pre}}\left( i \right)}} + {u_{{\rm{a}}, {\rm{pre}}(i)}}\frac{{\Delta T}}{{3{\rm{ }}600}}, \\ {d_i} = {d_{{\rm{pre}}\left( i \right)}} + {u_{{\rm{a}}, {\rm{pre}}(i)}}\frac{{\Delta T}}{{3{\rm{ }}600}}, \\ {E_i} = {E_{{\rm{pre}}\left( i \right)}} + {P_{{\rm{req}}, {\rm{pre}}(i)}}\frac{{\Delta T}}{{3{\rm{ }}600}}, \\ {P_{{\rm{bat}}, i}} = {P_{{\rm{req}}, i}} -{P_{{\rm{APU}}, i}}, \;\forall i \in \mathcal{J}\backslash \mathcal{S}, \\ {\mathit{\boldsymbol{x}}_i} \in \mathit{\boldsymbol{X}}, \forall i \in \mathcal{J}\backslash \{ {T_1}\}, \\ {P_{{\rm{APU}}, i}} \in \mathit{\boldsymbol{U}}, \forall i \in \mathcal{J}\backslash \mathcal{S}, \\ {P_{{\rm{bat}}, i}} \in \mathit{\boldsymbol{Y}},\forall i \in \mathcal{J}\backslash \mathcal{S}. \end{array} $

其中：α为充电价格，α=1.194 6元/(kWh)；β为柴油价格，β=5.67元/L；γ为SOC惩罚因子，γ=1 000；ρ_fuel为燃油密度；SOC_{ref, k}ⁱ为SOC消耗参考值SOC_refⁱ(k)。目标函数主要由3部分组成：电耗产生的费用、油耗产生的费用以及约束SOC跟随其消耗参考曲线的罚函数项。状态量中位移d_i和能耗E_i用于构建SOC消耗参考曲线，在第k个控制周期内的N步预测控制中，仅需获取d₁和E₁的值，可以分别基于d(k-1)、u_a(k-1) 和E(k-1)、P_req(k-1) 来计算。

5 仿真计算和结果分析

利用模型预测控制方法对由式(21) 给定的优化问题进行滚动优化，预测时域长度为N=20。在每一个控制周期内，利用DDP算法来求解式(21) 所示的优化问题。为了降低求解的难度，本文对上述优化问题进行了如下简化：当APU需要升功率时，只允许其向相邻高功率点调节；当APU需要降功率时，只允许其向相邻低功率点调节或者进入怠速状态，以此达到限制APU动态调节的幅度来降低发动机动态调节过程的油耗，同时减少DDP中遍历状态的个数以降低算法的复杂度。针对由中国城市公交工况构造的200 km循环工况，计算得到了12 m增程式城市客车的纯油耗为21.8 L/(100 km)，电耗为25.4 kWh/(100 km)，仿真结束时，动力电池组SOC值为0.250 9。图 8为前1 500 s的需求功率和优化得到的APU输出功率。基于改进的DDP算法，APU的输出功率没有出现频繁的调节。图 9中展示了此过程中动力电池组SOC曲线的变化趋势，电池组SOC曲线较好地跟随了SOC消耗参考曲线，虚线由每个控制周期内SOC消耗参考曲线上第1个点的值组成，即

$ \left\{ {{\rm{SOC}}_{_{{\rm{ref}}}}^{^1}\left( k \right)|k = 1, {\rm{ }}2, \cdots, 1{\rm{ }}500} \right\}. $

同样条件下，采用文[10]中的CDCS方法，仿真计算结果为纯油耗为24.1 L/(100 km)，电耗为25.4 kWh/(100 km)，可见本文的方法提高了燃油经济性。

6 结论

针对车辆在行驶过程中未来需求功率未知的情况，为了既保证增程式城市客车的动力性又提高燃油经济性，本文提出了一种基于电池SOC消耗管理和功率分配的整车能量分段跟踪优化方法。该方法通过在每个控制周期内构建一个需求功率预测序列和设计一条SOC消耗参考曲线，将全程的能量优化问题转化为一系列的阶段性优化问题，然后利用模型预测控制方法对功率分配策略进行滚动优化，并在中国城市公交工况下对该方法的燃油经济性进行了仿真验证，得到结论如下：

1) 利用模型预测控制方法对功率分配策略进行滚动优化，可以减少需求功率预测误差的影响，能够使得实际SOC曲线较好地跟随SOC消耗参考曲线。

2) 通过防止电池SOC在行程中被过早耗尽，并使其在行程结束时降到最低，能够在保证动力性的前提下提高增程式城市客车的燃油经济性。

为进一步改善车辆燃油经济性，未来将研究利用工况辨识技术来提高需求功率预测准确度的方法。