何朕， 姜曉明， 孟范偉， 周荻

（哈爾濱工業(yè)大學(xué)控制科學(xué)與工程系，黑龍江哈爾濱 150001）

0 引言

對(duì)于結(jié)構(gòu)不確定性，需要用結(jié)構(gòu)奇異值μ來描述系統(tǒng)的性能，而μ綜合法［1－2］是處理參數(shù)不確定性系統(tǒng)魯棒設(shè)計(jì)的一種有效的手段［3－4］。μ綜合中控制器和μ值的求取需要采用迭代的方法來進(jìn)行。雖然D-K迭代在Matlab中有現(xiàn)成的命令可供使用［5］，但是這些命令的使用需要設(shè)計(jì)者與之互動(dòng)，例如初始值的設(shè)置是否合理，迭代過程是否正常，是否還要進(jìn)行下去等。為此需要對(duì)這些命令在迭代過程中所產(chǎn)生的數(shù)據(jù)有充分的了解。可是一般的手冊(cè)或是非專業(yè)性的Matlab介紹都做不到這一點(diǎn)。本文從D-K迭代的理論出發(fā)，詳細(xì)分析了迭代過程中所產(chǎn)生的數(shù)據(jù)，使之排除設(shè)計(jì)過程中的一些似是而非的結(jié)論，引導(dǎo)設(shè)計(jì)者正確進(jìn)行D-K迭代來得出一個(gè)滿意的魯棒設(shè)計(jì)結(jié)果。本文將結(jié)合一個(gè)具體的μ綜合例題來進(jìn)行介紹，使所討論的問題更為具體，以便于了解和掌握。

1 結(jié)構(gòu)奇異值和μ綜合

本文主要討論μ綜合的算法問題。為節(jié)省篇幅，一些基本概念和定義將隨例子作些必要的說明而不單獨(dú)進(jìn)行討論。這里的例子是一個(gè)衛(wèi)星姿態(tài)控制系統(tǒng)，要求在參數(shù)攝動(dòng)10%下保證系統(tǒng)的性能，即魯棒性能問題。

設(shè)系統(tǒng)的運(yùn)動(dòng)方程［6］

式中:J1代表衛(wèi)星本體，J2為衛(wèi)星上的光學(xué)系統(tǒng)，二者是由連桿系統(tǒng)相連，J1=1，J2=0.1;k和b的名義值為k0=0.091，b0=0.000 36。本例取自文獻(xiàn)［6］，不過為了能突出本文的問題，b0取了一個(gè)更小的值。Tc為控制力矩，設(shè)輸出量為 θ2，若選擇 x=［θ2θ1］T作為狀態(tài)變量，并設(shè) Tc=u，則由式（1）和式（2）可寫得狀態(tài)空間表達(dá)式為

在下面的計(jì)算中按統(tǒng)一編號(hào)，輸出y就是圖1中的y2，而圖1中的u2就是式（1）中的u。設(shè)連桿的結(jié)構(gòu)參數(shù)隨溫度而有變化，即

式中k0為名義值;δk是界為1的不確定性，|δk|≤1;wk為相對(duì)變化范圍，設(shè)wk=10%。

本例中b的值與k有關(guān)，按泰勒級(jí)數(shù)展開得

式中b0為名義值。

將上述的參數(shù)攝動(dòng)代入式（3）可得

這里采用統(tǒng)一的符號(hào)Δ來表示不確定性，Δ=δk。式中C0x表示作用到不確定性Δ的信號(hào)，z3=C0x（見圖1）。Δ的輸出u3通過B0與狀態(tài)陣A相加。這樣，當(dāng)將不確定性單獨(dú)列出時(shí)，系統(tǒng)的狀態(tài)方程式就可以用矩陣形式表示為

而不確定性的方程式為

式（8）表明，當(dāng)考慮參數(shù)攝動(dòng)時(shí)，衛(wèi)星對(duì)象的輸入將由u2和u3兩部分組成，輸出也可分為y2和z3兩部分，而表示參數(shù)攝動(dòng)的不確定性Δ則和對(duì)象形成一種上線性分式變換的結(jié)構(gòu)，如圖1所示。

