当前社会经济有了巨大的发展，但人们出行时遇到的交通拥堵问题却日益严峻。政府、社会投入了大量资金，期望缓解交通拥堵，然而社会汽车保有量的增长速度远远快于道路建设速度，如何采用有效、科学的方案缓解交通拥堵，是当前人们面临的一个迫切问题^[1]。在过去几十年里，研究者们从不同的角度提出了各种交通模型，如格子流体力学模型、元胞自动机模型、跟驰模型和耦合映射模型等^[2-14]。在这些交通模型中，格子流体力学模型由于具有宏观离散性和计算的高效性，被人们广泛研究。1995年Bando等^[15]提出了最优速度模型(OV模型)，此后，许多学者对OV模型进行了改进研究，在原有模型中加入一些物理效应或交通因素等^{[3, 16]}。Nagatani^[17]在1998年首次提出了格子流体力学模型，该模型基于交通流流体动力学模型，从宏观上描述了连续车流之间的相互作用。运用非线性分析，得到在临界点附近的用mKdV方程描述的交通流拥堵和交通密度波的传播。1999年，Konishi等^[18]在离散时间的OV模型的基础上，首次提出了耦合映射(CM)模型中的延迟反馈控制(DFC)方法，通过理论分析和数值模拟，说明上述方法对控制交通拥堵有抑制作用。基于Nagatani和Konishi等的工作，许多研究者考虑不同的重要交通情况和控制参数，提出了各种控制模型^[19-25]，如考虑向后看效应^[26]、司机的延迟反应^{[27, 28]}、密度差效应^[29]、中断概率影响^[11]、上坡和下坡道路^[30]、基于非车道格子流体力学模型^[31]等。此外，唐铁桥等^[32]提出了改进的双车道格子流体力学模型。Ge等^[27]考虑驾驶员的反应延迟，得到了格子流体力学模型中的依赖于时间的Ginzburg-Landau (TDGL)方程，揭示了交通拥挤与交通流相变之间的关系。Jin等^[33]考虑延迟反馈控制，研究了OV模型中交通流的稳定性。此外，由于行人流和交通流具有许多共同的特征，Tian等^[34]在Nagatani的二维交通流模型的基础上提出了一种二维双向行人流。Cen等^[35]提出了考虑延迟时间连续累积效应的双向行人流的二维格子流体力学模型。2015年，Ge等^[36]在格子流体力学模型中引入了t时刻格子位置j与j+1的流量差作为反馈增益。Redhu等^[37]考虑向前看效应，提出了多相的格子流体力学模型。魏麒等^[38]研究了斜坡道路考虑平均流量差预期效应的格子流体力学模型。本研究根据Nagatani的格子流体力学模型，考虑把下游平均优化车流量与当前延迟车流量之间的差作为抑制交通阻塞的反馈控制项，建立控制方程。利用线性稳定性分析和非线性分析，推导出稳定条件以及描述密度波的mKdV方程的解，通过数值模拟验证理论分析结果。通过加入反馈控制项，研究反馈增益大小与反应时间长短对交通流拥堵的抑制作用。

1 模型建立

1998年Nagatani^[17]首次提出了格子流体力学模型，方程描述如下。

连续性方程：

$ \partial_t \rho_j+\rho_0\left(q_j-q_{j-1}\right)=0, $

(1)

动力学方程：

$ \partial_t q_j=a \rho_0 V\left(\rho_{j+1}(t)\right)-a q_j, $

(2)

其中，下标j表示格子位置；q_j=ρ_jv_j代表j格子位置时的交通流量；ρ_j和v_j分别表示t时刻在j格子位置处的交通密度和瞬时速度；ρ₀为总的平均密度；参量a表示交通流敏感系数，是延迟时间τ的倒数；取交通流优化速度函数V(ρ)如下：

$ V(\rho)=\frac{V_{\max }}{2}\left[\tanh \left(\frac{2}{\rho_0}-\frac{\rho}{\rho_0^2}-\frac{1}{\rho_c}\right)+\right.\left.\\ \tanh \left(\frac{1}{\rho_c}\right)\right]，$

