一類雙層薄膜結構振動能量采集器的數據驅動建模方法及應用1)

邱宏蘊 王志霞 丁 北 王 煒

(天津大學機械工程學院,天津 300350)

引言

在人們觀察世界,開展生產實踐,預測實踐結果的過程中,建模方法發揮著至關重要的作用.傳統建模方法通過對基礎理論的梳理歸納以及數學演繹,得到適用于當前問題的模型.然而,該方法極度依賴先驗知識,對于線性、直觀問題的分析較為有效,而對于包含強非線性因素的系統則適用性有限,也因此為基于實驗數據的建模方法提供了成長空間.

數據驅動建模方法,以工程實踐過程中獲取的物理實驗數據及有限元仿真結果為基礎,挖掘其中隱含的規律性,從而獲取真實反映系統運動狀態的動力學方程,現階段已逐步成為國內外學者關注的熱點問題.在相關領域,Frank 等[1]提出嶺回歸,實現了較高精度的模型預測功能.Tibshirani[2]基于最小二乘法(ordinary least square,OLS)和嶺回歸提出Lasso,該方法能夠充分壓縮模型,獲得復雜度更低、解釋性更好的模型.然而,Lasso 僅對于維度較低、數據量較大的問題具有較好的效果[3],且對特征值相差過大的數據模型存在過度壓縮現象[4],這無疑降低了其廣泛應用于實際場景的可能性.為了解決此類問題,Zou 等[5]提出彈性網絡(elasticnet,EN),該方法補充了Lasso 在高線性相關數據處理方面的不足,具有更好的模型壓縮效果.此后,諸多學者[6-8]建立了基于實驗數據的模型辨識方法,以識別動力系統的輸入、輸出和噪聲等信號在系統模型中的表現形式.上述方法主要識別非線性振動方程中的具體參數,并不能直接建立系統動力學模型,因此仍屬于參數識別的范疇.此后,Donoho[9]利用壓縮感知中關于近似最優子空間的研究,證明了Banach 空間理論在稀疏識別中應用的可行性.Wang 等[10]采用稀疏識別和壓縮感知理論構建非線性動力學系統,以實現對故障信號的預測和診斷,但是由于模型復雜度高以及頻繁出現的模型壓縮過度等問題[11],仍不能直接建立完整的動力學模型.直至2016 年,Brunton 等[12]提出了稀疏非線性動力系統辨識算法(sparse identification of nonlinear dynamical systems,SINDy),實現了從實驗數據到非線性控制方程的過渡,并且極大地簡化了模型復雜度,提高了算法魯棒性.

此后,Rudy 等[13]基于SINDy 提出了偏微分方程函數識別(partial differential equations functional identification of nonlinear dynamics,PDE-FIND),江昊等[14]、劉吉悅[15]通過對PDE-FIND 的研究,證明其對于強非線性偏微分方程的識別也具有較好精度,彌補了SINDy 難以識別非線性偏微分方程的不足,進一步提升了方法的實用性.聶滋森等[16]將稀疏正則化方法應用于懸臂梁損傷識別中,僅使用少量數據即可識別損傷的位置和強度,實現了正則化方法在小訓練集問題上的實際應用.與另一種廣泛使用的數據驅動方法動態模式分解[17](dynamic mode decomposition,DMD)相比,SINDy 無需假設基本形式,可直接通過實驗或者仿真數據獲取系統的非線性控制方程,因此實現了模型復雜性與計算精度之間的平衡,并且在機械和生物等諸多領域具備應用價值[18-19].

