道路交通运输的能源消耗和温室气体排放会对环境产生负面影响^[1]. 车辆编队的各成员车辆以队列形式在较短的车间距下稳定行驶，能够通过减少成员间空气阻力来降低能源消耗和碳排放. 乘客舒适度是评估队列控制性能的指标，它是指乘客（包括驾驶员）在车辆行驶过程中的舒适程度. 车辆剧烈变速可能对驾驶安全和乘客健康产生不良影响，因此乘客舒适度分析的重要性日益突出^[2]. 乘客舒适度、能耗以及跟踪性能在控制方案中存在耦合关系，这给队列控制带来了新的难点^[3]. 随着网联自动驾驶车辆(connected autonomous vehicle, CAV)市场渗透率的持续增长以及车路协同基础设施的推广部署，道路交通即将长期处于人工驾驶车辆(human-driven vehicle, HDV)和CAV共存的局面^[4]，不同驾驶模式的跟驰特性差异将加剧队列控制的复杂性. 如何在考虑乘客舒适度的前提下，设计面向CAV与HDV混合队列的协作控制方法以优化队列能耗是待解决的重要问题.

按照混合车辆队列控制目标是否促使各成员合作将已有研究工作一分为二^[5]，其中合作控制目标包括队列稳定性与跟踪性能，非合作目标主要包括能耗与乘客舒适度. 1）在研究混合车辆队列稳定性方面，边有钢等^[6]提出适用于双向多车跟随式拓扑的混合队列控制器，将前后车信息引入CAV的控制设计，分析前后车信息对混合队列稳定性的影响. 李淑庆等^[7]针对由协作式巡航控制(cooperative adaptive cruise control，CACC)车辆、自适应巡航控制(adaptive cruise control，ACC)车辆以及HDV组成的混合车队，考虑车辆之间的退化机制，应用传递函数理论，推导队列稳定性判别准则，分析不同队列规模下的稳定性. Zhan等^[8]将数据驱动方法与模型预测控制(model predictive control，MPC)相结合，分别设计集中式MPC和分布式MPC控制器，以消除停起波动. 2）为了优化混合车辆队列的能耗，Hu等^[9]提出灵活的间距策略，确保所有CAV编队成员都保持节能驾驶模式. Yao等^[10]通过引入燃料消耗和交通排放模型，研究CAV渗透率对混合交通流的燃油消耗和排放的影响. Zhao等^[11]提出递进视野的模型预测控制方法以降低能耗，还将该方法扩展到动态编队分割和合并规则. 3）车辆行驶平稳性将直接影响乘客舒适度，车辆的加速度过大会严重影响编队的乘客舒适度^[12]. 秦严严等^[13]分析CACC车辆、ACC车辆与HDV混合交通流的驾驶舒适性. Li等^[14]提出考虑乘客舒适度的CACC系统，在保证混合交通稳定性的前提下，提高驾驶员的身心舒适度. Tian等^[15]使用智能驾驶员模型(intelligent driver model，IDM)作为HDV的跟驰模型，将加速度导数作为控制变量，采用加速度描述车辆能耗，将加速度导数和加速度分别作为加权项，设计MPC控制律，在降低能耗的同时提高驾驶舒适性. 文献[13]、[14]未涉及能耗优化；文献[15]涉及的车辆能耗建模参量单一，描述混合队列车辆实时能耗的准确性欠佳.

本研究在队列稳定控制的基础上，以乘客舒适度为约束，以降低能耗为优化目标，提出混合车辆队列协作控制新算法. 1）针对队列中的CAV，建立乘客舒适度约束，采用实时优化(real-time optimization，RTO)的分布式模型预测控制(distributed model predictive control，DMPC)方法（RTODMPC）在确保车辆稳定协作的基础上降低实时能耗. 2）针对队列中的HDV，采用乘客舒适度和能耗性能良好的智能驾驶员跟驰模型，考虑前车扰动推理得到稳定跟驰条件. 3）引入贴合实际的基于动力的综合油耗模型(Virginia tech comprehensive power based fuel consumption model, VT-CPFM)，通过对瞬时油耗的计算提升能耗优化性能. 4）对所提算法进行仿真分析，验证算法的跟踪性能，将从初始状态到稳态的平均发动机功率作为能耗性能指标，对比不同控制方案的能耗优化性能.

由$ n+q+1 $辆车组成的混合车辆队列在平直道路上行驶，领航车辆为CAV，记为$ 0 $号车. 跟随车辆序列记为集合$ {Z} $，其中前$ n $辆车为CAV，记为集合$ {{C}_n} $，随后的$ q $辆为HDV，记为集合$ {H} $，有$ {Z} = {{C}_n} \cup {H} $. 将队列中所有CAV记为集合$ {{C}_{n+1}} $，因此，$ {{C}_{n+1}} = \{ 0\} \cup {{C}_n} $. 车辆编号$ i $为车辆实际序号. 如图1所示，CAV之间通过车车通信(vehicle to vehicle, V2V)技术获得信息.

混合车辆队列中的第$ i $辆CAV（$ i \in {{{C}}_{n+1}} $）在$ t $时刻的非线性纵向车辆动力学模型^[16]为

(1)$ \left. \begin{gathered} {{\dot p}_i}(t) = {v_i}(t), \\ {m_i}{a_i}(t) = \frac{{{T_i}(t)}}{{{r_{{\text{W}},i}}}} - \frac{1}{2}\rho {A_{\text{f}}}C_i^{\text{A}}v_i^2(t) - {m_i}gf_i^{\text{R}}, \\ {\tau _i}{{\dot T}_i}(t)+{T_i}(t) = {T_{{\text{des}},i}}(t). \\ \end{gathered} \right\} $