(3)

其中，V_max和ρ_c分别表示最大速度和安全临界密度。该函数在ρ=ρ_c=ρ₀处有一个上界和一个拐点。

在方程(2)中，j+1位置的最优流量与j位置的瞬时流量之差决定了j位置处交通流量的变化，表明交通流的流动变化由下游交通流量与当前位置流量来决定。所以，为了揭示下游车流对当前位置交通流演化的影响，在格子流体力学方程的基础上，把在反应时间内下游平均优化车流量与当前延迟车流量的差作为延迟反馈控制项，得到以下控制方程：

$ \partial_t \rho_j+\rho_0\left(q_j-q_{j-1}\right)=0, $

(4)

$ \begin{array}{r} \partial_t q_j(t)=a \rho_0 V\left(\rho_{j+1}(t)\right)-a q_j(t)+ \\ a \lambda\left[\frac{1}{t_d} \int t_{t-t_d} \rho_0 V\left(\rho_{j+1}(t)\right) d t-q_j\left(t-t_d\right)\right], \end{array} $

(5)

其中，t_d为反应时间，参数λ为延迟反馈增益。为了抑制交通拥堵的产生，在式(5)中考虑把下游最优车流量的平均累积效应与当前延迟车流量的差作为延迟反馈控制。当λ=0时，该模型与Nagatani的模型一致。

利用梯度积分公式对式(5)中的积分进行简化，并结合式(4)和式(5)消去瞬时速度，得到密度方程：

$ \begin{array}{*{20}{l}} \;\;\;\;\;\;\;\;\partial_t^2 \rho_j(t)+a \partial_t \rho_j(t)+a \lambda \partial_t \rho_j\left(t-t_d\right)+ \\ a \rho_0^2 V\left[\rho_{j+1}(t)-\rho_j(t)\right]+\frac{1}{2} a \lambda \rho_0^2 V\left[\rho_{j+1}(t)-\rho_j(t)+\right. \\ \left.\rho_{j+1}\left(t-t_d\right)-\rho_j\left(t-t_d\right)\right]=0。\end{array} $

(6)

2 线性稳定性分析

应用线性稳定性分析方法来研究稳态交通流在小扰动作用下的演化情况。在稳态时，交通流处于均匀的密度和流量，在小扰动[ρ^*, q^*]作用下，在均匀稳态附近增加一个微扰，即ρ_j=ρ^*+ρ_j⁰; q_j=q^*+q_j⁰，因此，方程(4)和(5)分别可改写为

$ \partial_t \rho_{j+1}^0+\rho_0\left(q_{j+1}^0-q_j^0\right)=0 , $

(7)

$ \begin{gathered} \partial_t q_j^0(t)=a \rho_0 V^{\prime} \rho_{j+1}^0(t)-a q_j^0(t)+ \\ \frac{1}{2} a \lambda \rho_0 V^{\prime}\left[\rho_{j+1}^0(t)-\rho_{j+1}^0\left(t-t_d\right)\right]-a \lambda q_j^0\left(t-t_d\right), \end{gathered} $

(8)

其中，$V^{\prime}=\left.\left(\frac{\partial V(\rho)}{\partial \rho}\right)\right|_{\rho=\rho *}$是优化速度函数的一阶导数，且小于零。对方程(7)和(8)进行拉普拉斯变换，可得：

$ s P_{j+1}(s)-\rho_{j+1}(0)+\rho_0\left(Q_{j+1}(s)-Q_j(s)\right)=0, $

(9)

$ \begin{gathered} s Q_{j+1}(s)-q_j(s)=a \rho_0 V^{\prime} P_{j+1}(s)-a Q_j(s)+ \\ \frac{1}{2} a \lambda \rho_0 V^{\prime} P_{j+1}(s)\left(1+e^{-s t_d}\right)-a \lambda Q_j(s) e^{-s t_d}, \end{gathered} $

