分?jǐn)?shù)階混沌系統(tǒng)動(dòng)力學(xué)特性分析與DSP實(shí)現(xiàn)

劉天明，閻慧臻，馬晨光，楊飛飛，曹穎鴻

(大連工業(yè)大學(xué) 信息科學(xué)與工程學(xué)院,遼寧 大連 116034)

0 引 言

近年來(lái)的分?jǐn)?shù)階微積分的相關(guān)研究表明，相比于整數(shù)階微積分，分?jǐn)?shù)階微積分可以更好地描述客觀物理世界[1]，因此針對(duì)分?jǐn)?shù)階微積分的研究成為熱點(diǎn)。針對(duì)Chen系統(tǒng)[2]、Liu系統(tǒng)[3]和簡(jiǎn)化Lorenz系統(tǒng)[4]等非線性混沌系統(tǒng)的研究發(fā)現(xiàn)，相比于整數(shù)階混沌系統(tǒng)，分?jǐn)?shù)階混沌系統(tǒng)具有更為復(fù)雜的動(dòng)力學(xué)特性[5-6]。

目前，分?jǐn)?shù)階混沌系統(tǒng)的求解算法主要有頻域法(FDM)[7]、預(yù)估校正法(ABM)[8]和Adomian分解法(ADM)[9]。這些算法普遍采用的是黎曼-劉維爾定義和卡普托定義，但它們都存在一些問(wèn)題：首先是其不能滿足整數(shù)階微積分所滿足的一些重要性質(zhì)[10]，例如乘積法則和鏈?zhǔn)椒▌t；其次，計(jì)算過(guò)程非常復(fù)雜[11-12]。因此，Khalil等[13]提出了一種新的分?jǐn)?shù)微積分定義——可整合分?jǐn)?shù)微積分。該分?jǐn)?shù)階微積分定義和性質(zhì)與ADM算法相結(jié)合可以很好彌補(bǔ)現(xiàn)有算法求解分?jǐn)?shù)階微分方程的不足[14-15]。CADM算法改進(jìn)ADM算法，降低了計(jì)算的復(fù)雜程度，并具備較快的收斂速度、計(jì)算速度和較小的資源消耗等優(yōu)點(diǎn)[10]。CADM算法求解分?jǐn)?shù)階混沌系統(tǒng)正逐漸成為研究的熱點(diǎn)[16-18]。

分?jǐn)?shù)階混沌系統(tǒng)的應(yīng)用[19-21]實(shí)際依賴信號(hào)處理的軟硬件技術(shù)。數(shù)字信號(hào)處理器(DSP)以其性能優(yōu)越、處理方便等優(yōu)點(diǎn)在工程上得到了廣泛的應(yīng)用?；诖?，本研究利用DSP技術(shù)，對(duì)提出的四維分?jǐn)?shù)階混沌系統(tǒng)進(jìn)行了硬件實(shí)現(xiàn)。

本研究在一個(gè)四維Sprott-B混沌系統(tǒng)的基礎(chǔ)上，利用可整合微積分定義構(gòu)造該四維混沌系統(tǒng)的分?jǐn)?shù)階形式。基于CADM算法求取該分?jǐn)?shù)階系統(tǒng)的數(shù)值解，對(duì)其動(dòng)力學(xué)行為進(jìn)行了分析，利用0-1測(cè)試驗(yàn)證產(chǎn)生混沌最小階數(shù)，同時(shí)用SE和C0復(fù)雜度算法分析了該系統(tǒng)隨機(jī)性。最后運(yùn)用DSP技術(shù)對(duì)該系統(tǒng)進(jìn)行了硬件實(shí)現(xiàn)。

1 分?jǐn)?shù)階混沌系統(tǒng)數(shù)值分析

1.1 分?jǐn)?shù)階混沌系統(tǒng)CADM求解算法

設(shè)一個(gè)分?jǐn)?shù)階系統(tǒng)方程為

(1)

(2)

(3)

將系統(tǒng)非線性項(xiàng)進(jìn)行分解：

(4)

式中：i=0，1，2，…；j=1，2，…。則方程數(shù)值解為

(5)

式中：

(6)

1.2 分?jǐn)?shù)階混沌系統(tǒng)求解

在文獻(xiàn)[22]提出的Sprott-B系統(tǒng)的基礎(chǔ)上設(shè)計(jì)了一個(gè)新的四維自治混沌系統(tǒng)如式(7)所示。

