考慮阻力約束的列車能量最優駕駛問題建模及分離迭代求解策略

劉良杰，馮江華，王 斌，胡云卿，黎向宇

(中車株洲電力機車研究所有限公司, 湖南 株洲 412001)

軌道交通車輛智能控制已成為當前的研究熱點，引起了國內外廣大學者的關注并取得了豐碩的成果[1-13]。文獻[1]從列車驅動策略節能和運行時刻表節能兩方面進行建模，采用Pontryagin極大值原理求解最優控制問題。文獻[2]提出了一種兩階段線性規劃算法優化地鐵列車運行時刻表，通過最小化列車運行過程中的能量消耗以及最大化利用列車制動時的再生能量，可實現節能19.27%以上。文獻[3-4]通過對列車運行時刻表進行優化求解實現節能。文獻[5-8]根據列車的特性、線路與速度約束等條件建立能量最優控制問題模型進行求解，實現節能。GE公司將列車的能量最優控制問題歸結為一個帶約束的非線性最優控制問題進行求解，并成功研發了優化節能裝置，并在北美內燃機車上應用[5]。西南交通大學馮曉云教授團隊借鑒優秀司機的駕駛經驗，建立能量最優問題模型，通過Pontryagin極大值原理進行求解[9-10]。北京交通大學徐洪澤教授團隊和毛保華教授團隊分別建立了列車的節能操作最優控制問題模型，徐洪澤教授團隊通過運行區間分段的方法將最優控制問題轉化為非線性規劃問題[11]；毛保華教授團隊通過遺傳算法求解節能操作最優控制問題[12]。文獻[13]建立了列車運行過程中能量最優的線性二次型最優模型，通過離散化策略將該問題轉化為凸二次規劃問題進行求解。

本文綜合考慮了列車的動力學模型，牽引、制動特性，列車阻力(坡道阻力和運行阻力)，線路限速等條件，在滿足安全、準點、平穩的條件下，建立了列車運行過程中的能量最優駕駛問題模型。列車運行過程的動力學方程(不考慮列車阻力)為常微分方程組(Ordinary Differential Equations，ODEs)，由于考慮了列車運行過程中的坡道阻力和運行阻力，列車運行過程的動力學方程增廣為微分代數方程組(Differential-Algebraic Equations，DAEs)，使得問題難以求解。為求解該模型，在時間域內將狀態變量和控制變量離散化[14-16]，將問題轉化為一般非線性規劃問題；然后提出了一種分離迭代策略，將該一般非線性規劃問題處理為一系列凸二次規劃問題，最后采用原-對偶預測校正內點算法求解[17]。

1 列車能量最優駕駛問題

列車運行過程的數學模型為

(1)

式中：t為時間；s(t)為列車的運行距離；v(t)為列車的運行速度；a(t)為列車的運行加速度；F(t)、r[v(t)]、gs[s(t)]分別為作用在列車單位質量上的牽引、制動力，運行阻力，坡道阻力。

r[v(t)]、gs[s(t)]分別為[18]

(2)

式中：α0,α1,α2為列車運行阻力表達式的系數；θ[s(t)]為列車運行過程中的坡道角。θ[s(t)]>0為上坡道，θ[s(t)]=0為平直道，θ[s(t)]<0為下坡道。g為重力加速度常數。

假設列車從A站運行到B站過程中，運行到地點C為t0(t0≥0)時刻，運行到地點D為tf(tf>t0)時刻，運行區間限速值為VLimit，最大加、減速度變化率為aRate(aRate>0)。依據列車的牽引、制動特性，設列車能夠達到的最大加速度值為amax，最大減速度值為-amin(amin<0)；設列車的最大牽引力為Fmax，最大電制動力為-Fmin(Fmin<0)，則列車在運行過程中受到的約束條件為

v(t)≤VLimit

(3)

(4)

amin≤a(t)≤amax

(5)

Fmin≤F≤Fmax

(6)

s(t0)=S0v(t0)=V0a(t0)=A0

(7)

s(tf)=S0+S=Sf,v(tf)=Vf,a(tf)=Af

(8)

式中：S0、V0、A0分別為t0時刻列車的運行位置、速度、加速度；S為從t0時刻起至tf時刻列車的運行距離；Sf、Vf、Af為tf時刻列車的運行位置、速度、加速度；t0,S0,V0,A0,tf,Sf,Vf,Af均為已知常數。

a(t)=F(t)-r[v(t)]-gs[s(t)]

r[v(t)]={α0+α1[3.6v(t)]+α2[3.6v(t)]2}×10-3