(10)

其中，L(ρ_j⁰)=P(s), L(q_j⁰)=Q(s), $L\left(\partial_t \rho_j^0\right)=$ $s P_j(s)-p_j(0), L\left(\partial_t q_j^0\right)=s Q_j(s)-q_j(0)$。L和s分别表示拉普拉斯变换和复数频率。消去方程(9)和(10)中的P_j+1(s)得到方程如下：

$ \begin{array}{*{20}{l}} \;\;\;\;\;Q_j(s)=\left(\frac{-a \rho_0^2 V^{\prime}-\lambda a \rho_0^2 V^{\prime}+\frac{\lambda a \rho_0^2 V^{\prime} s t _d}{2}}{s^2-a \rho_0^2 V^{\prime}+a s-\frac{a \lambda \rho_0^2 V^{\prime}}{2}-\frac{\lambda a \rho_0^2 V^{\prime} e^{{-s t}_d}}{2}+a \lambda s e^{{-s t}_d}}\right) \\ Q_{j+1}(s)+\frac{a \rho_0^2 V^{\prime}+\frac{a \lambda \rho_0 V^{\prime}\left(2-s t_d\right)}{2}}{d s} \rho_{j+1}(0)+\frac{s}{d s} q_j(0) 。\end{array} $

(11)

因此，传递函数为

$ G(s)=\frac{-a \rho_0^2 V^{\prime}-\lambda a \rho_0^2 V^{\prime}+\lambda a \rho_0^2 V^{\prime} s t_d / 2}{d s}, $

(12)

其中ds是特征多项式，

$ \begin{aligned} & d s=s^2-a \rho_0^2 V^{\prime}+a s-a \lambda \rho_0^2 V^{\prime} / 2- \\ & \lambda a \rho_0^2 V^{\prime} e^{-s t_d} / 2+a \lambda s e^{-s t_d} 。\end{aligned} $

(13)

Konishi等^{[18, 39]}基于控制理论得出关于交通拥堵的结论：如果满足以下两个条件，交通系统则不会发生拥堵, ①特征多项式ds是稳定的；②传递函数的H_∞范数小于或等于1，即$\|G(s)\|_{\infty}=\sup |G(j w)| \leqslant 1$。

根据Hurwitz稳定性准则，当ds大于零，特征多项式ds是稳定的^[37]。对于方程(13)，式中的每个系数都为正$\left[V^{\prime}=\left(\frac{\partial V(\rho)}{\partial \rho}\right) \leqslant 0\right]$，所以ds大于零，满足交通控制理论的第一个条件。对于第二个条件，分析$\|G(s)\|_{\infty}=\sup |G(j w)| \leqslant 1$。指数项e^－st_d≈1－st_d，化简后得到以下方程：

$ \begin{aligned} & |G(j w)|= \\ & \sqrt{\frac{\left(a \rho_0^2 V^{\prime}+a \lambda \rho_0^2 V^{\prime}\right)^2+\left(a \lambda \rho_0^2 V^{\prime} t_d w / 2\right)^2}{\left(w^2-a \lambda t_d w^2+a \rho_0^2 V^{\prime}+a \lambda \rho_0^2 V^{\prime}\right)^2+\left(a w+a \lambda w+a \lambda \rho_0^2 V^{\prime} t_d w / 2\right)^2}} \leqslant 1, \\ & w \in[0, \infty] 。\end{aligned} $

(14)

通过解方程(14)，得稳定条件：

$ a \geqslant-\frac{2 \rho_0^2 V^{\prime}}{1+\lambda-\lambda \rho_0^2 V^{\prime} t_d} 。$

(15)

显然，当λ=0，t_d=0，稳定条件与Nagatani模型的稳定条件(a≥－2ρ₀²V′)一致。

图 1反映了传递函数在不同反馈增益λ和反应时间t_d值的振幅变化。在图 1(a)中，反应时间t_d的值保持不变，改变反馈增益λ的值；而在图 1(b)中，反馈增益λ的值保持不变，改变反应时间t_d的值。从图 1可以看出，传递函数的幅值基本都小于1，说明增加反馈增益λ和反应时间t_d可以有效提高交通系统的稳定性。