(7)

式中：x，y，z，w為狀態(tài)變量；a，b，c為系統(tǒng)參數(shù)。

根據(jù)分?jǐn)?shù)階微積分定義，式(7)對(duì)應(yīng)的分?jǐn)?shù)階系統(tǒng)為

(8)

式中：q為系統(tǒng)階數(shù)。分解式(8)可得其中的線性項(xiàng)、非線性項(xiàng)和常數(shù)項(xiàng)分別為

(9)

(10)

(11)

根據(jù)式(4)對(duì)系統(tǒng)非線性項(xiàng)進(jìn)行分解，在保證精度的基礎(chǔ)上取前6項(xiàng)可得

(12)

(13)

(14)

設(shè)該分?jǐn)?shù)階混沌系統(tǒng)的初值為x0=[x1(0)，x2(0)，x3(0)，x4(0)]，展開(kāi)式第一項(xiàng)為

(15)

同時(shí)令

(16)

(17)

令h=t-t0為步長(zhǎng)，同時(shí)

(18)

則第二項(xiàng)可表示為

(19)

類似地，可得其他5項(xiàng)展開(kāi)式分別為

(20)

(21)

(22)

(23)

(24)

由此可得此時(shí)分?jǐn)?shù)階混沌系統(tǒng)解為

(25)

2 分?jǐn)?shù)階混沌系統(tǒng)動(dòng)力學(xué)特性分析

2.1 系統(tǒng)仿真

取式(7)系統(tǒng)參數(shù)a=4、b=9、c=5和q=0.8，步長(zhǎng)h=0.01，x0=[1，1，1，1]，系統(tǒng)相圖如圖1所示。此時(shí)李雅普諾夫指數(shù)為L(zhǎng)1=1.017 0、L2=0、L3=-12.547 9和L4=-16.720 6，系統(tǒng)維度DL=2.09。系統(tǒng)有一個(gè)正的李雅普諾夫指數(shù)，表明此時(shí)系統(tǒng)為混沌態(tài)。

此時(shí)系統(tǒng)Poincaré截面如圖2所示，圖中Poincaré截面既不是有限點(diǎn)集，也不是封閉曲線，是一些成片的具有分形結(jié)構(gòu)的密集點(diǎn)，這種結(jié)構(gòu)具備混沌系統(tǒng)的典型特征。

圖2 系統(tǒng)y-x平面的Poincaré截面(q=0.8)Fig.2 Poincaré section of the system on the y-x plane (q=0.8)

2.2 系統(tǒng)穩(wěn)定性分析

(26)

取a=4、b=9和c=5，此時(shí)系統(tǒng)平衡點(diǎn)為S1，2=(±3，±3，0，0)，系統(tǒng)特征方程為

λ32+9λ24+29λ16+126λ8+360=0

(27)

求解特征方程特征值，若其滿足

(28)

則系統(tǒng)在平衡點(diǎn)是漸進(jìn)穩(wěn)定的[23]。其中M是分?jǐn)?shù)階分母的最小公倍數(shù)。將特征值代入式(28)得

(29)

表明此時(shí)系統(tǒng)特征值滿足式(28)，所以系統(tǒng)在平衡點(diǎn)S1，2處是穩(wěn)定的。

2.3 系統(tǒng)參數(shù)變化對(duì)動(dòng)力學(xué)特性的影響

2.3.1 參數(shù)a對(duì)系統(tǒng)動(dòng)力學(xué)特性影響

取參數(shù)b=9、c=5和q=0.8，仿真步長(zhǎng)h=0.01，系統(tǒng)初值x0=[1，1，1，1]，當(dāng)a∈[1，5]時(shí)，系統(tǒng)的李雅普諾夫指數(shù)譜和分岔圖如圖3所示。為了更好觀察系統(tǒng)特性變化，在圖3(a)中舍去了最小的兩條李雅普指數(shù)曲線。當(dāng)a在[2.98，3.03]，[3.46，3.5]，[4.27，4.33]和[4.34，5]范圍內(nèi)取值時(shí)，系統(tǒng)最大的李雅普諾夫指數(shù)為零，所以a在這些區(qū)間內(nèi)取值系統(tǒng)是周期態(tài)。在其他范圍內(nèi)，最大的李雅普諾夫指數(shù)為正，系統(tǒng)是處于混沌態(tài)。表1總結(jié)了當(dāng)系統(tǒng)隨參數(shù)a變化時(shí)所具有的復(fù)雜動(dòng)力學(xué)行為。

