張宇琪，孫立山

(哈爾濱工業(yè)大學(xué)1.控制科學(xué)與工程系，2.電氣工程系，黑龍江哈爾濱 150001)

1971年，華裔科學(xué)家蔡少棠提出了憶阻器的概念［1］。2008年HP公司發(fā)現(xiàn)了一種用TiO2摻雜后得到的納米級(jí)電子元件，而這種元件的性質(zhì)恰好符合憶阻器的定義，于是制成了這種憶阻器［2］。

HP憶阻器發(fā)現(xiàn)后不久，美籍學(xué)者Itoh和蔡少棠提出了一系列憶阻器電路用來(lái)替換蔡氏二極管，得到了憶阻器混沌電路［3］。不過(guò)，人們重點(diǎn)研究的是一種光滑連續(xù)的三次函數(shù)磁控型憶阻器，物理上并未實(shí)現(xiàn)。雖然憶阻元件尚處在探索階段，但HP憶阻器已經(jīng)物理實(shí)現(xiàn)，只是尚未走進(jìn)普通實(shí)驗(yàn)室。因此，使用仿真軟件進(jìn)行仿真是有意義的。本文根據(jù)HP實(shí)驗(yàn)室提供的實(shí)物測(cè)試數(shù)據(jù)及相關(guān)理論，以及文獻(xiàn)［4］給出的改進(jìn)數(shù)學(xué)模型，將HP憶阻器應(yīng)用于蔡氏電路，得到了基于憶阻器的混沌電路。我們通過(guò)調(diào)整參數(shù)，并用Matlab和Orcad進(jìn)行了仿真，得到了不同于以往的混沌圖形。

1 HP憶阻器概述

蔡少棠給出的荷控型憶阻元件的數(shù)學(xué)模型為

式中，q為電荷量，φ為磁通量，M(q)為憶阻器阻值。

HP實(shí)驗(yàn)室給出的憶阻元件的基本模型如圖1所示。其中D為憶阻器的長(zhǎng)度，w表示元件的摻雜寬度，總阻值R為摻雜部分阻值與非摻雜部分阻值之和。w=0和w=D時(shí)憶阻器的極限阻值分別用Roff和Ron表示。摻雜寬度w隨外電場(chǎng)改變，u＞0時(shí)有向右的電流通過(guò)元件內(nèi)部，w增大，R變小;u＜0時(shí)有向左的電流通過(guò)元件內(nèi)部，w減小，R變大。

圖1 憶阻元件的基本模型

據(jù)此，HP實(shí)驗(yàn)室給出的流控型憶阻元件的數(shù)學(xué)模型為

其中，μv是表示離子在均勻場(chǎng)中移動(dòng)速率的常數(shù)，約為 10-14m2s-1V-1，通常 Roff/Ron=102-103。

當(dāng)Roff?Ron時(shí)，可近似用式(5)計(jì)算憶阻器阻值。正如文獻(xiàn)［2］所說(shuō)，人們?cè)趯ふ覒涀杵鞯臅r(shí)候總是想著磁場(chǎng)對(duì)于電阻的影響，而實(shí)際上，只要M(q)隨q改變而改變即可。HP憶阻器的特征為:磁通量φ并沒(méi)有起明確的作用;摻雜寬度w是有界的，在0至D之間。以往的研究都沒(méi)有同時(shí)具備這兩個(gè)特征，因此并非典型的HP憶阻器。

憶阻器中的狀態(tài)變量選擇值得思考。按照文獻(xiàn)［3］，蔡少棠將電荷q或磁通φ視為狀態(tài)變量，列出了狀態(tài)方程。而HP實(shí)驗(yàn)室認(rèn)為在不考慮邊界條件的情況下，狀態(tài)變量為摻雜寬度w(t)，正比于q(t)。由電路知識(shí)得，控制量為狀態(tài)變量，因此將狀態(tài)變量取為r=w/D，這與取為w是一樣的。則可用下面各式表示HP憶阻器的特性:

其中，Rinit為初始時(shí)憶阻器的電阻值。由HP實(shí)驗(yàn)室試驗(yàn)數(shù)據(jù)粗略計(jì)算，可取長(zhǎng)度為10nm的憶阻器，其k=104s－1A－1，Roff=40kΩ，Ron=0.1kΩ。

