一類雙層薄膜結(jié)構(gòu)振動(dòng)能量采集器的數(shù)據(jù)驅(qū)動(dòng)建模方法及應(yīng)用1)

邱宏蘊(yùn) 王志霞 丁 北 王 煒

(天津大學(xué)機(jī)械工程學(xué)院,天津 300350)

引言

在人們觀察世界,開展生產(chǎn)實(shí)踐,預(yù)測(cè)實(shí)踐結(jié)果的過程中,建模方法發(fā)揮著至關(guān)重要的作用.傳統(tǒng)建模方法通過對(duì)基礎(chǔ)理論的梳理歸納以及數(shù)學(xué)演繹,得到適用于當(dāng)前問題的模型.然而,該方法極度依賴先驗(yàn)知識(shí),對(duì)于線性、直觀問題的分析較為有效,而對(duì)于包含強(qiáng)非線性因素的系統(tǒng)則適用性有限,也因此為基于實(shí)驗(yàn)數(shù)據(jù)的建模方法提供了成長(zhǎng)空間.

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

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

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

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

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

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

1 識(shí)別算法

1.1 SINDy 算法及數(shù)據(jù)濾波算法

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

(1)數(shù)據(jù)濾波算法

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

傳統(tǒng)正則化方法的正則化項(xiàng)能量會(huì)隨著優(yōu)化進(jìn)程的進(jìn)行而逐漸降低,而正則化項(xiàng)作為優(yōu)化主體,能量的降低會(huì)導(dǎo)致優(yōu)化效果變差,曲率濾波使用已知的曲率來優(yōu)化其對(duì)應(yīng)的正則項(xiàng),從而達(dá)到更好的優(yōu)化效果.對(duì)一維含噪信號(hào)x(t),曲率可表示為

其中,Ix和Ixx分別表示x(t)對(duì)t的1 階、2 階導(dǎo)數(shù).對(duì)一維含噪信號(hào),曲率表示的正則化項(xiàng)可表示為

其中,λ為正則化參數(shù).在dt的步長(zhǎng)下,曲率優(yōu)化更新過程可表示為

正則化方法依賴正則化參數(shù)的選擇[27],然而在實(shí)際問題研究過程中,通常需要先進(jìn)行數(shù)據(jù)濾波,隨后使用濾波后的數(shù)據(jù)進(jìn)行SINDy 數(shù)據(jù)驅(qū)動(dòng).這就導(dǎo)致正則化參數(shù)不能隨著數(shù)據(jù)驅(qū)動(dòng)問題的變化而自動(dòng)選取,限制了SINDy 的工程場(chǎng)景應(yīng)用[28].

(2)SINDy 算法

假設(shè)振動(dòng)系統(tǒng)的方程可以表示為如下形式

其中,x(t)為系統(tǒng)在t時(shí)刻的狀態(tài),f(x(t))為系統(tǒng)控制方程.矩陣 Θ 為由給定非線性函數(shù)構(gòu)成的基函數(shù)庫

根據(jù)基函數(shù)庫,構(gòu)造狀態(tài)矩陣A

假定控制方程f(x(t))可由基函數(shù)庫中的部分項(xiàng)構(gòu)成,Ξ 為基函數(shù)庫中各成分系數(shù)構(gòu)成的矩陣,則振動(dòng)系統(tǒng)方程可表示為

因此,在已知導(dǎo)數(shù)矩陣與狀態(tài)矩陣A的前提下,對(duì)振動(dòng)系統(tǒng)方程的識(shí)別問題可轉(zhuǎn)化為對(duì)的識(shí)別問題.該問題可通過最小二乘法(OLS)解決,目標(biāo)函數(shù)Target可以表示為均方誤差形式

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

Lasso 法[2]通過在最小二乘法目標(biāo)函數(shù)后增加L1 正則化項(xiàng)的方法,達(dá)到稀疏效果.稀疏方法可以簡(jiǎn)化復(fù)雜的數(shù)學(xué)模型,得到項(xiàng)數(shù)較少的簡(jiǎn)單模型來表征系統(tǒng)的動(dòng)力學(xué)特性,并綜合考慮模型壓縮程度和模型精確性之間的矛盾.Lasso 法的求解目標(biāo)可表示為

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

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

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

1.2 ESINDy 算法

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

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

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

