廣義概率密度演化方程的相空間重構解法

蔣仲銘，劉毅，胥鵬，王成龍，王思琦，洪曉林

西南大學 工程技術學院，重慶 400715

由于在工程實踐中隨機性廣泛地存在于材料特性、外部荷載和邊界條件等因素中，因此隨機動力系統分析作為現代力學研究的重要分支，一直以來是人們研究的熱點．文獻[1]從概率守恒的基本原理出發，發展了概率密度演化分析理論，建立了廣義概率密度演化方程(GDEE)，為復雜非線性隨機動力系統的分析提供了一條可行的途徑[2-4]．已有的工作證明：這一新的分析理論在線性與非線性系統的隨機動力反應分析[5]、結構動力可靠度計算[6-7]、結構隨機最優控制等[8]方面均可獲得高效、準確的分析結果．

廣義概率密度演化方程作為一類偏微分方程，近年來眾多學者發展了一系列解析和半解析的求解方法[9-10]對其進行求解．本文從隨機動力系統的相空間出發，提出了一種可用于精確、高效求解非線性系統的GDEE求解方法——相空間重構法(PSRM)．本文利用該方法，得到了若干典型非線性振子的概率密度解，并將其與Monte Carlo模擬方法[11]、有限差分法[9]等數值求解方法進行了比較，驗證了該方法的準確性、有效性和便捷性．

1 廣義密度演化方程

文獻[5]發展了一類用于隨機系統分析的廣義概率密度演化方程，經過十數年的發展，這一方法廣泛地應用于結構的可靠性分析[5]、鐵路軌道設計[12]以及城市基礎設施防震減災設計[13-14]等領域之中．

不失一般性，設隨機動力系統為

(1)

初始條件為

x0=(x1，x2，…，xn)T

(2)

其中：X0=(X1，X2，…，Xn)T，n是物理系統的維數；G(·)為一般非線性系統，Θ=(Θ1，Θ2，…，Θs)T，其聯合概率密度函數為pΘ(θ)，s是系統中隨機變量的總數．顯然，對于特定的隨機動力系統，方程(1)，(2)的解存在且唯一，并連續地依賴于初始條件．即對于給定的初始條件，式(1)的解可寫為

X=H(Θ，t)

(3)

由概率守恒原理，存在廣義概率密度演化方程[2-3，15]

(4)

初始條件為

pXΘ(x，θ，t)=δ(x-x0)pΘ(θ)

(5)

由于廣義密度演化方程的解析解很難得到，工程中常用其數值解．

2 相空間重構法

在大多數非線性問題乃至工程問題中，式(4)和式(5)的解析求解都是相當困難的．由此，大量學者發展了一系列數值和近似的解析算法用于求解廣義概率密度演化方程[16-17]．本文通過深入研究廣義概率密度演化理論的物理機制，在此基礎上提出了一類基于非線性物理系統分解的相空間重構方法，可以有效避免采用傳統點演化方法求解廣義概率密度演化方程時所遇到的網格敏感性問題和數值耗散問題．

使用特征線法[18]得方程(4)和(5)的形式解為[19]

pXΘ(x，θ，t)=δ(x-H(θ，t))pΘ(θ)

(6)

不妨令

G(θ，x，t)=x-H(θ，t)

(7)

則方程(6)可化為

pXΘ(x，θ，t)=δ[G(θ，x，t)]pΘ(θ)

(8)

當a≠0時，存在

(9)

對一般的情況而言，關于g(x)的δ函數已在文獻[20]中給出

(10)

其中xi是g(x)的一個根．將式(10)代入式(8)的形式解可得到

(11)

(12)

聯立式(8)，(11)和(12)，對θ求積分可得

(13)

因此，在相空間X-θ的角度處理廣義概率密度演化方程(GDEE)的求解問題避免了反函數求解這一數學上的難點；同時，由于在實際中我們往往主要關注一個物理量(如位移)或是兩個物理量(位移和速度)的概率密度函數，這也大大簡化了X-θ空間中的G(θ，x，t)曲線，使得對物理系統G(θ，x，t)的分解成為可能．

基于物理系統分解的數值求解方法步驟如下：

2)對X-θ空間中的物理系統進行分類，將分布空間ΩΘ分為ΩΘ1，ΩΘ2，…，ΩΘNcla；Ncla為物理系統的分類個數；

3 數值算例

3.1 SDOF振子

單自由度振子(Single Degree of Freedom Oscillator)的無阻尼自由振動是一種簡諧振動，其固有頻率是系統本身的性質，與初始條件無關，且它的速度與加速度也是簡諧的．其方程可以寫為

(14)

其初始條件為

(15)

隨機變量ω在區間[0.1，0.4]內均勻分布，由于該方程為線性微分方程，其解析解易得為

X(ω，t)=0.1cos(ωt)

(16)

圖1為使用PSRM分解物理系統后的各個概率密度函數，將各物理系統的概率密度函數疊加得到圖2．

圖1 SDOF振子在各物理系統下的概率密度函數

圖2 SDOF振子各物理系統下的概率密度函數疊加

通過10 000次Monte Carlo模擬與本文提出的相空間重構法進行比較(圖3)，可以看出本文所提出的PSRM與精確解的吻合程度較好．圖4給出了代表點個數相同(200個)的情況下，使用Monte Carlo模擬方法、PDEM和 PSRM在典型時刻(t=20)響應的概率密度函數．由表1失效概率數據可以看出，相比傳統的PDEM，PSRM使用較少的代表點個數即可獲得較高的計算精度．