3 非线性分析

为了研究在有限扰动作用下产生交通拥堵形成的密度波，采用非线性分析的约化摄动法研究在临界点附近交通密度波所满足的方程，并对方程进行求解。在临界点(ρ_c, a_c)附近，对空间慢变量X和时间慢变量T作变换如下：

$ X=\varepsilon(j+b t), T=\varepsilon^3 t, $

(16)

其中, j和t分别表示空间变量和时间变量；参数ε的取值范围为$0<\varepsilon \ll 1$；b是待定系数；在j位置处的密度为

$ \rho_j(t)=\rho_c+\varepsilon R(X, T) 。$

(17)

将方程(16)和(17)代入密度方程(6)中，利用泰勒级数展开，展开到ε的五阶，得到非线性偏微分方程如下：

$ \begin{aligned} & \quad \varepsilon^2\left\{a b(1+\lambda)+a \rho_c^2 V^{\prime}(1+\lambda)\right\} \partial_X R+\varepsilon^3\left\{b^2+\right. \\ & \left.\frac{1}{2} a \rho_c^2 V^{\prime}+\frac{1}{2} a \lambda \rho_c^2 V^{\prime}\left(1-b t_d\right)-a b^2 \lambda t_d\right\} \partial_X^2 R+ \\ & \varepsilon^4\left\{a(1+\lambda) \partial_T R+\frac{1}{6} a(1+\lambda) \rho_c^2 V^{\prime \prime \prime} \partial_X R^3+\left[\frac{1}{6} a(1+\right.\right. \\ & \left.\left.\lambda) \rho_c^2 V^{\prime}-\frac{1}{4} a \lambda \rho_c^2 V^{\prime} b t_d\right] \partial_X^3 R\right\}+\varepsilon^5\left\{\left[2 b-2 a \lambda b t_d-\right.\right. \\ & \left.\frac{1}{2} a \lambda \rho_c^2 V^{\prime} t_d\right] \partial_X \partial_T R+\left[\frac{1}{12} a(1+\lambda) \rho_c^2 V^{\prime \prime \prime}-\right. \\ & \left.\frac{1}{12} a \lambda \rho_c^2 V^{\prime \prime \prime} b t_d\right] \partial_X^2 R^3+\left[\frac{1}{24} a(1+\lambda) \rho_c^2 V^{\prime}-\right. \\ & \left.\left.\frac{1}{12} a \lambda \rho_c^2 V^{\prime} b t_d\right] \partial_X^4 R\right\}=0, \end{aligned} $

(18)

其中，$V^{\prime}=\left.\frac{\partial V(\rho)}{\partial \rho}\right|_{\rho=\rho_c}, V^{\prime \prime \prime}=\left.\frac{\partial^3 V(\rho)}{\partial \rho^3}\right|_{\rho=\rho_c}$。在临界点(ρ_c, a_c)附近，a= a_c(1-ε²)。令b=－ρ_c²V′，消去方程(18)式中ε的二阶项以及三阶项，化简得以下形式：

$ \begin{aligned} & \;\;\;\;\varepsilon^4\left(\partial_T R-g_1 \partial_X^3 R+g_2 \partial_X R^3\right)+\varepsilon^5\left(g_3 \partial_X^2 R+\right. \\ & \left.g_4 \partial_X^4 R+g_5 \partial_X^2 R^3\right)=0, \end{aligned} $

(19)

