Optimization Operation of Electric Locomotive Based on Two-stage AdaptiveGauss Re-Collocation Pseudospectral Approach
-
摘要: 针对中车株洲电力机车有限公司设计的HXD1-C64电力机车, 提出一种基于两阶段自适应Gauss配点重构的伪谱法, 用于列车优化操纵问题高效快速求解. 首先, 建立了HXD1-C64电力机车优化操纵数学模型; 然后, 在推导Legendre-Gauss配点公式的基础上, 给出控制变量两阶段自适应Gauss配点策略, 第1阶段采用对分配点, 第2阶段引入斜率变化分析对控制变量配点进行自适应细分和合并; 最后, 在运行时间最短目标下对HXD1-C64电力机车优化操纵进行仿真实验. 结果显示, 相较于控制变量参数化方法和传统高斯伪谱方法(Gauss pseudospectral method, GPM), 改进方法获得了更优性能指标和牵引力控制品质, 计算时间分别减少91.93 %和33.88 %, 表明了所提方法的有效性.
-
关键词:
- HXD1-C64电力机车 /
- 优化操纵 /
- 两阶段Gauss配点重构 /
- Gauss伪谱法
Abstract: A two-stage adaptive Gauss re-collocation pseudospectral approach is proposed for solving optimization operation problems of HXD1-C64 electric locomotive designed by CRRC (China Railway Rolling Stock Corporation) Zhuzhou Institute Co., LTD. Firstly, the optimization operation mathematical model of HXD1-C64 electric locomotive is established. Next, the Legendre-Gauss collocation formula is derived and then a two-stage adaptive Gauss collocation strategy is given, where a dichotomy collocation is adopted in the first stage and a slope change rate analysis is introduced in the second stage to adaptively subdivide or eliminate control variable collocation points. Finally, the simulation tests are carried out on HXD1-C64 electric locomotive to achieve time minimum operation. Simulation results show that the improved method obtains better performance index and traction control quality, compared with the control variable parameterization algorithm and the traditional Gauss pseudospectral method (GPM), while the computational time decreases by 91.93 % and 33.88 %, respectively, revealing the effectiveness of the proposed approach. -
电力机车对于我国铁路运输至关重要, 而机车优化操纵一直是铁路运营部门和国内外学者长期关注的热点问题[1], 其不仅可以在安全、正点前提下为机车提供操控策略以最大限度节约能源[2], 还可以建立司机模拟操纵指导系统, 寻找机车最优操作方法[3]. 迄今为止, 国内外铁路科研工作者针对机车优化操纵问题相继开展了一系列的研究[4-5].
电力机车优化操纵, 其核心可以归结为一个最优控制问题[6], 即通过建立性能指标函数, 选择不同的自变量, 在一定约束条件下寻找使性能指标最优的控制策略. 机车优化操纵研究最早可以追溯到20世纪60 年代, Ichikawa[7]基于Pontryagin极大值原理解决了简单的列车优化控制问题. 近年来, 随着计算机技术的不断发展, 采用数值方法实现列车优化操纵成为研究的热点和前沿[6, 8]. 如文献[2]基于控制参数化研究了列车节能操纵最优控制, 通过对列车牵引特性曲线光滑化近似, 有效获取了列车节能最优操纵序列及相应的参考速度曲线; 文献[9]针对具有时间需求地铁列车开展节能操纵控制研究, 设计了一种近似动态规划方法, 有效提升了优化问题求解效率; 文献[10]基于多目标粒子群优化算法研究了最少能耗下的机车自动驾驶速度曲线优化.
通常来讲, 最优控制问题的数值算法可以归纳为四类: 迭代动态规划算法、间接法、直接法和智能优化算法. 其中, 迭代动态规划算法源自于Luus[11]提出的动态优化理论; 间接法主要通过Pontryagin极大值原理将最优控制问题转化为两点边界问题求解[12]; 直接法则采用离散化手段将最优控制问题近似为非线性规划(Nonlinear programming, NLP)问题直接求解[13]; 智能优化算法的特点是通过智能化手段求解离散优化问题[14]. 与迭代动态规划算法和智能优化方法相比, 直接法具有计算量较小、计算效率高的优势, 而相较于间接法其收敛域更宽广, 且不需要高精度的初值[15], 因此, 直接法也成为了当前最优控制问题数值求解法中的主流方法.
近几年, 作为直接法中的一种, 高斯伪谱法(Gauss pseudospectral method, GPM)因具有求解精度高、计算量少、结构简单、求解效率高等优点在航空航天、化工领域最优控制问题求解上得到众多学者的青睐[16-17], 同时也在机车轨迹规划中得到了应用[18]. 该方法的思想是: 在正交配点处将连续最优控制问题离散化, 并通过全局插值多项式逼近状态变量和控制变量, 从而将最优控制问题转化为非线性规划NLP问题进行求解[19]. 随着电力机车智能化控制需求的不断提升, 列车优化操纵对最优控制问题求解效率提出了更高的要求. 然而, 传统GPM方法全局统一配点离散状态变量和控制变量, 为了得到精确结果, 通常需要大量增加配点数, 然而却降低了计算效率[20], 影响其在列车高效优化操纵计算上的应用; 另一方面, 统一配点容易引起控制曲线的震荡, 影响控制品质[21]. 因此, 如何有效解决GPM方法离散配点和计算效率之间的矛盾, 进一步提升GPM 方法在机车优化操纵中的应用已经成为国内外众多学者探索的课题[18, 22].
为了实现列车优化操纵问题高效数值优化求解, 本文以中车株洲电力机车有限公司设计的HXD1-C64电力机车为对象在Gauss伪谱法框架下提出两阶段自适应Gauss配点重构伪谱法. 首先, 结合HXD1-C64电力机车动力学方程, 建立Bolza 形式优化操纵数学模型; 其次, 结合Legendre多项式, 推导Legendre-Gauss配点求解公式; 进一步, 针对控制变量离散配点引入斜率变化提出两阶段自适应Gauss配点重构策略, 通过自适应配点细分、合并, 优化控制变量配点数进而提升求解效率; 最后, 在算法实现的基础上, 对HXD1-C64电力机车单阶段运行时间最短优化操纵问题进行仿真实验, 并与控制向量参数化(Control vector parameterization, CVP)方法和传统GPM方法对比以验证所提方法的有效性.
本文的组织结构安排如下: 第1节是对HXD1-C64电力机车优化操纵数学模型的描述; 第2节介绍了Gauss伪谱法基本框架; 第3节给出了两阶段自适应Gauss配点重构策略过程; 第4节分析了本文算法实现过程; 第5节进行仿真实验以测试所提算法的性能; 第6节对本文工作进行总结.
1. 机车优化操纵数学模型
由HXD1电力机车头和C64重车构成的电力机车(简称HXD1-C64电力机车)如图1所示, 其在实际运行过程中, 一般包括起动加速、保速、惰行和制动几个阶段[3], 各阶段速度曲线用图2描述.
图 2 机车牵引、惰行、制动组合运行示意图Fig. 2 Traction, coasting, brake composite sketch map of electric locomotive根据图2运行过程, 建立HXD1-C64电力机车优化操纵数学模型. 首先, 作如下假设:
1)牵引力为理想牵引力, 不考虑牵引力曲线;
2)机车在平直道或具有一定坡道的轨道上运行;
3)只考虑机车头、重车基本阻力以及坡道阻力.
结合牛顿力学定理, 分析得到HXD1-C64电力机车动力学方程为
$ \begin{split} &\dot s(t) = v(t), \;\;\dot v(t) = a(t) \\ &F(t) = Ma(t), \;\;M = {M_1} + {M_2} \\ &F(t) = {F_u}(t) - {F_o}(t) \\ &{F_o}(t) = {F_{{o_1}(t)}} + {F_{{o_2}(t)}} + {F_s}(t) \\ &{F_s}(t){\rm{ = }}Mg\sin \theta \end{split} $
(1) 其中,
$ s(t) $ 表示机车行驶路程,$ v(t) $ 表示机车行驶速度,$ a(t) $ 是机车的加速度,$ {M_1} $ 表示HDX1机车头的质量,$ {M_2} $ 表示C64重车的质量,$ \theta $ 表示轨道坡度,$ g $ 是重力加速度. 机车所受外力之和$ F(t) $ 为粘着牵引力$ {F_u}(t) $ 与阻力$ {F_o}(t) $ 之差, 阻力由机车头基本阻力$ {F_{{o_1}(t)}} $ 、重车基本阻力$ {F_{{o_2}(t)}} $ 和坡道阻力$ {F_s}(t) $ 构成, 其中机车头基本阻力$ {F_{{o_1}(t)}} $ 计算式为$ \begin{array}{l} {M_1} \times (2.23 + 0.0053v(t) + 0.000675{v^2}(t))g \end{array} $
(2) 重车基本阻力
$ {F_{{o_2}(t)}} $ 由下式计算$ \begin{array}{l} {M_2} \times (0.92 + 0.0048v(t) + 0.000125{v^2}(t))g \end{array} $
(3) 定义
$ {x}(t): = {[s(t)\;v(t)]^{\rm{T}}} $ ,$ u(t): = {F_u}(t) $ , 结合HXD1-C64电力机车动力学方程和运行约束条件建立Bolza形式电力机车优化操纵数学模型为$ \begin{split} & {\rm{Min}}\;J = \phi ({{x}_0}, {t_0}, {{x}_f}, {t_f}) + \int_{{t_0}}^{{t_f}} {g({x}(t), u(t), t){\rm d}t}\\ &{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\dot{{x}}(t) = f({x}(t), u(t), t)\\ &\;\;\;\;\;\;E({{x}}({t_0}), {t_0}, {{x}}({t_f}), {t_f}) = 0\\ &\;\;\;\;\;\;C({{x}}(t), u(t), t) \le0, \;\;\;t \in [{t_0}, {t_f}] \end{split} $
(4) 其中,
$ J $ 表示机车优化操纵问题的性能指标, 通常为运行时间最短、牵引能耗最少、牵引力与制动力变化频率最小以及速度曲线最优等, 可以由终值项和积分项构成;$ f({{x}}(t), u(t), t) $ 表示HXD1-C64电力机车动力学方程;$ E({{x}}({t_0}), {t_0}, {{x}}({t_f}), {t_f}) = 0 $ 表示机车运行过程中初值和终值约束条件;$ C({{x}}(t), $ $ u(t), t) \le 0 $ 表示机车运行过程中定时、变坡道、任意线路限速、牵引力和制动力限制等约束条件;$ t \in [{t_0}, {t_f}] $ 表示机车运行时间区间.2. Gauss伪谱法
采用Gauss伪谱法对HXD1-C64电力机车优化操纵问题进行求解, 其主要过程为: 将状态和控制变量在Legendre-Gauss 点处离散, 然后通过Lagrange插值多项式逼近状态变量和控制变量, 从而将系统状态方程和性能指标函数积分项转化为代数运算, 最后将优化操纵问题转化为以节点处状态、控制变量为待优化参量的NLP问题.
2.1 Lagrange插值多项式逼近
本文中, 状态变量和控制变量的离散逼近采用Lagrange插值方法完成. 已知
$ N+1 $ 组离散点及相应的函数值为$ \left\{ {({\tau _1}, {y_1}), ({\tau _2}, {y_2}), \cdots\ , ({\tau _{N + 1}}, {y_{N{\rm{ + }}1}})} \right\} $ , 则变量$ y(\tau ) $ 经$N $ 次多项式插值后可以得到:$ \begin{array}{l} y(\tau ) \approx Y(\tau ) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{{\hat L}_i}(\tau )} {y_i} \end{array} $
(5) 其中,
$ Y(\tau ) $ 是$ N $ 次多项式近似,$ {\hat L_i}(\tau ) $ 是$ N $ 次Lagrange插值多项式, 用如下公式求取:$ \begin{array}{l} {\hat L_i}(\tau ){\rm{ = }}\displaystyle\prod\limits_{j = 1, j \ne i}^{N + 1} {{{\tau - {\tau _j}} \over {{\tau _i} - {\tau_j}}}} \end{array} $
(6) 定义
$ {\hat L_i}({\tau _j}){\rm{ = }}{\delta _{ij}} $ , 分析可知Lagrange插值多项式具有以下特点:$ \begin{array}{l} { \hat L_i}({\tau _j}) = {\delta _{ij}} = \left\{ \begin{aligned} &1, \quad i = j \\ &0, \quad i \ne j \\ \end{aligned} \right. \end{array} $
(7) 结合式(5)和(7)可知, 插值多项式在离散点上满足
$ Y({\tau _i}) = {y_i} $ ,$ i = 1, 2, \cdots , N+1 $ , 因而需要选择合适的离散点来求取$ y_i $ 的值. 然而, 传统的方法一般将时间区间等分为$ N $ 段, 然后采用Lagrange 插值多项式的线性组合来近似, 但是这种等距离分隔容易造成龙格(Runge)现象. 因此, 本文选择逼近效果更好的Legendre多项式在$ [ - 1, 1] $ 区间进行离散配点的选择, 其$ N $ 次Legendre多项式为$ \begin{split} &{P_{n + 1}}\left( \tau \right) = (\tau - {\alpha _n}){P_n}\left( \tau \right) - \beta _n^2{P_{n - 1}}\left( \tau \right) \\ &{P_0}\left( \tau \right) = 1,\; {P_{ - 1}}\left( \tau \right) = 0, \;n = 0, 1,\cdots , N \end{split} $
(8) 进而, 通过求解式(8)的根则可以得到
$ N $ 个位于开区间$ ( - 1, 1) $ 内的离散点, 记为$ \tau_i $ ,$ i = {\rm{1}}, {\rm{2}}, \cdots , N $ , 称为Legendre Gauss (LG)配点. 分析可知,$ \tau_i $ 为式(8)的根, 其求解式可由定理1给出.定理 1. 已知
$ N $ 次Legendre多项式为$ \begin{array}{l} {P_{n + 1}}\left( \tau \right) = (\tau - {\alpha _n}){P_n}\left( \tau \right) - \beta _n^2{P_{n - 1}}\left( \tau \right) \\ {P_0}\left( \tau \right) = 1, \;{P_{ - 1}}\left( \tau \right) = 0, \;n = 0, 1, \cdots , N \end{array} $
该多项式的根可由矩阵
$ H $ 的特征值求解, 其中,$ \begin{array}{l} {{H}} = { \left[ \begin{array}{ccccc} {\alpha _0} & {\beta _1} & \quad & \quad & \quad\\ {\beta _1} & {\alpha _1} & {\beta _2} & \quad & \quad\\ \quad & \quad & \ddots & \quad & \quad\\ \quad & \quad & \quad & {\alpha _{n - 2}} & {\beta _{n - 1}}\\ \quad & \quad & \quad & {\beta _{n - 1}} & {\alpha _{n - 1}}\\ \end{array} \right ]} \end{array} $
证明. 根据定义,
$ N $ 次Legendre多项式的根为以下方程的解$ \begin{array}{l} {P_{n + 1}}\left( \tau \right) = (\tau - {\alpha _n}){P_n}\left( \tau \right) - \beta _n^2{P_{n - 1}}\left( \tau \right) = 0 \end{array} $
(9) 当
$ n = 0 $ 时, 有$ \begin{array}{l} {P_1}\left( \tau \right) = (\tau - {\alpha _0}){P_0}\left( \tau \right) - \beta _0^2{P_{ - 1}}\left( \tau \right) \end{array} $
(10) 又因为
$ {P_{ - 1}}\left( \tau \right) = 0 $ ,$ {P_{0}}\left( \tau \right) = 1 $ , 所以$ {P_1}\left( \tau \right) = (\tau - $ ${\alpha _0}) $ .当
$ n = 1 $ 时, 有$ \begin{array}{l} {P_2}\left( \tau \right) = (\tau - {\alpha _1}){P_1}\left( \tau \right) - \beta _1^2{P_0}\left( \tau \right) \end{array} $
(11) 根据
$ {P_0}\left( \tau \right) = 1 $ , 得$ {P_2}\left( \tau \right) = (\tau - {\alpha _1}){P_1}\left( \tau \right) - \beta _1^2 $ , 则$ {P_2}\left( \tau \right) = 0 $ 的根为$ (\tau - {\alpha _1})(\tau - {\alpha _0}) - \beta _1^2 = 0 $ 的解. 令$ {A_0} = $ $ {\alpha _0},$ $ {A_1} = {\left[ {array}{cc} {A_0} & {\beta _1} \\ {\beta _1} & {\alpha _1}\\ {array} \right ]} $ , 则$ {A_1} $ 的特征值为$ \begin{array}{l} \left| {\lambda I - {A_1}} \right| = (\lambda - {A_0})(\lambda - {\alpha _1}) - \beta _1^2 = 0\end{array} $
(12) 分析可知,
$ {P_2}\left( \tau \right) = 0 $ 的根为矩阵$ {A_1} $ 的特征值.同理,
$ {P_3}\left( \tau \right) = 0 $ 的根为$ \begin{array}{l} (\tau - {\alpha _2})\left| {\tau I - {A_1}} \right| - \beta _2^2\left| {\tau I - {A_0}} \right| = 0 \end{array} $
(13) 的解. 令
$ {A_2} = {\left[ {array}{cc} {A_1} & {\beta _2} \\ {\beta _2} & {\alpha _2}\\ {array}\right ]} $ , 则$ {A_2} $ 的特征值为$ \begin{array}{l} \left| {\lambda I - {A_2}} \right| = \left| \begin{array}{cc} {\lambda I - {A_1}} &- { {\beta _2}}\\ { - {\beta _2}} & {\lambda - {\alpha _2}} \end{array} \right|=\\ \qquad {\left| \begin{array}{ccc} \lambda - {\alpha _0}& - {\beta _1} \quad\\ - {\beta _1} & \lambda - {\alpha _1} & - {\beta _2}\\ & - {\beta _2} & \lambda - {\alpha _2}\\ \end{array}\right |}=\\ \qquad (\lambda - {\alpha _2})\left| {\lambda I - {A_1}} \right| + {\beta _2}( - {\beta _2})\left| {\lambda I - {A_0}} \right| = 0 \end{array} $
(14) 因此,
$ {P_3}\left( \tau \right) = 0 $ 的根为矩阵$ {A_2} $ 的特征值. 递推可得,$ N $ 次Legendre多项式的根为矩阵$ H $ 的特征值. □2.2 状态与控制变量离散化
结合第2.1节Lagrange插值方法, 采用Legendre多项式对状态变量和控制变量进行离散逼近. 首先, 引入时间变量
$ \tau $ , 结合下式对问题(4)中的时间区间进行尺度变换$ \begin{array}{l} t = \dfrac{{t_f} \;-\; {t_0}} {2}\tau + \dfrac{{t_f} \;+\; {t_0}}{2} \end{array} $
(15) 这样, 时间区间
$ [{t_0}, {t_f}] $ 就转化到$ [ - 1, 1] $ 区间上. 将式(15)代入问题(4), 得到转化后机车优化操纵问题为$ \begin{split} &{\rm{min}}\;\;J = \phi ({{{x}}_0}, {t_0}, {{{x}}_f}, {t_f})\;+\\ &\qquad{{{t_f} - {t_0}} \over 2}\int_{ - 1}^1 {g({{x}}(\tau ), u(\tau ), \tau ){\rm d}\tau } \\ &{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{{{\rm d}{{x}}(\tau )} \over {{\rm d}\tau }} = {{{t_f} - {t_0}} \over 2}f({{x}}(\tau ), u(\tau ), \tau ) \\ &\qquad E({{x}}( - 1), {t_0}, {{x}}(1), {t_f}) = 0 \\ &\qquad C({{x}}(\tau ), u(\tau ), \tau ) \le 0, \;\;\;\tau \in [ - 1, 1] \end{split} $
(16) 进一步, 将状态变量在LG配点上进行Lagrange插值近似
$ \begin{array}{l} {{x}}(\tau ) \approx {{X}}(\tau ) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{L_i}(\tau )} {{X}}({\tau _i}) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{L_i}(\tau )} {{{X}}_i} \\ {L_i}(\tau ){\rm{ = }}\displaystyle\prod\limits_{j = 1, j \ne i}^{N + 1} {{{\tau - {\tau _j}} \over {{\tau _i} - {\tau _j}}}} {\rm{ = }}{{b(\tau )} \over {(\tau - {\tau _i})\dot b({\tau _i})}} \end{array} $
(17) 其中,
$ b(\tau ){\rm{ = }}\displaystyle\prod\nolimits_{i{\rm{ = 1}}}^{N + 1} {(\tau - {\tau _i})} $ . 同时, 将状态变量导数也进行离散化, 得到$ \begin{array}{l} \dot {{x}}(\tau ) \approx \dot {{X}}(\tau ) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{{\dot L}_i}(\tau )} {{{X}}_i} \end{array} $
(18) 只考虑LG配点处状态变量的导数值, 令
$ \tau = {\tau _k} $ ,$ k = 1, 2, \cdots , N+1 $ , 代入式(18)可得$ \begin{array}{l} \dot {{x}}({\tau _k}) \approx \dot {{X}}({\tau _k}) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{{\dot L}_i}({\tau _k})} {{{X}}_i} = \displaystyle\sum\limits_{i = 1}^{N + 1} {{D_{k, i}}} {{{X}}_i} \end{array} $
(19) 其中,
$ \begin{array}{l} {D_{k, i}}\;{\rm{ = }}\;{\dot L_i}({\tau _k}) = \left\{ \begin{aligned} &{{\dot b({\tau _k})} \over {({\tau _k} - {\tau _i})\dot b({\tau _i})}}, & k \ne i \\ & {{\ddot b({\tau _k})} \over {2\dot b({\tau _k})}}, & k = i \\ \end{aligned} \right. \end{array} $
(20) 同理, 将控制变量在时间区间
$ \left[ { - 1, 1} \right] $ 上采用LG配点Lagrange插值近似, 即$ \begin{array}{l} u(\tau ) \approx U(\tau ) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{L_i}(\tau )} U({\tau _i}) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{L_i}(\tau )} {U_i} \end{array} $
(21) 进而, 状态方程可由如下约束条件代替
$ \begin{array}{l} \displaystyle\sum\limits_{i = 1}^{N + 1} {{D_{k, i}}} {{{X}}_i} = {{{t_f} - {t_0}} \over 2}f({{{X}}_k}, {U_k}, {\tau _k}) \end{array} $
(22) 同时, Lagrange项的积分项可以近似表示为
$ \begin{array}{l} \int_{ - 1}^1 {g( {{x}}(\tau ), u(\tau ), \tau ){\rm d}\tau } \approx \displaystyle\sum\limits_{k = 1}^{N + 1} {{w_k}g({{x}}({\tau _k}), u({\tau _k}), {\tau _k})} \end{array} $
(23) 其中,
$ {w_k} $ 是Gauss积分公式中的积分权重.最终, 问题(16)转化为如下形式NLP问题
$ \begin{split} &{\rm{min}}\;\;J = \phi ({{{X}}_0}, {t_0}, {{{X}}_f}, {t_f})\;+\\ &\;\;\;\;\;\;\;\ {{{t_f} - {t_0}} \over 2}\sum\limits_{k = 1}^{N + 1} {{w_k}g({{{X}}_k}, {U_k}, {\tau _k})} \\ &{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\sum\limits_{i = 1}^{N + 1} {{D_{k, i}}} {{{X}}_i} - {{{t_f} - {t_0}} \over 2}f({{{X}}_k}, {U_k}, {\tau _k}) = 0 \\ &\;\;\;\;\;\;\;E({{{X}}_1}, {t_0}, {{{X}}_f}, {t_f}) = 0 \\ &\;\;\;\;\;\;\;C({{{X}}_k}, {U_k}, {\tau _k}) \le 0 \\ &\;\;\;\;\;\;\;{{{X}}_f} - {{{X}}_1} - \sum\limits_{i = 1}^{N + 1} {{{{X}}_i}} \sum\limits_{i = 2}^{N + 1} {{w_i}{D_{k, i}}} = 0 \end{split} $
(24) 分析可知, 经过LG配点离散化后, 问题(24)是一个有限维数NLP问题. 因此, 基于梯度的非线性优化求解方法, 比如序列二次规划(Sequential quadratic programming, SQP)方法、内点法等可以用于该问题的求解.
3. 两阶段自适应Gauss配点重构
在传统GPM方法中, 状态变量和控制变量均采用相同LG配点进行多项式插值, 因此NLP问题(24)中需要优化的参数个数为
$ \begin{array}{l} NUM = ({r_u} + {r_x}) \times N \end{array} $
(25) 其中,
$ {r_x} $ ,$ {r_u} $ 分别表示状态变量、控制变量维数. 为了保证优化问题的求解精度, 通常需要增加LG 配点数. 然而大量增加配点会成倍增加优化参数$ NUM $ 的个数, 降低计算效率[20], 从而影响GPM方法在机车优化操纵中的应用. 另一方面, Xiao等[21]也指出统一配点容易引起控制曲线震荡, 特别是含有跳变节点的优化问题. 因此, 为了解决上述问题, 本文提出两阶段自适应Gauss配点重构策略. 第1阶段, 控制变量采用对分LG 配点降低优化参数个数并全局优化计算; 第2阶段引入斜率变化自适应重构控制变量配点. 该两阶段自适应Gauss 配点重构过程如图3所示.第1阶段: 状态变量采用式(17)进行LG配点离散化, 控制变量在时间区间
$ \left[ { - 1, 1} \right] $ 上采用对分LG 配点离散, 即$ \begin{array}{l} {U_{2l - 1}} = {U_{2l}}, \;\;l = 1, 2, \cdots , \dfrac{N + 1}{2} \end{array} $
(26) 其中,
$ N $ 为奇数. 从而控制变量配点数减少为$ (N + 1)/2 $ 个, 其逼近公式为$ \begin{array}{l} u(\tau ) \approx U(\tau ) = \displaystyle\sum\limits_{i = 1}^{N + 1} {{L_i}(\tau )} U({\tau _i}) = \displaystyle\sum\limits_{l = 1}^{\frac{N + 1}{2}} {\displaystyle\sum\limits_{i = 1}^{2l} {{L_i}(\tau )} {U_l}} \end{array} $
(27) 将对分配点得到的控制变量代入式(24)后分析可知, 此时优化参数个数降低为
$ ({r_u}/2 + {r_x}) \times N $ 个.第2阶段: 在第1阶段求解的基础上, 对控制变量配点重构. 首先, 定义第
$ r $ 个控制变量在$ {\tau _{2l + 1}} $ 配点处斜率变化为$ \begin{array}{l} {S_{r, l{\rm{ + }}1}} = \dfrac{U_{l + 1}^r - U_l^r} {{\tau _{2l + 1}} - {\tau _{2l - 1}}}, \;\;l = 1, 2, \cdots , \dfrac{N + 1}{2} - 1 \end{array} $
(28) 其中,
$ U_{l + 1}^r $ ,$ U_l^r $ 分别表示第$ r $ 个控制变量在$ {\tau _{2l + 1}} $ 和$ {\tau _{2l - 1}} $ 配点处的值. 设斜率变化上下限阈值参数为$ {\varepsilon _{\rm{U}}} $ 和$ {\varepsilon _{\rm{L}}} $ , 分别进行对分细化和配点合并判断. 如果$ \begin{array}{l} {S_{r, l{\rm{ + }}1}} \ge {\varepsilon _{\rm{U}}} \end{array} $
(29) 则控制变量
$ U_{l + 1}^r $ 和$ U_l^r $ 形成突变, 需要对分细化, 在区间$ [{\tau _{2l - 1}}, {\tau _{2l + 1}}] $ 增加配点个数. 如果$ \begin{array}{l} {S_{r, l{\rm{ + }}1}}\le {\varepsilon _{\rm{L}}} \end{array} $
(30) 则控制变量
$ U_{l + 1}^r $ 和$ U_l^r $ 在区间$ [{\tau _{2l - 1}}, {\tau _{2l + 1}}] $ 光滑, 设置$ \begin{array}{l} U_{l + 1}^r = U_l^r \end{array} $
(31) 进而, 适当减少光滑区域配点数以提高计算效率. 最终, 通过自适应迭代得到配点对分细化和合并后的控制变量配点个数. 进一步, 两阶段自适应Gauss配点重构后可以分别得到状态变量和控制变量在时间区间
$ \left[ { - 1, 1} \right] $ 上的LG配点, 从而实现配点的自适应调整.4. 算法实现
基于上述算法推导, 给出HXD1-C64电力机车优化操纵问题两阶段自适应Gauss配点伪谱法实现过程, 其算法流程如图4所示, 主要包括以下几个步骤:
图 4 两阶段自适应Gauss配点伪谱法算法流程图Fig. 4 Flow chart of two-stage adaptive Gauss re-collocation pseudospectral approach步骤 1. HXD1-C64电力机车优化操纵问题数学模型导入;
步骤 2. 设置算法优化参数;
步骤 3. 第1阶段状态变量LG配点、控制变量对分LG配点, 得到有限维NLP问题;
步骤 4. 求解NLP问题, 如果求解成功则进入步骤5, 否则判断是否需要LG配点重选, 如果需要则进入步骤3, 否则结束;
步骤 5. 第2阶段配点重构, 根据第1阶段求解结果计算斜率变化, 重构控制变量LG配点并得到新的有限维NLP问题, 进入步骤6;
步骤 6. 求解有限维NLP问题, 如果求解成功则结束, 否则判断是否需要LG配点重选, 如果需要选则进入步骤3, 否则结束.
在算法实施过程中, Legendre多项式参数
$ {\alpha _n} $ 和$ {\beta _n} $ 的取值为$ \begin{array}{l} {\alpha _n} = 0, \;\;{\beta _n} =\frac{ n}{\sqrt {4{n^2} - 1}} \end{array} $
(32) Gauss积分权重因子
$ {w_k} $ 由拉格朗日多项式插值公式的积分求取, 采用Abramowitz提出的方法[23]进行求解$ \begin{array}{l} {w_k}\;{\rm{ = }}\;\dfrac{2}{ {(1 - \tau _k^2){{({{\dot P}_N}({\tau _k}))}^2}}}, \;\;k = 1, 2, \cdots , N+1 \end{array} $
(33) 此外, 为了得到精确的结果, 本文设置斜率变化阈值参数
$ {\varepsilon _{\rm{U}}} $ 和$ {\varepsilon _{\rm{L}}} $ 相等.5. 优化仿真实验
以HXD1-C64电力机车为对象, 考虑单阶段(即运行固定路程)时间最短优化操纵进行仿真实验. 分析可知, 时间最短优化操纵即运行时间
$ {t_f} $ 最小, 因此可只采用终值项表示, 因此问题(4)中优化性能指标函数为$ \begin{array}{l} {\rm{min}}J = \phi ({{{x}}_0}, {t_0}, {{{x}}_f}, {t_f}) = {t_f} \end{array} $
(34) 同时, HXD1-C64电力机车相关参数取值如表1所示. 本文选取胡云卿等[24-25]实现的CVP方法、传统GPM法同本文两阶段自适应Gauss配点重构伪谱法进行对比测试, 其中CVP方法采用统一网格控制变量离散, 状态变量通过四级五阶龙格库塔方法进行求解; 传统GPM方法控制变量、状态变量采用全局LG配点. 仿真实验在Intel CORE i5/2.3 GHz CPU和DDR3/1600 MHz 8 GB内存组成的计算机上进行.
表 1 HXD1-C64电力机车参数Table 1 Parameters of HXD1-C64 electric locomotive参数项 参数值 HXD1机车头质量 200 $\times10^3\;{\rm{kg} }$ C64 重车质量 84 $\times10^3\;{\rm{kg} }$ 粘着牵引力约束 $- 461\;{\rm{kN} } \leq {F_u}(t) \le 532\;{\rm{kN} }$ 路程约束 $s(t) = 1\;500\;{\rm{m} }$ 速度约束 $0 \leq v(t) \leq 11.1\;{\rm{m/s} }$ 5.1 仿真实验 1
本实验在平直道(即坡道为零)且只考虑基本阻力轨道上进行. 初始条件为
$ {{x}}(0)= {\rm{ }}{\left[ {0, 0} \right]^{\rm{T}}} $ , 初始控制参数统一设置为0, 重力加速度取$ {\rm g} = 9.8\; {{\rm m}^2}/{\rm s} $ , 斜率变化阈值取$ {\varepsilon _{\rm{U}}} = {\varepsilon _{\rm{L}}} = 0.04 $ . 实验过程中, CVP方法、传统GPM法和本文两阶段自适应Gauss配点伪谱法均在两种不同控制变量配点数下进行, 其中CVP方法为20、40, 传统GPM法为34、48, 本文方法为17、24, 三种方法优化结果对比如表2所示.表 2 HXD1-C64电力机车仿真实验1结果对比Table 2 Numerical results comparison of Test 1 for HXD1-C64 electric locomotive方法 测试 u(t) 配点数 x(t) 配点数 离散方法 NLP 求解器 性能指标 (s) 求解时间 (s) CVP 1 20 − 统一网格 内点法 142.25 50.41 2 40 − 统一网格 内点法 141.81 210.63 传统 GPM 1 34 34 LG配点 内点法 141.81 8.12 2 48 48 LG配点 内点法 141.62 19.19 本文方法 阶段 1 17 34 两阶段自适应 LG 配点 内点法 141.74 5.77 阶段 2 15 34 阶段 1 24 48 两阶段自适应 LG 配点 内点法 141.59 12.09 阶段 2 17 48 结果显示, 三种方法均获得了相近优化性能指标. 但是相比于其他两种方法, 提出方法在17个控制变量配点数下, 性能指标更优, 为141.59 s. 此外, 相较于CVP方法, 两次测试下本文方法分别减少88.55 %和94.26 %计算时间. 同时, 与传统GPM法相比, 采用两阶段自适应Gauss配点后, 控制变量配点数明显减少, 求解时间也得到有效降低, 两次测试分别节省28.94 %和37 %计算耗时, 显示了本文方法在提升求解效率方面的有效性.
图5给出了三种方法各自在20、34和15个控制变量配点数下得到的控制曲线, 图中控制变量采用分段线性参数化方式逼近. 由图可知, 传统GPM方法控制曲线在零值附近出现了振荡现象, 而本文方法控制曲线更加平稳, 显示本文方法解决传统GPM方法振荡问题的有效性. 同时, 通过配点重构, 本文方法在跳变点附近控制也比另外两种方法更加精细, 实现了在更少配点数和求解时间下取得更优控制品质和性能指标. 图6给出了本文方法在图5控制策略下得到的系统状态曲线, 图中路程和速度状态变量均满足表1中的约束条件, 表明了所得控制策略的有效性. 因而, 本文算法在平直道上得到的控制策略可以为列车司机牵引力操纵提供有效指导, 进一步还可以集成到列车操纵系统中为自动驾驶提供算法支撑.
图 6 仿真实验1两阶段自适应Gauss配点法状态曲Fig. 6 State profiles of two-stage adaptive Gauss collocation method for Test 15.2 仿真实验 2
结合我国电气化铁路区间线路最大坡度限制要求, 本实验在30 ‰ 坡道, 考虑基本阻力与坡道阻力的轨道上进行, 参数设置与仿真实验1一致. 三种方法均在两种不同控制变量配点数下进行测试, 结果如表3所示. 结果显示, 本文方法在两次测试中取得了141.75 s和141.67 s的性能指标, 均优于CVP方法和传统GPM方法, 同时求解时间也大大减少, 与CVP方法相比平均减少92.98 %, 相比于与传统GPM方法平均降低35.71 %, 进一步表明了提出方法在HXD1-C64电力机车优化操纵问题中的优势和改进算法的有效性.
表 3 HXD1-C64电力机车仿真实验2结果对比Table 3 Numerical results comparison of Test 2 for HXD1-C64 electric locomotive方法 测试 u(t) 配点数 x(t) 配点数 离散方法 NLP 求解器 性能指标 (s) 求解时间 (s) CVP 1 20 − 统一网格 内点法 142.26 63.01 2 40 − 统一网格 内点法 141.83 238.79 传统GPM 1 34 34 LG 配点 内点法 141.97 8.10 2 48 48 LG 配点 内点法 141.76 21.47 本文方法 阶段 1 17 34 两阶段自适应 LG 配点 内点法 141.75 5.66 阶段 2 15 34 阶段 1 24 48 两阶段自适应 LG 配点 内点法 141.67 12.60 阶段 2 15 48 图7给出了三种方法各自在20、34和15个控制变量配点数下得到的控制曲线. 由图可知, 坡道的存在, 使得机车在保速阶段需要提供90 KN左右的牵引力来克服坡道阻力, 而本文方法在该阶段由于自适应Gauss配点重构作用, 细化了控制变量在跳变点处的配点, 有效消除了振荡, 改善了控制品质进而提升了性能指标. 图8给出了本文方法在图7控制曲线下得到的路程和速度状态曲线, 比较图6 (b)和图8 (b)可知, 相比于平直道, 坡道的存在使得加速阶段速度 有所降低, 基于此, 本文算法还可以为列车坡道行驶提供精准的机车牵引力控制方案, 提升列车优化控制水平.
图 8 仿真实验2两阶段自适应Gauss配点法状态曲线Fig. 8 State profiles of two-stage adaptive Gauss collocation method for Test 26. 结论
本文提出一种用于HXD1-C64电力机车优化操纵问题计算的两阶段自适应Gauss配点重构伪谱法, 在Gauss伪谱法框架下, 给出了Legendre Gauss配点求解公式, 提出了控制变量两阶段自适应Gauss配点重构策略, 可以有效减少优化配点个数, 提升求解效率. 通过对HXD1-C64电力机车在平直道和30 ‰ 坡道上进行优化操纵仿真测试发现, 与CVP方法与传统GPM法相比, 改进方法能够在更少优化配点数下求得更优性能指标, 同时求解时间平均降低91.93 %和33.88 %, 显示了本文方法在机车优化操纵问题计算上的有效性和实际应用价值. 需要说明的是, 本文研究中牵引力假设为理想牵引力, 下一步将引入牵引力曲线在多阶段多性能指标要求下开展HXD1-C64电力机车优化操纵算法研究.
-
图 2 机车牵引、惰行、制动组合运行示意图
Fig. 2 Traction, coasting, brake composite sketch map of electric locomotive
图 4 两阶段自适应Gauss配点伪谱法算法流程图
Fig. 4 Flow chart of two-stage adaptive Gauss re-collocation pseudospectral approach
图 6 仿真实验1两阶段自适应Gauss配点法状态曲
Fig. 6 State profiles of two-stage adaptive Gauss collocation method for Test 1
图 8 仿真实验2两阶段自适应Gauss配点法状态曲线
Fig. 8 State profiles of two-stage adaptive Gauss collocation method for Test 2
表 1 HXD1-C64电力机车参数
Table 1 Parameters of HXD1-C64 electric locomotive
参数项 参数值 HXD1机车头质量 200 $\times10^3\;{\rm{kg} }$ C64 重车质量 84 $\times10^3\;{\rm{kg} }$ 粘着牵引力约束 $- 461\;{\rm{kN} } \leq {F_u}(t) \le 532\;{\rm{kN} }$ 路程约束 $s(t) = 1\;500\;{\rm{m} }$ 速度约束 $0 \leq v(t) \leq 11.1\;{\rm{m/s} }$ 
下载: 导出CSV
表 2 HXD1-C64电力机车仿真实验1结果对比
Table 2 Numerical results comparison of Test 1 for HXD1-C64 electric locomotive
方法 测试 u(t) 配点数 x(t) 配点数 离散方法 NLP 求解器 性能指标 (s) 求解时间 (s) CVP 1 20 − 统一网格 内点法 142.25 50.41 2 40 − 统一网格 内点法 141.81 210.63 传统 GPM 1 34 34 LG配点 内点法 141.81 8.12 2 48 48 LG配点 内点法 141.62 19.19 本文方法 阶段 1 17 34 两阶段自适应 LG 配点 内点法 141.74 5.77 阶段 2 15 34 阶段 1 24 48 两阶段自适应 LG 配点 内点法 141.59 12.09 阶段 2 17 48
下载: 导出CSV
表 3 HXD1-C64电力机车仿真实验2结果对比
Table 3 Numerical results comparison of Test 2 for HXD1-C64 electric locomotive
方法 测试 u(t) 配点数 x(t) 配点数 离散方法 NLP 求解器 性能指标 (s) 求解时间 (s) CVP 1 20 − 统一网格 内点法 142.26 63.01 2 40 − 统一网格 内点法 141.83 238.79 传统GPM 1 34 34 LG 配点 内点法 141.97 8.10 2 48 48 LG 配点 内点法 141.76 21.47 本文方法 阶段 1 17 34 两阶段自适应 LG 配点 内点法 141.75 5.66 阶段 2 15 34 阶段 1 24 48 两阶段自适应 LG 配点 内点法 141.67 12.60 阶段 2 15 48
下载: 导出CSV
-
[1] 1 González-Gil A, Palacin R, Batty P, Powell J P. A systems approach to reduce urban rail energy consumption. Energy Conversion and Management, 2014, 80: 509−524 [2] 仲维锋, 徐洪泽. 基于控制参数化方法的列车节能操纵优化研究. 铁道学报, 2017, 39(6): 72−79 doi: 10.3969/j.issn.1001-8360.2017.06.0102 Zhong Wei-Feng, Xu Hong-Ze. Energy-efficient train operation optimization using control parameterization Method. Journal of the China Railway Society, 2017, 39(6): 72−79 doi: 10.3969/j.issn.1001-8360.2017.06.010 [3] 周继续, 贺德强, 姚晓阳, 向伟彬, 陈二恒. 电力机车优化操纵国内外研究算法综述. 铁道科学与工程学报, 2015, 12(1): 177−183 doi: 10.3969/j.issn.1672-7029.2015.01.0293 Zhou Ji-Xu, He De-Qiang, Yao Xiao-Yang, Xiang Wei-Bin, Chen Er-Heng. Domestic and international algorithm analysis of electric locomotive operation optimization. Journal of Railway Science and Engineering, 2015, 12(1): 177−183 doi: 10.3969/j.issn.1672-7029.2015.01.029 [4] 4 Zhao X H, Ke B R, Lian K L. Optimization of train speed curve for energy saving using efficient and accurate electric traction models on the mass rapid transit system. IEEE Transactions on Transportation Electrification, 2018, 4(4): 922−935 doi: 10.1109/TTE.2018.2851785 [5] 5 Azpilicueta L, Astrain J J, Lopez-Iturri P, Granda F, Vargas-Rosales C, Villadangos J, et al. Optimization and design of wireless systems for the implementation of context aware scenarios in railway passenger vehicles. IEEE Transactions on Intelligent Transportation Systems, 2017, 18(10): 2838−2850 doi: 10.1109/TITS.2017.2688858 [6] 荀径, 杨欣, 宁滨, 王义惠, 李坤妃. 列车节能操纵优化求解方法综述. 铁道学报, 2014, 36(4): 14−20 doi: 10.3969/j.issn.1001-8360.2014.04.0036 Xun Jing, Yang Xin, Ning Bin, Wang Yi-Hui, Li Kun-Fei. Survey on trajectory optimization for train operation. Journal of the China Railway Society, 2014, 36(4): 14−20 doi: 10.3969/j.issn.1001-8360.2014.04.003 [7] 7 Ichikawa K. Application of optimization theory for bounded state variable problems to the operation of train. Bulletin of JSME, 1968, 11(47): 857−865 doi: 10.1299/jsme1958.11.857 [8] 8 Ying W, Yuan G, Yu X, Hansen I A, Miao J. Optimization models for high-speed train unit routing problems. Computers and Industrial Engineering, 2019, 127: 1273−1281 [9] Liu R M, Li S K, Yang L X, Yin J T. Energy-efficient subway train scheduling design with time-dependent demand based on an approximate dynamic programming approach. IEEE Transactions on Systems, Man, and Cybernetics, Systems, 2019, DOI: 10.1109/TSMC.2018.2818263 [10] 10 Domínguez M, Fernández-Cardador A, Cucala A P, Gonsalves T, Fernández A. Multi objective particle swarm optimization algorithm for the design of efficient ATO speed profiles in metro lines. Engineering Applications of Artificial Intelligence, 2014, 29(3): 43−53 [11] 11 Luus R. Time-optimal control of state constrained linear systems. AIChE Journal, 1973, 19(6): 1262−1263 doi: 10.1002/aic.690190631 [12] 12 Liu P, Li X Y, Liu X G, Hu Y Q. An improved smoothing technique-based control vector parameterization method for optimal control problems with inequality path constraints. Optimal Control Applications and Methods, 2017, 38(4): 586−600 [13] 13 Lin Q, Loxton R, Teo K L. The control parameterization method for nonlinear optimal control: a survey. Journal of Industrial and Management Optimization, 2014, 10(1): 275−309 [14] 14 Zhang P, Chen H, Liu X, Zhang Z. An iterative multiobjective particle swarm optimization-based control vector parameterization for state constrained chemical and biochemical engineering problems. Biochemical Engineering Journal, 2015, 103: 138−151 doi: 10.1016/j.bej.2015.07.004 [15] 雍恩米, 陈磊, 唐国金. 飞行器轨迹优化数值方法综述. 宇航学报, 2008, 29(2): 397−406 doi: 10.3873/j.issn.1000-1328.2008.02.00215 Yong En-Mi, Chen Lei, Tang Guo-Jin. A survey of numerical methods for trajectory optimization of spacecraft. Journal of Astronautics, 2008, 29(2): 397−406 doi: 10.3873/j.issn.1000-1328.2008.02.002 [16] 16 Xiao L, Liu X, Ma L, Zhang Z. An effective pseudospectral method for constraint dynamic optimisation problems with characteristic times. International Journal of Control, 2017, 91(3): 1−11 [17] Han F, Wang Z L, He L, Wu H L, Yang G, Duan G R. Trajectory plan for an ultra-short distance on-orbit service based on the Gaussian pseudo-spectral method. IEEE/CAA Journal of Automatica Sinica, 2018, DOI: 10.1109/JAS.2017.7510892 [18] 18 Wang Y, Schutter B D, Boom T, Ning B. Optimal trajectory planning for trains-A pseudospectral method and a mixed integer linear programming approach. Transportation Research Part C Emerging Technologies, 2013, 29: 97−114 doi: 10.1016/j.trc.2013.01.007 [19] 19 Darby C L, Hager W W, Rao A V. An hp-adaptive pseudospectral method for solving optimal control problems. Optimal Control Applications and Methods, 2011, 32(4): 476−502 [20] 徐少兵, 李升波, 成波. 最优控制问题的 Legendre 伪谱法求解及其 应用. 控制与决策, 2014, 29(12): 2113−212020 Xu Shao-Bing, Li Sheng-Bo, Cheng Bo. Theory and application of Legendre pseudo-spectral method for solving optimal control problem. Control and Decision, 2014, 29(12): 2113−2120 [21] 21 Xiao L, Liu P, Liu X, Zhang Z, Wang Y, Yang C, et al. Sensitivity-based adaptive mesh refinement collocation method for dynamic optimization of chemical and biochemical processes. Bioprocess and Biosystems Engineering, 2017, 40(9): 1375−1389 [22] 22 Wang P, Goverde R M. Multi-train trajectory optimization for energy-efficient timetabling. European Journal of Operational Research, 2019, 272(2): 621−635 doi: 10.1016/j.ejor.2018.06.034 [23] Lide D R. Handbook of Mathematical Functions. Florida: CRC Press, 2018. [24] 24 Liu X, Hu Y, Feng J, Liu K. A novel penalty approach for nonlinear dynamic optimization problems with inequality path constraints. IEEE Transactions on Automatic Control, 2014, 22(2): 141−145 [25] 胡云卿, 刘兴高, 薛安克. 带不等式路径约束最优控制问题的惩罚 函数法. 自动化学报, 2013, 39(12): 1996−200125 Hu Yun-Qing, Liu Xing-Gao, Xue An-Ke. A penalty method for solving inequality path constrained optimal control problems. Acta Automatica Sinica, 2013, 39(12): 1996−2001 期刊类型引用(5)
1. 胡云卿，江帆，张征方. 机车自动驾驶关键技术与问题研究. 控制与信息技术. 2024(04): 1-10 . 
百度学术2. 刘俊彤，陈征，张泽. 考虑自动驾驶仪延迟的非线性最优末制导方法. 航空兵器. 2024(04): 64-71 . 百度学术3. 刘平，陈卓，刘明杰，朴昌浩，Jang Soohyun，万凯林. 阶段约束下Gauss配点离散化平行车位自动泊车轨迹规划. 汽车工程. 2023(07): 1163-1173 . 百度学术4. 刘平，刘航，仇国庆，刘兴高. 热率约束下高超声速飞行器Gauss时间网格参数化轨迹规划. 控制理论与应用. 2022(12): 2283-2292 . 百度学术5. 孙昊，孙青林，滕海山，周朋，陈增强. 复杂环境下考虑动力学约束的翼伞轨迹规划. 航空学报. 2021(03): 372-381 . 百度学术其他类型引用(9)
-
下载:
计量
- 文章访问数: 2113
- HTML全文浏览量: 355
- PDF下载量: 143
- 被引次数: 14