圖3 SDOF振子PSRM與10 000次Mento Carlo效果對比

圖4 SDOF振子PSRM、200次Mento Carlo方法與PDEM比較

表1 SDOF振子在4種方法中不同閾值下的失效概率

3.2 Riccati振子

Riccati振子是控制理論中的重要方程，涉及Kalman濾波、信號頻率跟蹤[23]、線性二次調節器問題及模型簡化等諸多問題[24-25]．設一階Riccati方程為

(17)

其初始條件為

X(0)=1

(18)

圖5為使用PSRM分解Riccati振子得到的概率密度函數，疊加概率密度函數得到圖6．

圖5 Riccati振子在各物理系統下的概率密度函數

圖6 Riccati振子各物理系統的概率密度函數疊加

隨機變量θ在區間[1，10]均勻分布．用10 000次Monte Carlo模擬與本文提出的相空間重構法進行比較(圖7)，可以看出本文所提出的PSRM與精確解的吻合程度較好．圖8給出了代表點個數相同(200個)的情況下，使用Monte Carlo模擬、有限差分方法PDEM和 PSRM在典型時刻(t=20)響應的概率密度函數．由表2失效概率數據可以看出，相比傳統的PDEM，PSRM使用較少的代表點個數即可獲得較高的計算精度．

圖7 PSRM與10 000次Mento Carlo方法在Riccati方程上的對比實驗

圖8 Riccati方程PSRM與200次Mento Carlo 及PDEM比較

表2 Riccati方程在4種方法中不同閾值下的失效概率

3.3 Van der Pol振子

Van der Pol振子是非線性系統的經典模型之一，它起源于范德波爾對電子電路中三極管的振蕩效應的研究[26]．這一模型在物理學、生物學、神經學甚至經濟學中，都有著廣泛的應用[27-28]．設一維Van der Pol方程為

(19)

初始條件為

(20)

使用PSRM分解Van der pol振子的物理系統(圖9)，疊加各個物理系統的概率密度函數得到圖10．

圖9 Van der pol振子在各物理系統下的概率密度函數

圖10 Van der pol振子各物理系統的概率密度函數疊加

隨機變量θ在區間[0.1，5]均勻分布．用104次Monte Carlo模擬與本文提出的相空間重構法進行比較(圖11)，可以看出，本文所提出的PSRM與精確解的吻合程度較好．圖12給出了代表點個數相同(200個)的情況下，使用Monte Carlo模擬、PDEM和 PSRM在典型時刻(t=20)響應的概率密度函數．由表3失效概率數據可以看出，相比傳統的PDEM，PSRM使用較少的代表點個數即可獲得較高的計算精度．

圖11 Van der Pol振子PSRM與10 000次Mento Carlo比較

圖12 Van der pol振子PSRM與200次Mento Carlo及PDEM比較

表3 Van der pol振子在4種方法中不同閾值下的失效概率

3.4 Duffing振子

Duffing振子是一個典型的非線性振動系統，盡管是從簡單物理模型中得出來的非線性振動模型，但是其模型具有代表性．工程實際中的許多非線性振動問題的數學模型都可以轉化為該方程來研究，如腦電識別[29]、機械臂振動控制[30]、微弱信號檢測[31]等，Duffing系統也非常廣泛地被應用到實際工程中，例如尖銳碰摩轉子的故障檢測、微弱周期信號檢測、電力系統周期振蕩分析、周期電路系統的模擬與控制等[32-34]．設一維Duffing方程為

(21)

初始條件為

(22)

用PSRM分解Duffing振子的14個物理系統，得到各物理系統的概率密度函數(圖13)，并疊加各系統的概率密度函數(圖14)．

圖13 Duffing振子在各物理系統下的概率密度函數

圖14 Duffing振子各物理系統的概率密度疊加函數

隨機變量a在區間[0.1，5]均勻分布．用Monte Carlo模擬與本文提出的相空間重構法進行比較(圖15)，可以看出本文所提出的PSRM與精確解的吻合程度較好．圖16給出了代表點個數相同(600個)的情況下，使用Monte Carlo模擬、PDEM和 PSRM在典型時刻(t=20)響應的概率密度函數．由表4失效概率數據可以看出，相比傳統的PDEM，PSRM使用較少的代表點個數即可獲得較高的計算精度．

圖15 Duffing振子PSRM與10 000次Mento Carlo比較

圖16 Duffing振子PSRM與600次Mento Carlo及PDEM比較

表4 Duffing方程在4種方法中不同閾值下的失效概率

4 總結

廣義概率密度演化方程作為求解隨機動力系統的一種通用方法，近年來獲得了廣泛的關注．傳統的廣義概率密度演化方程數值算法，往往需要提前對概率空間進行剖分，通過獲取廣義概率密度演化方法的點演化方程來求解．本文通過深入研究廣義概率密度演化理論的物理機制，對代表物理系統的相軌跡進行分解，然后對分解后的廣義概率密度演化方程分別進行求解，疊加后即可得到原非線性系統的概率密度解．相比于傳統的數值求解方法，本方法不需要求解代表點的賦得概率，且所需的代表點個數較少．最后，通過4個數值算例，證明了該方法可以有效地計算強非線性系統隨機響應．

本文結果不僅可以作為求解廣義概率密度演化方程的一種半解析方法，同時也證明了廣義概率密度演化方程既可以從概率空間剖分角度進行求解，也可以從物理系統分解角度進行求解，有助于進一步深入研究和剖析隨機動力系統中隨機性與非線性的耦合機制．