離散憶阻混沌系統(tǒng)的Simulink 建模及其動力學特性分析*

扶龍香 賀少波 王會海 孫克輝

(中南大學物理與電子學院,長沙 410083)

為拓廣離散憶阻器的研究與應用,基于差分算子,構(gòu)建了具有平方非線性的離散憶阻模型,并實現(xiàn)了Simulink仿真.仿真結(jié)果表明,設(shè)計的憶阻器滿足廣義憶阻定義.將得到的離散憶阻引入三維混沌映射中,設(shè)計了一種新型四維憶阻混沌映射,并建立了該混沌映射的Simulink 模型.通過平衡點、分岔圖、Lyapunov 指數(shù)譜、復雜度、多穩(wěn)態(tài)分析了系統(tǒng)復雜動力學特性.本文從系統(tǒng)建模角度出發(fā),構(gòu)建離散憶阻與離散憶阻混沌映射,進一步驗證了離散憶阻的可實現(xiàn)性,為離散憶阻應用研究奠定了基礎(chǔ).

1 引言

1971 年,蔡少棠教授[1]在研究基本電學物理量即電壓V、電流I、電荷q、磁通φ之間的相互關(guān)系時,基于電路完備性考慮,預測存在除電阻、電感和電容外的第4 種無源非線性基本電路元件,用以描述電荷與磁通之間的關(guān)系,稱其為憶阻器.2008 年,惠普(Hewlett-Packard,HP)實驗室[2]利用TiO2和金屬鉑首次實現(xiàn)了世界上納米級憶阻器,從此掀起了憶阻器研究熱潮.近年來,憶阻器在許多應用領(lǐng)域都有深入研究,如圖像處理[3]、存儲器[4-6]、神經(jīng)網(wǎng)絡[7-9]、神經(jīng)動力學[10,11]等.但由于HP 憶阻器納米級工藝以及嚴格實驗條件的限制,迄今為止,憶阻器尚未得到商用.為了推進對憶阻器等記憶元件的研究,建立合適的仿真模型至關(guān)重要.近年來,很多基于SPICE (simulation pro gram with integrated circuit emphasis)的各類憶阻仿真模型相繼被提出[12-14].同時,Simulink 作為Matlab 中一種可視化建模工具,能對大多數(shù)非線性系統(tǒng)進行分析,具有流程清晰、易觀察、效率高等諸多優(yōu)點,便于更好地研究與分析非線性系統(tǒng).因此,采用Matlab/Simulink 仿真工具對憶阻器進行建模值得進一步研究[15-17].

憶阻器是一種非線性元件,其在混沌系統(tǒng)方面的應用也是憶阻器的研究熱點之一.在王春華等[18]對憶阻器相關(guān)應用總結(jié)中提到,可將憶阻器應用到混沌電路分為4 種方式,即用憶阻器替換蔡氏電路中非線性元件[19-21]、新增憶阻器[22-24]、憶阻器替換非Chua 系統(tǒng)中線性/非線性項的方法[25,26]、利用不同數(shù)學模型的憶阻器進行電路設(shè)計[27,28].本文將利用新增憶阻的方式,將離散憶阻引入到三維混沌映射中,構(gòu)建一種新型四維憶阻混沌映射,并研究離散憶阻混沌系統(tǒng)建模方法.

本文內(nèi)容安排如下.第2 節(jié),基于離散憶阻數(shù)學模型建立離散憶阻Simulink 仿真模型;第3 節(jié),將離散憶阻模型加入到三維Lorenz 混沌映射中,設(shè)計一新型混沌映射,并構(gòu)建其Simulink 模型;第4 節(jié),對該離散憶阻器及混沌系統(tǒng)Simulink 模型的應用進行分析,最后,對本文進行總結(jié).

2 離散憶阻建模與Simulink 仿真

2.1 離散憶阻數(shù)學模型

類似于連續(xù)憶阻器模型,離散憶阻器的定義為[29]