然而,工程實踐中的數字信號通常會受到噪聲干擾,即使微小的噪聲信號,都可能導致誤差的急劇放大,進而在一定程度上影響以SINDy 為代表的數據驅動建模算法的分析精度[20].因此,學者們將濾波模塊與數據驅動方法進行前后疊加,形成了提升建模方法可靠性的噪聲數據處理方法.如Brunton 采用的全變分正則化(total-variation regularization,TVR)[21]處理,以及Zou 等[22]提出的稀疏主成分分析(SPCA),在降維的同時也達到了降噪的目的.與前后疊加的噪聲處理方法類似,這些方法依次進行數據濾波與數據驅動,先后達到降噪與建模的目的,其降噪效果雖然顯著,卻未必能達到最佳的系統識別精度.因此,處理好數據濾波與數據驅動之間的關系,實現二者的有機結合,而并非簡單疊加,才是減小噪聲影響,提高識別精度的可能途徑,值得進一步研究.

另一方面,在數據處理過程中,當非線性函數庫構成的狀態矩陣特征值差異較大,各分量間線性相關程度較高時,矩陣將呈現病態矩陣的性質.SINDy在處理此類病態矩陣時,容易表現出模型壓縮過度[20]和多解[23]等問題,即所謂奇異性問題.面對奇異性問題時,SINDy 可能表現出較低的魯棒性,并影響其識別結果的準確性[24-25].因此,解決奇異性問題對于提升以SINDy 為代表的數據驅動方法的實用價值具有重要的現實意義.

作為應用,本文討論的電磁式薄膜振動能量采集器(electromagnetic vibration energy harvester,EMH),由環境振動作為激勵源,帶動薄膜上的永磁體振子做低頻振動,并導致相應的線圈內部的磁通量發生變化,產生感生電動勢.為深入研究該EMH的復雜振動特性,需要構建精確的動力學模型.然而,對于振幅較大的薄膜結構而言,強幾何非線性以及阻尼因素的影響較為明顯,這為直接運用力學理論建模增加了難度.在包含強非線性因素的控制方程中,線性項與非線性項系數數量級差距很大,狀態矩陣的構造方式更加復雜,不僅容易產生模型壓縮過度[20]和多解[23]等問題,對數據驅動算法的魯棒性以及建模精度也提出了更高要求.因此,嘗試數據驅動方法與非線性理論相結合的建模方式,更適合精確構建此類薄膜結構EMH 數學模型,進而開展復雜的振動行為分析.

綜上所述,針對數據驅動方法難以處理的強噪聲、強奇異性問題,本文將SINDy 算法進行相應的改進,提出了增強型稀疏非線性系統辨識算法(enhansed sparse identification of nonlinear dynamical systems,ESINDy).結合傳統非線性理論研究方法,通過研究一類雙薄膜式雙穩態振動能量采集器的動力學特性,對ESINDy 算法精度及稀疏效果進行分析和驗證,以期為數據驅動方法的實際應用提供參考.

1 識別算法

1.1 SINDy 算法及數據濾波算法

由于傳統數據驅動建模方法存在耗時長、模型復雜等問題,Brunton 等[12]提出了SINDy 方法,并在其中加入了TV-R 數據濾波算法.該方法是一種利用實驗數據進行數據驅動建模的方法,主要原理是將給定的一系列非線性函數組成基函數庫,并將基函數庫構建的狀態矩陣與TV-R 算法構造的導數矩陣形成線性殘差優化問題,求解該優化問題得到各個基函數分量的系數,最終獲得較為精確的控制方程.

(1)數據濾波算法

總變分正則化(TV-R)是最常用的數據濾波手段之一,但求解較為復雜.Gong[26]提出的曲率濾波(curvature filters,CF)在TV-R 的基礎上,使用更簡結的算法達到更好的濾波效果.

傳統正則化方法的正則化項能量會隨著優化進程的進行而逐漸降低,而正則化項作為優化主體,能量的降低會導致優化效果變差,曲率濾波使用已知的曲率來優化其對應的正則項,從而達到更好的優化效果.對一維含噪信號x(t),曲率可表示為

其中,Ix和Ixx分別表示x(t)對t的1 階、2 階導數.對一維含噪信號,曲率表示的正則化項可表示為