圖1中的W1、W2和W3均為常規(guī)的H∞設(shè)計(jì)中的權(quán)函數(shù)，為

式中的ρ應(yīng)在μ綜合設(shè)計(jì)中取最大值，舉例中的數(shù)據(jù)采用的是優(yōu)化設(shè)計(jì)最后一步的數(shù)據(jù)，此時(shí)ρ=0.11。

為了說明魯棒性能問題，先將對(duì)象P與權(quán)函數(shù)合寫成廣義對(duì)象G，將圖1整理成圖2的形式。圖2所示，即魯棒性能問題的框圖，要在存在參數(shù)攝動(dòng)Δ時(shí)，要求從輸入w到輸出z的H∞范數(shù)小于1，即

圖1 系統(tǒng)的框圖Fig．1 System block diagram

圖2 魯棒性能問題Fig．2 Robust performance problem

注意到這個(gè)魯棒性能要求也相當(dāng)于用一個(gè)攝動(dòng)塊ΔP跨接在w、z端上的穩(wěn)定性要求，如圖中虛線所示，‖ΔP‖∞≤1。設(shè)用M來表示控制器K回路閉合后的對(duì)象，M=Fl（G，K），則這個(gè)魯棒性能要求就完全等價(jià)于M在塊對(duì)角結(jié)構(gòu)的攝動(dòng)Δ和ΔP下的魯棒穩(wěn)定性問題，如圖3所示，或者寫成一個(gè)更為緊湊的形式，即

圖3 結(jié)構(gòu)不確定性問題Fig．3 Structured uncertainty problem

這種塊對(duì)角結(jié)構(gòu)的不確定性稱為結(jié)構(gòu)（化）不確定性。如果~Δ不是塊對(duì)角結(jié)構(gòu)，則根據(jù)小增益定理，圖3系統(tǒng)穩(wěn)定的條件就是ˉσ［M］≤1或‖M‖∞≤1。但是在現(xiàn)在的結(jié)構(gòu)不確定性下，這個(gè)奇異值（或H∞范數(shù)）條件就有保守性了，為此需要定義一個(gè)新的量，稱為結(jié)構(gòu)奇異值μΔ（M）。這里下角標(biāo)Δ表明結(jié)構(gòu)奇異值與具體的不確定性的結(jié)構(gòu)形式有關(guān)，不過一般書寫時(shí)也常略去這個(gè)下角標(biāo)。與常規(guī)的H∞設(shè)計(jì)中要求系統(tǒng)的最大奇異值的最小值小于1的概念類似，圖3系統(tǒng)穩(wěn)定的充要條件為

根據(jù)式（12）的要求來設(shè)計(jì)控制器K的方法稱為μ綜合。

結(jié)構(gòu)奇異值μ一般是通過上界來進(jìn)行計(jì)算的。而且是利用一個(gè)標(biāo)定陣D來獲得一個(gè)更緊密的上界［2］，即

這樣，當(dāng)用μ綜合法時(shí)就是要求解下列的優(yōu)化問題

注意到這里的D也是一個(gè)頻率函數(shù)D（ω），所以這個(gè)優(yōu)化問題是一個(gè)二參數(shù)的極小化問題

式（15）表明，當(dāng)D（ω）確定時(shí)，這是一個(gè)H∞優(yōu)化設(shè)計(jì)問題，可求得K。當(dāng)K確定時(shí)，每一個(gè)ω下求D，使 ˉσ［DMD－1］最小，D 代入后再重復(fù)第一步，求解H∞優(yōu)化設(shè)計(jì)問題。這樣交替進(jìn)行的求解過程稱D-K迭代。

2 D-K迭代的計(jì)算

進(jìn)行D-K迭代時(shí)，Matlab的 μ-Tools中有2個(gè)D-K迭代命令可供使用:

1）利用對(duì)話框命令dkitgui;