(a) x-y平面相圖

(a) a=5.00

(a) 李雅普諾夫指數(shù)譜

表1 參數(shù)a變化的系統(tǒng)狀態(tài)Tab.1 The system status with a

由表1可知，參數(shù)a∈[1，5]時(shí)，系統(tǒng)中出現(xiàn)了一種典型的混沌吸引子和6種不同類型的周期態(tài)。圖4中分別取了周期1(a=5)、周期2(a=4.4)、周期5(a=4.16)和混沌態(tài)(a=4)時(shí)的相圖和0-1測(cè)試，圖中可以看到當(dāng)系統(tǒng)處于周期態(tài)時(shí)，p-s平面上是有界的規(guī)則的運(yùn)動(dòng)；當(dāng)系統(tǒng)處于混沌態(tài)時(shí)，p-s平面上是類似于布朗運(yùn)動(dòng)的無(wú)界運(yùn)動(dòng)。

2.3.2 參數(shù)b對(duì)系統(tǒng)動(dòng)力學(xué)特性影響

取參數(shù)a=4、c=5和q=0.8，仿真步長(zhǎng)h=0.01，系統(tǒng)初值x0=[1，1，1，1]，圖5是系統(tǒng)參數(shù)b∈[7，11]變化時(shí)的李雅普諾夫指數(shù)譜和分岔圖。圖5(a)同樣舍去最小的兩條李雅普諾夫指數(shù)曲線。從分岔圖中可以看到，隨著參數(shù)b增加，系統(tǒng)通過(guò)倍周期分岔從周期態(tài)進(jìn)入混沌狀態(tài)，同時(shí)當(dāng)參數(shù)b在[8.53，8.59]和[10.21，10.25]處時(shí)出現(xiàn)兩個(gè)較為明顯的周期窗口，這與李雅普諾夫指數(shù)譜是相對(duì)應(yīng)的。

(a) 李雅普諾夫指數(shù)譜

如表2所示，系統(tǒng)參數(shù)b∈(7，11)時(shí)，系統(tǒng)出現(xiàn)了各種不同的動(dòng)力學(xué)行為，有多種類型的周期態(tài)和混沌態(tài)出現(xiàn)，說(shuō)明參數(shù)b的變化對(duì)系統(tǒng)動(dòng)力學(xué)特性有較大的影響。

表2 參數(shù)b變化的系統(tǒng)狀態(tài)Tab.2 The system status with b

2.3.3 階數(shù)q對(duì)系統(tǒng)動(dòng)力學(xué)特性的研究

取a=4、b=9和c=5，步長(zhǎng)h=0.01，系統(tǒng)初值x0=[1，1，1，1]，q∈[0.5，1]。系統(tǒng)隨階數(shù)q變化的指數(shù)譜和分岔圖如圖6所示。從圖6(a)李雅普諾夫指數(shù)譜中可見(jiàn)，當(dāng)q∈[0.53，1]，系統(tǒng)處于混沌狀態(tài)。當(dāng)q<0.52時(shí)，在分岔圖和李雅普諾夫指數(shù)譜中均沒(méi)有數(shù)值，系統(tǒng)處于發(fā)散狀態(tài)。

(a) q=0.53

(a) SE復(fù)雜度

(a) 李雅普諾夫指數(shù)譜

為了進(jìn)一步分析系統(tǒng)適用于保密通信等領(lǐng)域的參數(shù)選擇范圍，更好的觀測(cè)混沌系統(tǒng)的動(dòng)力學(xué)特性[24-27]。取系統(tǒng)階數(shù)q∈(0.5，1)，序列長(zhǎng)度N=50 000，此時(shí)系統(tǒng)SE復(fù)雜度和C0復(fù)雜度如圖7所示。由圖(7)知，當(dāng)q大于0.52時(shí)，SE復(fù)雜度和C0復(fù)雜度的變化趨勢(shì)是逐漸減小的，當(dāng)q=1時(shí)，復(fù)雜度的值達(dá)到最小?？梢?jiàn)，分?jǐn)?shù)階混沌系統(tǒng)是比整數(shù)階混沌系統(tǒng)更為復(fù)雜，隨著階數(shù)減小，系統(tǒng)復(fù)雜度會(huì)相應(yīng)變大，說(shuō)明該系統(tǒng)的分?jǐn)?shù)階狀態(tài)比整數(shù)階狀態(tài)具有高的應(yīng)用價(jià)值。