其中,λ為正則化參數.在dt的步長下,曲率優化更新過程可表示為

正則化方法依賴正則化參數的選擇[27],然而在實際問題研究過程中,通常需要先進行數據濾波,隨后使用濾波后的數據進行SINDy 數據驅動.這就導致正則化參數不能隨著數據驅動問題的變化而自動選取,限制了SINDy 的工程場景應用[28].

(2)SINDy 算法

假設振動系統的方程可以表示為如下形式

其中,x(t)為系統在t時刻的狀態,f(x(t))為系統控制方程.矩陣 Θ 為由給定非線性函數構成的基函數庫

根據基函數庫,構造狀態矩陣A

假定控制方程f(x(t))可由基函數庫中的部分項構成,Ξ 為基函數庫中各成分系數構成的矩陣,則振動系統方程可表示為

因此,在已知導數矩陣與狀態矩陣A的前提下,對振動系統方程的識別問題可轉化為對的識別問題.該問題可通過最小二乘法(OLS)解決,目標函數Target可以表示為均方誤差形式

化簡式(9),最小二乘法解可表示為

Lasso 法[2]通過在最小二乘法目標函數后增加L1 正則化項的方法,達到稀疏效果.稀疏方法可以簡化復雜的數學模型,得到項數較少的簡單模型來表征系統的動力學特性,并綜合考慮模型壓縮程度和模型精確性之間的矛盾.Lasso 法的求解目標可表示為

求解駐點,并將式(10)代入

SINDy 以Lasso 解作為循環主體函數,將式(12)不斷迭代,使得均值誤差處在一定容許范圍內,從而求得稀疏矩陣的局部最優解,進而完成控制方程的重建.

雖然SINDy 具有不依賴先驗知識,能夠直接從數據中建立控制方程的算法優勢,但由于Lasso 方法本身對特征值差距較大的系數識別效果較差[4],這就使得SINDy 面對病態矩陣時也產生了類似的諸如模型壓縮過度等奇異性問題[20,23],需要進一步改進.

1.2 ESINDy 算法

具體而言,關于SINDy 存在受噪聲數據影響顯著和難以識別病態矩陣等不足,主要表現為: 一方面,SINDy 雖包括基礎的抗噪功能,但面對強噪聲就需要引入額外的降噪手段,這其中數據處理和識別過程會放大噪聲,需要加入額外的數據濾波模塊減小噪聲影響.現有數據驅動方法針對噪聲的處理通常是在數據驅動前對數據信號進行降噪預處理.這種方法雖然簡單便捷,但數據濾波與識別作為相對獨立的兩個過程被割裂開來,因此從算法的整體性和完整性角度考慮,還存在提升空間[28].另一方面,SINDy 通過構造時間序列矩陣,從而將控制方程的識別問題轉化為線性優化問題,并使用Lasso 法作為循環主體函數對優化問題進行分析.時間序列矩陣構造方法的局限性會加大其內部低階成分與高階成分間線性相關程度和特征值差異,導致以Lasso作為循環函數體的SINDy 在處理相關數據時存在明顯的局限性.

有鑒于此,本節從數據濾波模塊和稀疏識別循環函數體兩個方向入手,提出改進型SINDy 分析方法(ESINDy),如圖1 所示.該方法將曲率濾波正則化參數與識別結果的稀疏性(sparsity,SP)進行關聯,實現數據濾波與數據驅動的有機結合,以改善數據濾波效果;同時,使用彈性網絡(EN)代替Lasso 作為循環主體函數,使其在面對線性較強、特征值差異較大的數據時有更好的處理效果;最后,借助帕累托前沿分析(Pareto front analysis)[29]確定最小均值誤差,從而自動篩選最優解.

圖1 SINDy 分析模式與ESINDy 分析模式對比Fig.1 Comparison between SINDy and ESINDy

(1)數據驅動部分

