彭長樂，陳 城，侯和濤

(山東大學土建與水利學院，山東，濟南 250000)

結構控制的概念在土木工程領域提出以來，經過多年來的研究發(fā)展，在減小結構的動力荷載響應方面顯示出了巨大的潛力[1-2]，其中磁流變阻尼器因為其低能耗、出力大、響應快、性能可靠等優(yōu)點備受研究者們關注[3-4]。特別是自Lord公司研制出最大出力200 kN可用于土木工程領域的足尺阻尼器以來，磁流變阻尼器更成為結構控制領域的熱點。要使磁流變阻尼器的控制性能得到充分的發(fā)揮，獲取一個準確的數值模型是十分必要的。迄今為止，許多磁流變阻尼器的模型相繼被提出來，其中包括Bouc-Wen 模型[5]、Viscous-Dahl模型[6]，Hyperbolic Tangent模型[7]和Maxwell Nonlinear Slider (MNS)模型[8]等。這些模型的參數多以進化算法設計目標函數，基于一系列的正弦位移識別試驗，通過優(yōu)化計算出最優(yōu)參數[9]。

然而由于在建模情況中所做的簡化假設，以及有限的觀測數據，傳統(tǒng)的確定性參數現(xiàn)象學模型很難考慮由于環(huán)境因素、模型誤差以及測量誤差帶來的不確定性，從而導致數值模型的出力結果與實驗結果存在不同程度的差異。這些差異可能使得數值模擬中的結構在地震響應下響應結果與實際實驗不符合，也會一定程度影響磁流變阻尼器的控制算法效果。Caicedo等[10]證明使用模型參數的概率表征方式有助于更真實的模擬預測。

MNS模型利用Hershel-Bulkley粘塑性來描述阻尼器中發(fā)生的MR流體的后屈服非牛頓流體行為，即剪切稀化和增稠行為[8]。MNS模型能比較準確描述阻尼器的位移與出力、速度與出力的非線性滯回行為。但由于MNS模型本身屬于現(xiàn)象學模型的一種，阻尼器的復雜非線性行為使得建模仍存在偏差，因此單一的確定性模型并不能完全滿足預測的要求。而不確定性參數能幫助克服參數優(yōu)化時出現(xiàn)“異參同效”的影響，最為有效地利用阻尼器模型的預測能力。

本文運用馬爾可夫蒙特卡洛方法對磁流變阻尼器的MNS模型進行了不確定性分析。通過選擇合適的似然函數生成模型參數的后驗分布，并對參數分布形式，參數之間的內在相關性進行分析。在此基礎上進一步與確定性參數在正弦位移曲線實驗和實時混合模擬實驗下的出力誤差、耗能預測方面進行對比總結，分析不確定性模型的特點。

1 MNS確定性模型介紹

磁流變阻尼器的MNS模型由屈服前和屈服后兩部分組成[9], 如圖1所示。

1.1 屈服前模式

屈服前的阻尼器行為由Maxwell單元描述，由一個系數為c的阻尼元件與一個剛度為k的彈簧元件串聯(lián)而成。其阻尼力f為：

圖1 MNS模型示意圖Fig.1 Schematic diagram of MNS model

式中：y表示彈簧元件與阻尼元件的總位移；z表示阻尼元件的位移。當阻尼器處于屈服前模式，阻尼元件與彈簧元件的總速度與實際阻尼器速度相等。在阻尼器初始狀態(tài)時y與x相等，直到阻尼器出現(xiàn)第1次屈服。在屈服后模式下，y與z需要通過式(1)的變式，并利用式(2)計算的屈服后出力f值持續(xù)更新。

當最終模式從屈服后返回屈服前，這時的y=x+σ，并且在此后的屈服前計算中σ恒定不變。

1.2 屈服后模式

當阻尼器行為模式從預屈服變?yōu)榍?，阻尼力的斜率發(fā)生改變，阻尼力通常隨著速度的增加而增加。在屈服后模式下，MNS模型通過運用非線性滑塊(nonlinear slider)來描述阻尼器出力與速度之間的關系，將正負向出力分開兩組表示，從而有利于提高阻尼器不對稱下的預測。正向阻尼器出力的數學表達式如下：

第1部分基于Herschel-Bulkley粘塑性理論，包括、a、b、n這4個參數。第2部分是與Herschel-Bulkley粘塑性理論定義的切線。負向阻尼器出力可以類似公式表達，因此同樣也包含4個負向出力參數。從加載、卸載路徑上來看，磁流變阻尼器的出力也存在差異，可以利用增加式(6)的方式模擬不同階段下的出力。