$ \begin{aligned} 其中, & g_1=-\frac{1}{1+\lambda}\left[\frac{1}{6}(1+\lambda) \rho_c^2 V^{\prime}-\frac{1}{4} \lambda \rho_c^2 V^{\prime} b t_d\right], \\ & g_2=\frac{1}{6} \rho_c^2 V^{\prime \prime \prime}, g_3=\frac{1}{a_c(1+\lambda)}\left(\rho_c^2 V^{\prime}\right)^2, \\ & g_4=\frac{1}{1+\lambda}\left[\frac{1}{24}(1+\lambda) \rho_c^2 V^{\prime}-\frac{1}{12} \lambda \rho_c^2 V^{\prime} b t_d\right], \\ & g_5=\frac{1}{1+\lambda}\left[\frac{1}{12}(1+\lambda) \rho_c^2 V-\frac{1}{12} \lambda \rho_c^2 V^{\prime \prime \prime} b t_d\right]。\end{aligned} $

为了得到标准的mKdV方程，采用以下变换：

$ T^{\prime}=g_1 T, R=\sqrt{\frac{g_1}{g_2}} R^{\prime} 。$

(20)

因此，方程(19)可写为

$ \partial_{T^{\prime}} R^{\prime}-\partial_X^3 R^{\prime}+\partial_X R^{\prime 3}+\varepsilon M\left[R^{\prime}\right]=0, $

(21)

其中，$M\left[R^{\prime}\right]=\sqrt{\frac{1}{g_1}}\left[g_3 \partial_X^2 R^{\prime}+\frac{g_1 g_5}{g_2} \partial_X^2 R^{\prime 3}+\right.$$\left.g_4 \partial_X^4 R^{\prime}\right] ; \varepsilon M\left[R^{\prime}\right]=O(\varepsilon)$为校正项。忽略O(ε)项，得到扭结-反扭结波的解为

$ R_0^{\prime}\left(X, T^{\prime}\right)=\sqrt{c} \tanh \sqrt{\frac{c}{2}}\left(X-c T^{\prime}\right), $

(22)

其中R′₀(X, T′)满足可解性条件：

$ \begin{aligned} & \quad\left(R_0^{\prime}, M\left[R^{\prime}\right]\right)=\int_{-\infty}^{+\infty} \mathrm{d} X R_0^{\prime}\left(X, T^{\prime}\right) M[(X, \\ & \left.\left.T^{\prime}\right)\right] \text { 。} \end{aligned} $

(23)

通过积分化简，可得传播速度c的值：

$ c=\frac{5 g_2 g_3}{2 g_2 g_4-3 g_1 g_5} 。$

(24)

因此，得到mKdV方程(20)的扭结-反扭结波的解为

$ \rho_j=\rho_c+\varepsilon \sqrt{\frac{g_1 c}{g_2}} \tanh \left(\sqrt{\frac{c}{2}}\left(X-\operatorname{cg} g_1 T\right)\right) 。$

(25)

扭结-反扭结波的振幅A为

$ A=\sqrt{\frac{g_1 c}{g_2}\left(\frac{a_c}{a}-1\right)} 。$

(26)

扭结-反扭结波的解代表低密度下的自由流和高密度下的堵塞相组成的共存相，其共存曲线可表示为

$ \rho_j=\rho_c \pm A \text {, } $

(27)

其中, ρ_j=ρ_c－A表示交通阶段中的自由流；ρ_j=ρ_c+A表示交通阶段中的堵塞相。

图 2为交通控制方程中密度敏感度(ρ, a)相图。每个图中都有3对曲线，实线为中性稳定线，虚线为共存曲线。根据交通流稳定性条件，将相图划分为3个区域：共存曲线以上的稳定区，中性稳定线以下的不稳定区，以及共存曲线与中性稳定线之间的亚稳定区。由图 2可以看出，中性稳定线和共存曲线随着反馈增益λ和反应时间t_d值的增加而下降，表明增加反馈增益λ和反应时间t_d可以有效地抑制交通拥堵产生。

4 数值模拟

采用周期边界条件进行数值模拟，验证线性稳定性分析和非线性分析的结论。为了便于模拟，对方程(6)应用差分方法进行离散，得到差分方程如下：