使用L1 正則化項的Lasso 能夠實現良好的稀疏效果,但對病態數據的識別能力存在顯著不足,導致SINDy 在處理此類數據時出現精度降低、壓縮過度等問題.使用L2 正則化項的嶺回歸,雖然無法得到精簡模型,卻能在處理上述數據時具有更好的魯棒能力,實現較高精度的識別.EN 同時使用L1 和L2 正則化項,巧妙地結合了L1 的稀疏能力與L2 的魯棒能力,因而具有更好的識別效果.EN 的等價優化問題為

經推導可知閉式解可表示為

(2)數據濾波部分

TV-R 和曲率濾波等正則化濾波算法具有處理頻帶寬、降噪效果好等優點,然而正則化參數通常使用經驗選擇,極大地限制了其通用性[27].在所有可能解中,稀疏性意味著能夠使用較少的項描述模型的本質,從而制定較為精簡快速的控制方案.因此,篩選適合的正則化參數,使得識別結果稀疏度最大,并使用該參數執行曲率濾波,即可識別到更為稀疏的模型.

ESINDy 算法流程如表1 所示,其中:X為待處理數據,λ1,λ2為正則化參數,iter為循環次數,Sp為通過L1 與L2 范數間差異定義的稀疏度,表達式為

表1 ESINDy 算法流程Table 1 ESINDy algorithm

上述方法不僅解決了正則化參數篩選困難的問題,還將曲率濾波與SINDy 有機結合,賦予了正則化參數新的物理意義,從而提高數據濾波效果,間接提高了稀疏識別的準確性.接下來將針對R?sler 系統進行研究,以驗證ESINDy 方法處理奇異性問題的有效性.

1.3 R?sler 系統識別

R?sler 系統作為混沌系統,識別結果的微小差異會顯著反映在系統中,并且標準R?sler 系統的計算集具有一定的奇異性,這無疑對SINDy 算法的魯棒性和奇異問題的識別能力提出了新的挑戰.

數學上常使用條件數判定矩陣奇異及病態程度,矩陣奇異性越強,條件數越大.R?sler 系統數值解狀態矩陣條件數為 3.227×103,具有一定病態矩陣的性質.分別使用SINDy 與ESINDy 識別該系統,x分量識別結果如表2 所示,關鍵參數識別結果如表3 所示.

表2 x 分量識別結果Table 2 x component identification result

表3 部分參數識別結果Table 3 Partial parameter identification results

由表3 可知,ESINDy 較SINDy 具有更高的識別精度.就識別結果中稀疏項的正確率而言,ESINDy(8 項,100%)較SINDy (6 項,75%)更高,具有更好的稀疏效果.

如圖2 所示,使用兩種方法分別重構系統并計算x分量數值解,ESINDy (90.71%)的準確率遠高于SINDy (44.59%),證明ESINDy 在面對奇異問題時擁有比SINDy 更優越的性能.

圖2 重建系統x 分量數值解Fig.2 Numerical solution of x component of reconstruction system

2 實驗模型

2.1 非線性結構設計

選用一類雙薄膜式雙穩態振動系統作為算法的實際應用范例,結構示意圖如圖3 所示.該結構主要由薄膜、中間質量體、外殼和外部磁鐵組成,其中:質量體攜帶磁鐵通過切割鑲嵌在外殼中的線圈產生電流,兩端通過磁鐵固定在薄膜上,并與外部磁鐵相吸引;上下端蓋可通過旋轉改變外部磁鐵與薄膜上磁鐵間距,從而改變磁力大小.實驗與仿真結果都表明,對常見的單層薄膜結構而言,雙薄膜的設計可以有效地避免由于質量體周向旋轉而產生的無關振型.

圖3 雙薄膜式雙穩態振動系統結構Fig.3 Double-film-type bistable vibration system

質量體的振動微分方程可表示為

其中,c()表示阻尼項,k(x)表示彈性項,e(x)表示磁力項,h(t)表示激勵項.

2.2 非線性項分析