其中,n=0,1,2,···,Δq(tn)=q(tn+1)-q(tn)為前向差分.因為

等號兩邊相加,則憶阻“內(nèi)部電荷”計算式為

顯然,q(t0)為離散憶阻內(nèi)部初始電荷.在實際計算過程中,該離散憶阻器為一離散累加數(shù)學模型.一般地,可以將該離散憶阻改寫為迭代式,即:

式中,i(tn)為離散輸入電流信號,V(tn)為離散輸出電壓信號,qn為憶阻內(nèi)部電荷,M(qn)為憶阻值.其中M(·)函數(shù)將采用平方函數(shù)[30],因此,該離散憶阻器模型可描述為

式中,α和β為憶阻器參數(shù),令k=1.

2.2 Simulink 仿真實現(xiàn)

根據(jù)離散憶阻器的數(shù)學模型,構(gòu)建離散憶阻器的Simulink 模型,設(shè)計的模型如圖1 所示.該模型將離散憶阻數(shù)學模型中的加減乘除分類整合,且由離散時間求和模塊結(jié)合前向差分算子實現(xiàn)電荷q的迭代,其中可通過改變Constant 與Gain1 的取值來改變α與β的值,Scope1 顯示的是憶阻器輸入電壓vn與流經(jīng)其電流in之間的關(guān)系.首先,對該離散憶阻Simulink 模型進行驗證,這里取α=1,β=0.000002.發(fā)現(xiàn)在正弦電流的作用下,憶阻器電流-電壓特性曲線見圖2(a),具有遲滯環(huán)特性.其次,分別改變輸入正弦信號I(tn)=A0sin(2πwtn)的幅值和頻率以討論輸入信號幅值和頻率對該憶阻模型的影響.取不同的幅值A(chǔ)0=1,1.5,2 A,其中ω=0.005 rad/s,結(jié)果如圖2(b)所示.可見,在輸入信號頻率保持不變的情況下,該憶阻模型的滯回曲線隨著輸入信號幅值的增大而增大.保持輸入信號幅值A(chǔ)0=1.5 A 不變,改變輸入信號頻率,即ω分別取0.1,0.01,0.003 rad/s,結(jié)果如圖2(c)所示,可見在輸入信號幅值一定的情況下,該憶阻模型的滯回曲線隨著輸入信號頻率增大而減小,最終退化為線性器件.經(jīng)以上分析,根據(jù)憶阻器特性判定原則可知,所構(gòu)建的離散憶阻Simulink 模型滿足憶阻定義[31],可用于進一步的應用研究.

圖1 離散憶阻Simulink 模型Fig.1.Discrete memristor Simulink model.

圖2 離散憶阻器的電流-電壓特性曲線 (a)電流-電壓特性曲線;(b)幅值變化;(c)頻率變化Fig.2.Current-voltage characteristic curves of discrete memristor:(a)Current-voltage characteristic curve;(b)amplitude change;(c)frequency change.

3 離散憶阻混沌系統(tǒng)及其Simulink仿真

3.1 離散憶阻混沌系統(tǒng)

Sprott[32]報道了一種三維混沌映射,并稱之為Lorenz 混沌映射,其動力學方程為

采用結(jié)構(gòu)化設(shè)計方法,該映射可以被設(shè)計為圖3 所示的結(jié)構(gòu)圖.可見,該系統(tǒng)存在3 個延時模塊,1 個加法器和1 個乘法器.

圖3 三維Lorenz 混沌映射(6)系統(tǒng)框圖Fig.3.Structure diagram of the three-dimensional Lorenz chaotic map (6).

引入離散憶阻,得到一新型三維憶阻混沌映射,其定義為

其中,a為系統(tǒng)參數(shù),用于控制憶阻的輸出值對系統(tǒng)動力學行為的影響.根據(jù)上述系統(tǒng)方程(7),繪制引入憶阻器后的混沌系統(tǒng)框圖,如圖4 所示.可見,離散憶阻增加在xn與yn的乘積項上.