1.3 模式改變判別

當阻尼器出力從Maxwell單元達到屈服后曲線，模型從屈服前轉換屈服后。數學上這種狀態(tài)可以表達為：

當模型從屈服后回到屈服前，只需要速度滿足:式(8)就會出現(xiàn)：

式中，y可以利用式(1)的變形計算：

2 試驗數據

本文使用的試驗數據基于科羅拉多大學智能結構技術實驗室完成的足尺磁流變阻尼器性能測試試驗[10]。由于篇幅的限制，本文主要考慮阻尼器在最大電流2.5A下，即為Passive-On模式下的不確定性研究。表1列出了本文選取的6組正弦位移加載實驗的頻率幅值，每組都在2.5 A的電流輸入下進行10次往復加載。

表1 正弦試驗參數Table 1 Sinusoidal test parameters

通過利用粒子群算法[9]，設計最小均方根誤差為目標函數，經過多次迭代，計算出模型與試驗相對誤差最小的確定性MNS模型參數值列于表2中。

表2 MNS模型參數信息統(tǒng)計結果Table 2 Statistics of MNS model parameter

3 馬爾可夫鏈蒙特卡洛方法

馬爾科夫鏈蒙特卡洛算法可以獲得參數后驗分布的一系列樣本，適用于非標準分布不獨立的多變量分布模擬[11]。其基本步驟如下：

1) 建立先驗分布：由于對于本文模型沒有任何先驗信息。因此本文的所有參數的先驗分布Prior都假設為在區(qū)間內的均勻分布。

2) 初始點：合適的參數初始值可以減少達到收斂所需的計算時間，提高運算效率。本文中選擇已有文獻中的確定性參數作為初始點。

3) 計算似然值：每組參數值會生成對應的阻尼力模型輸出，通過計算與試驗輸出的相對誤差，得到模型擬合程度的函數如下：

假設似然函數符合均值為0，標準差為ε的高斯分布。ε的選值直接影響后驗參數的范圍，過小的ε會使參數收斂到一個狹窄的區(qū)域，不確定性太小。過大的ε則會使得參數范圍太大，導致不確定分析的結果過度偏離試驗結果。

4) 計算后驗概率：

5) 生成新樣本：加上隨機增量會新生成一組參數，新參數計算的后驗如果滿足式(13)則認為可以接受：

其中，γ服從[0，1]的均勻分布。重復步驟3)~步驟5)，生成長度為N的馬爾科夫鏈。選擇收斂后的鏈長M為參數樣本進行不確定性分析。

4 參數分析結果

通過選取ε值為0.5%，設定抽樣長度為十萬，得到了相應的后驗分布參數樣本。本文利用Gelman-Rubin原則[12]，通過分析多個馬爾可夫鏈之間的差異來評估MCMC收斂。通過比較每個模型參數的估計的鏈間和鏈內方差，確定MCMC的結果在50000次后基本達到收斂。最終選擇收斂后的馬爾科夫鏈作為合理近似參數樣本。

圖2做出不確定參數分布與粒子群算法得到最優(yōu)確定性參數的對比?？梢园l(fā)現(xiàn)所有參數都為單峰分布。其中c與k呈現(xiàn)正態(tài)分布，其他參數則表現(xiàn)出對數正態(tài)分布的趨勢。從圖2中還可以發(fā)現(xiàn)MNS模型的確定性參數值都處在不確定性參數的分布之中，并且大部分都位于概率密度函數的峰值附近，卻不一定恰好出現(xiàn)在頂峰位置。這主要是因為參數之間出現(xiàn)一些相關性，導致最優(yōu)的參數組合不一定都處于各自參數的最可能出現(xiàn)的點，這也從側面反映了參數之間的相互影響會導致確定性參數的取值，單組參數并不能完全發(fā)揮出模型的預測能力。表2列出了對各參數的統(tǒng)計信息進行分析得到的前3階矩。

圖2 不確定性參數分布與確定性參數對比Fig.2 Comparison of probabilistic parameter distribution and deterministic parameters

圖3 MNS模型參數變異系數Fig.3 Coefficient of variation for MNS model parameters

圖3使用變異系數來具體定量分析參數本身的不確定性大小。從圖3中可以觀測到，在由馬爾可夫鏈蒙特卡洛方法得到的不確定性參數中，不同參數分布的變異系數不盡相同。對于有明顯物理意義的參數，描述阻尼器在未屈服狀態(tài)下的阻尼系數c與剛度系數k的變異系數明顯要小于其他參數。而用于擬合阻尼器在屈服后Hershel-Bulkley曲線的模型參數，則顯示出較強的不確定性。這說明對于有較強客觀物理意義的參數，其識別能力更強，不確定性更小。對于現(xiàn)象學模型中的參數，由于其本身沒有明顯物理意義，一般為擬合曲線而設計，所存在的參數組合可能性更多，不確定性更大。