對強非線性耦合系統而言,非線性項會極大影響系統的非線性響應及動力學行為[30],因此識別非線性項是識別系統方程的重中之重.系統控制方程中非線性成分主要由彈力、磁力提供,使用板殼理論分析彈力項.撓度方程可表示為

其中,ω 為撓度,q(ρ,θ)為極坐標表示下的載荷,D為抗彎剛度.實驗模型可簡化為中間受集中載荷環形模型,積分求解式(18),最大撓度可表示為

其中,P為作用于中間區域的集中載荷,R為圓環半徑,孔徑比 ξ=r/R.該模型最大撓度為

通過對比式(19)和式(20),發現 ωm與 ωp成正比.因此,可由中間孔的尺寸確定撓度改變的程度ωm=G(ξ)ωp,定義無量綱偏移系數

薄膜恢復力與撓度的關系為

由于薄膜為大變形,力-變關系實際滿足如下關系

當 ω <ωL時,式(23)所述關系可以認為是小變形.3 次擬合關系近似斜率與線性關系斜率應滿足下述條件

撓度小于半徑的1% 可以認為是小變形,即ωL=0.01R,則恢復力與薄膜變形關系可表示為

振子質量為M,故彈力項可表示為

對于磁力項,由磁偶極子理論[31-32],一對圓形磁鐵間磁力Fe(l)可表示為

其中,C為與永磁體剩余磁通密度、空間磁導率、永磁體體積有關的常數,a和b為與永磁體半徑和厚度有關的常數.EMH 具有兩對距離相同的磁鐵,假設在振子位移為x時,磁鐵間磁力可表示為

其中,l為兩端磁鐵間距離.磁力項e(x)通過泰勒展開可表示為3 次多項式形式

其中,M為振子質量,c1和c3為未知定常數.磁鐵力矩關系可表示為

其中,f(x)為不包含關于x的1 次和3 次項的未知函數,a1和a3為待擬合系數且有如下關系

通過測量并擬合單對磁鐵力-距關系的a1和a3并代入式(31),即可獲得c1和c3,并由式(29)得到磁力項.

2.3 參數測量及仿真

在對EMH 系統進行識別前,需要先對其動力學建模后的結果進行實驗與仿真,以驗證動力學建模結果的正確性,便于與數據驅動建模結果進行比較.

針對磁力項,使用COMSOL Multiphysics 5.X進行仿真計算,并搭建如圖4(a)所示實驗系統,測量磁力與磁鐵間距離的關系.

圖4 磁力及彈力參數測量裝置 (1.測力計,2.手輪,3.磁鐵,4.EMH)Fig.4 Magnetic force and elastic force parameters measuring device(1.dynamometer,2.rocker,3.magnet,4.EMH)

實驗與理論、仿真結果比對結果如圖5 所示.實驗中兩磁鐵距離在 5 ～10 mm 間變化,該區域實驗和仿真數據較為接近,因此對距離l>5 mm 的仿真數據進行3 次多項式擬合

圖5 磁力的實驗、仿真與理論結果Fig.5 Experimental,simulation and theoretical results of magnetic force

將式(32)代入式(29),磁力項可表示為

針對彈力項,使用Ansys Workbench 進行有限元仿真,幾何參數與仿真參數見表4,仿真模型為周邊夾支邊界條件的圓環模型.

表4 Ansys 仿真采用的模型參數Table 4 Model parameters used in Ansys simulation

EMH 可視為孔徑比 ξ=2:13 的圓環,根據式(25),恢復力與薄膜中心點位移應滿足

此處考慮內部氣體壓強影響,則導致力-變關系存在一定偏移

其中,假設內部空氣壓強比大氣壓高1%,則中心偏移量b≈2 mm,此時力-變關系可表示為

