0 引言

很多专家对随机用户均衡(stochastic user equilibrium, SUE)问题的求解进行了研究, 主要包括基于路段和基于路径两类算法, 前者主要有连续平均法^[1-2]、Dial′s算法及变形算法^[3-4]等, 后者主要有Frank-Wolfe算法^[5]、梯度投影法^[6]、信赖域法^[7]、截断牛顿法^[8]、智能算法中的遗传算法^[9]、灵敏度分析法^[10-11]及其它算法^[12-20]。本研究主要探讨基于路径的算法, 总的来说比较有效的是梯度投影法和截断牛顿法。但是梯度投影法仅具有线性收敛速度, 为此Bojian Zhou在文献[8]中提出了改进的截断牛顿法, 并通过算例验证了方法的有效性。截断牛顿法实际上是牛顿法的一种改进, 因为牛顿法只能用来求解无约束优化问题, 但是截断牛顿法可以用来求解带线性等式约束的规划问题。使用牛顿法不可避免地要计算目标函数的Hessian矩阵, 如果变量个数较多, Hessian矩阵的计算和存储是一项计算量非常大的工作, 而且对于某些函数来说, 还有可能会出现Hessian矩阵不正定的情形。为此尝试用拟牛顿方程替换牛顿方程, 从而构造截断拟牛顿法。本研究中构建的改进截断拟牛顿法(modified truncated quasi-Newton, MTQN)选取目前最流行也是最有效的BFGS秩2拟牛顿校正公式, 既避免了Hessian矩阵的计算, 也保证了拟牛顿方程解的存在唯一性。线搜索采用的是Armijo非精确线搜索, 但这种搜索方法不能保证搜索方向的可行性, 为克服Armijo线搜索的不足, 选取改进的BFGS校正公式。最后将本研究所提算法与文献[8]中的截断牛顿法进行比较, 验证了结果的正确性和算法的适用性。

1 基于BFGS公式的改进截断拟牛顿法

考虑带线性等式约束的非线性规划问题$\text{P}_{1}$如下:

$ \begin{array}{l} \mathop {\min }\limits_{x \in {{\bf{R}}^n}} \;\;f\left( \mathit{\boldsymbol{x}} \right),\\ {\rm{s}}.\;{\rm{t}}.\;\mathit{\boldsymbol{Ax}} = \mathit{\boldsymbol{b}}, \end{array} $

式中: $f$: $\mathbf{R}^{n}$→$\mathbf{R}$二次连续可微; $\boldsymbol{A}$为$m$×$n$阶行满秩矩阵; $\boldsymbol{b}$∈$\mathbf{R}^{m}$。为保证全局解的存在, 此处假设$f$为严格凸函数。

假设$\boldsymbol{x}_{0}$为$\text{P}_{1}$的可行点, $\boldsymbol{Z}_{k}$为矩阵$\boldsymbol{A}$的零空间在第$k$步对应的$n$×($n$-$m$)维基矩阵, $\boldsymbol{Z}_{k}$依据基变量的选取来确定。则$\text{P}_{1}$的任何一个可行点$\boldsymbol{x}$均可表示为:

$ \mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{x}}_0} + {\mathit{\boldsymbol{Z}}_k}\mathit{\boldsymbol{y}}, $

式中: $\boldsymbol{y}$为$n$-$m$维列向量, 称为简化变量, 对应的$\boldsymbol{x}$称为原始变量。由此模型$\text{P}_{1}$可转化为无约束优化问题$\text{P}_{2}$如下:

$ \mathop {\min }\limits_{y \in {{\bf{R}}^{n - m}}} \;\;\;{\varphi _k}\left( \mathit{\boldsymbol{y}} \right) = f\left( {{\mathit{\boldsymbol{x}}_0} + {\mathit{\boldsymbol{Z}}_k}\mathit{\boldsymbol{y}}} \right), $