為考慮參數相關性的影響，本文對參數組合間的相關系數矩陣進行分析計算，參數之間的某些較強相關性說明參數相互間是多余的，可以為模型簡化提供指導。圖4可視化出參數之間的相關性大小，從中明顯可以看出：描述屈服前的參數阻尼系數c與剛度系數k存在較強的相關性，為0.5左右；描述屈服后的模式a，b，n之間同樣存在很強相關性，大致為0.6~0.9，并且僅僅出現(xiàn)在相同正負的方向中；其余的參數之間相關性不強，幾乎在-0.2到0.2之間，可以基本忽略相關性而認為這些參數是相互獨立的。

圖4 MNS模型參數相關系數Fig.4 Correlation between MNS model parameters

5 基于不確定性模型的試驗預測

模型參數不確定性分析的最終目的是利用模型參數對真實物理模型的響應進行有效的預測。本節(jié)分別從試驗的阻尼力輸出和耗能預測兩個指標，對不確定性模型和確定性模型進行對比評估。以表1中25.4 mm，1 Hz正弦位移加載試驗為例，圖5給出了磁流變阻尼器阻尼力輸出時程比較，其中實線為不確定性參數模擬的95%預測區(qū)間以及預測均值，虛線為確定性模型預測曲線，點劃線為試驗值?？梢钥吹酱_定性參數預測曲線幾乎與不確定性曲線的均值預測重合。另外試驗觀測值大部分都處于95%的預測區(qū)間之內，表明不確定性模型具有良好的預測能力。值得指出，試驗值仍有處于預測之外的值，這主要受限制于模型本身的精確程度，參數的不確定性分析，并不能完全補償模型誤差的影響。

圖5 磁流變阻尼器輸出力比較Fig.5 Comparison of force outputs

為定量分析，本文使用歸一化均方根(Root Mean Square, RMS)誤差作為評價模型參數優(yōu)劣的目標函數。

式中：ERMS越小，模型預測與實際輸入的差距越??；當ERMS=0時，表示模型與試驗結果完全一致。對于表1中的六組模型識別試驗，圖6比較了MNS 模型對應確定性和不確定性參數的ERMS。

圖6 模型誤差分布比較Fig.6 Comparison of RMS error

從圖6可以發(fā)現(xiàn)，本文中得到的不確定性參數的誤差基本處在3%~4%之間，而確定性算法得到一組參數值的誤差接近不確定性模型的誤差最小值處。在參數識別過程中，不確定性模型適當放寬了接受參數的范圍，而非僅僅接受最優(yōu)的一組參數。使用這些參數的組合可以在合理范圍內拓寬模型的預測區(qū)間，從而提高模型的預測性能。

為了進一步檢驗對MNS模型參數不確定性分析的預測效果，本文采用了Chen等[13]進行的實時混合模擬試驗結果。該系列混合模擬試驗選用了1999年Chi-Chi地震在TCU105站點觀測到的地震運動加速度并縮放到Design Basis Earthquake (DBE)和Maximum Considered Earthquake (MCE)兩種情況。磁流變阻尼器在DBE和MCE下實時混合模擬試驗(RTHS)的最大位移分別在50 mm和70 mm左右，對應的磁流變阻尼器速度-出力曲線如圖7和圖8所示。

對應不同情況下阻尼器的出力情況，計算出的均方根誤差ERMS分布如圖9和圖10所示。對于這兩組混合模擬試驗結果，確定性參數的阻尼力輸出誤差不再處于不確定性參數模型的誤差最小點。這說明對于不同的實驗，最優(yōu)參數并不是唯一的，這也同時說明了不確定性分析的重要性。完全靠一組參數值是很難做到在實驗情況出現(xiàn)變化時仍保證較高預測精度的。而在控制不確定性大小合適的情況下，可以彌補確定性參數的不足，從而提高MNS模型對不同試驗的預測能力。

圖7 實時混合模擬試驗DBE下速度與出力的關系Fig.7 Velocity-force curve of RTHS with DBE

圖8 實時混合模擬試驗MCE下速度與出力的關系Fig.8 Velocity-force curve of RTHS with MCE

圖9 實時混合模擬試驗中震下ERMSFig.9 ERMS for real time hybrid simulation with DBE

圖10 實時混合模擬試驗大震下ERMSFig.10 ERMS of real time hybrid simulation test with MCE

