王永坤,孫明瑋,劉忠信,陳增強
(南開大學 自動化系,300071天津)
存在撓性結構或包含撓性部件的一類系統(tǒng)稱為撓性系統(tǒng).撓性系統(tǒng)中的撓性部件分為很多種,例如衛(wèi)星的太陽能帆板[1-3]、機器人的機械臂[4]、大型運載火箭或具有一定范圍長細比的有翼飛行器等.此類系統(tǒng)的一個主要特點是系統(tǒng)中存在弱阻尼的諧振模態(tài),在控制系統(tǒng)設計時需要特殊考慮,以保證閉環(huán)系統(tǒng)穩(wěn)定[5].撓性系統(tǒng)振動模態(tài)諧振頻率的理想值或者是通過理論方法計算出來的,或者是通過地面振動試驗測量出來的.理論計算中分布質量和結構特性存在的誤差,地面試驗中空天的不一致性以及測量誤差等,都會導致諧振頻率誤差.飛行器在飛行過程中,振動模態(tài)的一些特征參數(shù)會隨著溫度和質量的變化而變化,因此通過理論計算或者模態(tài)振動試驗得到的諧振頻率具有一定的不確定性,在控制系統(tǒng)的設計過程中必須考慮其對于閉環(huán)魯棒穩(wěn)定性的影響.文獻[6-7]指出采用H∞回路成形法[8-10]設計撓性控制系統(tǒng)時,由于系統(tǒng)中弱阻尼模態(tài)的存在,會減小諧振頻率的實際攝動范圍,并不能保證系統(tǒng)具有期望的魯棒性,但是并沒有給出計算諧振頻率實際攝動范圍的方法,同時文獻[7]計算出的系統(tǒng)H∞范數(shù)的逆并不代表系統(tǒng)中諧振頻率攝動值的百分比,它實際是系統(tǒng)綜合攝動范圍的一個指標.由于諧振頻率分別出現(xiàn)在一次項和二次項中,導致其非線性的耦合攝動形式,傳統(tǒng)的線性或雙線性分析方法不再適用[11-12].目前還沒有發(fā)現(xiàn)文獻報道存在系統(tǒng)性的計算方法可以直接確定撓性系統(tǒng)諧振頻率的攝動范圍.針對這個問題,工程上只能采用攝動試驗法,通過試湊的方式摸索得到諧振頻率的穩(wěn)定性邊界,該方法不僅費時費力,而且精度難以得到保證.因此,本文給出一種基于D分割法[13]的系統(tǒng)魯棒穩(wěn)定性分析方法,給出弱阻尼模態(tài)下諧振頻率攝動范圍的幾何求法,并通過一個算例驗證了該方法的有效性,同時通過另一算例說明了本方法可直接推廣應用到具有多個模態(tài)的撓性系統(tǒng)分析中.
本文以具有太陽能帆板的衛(wèi)星為例來分析撓性系統(tǒng)諧振頻率的攝動范圍,但文中方法也可直接推廣到有翼飛行器等撓性對象的分析中.一般來說,撓性結構振動具有無窮多個模態(tài),但是振動能量主要集中在低階項上,因此撓性系統(tǒng)設計時常以一階模態(tài)為主,將高階模態(tài)作為未建模動態(tài)來處理.在控制系統(tǒng)設計時要求控制器的增益在一階模態(tài)頻率后迅速衰減.圖1給出了衛(wèi)星示意圖.在圖1(a)中,θ1是衛(wèi)星主體與星體的夾角,θ2是衛(wèi)星姿態(tài)角,代表星體傳感器與儀表設備的夾角,圖1(b)所示為將傳感器安裝到與θ2相連的圓盤上的衛(wèi)星等效機械系統(tǒng)圖.對圖1所示模型,假設用彈簧鏈接的兩個質體的扭矩常數(shù)為k,粘性阻尼系數(shù)為b,建立線性模型如下[4]:
其中Tc為作用在主體上的控制力矩,J1、J2分別是圖1(b)中兩個質體的轉動慣量.由式(1)可得含有一維振動模態(tài)的撓性系統(tǒng)模型為
其中
式中c0、c1是模態(tài)幅值,只考慮一階模態(tài)時c0=-c1,ωn和ξ分別是系統(tǒng)的諧振頻率和阻尼比.在衛(wèi)星飛行過程中,參數(shù)k會隨溫度的波動以及外形結構特性而變化,因此由式(3)可知諧振頻率ωn會發(fā)生攝動.從式(2)可以看出ωn分別出現(xiàn)在一次項和二次項中,導致其攝動形式是非線性的,其對于閉環(huán)控制系統(tǒng)的影響不是簡單的線性或者雙線性形式,傳統(tǒng)的分析方法無法直接應用[11-12].
圖1 衛(wèi)星及其雙體
考慮如下一類反饋控制系統(tǒng)的閉環(huán)特征多項式:
其中ak(z)和bk(z)分別為參數(shù)向量z的連續(xù)函數(shù),z=[z1,z2,…,zl]Τ.L為系統(tǒng)的時滯常數(shù),L∈Η?R,zi∈Q?RL,R為實數(shù)集.
對于給定l階向量z和L,若特征多項式δ(s,L,z)無非負實部零點,則稱其漸近穩(wěn)定.由于δ(s,L,z)的零點為z和L的連續(xù)函數(shù),空間Η×Q被超曲面分割成若干個區(qū)域,超曲面上的點對應的δ(s,L,z)至少有一個純虛根或s=∞,這種分解法稱為Η×Q的D分割法.D分割法的實質是將s平面上的虛軸映射為參數(shù)空間Η×Q的超曲面,根據(jù)實數(shù)根穿越原點和穿越無窮,復根穿越虛軸得出D分割邊界[14].
其中
顯然,隨著z和L的連續(xù)變化,穩(wěn)定系統(tǒng)變化到不穩(wěn)定系統(tǒng)中間必然穿越?D.換句話說,研究魯棒穩(wěn)定邊界只需要通過分析?D就可以實現(xiàn).
由于本文所研究問題的時滯特性可以忽略,因此考慮如下特征多項式:
其中ai(i=0,1,…,n)是多項式系數(shù),假設ai(i=0,1,…,n)均可由特定參數(shù)d的線性函數(shù)表示,d是式(5)中的不確定參數(shù),則上式可寫為
其中p(s)和q(s)是兩個互質的多項式.假設當式(6)中的d=0時,Pc(s)是穩(wěn)定的.
此時,基于D分割法,魯棒穩(wěn)定邊界?D上的p(s)和q(s)分別是
其中f(ω)、e(ω)、h(ω)、g(ω)均為實數(shù)多項式.由Pc(jω)=0有
消去上式中d,可得
求解式(8)得到的正實根ωi,表示系統(tǒng)的截止頻率.將ωi代入式(7)可得
定義
其中d+,d-即為在滿足系統(tǒng)魯棒穩(wěn)定前提下,不確定參數(shù)d的正負攝動邊界.
當式(5)的某些系數(shù)中存在d1,d2兩個不確定參數(shù),并且這些系數(shù)均可表示為d1,d2的雙線性函數(shù)形式時,就成為一個二維線性不確定問題.假設當d1=d2=0時,Pc(s)是穩(wěn)定的.為了確定d1、d2的魯棒穩(wěn)定區(qū)域,首先,令d2=0,采用2.2節(jié)的方法求得此時參數(shù)d1的魯棒穩(wěn)定區(qū)域為然后,在區(qū)域內均勻采樣N個點,對于其中每一個d1,i(i=1,2,…,N),求解此時關于不確定參數(shù)d2的魯棒穩(wěn)定區(qū)域;最后,分別連接點和和其中i=1,2,…,N-1,不確定參數(shù)的魯棒穩(wěn)定區(qū)域如圖2陰影部分所示.對于多維線性不確定系統(tǒng),其魯棒性分析方法類似.對于存在二維或者多維不確定參數(shù)的系統(tǒng),當各個不確定參數(shù)之間存在一定的關系,而參數(shù)之間的關系曲線與上述得到的魯棒穩(wěn)定區(qū)域的交點即為各個參數(shù)的攝動邊界.而攝動點一般在攝動邊界內選取是有意義的.
圖2 不確定參數(shù)魯棒穩(wěn)定區(qū)域示意圖
考慮式(2)所表示的撓性系統(tǒng),系統(tǒng)的閉環(huán)特征多項式為
其中C(s)是控制器,諧振頻率ωn是攝動因子,分別出現(xiàn)在一次項和二次項中.
現(xiàn)在考慮一個更廣義的二維不確定性系統(tǒng).將式(9)中ωn的一次項和二次項分別表示為
其中ωn和ωn2分別表示諧振頻率攝動后的一次項和二次項,n是諧振頻率的標稱值,d1和d2分別表示攝動的百分比,此時存在d1和d2兩個不確定參數(shù),轉化為二維不確定性問題.采用2.3節(jié)的方法,假設在給定n、ξ、c0、c1和C(s)的情況下,繪制出不確定參數(shù)魯棒穩(wěn)定區(qū)域即圖3中的陰影部分.
圖3 魯棒穩(wěn)定參數(shù)區(qū)域
對于包含二階振動模態(tài)的撓性系統(tǒng),即存在ωn1和ωn2兩個不確定參數(shù).首先令ωn2固定為標稱值,通過上述方法求得ωn1的攝動區(qū)間.然后在ωn1的攝動區(qū)間內均勻采樣,對于每一個采樣值,可以求出ωn2的攝動邊界值,最終得到[ωn1,ωn2]的整體攝動區(qū)域.對于包含多階振動模態(tài)的撓性系統(tǒng),分析方法類似.
例1 考慮具有兩個太陽能帆板的衛(wèi)星模型[7],其中c0=1.731 9×10-5,c1=3.785 9×10-4,ωn1=1.539 rad/s,ξ1=0.003,此類系統(tǒng)具有明顯的撓性,一階模態(tài)的c1大于剛性模態(tài)的c0.
這里直接引用文獻[7]中的控制器設計結果,通過回路成形法得到H∞控制器.首先給被控對象加入串聯(lián)補償環(huán)節(jié),使其具有期望的開環(huán)特性,即回路成形.然后針對補償后的廣義對象,設計H∞控制器,保證閉環(huán)系統(tǒng)的穩(wěn)定性.其中補償函數(shù)選擇為
由式(2)、(10)和(11)可得廣義的系統(tǒng)傳遞函數(shù)為
當d1=d2=0時,針對標稱系統(tǒng)設計的H∞控制器為
因此可得閉環(huán)特征多項式
令d2=0,根據(jù)2.2節(jié)一維不確定參數(shù)分析方法可以計算出,d1的允許攝動區(qū)間為[-0.611 4,0.312 3].由于直線d1=d2只通過第一、三象限,因此針對每一個采樣的d1,i;當d1,i≥0時,只需求出此時對應的d2攝動的上邊界值當d1,i<0時,只需求出此時對應的d2攝動的下邊界值,結果如圖4所示.直線d1=d2與一、三象限邊界的交點分別為(0.138 6,0.138 6)和(-0.351 4,-0.351 4),表明諧振頻率ωn1增加13.86%或減小35.14%后,系統(tǒng)處于臨界穩(wěn)定狀態(tài),ωn1的攝動范圍即為[-35.14%,13.86%].
圖4 諧振頻率ωn1的攝動范圍
攝動系統(tǒng)的階躍響應曲線如圖5所示,其中實線表示標稱系統(tǒng)的響應曲線,短虛線和長虛線分別代表諧振頻率攝動-35.14%和13.86%時的響應曲線.可以看出當ωn1的攝動達到所求的正負邊界時,系統(tǒng)出現(xiàn)等幅振蕩,驗證了本文提出方法的正確性.
圖5 系統(tǒng)的階躍響應曲線
而采用文獻[7]的分析方法,求得系統(tǒng)的H∞范數(shù)γ=2.486,則γ-1≈0.402 3,根據(jù)小增益定理可得只要對象的攝動小于40.23%,就可以保證閉環(huán)魯棒穩(wěn)定.而實際情況是當諧振頻率正攝動14%時就達到了穩(wěn)定邊界,文獻[7]得出的結論是:此情況的出現(xiàn)是由于對象存在弱阻尼模態(tài),壓縮了實際允許的攝動范圍,因此導致H∞回路成形法設計的系統(tǒng)并不一定具有魯棒性.
本文經(jīng)過分析,γ-1≈0.402 3表示當系統(tǒng)中的不確定參數(shù)ωn1攝動值達到最大時,對應整體不確定項的H∞范數(shù)為0.402 3.定義
其中G和G0分別表示實際對象和標稱對象,ΔG表示對象的整體不確定項,其H∞范數(shù)與Δωn1的關系如圖6所示,其中橫坐標表示諧振頻率的攝動值,縱坐標表示‖ΔG‖∞.從圖中可以看出當諧振頻率ωn1攝動達到正向邊界即13.86%時,‖ΔG‖∞為0.407 5,幾乎等于文獻[7]給出的0.402 3.因此,系統(tǒng)H∞范數(shù)的逆并不表示某一參數(shù)的攝動范圍,而是一種結構參數(shù)攝動引起的非結構不確定性變化范圍,二者本質上是不同的.
圖6 整體不確定項的H∞范數(shù)
例2 仍采用算例1的模型以及控制器設計結果,并在此基礎上加入一個二階振動模態(tài),得到如下包含兩個振動模態(tài)的撓性系統(tǒng):
其中c0、c1、ωn1、ξ1保持例1中的數(shù)值不變,c2=1.02×10-4,ωn2=5.2,ξ2=0.002,諧振頻率ωn1和ωn2是系統(tǒng)中的不確定參數(shù).
首先令ωn2固定為標稱值,對ωn1進行攝動分析.通過本文提出的方法可得ωn1的攝動比例為d1∈[-0.36,0.163].對比例1可見,由于二階模態(tài)的能量比較低,對于一階模態(tài)振動頻率的攝動范圍影響不是很大,因此一般在對撓性系統(tǒng)進行分析時,通常只需考慮能量較大的低階振動模態(tài)即可.當ωn1的攝動達到臨界值時,系統(tǒng)的階躍響應曲線如圖7所示出現(xiàn)等幅振蕩.
在區(qū)間[-0.36,0.163]內均勻采樣,對于每一個對應的ωn1值,求得ωn2攝動比例的邊界值,最終可得ωn1、ωn2的攝動魯棒穩(wěn)定區(qū)域,即圖8中曲線與坐標軸所圍區(qū)域,曲線表示攝動的邊界.需要說明的是二階模態(tài)諧振頻率ωn2的攝動沒有上界,反映了諧振頻率越高對于整個剛體的影響越小,符合實際經(jīng)驗.在求得區(qū)域選取內點(-0.2,-0.5),階躍響應如圖9(a)所示,可以看出此時系統(tǒng)處于收斂的趨勢;在邊界上選取(-0.32,-0.713)和(0.07,-0.793),階躍響應如圖9(b)所示,輸出產(chǎn)生等幅振蕩.
圖7 ωn 攝動時系統(tǒng)的階躍響應曲線
圖8 諧振頻率ωn1、ωn2的攝動魯棒穩(wěn)定范圍
圖9 ωn1、ωn2攝動時系統(tǒng)的階躍響應曲線
本文針對撓性系統(tǒng)中不確定諧振頻率的魯棒穩(wěn)定性分析問題,設計了一種基于D分割法的幾何分析方法,給出了振動模態(tài)諧振頻率的精確攝動范圍.通過數(shù)值算例驗證了該方法的有效性,對于實際工程具有一定參考價值.該方法由于突破了傳統(tǒng)魯棒穩(wěn)定性分析要求攝動參數(shù)必須是線性或者雙線性的形式限制,也可推廣到其他含參數(shù)非線性攝動形式的場合,通過繪制相應的幾何曲線求取特定參數(shù)的魯棒穩(wěn)定性范圍.
[1]MCFARLANE D C,GLOVER K.Robust controller design using normalized coprime factor plant descriptions[M].New York:Springer Verlag,1990:205-208.
[2]劉暾,楊大明.撓性衛(wèi)星動力學及姿態(tài)控制模型的建立[J].哈爾濱工業(yè)大學學報,1985,17(A8):1-17.
[3]崔乃剛,劉家夫,榮思遠.柔性太陽帆航天器動力學建模與姿態(tài)控制[J].哈爾濱工業(yè)大學學報,2011,43(7):1-5.
[4]FRANKLIN G F,POWELL J D,ABBAS E N.Feedback control of dynamic systems[M].4th Ed.Reading,MA:Addison-Wesley,1986:488-491.
[5]游江,趙國良,羅耀華.H∞混合靈敏度的魯棒性設計及在撓性系統(tǒng)控制中的應用[J].儀器儀表學報,2006,27(6):556-558.
[6]何朕,孟范偉,劉偉,等.H∞回路成形設計的魯棒性[J].自動化學報,2010,36(6):890-893.
[7]何朕,孟范偉,王廣雄,等.撓性系統(tǒng)的控制設計和H∞回路成形法[J].控制理論與應用,2013,30(8):1014-1020.
[8]MCFARLANE D C,GLOVER K.A loop shaping design procedure usingH∞synthesis[J].IEEE Transactions on Automatic Control,1992,37(6):759-769.
[9]KAITWANIDVILAI S,OLRANTHICHACHAT P.Robust loop shaping—fuzzy gain scheduling control of a servopneumatic system using particle swarm optimization approach[J].Mechatronics,2011,21(1):11-21.
[10]OSINUGA M,PATRA S,LANZON A.Smooth weight optimization inH∞loop-shaping design[J].Systems&Control Letters,2010,59(11):663-670.
[11]WEDELL E,CHUANG C H,WIE B.Stability robustness margin computation for structured realparameter perturbations[J].Journal of Guidance,Control,and Dynamics,1991,14(3):607-614.
[12]WIE B,LU J B.Two real critical constraints for real parametermargin computation [J]. Journal of Guidance,Control,and Guidance,1994,17(3):561-569.
[13]GRYAZINA E N,POLYAK B T.Stability regions in the parameterspace:D-decomposition revisited[J].Automatica,2006,42(1):13-26.
[14]GRYZINA E N.The D-decomposition theory[J].Automation and Remote Control,2004,65(12):1872-1884.