截取$φ_{k}$($\boldsymbol{y}$)的二阶Taylor展开式为

$ {\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k} + \mathit{\boldsymbol{l}}} \right) \approx {\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k}} \right) + {\mathit{\boldsymbol{l}}^{\rm{T}}}\nabla {\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k}} \right) + \frac{1}{2}{\mathit{\boldsymbol{l}}^{\rm{T}}}{\nabla ^2}{\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k}} \right)\mathit{\boldsymbol{l}}, $

对应的简化牛顿方程为

$ {{\mathit{\boldsymbol{\tilde G}}}_k}\mathit{\boldsymbol{l}} = - {{\mathit{\boldsymbol{\tilde g}}}_k}, $

式中:

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\tilde G}}}_k} = {\nabla ^2}{\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k}} \right) = \mathit{\boldsymbol{Z}}_k^{\rm{T}}{\nabla ^2}f\left( {{\mathit{\boldsymbol{x}}_k}} \right){\mathit{\boldsymbol{Z}}_k} = \mathit{\boldsymbol{Z}}_k^{\rm{T}}{\mathit{\boldsymbol{G}}_k}{\mathit{\boldsymbol{Z}}_k},}\\ {{{\mathit{\boldsymbol{\tilde g}}}_k} = \nabla {\varphi _k}\left( {{\mathit{\boldsymbol{y}}_k}} \right) = \mathit{\boldsymbol{Z}}_k^{\rm{T}}\nabla f\left( {{\mathit{\boldsymbol{x}}_k}} \right) = \mathit{\boldsymbol{Z}}_k^{\rm{T}}{\mathit{\boldsymbol{g}}_k},} \end{array} $

(1)

式中: $\boldsymbol{g}_{k}$和$\boldsymbol{G}_{k}$分别为$f$的梯度和Hessian矩阵; $\tilde{\boldsymbol{g}}_{k}$和$\tilde{\boldsymbol{G}}_{k}$分别称为$f$的简化梯度和简化Hessian矩阵。

为克服$\boldsymbol{G}_{k}$计算量大和可能不正定的缺点, 本研究中利用$\boldsymbol{H}_{k}$来近似$\boldsymbol{G}_{k}$, 并利用BFGS校正公式进行更新如下:

$ {\mathit{\boldsymbol{H}}_{k + 1}} = {\mathit{\boldsymbol{H}}_k} - \frac{{{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{s}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}}}{{\mathit{\boldsymbol{s}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{s}}_k}}} + \frac{{{\mathit{\boldsymbol{y}}_k}\mathit{\boldsymbol{y}}_k^{\rm{T}}}}{{\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k}}}, $

式中:$\boldsymbol{s}_{k}$=$\boldsymbol{x}_{k+1}$-$\boldsymbol{x}_{k}$为位移差; $\boldsymbol{y}_{k}$=$\boldsymbol{g}_{k+1}$-$\boldsymbol{g}_{k}$为梯度差。对应的近似简化Hessian矩阵:

$ {{\mathit{\boldsymbol{\tilde H}}}_k} = \mathit{\boldsymbol{Z}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{Z}}_k}, $

(2)

进而得到基于BFGS校正公式的简化拟牛顿方程为

$ {{\mathit{\boldsymbol{\tilde H}}}_k}\mathit{\boldsymbol{l}} = - {{\mathit{\boldsymbol{\tilde g}}}_k}。$

(3)

假设$\boldsymbol{l}_{k}$为式(3)的解, 则可以得到原始变量空间中的搜索方向

$ {\mathit{\boldsymbol{d}}_k} = {\mathit{\boldsymbol{Z}}_k}{\mathit{\boldsymbol{l}}_k}, $

(4)

值得注意的是, BFGS校正公式在使用精确线搜索或Wolfe搜索准则时可以保持$\boldsymbol{H}_{k}$的正定性, 进而保持$\tilde{\boldsymbol{H}}_{k}$的正定性, 从而使得式(3)存在唯一解$\boldsymbol{l}_{k}$, 此时有