圖4 加入憶阻器后的三維混沌映射(7)系統(tǒng)框圖Fig.4.Structure diagram of the three-dimensional chaotic map with memristor (7).

結(jié)合離散憶阻數(shù)學模型,令:

進而可以得到系統(tǒng)的四維表達式為

根據(jù)(9)式繪制設(shè)計的憶阻混沌系統(tǒng)框圖,如圖5 所示,其中陰影部分描述的是離散憶阻的結(jié)構(gòu).

圖5 四維憶阻混沌映射(9)系統(tǒng)框圖Fig.5.Structure diagram of the four-dimensional memristor chaotic map (9).

3.2 離散憶阻混沌系統(tǒng)Simulink 仿真

根據(jù)四維憶阻混沌系統(tǒng),建立該系統(tǒng)Simulink仿真模型,如圖6 所示,其中憶阻模塊采用圖1 所示的模型.模型中增益Gain2 表示系統(tǒng)參數(shù)a的值,Scope2 顯示的是系統(tǒng)xn,yn,zn時間序列,并可改變3 個延時器以及憶阻模塊中離散時間求和模塊中的初值,從而分別改變系統(tǒng)初始狀態(tài)x0,y0,z0,w0值.

圖6 四維憶阻混沌映射(9)Simulink 模型Fig.6.Simulink model of thefour-dimensional memristor chaotic map (9).

設(shè)置該模型系統(tǒng)參數(shù)a=0.9,初始狀態(tài)為x0=0.5,y0=0.5,z0=0.5,w0=0.4,憶阻參數(shù)α=1,β取值分別為—0.1,—0.02,—0.000002.將示波器中的圖像保存到WorkSpace,再利用Matlab 繪制出圖像,得到x-y平面的吸引子相圖,如圖7(a)—(c)所示.可見,隨著β的取值的不同,吸引子大小形狀產(chǎn)生較大變化.設(shè)置系統(tǒng)初始狀態(tài)為x0=0.5,y0=0.5,z0=0.5,w0=0.4,憶阻參數(shù)α=1,β=—0.001,改變系統(tǒng)參數(shù)a,分別取值為0.25,0.5,0.9,得到的xn-yn平面吸引子相圖,如圖7(d)—(f)所示,表明參數(shù)a也能夠改變吸引子大小與形狀.繪制圖7 中各個狀態(tài)下對應的Lyapunov 指數(shù),如圖8所示.可見,系統(tǒng)參數(shù)以及憶阻參數(shù)的變動能夠改變系統(tǒng)的動力學行為.

圖7 四維憶阻混沌映射(9)的Simulink 仿真吸引子圖 (a)β=—0.1,超混沌狀態(tài);(b)β=—0.02,超混沌狀態(tài);(c)β=—0.000002,混沌狀態(tài);(d)a=0.25,混沌狀態(tài);(e)a=0.5,混沌狀態(tài);(f)a=0.9,混沌狀態(tài)Fig.7.Simulink simulation results of the four-dimensional memristor chaotic map (9):(a)β=—0.1,hyperchaotic;(b)β=—0.02,hyperchaotic;(c)β=—0.000002,chaos;(d)a=0.25,chaos;(e)a=0.5,chaos;(f)a=0.9,chaos.

圖8 四維憶阻混沌映射(9)的Simulink 仿真吸引子圖對應Lyapunov 指數(shù)譜 (a)β=—0.1,超混沌狀態(tài);(b)β=—0.02,超混沌狀態(tài);(c)β=—0.000002,混沌狀態(tài);(d)a=0.25,混沌狀態(tài);(e)a=0.5,混沌狀態(tài);(f)a= 0.9,混沌狀態(tài)Fig.8.Simulink simulation attractor diagram of the four-dimensional memristor chaotic map (9)corresponds to Lyapunov Exponent spectra:(a)β=—0.1,hyperchaotic;(b)β=—0.02,hyperchaotic;(c) β=—0.000002,chaos;(d)a=0.25,chaos;(e)a=0.5,chaos;(f)a=0.9,chaos.