$ \begin{array}{*{20}{l}} \;\;\;\;\;\rho_j(t+2 \Delta t)=2 \rho_j(t+\Delta t)-\rho_j(t)- \\ a \Delta t\left[\rho_j(t+\Delta t)-\rho_j(t)\right]-a \rho_0^2(\Delta t)^2\left[V\left(\rho_{j+1}(t)\right)-\right. \\ \left.V\left(\rho_j(t)\right)\right]-a \lambda \Delta t\left[\rho_j\left(t-t_d+\Delta t\right)-\rho_j\left(t-t_d\right)\right]- \\ \frac{1}{2} a \lambda \rho_0^2(\Delta t)^2\left[V\left(\rho_{j+1}(t)\right)-V\left(\rho_j(t)\right)+V\left(\rho_{j+1}(t-\right.\right. \\ \left.\left.\left.t_d\right)\right)-V\left(\rho_j\left(t-t_d\right)\right)\right] 。\end{array} $

(28)

数值模拟考虑道路离散格子单元数N=100，交通流初始状态均匀，密度为ρ₀=ρ_c=0.25，敏感度系数a=1.65，最大速度V_max=2，模拟的时间步长Δt=0.1，交通流的初始扰动设置为

$ \rho_j(w)=\left\{\begin{array}{l} \rho_0 ; j \neq 50, 51 \\ \rho_0+0.1 ; j=50, \\ \rho_0-0.1 ; j=51 \end{array}\right. $

其中, w=1, 2, …, 5表示初始的5个时间步长。

数值模拟结果如图 3、图 4和图 5所示。图 3显示了密度随时间变化的情况，从图中可以清楚地观察到扭结-反扭结密度波。在图 3(a)中，当反馈增益λ=0.0和反应时间t_d=0.0时，密度波的振幅急剧振荡。对比图 3(b)、(c)和对比图 3(c)、(d)、(e)可知，密度波振幅波动随反馈增益λ和反应时间t_d的增加而减小，有利于交通流的稳定。

图 4为密度波的时空演化图。当不考虑反馈增益λ和反应时间t_d时，出现了扭结-反扭结的交通波现象，且波动幅度较大。当引入反馈增益λ和反应时间t_d时，密度波的波动减小。由于交通振荡的幅度与车辆系统的稳定性成正比，因此可知反馈增益λ和反应时间t_d的增加有利于车辆交通系统的稳定。

在不同反馈增益λ和反应时间t_d下，纵坐标为ρ₂₅(t)－ρ₂₅(t－1)与横坐标为ρ₂₅(t)相空间密度差的散点图如图 5所示。在不考虑反馈增益λ和反应时间t_d的情况下，图 5(a)表现出明显的滞后现象和混沌行为，可以看出交通系统趋于不稳定。对比图 5(b)、(c)、(d)和(e)可知，随着反馈增益λ和反应时间t_d的增加，滞后曲线逐渐缩小。在图 5(e)中，滞后曲线缩小成一个点，说明交通拥堵得到了有效缓解。

5 结论

本研究引入下游平均优化车流量与当前延迟车流量的差作为延迟反馈控制，提出了具有延迟反馈控制的交通流格子流体力学模型。通过线性稳定性分析，得到了传递函数和稳定条件；通过非线性分析得到了mKdV方程的扭结-反扭结密度波解。结果表明，反馈增益λ和反应时间t_d能有效缓解交通拥堵。为了验证理论分析，进行数值模拟，详细研究了在反馈控制作用下的密度-时间图和密度波的时空演化图。交通密度波的振荡减小表明反馈控制实现了抑制交通堵塞的控制策略。通过分析瞬时的交通密度差变化的散点图，发现在反馈增益λ=0时存在交通滞后和混沌现象。当加入反馈控制时，密度与密度差的空间散点图缩小为一个狭窄的有限区域，形成一个具有一定周期的极限环的稳态。交通滞后曲线的收缩表明交通系统的稳定性得到增强，有效地抑制了交通拥堵。仿真结果与理论分析结果一致，说明考虑下游平均优化车流量与当前延迟车流量的差作为反馈控制，有助于提高对交通系统的控制。