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

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

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

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

1 列車能量最優(yōu)駕駛問題

列車運(yùn)行過程的數(shù)學(xué)模型為

(1)

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

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

(2)

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

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

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時(shí)刻列車的運(yùn)行位置、速度、加速度；S為從t0時(shí)刻起至tf時(shí)刻列車的運(yùn)行距離；Sf、Vf、Af為tf時(shí)刻列車的運(yùn)行位置、速度、加速度；t0,S0,V0,A0,tf,Sf,Vf,Af均為已知常數(shù)。

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 離散化策略

對時(shí)間區(qū)間t∈[t0,tf]進(jìn)行離散化，等分為N個(gè)子時(shí)間區(qū)間，每個(gè)子時(shí)間區(qū)間的長度為T，則有

(10)

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

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

(11)

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

s0=S0v0=V0a0=A0

(12)

sN=SfvN=VfaN=Af

(13)

(14)

設(shè)在[ti-1,ti],i=1,…,N子區(qū)間列車的運(yùn)行距離為Δsi，Δsi可計(jì)算為

?

i=3,4,…,N

(15)

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

i=2,3,…,N

(16)

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

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

(17)

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

(18)

(19)

sN=Sf≈S0+NV0T+

(20)

針對式(9)中的不等式約束，在任意時(shí)刻加速度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可轉(zhuǎn)化為

i=1,…,N-1

(24)

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

(25)

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

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)

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

2.2 一般非線性規(guī)劃問題求解

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

‖F(xiàn)(k)-F(k-1)‖≤ε

(27)

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

圖1 分離迭代策略流程圖

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

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;常數(shù)項(xiàng)q=0；AE∈R4×(N+1)，AI∈R(7N-1)×(N+1)分別為等式約束和不等式約束的Jacobi矩陣；bE∈R4,bI∈R7N-1。

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

(29)

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

Step6計(jì)算下一個(gè)迭代點(diǎn)為

Step7計(jì)算新的對偶間隔

3 算例計(jì)算與仿真

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

當(dāng)不考慮列車運(yùn)行過程中的阻力時(shí)，列車從地點(diǎn)C運(yùn)行至tf=200 s時(shí)刻到地點(diǎn)D過程中的能量最優(yōu)駕駛問題數(shù)學(xué)模型為

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 加速度與時(shí)間曲線

圖3 速度與時(shí)間曲線

圖4 運(yùn)行距離與時(shí)間曲線

圖5 速度與運(yùn)行距離曲線

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

圖6 加速度與時(shí)間曲線

圖7 速度與時(shí)間曲線

圖8 運(yùn)行距離與時(shí)間曲線

圖9 速度與運(yùn)行距離曲線

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

考慮列車運(yùn)行過程中的阻力時(shí)，假設(shè)列車運(yùn)行過程中的坡道信息為

(31)

列車的運(yùn)行阻力計(jì)算公式為

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)

列車從地點(diǎn)C運(yùn)行至tf=115 s時(shí)刻到達(dá)地點(diǎn)D過程中的能量最優(yōu)駕駛問題數(shù)學(xué)模型為

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 加速度與時(shí)間曲線

圖11 速度與時(shí)間曲線

圖12 運(yùn)行距離與時(shí)間曲線

圖13 速度與運(yùn)行距離曲線

圖14 加速度、牽引力、阻力與運(yùn)行距離曲線

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

3.2 算法性能分析

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

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

圖15 迭代過程中的牽引/制動(dòng)力(單位質(zhì)量下)與時(shí)間關(guān)系

圖16 迭代過程中的速度與時(shí)間關(guān)系

圖17 迭代過程中的運(yùn)行距離與時(shí)間關(guān)系

4 結(jié)束語

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