gs[s(t)]=gsin{θ[s(t)]}≈gθ[s(t)]

s(t0)=S0v(t0)=V0a(t0)=A0

s(tf)=Sfv(tf)=Vfa(tf)=Af

v(t)≤VLimit

amin≤a(t)≤amax

Fmin≤F(t)≤Fmax

(9)

2 離散化策略及問題求解

2.1 離散化策略

對時間區間t∈[t0,tf]進行離散化，等分為N個子時間區間，每個子時間區間的長度為T，則有

(10)

在離散的時間點序列t0,t1,…,tN-1,tN上，對加速度a(t)離散化得到離散的加速度序列a0,a1,…,aN-1,aN；對速度v(t)離散化得到離散的速度序列v0,v1,…,vN-1,vN；對運行距離s(t)離散化得到離散的運行距離序列s0,s1,…,sN-1,sN；對牽引力F(t)離散化得到離散的牽引力序列F0,F1,…,FN-1,FN；對列車的運行阻力r[v(t)]離散化得到離散的運行阻力r(v0),r(v1),…,r(vN-1),r(vN)；對列車的坡道阻力gs[s(t)]離散化得到離散的坡道阻力gs(s0),gs(s1),…,gs(sN-1),gs(sN)。那么，在離散點處的加速度值可表示為

ai=Fi-r(vi)-gs(si)i=0,1,…,N

(11)

在t0和tf時刻，由式(7)、式(8)可得

s0=S0v0=V0a0=A0

(12)

sN=SfvN=VfaN=Af

(13)

(14)

設在[ti-1,ti],i=1,…,N子區間列車的運行距離為Δsi，Δsi可計算為

?

i=3,4,…,N

(15)

由于si=si-1+Δsi=S0+Δs1+…+Δsi，可得

i=2,3,…,N

(16)

則式(9)中的等式約束可轉化為

a0=A0=F0-r(v0)-gs(s0)

(17)

aN=Af=FN-r(vN)-gs(sN)

(18)

(19)

sN=Sf≈S0+NV0T+

(20)

針對式(9)中的不等式約束，在任意時刻加速度ai需要滿足

amin≤Fi-r(vi)-gs(si)≤amaxi=1,…,N-1

(21)

牽引力Fi需要滿足：

Fmin≤Fi≤Fmaxi=0,…,N

(22)

aRatei=1,…,N

(23)

限速條件v(t)≤VLimit可轉化為

i=1,…,N-1

(24)

式(9)的性能指標函數通過離散化后可近似表示為

(25)

令F=[F0,F1,…,FN]T，式(9)可轉化為一般非線性規劃問題

s.t.F0-r(v0)-gs(s0)=A0

FN-r(vN)-gs(sN)=Af

Fi-r(vi)-gs(si)≥amini=1,…,N-1

-Fi+r(vi)+gs(si)≥-amaxi=1,…,N-1

Fi≥Fmini=0,1,…,N

-Fi≥-Fmaxi=0,1,…,N

Fi-r(vi)-gs(si)-

[Fi-1-r(vi-1)-gs(si-1)]≥-aRateT

i=1,…,N

Fi-1-r(vi-1)-gs(si-1)-

[Fi-r(vi)-gs(si)]≥-aRateT

i=1,…,N

i=1,…,N-1

(26)

在一般非線性規劃問題中(式(26))，有N+1個決策變量，4個等式約束，7N-1個不等式約束。

2.2 一般非線性規劃問題求解

當不考慮列車的阻力時，式(26)中不包含r(vi),gs(si),i=0,1,…,N，此時式(26)為凸二次規劃問題，直接采用原-對偶預測校正內點算法[17]進行求解。當考慮了列車運行過程中的坡道阻力和運行阻力時，式(26)為一般非線性規劃問題，增加了求解難度。在式(1)中，r[v(t)]為v(t)的函數，gs[s(t)]為s(t)的函數，不容易直接通過求解式(1)DAEs得到v(t),s(t)表達式，因此，對一般非線性規劃問題(式(26))直接進行求解非常困難。可以采用分離迭代策略：令v=[v0,v1,…,vN]T，s=[s0,s1,…,sN]T，假設初始估計值F(0)已知，由式(14)、式(16)可求得v(0),s(0)，代入式(26)中，則問題(26)變為凸二次規劃問題，通過求解可得到式(26)的一個近似解F(1)，再由式(14)和式(16)求得v(1),s(1)，代入式(26)繼續迭代求解凸二次規劃問題。假設迭代求解凸二次規劃問題k次后，得到序列{F(k)}，若存在一個充分小的正數ε>0，滿足