4 動力學特性分析

4.1 平衡點及其穩(wěn)定性分析

令(9)式左右兩邊相等,即可求出系統(tǒng)平衡點為

其中,C可以為任意值,也就是說w(e)可以取任意值,即系統(tǒng)具有無限多的平衡點.

對系統(tǒng)進行線性化分析,得其Jacobian 矩陣為

代入平衡點后矩陣為

經(jīng)計算得其特征值為λi=[1,—1,0.5+0.866i,0.5—0.866i].可見,系統(tǒng)存在不穩(wěn)定平衡點,可產(chǎn)生混沌,即驗證了前述Simulink 仿真結(jié)果.

4.2 分岔圖與Lyapunov 指數(shù)譜分析

系統(tǒng)(9)具有無限多不穩(wěn)定的平衡點,表明此離散憶阻混沌系統(tǒng)具有豐富的動力學行為,為了進一步研究系統(tǒng)的復雜動力學特性,分析了當a在[0,1]范圍變化時,Lyapunov 指數(shù)譜和系統(tǒng)的分岔圖.

Lyapunov 指數(shù)是相空間中鄰近軌跡發(fā)散或收斂的平均指數(shù)率,出現(xiàn)正的Lyapunov 指數(shù)時,表明系統(tǒng)是混沌的.初始狀態(tài)設(shè)置為x0=0.5,y0=0.5,z0=0.5,w0=0.4,憶阻參數(shù)α=1,β=—0.001,迭代次數(shù)為N=10000,通過Matlab 程序仿真繪制xn隨參數(shù)a變化的分岔圖和Lyapunov 指數(shù)譜,如圖9 所示.可見,Lyapunov 指數(shù)譜與分岔圖吻合,兩圖都表明了系統(tǒng)能產(chǎn)生復雜的動力學行為.當a∈[0,0.26)時,Lyapunov 指數(shù)有兩個大于零的值,但其值較小,且分岔圖形狀與a∈[0.26,1]時有明顯差別.因此,進一步分析a的取值與該區(qū)間xn序列值及對應的Lyapunov 指數(shù),即a=0.1,憶阻參數(shù)和初值保持不變,繪制xn序列值及對應Lyapunov 指數(shù)如圖10 所示.通過xn序列發(fā)現(xiàn)此時系統(tǒng)處于非周期狀態(tài),存在混沌成分,且Lyapunov 指數(shù)值很小,處于弱混沌狀態(tài).當a=0.9 時,Lyapunov 指數(shù)中有兩個正數(shù),且值相對較大,結(jié)合分岔圖可以發(fā)現(xiàn)該系統(tǒng)處于超混沌狀態(tài).通過分岔圖還可以發(fā)現(xiàn)在超混沌區(qū)間中出現(xiàn)了周期窗口,表明系統(tǒng)的動力學行為對參數(shù)a的取值比較敏感.

圖9 四維憶阻混沌映射(9)分岔圖與Lyapunov 指數(shù)譜 (a)分岔圖;(b)Lyapunov 指數(shù)譜Fig.9.Bifurcation and Lyapunov exponent spectra of the four-dimensional memristor chaotic map (9):(a)Bifurcation diagram;(b)Lyapunov exponent (LE)spectra.

圖10 弱混沌系統(tǒng)動力學分析 (a)xn 序列圖;(b)Lyapunov 指數(shù)譜Fig.10.Dynamic analysis of weakly chaotic system:(a)xn sequence diagram;(b)Lyapunov exponent spectra.