2）利用腳本文件dkit。這個(gè)命令實(shí)際上是調(diào)用一個(gè)已編好的程序，只要使用者對(duì)變量進(jìn)行初始化即可。

這兩個(gè)命令實(shí)質(zhì)上是一樣的，利用對(duì)話框命令時(shí)，設(shè)計(jì)時(shí)的互動(dòng)性更強(qiáng)一些，設(shè)計(jì)者可一步一步地了解迭代的進(jìn)程和選擇必要的參數(shù)。如果不需要更多的互動(dòng)，則可以選用dkit，可以較快地得出結(jié)果，甚至可以選用自動(dòng)迭代的方式直接給出結(jié)果。

下面主要是結(jié)合dkit的結(jié)果來進(jìn)行說明。在進(jìn)行D-K迭代前首先需要建立廣義對(duì)象的系統(tǒng)陣。第一步需要將圖1中的各部分用pck或nd2sys命令轉(zhuǎn)換成系統(tǒng)陣的形式，因?yàn)樵贖∞操作中是不支持傳遞函數(shù)的，都采用由狀態(tài)空間表達(dá)式所構(gòu)成的系統(tǒng)陣的形式。接下來將這些環(huán)節(jié)通過sysic命令連接成廣義對(duì)象，具體的操作步驟可以參見文獻(xiàn)［5］。廣義對(duì)象如圖4所示。需要強(qiáng)調(diào)的是，利用sysic建立廣義對(duì)象G時(shí)的輸入輸出變量順序都應(yīng)該按圖2所示的從上到下的順序進(jìn)行排列，否則可能得到的是另外一個(gè)廣義對(duì)象。

圖4 本例中的廣義對(duì)象Fig．4 Interconnection of the generalized plant

要將這個(gè)廣義對(duì)象接入到dkit的M文件中，還需要定義一組用戶變量，包括廣義對(duì)象的名稱G、對(duì)象測(cè)量輸出y的維數(shù)、控制輸入u的維數(shù)、攝動(dòng)塊的結(jié)構(gòu)參數(shù)以及頻率范圍。在定義了用戶變量后，若將這個(gè)文件保存為并命名為satellite－dk．m。在命令窗口輸入下述命令，即

這樣，dkit中就有了用戶所定義的變量，就可以運(yùn)行了。

作為例子，在進(jìn)行相應(yīng)的變量定義后，通過dkit進(jìn)行D－K迭代，第一次的迭代過程如表1所示，得γ值為65.625。

表1 第一步γ迭代的結(jié)果Table 1 First γ-iteration results

表1中λ1、λ3分別代表Hamilton陣H∞和J∞特征值的最小實(shí)部，而λ2、λ4分別表示Riccati方程解X∞和Y∞特征值的最小實(shí)部，ρ表示 X∞Y∞的譜半徑。這第一步中的標(biāo)定陣D=I，在綜合控制器K（s）的γ迭代中求解的就是H∞優(yōu)化問題

這些離散數(shù)據(jù)ˉσ［D（ωi）M（G，K）D－1（ωi）］中的最大值便是這第一次迭代所得出的μ值，本例中為2.987（見表2）。

求得K后還要進(jìn)行μ分析，即將K（s）代入系統(tǒng)中，并在各頻率點(diǎn)上求極小，為

表2 迭代結(jié)果Table 2 Iteration results

第二步D-K迭代時(shí)，利用第一步所得的K（s），并對(duì)D（ωi）的數(shù)據(jù)用有理函數(shù)陣擬合得D（s），再依次求解式（16）、式（17）的優(yōu)化問題，表2所示就是本例的各次迭代結(jié)果。第3次迭代后μ的最大值為0.997，就可以停止迭代了。圖5所示即為3次迭代后所得的最終的控制器K奇異值bode圖。

圖5 三次迭代后的控制器奇異值bode圖Fig．5 Bode plot of controller after 3 iterations

3 算法分析

