高超聲速飛行器自適應(yīng)高階終端滑?？刂?/h1>
2014-06-05 15:30:42田栢苓
宗 群，蘇 芮，王 婕，田栢苓
(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)
高超聲速飛行器自適應(yīng)高階終端滑?？刂?/p>
宗 群，蘇 芮，王 婕，田栢苓
(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)
針對高超聲速飛行器的縱向運(yùn)動(dòng)模型，研究了飛行器輸出跟蹤控制問題，提出了一種將動(dòng)態(tài)逆方法與高階終端滑模控制相結(jié)合的魯棒自適應(yīng)控制方法. 首先，利用反饋線性化方法對高超聲速飛行器縱向模型輸入輸出線性化；通過設(shè)計(jì)具有全局魯棒性的終端滑模面，提高系統(tǒng)的輸出收斂速度；同時(shí)，采用自適應(yīng)高階滑模控制律，在不確定上界未知條件下對其進(jìn)行自適應(yīng)估計(jì)，從而實(shí)現(xiàn)控制器增益的實(shí)時(shí)在線調(diào)整，減少系統(tǒng)抖振；最后，基于Lyapunov理論證明了此控制策略可以保證閉環(huán)系統(tǒng)穩(wěn)定. 仿真結(jié)果表明，所設(shè)計(jì)的控制器能夠?qū)崿F(xiàn)高超聲速飛行器縱向爬升機(jī)動(dòng)中速度和高度的穩(wěn)定跟蹤控制.
魯棒自適應(yīng)控制；飛行控制；高階滑模；終端滑模
高超聲速飛行器飛行環(huán)境復(fù)雜、飛行包線跨度大、動(dòng)態(tài)特性易變且采用機(jī)體/發(fā)動(dòng)機(jī)一體化設(shè)計(jì)技術(shù)，導(dǎo)致氣動(dòng)參數(shù)與發(fā)動(dòng)機(jī)推力之間存在強(qiáng)烈耦合．這些特點(diǎn)對飛行控制系統(tǒng)研究設(shè)計(jì)提出了許多新的挑戰(zhàn)，使高超聲速飛行控制成為當(dāng)前飛行器控制研究領(lǐng)域的前沿問題．
目前在控制領(lǐng)域中受到廣泛采用和認(rèn)可的模型是美國空軍實(shí)驗(yàn)室學(xué)者Bolender和Doman[1]研究的一體化解析式模型，該模型全面地刻畫了吸氣式高超聲速飛行器的縱向動(dòng)力學(xué)行為．在此基礎(chǔ)上Parker等[2]考慮高超聲速飛行器的氣推耦合和彈性振動(dòng)，建立了曲線擬合模型，利用線性二次型調(diào)節(jié)器(LQR)理論設(shè)計(jì)了一種多輸入多輸出的最優(yōu)控制器；隨后，文獻(xiàn)[3]提出增加鴨翼變量，消除控制舵面對升力的氣動(dòng)耦合項(xiàng)所導(dǎo)致的非最小相位行為，設(shè)計(jì)了基于反饋線性化方法的LQR控制器．Rehman等[4]針對反饋線性化方法應(yīng)用中存在的輸入耦合和模型不確定性問題，提出了一種由魯棒反饋線性化和最小最大LQR相結(jié)合的控制方法，但該方法無法滿足飛行器控制的實(shí)時(shí)性．
為了獲得更好的控制性能，一些學(xué)者提出Terminal滑?？刂撇呗訹5]，確?；Ｃ嫔系母櫿`差在有限時(shí)間內(nèi)收斂到零．進(jìn)一步，文獻(xiàn)[6]提出了適用于高階非線性系統(tǒng)的Terminal滑模面，克服了滑模面函數(shù)導(dǎo)數(shù)不連續(xù)的缺點(diǎn)，消除了滑?？刂频牡竭_(dá)階段，保證了系統(tǒng)的全局魯棒性．雖然Terminal滑?？刂瓶梢允瓜到y(tǒng)狀態(tài)在有限時(shí)間內(nèi)收斂，提高收斂速度和跟蹤精度，但是仍然沒有克服抖振問題．近年提出的高階滑?？刂品椒▽⒏哳l切換控制轉(zhuǎn)移到滑模變量的高階導(dǎo)數(shù)上，有效消除了抖振，同時(shí)保留著傳統(tǒng)滑模的良好特性[7-8]．此外，Su等[9]采用自適應(yīng)律對有界不確定性進(jìn)行估計(jì)，提出了自適應(yīng)滑模控制，并應(yīng)用于機(jī)械臂的跟蹤控制．
為了針對高超聲速飛行器直接設(shè)計(jì)高性能的自適應(yīng)魯棒非線性控制器，本文將高次冪多項(xiàng)式型非奇異終端滑模與高階滑?？刂品椒ㄏ嘟Y(jié)合，用于消除抖振并提高系統(tǒng)控制性能．針對飛行器輸入輸出線性化模型設(shè)計(jì)非奇異終端滑模面，同時(shí)利用自適應(yīng)方法估計(jì)有界的不確定上界，進(jìn)一步利用自適應(yīng)高階滑模控制策略，實(shí)時(shí)在線調(diào)整控制增益，使系統(tǒng)的狀態(tài)收斂至平衡點(diǎn)，消除控制抖振，對高超聲速飛行器速度與高度的階躍參考信號(hào)實(shí)現(xiàn)穩(wěn)定跟蹤.
1 模型描述
1.1 吸氣式高超聲速飛行器模型
本文考慮的吸氣式高超聲速飛行器縱向動(dòng)力學(xué)模型來源于文獻(xiàn)[1]的研究成果．基于平面地球假設(shè)，將飛行器建模為一個(gè)彈性結(jié)構(gòu)，推導(dǎo)的運(yùn)動(dòng)方程中包含彈性效應(yīng)，同時(shí)剛體動(dòng)力學(xué)和彈性動(dòng)力學(xué)之間通過氣動(dòng)力產(chǎn)生耦合．運(yùn)動(dòng)方程表達(dá)式為
式中：v為速度；h為高度；γ為航跡角；α為攻角；Q為俯仰角速率；它們?yōu)橄到y(tǒng)的5個(gè)剛體狀態(tài)；iη和 η˙i(i=1,2,3)為系統(tǒng)的6個(gè)彈性狀態(tài)．選擇速度v和高度h為系統(tǒng)輸出量，控制輸入為油門開度φ、升降舵偏轉(zhuǎn)量δe、鴨翼偏轉(zhuǎn)量δc和擴(kuò)散面積比Ad．控制輸入通過氣動(dòng)力(升力L和阻力D)、俯仰力矩M、廣義力Ni和發(fā)動(dòng)機(jī)推力T影響飛行器動(dòng)態(tài)．在標(biāo)準(zhǔn)質(zhì)量下，彈性結(jié)構(gòu)的阻尼和自然頻率分別為?m,i和ωm,i．其他變量此處不再說明，其具體含義可參見文獻(xiàn)[3，10]．
將前體轉(zhuǎn)角τ1和后體轉(zhuǎn)角τ2的變化量分別記為Δτ1和Δτ2，廣義坐標(biāo)與Δτ1、Δτ2滿足線性映射關(guān)系，該線性映射與飛行器質(zhì)量m相關(guān)．俯仰力矩M、氣動(dòng)力(升力L和阻力D)、廣義力Ni及發(fā)動(dòng)機(jī)推力T的表達(dá)式分別為
式中：CM為俯仰力矩系數(shù)；CL和CD分別為氣動(dòng)力升力系數(shù)和阻力系數(shù)；CNi為廣義力系數(shù)；CT,φ和CT為推力系數(shù)．這些系數(shù)的計(jì)算式分別為
此外，動(dòng)壓q和自由流馬赫數(shù)Ma的表達(dá)式為
式中：c為聲速；k0為比熱系數(shù)，k0=1.4；R為氣體常數(shù)，R=286.9,J/(kg·K)；T∞(h)為自由流溫度；ρ(h)為與高度相關(guān)的空氣密度．T∞(h)和ρ(h)由1976年美國標(biāo)準(zhǔn)大氣中以高度為參數(shù)查詢獲得．文中選取的物理參數(shù)如下所示：參考面積Sref=1.579,m2，平均氣動(dòng)弦長c=5.181,6,m，慣性力矩Iyy=677,909,kg/m，彈性結(jié)構(gòu)的阻尼為?m,i=0.02，重力加速度g= 9.81,m/s，自然頻率為ωm,1=21.17,rad/s，ωm,2= 53.92,rad/s，ωm,3=109.1,rad/s[11]．
1.2 高超聲速飛行器模型輸入輸出線性化
針對高超聲速飛行器模型(1)進(jìn)行輸入輸出線性化．首先，忽略彈性狀態(tài)ηi和η˙i的影響，并將鴨翼偏轉(zhuǎn)量δc設(shè)為舵偏轉(zhuǎn)量δe的函數(shù)，具體表達(dá)式為
式中kec為比例系數(shù)，k=-()，這樣即可消除
ec控制舵面對升力的氣動(dòng)耦合項(xiàng)所導(dǎo)致的非最小相位行為[4]．在此基礎(chǔ)上，進(jìn)一步忽略力和力矩系數(shù)式(3)中控制舵面對升力、阻力系數(shù)的耦合項(xiàng)以及彈性模態(tài)引起的不確定性，將式(3)簡化為
簡化后的模型只含有5個(gè)剛體狀態(tài)，即系統(tǒng)階數(shù)為5．考慮到動(dòng)態(tài)逆條件[12]，即系統(tǒng)的階數(shù)與系統(tǒng)的相對階相等，引入二階執(zhí)行器動(dòng)態(tài)使系統(tǒng)階數(shù)升為7，其具體形式為
式中φc為燃油當(dāng)量比．
經(jīng)過上述對模型的簡化與處理，下面將給出系統(tǒng)相對階的計(jì)算過程．
假設(shè)在實(shí)際飛行過程中油門全開，即擴(kuò)散面積比Ad=1，對系統(tǒng)輸出y=[vh]Τ分別求3次與4次導(dǎo)數(shù)，直至出現(xiàn)控制輸入u =[φcδe]Τ，得到
其中vf、hf、11b、12b、21b、22b的具體表達(dá)式分別為
由式(8)可知，對速度v和高度h分別微分3次和4次，控制輸入u =[φcδe]Τ出現(xiàn)在最后的微分式子中，則系統(tǒng)的相對階為7，與系統(tǒng)的階數(shù)相等，滿足動(dòng)態(tài)逆條件，系統(tǒng)的零動(dòng)態(tài)是穩(wěn)定的，因此高超聲速飛行器縱向模型可以實(shí)現(xiàn)輸入輸出線性化．令，考慮彈性飛行器模型的不確定性，相應(yīng)地方程(8)可寫為積分鏈形式，即
2 自適應(yīng)高階終端滑?？刂破髟O(shè)計(jì)
本文討論的高超聲速飛行器跟蹤控制問題是針對彈性飛行器模型(1)存在不確定的情形下，設(shè)計(jì)高階終端滑?？刂破?，使得速度v和高度h實(shí)現(xiàn)對期望速度dv和期望高度dh的穩(wěn)定跟蹤．
2.1 終端滑模面設(shè)計(jì)
定義飛行速度與高度的跟蹤誤差分別為
基于跟蹤誤差，分別設(shè)計(jì)速度滑模面vS與高度滑模面hS[12]為
式中：cv和ch為正系數(shù)向量，cv=[c11,c12,c13]，ch=
分析可知，式(19)與式(20)及后續(xù)穩(wěn)定性證明過程中需要用到ve以及he的高階導(dǎo)數(shù)，而其高階導(dǎo)數(shù)無法直接測量．本文利用具有魯棒性強(qiáng)、精度高、收斂快等優(yōu)點(diǎn)的高階滑模微分器[14]，獲取狀態(tài)ve與he的微分估計(jì)值．由于篇幅所限，這里只給出ve的高階滑模微分器的表達(dá)式為
式中：ev為待微分信號(hào)；vi(i=0,1)為ev的第i+1階微分估計(jì)值；zk(k=1,2)為系統(tǒng)(21)中的內(nèi)部狀態(tài)；λ0,λ1,λ2為微分器的設(shè)計(jì)參數(shù)．
根據(jù)文獻(xiàn)[6]的研究，Pv(t)=[pv(t),p˙v(t),˙p˙v(t)]Τ與Ph(t)=[ph(t),p˙h(t),˙p˙h(t),˙p˙˙h(t)]Τ滿足下列假設(shè).
假設(shè)1：P(t)∈C3[0,∞)，P(t)∈C4[0,∞)且
vh當(dāng)t≥t時(shí)P(t)和P(t)恒為零
0
vh (i=1,2,3;j =1,2,3,4)．
設(shè)計(jì)Pv(t)和Ph(t )的表達(dá)式分別為
根據(jù)假設(shè)1，Pv(t)和Ph(t)應(yīng)為連續(xù)光滑函數(shù)，因而當(dāng)t=t0時(shí)，Pv(t)和Ph(t)及其導(dǎo)數(shù)為零．相應(yīng)地，可求出系數(shù)aij(i=1,2,3;j =1,2,3)與bij(i=1,…,4; j=1,…,4)的具體值，具體求解過程可參考文獻(xiàn)[6]．
假設(shè)2：假設(shè)Δ1和Δ2是有界的，即max{Δ1,Δ2}≤M0，其中M0為不確定上界．
2.2 終端滑?？刂破髟O(shè)計(jì)
通過式(19)和式(20)可知，系統(tǒng)的初始狀態(tài)在滑模面上，消除了滑模到達(dá)過程，控制器設(shè)計(jì)具有全局魯棒性．
定理1：針對非線性不確定系統(tǒng)式(15)和式(16)，若假設(shè)1與假設(shè)2成立，則控制律
能夠確保Sv、Sh及誤差變量Ev與Eh快速收斂到零，其中B=[b1,b2]Τ、b1=[b11b12]、b2=[b21b22]，切換增益系數(shù)K滿足K≥max{c13M0,c24M0}．
證明：記S =[SvSh]Τ，取Lyapunov函數(shù)為
由式(19)和式(20)可知，分別對速度滑模變量Sv與高度滑模變量Sh求一階導(dǎo)數(shù)可出現(xiàn)控制量u，即
令μ1=K-c13M0，μ2=K-c24M0，則μ1和μ2均為非負(fù)值．對求取一階導(dǎo)數(shù)，并將定理1中控制律u代入，得到
設(shè)計(jì)上述的高階滑?？刂坡尚枰啦淮_定的上界，即要求M0的值已知，但實(shí)際應(yīng)用中很難預(yù)先獲得M0的確定值．因此，為滿足滑模面函數(shù)S=0的可達(dá)性，控制器的增益往往取得過大，帶有很大的保守性，容易導(dǎo)致抖振．下文將提出一種自適應(yīng)控制律的設(shè)計(jì)方法，控制增益能夠根據(jù)系統(tǒng)狀態(tài)自動(dòng)調(diào)節(jié)，不需事先知道系統(tǒng)不確定的上界．
根據(jù)式(25)設(shè)計(jì)控制律
式中?K為自適應(yīng)增益系數(shù)．
根據(jù)式(25)將集中的不確定記為
假設(shè)存在正數(shù)Kd為系統(tǒng)穩(wěn)定后K?的最終值，則Kd必須滿足Kd＞Δ．設(shè)計(jì)自適應(yīng)律為
式中：μ為自適應(yīng)調(diào)整參數(shù)，μ＞0；K?自動(dòng)調(diào)節(jié)的速率取決于μ的大?。旅鎸⑼ㄟ^Lyapunov方法證明控制律(26)的有效性．令σ=Kd-Δ，則α＞0．此外，令K?=K?-Kd．
選取Lyapunov函數(shù)為
對其求取一階導(dǎo)數(shù)為
證畢．
3 仿真分析
將上述控制器(26)在Matlab/Simulink環(huán)境中進(jìn)行仿真以驗(yàn)證其有效性．針對飛行器縱向爬升機(jī)動(dòng)過程，分別考慮飛行速度階躍變化量vΔ=30.5,m/s及飛行高度階躍變化量hΔ=610,m.其中，飛行器的初值設(shè)置及兩種測試情形下的穩(wěn)態(tài)值見表1．
表1 飛行器初始狀態(tài)及兩種情形下的穩(wěn)態(tài)值Tab.1 Initial condition of the flight vehicle and steady-state values for two different cases
選取自適應(yīng)調(diào)整參數(shù)μ=0.2，滑模面增益設(shè)為cv=[662]與ch=[6664]．跟蹤速度階躍參考信號(hào)時(shí)，飛行器的狀態(tài)(速度和高度)變化情況見圖1，控制輸入燃油當(dāng)量比cφ與升降舵偏轉(zhuǎn)量eδ的變化見圖2，自適應(yīng)增益?K的變化如圖3所示．同理，跟蹤高度階躍參考信號(hào)時(shí)，飛行器的狀態(tài)、控制輸入及自適應(yīng)增益?K的變化如圖4～圖6所示．
圖1 飛行器狀態(tài)變化情況(跟蹤30.5,m/s速度階躍指令)Fig.1 Change of vehicle state(tracking a 30.5,m/s stepvelocity command)
圖2 控制量變化情況(跟蹤30.5,m/s速度階躍指令)Fig.2 Change of the control inputs(tracking a 30.5,m/s step-velocity command)
圖3 跟蹤速度階躍指令(30.5,m/s)時(shí)的自適應(yīng)增益?KFig.3 Adaptive gain parameter ?K of tracking a 30.5,m/s step-velocity command
圖4 飛行器狀態(tài)變化情況(跟蹤610,m高度階躍指令)Fig.4 Change of the vehicle state(tracking a 610,m stepaltitude command)
圖5 跟蹤610,m高度階躍指令時(shí)的自適應(yīng)增益?KFig.5 Adaptive gain parameter ?K of tracking a 610,m step-altitude command
從仿真結(jié)果可以看出，自適應(yīng)高階終端滑?？刂票ＷC高超聲速飛行器飛行速度和高度能快速跟蹤上各自的參考指令值，并且沒有穩(wěn)態(tài)誤差，控制輸入cφ與eδ在飛行器進(jìn)入穩(wěn)定狀態(tài)后都各自趨于穩(wěn)態(tài)值．
圖6 控制量變化情況(跟蹤610,m高度階躍指令)Fig.6 Change of the control inputs(tracking a 610,m stepaltitude command)
4 結(jié) 語
本文針對高超聲速飛行器縱向運(yùn)動(dòng)學(xué)模型具有高度非線性、多變量耦合、多不確定及易受彈性影響等特點(diǎn)，提出了一種基于非線性動(dòng)態(tài)逆的自適應(yīng)高階終端滑模魯棒控制器設(shè)計(jì)方法．采用非線性動(dòng)態(tài)逆方法設(shè)計(jì)飛行器速度與高度跟蹤的高階終端滑?？刂破鳎⒉捎米赃m應(yīng)策略降低了系統(tǒng)對模型有界不確定的保守性．仿真結(jié)果表明，所提出的控制方法具有較好的跟蹤性能與魯棒性．
［1］ Bolender M A，Doman D B. Nonlinear longitudinal dynamical model of an air-breathing hypersonic vehicle [J]. Journal of Spacecraft and Rockets，2007，44(2)：374-387.
［2］ Parker J T，Serrani A，Yurkovich S，et al. Approximate feedback linearization of an air-breathing hypersonic vehicle[C]//AIAA Guidance，Navigation，and Control Conference and Exhibit. Keystone，Colorado：American Institute of Aeronautics and Astronautics，2006：1-16.
［3］ Parker J T，Serrani A，Yurkovich S，et al. Controloriented modeling of an air-breathing hypersonic vehicle[J]. Journal of Guidance，Control and Dynamics，2007，30(3)：856-869.
［4］ Rehman O U，F(xiàn)idan B，Petersen I. Uncertainty modeling and robust minimax LQR control of hypersonic flight vehicles[C]//AIAA Guidance，Navigation，and Control Conference. Toronto，Ontario，Canada：American Institute of Aeronautics and Astronautics，2010：1-15.
［5］ Wu Y Q，Yu X，Man Z H. Terminal sliding mode control design for uncertain dynamic system [J]. Systems and Control Letters，1998，34(2)：281-287.
［6］ Zhuang Kaiyu，Su Hongye，Zhang Keqin，et al. Adaptive terminal sliding mode control for high order nonlinear dynamic systems [J]. Journal of Zhejiang University Science，2003，4(1)：58-63.
［7］ Levant A. High order sliding modes，differentiation and output feedback control [J]. International Journal of Control，2003，76(9)：924-941.
［8］ Defoort M，F(xiàn)loquet T，Kokosy A，et al. A novel higher order sliding mode control scheme [J]. System & Control Letters，2009，58(2)：102-108.
［9］ Su C Y，Leung T P. A sliding mode controller with bound estimation for robot manipulators [J]. IEEE Transactions on Robotics and Automation，1993，9(2)：208-214.
［10］Wang Jie，Zong Qun，Tian Bailing，et al. Flight control for flexible hypersonic vehicle based on quasicontinuous high-order sliding mode[J]. Journal of Systems Engineering and Electronics，2013，24(2)：288-295.
［11］Sigthorsson D O，Serrani A. Development of linear parameter-varying models of hypersonic air-breathing vehicles[C]//AIAA Guidance，Navigation，and Control Conference. Chicago，Illinois：American Institute of Aeronautics and Astronautics，2009：1-21.
［12］Isidori A. Nonlinear Control Systems [M]. Berlin：Springer-Verlag，1995.
［13］Xu H J，Mirmirani M D，Ioannou P A. Adaptive sliding mode control design for hypersonic flight vehicle [J]. Journal of Guidance，Control，and Dynamics，2004，27(5)：829-838.
［14］Levant A. High-order sliding modes，differentiation and output-feedback control [J]. International Journal of Control，2003，76(9/10)：924-941.
(責(zé)任編輯：孫立華)
Adaptive High-Order Terminal Sliding Mode Control for Hypersonic Flight Vehicle
Zong Qun，Su Rui，Wang Jie，Tian Bailing
(School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China)
A robust control strategy combining adaptive high order terminal sliding mode with dynamic inversion was designed for the uncertain longitudinal dynamics of hypersonic vehicles about the output tracking problem. The feedback linearization method was adopted in this paper to linearize the input and output of a hypersonic vehicle longitudinal model. Then，the terminal sliding surface which possessed the global robustness was designed to improve the convergence rate of the output vector. Afterwards，the gain of the controller achieved real-time online adjustment by the adaptive high order sliding mode control strategy，through which the chattering was reduced. Finally，the finitetime stability of the closed-loop system was guaranteed based on the Lyapunov theory. Simulation results show that the designed controller is effective in climbing maneuvers of velocity and altitude reference commands respectively.
robust adaptive control；flight control；high-order sliding mode；terminal sliding mode
TP273
A
0493-2137(2014)11-1031-07
10.11784/tdxbz201309014
2013-09-03；修回日期：2013-11-01.
國家自然科學(xué)基金資助項(xiàng)目(61273092，61203012)；教育部科學(xué)技術(shù)研究重大資助項(xiàng)目(311012)；天津市自然科學(xué)基金資助項(xiàng)目(12JCZDJC30300)；中國航空科學(xué)基金資助項(xiàng)目(20125848004)；天津市過程檢測與控制重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(TKLPMC-201315)；天津大學(xué)自主創(chuàng)新基金資助項(xiàng)目(2013XQ-0022).
宗 群（1961— ），男，博士，教授．
宗 群，zongqun@tju.edu.cn
時(shí)間：2014-03-25.
http：//www.cnki.net/kcms/doi/10.11784/tdxbz201309014.html.
猜你喜歡
一半模型童話王國·奇妙邏輯推理(2024年5期)2024-06-19 16:03:38 Smartflower POP 一體式光伏系統(tǒng)工業(yè)設(shè)計(jì)(2022年8期)2022-09-09 07:43:20 WJ-700無人機(jī)系統(tǒng)軍民兩用技術(shù)與產(chǎn)品(2021年10期)2021-03-16 06:05:30 ZC系列無人機(jī)遙感系統(tǒng)北京測繪(2020年12期)2020-12-29 01:33:58 重要模型『一線三等角』中學(xué)生數(shù)理化·七年級(jí)數(shù)學(xué)人教版(2020年10期)2020-11-26 08:24:50 重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布數(shù)學(xué)物理學(xué)報(bào)(2020年2期)2020-06-02 11:29:24 瞞天過?！律O(shè)計(jì)萌到家藝術(shù)啟蒙(2018年7期)2018-08-23 09:14:18 連通與提升系統(tǒng)的最后一塊拼圖 Audiolab 傲立 M-DAC mini家庭影院技術(shù)(2017年9期)2017-09-26 03:41:45 設(shè)計(jì)秀海峽姐妹(2017年7期)2017-07-31 19:08:17 有種設(shè)計(jì)叫而專Coco薇(2017年5期)2017-06-05 08:53:16
宗 群，蘇 芮，王 婕，田栢苓
(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)
高超聲速飛行器自適應(yīng)高階終端滑?？刂?/p>
宗 群，蘇 芮，王 婕，田栢苓
(天津大學(xué)電氣與自動(dòng)化工程學(xué)院，天津 300072)
針對高超聲速飛行器的縱向運(yùn)動(dòng)模型，研究了飛行器輸出跟蹤控制問題，提出了一種將動(dòng)態(tài)逆方法與高階終端滑模控制相結(jié)合的魯棒自適應(yīng)控制方法. 首先，利用反饋線性化方法對高超聲速飛行器縱向模型輸入輸出線性化；通過設(shè)計(jì)具有全局魯棒性的終端滑模面，提高系統(tǒng)的輸出收斂速度；同時(shí)，采用自適應(yīng)高階滑模控制律，在不確定上界未知條件下對其進(jìn)行自適應(yīng)估計(jì)，從而實(shí)現(xiàn)控制器增益的實(shí)時(shí)在線調(diào)整，減少系統(tǒng)抖振；最后，基于Lyapunov理論證明了此控制策略可以保證閉環(huán)系統(tǒng)穩(wěn)定. 仿真結(jié)果表明，所設(shè)計(jì)的控制器能夠?qū)崿F(xiàn)高超聲速飛行器縱向爬升機(jī)動(dòng)中速度和高度的穩(wěn)定跟蹤控制.
魯棒自適應(yīng)控制；飛行控制；高階滑模；終端滑模
高超聲速飛行器飛行環(huán)境復(fù)雜、飛行包線跨度大、動(dòng)態(tài)特性易變且采用機(jī)體/發(fā)動(dòng)機(jī)一體化設(shè)計(jì)技術(shù)，導(dǎo)致氣動(dòng)參數(shù)與發(fā)動(dòng)機(jī)推力之間存在強(qiáng)烈耦合．這些特點(diǎn)對飛行控制系統(tǒng)研究設(shè)計(jì)提出了許多新的挑戰(zhàn)，使高超聲速飛行控制成為當(dāng)前飛行器控制研究領(lǐng)域的前沿問題．
目前在控制領(lǐng)域中受到廣泛采用和認(rèn)可的模型是美國空軍實(shí)驗(yàn)室學(xué)者Bolender和Doman[1]研究的一體化解析式模型，該模型全面地刻畫了吸氣式高超聲速飛行器的縱向動(dòng)力學(xué)行為．在此基礎(chǔ)上Parker等[2]考慮高超聲速飛行器的氣推耦合和彈性振動(dòng)，建立了曲線擬合模型，利用線性二次型調(diào)節(jié)器(LQR)理論設(shè)計(jì)了一種多輸入多輸出的最優(yōu)控制器；隨后，文獻(xiàn)[3]提出增加鴨翼變量，消除控制舵面對升力的氣動(dòng)耦合項(xiàng)所導(dǎo)致的非最小相位行為，設(shè)計(jì)了基于反饋線性化方法的LQR控制器．Rehman等[4]針對反饋線性化方法應(yīng)用中存在的輸入耦合和模型不確定性問題，提出了一種由魯棒反饋線性化和最小最大LQR相結(jié)合的控制方法，但該方法無法滿足飛行器控制的實(shí)時(shí)性．
為了獲得更好的控制性能，一些學(xué)者提出Terminal滑?？刂撇呗訹5]，確?；Ｃ嫔系母櫿`差在有限時(shí)間內(nèi)收斂到零．進(jìn)一步，文獻(xiàn)[6]提出了適用于高階非線性系統(tǒng)的Terminal滑模面，克服了滑模面函數(shù)導(dǎo)數(shù)不連續(xù)的缺點(diǎn)，消除了滑?？刂频牡竭_(dá)階段，保證了系統(tǒng)的全局魯棒性．雖然Terminal滑?？刂瓶梢允瓜到y(tǒng)狀態(tài)在有限時(shí)間內(nèi)收斂，提高收斂速度和跟蹤精度，但是仍然沒有克服抖振問題．近年提出的高階滑?？刂品椒▽⒏哳l切換控制轉(zhuǎn)移到滑模變量的高階導(dǎo)數(shù)上，有效消除了抖振，同時(shí)保留著傳統(tǒng)滑模的良好特性[7-8]．此外，Su等[9]采用自適應(yīng)律對有界不確定性進(jìn)行估計(jì)，提出了自適應(yīng)滑模控制，并應(yīng)用于機(jī)械臂的跟蹤控制．
為了針對高超聲速飛行器直接設(shè)計(jì)高性能的自適應(yīng)魯棒非線性控制器，本文將高次冪多項(xiàng)式型非奇異終端滑模與高階滑?？刂品椒ㄏ嘟Y(jié)合，用于消除抖振并提高系統(tǒng)控制性能．針對飛行器輸入輸出線性化模型設(shè)計(jì)非奇異終端滑模面，同時(shí)利用自適應(yīng)方法估計(jì)有界的不確定上界，進(jìn)一步利用自適應(yīng)高階滑模控制策略，實(shí)時(shí)在線調(diào)整控制增益，使系統(tǒng)的狀態(tài)收斂至平衡點(diǎn)，消除控制抖振，對高超聲速飛行器速度與高度的階躍參考信號(hào)實(shí)現(xiàn)穩(wěn)定跟蹤.
1 模型描述
1.1 吸氣式高超聲速飛行器模型
本文考慮的吸氣式高超聲速飛行器縱向動(dòng)力學(xué)模型來源于文獻(xiàn)[1]的研究成果．基于平面地球假設(shè)，將飛行器建模為一個(gè)彈性結(jié)構(gòu)，推導(dǎo)的運(yùn)動(dòng)方程中包含彈性效應(yīng)，同時(shí)剛體動(dòng)力學(xué)和彈性動(dòng)力學(xué)之間通過氣動(dòng)力產(chǎn)生耦合．運(yùn)動(dòng)方程表達(dá)式為
式中：v為速度；h為高度；γ為航跡角；α為攻角；Q為俯仰角速率；它們?yōu)橄到y(tǒng)的5個(gè)剛體狀態(tài)；iη和 η˙i(i=1,2,3)為系統(tǒng)的6個(gè)彈性狀態(tài)．選擇速度v和高度h為系統(tǒng)輸出量，控制輸入為油門開度φ、升降舵偏轉(zhuǎn)量δe、鴨翼偏轉(zhuǎn)量δc和擴(kuò)散面積比Ad．控制輸入通過氣動(dòng)力(升力L和阻力D)、俯仰力矩M、廣義力Ni和發(fā)動(dòng)機(jī)推力T影響飛行器動(dòng)態(tài)．在標(biāo)準(zhǔn)質(zhì)量下，彈性結(jié)構(gòu)的阻尼和自然頻率分別為?m,i和ωm,i．其他變量此處不再說明，其具體含義可參見文獻(xiàn)[3，10]．
將前體轉(zhuǎn)角τ1和后體轉(zhuǎn)角τ2的變化量分別記為Δτ1和Δτ2，廣義坐標(biāo)與Δτ1、Δτ2滿足線性映射關(guān)系，該線性映射與飛行器質(zhì)量m相關(guān)．俯仰力矩M、氣動(dòng)力(升力L和阻力D)、廣義力Ni及發(fā)動(dòng)機(jī)推力T的表達(dá)式分別為
式中：CM為俯仰力矩系數(shù)；CL和CD分別為氣動(dòng)力升力系數(shù)和阻力系數(shù)；CNi為廣義力系數(shù)；CT,φ和CT為推力系數(shù)．這些系數(shù)的計(jì)算式分別為
此外，動(dòng)壓q和自由流馬赫數(shù)Ma的表達(dá)式為
式中：c為聲速；k0為比熱系數(shù)，k0=1.4；R為氣體常數(shù)，R=286.9,J/(kg·K)；T∞(h)為自由流溫度；ρ(h)為與高度相關(guān)的空氣密度．T∞(h)和ρ(h)由1976年美國標(biāo)準(zhǔn)大氣中以高度為參數(shù)查詢獲得．文中選取的物理參數(shù)如下所示：參考面積Sref=1.579,m2，平均氣動(dòng)弦長c=5.181,6,m，慣性力矩Iyy=677,909,kg/m，彈性結(jié)構(gòu)的阻尼為?m,i=0.02，重力加速度g= 9.81,m/s，自然頻率為ωm,1=21.17,rad/s，ωm,2= 53.92,rad/s，ωm,3=109.1,rad/s[11]．
1.2 高超聲速飛行器模型輸入輸出線性化
針對高超聲速飛行器模型(1)進(jìn)行輸入輸出線性化．首先，忽略彈性狀態(tài)ηi和η˙i的影響，并將鴨翼偏轉(zhuǎn)量δc設(shè)為舵偏轉(zhuǎn)量δe的函數(shù)，具體表達(dá)式為
式中kec為比例系數(shù)，k=-()，這樣即可消除
ec控制舵面對升力的氣動(dòng)耦合項(xiàng)所導(dǎo)致的非最小相位行為[4]．在此基礎(chǔ)上，進(jìn)一步忽略力和力矩系數(shù)式(3)中控制舵面對升力、阻力系數(shù)的耦合項(xiàng)以及彈性模態(tài)引起的不確定性，將式(3)簡化為
簡化后的模型只含有5個(gè)剛體狀態(tài)，即系統(tǒng)階數(shù)為5．考慮到動(dòng)態(tài)逆條件[12]，即系統(tǒng)的階數(shù)與系統(tǒng)的相對階相等，引入二階執(zhí)行器動(dòng)態(tài)使系統(tǒng)階數(shù)升為7，其具體形式為
式中φc為燃油當(dāng)量比．
經(jīng)過上述對模型的簡化與處理，下面將給出系統(tǒng)相對階的計(jì)算過程．
假設(shè)在實(shí)際飛行過程中油門全開，即擴(kuò)散面積比Ad=1，對系統(tǒng)輸出y=[vh]Τ分別求3次與4次導(dǎo)數(shù)，直至出現(xiàn)控制輸入u =[φcδe]Τ，得到
其中vf、hf、11b、12b、21b、22b的具體表達(dá)式分別為
由式(8)可知，對速度v和高度h分別微分3次和4次，控制輸入u =[φcδe]Τ出現(xiàn)在最后的微分式子中，則系統(tǒng)的相對階為7，與系統(tǒng)的階數(shù)相等，滿足動(dòng)態(tài)逆條件，系統(tǒng)的零動(dòng)態(tài)是穩(wěn)定的，因此高超聲速飛行器縱向模型可以實(shí)現(xiàn)輸入輸出線性化．令，考慮彈性飛行器模型的不確定性，相應(yīng)地方程(8)可寫為積分鏈形式，即
2 自適應(yīng)高階終端滑?？刂破髟O(shè)計(jì)
本文討論的高超聲速飛行器跟蹤控制問題是針對彈性飛行器模型(1)存在不確定的情形下，設(shè)計(jì)高階終端滑?？刂破?，使得速度v和高度h實(shí)現(xiàn)對期望速度dv和期望高度dh的穩(wěn)定跟蹤．
2.1 終端滑模面設(shè)計(jì)
定義飛行速度與高度的跟蹤誤差分別為
基于跟蹤誤差，分別設(shè)計(jì)速度滑模面vS與高度滑模面hS[12]為
式中：cv和ch為正系數(shù)向量，cv=[c11,c12,c13]，ch=
分析可知，式(19)與式(20)及后續(xù)穩(wěn)定性證明過程中需要用到ve以及he的高階導(dǎo)數(shù)，而其高階導(dǎo)數(shù)無法直接測量．本文利用具有魯棒性強(qiáng)、精度高、收斂快等優(yōu)點(diǎn)的高階滑模微分器[14]，獲取狀態(tài)ve與he的微分估計(jì)值．由于篇幅所限，這里只給出ve的高階滑模微分器的表達(dá)式為
式中：ev為待微分信號(hào)；vi(i=0,1)為ev的第i+1階微分估計(jì)值；zk(k=1,2)為系統(tǒng)(21)中的內(nèi)部狀態(tài)；λ0,λ1,λ2為微分器的設(shè)計(jì)參數(shù)．
根據(jù)文獻(xiàn)[6]的研究，Pv(t)=[pv(t),p˙v(t),˙p˙v(t)]Τ與Ph(t)=[ph(t),p˙h(t),˙p˙h(t),˙p˙˙h(t)]Τ滿足下列假設(shè).
假設(shè)1：P(t)∈C3[0,∞)，P(t)∈C4[0,∞)且
vh當(dāng)t≥t時(shí)P(t)和P(t)恒為零
0
vh (i=1,2,3;j =1,2,3,4)．
設(shè)計(jì)Pv(t)和Ph(t )的表達(dá)式分別為
根據(jù)假設(shè)1，Pv(t)和Ph(t)應(yīng)為連續(xù)光滑函數(shù)，因而當(dāng)t=t0時(shí)，Pv(t)和Ph(t)及其導(dǎo)數(shù)為零．相應(yīng)地，可求出系數(shù)aij(i=1,2,3;j =1,2,3)與bij(i=1,…,4; j=1,…,4)的具體值，具體求解過程可參考文獻(xiàn)[6]．
假設(shè)2：假設(shè)Δ1和Δ2是有界的，即max{Δ1,Δ2}≤M0，其中M0為不確定上界．
2.2 終端滑?？刂破髟O(shè)計(jì)
通過式(19)和式(20)可知，系統(tǒng)的初始狀態(tài)在滑模面上，消除了滑模到達(dá)過程，控制器設(shè)計(jì)具有全局魯棒性．
定理1：針對非線性不確定系統(tǒng)式(15)和式(16)，若假設(shè)1與假設(shè)2成立，則控制律
能夠確保Sv、Sh及誤差變量Ev與Eh快速收斂到零，其中B=[b1,b2]Τ、b1=[b11b12]、b2=[b21b22]，切換增益系數(shù)K滿足K≥max{c13M0,c24M0}．
證明：記S =[SvSh]Τ，取Lyapunov函數(shù)為
由式(19)和式(20)可知，分別對速度滑模變量Sv與高度滑模變量Sh求一階導(dǎo)數(shù)可出現(xiàn)控制量u，即
令μ1=K-c13M0，μ2=K-c24M0，則μ1和μ2均為非負(fù)值．對求取一階導(dǎo)數(shù)，并將定理1中控制律u代入，得到
設(shè)計(jì)上述的高階滑?？刂坡尚枰啦淮_定的上界，即要求M0的值已知，但實(shí)際應(yīng)用中很難預(yù)先獲得M0的確定值．因此，為滿足滑模面函數(shù)S=0的可達(dá)性，控制器的增益往往取得過大，帶有很大的保守性，容易導(dǎo)致抖振．下文將提出一種自適應(yīng)控制律的設(shè)計(jì)方法，控制增益能夠根據(jù)系統(tǒng)狀態(tài)自動(dòng)調(diào)節(jié)，不需事先知道系統(tǒng)不確定的上界．
根據(jù)式(25)設(shè)計(jì)控制律
式中?K為自適應(yīng)增益系數(shù)．
根據(jù)式(25)將集中的不確定記為
假設(shè)存在正數(shù)Kd為系統(tǒng)穩(wěn)定后K?的最終值，則Kd必須滿足Kd＞Δ．設(shè)計(jì)自適應(yīng)律為
式中：μ為自適應(yīng)調(diào)整參數(shù)，μ＞0；K?自動(dòng)調(diào)節(jié)的速率取決于μ的大?。旅鎸⑼ㄟ^Lyapunov方法證明控制律(26)的有效性．令σ=Kd-Δ，則α＞0．此外，令K?=K?-Kd．
選取Lyapunov函數(shù)為
對其求取一階導(dǎo)數(shù)為
證畢．
3 仿真分析
將上述控制器(26)在Matlab/Simulink環(huán)境中進(jìn)行仿真以驗(yàn)證其有效性．針對飛行器縱向爬升機(jī)動(dòng)過程，分別考慮飛行速度階躍變化量vΔ=30.5,m/s及飛行高度階躍變化量hΔ=610,m.其中，飛行器的初值設(shè)置及兩種測試情形下的穩(wěn)態(tài)值見表1．
表1 飛行器初始狀態(tài)及兩種情形下的穩(wěn)態(tài)值Tab.1 Initial condition of the flight vehicle and steady-state values for two different cases
選取自適應(yīng)調(diào)整參數(shù)μ=0.2，滑模面增益設(shè)為cv=[662]與ch=[6664]．跟蹤速度階躍參考信號(hào)時(shí)，飛行器的狀態(tài)(速度和高度)變化情況見圖1，控制輸入燃油當(dāng)量比cφ與升降舵偏轉(zhuǎn)量eδ的變化見圖2，自適應(yīng)增益?K的變化如圖3所示．同理，跟蹤高度階躍參考信號(hào)時(shí)，飛行器的狀態(tài)、控制輸入及自適應(yīng)增益?K的變化如圖4～圖6所示．
圖1 飛行器狀態(tài)變化情況(跟蹤30.5,m/s速度階躍指令)Fig.1 Change of vehicle state(tracking a 30.5,m/s stepvelocity command)
圖2 控制量變化情況(跟蹤30.5,m/s速度階躍指令)Fig.2 Change of the control inputs(tracking a 30.5,m/s step-velocity command)
圖3 跟蹤速度階躍指令(30.5,m/s)時(shí)的自適應(yīng)增益?KFig.3 Adaptive gain parameter ?K of tracking a 30.5,m/s step-velocity command
圖4 飛行器狀態(tài)變化情況(跟蹤610,m高度階躍指令)Fig.4 Change of the vehicle state(tracking a 610,m stepaltitude command)
圖5 跟蹤610,m高度階躍指令時(shí)的自適應(yīng)增益?KFig.5 Adaptive gain parameter ?K of tracking a 610,m step-altitude command
從仿真結(jié)果可以看出，自適應(yīng)高階終端滑?？刂票ＷC高超聲速飛行器飛行速度和高度能快速跟蹤上各自的參考指令值，并且沒有穩(wěn)態(tài)誤差，控制輸入cφ與eδ在飛行器進(jìn)入穩(wěn)定狀態(tài)后都各自趨于穩(wěn)態(tài)值．
圖6 控制量變化情況(跟蹤610,m高度階躍指令)Fig.6 Change of the control inputs(tracking a 610,m stepaltitude command)
4 結(jié) 語
本文針對高超聲速飛行器縱向運(yùn)動(dòng)學(xué)模型具有高度非線性、多變量耦合、多不確定及易受彈性影響等特點(diǎn)，提出了一種基于非線性動(dòng)態(tài)逆的自適應(yīng)高階終端滑模魯棒控制器設(shè)計(jì)方法．采用非線性動(dòng)態(tài)逆方法設(shè)計(jì)飛行器速度與高度跟蹤的高階終端滑?？刂破鳎⒉捎米赃m應(yīng)策略降低了系統(tǒng)對模型有界不確定的保守性．仿真結(jié)果表明，所提出的控制方法具有較好的跟蹤性能與魯棒性．
［1］ Bolender M A，Doman D B. Nonlinear longitudinal dynamical model of an air-breathing hypersonic vehicle [J]. Journal of Spacecraft and Rockets，2007，44(2)：374-387.
［2］ Parker J T，Serrani A，Yurkovich S，et al. Approximate feedback linearization of an air-breathing hypersonic vehicle[C]//AIAA Guidance，Navigation，and Control Conference and Exhibit. Keystone，Colorado：American Institute of Aeronautics and Astronautics，2006：1-16.
［3］ Parker J T，Serrani A，Yurkovich S，et al. Controloriented modeling of an air-breathing hypersonic vehicle[J]. Journal of Guidance，Control and Dynamics，2007，30(3)：856-869.
［4］ Rehman O U，F(xiàn)idan B，Petersen I. Uncertainty modeling and robust minimax LQR control of hypersonic flight vehicles[C]//AIAA Guidance，Navigation，and Control Conference. Toronto，Ontario，Canada：American Institute of Aeronautics and Astronautics，2010：1-15.
［5］ Wu Y Q，Yu X，Man Z H. Terminal sliding mode control design for uncertain dynamic system [J]. Systems and Control Letters，1998，34(2)：281-287.
［6］ Zhuang Kaiyu，Su Hongye，Zhang Keqin，et al. Adaptive terminal sliding mode control for high order nonlinear dynamic systems [J]. Journal of Zhejiang University Science，2003，4(1)：58-63.
［7］ Levant A. High order sliding modes，differentiation and output feedback control [J]. International Journal of Control，2003，76(9)：924-941.
［8］ Defoort M，F(xiàn)loquet T，Kokosy A，et al. A novel higher order sliding mode control scheme [J]. System & Control Letters，2009，58(2)：102-108.
［9］ Su C Y，Leung T P. A sliding mode controller with bound estimation for robot manipulators [J]. IEEE Transactions on Robotics and Automation，1993，9(2)：208-214.
［10］Wang Jie，Zong Qun，Tian Bailing，et al. Flight control for flexible hypersonic vehicle based on quasicontinuous high-order sliding mode[J]. Journal of Systems Engineering and Electronics，2013，24(2)：288-295.
［11］Sigthorsson D O，Serrani A. Development of linear parameter-varying models of hypersonic air-breathing vehicles[C]//AIAA Guidance，Navigation，and Control Conference. Chicago，Illinois：American Institute of Aeronautics and Astronautics，2009：1-21.
［12］Isidori A. Nonlinear Control Systems [M]. Berlin：Springer-Verlag，1995.
［13］Xu H J，Mirmirani M D，Ioannou P A. Adaptive sliding mode control design for hypersonic flight vehicle [J]. Journal of Guidance，Control，and Dynamics，2004，27(5)：829-838.
［14］Levant A. High-order sliding modes，differentiation and output-feedback control [J]. International Journal of Control，2003，76(9/10)：924-941.
(責(zé)任編輯：孫立華)
Adaptive High-Order Terminal Sliding Mode Control for Hypersonic Flight Vehicle
Zong Qun，Su Rui，Wang Jie，Tian Bailing
(School of Electrical Engineering and Automation，Tianjin University，Tianjin 300072，China)
A robust control strategy combining adaptive high order terminal sliding mode with dynamic inversion was designed for the uncertain longitudinal dynamics of hypersonic vehicles about the output tracking problem. The feedback linearization method was adopted in this paper to linearize the input and output of a hypersonic vehicle longitudinal model. Then，the terminal sliding surface which possessed the global robustness was designed to improve the convergence rate of the output vector. Afterwards，the gain of the controller achieved real-time online adjustment by the adaptive high order sliding mode control strategy，through which the chattering was reduced. Finally，the finitetime stability of the closed-loop system was guaranteed based on the Lyapunov theory. Simulation results show that the designed controller is effective in climbing maneuvers of velocity and altitude reference commands respectively.
robust adaptive control；flight control；high-order sliding mode；terminal sliding mode
TP273
A
0493-2137(2014)11-1031-07
10.11784/tdxbz201309014
2013-09-03；修回日期：2013-11-01.
國家自然科學(xué)基金資助項(xiàng)目(61273092，61203012)；教育部科學(xué)技術(shù)研究重大資助項(xiàng)目(311012)；天津市自然科學(xué)基金資助項(xiàng)目(12JCZDJC30300)；中國航空科學(xué)基金資助項(xiàng)目(20125848004)；天津市過程檢測與控制重點(diǎn)實(shí)驗(yàn)室開放基金資助項(xiàng)目(TKLPMC-201315)；天津大學(xué)自主創(chuàng)新基金資助項(xiàng)目(2013XQ-0022).
宗 群（1961— ），男，博士，教授．
宗 群，zongqun@tju.edu.cn
時(shí)間：2014-03-25.
http：//www.cnki.net/kcms/doi/10.11784/tdxbz201309014.html.
猜你喜歡
童話王國·奇妙邏輯推理(2024年5期)2024-06-19 16:03:38
工業(yè)設(shè)計(jì)(2022年8期)2022-09-09 07:43:20
軍民兩用技術(shù)與產(chǎn)品(2021年10期)2021-03-16 06:05:30
北京測繪(2020年12期)2020-12-29 01:33:58
中學(xué)生數(shù)理化·七年級(jí)數(shù)學(xué)人教版(2020年10期)2020-11-26 08:24:50
數(shù)學(xué)物理學(xué)報(bào)(2020年2期)2020-06-02 11:29:24
藝術(shù)啟蒙(2018年7期)2018-08-23 09:14:18
家庭影院技術(shù)(2017年9期)2017-09-26 03:41:45
海峽姐妹(2017年7期)2017-07-31 19:08:17
Coco薇(2017年5期)2017-06-05 08:53:16