初始狀態(tài)設(shè)置為x0=0.5,y0=0.5,z0=0.5,w0=0.4,N=10000,憶阻參數(shù)α=1,β=—0.001,系統(tǒng)參數(shù)a分別取0.25,0.5,繪制四維憶阻混沌映射(9)隨初值z0變化的Lyapunov 指數(shù)譜,結(jié)果如圖11(a)和圖11(b)所示,原三維Lorenz 混沌映射(6)的Lyapunov 指數(shù)譜,如圖11(c)所示,對比可見,加入憶阻器后系統(tǒng)動力學行為變得更加復雜了,由混沌狀態(tài)進入了超混沌狀態(tài).

圖11 隨初值z0 變化的Lyapunov 指數(shù)譜 (a)四維憶阻混沌映射(9),a=0.25;(b)四維憶阻混沌映射(9),a=0.5;(c)三維Lorenz 混沌映射(6)Fig.11.Lyapunov exponent spectra with initial value z0:(a)Four-dimensional memristor chaotic map (9),a=0.5;(b)four-dimensional memristor chaotic map (9),a=0.5;(c)three-dimensional Lorenz chaotic map (6).

4.3 復雜度分析

復雜度是衡量一個混沌序列與隨機序列接近程度的指標,復雜度越大意味著混沌序列越隨機[33].繪制當初始狀態(tài)設(shè)置為x0=0.5,y0=0.5,z0=0.5,w0=0.4,憶阻參數(shù)α=1,β=—0.001,迭代次數(shù)N=10000,a在[0,1]之間的系統(tǒng)樣本熵復雜度[34],如圖12 所示.可見系統(tǒng)樣本熵復雜度呈上升趨勢,且在a=0.26 后,系統(tǒng)進入高復雜度狀態(tài),系統(tǒng)中的混沌序列越來越接近隨機序列,再次表明a值約為0.26,系統(tǒng)處于混沌狀態(tài).同時,可以發(fā)現(xiàn)該混沌系統(tǒng)高復雜度區(qū)域較寬,且與對應的分岔圖和Lyapunov 指數(shù)譜匹配.

圖12 四維憶阻混沌映射(9)樣本熵復雜度Fig.12.Sample entropy (SampEn)complexity of the fourdimensional memristor chaotic map (9).

4.4 多穩(wěn)態(tài)分析

設(shè)系統(tǒng)初始狀態(tài)為x0=0.5,z0=0.5,w0=0.4,N=10000,系統(tǒng)參數(shù)a=0.9,憶阻參數(shù)α=1,β=—0.001,憶阻初值y0分別取值為1,0.5,0.25,繪制xn-yn平面的吸引子共存圖,如圖13(a)所示.取系統(tǒng)參數(shù)a=0.9,初始狀態(tài)設(shè)置為x0=0.5,y0=0.5,w0=0.4,N=10000,憶阻參數(shù)α=1,β=—0.001.憶阻初值z0分別取值為1,0.5,0.25,繪制xn-yn平面的吸引子共存圖,如圖13(b)所示.設(shè)置系統(tǒng)初始狀態(tài)為x0=0.5,y0=0.5,z0=0.5,N=10000,系統(tǒng)參數(shù)a=0.9,憶阻參數(shù)α=1,β=—0.001,憶阻初值w0分別取值為1,0.5,0.25,繪制x-y平面的吸引子共存圖,如圖13(c)所示.可見設(shè)計的離散憶阻系統(tǒng)存在吸引子共存現(xiàn)象.設(shè)系統(tǒng)參數(shù)a=0.9,圖14 分別表示隨初值y0,z0,w0變化的分岔圖.可見,隨著系統(tǒng)初值的變化,各分岔圖都出現(xiàn)了較多的周期窗口且吸引子大小也在不斷變化,進一步表明系統(tǒng)存在無限多共存吸引子.

圖13 xn-yn 平面上共存吸引子 (a)初值y0=1,0.5,0.25;(b)初值z0=1,0.5,0.25;(c)初值w0=1,0.5,0.25Fig.13.Coexisting attractors in the xn-yn plane:(a)y0=1,0.5,0.25;(b)z0=1,0.5,0.25;(c)w0=1,0.5,0.25.