由圖6可知q<0.52時(shí)，系統(tǒng)是發(fā)散的，為了驗(yàn)證系統(tǒng)產(chǎn)生混沌的最小階數(shù)，引入了0-1測(cè)試。0-1測(cè)試的結(jié)果如圖8所示，當(dāng)參數(shù)q=0.53時(shí)，測(cè)試結(jié)果為類布朗運(yùn)動(dòng)，說(shuō)明系統(tǒng)是混沌的。如圖8(b)所示當(dāng)參數(shù)q=0.52時(shí)，測(cè)試結(jié)果為穩(wěn)定點(diǎn)，說(shuō)明系統(tǒng)是非混沌的。由此可得系統(tǒng)產(chǎn)生混沌的最小階數(shù)0.53×4=2.12。

3 分?jǐn)?shù)階混沌系統(tǒng)DSP實(shí)現(xiàn)

完成分?jǐn)?shù)階混沌系統(tǒng)的硬件實(shí)現(xiàn)是其實(shí)際應(yīng)用的基礎(chǔ)。首先對(duì)DSP進(jìn)行初始化，并進(jìn)行系統(tǒng)初始設(shè)置。然后對(duì)系統(tǒng)式(8)進(jìn)行迭代，將迭代結(jié)果數(shù)據(jù)推送入硬件系統(tǒng)后處理并將結(jié)果從DSP中輸出，最后使用數(shù)據(jù)進(jìn)行初始值的替換并進(jìn)行迭代，直到結(jié)束為止。在數(shù)據(jù)處理環(huán)節(jié)，DSP產(chǎn)生的數(shù)據(jù)被傳送入DAC8552中轉(zhuǎn)換成模擬信號(hào)傳送給示波器。

DSP硬件平臺(tái)中采用的DSP-芯片為T(mén)MS320F28335，采用16位雙通道D/A轉(zhuǎn)換器DAC8552對(duì)由DSP產(chǎn)生的時(shí)間序列進(jìn)行轉(zhuǎn)換。D/A轉(zhuǎn)換器由DSP通過(guò)SPI(串行外圍接口)控制，把數(shù)字信號(hào)轉(zhuǎn)化為模擬信號(hào)傳送到示波器上顯示。

圖9為通過(guò)DSP平臺(tái)得到的混沌系統(tǒng)信號(hào)相圖。此處系統(tǒng)參數(shù)a=4、b=9和c=5，系統(tǒng)階數(shù)q=0.8，步長(zhǎng)h=0.01，將其與計(jì)算機(jī)仿真得到的系統(tǒng)相圖進(jìn)行比較，二者結(jié)果是一致的。

圖9 DSP實(shí)現(xiàn)系統(tǒng)相圖Fig.9 The system phase diagram of DSP implementation

4 結(jié) 論

設(shè)計(jì)了一種四維分?jǐn)?shù)階系統(tǒng)，采用CADM分解算法求解了該四維分?jǐn)?shù)階混沌系統(tǒng)的數(shù)值解并分析了該系統(tǒng)動(dòng)力學(xué)特性。結(jié)果表明，該分?jǐn)?shù)階混沌系統(tǒng)具有復(fù)雜的動(dòng)力學(xué)行為。系統(tǒng)產(chǎn)生混沌的最小階數(shù)為2.12，而且通過(guò)復(fù)雜度分析可知，當(dāng)q=0.53時(shí)，系統(tǒng)復(fù)雜度最高，此時(shí)混沌序列隨機(jī)性最好，安全性能最高。最后在DSP平臺(tái)上完成了該系統(tǒng)的硬件實(shí)現(xiàn)，結(jié)果體現(xiàn)了CADM算法的正確性以及分?jǐn)?shù)階混沌系統(tǒng)的物理可實(shí)現(xiàn)性。研究結(jié)果為分?jǐn)?shù)階混沌系統(tǒng)運(yùn)用于保密通信等領(lǐng)域，提供了理論和實(shí)際應(yīng)用基礎(chǔ)。