式中：$ {p_i}(t) $、$ {v_i}(t) $和$ {a_i}(t) $分别为车辆$ i $在$ t $时刻的位置、速度和加速度，$ {T_i} $和$ {T_{{\text{des}},i}} $分别为车辆$ i $车轮上的实际和期望驱动/制动扭矩，$ {m_i} $和$ {r_{{\text{W}},i}} $分别为车辆$ i $的质量和轮胎半径；$ {\tau _i} $为车辆$ i $的传动系统时滞，$ \rho $为空气密度，$ {A_{\text{f}}} $为车辆正向截面面积，$ C_i^{\text{A}} $为空气阻力系数，$ f_i^{\text{R}} $和$ g $分别为滚动阻力和重力加速度系数. 定义$ \Delta t $为采样间隔，$ {u_i}(t) $为车辆$ i $的控制输入，式(1)经变换、离散化的表达式^[17]为

(2)$ \left. \begin{gathered} {p_i}(t+1) = {p_i}(t)+{v_i}(t) \cdot \Delta t, \\ {v_i}(t+1) = {v_i}(t)+{a_i}(t) \cdot \Delta t, \\ {a_i}(t+1) = {a_i}(t)+{u_i}(t) \cdot \Delta t. \\ \end{gathered} \right\} $

(3)$ \begin{split} {u_i}(t) =& \frac{1}{{{\tau _i}}}\left\{ {\frac{1}{{{m_i}}}\left[ {\frac{{{T_{{\text{des}},i}}(t)}}{{{r_{{\mathrm{W}},i}}}} - {m_i}gf_i^{\text{R}} - } \right.} \right. \\& \left. {\left. {\frac{1}{2}\rho {A_{\text{f}}}C_i^{\text{A}}{v_i}(t)(2{\tau _i}{a_i}(t)+{v_i}(t))} \right] - {a_i}(t)} \right\}. \end{split} $

令车辆$ i $在时刻$ t $的状态向量$ {{\boldsymbol{x}}_i}(t) = {[{p_i}(t),\;\;{v_i}(t),\;\;{a_i}(t)]^{\text{T}}} $，整理式(2)为

(4)$ \begin{aligned} \left.\begin{array}{*{20}{l}}{{\boldsymbol{x}}_i}(t+1) = {\boldsymbol{A}}{{\boldsymbol{x}}_i}(t)+{\boldsymbol{B}}{u_i}(t);\\{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} 1&{\Delta t}&0 \\ 0&1&{\Delta t} \\ 0&0&1 \end{array}} \right],\qquad {\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ {\Delta t} \end{array}} \right]. \end{array}\right\}\end{aligned}$

采用图论模型，将混合车辆队列跟随序列中的CAV间的通信拓扑记为有向图$ {{G}_n} = \{ {{C}_n},{\text{ }}{{\text{E}}_n}\} $，将图$ {{G}_n} $的邻接矩阵记为$ {\boldsymbol{S}} = [{s_{ij}}] \in {{{\bf{R}}}^{n \times n}} $，如果跟随车辆$ i $可以获取跟随车辆$ j $的信息，则对应$ {s_{ij}} = 1 $，否则$ {s_{ij}} = 0 $. 假设有向图$ {{G}_n} $不存在自环，即$ {s_{ii}} = 0 $；车辆$ i $的内邻集$ {{N}_i} $定义为车辆$ i $所有入边邻接点组成的集合，记为$ {{N}_i} = \{ j|j \in {{C}_n},{\text{ }}j \ne i,{\text{ }}{s_{ij}} = 1\} $，车辆$ i $的外邻集$ {{O}_i} $定义为车辆$ i $所有出边邻接点组成的集合，记为$ {{O}_i} = \{ j|j \in {{C}_n},{\text{ }}j \ne i,{\text{ }}{s_{ji}} = 1\} $；$ {{G}_n} $的入度矩阵记为$ {\boldsymbol{D}} = {\text{diag}}\;({d_{ii}}) \in {{{\bf{R}}}^{n \times n}} $，其中$ {d_{ii}} = \displaystyle\sum\nolimits_{j \in {{N}_i}} {{s_{ij}}} $；图$ {{G}_n} $的拉普拉斯矩阵记为$ {\boldsymbol{Y}} = [{y_{ij}}] \in {{{\bf{R}}}^{n \times n}} $，即$ {\boldsymbol{Y}} = {\boldsymbol{D}} - {\boldsymbol{S}} $. 将包含领航车辆在内的CAV集合$ {{C}_{n+1}} $中车辆间的通信拓扑记为有向图$ {{G}_{n+1}} = \{ {{C}_{n+1}},{\text{ }}{{{E}}_{n+1}}\} $. 将领航车辆与其他$ n $辆跟随车辆间的通信连接记为连接矩阵$ {\boldsymbol{W}} = {\text{diag}}\;({w_{ii}}) \in {{{\bf{R}}}^{n \times n}} $，如果跟随车辆$ i $可以获得领航车辆的信息，则$ {w_{ii}} = 1 $，否则$ {w_{ii}} = 0 $；考虑领航车辆，将跟随车辆$ i $的内邻集记为$ {{I}_i} = \{ j|j \in {{C}_{n+1}},{\text{ }}j \ne i,{\text{ }}{s_{ij}} = 1\} $. CAV队列不同通信节点连接决定了CAV队列中不同的通信拓扑结构. CAV队列有4种典型的通信拓扑：前驱跟随(predecessor-following, PF)、前驱领航跟随(predecessor-leader-following, PLF)、两前驱跟随(two-predecessor-following, TPF)和两前驱领航跟随(two-predecessor- leader-following, TPLF). 假设CAV车间通信没有通信延迟和丢包. 在本研究中，HDV没有配备车间通信设备，因此不存在V2V链路，不考虑HDV通信拓扑.

