王 青,溫李慶,2,李江雄,柯映林,李 濤,張世炯
(1.浙江大學機械工程學院,浙江杭州310027;2.中國人民解放軍94982部隊,安徽安慶246005;3.中航工業(yè)成都飛機工業(yè)(集團)有限責任公司總裝廠,四川成都610092)
基于Petri網(wǎng)的飛機總裝配生產線建模及優(yōu)化方法
王 青1,溫李慶1,2,李江雄1,柯映林1,李 濤3,張世炯3
(1.浙江大學機械工程學院,浙江杭州310027;2.中國人民解放軍94982部隊,安徽安慶246005;3.中航工業(yè)成都飛機工業(yè)(集團)有限責任公司總裝廠,四川成都610092)
針對飛機總裝配生產線最優(yōu)化調度問題,建立能夠描述飛機脈動式總裝配生產線工藝流程間復雜關系的賦時庫所Petri網(wǎng)(TPPN)仿真模型.在飛機總裝配工序的前后繼關系和有限的車間資源約束下,以工期最短、工作組最少和工作組的總效率最高為優(yōu)化目標,提出生產線調度問題的多目標優(yōu)化模型.采用逐層法進行最優(yōu)化求解,根據(jù)關鍵路徑推算公式確定生產線的最短工期和最少工作組數(shù),采用貪心-匈牙利兩級遞階算法獲得生產線的最佳人工分配方案和最優(yōu)作業(yè)排序.通過應用實例,驗證了優(yōu)化方法的準確性和有效性.
Petri網(wǎng);飛機總裝配;生產線;最優(yōu)化調度
飛機裝配是根據(jù)飛機設計要求將飛機零件、組件、部件進行定位、連接,形成高層次裝配體或整機的復雜過程.由于飛機結構形式及系統(tǒng)組成復雜,裝配工作量巨大,約占整個飛機制造勞動量的40%~50%[1].波音公司上世紀90年代在B717總裝配中基于精益思想建立連續(xù)移動式總裝配生產線,后續(xù)推廣到B737、B777等飛機;洛克希德馬丁公司在F-22、F-35戰(zhàn)斗機總裝配中建立了移動生產線;空客公司主要以站位移動式(簡稱脈動式)總裝配生產線為主[2].
飛機總裝配生產線技術在我國飛機制造企業(yè)中的應用尚處于起步階段,在質量控制、作業(yè)安排、人員管理、物流配送等方面都不能滿足新型飛機裝配的要求.如何實現(xiàn)脈動式或連續(xù)移動式飛機總裝配生產線的合理規(guī)劃和工藝調度,是目前我國飛機制造企業(yè)迫切需要解決的問題.
飛機裝配生產線仿真建模及最優(yōu)化調度是實現(xiàn)生產線合理規(guī)劃的基礎.Ramchandani[3]在基本Petri網(wǎng)的基礎上率先提出賦時Petri網(wǎng)(timed Petri Net,TPN)理論.Rivera等[4]探討了離散制造系統(tǒng)簡化Petri網(wǎng)模型的建立方法.針對柔性制造系統(tǒng)(flexible manufacturing system,FMS)調度問題,Li等[5]提出基于TPN的同步激發(fā)策略,使用該策略能夠快速得到FMS最優(yōu)調度解.Zhang等[6]應用TPN模型解決了FMS動態(tài)調度難題.在汽車制造領域,Young等[7]采用基于TPN模型的兩級分層調度方法,獲得了汽車生產線的最優(yōu)調度方案,適用于大規(guī)模、混線生產的流水線作業(yè)模式.在飛機裝配領域,Liu等[8-9]提出面向對象分層次隨機Petri網(wǎng)(HOSPN)的建模方法,應用該方法建立了客機裝配的工藝流程模型.該方法的建模過程復雜,主要用于生產線布局的改善,未涉及資源調度的優(yōu)化問題.
本文根據(jù)某型飛機脈動式總裝配生產線的建設要求,針對生產線最優(yōu)化調度問題,建立該型飛機脈動式總裝配生產線的賦時庫所Petri網(wǎng)(timed place Petri net,TPPN)仿真模型.通過綜合考慮并合理簡化生產線的約束條件及優(yōu)化目標,提出生產線調度問題的多目標優(yōu)化模型.采用逐層求解的方法得到生產線的最優(yōu)化調度結果,將求解結果應用于飛機生產線調度.
1.1 飛機脈動式總裝配生產線簡介
某型飛機脈動式總裝配生產線總體布局如圖1所示。依據(jù)生產均衡性原則,將生產線劃分為4個站位,分別為:大部件對接站位、導管電纜安裝站位、系統(tǒng)檢測站位、飛機交付站位。飛機大部件從大部件對接站位進入生產線,最終在飛機交付站位完成所有總裝配工作。在連續(xù)生產情況下,生產線的生產周期等于各站位工作節(jié)拍的最大值.生產線站位使用工藝流程圖將站位工序及流程以圖表的形式進行說明,作為各站位調度工作的參考.大部件對接站位工藝的流程圖如圖2所示.
利用工藝流程圖可以簡潔地說明生產線需要完成的工序以及完成這些工序的大致流程,但缺乏對生產流程中工序間的前后繼約束關系及工序的并發(fā)、同步、異步特征的描述能力.
1.2 生產線TPPN模型的構建方法
飛機脈動式總裝配生產線是典型的離散事件系統(tǒng).采用賦時庫所Petri網(wǎng),即TPPN理論,可以建立考慮時間因素的離散事件系統(tǒng)的圖形化仿真模型[10].仿真模型對單項工序的描述如圖3所示.圖中,變遷ti1表示工序i的開工,ti2表示工序i的完工.對庫所pi賦予時間權值di,表示從ti1發(fā)生(pi獲得標識)開始,至少要經(jīng)過di個單位時間ti2才可以發(fā)生.
圖2 大部件對接站位工藝流程Fig.2 Process of docking station
圖3 單項工序的TPPN模型Fig.3 TPPN model for individual procedure
針對飛機脈動式總裝配生產線,可以通過如下步驟建立TPPN模型.
1)若工序i是工序j的前提工序,則在ti2和tj1之間加入庫所pij,如圖4所示,并對pij賦予時間權值dij=0.
2)對于那些無前提工序的工序,將代表它們各自開工的變遷合并成為一個,并引入初始庫所po,對po賦予時間權值do=0.
3)對于那些無后續(xù)工序的工序,把代表它們各自完工的變遷合并成一個,引入結束庫所pe,對pe賦予時間權de=0.
4)設置初始標識M0,使得M0(po)=1, M0(p)=0( p ≠po).
5)在不改變工序之間銜接關系的前提下,消去某些零權庫所(po、pe不能消去),把它們所連接的前一工序的完工變遷和后一工序的開工變遷合并成一個.最大限度消去零權庫所后,生產線每個工序由一個庫所來表示.對于化簡過程中為了不改變工序間的銜接關系而不得不保留的那些零權庫所,稱為虛擬工序庫所.
圖4 銜接工序TPPN模型Fig.4 TPPN model for linking procedure
1.3 生產線TPPN模型的運行規(guī)則
飛機脈動式總裝配生產線中的某道工序必須在所有前驅工序均已完工的情況下才能進入開工狀態(tài),生產線的運行過程是生產線工序在等待開工、加工中、完工等狀態(tài)間不斷轉變的過程.為了對生產線的運行進行仿真分析,采用生產線TPPN模型的運行規(guī)則,定義如下.
定義1 TPPN模型的使能規(guī)則[11]如下:
變遷ti∈T在標識M下使能(數(shù)學表達式為:M[ti>),當且僅當
式中:·t表示t的所有輸入庫所.
定義2 在標志M下使能的變遷t的激發(fā)將產生新標識M′[11]:
1.4 飛機脈動式總裝配生產線TPPN模型的建立
在建立飛機脈動式總裝配生產線的TPPN仿真模型時,由于生產線過于復雜且各站位的建模過程類似,本文僅詳細介紹大部件對接站位的建模過程.大部件對接站位工藝流程如圖2所示,其中0140、0150號工序的并行工序多,生產安排中存在工藝路線約束復雜、可并行作業(yè)工序多、資源調度困難等問題.0140、0150號工序的前后工序均為簡單串行工序,不存在調度難題,因此建立大部件對接站位0140、0150號工序及并行工序的TPPN仿真模型.模型中,工序編號、工序內容、工期及前驅工序如表1所示,根據(jù)表1建立的TPPN仿真模型如圖5所示.
表1 大部件對接站位TPPN模型工序信息Tab.1 Information of TPPN model for docking station
圖5 大部件對接站位重點工序TPPN仿真模型Fig.5 TPPN model of key procedure in docking station
2.1 生產線調度問題描述
飛機的總裝配從飛機大部件進入裝配線開始,需要在計劃期內完成飛機的部件裝配及功能試驗等所有工序,這些工序分屬于不同的站位.生產線調度問題主要集中在站位內工序如何安排和資源如何分配,站位間的調度更依賴經(jīng)驗豐富的工藝員進行全局規(guī)劃,因此將生產線調度問題簡化為生產線站位內調度問題.生產線站位內調度問題描述如下:m道工序需要指派給n個工作組完成;工序間存在前后繼關系約束;每道工序可由多個但僅能選其中一個工作組完成;不同的工作組完成同一道工序的效率不同;同一個工作組在某一時段只能從事一道工序;工序一旦開工便不能中斷.
調度目標為:工期最短;工作組最少;工作組的總效率最高,且逐級遞進地滿足上述3個目標.
2.2 生產線調度問題的多目標優(yōu)化模型
由問題描述可知,生產線最優(yōu)化調度問題是完全分層多目標優(yōu)化問題,數(shù)學模型如下:
式中:f1(X) 為第1優(yōu)先層,f2(X) 為第2優(yōu)先層,-f3(X) 為第3優(yōu)先層.
式中:Cmax為工期;n為工作組數(shù),其中i=1,2,…,n為工作組編號;m為工序數(shù),j=1,2,…,m為工序編號;w為工作組總效率;sj為工序j開工時間;cj為工序j收工時間;Qj為工序j的緊前工序集;pij為工作組i完成工序j所需時間;d=0,1,…,Cmax為離散時間節(jié)點,定義決策變量:
式(4)中的約束條件依次如下.工序前后繼關系約束;某時間節(jié)點一道工序只能由一個工作組開工;某時間節(jié)點一個工作組只能開工一道工序;工序加工過程連續(xù)性約束;工序的完工時間不能超過總工期.
針對式(4)的多目標優(yōu)化模型,根據(jù)模型的優(yōu)先層次逐層求解,上一優(yōu)先層的最優(yōu)解即為下一優(yōu)先層的可行域.首先,求解第1優(yōu)先層,確定生產線的最短工期;然后,基于最短工期的計算結果求解第2優(yōu)先層,確定生產線需要的最少工作組數(shù);最后,采用貪心-匈牙利兩級遞階算法求解第3優(yōu)先層.貪心算法通過合理地選取決策點,將最優(yōu)化調度問題轉換為分階段指派子問題,指派子問題利用匈牙利算法求出局部最優(yōu)解,最終組合子問題的局部最優(yōu)解得到工作組總效率最高的最優(yōu)化調度結果.
3.1 最短工期的確定
飛機脈動式總裝配生產線的TPPN模型描述了實際生產線中工序的前后繼關系及各工序耗時.根據(jù)該模型,采用關鍵路徑推算公式可以計算出各工序的最早可開工時間、完成整個生產線流程的最短工期、保證最短工期的各工序pi最晚必須開工時間.推算公式如下.
定理1 E(pi)表示工序pi的最早可開工時間,TE表示完成整個生產線流程的最短工期, L(pi)表示為保證生產線在最短工期TE內完工、工序pi最晚必須的開工時間[11].
式中:·p和p·分別表示庫所p的所有輸入與輸出變遷,pe為Petri網(wǎng)中的結束工序.
3.2 最少工作組的確定
基于最短工期的計算結果,利用定理2可以得到生產線的關鍵路徑;利用定理3可以得到最少并行作業(yè)工作組數(shù).
定理2 生產線的關鍵路徑集{pi}為滿足下列條件的所有工序[11]:
1)位于po指向pe的一條有向路上;
2)該有向路上所有的工序pi均滿足E(pi)=L(pi);
3)非開始、結束及虛擬工序.
定理3 設∑為某裝配生產線的TPPN模型.P(ti,tj)是∑中的一個并發(fā)段,MTL(∑ )是∑的一條關鍵路徑,P(ti,tj)表示并發(fā)段P(ti,tj)中全部庫所的集合,則可得下式[11].
1)并發(fā)段P(ti,tj)的最短執(zhí)行時間為
2)為了保證執(zhí)行時間最短,該并發(fā)段所需的工作組的最小數(shù)目為
3.3 最優(yōu)化調度問題的分解和轉換
求解最優(yōu)化問題的算法通常需要經(jīng)過一系列的步驟,在每個步驟都面臨多種選擇.貪心算法在問題的每一步都作出當前最佳的選擇,即得到問題的局部最優(yōu)解,并組合局部最優(yōu)解得到問題的全局最優(yōu)解[12].
深入分析飛機脈動式總裝配生產線可知,調度過程存在以下特點:隨著時間的推移,生產線不斷有工序完工,從而不斷有工作組進入空閑等待派工、可開工工序進入等待指派狀態(tài).最優(yōu)化調度是在時間推移中做出最優(yōu)指派的決策.貪心算法適用于飛機脈動式總裝配生產線最優(yōu)化調度問題,貪心-匈牙利兩級遞階算法描述如下.
算法1 貪心-匈牙利兩級遞階算法
1)初始化.初始化時間t=0,輸入矩陣I、輸出矩陣O,工序時延集D,已開工工序剩下時延集d,已完工工序標識集W,新完工工序標識集w,已開工工序標識集M,待開工工序標識集N,激發(fā)的變遷集δ,表示指派結果優(yōu)劣的全局誤工系數(shù)cost=0.
2)時間推移.由已開工工序集剩下時延d(pi)、待開工工序集最晚必須開工時間L(pj)和時間t得到變遷激發(fā)時延α=min{{d(pi)}∪{L(pj)-t}},時間推移結果為t=t+α.由時延公式可知,時間推移的結果為產生新完工工序或產生急需派工工序.
3)判斷新完工工序標識集w是否為空.若是,則表示未產生新完工工序,轉4)進一步判斷;若否,則表示產生新完工工序,轉5).
4)判斷待開工工序標識集N是否為空.若是,則無等待派工工序,轉5);若否,則表示有緊急工序等待指派,此為分階段指派決策點,轉10)進行派工.
5)調用激發(fā)函數(shù).
6)判斷激發(fā)函數(shù)返回的激發(fā)變遷集δ末位否為1.若是,則表示仿真到達結束工序,調度結束,輸出結果;若否,則轉7).
7)判斷δ是否為空.若是,則無變遷激發(fā),轉9);若否,變遷激發(fā)將產生新的待開工工序,更新已完工工序標識集W、待開工工序標識集N,轉8).
8)判斷待開工工序標識集N中是否有虛擬工序.虛擬工序D(di)=0,不參與指派.若是,則令W(di)=1,N(di)=1,轉5);若否,則轉9).
9)判斷已開工工序標識集M是否為空.若是,則工作組均處于空閑等待派工狀態(tài),此為分階段指派決策點,轉10)進行派工;若否,則轉2).
10)調用指派函數(shù).指派函數(shù)面向待開工工序標識集N對空閑工作組進行派工,并將分階段派工的局部誤工系數(shù)以迭加的形式計入全局誤工系數(shù)cost.求解指派問題的匈牙利算法在3.4節(jié)中闡述.
算法1的流程如圖6所示.
圖6 貪心算法流程圖Fig.6 Flow chart of greedy algorithm
3.4 指派問題的求解
貪心算法將最優(yōu)化調度問題轉換為分階段指派子問題,在分階段指派決策點遇到如下子問題的優(yōu)化:指派n個工作組完成m道工序,已知工作組i完成工序j的效率為cij,各工作組完成一道不同的工序,如何優(yōu)化安排工作組從事的工序才能使得整體工作效率最高?該問題為指派問題,也稱為最佳匹配問題.匈牙利算法是解決該類問題的有效算法,該算法的核心是尋找增廣路徑,用增廣路徑求取二分圖的最佳匹配[13].
匈牙利算法解決的是效率矩陣為方陣、即工作組數(shù)和工序數(shù)相等的指派問題,對于效率矩陣不是方陣的指派問題模型,存在以下幾種情況及對應的解決辦法.
1)工作組多工序少.可以添加一些虛擬的工序,各工作組完成這些工序的效率可以取為0.
2)工作組少工序多.可以添加一些虛擬的工作組,他們完成各工序的效率可以取為0.
通過上述方法添加虛擬工作組或虛擬工序,將效率矩陣為非方陣的指派問題轉化成效率矩陣為方陣的指派問題,且轉換后不影響指派結果.
匈牙利算法只能求解最小化優(yōu)化問題,因此引入誤工系數(shù),將最大效率對應最小誤工系數(shù),即由效率矩陣C得到誤工系數(shù)矩陣B,將最大化效率問題轉換為最小化誤工系數(shù)問題.轉換方法為在效率矩陣C中找出最大數(shù)m=max{cij},然后進行矩陣變換B=m I-C.通過轉換,將原最大化效率指派問題模型轉換為同解的最小化誤工系數(shù)指派問題模型:
再利用匈牙利算法求解.
以某型飛機總裝配生產線為例,由于站位工藝的分組及調度主要依賴個人經(jīng)驗,易于出現(xiàn)人員管理混亂、裝配作業(yè)交叉、站位計劃拖延等問題,具體而言:一、二站位經(jīng)常出現(xiàn)工作節(jié)拍超過計劃節(jié)拍,飛機交付拖延;三、四站位可以按計劃節(jié)拍生產,但工作組空閑等待情況嚴重,資源利用率不高.
本文生產線TPPN仿真模型的運行規(guī)則及多目標優(yōu)化模型的解法采用Matlab編程實現(xiàn),應用本文算法已較好地解決該生產線存在的問題.限于篇幅,本文僅詳述第一站位(大部件對接站位)的求解過程,其他站位的求解過程類似,僅列出求解結果.
4.1 第一站位求解結果與分析
由表1、圖5得到第一站位調度問題的初始條件,運行最短工期函數(shù),得到如表2所示的結果,取最少工作組數(shù)完成裝配任務.
表2 第一站位最短工期函數(shù)運行結果Tab.2 Results of shortest duration function for first station
所選裝配任務除去開始、結束及虛擬工序,總計有12道工序指派給3個工作組.12道工序按工作內容大致區(qū)分為如表3所示的3大類.從專業(yè)化生產的角度出發(fā),規(guī)定每個工作組優(yōu)先指派某一類工序,以效率來表明工作組對該類工序的熟悉程度.效率權重值為1,數(shù)值越大,熟悉度越高,如表4所示.
表3 第一站位工序分類Tab.3 Classification of procedures in first station
表4 第一站位工作組效率Tab.4 Efficiency of working groups in first station
根據(jù)表4、式(11),可得最小化指派問題的誤工系數(shù)矩陣.
根據(jù)以上初始條件,運用本文算法得到最優(yōu)化調度結果,優(yōu)化調度前、后的甘特圖如圖7所示.將調度結果用于指導第一站位調度,得到應用最優(yōu)化調度方法前、后該站位各參數(shù)的對比,如表5所示.
圖7 第一站位復雜工序優(yōu)化前、后對比圖Fig.7 Comparison of Gantt chart for complex process in first station before and after using optimal method
表5 第一站位應用最優(yōu)化方法前、后各參數(shù)對比Tab.5 Comparison of various parameters of first station before and after using optimal method
4.2 其他站位求解結果
分別建立導管電纜安裝站位、系統(tǒng)檢測站位、飛機交付站位重點工序的TPPN仿真模型、多目標優(yōu)化模型并調用多目標優(yōu)化模型的求解函數(shù),得到各站位重點工序的最優(yōu)化調度結果.將調度結果用于指導各站位的調度,得到應用最優(yōu)化調度方法前、后各站位相關參數(shù)對比如表6所示.
表6 各站位應用最優(yōu)化調度方法前、后相關參數(shù)對比Tab.6 Comparison of various parameters of other stations before and after using optimal method
由表5、6可知,通過應用本文提出的最優(yōu)化調度方法,使得某型飛機生產線的生產周期由41 d減少至30 d,工作組數(shù)由12組減少至10組,且工作組均被指派從事其最擅長的工作,生產效率得到顯著提高.
(1)本文所提出的生產線工藝流程TPPN仿真模型體現(xiàn)了飛機生產線離散狀態(tài)、事件驅動的特性,描述了生產線工序間的前后繼約束關系及工序的并發(fā)、同步、異步特征.
(2)飛機生產線調度問題的多目標優(yōu)化模型是對復雜生產線調度問題的合理簡化.通過逐層求解多目標優(yōu)化模型,能夠迅速地得到滿足需求的最優(yōu)解.
(3)通過應用上述建模及優(yōu)化方法,使飛機生產線的生產周期縮短,工作組數(shù)量減少且均被指派其最擅長的工作,生產效率得到顯著提高.
[1]王云渤,張關康,馮宗律,等.飛機裝配工藝學(修訂本)[M].北京:國防工業(yè)出版社,1990.
[2]范玉青,梅中義,陶劍,等.大型飛機數(shù)字化制造工程[M].北京:航空工業(yè)出版社,2011.
[3]RAMCHANDANI C.Analysis of asynchronous concurrent systems by timed Petri nets[D].Cambridge:MIT,1974.
[4]RIVERA R I,RAMIREZ T A,LOPEZ M E.Building reduced petri net models of discrete manufacturing systems[J].Manufacturing and Computer Modeling,2005, 41(8):923- 937.
[5]LI Cheng,WU Wei-min,RONG Gang.Heuristic search and concurrency strategy based on Petri net for FMS scheduling[J].Networking,Sensing and Control, 2014,64(5):80- 85.
[6]ZHANG Wei-jun,THEODOR F,YANG Hua-shu.Dy-namic scheduling in flexible assembly system based on timed Petri nets model[J].Robotics and Computer-Integrated Manufacturing.2005,21(6):550- 558.
[7]YOUNG W K,INABA A,SUZUKI T,et al.Hierarchical scheduling for large-scale production system based on continuous and timed Petri net model[C]//Proceedings of the 41st SICE Annual Conference.Piscataway:IEEE,2002:268- 271.
[8]LIU Xia,YE Wen-hua,WEI Bi-sheng,et al.Research on multi-level modeling method for aircraft assembly line[J].Advanced Materials Research,2012,490- 495:538- 542.
[9]LU Hu,LIU Xia,PANG Wei,et al.Modeling and simulation of aircraft assembly line based on quest[J].Advanced Materials Research,2012,569(2):666- 669.
[10]GRADISAR D,MUSIC G.Production-process modeling based on production-management data:a Petri-net approach[J].International Journal of Computer Integrated Manufacturing,2007,20(8):794- 810.
[11]吳哲暉.Petri網(wǎng)導論[M].北京:機械工業(yè)出版社, 2006.
[12]CORMEN T H,LEISERSON C E,RIVEST R L,et al.Introduction to algorithms[M].Cambridge:MIT, 2013.
[13]王樹禾.圖論[M].北京:科學出版社,2004.
Modeling and optimization for aircraft final assembly line based on Petri net
WANG Qing1,WEN Li-qing1,2,LI Jiang-xiong1,KE Ying-lin1,LI Tao3, ZHANG Shi-jiong3
(1.School of Mechanical Engineering,Zhejiang University,Hangzhou 310027,China;2.94982 Troops,PLA,An'qing 246005,China;3.Assembly Plant,AVIC Chengdu Aircraft Corporation,Chengdu 610092,China)
A timed place Petri Net(TPPN)model for pulse final assembly line was constructed to describe the complex logical relationship of the assembly process in order to optimize the scheduling of aircraft final assembly line.A multi-objective optimization model was proposed under the constraints of limited workshop resources and the successive relationship in aircraft assembly process in order to achieve the following goals:the shortest duration,the minimum working groups,and the highest overall efficiency of working groups.A hierarchical solving strategy was adopted to obtain the optimal solution.The formula of critical path method(CPM)was used to figure out the shortest duration and the minimum working groups.Then the greedy Hungary algorithm(GHA)was applied to obtain the reasonable arrangement of human resources and optimal job sorting.An example of aircraft final assembly line was provided to verify the accuracy and effectiveness of the simulation and optimization method.
Petri net;aircraft final assembly;production line;optimal scheduling
10.3785/j.issn.1008-973X.2015.07.004
V 260
A
1008- 973X(2015)07- 1224- 08
2014- 10- 17. 浙江大學學報(工學版)網(wǎng)址:www.journals.zju.edu.cn/eng
國家自然科學基金資助項目(51375442).
王青(1979-),男,副教授,從事飛機數(shù)字化裝配技術及系統(tǒng)集成的研究.ORCID:0000-0003-4318-7867.
E-mail:wqing@zju.edu.cn