圖14 隨初值變化的分岔圖 (a)初值y0 變化,x0=0.5,z0=0.5,w0=0.4;(b)初值z0 變化,x0=0.5,y0=0.5,w0=0.4;(c)初值w0 變化,x0=0.5,y0=0.5,z0=0.5Fig.14.Bifurcation diagrams with the initial values:(a)Initial value y0 changes,x0=0.5,z0=0.5,w0=0.4;(b)initial value z0 changes,x0=0.5,y0=0.5,w0=0.4;(c)initial value w0 changes,x0=0.5,y0=0.5,z0=0.5.

對于不同系統(tǒng)初值條件下,共存的無限多吸引子的吸引盆可以描繪吸引子與初始條件的關(guān)系[35].系統(tǒng)參數(shù)分別為a=0.25,0.5,0.9,初始值為x0=0.5,y0=0.5,N=10000,憶阻參數(shù)α=1,β=—0.001,z0-w0初值平面上的吸引盆如圖15,圖中不同顏色代表不同的吸引子類型.可見,總體上,a較大時吸引子類型變化更加復雜.根據(jù)圖15 中不同a值吸引盆,分析對應的隨初值z0變化的樣本熵復雜度如圖16 所示,發(fā)現(xiàn)a=0.9 時,復雜度值總體上明顯高于a=0.25 和a=0.5 時的值,且對初值的變化更加敏感,進一步驗證了圖15 中的結(jié)果.

圖15 初值平面上吸引盆 (a)a =0.25;(b)a =0.5;(c)a =0.9Fig.15.Attraction basins in the planes constructed by two different initial conditions:(a)a =0.25;(b)a =0.5;(c)a =0.9.

圖16 在不同參數(shù)a 下隨初值z0 變化的樣本熵復雜度 (a)a =0.25;(b)a =0.5;(c)a =0.9Fig.16.Sample entropy complexity with initial value z0 under different system parameter a:(a)a =0.25;(b)a =0.5;(c)a =0.9.

原系統(tǒng)隨初值z0變化的復雜度如圖17 所示,對比圖16,再次驗證加入憶阻器后系統(tǒng)復雜度變大,且高復雜度的區(qū)域更寬,表明設(shè)計的憶阻混沌系統(tǒng)具有較原系統(tǒng)更高的復雜度.

圖17 三維Lorenz 混沌映射(6)樣本熵復雜度Fig.17.Sample entropy complexity of the three-dimensional Lorenz chaotic map (6).

5 結(jié)論

本文在離散憶阻數(shù)學模型基礎(chǔ)上構(gòu)建了離散憶阻Simulink 模型.在Simulink 模型中借助離散時間求和模塊實現(xiàn)電荷的迭代,以體現(xiàn)離散憶阻器獨特的記憶特性.仿真結(jié)果發(fā)現(xiàn),設(shè)計的憶阻器滿足憶阻定義,表明本文從仿真的角度實現(xiàn)了憶阻器.同時,將此離散憶阻器引入到三維Lorenz 混沌映射中,實現(xiàn)了一種新型的四維憶阻混沌映射,并建立了該新型四維憶阻混沌映的Simulink 模型.通過分析該離散混沌映射的平衡點及穩(wěn)定性,發(fā)現(xiàn)其具有無窮多平衡點.通過分析該混沌映射分岔圖、Lyapunov 指數(shù)譜、系統(tǒng)復雜度、吸引子共存、吸引盆以及對應初值平面復雜度,研究了其復雜的動力學特性.數(shù)值仿真分析結(jié)果表明新型四維憶阻混沌映具有較原系統(tǒng)更豐富的動力學行為以及更高的復雜度.本文的研究進一步驗證了離散憶阻模型的可實現(xiàn)性及其潛在的應用價值,為離散憶阻在非線性系統(tǒng)中的應用研究奠定了基礎(chǔ).