柳子強(qiáng) 王曉旭 王大志
(中國第一汽車股份有限公司研發(fā)總院,長春 130013)
目前,可靠性試驗中采用的城市汽車行駛工況通常為歐洲城市循環(huán)工況,該工況加減速過程簡單,與實際車輛行駛狀態(tài)差距較大,已不能滿足車輛開發(fā)測試的需要[1]。
使用用戶大數(shù)據(jù)可以更準(zhǔn)確地描述用戶實際駕駛狀態(tài),有效提高試驗工況與用戶駕駛行為的關(guān)聯(lián)性。本文通過大數(shù)據(jù)的提取和處理,并結(jié)合主成分分析、K 均值聚類、馬爾可夫排序等數(shù)據(jù)處理方法,實現(xiàn)對用戶各種駕駛行為下典型工況的構(gòu)建,以進(jìn)一步提高試驗結(jié)果與用戶實際駕駛狀態(tài)的一致性。
本文采用的大數(shù)據(jù)主要來自某企業(yè)某車型用戶大數(shù)據(jù)平臺,基于用戶實際駕駛車輛定期反饋的數(shù)據(jù)進(jìn)行研究,用戶的選取主要參照以下原則:選取該車型用戶車輛總行駛里程由高到低排序的前10 個城市,下載相應(yīng)用戶某自然年全年的行駛數(shù)據(jù)作為分析對象。前10 個城市所有用戶的總行駛里程如表1所示。
表1 某車型用戶總行駛里程由高到低排序的前10個城市
進(jìn)行數(shù)據(jù)提取和聚合后,需要對大數(shù)據(jù)片段進(jìn)行劃分,從而形成多個完整的行駛工況。選取具有代表性的行駛工況,將這些工況的速度-時間數(shù)據(jù)按一定的時間周期劃分,可得到多個運行片段[2],片段劃分遵循如下原則:
a.將連續(xù)2 個速度為0 的時間段劃分為一個運行片段,某個片段的速度-時間曲線如圖1所示。
圖1 某單一片段速度-時間曲線
b.計算劃分后的運行片段的最大速度,若最大速度不超過5 km/h,則將此片段視為無效片段,并從片段集合中去除。
c.對剩余的有效運行片段進(jìn)行編號。
本文最終將獲得的數(shù)據(jù)劃分為104 357 個運行片段。
從對駕駛特點的描述角度,本文選取了18個特征參數(shù),如表2所示。
表2 工況特征參數(shù)
對每一個駕駛循環(huán)工況進(jìn)行18 維特征參數(shù)計算,最終形成帶有特征參數(shù)信息的片段集合,稱為特征參數(shù)矩陣。
本文采用主成分分析法將高維且具有一定相關(guān)性的復(fù)雜特征指標(biāo)轉(zhuǎn)化為低維的多個互不相關(guān)的主成分,并保留原始特征參數(shù)中的大量信息[3],為后續(xù)的分類和分析工作提供基礎(chǔ)。
假設(shè)共有n個運行片段,每個運行片段有p個特征參數(shù),將計算出的所有運行片段的特征參數(shù)矩陣記為X。設(shè)隨機(jī)變量Xi的均值為μi、協(xié)方差矩陣為Σi,則通過主成分分析,對p個特征參數(shù)進(jìn)行線性變換,生成的綜合指標(biāo)即為主成分,記為y1、y2、……yp。其中特征參數(shù)矩陣X可以表示為:
式中,xij(i=1,2,…,n,j=1,2,…,p)為特征參數(shù)矩陣中的特征參數(shù)。
經(jīng)過變換后主成分的表達(dá)式為:
式中,x1~xp為原始特征參數(shù);lij(i,j=1,2,…,p)為主成分變換后的參數(shù)矩陣;l1′~lp′為各主成分中特征參數(shù)的載荷系數(shù)。
根據(jù)相關(guān)系數(shù)矩陣或協(xié)方差矩陣求解各主成分yp,其中方差及協(xié)方差計算公式分別為:
式中,lj為協(xié)方差參數(shù)矩陣。
各特征參數(shù)的量綱不同,會引起各變量的分散程度差異較大。在通過協(xié)方差矩陣求解特征值與對應(yīng)的特征向量前,為消除量綱不同帶來的不合理影響,常對各原始變量進(jìn)行標(biāo)準(zhǔn)化處理。標(biāo)準(zhǔn)化處理后的變量為:
式中,E(xj)為X向量的平均值。
向量Z=(z1,z2,…,zp)T的協(xié)方差矩陣為相關(guān)系數(shù)矩陣ρ=[ρ(xi,xj)]p×p。主成分分析后主成分的協(xié)方差矩陣為對角矩陣Λ,其對角線元素為相關(guān)系數(shù)矩陣ρ的特征根λ1、λ2、……、λp。其中相關(guān)系數(shù)矩陣
各主成分的總累積貢獻(xiàn)量可以用r′表示:
根據(jù)式(6),第j個主成分的貢獻(xiàn)率為λj/r′,則前m個主成分的累積方差貢獻(xiàn)率為,當(dāng)累積貢獻(xiàn)率達(dá)到80%~90%時,提取前m個主成分,可以代替原始變量中大部分特征信息量。
當(dāng)特征參數(shù)變量x1~xp在某個主成分上的載荷系數(shù)近似時,對其主成分的解釋較為困難,可以通過因子分析中方差最大化正交旋轉(zhuǎn)法來實現(xiàn)對因子載荷矩陣的旋轉(zhuǎn),使每個變量在主成分上方差最大化,載荷矩陣每行或每列的元素平方向0 和1 兩級分化,因子載荷越接近1 說明相關(guān)性越強(qiáng)[4],以此實現(xiàn)對因子載荷系數(shù)的解釋。
假設(shè)共有n個特征參數(shù),主成分因子數(shù)量為m,若每次進(jìn)行2 個因子的旋轉(zhuǎn)計算,共有m(m-1)/2 種旋轉(zhuǎn)方法,將此視為一次循環(huán),直到因子載荷矩陣中各列的方差和最大且收斂。假設(shè)載荷矩陣為A:
令正交矩陣為Q:
式中,φ為矩陣正交后的旋轉(zhuǎn)角度。旋轉(zhuǎn)后的因子載荷矩陣B為:
為使旋轉(zhuǎn)后的因子載荷矩陣各列方差總和V=V1+V2最大,Vj的計算公式為:
以上計算過程均可通過SPSS軟件進(jìn)行,最終得到主成分分析的成分矩陣,將成分矩陣與原始的特征參數(shù)矩陣相乘,即可得到主成分分析結(jié)果。
通過主成分分析,可以得到降維后的結(jié)果,其總方差解釋如表3所示。
表3 總方差解釋
其中方差百分比即為此主成分所包含的原始信息的比例,累積百分比80%以上的成分即可覆蓋大部分的原始信息。本文中,前4 個主成分累積貢獻(xiàn)量為80.158%,因此可將原始數(shù)據(jù)的特征參數(shù)矩陣由18維降低到4維。
作為經(jīng)典的數(shù)據(jù)挖掘方法,K 均值聚類方法是一種無監(jiān)督的學(xué)習(xí)方法[5]。在不明確運行片段分為哪些典型工況時,運行K 均值聚類算法給定聚類數(shù)量K,即為K類典型工況。按照點與點之間的距離,將每個點分到距離最近的類簇中心所代表的類別中,所有樣本點分配完成后重新計算該類簇中所有樣本點的平均值,即為新的聚類中心點,之后繼續(xù)迭代分類,直至類簇中心點的變化很小或達(dá)到給定的迭代次數(shù)[6]。
K 均值聚類需要事先指定K值,因此需要一種可以確定最優(yōu)K值的方法,目前主要有4 種方法,即方差比準(zhǔn)則、大衛(wèi)-博丁指數(shù)、輪廓系數(shù)法、手肘法。
在實際分析中,分類數(shù)量不應(yīng)出現(xiàn)較大值,因此可在小范圍內(nèi)采用窮舉法,然后根據(jù)判定式進(jìn)行最優(yōu)解的確認(rèn),本文對分類數(shù)在2~10范圍內(nèi)進(jìn)行窮舉,綜合考慮4種方法的計算結(jié)果,確定K=4。
本文使用SPSS 軟件對主成分分析后形成的4維特征參數(shù)矩陣進(jìn)行聚類,取K=4,最大循環(huán)次數(shù)為100次,部分結(jié)果如表4所示。
表4 K均值聚類部分結(jié)果
在結(jié)果中,每個運行片段都有一個對應(yīng)的聚類編號,以及該片段到聚類中心的距離,與聚類中心距離越近,代表這一片段越能反映其所在分類的特征。因此,選取4 類中與聚類中心最近的4 個片段,即為這4類各自的代表片段,結(jié)果如表5所示。
表5 代表片段編號和距離
由于代表片段為與聚類中心距離最近的片段,因此在定義分類的含義時,可以直接觀察代表片段的特征對分類進(jìn)行解釋。綜合考慮各類代表片段的信息后,對各分類的特征解釋如下:
a. 分類1。城郊工況,主要典型特征為車速處于中高速,有較為頻繁的加減速,加速度較大,其典型片段如圖2所示。
圖2 城郊工況速度曲線
b. 分類2。高速工況,主要典型特征為車速長時處于高速段,減速不劇烈,加速度也較小,其典型片段如圖3所示。
圖3 高速工況速度曲線
c. 分類3。城市擁堵工況,主要典型特征為車速處于低速段,加速度較小,其典型片段如圖4 所示。
圖4 城市擁堵工況速度曲線
d. 分類4。城市快速路,主要典型特征為速度處于中速,有一定加減速,加速度較大,其典型片段如圖5所示。
圖5 城市快速路工況速度曲線
選取各分類中用戶行駛過程損傷較大的片段作為試驗工況,以便縮短試驗工況的時間,利用20 位用戶行駛大數(shù)據(jù)全過程的平均總損傷與試驗工況的損傷進(jìn)行協(xié)同計算,最終得出各片段循環(huán)次數(shù)如表6 所示,為了便于計算,循環(huán)次數(shù)進(jìn)行取整。
表6 各試驗片段及循環(huán)次數(shù)
試驗片段的轉(zhuǎn)速和轉(zhuǎn)矩變化情況如圖6所示。
圖6 試驗片段轉(zhuǎn)速變化情況
準(zhǔn)確復(fù)現(xiàn)用戶實際使用條件下的載荷作用效果是電驅(qū)動系統(tǒng)可靠性評價的核心,其本質(zhì)是通過合理的工況組合順序使得各部件性能同步退化到全生命周期設(shè)定目標(biāo)。本文根據(jù)馬爾可夫原理分析原始工況片段之間的順序關(guān)系。
馬爾可夫鏈實際上是一組離散隨機(jī)變量的集合[7],具體指對概率空間(Ω,F,P)內(nèi)以一維可數(shù)集為指數(shù)集的隨機(jī)變量集合X={Xn:n>0},假設(shè)隨機(jī)變量的取值范圍均在可數(shù)集內(nèi),X=si,si∈s,且隨機(jī)變量的條件概率滿足如下關(guān)系:
將式(1)中的X稱為馬爾可夫鏈,對于一個固定的馬爾可夫鏈模型,式(11)表明,隨著馬爾可夫鏈的增長,鏈中事件參數(shù)的分布不變?;隈R爾可夫鏈的這種性質(zhì),通過模型時間轉(zhuǎn)移矩陣構(gòu)建循環(huán)工況能夠代表整個原始數(shù)據(jù)中用戶實際行駛工況[8]。
采用最大似然估計法計算各類典型工況間的狀態(tài)轉(zhuǎn)移矩陣,馬爾可夫過程可以通過貝葉斯公式求得穩(wěn)態(tài)概率:
式中,P(Z0)為Z0事件的先驗概率;P(Zτ|Zτ-1)為在Zτ-1事件發(fā)生的前提下,Zτ事件發(fā)生的概率。重復(fù)上述過程,N次重復(fù)觀察的公式為:
通過極大似然函數(shù),假設(shè)從工況r轉(zhuǎn)移到工況s,可得出各工況間的狀態(tài)轉(zhuǎn)移概率方程:
式中,Prs為當(dāng)前狀態(tài)工況r轉(zhuǎn)移到下個時刻狀態(tài)工況s的概率,r,s=1,2,3,4;Nrs為工況r轉(zhuǎn)移到工況s的頻次[9]。
將聚類后的4種典型工況特征作為馬爾可夫過程的4個狀態(tài),構(gòu)成狀態(tài)空間,根據(jù)每個片段所屬的工況類別構(gòu)建馬爾可夫鏈模型。對于一個固定的馬爾可夫過程,通過各工況間的狀態(tài)轉(zhuǎn)移概率方程,即可計算各工況狀態(tài)轉(zhuǎn)移概率矩陣[10]。
本文使用nCode 進(jìn)行輔助計算,馬爾可夫轉(zhuǎn)移概率矩陣實際計算結(jié)果如表7所示。根據(jù)表7,最終的工況順序為工況3-工況4-工況2-工況1。
表7 馬爾可夫轉(zhuǎn)移概率矩陣
根據(jù)以上計算結(jié)果,本文構(gòu)建的工況為:試驗工況3 運行150 次,試驗工況4 運行100 次,試驗工況2 運行6 次,試驗工況1 運行8 次;以上工況循環(huán)100次。
按運行順序,各工況示意如圖7~圖10所示。
圖7 工況3
圖8 工況4
圖9 工況2
圖10 工況1
目前,電驅(qū)系統(tǒng)普遍使用的循環(huán)耐久試驗工況包括定轉(zhuǎn)速定轉(zhuǎn)矩運行,不同轉(zhuǎn)速、轉(zhuǎn)矩下的持續(xù)工況循環(huán)運行,以及交變工況的循環(huán)運行。采用本文提出的方法構(gòu)建的工況與當(dāng)前普遍采用的工況相比,具有以下優(yōu)點:
a.工況由用戶數(shù)據(jù)提取和處理獲得,且各工況有對應(yīng)的含義,能更準(zhǔn)確地反映用戶實際駕駛時的工況信息。
b. 由于采用了各工況中大數(shù)據(jù)計算得到的損傷值偏大的用戶數(shù)據(jù)作為代表片段,通過本文方法計算出的工況實際用時將少于目前工況用時,有利于試驗周期的壓縮。
本文通過對提取后的用戶大數(shù)據(jù)進(jìn)行主成分分析、K 均值聚類、馬爾可夫排序,構(gòu)建了基于用戶大數(shù)據(jù)的可靠性試驗工況,與當(dāng)前采用的試驗參考工況相比,更貼合用戶實際駕駛狀態(tài)。