$ \mathit{\boldsymbol{d}}_k^{\rm{T}}{\mathit{\boldsymbol{g}}_k} = \mathit{\boldsymbol{l}}_k^{\rm{T}}\mathit{\boldsymbol{Z}}_k^{\rm{T}}{\mathit{\boldsymbol{g}}_k} = - \mathit{\boldsymbol{\tilde g}}_k^{\rm{T}}\mathit{\boldsymbol{\tilde H}}_k^{ - 1}\mathit{\boldsymbol{Z}}_k^{\rm{T}}{\mathit{\boldsymbol{g}}_k} = - \mathit{\boldsymbol{\tilde g}}_k^{\rm{T}}\mathit{\boldsymbol{\tilde H}}_k^{ - 1}{{\mathit{\boldsymbol{\tilde g}}}_k} < 0, $

故对应的$\boldsymbol{d}_{k}$为$f$在$\boldsymbol{x}_{k}$处的下降方向, 但是在本研究的线搜索使用易于编程实现的Armijo准则, 为此参考文献[21], 采用校正方式

$ {\mathit{\boldsymbol{H}}_{k + 1}} = \left\{ \begin{array}{l} {\mathit{\boldsymbol{H}}_k},\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k} \le 0\\ {\mathit{\boldsymbol{H}}_k} - \frac{{{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{s}}_k}\mathit{\boldsymbol{s}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}}}{{\mathit{\boldsymbol{s}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}_k}{\mathit{\boldsymbol{s}}_k}}} + \frac{{{\mathit{\boldsymbol{y}}_k}\mathit{\boldsymbol{y}}_k^{\rm{T}}}}{{\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k}}},\mathit{\boldsymbol{y}}_k^{\rm{T}}{\mathit{\boldsymbol{s}}_k} > 0 \end{array} \right.。$

(5)

式(5)可以保证$\boldsymbol{H}_{k}$的正定型, 一般在初始步选取$\boldsymbol{H}_{0}$=$\boldsymbol{I}_{n}$或$\boldsymbol{H}_{0}$=$\boldsymbol{G}_{0}$=∇$^{2}f$($\boldsymbol{x}_{0}$)。

依据式(1)(2)(3)(5)及结论给出基于BFGS校正公式的改进截断拟牛顿算法(BFGS-MTQN)如下:

算法1 (BFGS-MTQN算法)

Step 0:给定初始点$\boldsymbol{x}_{0}$, 终止误差$ε$, 初始对称正定阵$\boldsymbol{H}_{0}$, 令$k$:=0;

Step 1:若$\boldsymbol{x}_{k}$满足终止误差$ε$, 停止迭代, 输出最优解$\boldsymbol{x}^{*}$≈$\boldsymbol{x}_{k}$, 否则执行Step 2;

Step 2:计算$\boldsymbol{g}_{k}$, 选取基变量并计算$\boldsymbol{A}$对应的基矩阵$\boldsymbol{Z}_{k}$, 利用公式(1)(2)计算得到$\tilde{\boldsymbol{g}}_{k}$和$\tilde{\boldsymbol{H}}_{k}$, 进而得到拟牛顿方程$\tilde{\boldsymbol{H}}_{k}\boldsymbol{l}$=-$\tilde{\boldsymbol{g}}_{k}$, 求得其解为$\boldsymbol{l}_{k}$;

Step 3:利用式(4)计算原变量空间中对应的搜索方向$\boldsymbol{d}_{k}$=$\boldsymbol{Z}_{k}\boldsymbol{l}_{k}$;

Step 4:沿着搜索方向$\boldsymbol{d}_{k}$利用Armijo准则(算法2)进行线搜索得到步长$λ_{k}$;