D-K迭代實(shí)際上是一種交互式的設(shè)計(jì)過程，中間會(huì)交替出現(xiàn)不同的γ值和μ值，究竟哪些數(shù)據(jù)會(huì)影響到對(duì)設(shè)計(jì)過程的判斷還需要分別從K和D的求解算法上來進(jìn)行分析。

3.1 二分法γ迭代和γ的上界

D－K迭代中當(dāng)D（ω）確定下來后就是一個(gè)H∞優(yōu)化設(shè)計(jì)問題，見式（15）。H∞優(yōu)化問題是指求解最小的H∞范數(shù)值γ，并給出H∞控制器K（s）。這個(gè)最小的γ值一般采用二分法來進(jìn)行。即先給定一個(gè)γ值的范圍，例如0～100，然后判斷系統(tǒng)的H∞范數(shù)是否小于這個(gè)給定的上限值100。如果滿足γ＜100的條件，再將這個(gè)上限折半，看是否小于50。如果不滿足，則在50～100之間找。這樣折半尋找，稱為γ迭代，一直到滿足誤差要求的最小值。這里的H∞范數(shù)是否小于某個(gè)給定的γ值所依據(jù)的是H∞理論中的2個(gè)Riccati方程法，即DGKF法［7］。DGKF法指出，H∞問題是否有解（即H∞范數(shù)＜ γ）的條件是Hamilton陣H∞和J∞是否屬于dom（Ric），且對(duì)應(yīng)的解X∞和Y∞是否是半正定和ρ（X∞Y∞）是否小于1，式中ρ為譜半徑。因此在γ迭代過程中每次都要計(jì)算H∞陣和J∞陣的特征值和廣義特征向量。用計(jì)算機(jī)求解矩陣的特征值時(shí)要考慮到算法的數(shù)值穩(wěn)定性問題［8］。為此在求解特征值時(shí)一般是將矩陣變換成Schur上三角形，Schur上三角形陣的對(duì)角元就是特征值。Schur變換的算法不是一種有限步的算法，而是一種迭代的QR算法［8］。QR算法是一個(gè)公認(rèn)的算法，一般都可很快收斂到真實(shí)值。但若遇到病態(tài)的特征值，則會(huì)出現(xiàn)誤差［8］。所謂病態(tài)特征值是指2個(gè)特征值非常接近的情況（例如有重根的情況）。在二分法γ迭代過程中，不斷減小的γ值代入到高階的H∞和J∞陣中難免會(huì)出現(xiàn)重根的情況，這時(shí)就可能會(huì)得出錯(cuò)誤的結(jié)果。具體來說，本例中dkit是按各默認(rèn)值來運(yùn)行的，γmax=100，第一次的判別是γ＜100（見表1），折半后要檢查是否小于50，表1第二行表明這時(shí)Y∞陣的一個(gè)最小的特征值已小于0，即Y∞已非半正定，表明系統(tǒng)的H∞范數(shù)不小于50，其實(shí)這是錯(cuò)判了。經(jīng)檢查，Y∞是Hamilton陣J∞的解，而γ=50時(shí)，這個(gè)J∞陣恰好存在兩對(duì)重根:－1 000、－1 000和1 000、1 000。這2對(duì)病態(tài)特征值導(dǎo)致γ=50時(shí)迭代失敗。這個(gè)判斷失敗，導(dǎo)致需要到50～100之間的二分法尋找，最終得這次γ迭代的結(jié)果是65.625。由此可見H∞優(yōu)化中先后會(huì)出現(xiàn)2個(gè)γ值，一個(gè)是根據(jù)是否有解的條件得到的65.625，另一個(gè)則是加上控制器后實(shí)際得到的H∞范數(shù)。這前一個(gè)值可以稱為是γ的上界。所以當(dāng)采用dkit進(jìn)行D-K迭代時(shí)，首先顯示的就是γ，如果見到這么大的γ值（65.625），就可能放棄這次數(shù)據(jù)而令外尋找別的方案。其實(shí)這只是一個(gè)假象，真正的γ值并沒有這么大，就本例而言，經(jīng)過驗(yàn)證其真實(shí)值約等于第一步迭代產(chǎn)生的μ值，見表2的設(shè)計(jì)結(jié)果。