文獻(xiàn)［2］還提到了邊界情況，即當(dāng)完全摻雜或完全未摻雜時(shí)的情況。文中給出了一種通俗的假設(shè):當(dāng)上述兩種情況發(fā)生時(shí)，電阻保持不變，直到離開邊界范圍。在納米級(jí)范圍，很小的電壓能產(chǎn)生巨大的電場(chǎng)，進(jìn)而在離子移動(dòng)過(guò)程中產(chǎn)生重要的非線性。這種非線性在邊界處表現(xiàn)得更加明顯，使得摻雜與非摻雜分界線移動(dòng)速度逐漸趨于零，這一現(xiàn)象稱為摻雜漂移，可以用窗函數(shù)f(r)表示。雖然能用來(lái)解釋HP憶阻器的確切的窗函數(shù)f(r)并沒(méi)找到，但已有一些模型供可我們使用。HP實(shí)驗(yàn)室在研究時(shí)給出了一種假設(shè):w(1－w)/D2。依據(jù)這種假設(shè)，給出了仿真的u-i關(guān)系圖，并與實(shí)際測(cè)量的關(guān)系圖進(jìn)行比較，兩者十分接近。但文獻(xiàn)［2］中同時(shí)說(shuō)，不考慮非線性漂移的情況也有報(bào)道。因此，文中指出，可以根據(jù)豐富的u-i關(guān)系圖來(lái)分析邊界問(wèn)題。

圖2 HP實(shí)驗(yàn)室滯后環(huán)的u-i關(guān)系

在文獻(xiàn)［4］中給出了f(r)的另一種形式:

其中，p是一個(gè)正整數(shù)。若p較小，則分界線快要達(dá)到邊界時(shí)速度逐漸減小到0;若p越大，這種控制速度的功能就越不明顯;p趨于無(wú)窮，則起不到邊界限制的作用，本文采用了這種模型。

事實(shí)上，非線性漂移很復(fù)雜，后續(xù)的學(xué)術(shù)報(bào)道中指出，即使經(jīng)過(guò)同一位置，分界線朝兩個(gè)方向的速度并不相同，且移動(dòng)速度還和電場(chǎng)強(qiáng)度有關(guān)。因此，上述模型都是簡(jiǎn)單模型，不可能完全表征HP憶阻器，但這對(duì)我們仿真分析已經(jīng)足夠了。

2 憶阻混沌電路的Matlab仿真

我們可用HP憶阻器模型與兩個(gè)電感、一個(gè)電容和一個(gè)負(fù)電阻構(gòu)成一個(gè)混沌電路。如圖3所示。

圖3 標(biāo)準(zhǔn)蔡氏憶阻器電路

列出圖示電路的狀態(tài)方程，共有四個(gè)狀態(tài)變量:

其中，M(r)=rRon+(1－r)Roff，窗函數(shù) f(r)=1－(2r－1)2p。

我們做變換使 x=i1，y=u，z=i2，e=r，令 a=L2/L1，b=1/C，c=1/L2，R/L2=1，k=2.5，M(e)=M(r)=0.025e+10(1－e)，可得到其無(wú)量綱方程:

取 a=2.7，b=3.6，c=0.79，初始條件為 x=0，y=0.05，z=0，e=0.3。代入 Matlab 中進(jìn)行仿真，得到的相空間或相平面上的投影如圖4所示。

圖4 混沌吸引子的投影

因方程為剛性的，最好用ode15s或ode23s而不宜用ode45算法，否則所用時(shí)間較長(zhǎng)。Lyapunov指數(shù)分別為 L1=0.172，L2=0.001，L3=-0.005，L4=12.46，Lyapunov維數(shù) dL=3.01。狀態(tài)變量 y的時(shí)域波形如圖5(a)所示，其特性貌似隨機(jī)非周期性。在y=0的截面上三維Poincaré映射軌線在z-e平面上的投影如圖5(b)所示。從相軌圖、Lyapunov指數(shù)和維數(shù)、Poincaré映像以及時(shí)域波形可以得出該電路是混沌電路。

圖5 狀態(tài)變量y的時(shí)域波形和Poincaré映射

以b為變量的分岔圖及Lyapunov指數(shù)譜如圖6所示。此處畫分岔圖所選取的Poincaré截面y=0是一個(gè)在四維相空間中垂直于y軸過(guò)原點(diǎn)的平面，Lyapunov指數(shù)譜省略了L4。由圖可知，分岔圖和Lyapunov指數(shù)譜是基本一致的:當(dāng)b＜4.7時(shí)，系統(tǒng)混沌;當(dāng)b＞6.7時(shí)，系統(tǒng)收斂到一固定點(diǎn);4.7＜b＜6.7時(shí)，最大Lyapunov指數(shù)時(shí)正時(shí)負(fù)，動(dòng)力學(xué)行為比較復(fù)雜，在此不做研究。

圖6 參數(shù)b變化的分岔圖和Lyapunov指數(shù)譜

該混沌系統(tǒng)的平衡點(diǎn)集為 A= {(x，y，z，e)|x=y=z=0，e=const}，即e坐標(biāo)上的點(diǎn)集均是平衡點(diǎn)。由e的定義知，這里的c是一個(gè)0到1的常數(shù)，這與一般的混沌系統(tǒng)是不同的。當(dāng)e=0.3時(shí)，有λ1，2=0.2379 ± j1.5297λ3=-18.3961，λ4=0，平衡點(diǎn)是不穩(wěn)定的焦點(diǎn);當(dāng) e=0.9 時(shí)，λ1，2=-0.9723 ±j3.0357，λ3=0.1839，λ4=0，平衡點(diǎn)是不穩(wěn)定的鞍點(diǎn);當(dāng) e=0.7 時(shí)，λ1，2=-0.1794 ± j1.3961，λ3=-6.7885，λ4=0，平衡點(diǎn)是穩(wěn)定點(diǎn)。