(1)數(shù)據(jù)驅(qū)動(dòng)部分

使用L1 正則化項(xiàng)的Lasso 能夠?qū)崿F(xiàn)良好的稀疏效果,但對(duì)病態(tài)數(shù)據(jù)的識(shí)別能力存在顯著不足,導(dǎo)致SINDy 在處理此類數(shù)據(jù)時(shí)出現(xiàn)精度降低、壓縮過度等問題.使用L2 正則化項(xiàng)的嶺回歸,雖然無法得到精簡(jiǎn)模型,卻能在處理上述數(shù)據(jù)時(shí)具有更好的魯棒能力,實(shí)現(xiàn)較高精度的識(shí)別.EN 同時(shí)使用L1 和L2 正則化項(xiàng),巧妙地結(jié)合了L1 的稀疏能力與L2 的魯棒能力,因而具有更好的識(shí)別效果.EN 的等價(jià)優(yōu)化問題為

經(jīng)推導(dǎo)可知閉式解可表示為

(2)數(shù)據(jù)濾波部分

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

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

表1 ESINDy 算法流程Table 1 ESINDy algorithm

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

1.3 R?sler 系統(tǒng)識(shí)別

R?sler 系統(tǒng)作為混沌系統(tǒng),識(shí)別結(jié)果的微小差異會(huì)顯著反映在系統(tǒng)中,并且標(biāo)準(zhǔn)R?sler 系統(tǒng)的計(jì)算集具有一定的奇異性,這無疑對(duì)SINDy 算法的魯棒性和奇異問題的識(shí)別能力提出了新的挑戰(zhàn).

數(shù)學(xué)上常使用條件數(shù)判定矩陣奇異及病態(tài)程度,矩陣奇異性越強(qiáng),條件數(shù)越大.R?sler 系統(tǒng)數(shù)值解狀態(tài)矩陣條件數(shù)為 3.227×103,具有一定病態(tài)矩陣的性質(zhì).分別使用SINDy 與ESINDy 識(shí)別該系統(tǒng),x分量識(shí)別結(jié)果如表2 所示,關(guān)鍵參數(shù)識(shí)別結(jié)果如表3 所示.

表2 x 分量識(shí)別結(jié)果Table 2 x component identification result

表3 部分參數(shù)識(shí)別結(jié)果Table 3 Partial parameter identification results

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

如圖2 所示,使用兩種方法分別重構(gòu)系統(tǒng)并計(jì)算x分量數(shù)值解,ESINDy (90.71%)的準(zhǔn)確率遠(yuǎn)高于SINDy (44.59%),證明ESINDy 在面對(duì)奇異問題時(shí)擁有比SINDy 更優(yōu)越的性能.

圖2 重建系統(tǒng)x 分量數(shù)值解Fig.2 Numerical solution of x component of reconstruction system

2 實(shí)驗(yàn)?zāi)Ｐ?/h2>

2.1 非線性結(jié)構(gòu)設(shè)計(jì)

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

圖3 雙薄膜式雙穩(wěn)態(tài)振動(dòng)系統(tǒng)結(jié)構(gòu)Fig.3 Double-film-type bistable vibration system

質(zhì)量體的振動(dòng)微分方程可表示為

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

2.2 非線性項(xiàng)分析

對(duì)強(qiáng)非線性耦合系統(tǒng)而言,非線性項(xiàng)會(huì)極大影響系統(tǒng)的非線性響應(yīng)及動(dòng)力學(xué)行為[30],因此識(shí)別非線性項(xiàng)是識(shí)別系統(tǒng)方程的重中之重.系統(tǒng)控制方程中非線性成分主要由彈力、磁力提供,使用板殼理論分析彈力項(xiàng).撓度方程可表示為

其中,ω 為撓度,q(ρ,θ)為極坐標(biāo)表示下的載荷,D為抗彎剛度.實(shí)驗(yàn)?zāi)Ｐ涂珊?jiǎn)化為中間受集中載荷環(huán)形模型,積分求解式(18),最大撓度可表示為

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

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

薄膜恢復(fù)力與撓度的關(guān)系為

由于薄膜為大變形,力-變關(guān)系實(shí)際滿足如下關(guān)系

當(dāng) ω <ωL時(shí),式(23)所述關(guān)系可以認(rèn)為是小變形.3 次擬合關(guān)系近似斜率與線性關(guān)系斜率應(yīng)滿足下述條件