3.2 D擬合的收斂問題

D-K迭代中當(dāng)上一步求得K（s）后，下一步就是要求解D（ω）使 ˉσ［D（ω）Fl（G，K）D－1（ω）］達(dá)到最小。這一般是在上一步先得出各離散頻率點(diǎn)上達(dá)到極小值時(shí)的D（ωi），然后再用有理函數(shù)陣D（ω）來進(jìn)行擬合，將擬合所得到的D（ω）代入式（16），再求解H∞優(yōu)化問題。所以這里也有兩套數(shù)據(jù)的問題，一個(gè)是逐點(diǎn)尋優(yōu)所得的范數(shù)，另一個(gè)則是用有理函數(shù)代替后求得的H∞優(yōu)化解。如果D（ω）擬合得不好，則這兩套數(shù)據(jù)就會(huì)有矛盾。一般來說，D-K迭代的過程總是收斂的，但如果迭代過程中發(fā)現(xiàn)μ值收斂不好（有波動(dòng)或變大），則很可能就是標(biāo)定陣D（ω）的擬合不好。一個(gè)可能的原因是求極小值時(shí)離散點(diǎn)不夠密，沒有完全覆蓋系統(tǒng)中的高頻分量，例如弱阻尼的諧振模態(tài)。

4 結(jié)論

D-K迭代中包含γ迭代和標(biāo)定陣D的求極小，在依次迭代中會(huì)出現(xiàn)多個(gè)范數(shù)值，弄清這些數(shù)據(jù)產(chǎn)生的背景，對(duì)正確進(jìn)行D-K迭代是至關(guān)重要的。其實(shí)是對(duì)于γ迭代來說，因?yàn)镠amilton陣的特征值隨每次的γ值而有變化，當(dāng)出現(xiàn)病態(tài)特征值時(shí)，γ迭代過程就會(huì)停止在一個(gè)較高的虛假值上，加上控制器后的系統(tǒng)的H∞范數(shù)才是真正的γ值。而且這里的討論也具有普遍意義，因?yàn)镠∞控制理論中有多種方法，包括線性矩陣不等式法，這些算法得出的并不是真實(shí)的γ值，而只是一個(gè)γ的上界，加上控制器后實(shí)際得到的γ值才是真正的H∞范數(shù)。

［1］ 周克敏，DOYLE J C，GLOVER K．魯棒與最優(yōu)控制［M］．毛劍琴，鐘宜生，林巖，等譯．北京:國(guó)防工業(yè)出版社，2002:306－338．

［2］STEIN G，DOYLE J C．Beyond singular values and loop shapes［J］．Journal of Guidance，Control and Dynamics，1991，14（1）:5－16．

［3］LANZON A，TSIOTRAS P．A combined application of H∞loop shaping and μ-synthesis to control high-speed flywheels ［J］．IEEE Transactions on Control Systems Technology，2005，13（5）:766－777．

［4］YIN Guodong，CHEN Nan，LI Pu．Improving handling stability performance of four-wheel steering vehicle via μ-synthesis robust control［J］．IEEE Transactions on Vehicular Technology，2007，56（5）:2432－2439．

［5］BALAS G J，DOYLE J C，GLOVER K，et al．μ-Analysis and Synthesis Toolbox［M］．Natic:MathWorks and MUSYN，1993．

［6］ FRANKLIN G F，POWELL J D，EMAMI-NAEINI A．動(dòng)態(tài)系統(tǒng)的反饋控制［M］．朱齊丹，張麗珂，原新，等譯．4版．北京:電子工業(yè)出版社，2004:488－498．

［7］DOYLE J C，GLOVER K，KHARGONEKAR P P，et al．State space solutions to standard H2and H∞control problems［J］．IEEE Transactions on Automatic Control，1989，34（8）:831－847．

［8］MACIEJOWSKI J M．Multivariable Feedback Design［M］．New York:Addison-Wesley Publishing Company，1989．

（編輯:張?jiān)婇w）