從圖9和圖10進一步可以看出不確定模型在正弦位移加載方式下的誤差基本處于3%~4%之間，而在實時混合模擬試驗的位移響應下的誤差增加至7%~11%。這說明了MNS模型在有規(guī)律正弦位移加載方式下的預測能力要強于模型對于實時混合模擬試驗的預測。這提醒我們在參數優(yōu)化過程中，過分依賴使用有規(guī)律的位移加載方法會低估模型在真實結構位移中的誤差。

磁流變阻尼器作為建筑結構中的控制裝置，主要起耗散能量的作用。因此在用模型預測與實際出力的誤差來衡量參數優(yōu)劣外，本文進一步用耗能預測來檢驗模型不確定性分析的預測效果。

首先對正弦位移識別實驗本身進行耗能預測，圖11對不確定性參數預測的耗能分布與正弦位移

圖11 不同頻率、幅值下耗能預測(直方圖為不確定性模型耗能預測，虛線為確定性模型預測耗能，點劃線為試驗耗能)Fig.11 Energy dissipation prediction under different frequencies and amplitudes (The histogram is the uncertainty prediction, the dashed line is the deterministic prediction, and the dotted line is the experimental energy consumption)

實驗實際耗能預測結果進行對比。從圖中可以觀測到：使用不確定性參數模型得到的耗能分布基本上都處于近似正態(tài)分布；確定性模型預測一般都接近于不確定性模型預測的均值處。并且試驗得到的實際耗能都處于不確定性模型預測區(qū)間內部，這說明使用這組不確定性參數組合預測這6組試驗是完全可行的。

圖12對圖7和圖8中實時混合模擬試驗進行耗能預測分析。確定性參數預測仍接近不確定性模型預測均值附近，并且不確定性參數的預測耗能區(qū)間依然將試驗耗能包絡其中。但是值得注意的是，相對于正弦位移識別試驗。實時混合模擬下的試驗耗能，處于耗能預測的更邊緣位置，這說明模型的耗能預測能力下降，但是仍處于5%~95%預測范圍區(qū)間之內。這與使用ERMS分析得到的結論類似，即正弦位移下MSN模型預測比實時混合模擬下的結構響應出力要更準確。

圖12 不同地震水平下磁流變阻尼器耗能預測(直方圖為不確定性模型耗能預測，虛線為確定性模型預測耗能，點劃線為試驗耗能)Fig.12 Energy dissipation prediction of MR damper for real-time hybrid simulation under difterent seismic levers (The histogram is the uncertainty prediction, the dashed line is the deterministic prediction, and the dotted line is the experimental energy consumption)

6 結論

本文對200 kN足尺磁流變阻尼器的MNS模型在2.5 A下的參數不確定性進行了分析。通過選擇合適的不確定性大小與似然函數，基于馬爾科夫鏈蒙特卡洛方法，生成了一系列參數組合值。并使用后驗參數分別對正弦位移試驗和實時混合模擬下的阻尼器出力進行了預測評估。根據分析結果得出以下結論。

(1) 從生成參數的分布與相關性中分析發(fā)現(xiàn)：最優(yōu)的一組參數值被包含在不確定性參數的分布中，并且大部分確定性參數都趨近于均值處；某些參數間存在強烈的相關性，這為去除多余參數，簡化模型提供了思路；參數的變異系數不盡相同，有實際物理意義的較小，更容易識別，而現(xiàn)象模型參數則存在更多可能性。

(2) 利用不確定性分析的參數分布預測試驗出力可以發(fā)現(xiàn)，大部分試驗出力處于不確定性模型的5%~95%預測區(qū)間以內，但是不足以完全包絡。這意味著不確定性模型可以提高預測能力，但是并不可能完全彌補模型本身的誤差。

(3) 在對歸一化均方根誤差的分析中，選用的正弦位移識別試驗時，確定性值的誤差最小，而在混合模擬驗證試驗中，這組確定性值的出力預測卻并不處于誤差最小值點，并且在實時混合模擬驗證試驗中誤差明顯變大了。說明單一最優(yōu)參數的適應性較弱，而使用不確定性參數則可以減小由于試驗情況改變引起的誤差增大。

(4) 最后通過對正弦識別試驗和實時混合模擬驗證試驗進行耗能預測。所有選取的試驗中，確定性模型預測處于不確定性模型預測區(qū)間的均值處，而阻尼器實際耗能都處于預測區(qū)間內部，這說明不確定性模型對于阻尼器耗能的預測是可行的，可以為進一步預測建筑結構在使用磁流變阻尼器時的位移響應預測提供參考。