VT-CPFM将油耗表征为车辆功率的二阶多项式函数，能够精确计算车辆的实时油耗^[18]. 采用VT-CPFM表示车辆的综合油耗. 车辆$ i $的燃油消耗$ {F_{{\text{C}},i}}(t) $的表达式为

(5)$ {F_{{\text{C}},i}}(t) = \left\{ \begin{array}{ll} {\alpha _{0,i}}+{\alpha _{1,i}}{P_i}(t)+{\alpha _{2,i}}{P_i}{(t)^2}{\text{ }},&{P_i}(t) \geqslant 0; \\ {\alpha _{0,i}}{\text{ , }}&{P_i}(t) < 0. \\ \end{array} \right. $

式中：$ {\alpha }_{0,i},\text{ }{\alpha }_{1,i},\text{ }{\alpha }_{2,i} $为车辆$ i $的VT-CPFM系数，$ {P_i}(t) $为车辆$ i $在时刻$ t $的瞬时功率，

(6)$ {P_i}(t) = \frac{{{T_i}(t){v_i}(t)}}{{{\eta _{{\text{T}},i}}{r_{{\text{W}},i}}}} = \dfrac{{\left[ {{m_i}{a_i}(t)+\dfrac{1}{2}\rho {A_{\text{f}}}C_i^{\text{A}}v_i^2(t)+{m_i}gf_i^{\text{R}}} \right]{v_i}(t)}}{{{\eta _{{\text{T}},i}}}}. $

其中$ {\eta _{{\text{T}},i}} $为车辆$ i $传动系统的机械效率.

跟踪性能控制目标要求队列中车辆与前车保持相同的车速，队列车辆之间还须保持期望的车间距. 对于CAV和HDV，车辆与前车的跟踪性能控制目标为

(7)$ \left. \begin{gathered} \mathop {\lim }\limits_{t \to +\infty } \left[ {{p_i}(t) - {p_{i - 1}}(t)+{d_0}} \right] = 0, \\ \mathop {\lim }\limits_{t \to +\infty } \left[ {{v_i}(t) - {v_{i - 1}}(t)} \right] = 0, \\ \mathop {\lim }\limits_{t \to +\infty } \left[ {{a_i}(t) - {a_{i - 1}}(t)} \right] = 0. \\ \end{gathered} \right\} $

式中：$ {d_0} $为期望的车间距. 对于CAV中的车辆$ i $与车辆$ j $，式(7)可写为

(8)$ \mathop {\lim }\limits_{t \to +\infty } \left\| {{{\boldsymbol{x}}_i}(t) - {{\boldsymbol{x}}_j}(t) - {{{\boldsymbol{\tilde d}}}_{ji}}} \right\| = 0. $

其中$ {{\boldsymbol{\tilde d}}_{ji}} = {[(j - i){d_0},\;\;0,\;\;0]^{\text{T}}} $.

为了在实现混合队列CAV与HDV协作的基础上降低能耗，同时保障乘客舒适度，本研究提出实时优化的分布式模型预测控制与智能驾驶员模型结合的协作控制方法（RTODMPC-IDM）. 1）针对CAV队列控制，DMPC能够在最优控制的框架下有效处理多约束和多目标优化问题，RTO能够平衡多目标间的冲突. RTODMPC为分层架构，由协作一致控制层（上层）和能耗优化控制层（下层）构成；上层保证车辆跟踪性能，下层降低车辆能耗. 在RTODMPC中，所有CAV通过V2V传输共享自身在预测范围内的假设轨迹信息，每辆 CAV 结合自车状态按序求解上层与下层优化问题. 通过引入控制输入和加速度约束保障乘客舒适度，采用VT-CPFM构建下层目标函数，以更充分降低能耗. 2）针对HDV，选择IDM. 原因是IDM跟驰模型具有良好的乘客舒适性，同时保持低能耗^[19]. 在IDM减速操作中，HDV在达到稳定的跟驰状态前或者到达停止线前，减速度逐渐增大至舒适值，然后平滑减小到0；加速操作被最大加速度约束. 推导前车扰动下跟驰行为的线性稳定性条件，在此基础上，每辆HDV通过IDM根据自车与前车的相对距离、相对速度及自车速度，确定自车加速度并保障乘客舒适度.

RTODMPC的上下层具有相同的时间尺度，能够实现双层实时求解. 考虑当前时刻$ t $，令$ {N_{\text{p}}} $为预测时域长度($ {N_{\text{p}}} > 0 $)，在预测时域$ [t,{\text{ }}t+{N_{\text{p}}} - 1] $内，定义CAV_i，$i \in {{{C}}_{n+1}} $在$ t $时刻的预测状态和假设状态分别为$ {\boldsymbol{x}}_i^{\text{p}}(t) $和$ {\boldsymbol{x}}_i^{\text{a}}(t) $. 将车辆$ i $在$ t+k $时刻预测状态及其控制输入分别记为$ {\boldsymbol{x}}_i^{\text{p}}(k|t) $和$ u_i^{\text{p}}(k|t) $，车辆$ i $在$ t+k $时刻的假设状态及其控制输入分别记为$ {\boldsymbol{x}}_i^{\text{a}}(k|t) $和$ u_i^{\text{a}}(k|t) $，其中$ k \in [0,{\text{ }}{N_{\text{p}}} - 1] $. 矩阵加权范数$ {\left\| {{\boldsymbol{x}}} \right\|_{P}} = \sqrt {{{{\boldsymbol{x}}}^{\rm{T}}}{{\boldsymbol{Px}}}} $，其中$ {{\boldsymbol{P}}} $为对称正定矩阵；$ \left| I_i \right| $为集合$ I_i$中元素个数.