‖F(k)-F(k-1)‖≤ε

(27)

則終止迭代，由F(k)可求得v(k),s(k)。此時由式(11)、式(14)、式(16)可知F(k-1)≈F(k),v(k-1)≈v(k),s(k-1)≈s(k)，取F(*)=F(k)作為式(26)的一個近似最優解，v(*)=v(k)為近似最優速度曲線，s(*)=s(k)為近似最優運行距離曲線。分離迭代策略的流程圖見圖1。

圖1 分離迭代策略流程圖

在分離迭代策略中，將r(vi),gs(si),i=0,1,…,N代入式(26)，則問題可簡化為凸二次規劃問題，共有N+1個決策變量，4個等式約束，7N-1個不等式約束，其規范化表達式為

s.t.AEx=bE

AIx≥bI

(28)

式中：決策變量x=[F0,F1,…,FN]T∈RN+1；H=2I∈R(N+1)×(N+1)為正定矩陣，I∈R(N+1)×(N+1)為單位矩陣；向量p=0∈RN+1;常數項q=0；AE∈R4×(N+1)，AI∈R(7N-1)×(N+1)分別為等式約束和不等式約束的Jacobi矩陣；bE∈R4,bI∈R7N-1。

針對式(28)，采用原-對偶預測校正內點算法[17]進行求解。對式(28)不等式約束條件添加松弛變量y≥0。令μ∈Rm2≥0為不等式約束y≥0的拉格朗日乘子向量，λ∈Rm1為等式約束AEx=bE的拉格朗日乘子向量。設點z(k)=(x(k),y(k),λ(k),μ(k))為當前迭代點，則下一個迭代點為

(29)

Step1輸入相關參數H、p、q、AE、bE、AI、bI,選擇初始點(x(0),y(0),λ(0),μ(0))，并設定精度要求tol，令k=0。

Step6計算下一個迭代點為

Step7計算新的對偶間隔

3 算例計算與仿真

假設列車從A站運行到B站過程中，運行到地點C為t0=15 s時刻，在t0時刻列車的加速度A0=0.7 m/s2，速度V0=40 km/h，運行距離S0=100 m；運行到地點D為tf時刻，列車的加速度Af=-0.5 m/s2，速度Vf=5 km/h，運行距離Sf=S+S0=1 900 m；假設列車在運行過程中最大牽引、制動能力為1.0 m/s2，能夠允許達到的最大加速度為0.9 m/s2，最大減速度為0.75 m/s2，站間限速值為80 km/h(22.22 m/s)，最大加、減速度變化率為0.75 m/s3。

當不考慮列車運行過程中的阻力時，列車從地點C運行至tf=200 s時刻到地點D過程中的能量最優駕駛問題數學模型為

s.t.F0=0.7

FN=-0.5

Fi≥-0.75i=1,…,N-1

-Fi≥-0.9i=1,…,N-1

Fi≥-1.0i=0,1,…,N

-Fi≥-1.0i=0,1,…,N

Fi-Fi-1≥-0.75Ti=1,…,N

Fi-1-Fi≥-0.75Ti=1,…,N

i=1,…,N-1

(30)

圖2 加速度與時間曲線

圖3 速度與時間曲線

圖4 運行距離與時間曲線

圖5 速度與運行距離曲線

圖2～圖5中，列車從t0=15 s時加速度逐漸減小(大于0)，并在60.6 s時加速度為0，表示該階段為加速階段；此后加速度小于0，表示此階段為制動階段。加速度大于0時速度不斷增加，并在60.6 s時達到最高運行速度12.62 m/s(小于限速值22.22 m/s，在運行時間富裕的情況下列車并非貼限速運行)；加速度小于0時速度不斷減小直到制動運行至地點D。運行距離隨著運行時間單調上升，且運行到地點D時運行距離恰為Sf=1 900 m。運行速度與運行距離關系呈現出光滑的圓弧形曲線。

圖6 加速度與時間曲線

圖7 速度與時間曲線

圖8 運行距離與時間曲線

圖9 速度與運行距離曲線

圖6～圖9中，列車在啟動階段的最大加速度為0.76 m/s2(小于允許的最大加速度0.9 m/s2)，表明不需要以最大牽引力起車就能滿足運行要求；在44.0 ～72.9 s時間段內加速度為0，表明列車處于惰行階段(不消耗能量)，隨后制動直到運行至地點D，其中在101.7 ～114.6 s時間段內以最大減速度0.75 m/s2減速制動運行。在44.0～72.9 s的惰行時間段內，速度恰為限速值22.22 m/s。運行距離隨著運行時間單調上升，且運行到地點D時運行距離恰為Sf=1 900m。運行速度與運行距離關系呈現出光滑的圓弧形曲線。