撓度小于半徑的1% 可以認(rèn)為是小變形,即ωL=0.01R,則恢復(fù)力與薄膜變形關(guān)系可表示為

振子質(zhì)量為M,故彈力項(xiàng)可表示為

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

其中,C為與永磁體剩余磁通密度、空間磁導(dǎo)率、永磁體體積有關(guān)的常數(shù),a和b為與永磁體半徑和厚度有關(guān)的常數(shù).EMH 具有兩對(duì)距離相同的磁鐵,假設(shè)在振子位移為x時(shí),磁鐵間磁力可表示為

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

其中,M為振子質(zhì)量,c1和c3為未知定常數(shù).磁鐵力矩關(guān)系可表示為

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

通過測(cè)量并擬合單對(duì)磁鐵力-距關(guān)系的a1和a3并代入式(31),即可獲得c1和c3,并由式(29)得到磁力項(xiàng).

2.3 參數(shù)測(cè)量及仿真

在對(duì)EMH 系統(tǒng)進(jìn)行識(shí)別前,需要先對(duì)其動(dòng)力學(xué)建模后的結(jié)果進(jìn)行實(shí)驗(yàn)與仿真,以驗(yàn)證動(dòng)力學(xué)建模結(jié)果的正確性,便于與數(shù)據(jù)驅(qū)動(dòng)建模結(jié)果進(jìn)行比較.

針對(duì)磁力項(xiàng),使用COMSOL Multiphysics 5.X進(jìn)行仿真計(jì)算,并搭建如圖4(a)所示實(shí)驗(yàn)系統(tǒng),測(cè)量磁力與磁鐵間距離的關(guān)系.

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

實(shí)驗(yàn)與理論、仿真結(jié)果比對(duì)結(jié)果如圖5 所示.實(shí)驗(yàn)中兩磁鐵距離在 5 ～10 mm 間變化,該區(qū)域?qū)嶒?yàn)和仿真數(shù)據(jù)較為接近,因此對(duì)距離l>5 mm 的仿真數(shù)據(jù)進(jìn)行3 次多項(xiàng)式擬合

圖5 磁力的實(shí)驗(yàn)、仿真與理論結(jié)果Fig.5 Experimental,simulation and theoretical results of magnetic force

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

針對(duì)彈力項(xiàng),使用Ansys Workbench 進(jìn)行有限元仿真,幾何參數(shù)與仿真參數(shù)見表4,仿真模型為周邊夾支邊界條件的圓環(huán)模型.

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

EMH 可視為孔徑比 ξ=2:13 的圓環(huán),根據(jù)式(25),恢復(fù)力與薄膜中心點(diǎn)位移應(yīng)滿足

此處考慮內(nèi)部氣體壓強(qiáng)影響,則導(dǎo)致力-變關(guān)系存在一定偏移

其中,假設(shè)內(nèi)部空氣壓強(qiáng)比大氣壓高1%,則中心偏移量b≈2 mm,此時(shí)力-變關(guān)系可表示為

搭建如圖4(b)所示的實(shí)驗(yàn)系統(tǒng),測(cè)量恢復(fù)力與中央磁鐵位移的關(guān)系,并將測(cè)試結(jié)果與理論、仿真結(jié)果比對(duì),如圖6 所示.首先,無偏移理論與仿真數(shù)據(jù)吻合,這說明從理論角度對(duì)于薄膜結(jié)構(gòu)彈力的分析較為準(zhǔn)確;其次,在考慮中心偏移的情況下,實(shí)驗(yàn)測(cè)試數(shù)據(jù)與式(36)的力-變形關(guān)系基本吻合,也證明了偏移假設(shè)的必要性.因此,后續(xù)分析以式(36)的彈力關(guān)系作為實(shí)際表達(dá)式.

圖6 恢復(fù)力實(shí)驗(yàn)、仿真與理論結(jié)果Fig.6 Results of experiment,simulation and theory of resilience

3 實(shí)驗(yàn)