Step 5:置$\boldsymbol{x}_{k+1}$=$\boldsymbol{x}_{k}$+$λ_{k}\boldsymbol{d}_{k}$, 利用式(5)更新得到$\boldsymbol{H}_{k+1}$, 置$k$:=$k$+1, 转Step 1。

算法2 (Armijo非精确线搜索)

Step 1:给定参数$σ$∈(0, 0.5)及步长调整因子$ω$∈(0, 1), 置$λ$=1;

Step 2:验证以下不等式是否成立:

$ f\left( {{\mathit{\boldsymbol{x}}_k} + \lambda {\mathit{\boldsymbol{d}}_k}} \right) - f\left( {{\mathit{\boldsymbol{x}}_k}} \right) \le \sigma \lambda \mathit{\boldsymbol{g}}_k^{\rm{T}}{\mathit{\boldsymbol{d}}_k}; $

Step 3:判断, 若Step 2中不等式成立, 则终止迭代, 输出$λ_{k}$=$λ$, 否则令$λ$=$ωλ$, 转Step 2。

2 SUE模型及算法实现 2.1 SUE模型及等价模型

符号说明: $L$为所有路段组成的集合, ∀$a$∈$L$; $O$为所有出发点的集合, ∀ $r$∈$O$; $D$为所有目的点的集合, ∀$s$∈$D$; $I^{r,s}$为OD(origin-destination)对($r$, $s$)间的所有路径组成的集合; $x_{a}$为路段$a$上的交通流量, 组成的向量记为$\boldsymbol{x}$; $h^{r,s}_{i}$为OD对($r$, $s$)间第$i$条路径的流量, 组成的向量记为$\boldsymbol{h}^{r,s}$, 所有$\boldsymbol{h}^{r,s}$组成向量记为$\boldsymbol{h}$; $t_{a}$($x$)为路段$a$的行驶时间函数, 只与本路段的流量$x_{a}$有关, 并且关于$x_{a}$连续可微且严格单增; $C_{a}$为路段$a$的最大通行能力; $q^{r,s}$为OD对($r$, $s$)间的交通需求; $δ^{r,s}_{a,i}$为路段-路径关联变量, 若OD对($r$, $s$)间的第$i$条路径经过路段$a$, 则$δ^{r,s}_{a,i}$=1, 否则为0。

一般的SUE问题的模型$\text{P}_{3}$为:

$ \begin{array}{l} \min \;\;Z\left( {\mathit{\boldsymbol{h}},\mathit{\boldsymbol{x}}} \right) = \sum\limits_{a \in L} {\int_0^{{x_a}} {{t_a}\left( w \right){\rm{d}}w} } + \frac{1}{\theta }\sum\limits_{r \in O} {\sum\limits_{s \in D} {\sum\limits_{i \in {I^{r,s}}} {h_i^{r,s}\ln h_i^{r,s}} } } ,\\ {\rm{s}}.\;{\rm{t}}.\left\{ \begin{array}{l} \sum\limits_{i \in {I^{r,s}}} {h_i^{r,s} = {q^{r,s}}} ,\forall r \in O,s \in D\\ {x_a} = \sum\limits_{r \in O} {\sum\limits_{s \in D} {\sum\limits_{i \in {I^{r,s}}} {h_i^{r,s}\delta _{a,i}^{r,s}} } } ,\forall a \in L\\ h_i^{r,s} \ge 0,\forall r \in O,s \in D,i \in {I^{r,s}} \end{array} \right., \end{array} $

式中:$θ$＞0为常数; 第一个约束条件代表OD对-路径间流量的守恒约束; 第二个约束条件为路段-路径的流量守恒约束; 第三个约束条件为路径流量非负约束。

根据目标函数的形式易得第三个约束条件自动成立, 另外为了使用BFGS-MTQN算法进行求解, 将第二个约束条件代入目标函数, 最终模型$\text{P}_{3}$转化为以路径流量作为决策变量的线性等式约束规划$\text{P}_{4}$为:

$ \begin{array}{l} \min \;\;\;f\left( \mathit{\boldsymbol{h}} \right) = \sum\limits_{a \in L} {\int_0^{\sum\limits_{r \in R} {\sum\limits_{s \in S} {\sum\limits_{i \in {I^{r,s}}} {h_i^{r,s}\delta _{a,i}^{r,s}} } } } {{t_a}\left( w \right){\rm{d}}w} } + \frac{1}{\theta }\sum\limits_{r \in O} {\sum\limits_{s \in D} {\sum\limits_{i \in {I^{r,s}}} {h_i^{r,s}\ln h_i^{r,s}} } } ,\\ {\rm{s}}.\;{\rm{t}}.\;\;\sum\limits_{i \in {I^{r,s}}} {h_i^{r,s}} = {q^{r,s}},\forall r \in O,s \in D。\end{array} $

关于目标函数的严格凸性和等式约束系数矩阵的行满秩性在文献[10]中已给出证明, 此处不再详细介绍。

2.2 SUE模型BFGS-MTQN算法的收敛性

关于MTQN算法的收敛性在文献[21]中已经给出, 通过对路段行驶时间函数的假设可以验证随机用户均衡模型中的函数均具有文献[21]中所给的假设条件, 因此BFGS-MTQN算法在求解SUE模型时是全局收敛的, 而且具有超线性收敛速度。

2.3 SUE模型BFGS-MTQN算法实现中的问题

(1) 基变量及基矩阵的计算

针对每个OD对, 分别选取其中的任意一条路径的流量作为基变量, 得到相应的基矩阵$\boldsymbol{Z}_{k}$(具体的计算方法参考下文中算例), 利用式(5)得到的$\boldsymbol{H}_{k}$, 再利用式(2)计算得到$\tilde{\boldsymbol{H}}_{k}$, 并计算$\tilde{\boldsymbol{H}}_{k}$的条件数, 由此得到每个OD对间所有路径对应的条件数, 选取最小值对应的路径流量作为最终的基变量并计算相应的基矩阵$\boldsymbol{Z}_{k}$。

(2) Armijo线搜索的改进

SUE模型中要求路径流量非负, 即$\boldsymbol{h}^{r,s}_{k+1}$=$\boldsymbol{h}^{r,s}_{k}$+$λ_{k}\boldsymbol{d}^{r,s}_{k}$≥0, 为此构造$λ_{\text{max}}$为:

$ {\lambda _{\max }} = \left\{ {\begin{array}{*{20}{c}} {\infty ,}&{\mathit{\boldsymbol{d}}_k^{r,s} \ge 0}\\ {\mathop {\min }\limits_{i \in {I^{r,s}}} \left\{ { - \frac{{h_i^{r,s}}}{{d_i^{r,s}}}\left| {d_i^{r,s} < 0} \right.} \right\},}&{否则} \end{array}} \right.。$

然后在算法2的Step 2中加入不等式条件$λ$≤$λ_{\text{max}}$即可保证每次迭代所得解的可行性。

3 算例分析

在此考察如图 1所示的交通路网^[8], 图中共9个节点、12条路段, 各条路段上的行驶时间采用美国公共道路署提出的费用-流量(Bureau of Public Road, BPR)函数

$ {t_a}\left( {{x_a}} \right) = {t_{a0}}\left[ {1 + \alpha {{\left( {\frac{{{x_a}}}{{{C_a}}}} \right)}^\beta }} \right], $

式中: $t_{a0}$为路段$a$的自由流行驶时间, 在本例假设路段7、8、9的自由流行驶时间均为1, 其余路段均为2;$C_{a}$为每个路段的最大容量, 均设为100;取$α$=0.6, $β$=4。

在此算例中假设只有一个OD对(1, 9), 由此可以枚举此OD对之间的所有路径如表 1所示。

假设OD对需求为150, 并将其平均分配到各条路径上即可得到初始路径流量为$\boldsymbol{h}_{0}$=(25, 25, 25, 25, 25, 25), 由此可以得到此时的Hessian矩阵$\boldsymbol{G}_{0}$, 取$\boldsymbol{H}_{0}$=$\boldsymbol{G}_{0}$, 若选取第一条路径的流量当作基变量, 则线性等式约束对应的齐次线性方程组可变形为$h_{1}=-∑\limits^{6}_{i=2}h_{i}$, 从而可以得到其基础解系组成的基矩阵$\boldsymbol{Z}_{0}$及Hessian矩阵的近似简化矩阵$\tilde{\boldsymbol{H}}_{0}$为:

$ {\mathit{\boldsymbol{Z}}_0} = \left( {\begin{array}{*{20}{c}} { - 1}&{ - 1}&{ - 1}&{ - 1}&{ - 1}\\ 1&0&0&0&0\\ 0&1&0&0&0\\ 0&0&1&0&0\\ 0&0&0&1&0\\ 0&0&0&0&1 \end{array}} \right),{{\mathit{\boldsymbol{\tilde H}}}_0} = \mathit{\boldsymbol{Z}}_0^{\rm{T}}{\mathit{\boldsymbol{H}}_0}{\mathit{\boldsymbol{Z}}_0}。$

利用2-范数得到此时$\tilde{\boldsymbol{H}}_{0}$的条件数为8.2638(可简称为路径1的条件数), 类似的方法可以得到所有路径对应条件数为(8.2638, 8.2840, 8.2625, 7.9758, 8.2701, 7.8691), 根据文献[8]可知条件数越小收敛速度越快, 因此在本次迭代中选取路径6作为基变量后再进行循环。

为验证BFGS-MTQN算法的适用性, 首先将模型$\text{P}_{3}$及算法2中的参数设置如下:$θ$=0.5, $σ$=0.25, $ω$=0.5, 然后将两种不同OD需求下的算法结果与MTN方法所得结果进行比较如表 2所示。

表 2中的准确解是通过Matlab自带的非线性规划命令求解得到, 对应的目标函数值略小于MTN和BFGS-MTQN两种方法所得到的结果。由表 2可知:(1)本研究中构建的BFGS-MTQN算法得到的解的误差和迭代次数均小于MTN算法; (2)BFGS-MTQN算法得到的解虽然与准确解存在一定误差, 但迭代次数要远小于准确解算法。通过多次变化OD对需求和参数的值进行试验, 发现以上两个结论都是成立的, 从而验证了BFGS-MTQN算法在求解SUE问题时的合理性和优越性。

值得说明的是:文献[8]中关于线搜索算法中的参数$σ$的范围是(0, 1), 而在实际计算中发现当$σ$＞0.5时, 会使得线搜索的迭代次数突然增多, 而且在某些情况下还会陷入死循环, 为此将$σ$的范围修正为(0, 0.5)。在(0, 0.5)范围内, $σ$的取值并不会影响线搜索的迭代次数, 因此可以在(0, 0.5)随机选取$σ$的取值。

4 结语

本研究根据随机用户均衡模型的特点, 以路径流量作为决策变量, 构造了一种基于BFGS校正公式的截断拟牛顿算法, 在文中的假设条件下具有超线性的收敛速度。与截断牛顿法相比较, 基于BFGS校正公式的截断拟牛顿算法可以避免二阶Hessian矩阵的计算, 而且通过算例结果分析得到, 此算法在迭代次数和误差两个方面均优于截断牛顿法。最后给出了算法实施过程中一些实用的结论。

在以后的工作中, 从以下几个方面展开:(1)针对本研究中所得到的一些结论, 尝试能否给出理论性的证明; (2)正在尝试将多种拟牛顿公式所得结果进行比较及改进, 并且已经得到了一些比较理想的结果; (3)利用此算法来求解更大规模和更多类型的交通配流及改进模型。