所以，不同位置的平衡點(diǎn)有不同的穩(wěn)定性。因此e在不同的初始值e時(shí)，其運(yùn)行軌道會(huì)在混沌行為、周期行為或穩(wěn)定的匯之間發(fā)生狀態(tài)轉(zhuǎn)移，在一定的電路參數(shù)下零特征值(λ4=0)對(duì)憶阻器混沌電路的穩(wěn)定性也有很大影響。

3 憶阻混沌電路的Orcad仿真

憶阻器與二極管等非線性元件存在差異，憶阻器有記憶功能，而其他元件無(wú)此功能。例如，可通過(guò)伏安法測(cè)得二極管的伏安特性曲線，而憶阻器只要加正電壓v+，電阻就會(huì)一直改變，從關(guān)狀態(tài)經(jīng)過(guò)一段時(shí)間約D2Roff/(2μvv+Ron)也會(huì)變?yōu)殚_狀態(tài)。因此，不能用伏安法測(cè)憶阻器的伏安特性，而需要應(yīng)用正弦電壓掃描。

根據(jù)HP憶阻器的非線性摻雜漂移模型建立一個(gè)器件模型，可以在Orcad中像電阻一樣使用。圖7為用正弦交流電源供電時(shí)的滯后環(huán)圖像，與HP實(shí)驗(yàn)室的結(jié)果相差不大，差異主要在于窗函數(shù)的構(gòu)建。由圖可見，當(dāng)交流電壓幅值為2V時(shí)，邊界效應(yīng)小，如圖7(a)所示;當(dāng)幅值為4V時(shí)，邊界效應(yīng)大，如圖7(b)所示。

圖7 憶阻器Orcad仿真結(jié)果

目前，憶阻器的Roff單位均為千歐級(jí)，而頻率超過(guò)10Hz就不會(huì)有滯后環(huán)。圖中Ron=0.1kΩ，Roff=40kΩ ，Rinit=28kΩ，p=1，f=1Hz，k=104s-1A-1。由于電路方程為剛性的，在Orcad仿真時(shí)要設(shè)置最大步長(zhǎng)TMAX為1ms，相對(duì)精度Reltol為1E-6，否則得不到結(jié)果，這也說(shuō)明了Orcad的數(shù)值仿真能力略遜于Matlab。

使用Orcad進(jìn)行仿真，電感可以用回轉(zhuǎn)器來(lái)等效替代，實(shí)現(xiàn)無(wú)內(nèi)阻的電感，負(fù)電阻可用負(fù)電阻變換電路來(lái)實(shí)現(xiàn)。仿真電路如圖8所示，此處用到的憶阻器中，Ron=0.025Ω，Roff=10Ω，Rinit=7Ω。仿真結(jié)果如9所示，其結(jié)果與Matlab仿真結(jié)果相同。

4 結(jié)語(yǔ)

圖8 Orcad仿真電路圖

圖9 憶阻器混沌電路的Orcad仿真結(jié)果

本文介紹了HP憶阻器非線性摻雜漂移特性，說(shuō)明了HP憶阻器的獨(dú)特性質(zhì)。將HP憶阻器與蔡氏電路結(jié)合起來(lái)，構(gòu)成了基于HP憶阻器的混沌電路。對(duì)憶阻器的Orcad模型進(jìn)行了特性曲線的測(cè)試，與文獻(xiàn)［2］中結(jié)果相同?；煦珉娐返腛rcad仿真結(jié)果與Matlab結(jié)果相同。Matlab可以進(jìn)行大量的數(shù)值計(jì)算，如計(jì)算Lyapunov指數(shù)等，而Orcad則可以設(shè)計(jì)電路元件參數(shù)。因此，要在不同場(chǎng)合使用不同軟件，發(fā)揮其優(yōu)勢(shì)，使二者相輔相成。

［1］ L.O.Chua.Memristor-the missing circuit element［J］.IEEE Transaction on Circuit Theory，1971，18(5):507-519

［2］ D.B .Strukov，G.S.Snider，D.R.Stewart and R.S.Williams.The missing memristor found［J］.Naturre，2008，453:80-83

［3］ M.ltoh，L.O.Chua.Memristor Oscillators［J］.International Journal of Bifurcation and Chaos，2008，18(11):3183-3206

［4］ Y.N.JOGLEKAR，S.J.WOLF，The elusive memristor:properties of basic electrical circuits.European Journal of Physics，2009.30:661-675