劉良杰,馮江華,王 斌,胡云卿,黎向宇
(中車株洲電力機(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]。
列車運(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)
對時(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è)不等式約束。
當(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ì)算新的對偶間隔
假設(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é)能。
在能量最優(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)系
本文根據(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)駕駛問題。