搭建如圖4(b)所示的實驗系統,測量恢復力與中央磁鐵位移的關系,并將測試結果與理論、仿真結果比對,如圖6 所示.首先,無偏移理論與仿真數據吻合,這說明從理論角度對于薄膜結構彈力的分析較為準確;其次,在考慮中心偏移的情況下,實驗測試數據與式(36)的力-變形關系基本吻合,也證明了偏移假設的必要性.因此,后續分析以式(36)的彈力關系作為實際表達式.

圖6 恢復力實驗、仿真與理論結果Fig.6 Results of experiment,simulation and theory of resilience

3 實驗

搭建如圖7 所示振動測量實驗系統,使用位移傳感器測量振子振幅,加速度計測量樣機外殼加速度.測量得到的振幅和加速度數據信噪比約為31 dB,為混雜了一定強度隨機噪聲的采集信號.使用振幅和加速度數據構造狀態矩陣,狀態矩陣條件數為 1.056 3×1013,表現出病態矩陣的性質.

圖7 振動實驗測試系統Fig.7 Vibration experimental system

分別使用SINDy,ESINDy 兩種分析方法進行數據驅動建模,兩種方法的識別結果和理論結果的對比如表5 所示.

表5 理論和識別結果對比Table 5 Comparison between theory and identification results

相比較傳統理論手段建模,數據驅動方法具有精度更高,能夠識別到難以識別的阻尼和非線性項等優勢;相比較單一使用數據驅動方法建模,理論與數據驅動相結合的方法能夠預先估計可能含有的項及部分項的系數,具有更高的精度和更寬的應用場景.

比較兩種數據驅動方法,從識別結果來看,ESINDy的識別準確率(83.07%)相較于SINDy (62.48%)有顯著提高.就稀疏項的識別準確率而言,ESINDy 正確識別稀疏項數(4 項)與SINDy (4 項)相仿,ESINDy 錯誤識別稀疏項數(0 項)低于SINDy(1 項),表明SINDy 面對病態矩陣確實存在過度壓縮現象,而ESINDy 則沒有此類問題.

從系統重構效果來看,SINDy 對系數篩選的錯誤導致相圖中重要信息的丟失.如圖8 所示,ESINDy重構相圖表現出與真實相圖類似的馬鞍形狀,且鞍點與穩定點的位置幾乎相同,而SINDy 重構相圖僅表現出一定的非線性特點.

圖8 振動系統真實相圖與SINDy,ESINDy 重構系統相圖Fig.8 Real phase diagram of vibration system,SINDy and ESINDy reconstruction system

4 結論

本文針對傳統SINDy 方法難以處理強噪聲數據和奇異系統的缺點進行改進,提出改進型ESINDy算法.改進主要集中在將曲率濾波與SINDy 進行結合,從而提供更適合的正則化參數,以獲得更好的降噪效果;使用ElasticNet 作為SINDy 的循環主函數體,從而大大增強了SINDy 對奇異系統的處理能力.使用ESINDy 算法識別R?sler 模型和EMH 模型,研究得出如下結論.

(1)通過對R?sler 系統的分析,證明相較于SINDy,ESINDy 在稀疏項篩選精度和有效項識別精度上均有顯著改善,精度提升約15%～20%,說明改進方法有效提高了針對奇異問題的分析能力與分析精度.

(2)通過對EMH 系統的分析,證明數據驅動方法具備重建控制方程的能力,其中: 基于對重構系統的數值計算與分析,系統的動力學特點得以更加深入而準確地呈現,即ESINDy 能夠準確識別出EMH 系統的鞍點、穩定點、周期、振幅等隱含信息,且避免了諸如過度壓縮等問題,表現出比SINDy更優越的性能.

(3)稀疏辨識算法既可以識別理論分析方法無法計算的參數,又可以完善系統的控制方程.隨著研究對象日益復雜,利用稀疏辨識算法結合理論方法對物理模型進行預測的手段具備潛在優勢.

在后續研究中,將發揮ESINDy 方法的優勢,對于EMH 的復雜振動問題進行進一步的深入分析.