考慮列車運行過程中的阻力時，假設列車運行過程中的坡道信息為

(31)

列車的運行阻力計算公式為

r[v(t)]=10-3gs{α0+α1[3.6v(t)]+α2[3.6v(t)]2}

α0=0.92α1=0.004 8α2=0.000 125

(32)

列車從地點C運行至tf=115 s時刻到達地點D過程中的能量最優駕駛問題數學模型為

s.t.F0-r(v0)-gs(s0)=0.7

FN-r(vN)-gs(sN)=-0.5

Fi-r(vi)-gs(si)≥-0.75i=1,…,N-1

-Fi+r(vi)+gs(si)≥-0.9i=1,…,N-1

Fi≥-1.0i=0,1,…,N

-Fi≥-1.0i=0,1,…,N

Fi-r(vi)-gs(si)-[Fi-1-r(vi-1)-gs(si-1)]≥

-0.75Ti=1,…,N

Fi-1-r(vi-1)-gs(si-1)-[Fi-r(vi)-gs(si)]≥

-0.75Ti=1,…,N

i=1,…,N-1

(33)

圖10 加速度與時間曲線

圖11 速度與時間曲線

圖12 運行距離與時間曲線

圖13 速度與運行距離曲線

圖14 加速度、牽引力、阻力與運行距離曲線

圖10中，列車從t0=15 s時加速度逐漸減小(大于0)直到在47.7 s時加速度為0，表示該階段為加速階段；并在47.7～60.0 s時間段內的大部分時間加速度為0(牽引力不為0)；在60.0～69.5 s加速度小于0，列車減速運行；在69.5～77.5 s時間段內由于受到下坡道阻力的影響加速度大于0；此后加速度小于0，表示此階段為減速制動階段。圖11中，速度不斷增加，并在47.7 s時達到限速值22.22 m/s直到60.0 s；在60.0～77.5 s時間內速度先減小后增大(小于限速值22.22 m/s)，此后開始減速制動運行至地點D。圖12中，運行距離隨著運行時間單調上升，且運行到地點D時運行距離恰為Sf=1 900 m。圖13中，運行速度與運行距離關系呈現出光滑的圓弧形曲線。圖14中，加速度是由牽引/制動力克服列車運行過程中的阻力后生成；列車在運行區間1 000～1 200 m時(30‰的大上坡道)，以較小的牽引力(小于列車阻力)沖坡運行，同時滿足列車運行過程中對速度的要求，因此在該運行區間能夠明顯體現出節能。

3.2 算法性能分析

在能量最優駕駛問題(式(33))的求解過程中，不斷迭代求解凸二次規劃問題(式(28))，隨著迭代次數的增加，式(28)的最優解將逐漸逼近式(33)的最優解。迭代過程中得到的牽引、制動力與時間，速度與時間，運行距離與時間的曲線見圖15～圖17。

圖15中，迭代過程中的牽引/制動力與時間曲線隨著迭代次數的增加逐漸趨近于一條曲線，第3次迭代和第4次迭代得到的牽引/制動力與時間曲線幾乎重合，則第4次迭代得到的牽引/制動力為式(33)的近似最優解。圖16中，迭代過程中第2次迭代到第4次迭代的速度與時間曲線幾乎重合為一條曲線。圖17中，迭代過程中的運行距離與時間曲線幾乎重合為一條曲線。

圖15 迭代過程中的牽引/制動力(單位質量下)與時間關系

圖16 迭代過程中的速度與時間關系

圖17 迭代過程中的運行距離與時間關系

4 結束語

本文根據列車的動力學模型，牽引、制動特性，線路限速，列車阻力，乘坐舒適性等約束條件，建立了列車在運行過程中能量最優控制問題的數學模型。由于考慮了列車運行過程中的阻力，約束條件中列車運行過程的動力學方程由ODEs增廣為DAEs，使得問題難以求解。為求解該模型，在時間域內將狀態變量和控制變量離散化可以將問題轉化為非線性規劃問題，采用提出的分離迭代策略，可以將該復雜的非線性規劃問題轉化為凸二次規劃問題進行迭代求解。通過對兩站間運行的列車能量最優駕駛問題進行求解計算和仿真，驗證了所提出方法的有效性。本文考慮的限速值為固定值，后續將進一步研究限速值隨線路實際位置變化的復雜情況下能量最優駕駛問題。