搭建如圖7 所示振動(dòng)測(cè)量實(shí)驗(yàn)系統(tǒng),使用位移傳感器測(cè)量振子振幅,加速度計(jì)測(cè)量樣機(jī)外殼加速度.測(cè)量得到的振幅和加速度數(shù)據(jù)信噪比約為31 dB,為混雜了一定強(qiáng)度隨機(jī)噪聲的采集信號(hào).使用振幅和加速度數(shù)據(jù)構(gòu)造狀態(tài)矩陣,狀態(tài)矩陣條件數(shù)為 1.056 3×1013,表現(xiàn)出病態(tài)矩陣的性質(zhì).

圖7 振動(dòng)實(shí)驗(yàn)測(cè)試系統(tǒng)Fig.7 Vibration experimental system

分別使用SINDy,ESINDy 兩種分析方法進(jìn)行數(shù)據(jù)驅(qū)動(dòng)建模,兩種方法的識(shí)別結(jié)果和理論結(jié)果的對(duì)比如表5 所示.

表5 理論和識(shí)別結(jié)果對(duì)比Table 5 Comparison between theory and identification results

相比較傳統(tǒng)理論手段建模,數(shù)據(jù)驅(qū)動(dòng)方法具有精度更高,能夠識(shí)別到難以識(shí)別的阻尼和非線性項(xiàng)等優(yōu)勢(shì);相比較單一使用數(shù)據(jù)驅(qū)動(dòng)方法建模,理論與數(shù)據(jù)驅(qū)動(dòng)相結(jié)合的方法能夠預(yù)先估計(jì)可能含有的項(xiàng)及部分項(xiàng)的系數(shù),具有更高的精度和更寬的應(yīng)用場(chǎng)景.

比較兩種數(shù)據(jù)驅(qū)動(dòng)方法,從識(shí)別結(jié)果來看,ESINDy的識(shí)別準(zhǔn)確率(83.07%)相較于SINDy (62.48%)有顯著提高.就稀疏項(xiàng)的識(shí)別準(zhǔn)確率而言,ESINDy 正確識(shí)別稀疏項(xiàng)數(shù)(4 項(xiàng))與SINDy (4 項(xiàng))相仿,ESINDy 錯(cuò)誤識(shí)別稀疏項(xiàng)數(shù)(0 項(xiàng))低于SINDy(1 項(xiàng)),表明SINDy 面對(duì)病態(tài)矩陣確實(shí)存在過度壓縮現(xiàn)象,而ESINDy 則沒有此類問題.

從系統(tǒng)重構(gòu)效果來看,SINDy 對(duì)系數(shù)篩選的錯(cuò)誤導(dǎo)致相圖中重要信息的丟失.如圖8 所示,ESINDy重構(gòu)相圖表現(xiàn)出與真實(shí)相圖類似的馬鞍形狀,且鞍點(diǎn)與穩(wěn)定點(diǎn)的位置幾乎相同,而SINDy 重構(gòu)相圖僅表現(xiàn)出一定的非線性特點(diǎn).

圖8 振動(dòng)系統(tǒng)真實(shí)相圖與SINDy,ESINDy 重構(gòu)系統(tǒng)相圖Fig.8 Real phase diagram of vibration system,SINDy and ESINDy reconstruction system

4 結(jié)論

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

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

(2)通過對(duì)EMH 系統(tǒng)的分析,證明數(shù)據(jù)驅(qū)動(dòng)方法具備重建控制方程的能力,其中: 基于對(duì)重構(gòu)系統(tǒng)的數(shù)值計(jì)算與分析,系統(tǒng)的動(dòng)力學(xué)特點(diǎn)得以更加深入而準(zhǔn)確地呈現(xiàn),即ESINDy 能夠準(zhǔn)確識(shí)別出EMH 系統(tǒng)的鞍點(diǎn)、穩(wěn)定點(diǎn)、周期、振幅等隱含信息,且避免了諸如過度壓縮等問題,表現(xiàn)出比SINDy更優(yōu)越的性能.

(3)稀疏辨識(shí)算法既可以識(shí)別理論分析方法無法計(jì)算的參數(shù),又可以完善系統(tǒng)的控制方程.隨著研究對(duì)象日益復(fù)雜,利用稀疏辨識(shí)算法結(jié)合理論方法對(duì)物理模型進(jìn)行預(yù)測(cè)的手段具備潛在優(yōu)勢(shì).

在后續(xù)研究中,將發(fā)揮ESINDy 方法的優(yōu)勢(shì),對(duì)于EMH 的復(fù)雜振動(dòng)問題進(jìn)行進(jìn)一步的深入分析.