考虑式(4)描述的CAV动力学系统，车辆$ i $根据自身及其内邻集$ {I_i} $中车辆的假设状态序列、预测状态序列及其控制输入序列，构建最优控制问题：

(9)$ \mathop {\min }\limits_{u_i^{\text{p}}(t)} {J_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(t),{\text{ }}u_i^{\text{p}}(t)} \right]. $

(10)$ \begin{split}& {J_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(t),{\text{ }}u_i^{\text{p}}(t)} \right] = \\ & \qquad \sum\limits_{k = 0}^{{\text{ }}{N_{\text{p}}} - 1} {{L_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(k|t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(k|t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(k|t),{\text{ }}u_i^{\text{p}}(k|t)} \right]} . \end{split}$

(11)$ \begin{split}& {L_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(k|t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(k|t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(k|t),{\text{ }}u_i^{\text{p}}(k|t)} \right] = \\&\qquad {\text{ }}{\left\| {{\boldsymbol{x}}_i^{\text{p}}(k|t) - {\boldsymbol{x}}_i^{\text{a}}(k|t)} \right\|_{{{\boldsymbol{F}}_i}}}+{\left\| {u_i^{\text{p}}(k|t)} \right\|_{{{\boldsymbol{R}}_i}}}+ \\& \qquad{\text{ }}\sum\limits_{j \in {{{I}}_i}} {{{\left\| {{\boldsymbol{x}}_i^{\text{p}}(k|t) - {\boldsymbol{x}}_j^{\text{a}}(k|t) - {{{\boldsymbol{\tilde d}}}_{ji}}} \right\|}_{{{\boldsymbol{G}}_i}}}} . \end{split} $

式(10)为协作一致控制层最优控制问题的目标函数，其中$ {{\boldsymbol{R}}_i} \in {{\bf{R}}}{\text{, }}{{\boldsymbol{F}}_i} \in {{{\bf{R}}}^{3 \times 3}},{\text{ }}{{\boldsymbol{G}}_i} \in {{{\bf{R}}}^{3 \times 3}} $是正定权重矩阵. 协作一致控制层最优控制问题的约束为

(12)$ {\boldsymbol{x}}_i^{\text{p}}(0|t) = {{\boldsymbol{x}}_i}(t), $

(13)$ {\boldsymbol{x}}_i^{\text{p}}(k+1|t) = {\boldsymbol{Ax}}_i^{\text{p}}(k|t)+{\boldsymbol{B}}u_i^{\text{p}}(k|t), $

(14)$ {\boldsymbol{x}}_i^{\text{p}}({N_{\text{p}}}|t) = {\boldsymbol{x}}_i^{\text{a}}({N_{\text{p}}}|t), $

(15)$ {\boldsymbol{x}}_i^{\text{a}}({N_{\text{p}}}|t) = {\boldsymbol{Ax}}_i^{\text{a}}({N_{\text{p}}}|t - 1)+{\boldsymbol{B}}u_i^ * ({N_{\text{p}}}|t - 1). $

式(12)为初始条件，式(13)为系统约束，式(14)保证预测状态和假设状态终端相等，式(15)为假设状态终端的迭代. 为了保证乘客舒适度，分别对加速度及其变化率进行约束：设置加速度约束为$ {a_{{\text{LB}}}} < {a_i} < {a_{{\text{UB}}}} $，其中$ {a_{{\text{LB}}}} $和$ {a_{{\text{UB}}}} $为满足乘客舒适度要求的加速度下限和上限，控制输入反映加速度变化率，设置控制输入约束为$ - 1 \leqslant {u_i} \leqslant 1 $. 式(15)中的$ u_i^ * ({N_{\text{p}}}|t - 1) $的表达式^[20]为

(16)$ \begin{split} u_i^ * ({N_{\text{p}}}|t - 1) = \frac{1}{{\left| {{{I}_i}} \right|}}{\boldsymbol{K}}\sum\limits_{j \in {{{I}}_i}} {\left[ {{\boldsymbol{x}}_j^{\text{a}}({N_{\text{p}}}|t - 1)+{{{\boldsymbol{\tilde d}}}_{ji}} - {\boldsymbol{x}}_i^{\text{a}}({N_{\text{p}}}|t - 1)} \right]} , \end{split} $

(17)$ {\boldsymbol{K}} = - {({{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{PB}}+{{\boldsymbol{I}}_3})^{ - 1}}{{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{PA}}. $

式(17)中$ {\boldsymbol{P}} > {\bf{0}} $，是矩阵代数Riccati方程（matrix algebraic Riccati equation, MARE）的解：

(18)$ {\boldsymbol{P}} = {{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{PA}}+{\boldsymbol{Q}} - (1 - {\delta ^2}){{\boldsymbol{A}}^{\text{T}}}{\boldsymbol{PB}}{({{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{PB}}+{{\boldsymbol{I}}_3})^{ - 1}}{{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{PA}}. $

式中：$ {\boldsymbol{Q}} > {\bf{0}} $为常数矩阵，参数$ \delta $满足：$ \mathop {\max }\limits_{i \in {\mathcal{C}_n}} \left| {{\lambda _i}({{\boldsymbol{D}}^{ - 1}}{\boldsymbol{S}})} \right| < \delta < 1 $. 求解式(9)，得到本层的最优控制输入：

(19)$ u_i^0(t) = \mathop {\arg \min }\limits_{u_i^{\text{p}}(t)} {J_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(t),{\text{ }}u_i^{\text{p}}(t)} \right]. $

本层最优状态通过车辆状态方程式(4)和最优控制输入式(19)完成更新迭代：

(20)$ {\boldsymbol{x}}_i^0(t) = {\boldsymbol{Ax}}_i^0(t - 1)+{\boldsymbol{B}}u_i^0(t - 1). $

根据式(19)、(20)，得到目标函数的最优值:

(21)$ J_{{\text{c}},i}^0(t) = {J_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^0(t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(t),{\text{ }}u_i^0(t)} \right]. $

上层的最优成本构建下层稳定性的上限约束，在保证跟踪性能和稳定性的基础上降低能耗. 引入VT-CPFM构建能耗优化控制层的目标函数，设计最优化问题式为

(22)$ \mathop {\min }\limits_{u_i^{\text{p}}(t)} {J_{{\text{e}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}u_i^{\text{p}}(t)} \right]{\text{ }},{\text{ }}i \in {{C}_{n+1}}. $

(23)$ {J_{{\text{e}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}u_i^{\text{p}}(t)} \right]{\text{ = }}\sum\limits_{k = 0}^{{N_{\text{p}}} - 1} {{L_{{\text{e}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(k|t),{\text{ }}u_i^{\text{p}}(k|t)} \right]} , $

(24)$ {L_{{\text{e}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(k|t),{\text{ }}u_i^{\text{p}}(k|t)} \right] = {F_{{\text{C}},i}}(k|t) \cdot \Delta t. $

式(23)为能耗优化控制层最优控制问题的目标函数，式(24)为车辆$ i $的瞬时能耗目标函数，$ {F_{{\text{C,i}}}}(t) $为VT-CPFM中车辆的燃油消耗率. Bian等^[20]证明了多层分布式模型预测控制下CAV队列的稳定性，给出下层分布式模型预测控制须遵循上层的迭代可行性及系统稳定性保证约束：

(25)$ {J_{{\text{c}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(t),{\text{ }}{\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(t),{\text{ }}u_i^{\text{p}}(t)} \right] \leqslant J_{{\text{c}},i}^0(t)+\beta {\varepsilon _i}(t - 1){\text{.}} $

(26)$ \begin{split}& \beta \in (0,{\text{ }}1) \text{，} {\varepsilon _i}(t - 1) = {L_{{\text{c}},i}}\left[ {\boldsymbol{x}}_i^ * (0|t - 1),{\text{ }}{\boldsymbol{x}}_i^{\text{a}}(0|t - 1),\right. \\ &\qquad\left. {\boldsymbol{x}}_{j \in {{{I}}_i}}^{\text{a}}(0|t - 1),{\text{ }}u_i^ * (0|t - 1) \right]+ \\& \qquad\sum\limits_{k = 1}^{{N_{\text{p}}} - 1} \left[ {{{\left\| {{\boldsymbol{x}}_i^ * (k|t - 1) - {\boldsymbol{x}}_i^{\text{a}}(k|t - 1)} \right\|}_{{{\boldsymbol{F}}_i}}} - } \right. \\ &\qquad \sum\limits_{j \in {{{O}}_i}} \left. {\left\| {{\boldsymbol{x}}_i^ * (k|t - 1) - {\boldsymbol{x}}_i^{\text{a}}(k|t - 1)} \right\|}_{{{\boldsymbol{G}}_j}} \right]. \end{split} $

为了保证上层优化的可行解对下层优化，对于优化问题式(22)，其约束包括协作一致控制层的所有约束，即式(12)~(15)和乘客舒适度约束；优化问题式(22)还须满足迭代可行性及系统稳定性保证约束，即式(25). 结合目标函数及约束条件，求得下层的最优控制输入：

(27)$ u_i^ * (t) = \mathop {\arg \min }\limits_{u_i^{\text{p}}(t)} {J_{{\text{e}},i}}\left[ {{\boldsymbol{x}}_i^{\text{p}}(t),{\text{ }}u_i^{\text{p}}(t)} \right]. $

通过车辆状态方程式(4)及$ u_i^ * (t) $求得能耗优化控制层最优状态：

(28)$ {\boldsymbol{x}}_i^ * (t) = {\boldsymbol{Ax}}_i^ * (t - 1)+{\boldsymbol{B}}u_i^ * (t - 1). $

将最优控制输入$ u_i^ * (t) $用于更新假设状态的控制输入，表达式为

(29)$ u_i^{\text{a}}(k|t+1) = \left\{ \begin{array}{ll} u_i^ * (k+1|t),&k \in [0,{N_{\text{p}}} - 2]; \\ u_i^ * ({N_{\text{p}}}|t){\text{ }},& k = {N_{\text{p}}} - 1. \\ \end{array} \right. $

其中$ u_i^{\text{a}}(k|0) \in [ - 1,1] $. 实时优化的分布式模型预测控制中的假设状态$ {\boldsymbol{x}}_i^{\text{a}} $通过式(28)和车辆状态方程式(4)完成迭代：

(30)$ {\boldsymbol{x}}_i^{\text{a}}(k|t+1) = \left\{ \begin{array}{ll} {\boldsymbol{x}}_i^ * (k+1|t){\text{ , }}&k \in [0,{N_{\text{p}}} - 1] ; \\ {\boldsymbol{Ax}}_i^{\text{a}}({N_{\text{p}}} - 1|t+1)+ & \\\qquad {\boldsymbol{B}}u_i^a({N_{\text{p}}} - 1|t+1){\text{ }},&k = {N_{\text{p}}}. \\ \end{array} \right. $

其中$ - 1 \leqslant u_i^a(k|t) \leqslant 1 $，初始时刻的假设状态为

(31)$ {\boldsymbol{x}}_i^{\text{a}}(k|0) = \left\{ \begin{array}{ll} {{\boldsymbol{x}}_i}(0){\text{ , }}&k = 0; \\ {\boldsymbol{Ax}}_i^{\text{a}}(k - 1|0)+ {\boldsymbol{B}}u_i^{\text{a}}(k - 1|0){\text{ }},&k \in [1,{N_{\text{p}}}] . \\ \end{array} \right. $

在下层最优控制输入序列中，$ u_i^ * (0|t) $用于$ t $时刻车辆$ i $的控制，$ {\boldsymbol{x}}_i^ * (k|t) $，$ k \in \left[ {1{\text{ }},{N_{\text{p}}}} \right] $和$ u_i^ * (k|t) $，$ k \in \left[ {1{\text{ }},{N_{\text{p}}}} \right] $用以生成在$ t+1 $时刻的假设轨迹.

根据IDM跟驰模型，将HDV_i，$ i \in {{H}} $的纵向动力学描述为

(32)$ \left. \begin{gathered} {p_i}(t+\Delta t) = {p_i}(t)+{v_i}(t) \cdot \Delta t+\frac{1}{2} \cdot {a_i}(t) \cdot {(\Delta t)^2}, \\ {v_i}(t+\Delta t) = {v_i}(t)+{a_i}(t) \cdot \Delta t, \\ {a_i}(t) = {f^{{\mathrm{CF}}}}\left[ {{s_i}(t),{\text{ }}\Delta {v_i}(t),{\text{ }}{v_i}(t)} \right]. \\ \end{gathered} \right\} $

(33)$ {s_i}(t) = {p_{i - 1}}(t) - {p_i}(t),\quad \Delta {v_i}(t) = {v_{i - 1}}(t) - {v_i}(t). $

式中：$ {p_i}(t) $、$ {v_i}(t) $和$ {a_i}(t) $分别为HDV_i在$ t $时刻的位置、速度和加速度. 加速度$ {a_i}(t) $由相对距离$ {s_i}(t) $、相对速度$ \Delta {v_i}(t) $和HDV速度$ {v_i}(t) $确定. 在平衡交通状态下，$ {a_i}(t) = 0 $，每辆车以相同的平衡速度$ {v^ * } $和相应的平衡间距$ {s^ * } $移动，在本研究中$ {s^ * } = {d_0} $. 平衡状态$ ({s^ * },{\text{ }}{v^ * }) $应满足：

(34)$ {f^{{\text{CF}}}}({s^ * },\; 0,\; {v^ * }) = 0. $

采用IDM车辆跟驰模型来描述HDV的加速度. HDV_i，$i \in {{H}} $的$ {a_i}(t) $为

(35)$ {a_i}(t) = {a_{\max }} \cdot \left\{ {1 - {{\left[ {\frac{{{v_i}(t)}}{{{v^0}}}} \right]}^4} - {{\left[ {\frac{{{s^ * }({v_i}(t),{\text{ }}\Delta {v_i}(t))}}{{{p_{i - 1}}(t) - {p_i}(t) - l}}} \right]}^2}} \right\}, $

(36)$ {s^ * }\left[ {{v_i}(t),{\text{ }}\Delta {v_i}(t)} \right] = {s^0}+T{v_i}(t)+\frac{{{v_i}(t) \cdot {{\left[ {{v_{i - 1}}(t) - {v_i}(t)} \right]}^2}}}{{2\sqrt {{a_{\max }}{b_{{\text{comf}}}}} }}. $

式中：$ {a_{\max }} $为最大加速度；$ {b_{{\text{comf}}}} $为舒适减速度，可用于保障乘客舒适性；$ {v^0} $为自由车流时的期望速度；$ T $为安全时间间隔；$ {s^ * }({v_i},{\text{ }}\Delta {v_i}) $、$ {s^0} $分别为期望的安全距离和静止安全距离. 忽略车辆长度不失一般性，记车辆长度$ l = 0 $.

在CAV满足控制系统稳定性的基础上，考虑混合车辆队列中CAV尾车对HDV跟驰稳定性的影响，分析得到跟驰稳定性条件. 假设CAV尾车$ n $对HDV_n+1带来扰动，使得车辆$ n+1 $在$ t $时刻偏离稳态位置$ p_{n+1}^ * (t) $. 将偏离记作$ {\xi _{n+1}}(t) $，

(37)$ {\xi _{n+1}}(t) = {p_{n+1}}(t) - p_{n+1}^ * (t). $

式中：$ {p_{n+1}}(t) $为车辆$ n+1 $在$ t $时刻受扰动影响的实际位置. 记物理延迟时延为$ T' $，由反应时间和机械延迟组成. 对式(37)两边二阶求导，得到

(38)$ {\ddot \xi _{n+1}}(t+T') = a(t+T'). $

结合式(32)，将式(38)线性化，得到

(39)$ \begin{split} {{\ddot \xi }_{n+1}}(t) =& f_{n+1}^{{\text{CF}}1}\left[ {{\xi _n}(t) - {\xi _{n+1}}(t)} \right]+ \\& f_{n+1}^{{\text{CF2}}}{{\ddot \xi }_{n+1}}(t)+{\text{ }}f_{n+1}^{{\text{CF}}3}\left[ {{{\dot \xi }_{n+1}}(t) - {{\dot \xi }_n}(t)} \right]. \end{split} $

式中：$ f_{n+1}^{{\text{CF1}}} $、$ f_{n+1}^{{\text{CF2}}} $、$ f_{n+1}^{{\text{CF3}}} $分别为跟驰模型对车间距、速度和速度差的偏导数. 根据IDM跟驰模型式(35)，$ f_{n+1}^{{\text{CF1}}} $、$ f_{n+1}^{{\text{CF2}}} $、$ f_{n+1}^{{\text{CF3}}} $的表达式分别为

(40)$ f_{n+1}^{{\text{CF1}}} = {\left. {\frac{{\partial f}}{{\partial s}}} \right|_{({s^ * },{v^ * },0)}} = \frac{{2{a_{\max }}}}{{{s^0}}}{\left(\frac{{{s^ * }}}{{{s^0}}}\right)^2} \geqslant 0, $

(41)$ f_{n+1}^{{\text{CF2}}} = {\left. {\frac{{\partial f}}{{\partial v}}} \right|_{({s^ * },{v^ * },0)}} = - \frac{{4{a_{\max }}}}{{{v^0}}}{\left(\frac{{4{a_{\max }}}}{{{v^0}}}\right)^3} - \frac{{2{a_{\max }}T}}{{{s^0}}}\left(\frac{{{s^ * }}}{{{s^0}}}\right) \leqslant 0, $

(42)$ f_{n+1}^{{\text{CF3}}} = {\left. {\frac{{\partial f}}{{\partial \Delta v}}} \right|_{({s^ * },{v^ * },0)}} = - \frac{{{a_{\max }}{v_i}(t)}}{{{s^0}\sqrt {{a_{\max }}{b_{c{\text{omf}}}}} }}\left(\frac{{{s^ * }}}{{{s^0}}}\right) \leqslant 0. $

由式(39)得到差分方程：

(43)$ \begin{split}& {{\dot \xi }_{n+1}}(t+2T') - {{\dot \xi }_{n+1}}(t+T') = T'f_{n+1}^{{\text{CF1}}}\left[ {{\xi _n}(t) - {\xi _{n+1}}(t)} \right]+ \\& {\text{ }}f_{n+1}^{{\text{CF2}}}\left[ {{\xi _{n+1}}(t+T') - {\xi _{n+1}}(t)} \right]+ \\& {\text{ }}f_{n+1}^{{\text{CF3}}}\left[ {{\xi _{n+1}}(t+T') - {\xi _{n+1}}(t) - {\xi _n}(t+T')+{\xi _n}(t)} \right]. \\[-1pt]\end{split} $

假设扰动$ \xi (t) = {\text{exp}}\;({\text{i}}n\omega +zt) $^[21]，则$ \dot \xi (t) = z \cdot \exp\;({\text{i}}n\omega + zt) $、$ \ddot \xi (t) = {z^2} \cdot \exp \;({\text{i}}n\omega + zt) $. 将$ \xi (t) $、$ \dot \xi (t) $、$ \ddot \xi (t) $代入式(43)，化简整理得到

(44)$ \begin{split} & ({{\mathrm{e}}^{T'z}} - 1)\left[ {{{\mathrm{e}}^{T'z}}z - f_{n+1}^{{\text{CF}}2}+f_{n+1}^{{\text{CF3}}}({{\mathrm{e}}^{ - {\text{i}}\omega }} - 1)} \right] =\\&\qquad T'f_{n+1}^{{\text{CF1}}}({{\mathrm{e}}^{ - {\text{i}}\omega }} - 1).\end{split} $

将$ z $展开为$ {z_1}({\text{i}}\omega )+{z_2}{({\text{i}}\omega )^2} $代入式(44)，得到系数$ z $的一阶$ {z_1} $和二阶项$ {z_2} $，

(45)$ {z_1} = {{f_{n+1}^{{\text{CF1}}}} / {f_{n+1}^{{\text{CF2}}}}}, $

(46)$ {z_2} = \frac{{(1 - f_{n+1}^{{\text{CF2}}}T')z_1^2 - f_{n+1}^{{\text{CF3}}}{z_1} - \dfrac{1}{2}f_{n+1}^{{\text{CF1}}}}}{{f_{n+1}^{{\text{CF1}}}}}. $

如果$ {z_2} < 0 $，则均匀稳态车流在受到扰动后，偏离将无法收敛到0；如果$ {z_2} > 0 $，则$ \xi (t) $产生的振荡幅值将呈衰减趋势，均匀稳态车流可以保持稳定^[21]. 将式(45)代入式(46)，根据稳定判别条件$ {z_2} > 0 $，化简整理得到HDV跟驰稳定性条件：

(47)$f_{n+1}^{{\text{CF2}}}f_{n+1}^{{\text{CF3}}}+\frac{1}{2}f_{n+1}^{{\text{CF2}}}f_{n+1}^{{\text{CF2}}}+f_{n+1}^{{\text{CF1}}}f_{n+1}^{{\text{CF2}}}{T^{'}} - f_{n+1}^{{\text{CF1}}} > 0.$

依次类推，跟随序列中HDV的IDM参数须遵循式(47)对应的稳定性条件，以保持跟驰稳定.

验证混合车辆队列协作控制算法的协作性能和能耗优化效果. 所有的仿真算例均在Matlab R2018a环境下运行，计算平台参数为Intel Core i7 CPU/3.40 GHZ. 仿真使用的CAV通信拓扑如图2所示. RTODMPC-IDM参数设置如表1所示，其中IDM参数根据跟驰模型稳态方程式(34)和跟驰稳定性条件式(47)选取，$ {U}[ - 200,{\text{ }}200] $为$ - 200 $~$ 200 $服从均匀分布的随机数. 仿真对车辆参数进行差异化设置，体现混合队列中车辆成员的异质性^[20,22].

设置2个典型仿真场景：恒速领航车辆场景和变速领航车辆场景，验证混合车辆队列协作中的跟踪性能以及加速度和控制输入是否满足乘客舒适度约束. 恒速领航车辆场景中领航车辆以恒定的速度行驶，是比较温和的条件；变速领航车辆场景中领航车辆速度存在调整.

设跟随车辆序列n = 3，q = 3. 仿真开始时，领航车辆的初始状态为$ {a_0}(0) = 0{\text{ }}{{\text{m}} / {{{\text{s}}^2}}} $、$ {v_0}(0) = 10{\text{ }}{{\text{m}} / {\text{s}}} $、$ {p_0}(0) = 0{\text{ m}} $；仿真中CAV间选取PF通信拓扑. 跟随车辆$ i $的初始状态为 $ {p_i}(0) = - i \cdot {d_0}+U{\rm{[ - 4, \;4]\; m}} $、$ {v_i}(0) = 10+U{\rm{[ - 2,\; 2] }}{{\text{ m}} / {\text{s}}} $、$ {a_i}(0) = $0 m/s². 采用RTODMPC-IDM的混合队列各车辆位置、速度、加速度随时间的变化情况如图3所示；车辆位置、速度、加速度的误差随时间的变化情况如图4所示，CAV的控制输入随时间的变化情况如图5所示. 结合图3、4可知，跟踪误差渐近收敛于零，这验证了所提控制方法的跟踪性能；由图3、5可知，车辆加速度始终位于$ \left[ { - 2,\;2} \right] $，控制输入始终位于$ \left[ { - 1,\;1} \right] $，因此，加速度及控制输入均满足乘客舒适度约束.

变速领航车辆场景的车辆个数和跟随车辆初始状态与恒速领航车辆场景设置相同，CAV间选取PLF通信拓扑. 本场景中领航车辆速度存在变化，加速度设置为

(48)$ {a_0}(t) = \left\{ \begin{array}{ll} 0,&{\text{ 0 s}} \leqslant t < 3{\text{ s}}; \\ t - 3,&{\text{ 3 s}} \leqslant t < 4{\text{ s}} ; \\ {\text{1}},&{\text{ 4 s}} \leqslant t < 5{\text{ s}}; \\ 6 - t,&{\text{ 5 s}} \leqslant t < 7{\text{ s}} ; \\ - 1,&{\text{ }}7{\text{ s}} \leqslant t < 8{\text{ s}} ; \\ t - 9,&{\text{ }}8{\text{ s}} \leqslant t < 9{\text{ s}}; \\ {\text{0, }}&{\text{ }}9{\text{ s}} \leqslant t {\text{.}} \\ \end{array} \right. $

采用RTODMPC-IDM的混合队列中各车辆位置、速度、加速度在控制下随时间的变化情况如图6所示；车辆位置、速度、加速度的误差随时间的变化情况如图7所示；CAV控制输入变化情况如图8所示. 结合图6、7可知，跟踪误差渐近收敛于零，这验证了所提控制方法在变速领航下仍能保证跟踪性能. 由图6、8可知，车辆加速度始终位于$ \left[ { - 2,\;2} \right] $，控制输入始终位于$ \left[ { - 1,\;1} \right] $，因此，在变速领航车辆场景下，加速度及控制输入依然均满足乘客舒适度约束.

以Human-Centered CACC^[14]为对比方法，仿真验证不同方法的能耗控制性能. 考虑到初始车辆状态分布会对编队协作过程中的能耗产生影响，选取3种不同初始场景，分别记为场景1、场景2、场景3. 设置混合车辆序列$ n=5、q=2 $；各场景领航车辆状态信息参考文献[14]. 3种场景的跟随车辆状态分布如表2所示. 为了衡量编队协作过程中混合队列成员单车能耗，定义从初始状态到稳态过程中的平均发动机功率为

(49)$ {P_{{\mathrm{a}},{\text{ }}i}} = \frac{1}{{{t_{\mathrm{s}}} - {t_0}}}{\displaystyle\sum\limits_{{t_0}}^{{t_{\mathrm{s}}}} {{P_i}(t)} }. $

式中：$ {t_0} $为起始时间，$ {t_{\mathrm{s}}} $为到达稳定状态的时间.如图9所示为2种方法在3种场景中的能耗对比. 可以看出，在随机的跟随车辆初始状态分布下，相比Human-Centered CACC，在RTODMPC-IDM控制下，3种场景的单车平均功率优化效果范围分别为16%~36%，16%~40%，6%~36%，单车平均功率最多降低了40%，分析原因：对于CAV，RTODMPC采用精确的VT-CPFM，能够降低车辆瞬时能耗；对于HDV，IDM跟驰模型具有良好的低能耗跟驰特性，因此车辆能耗得到优化.

本研究针对混合车辆队列最优能耗控制问题，以保障乘客舒适度为约束条件，提出结合实时优化分布式模型预测控制与智能驾驶员模型的协作控制方法. 对于CAV，设置分层控制架构，建立协作一致控制层（上层）和能耗优化控制层（下层），将精确的VT-CPFM引入能耗模型下层，从而在保证稳定协作的同时优化瞬时能耗. 对于HDV，采用IDM跟驰模型，分析得出跟驰模型的稳定性条件. 仿真结果表明：在所提方法控制下，车辆成员可以稳定实现协作，该方法能够满足乘客舒适度且跟踪性能良好. 3组不同车辆初始状态分布情况下的对比仿真结果表明，相较Human- Centered CACC，所提方法能够有效降低能耗. 未来研究将针对更普适的由CAV领航的CAV与HDV随机混合分布的队列，计划视本研究中的队列为子队列，以便将原混合队列从各领航CAV处拆分形成多个CAV+HDV子队列，在保障子队列间安全的前提下，通过控制每个子队列领航CAV完成整体队列控制.