胡煒 廖建彬 杜永乾
1) (福州大學物理與信息工程學院, 福州 350116)
2) (集美大學輪機工程學院, 福建省船舶與海洋工程重點實驗室, 廈門 361021)
3) (西北工業大學深圳研究院, 深圳 518057)
4) (西北工業大學電子信息學院, 西安 710072)
憶阻網絡是一種基于憶阻器單元的大規模非線性電路, 在下一代人工智能、生物電子、高性能存儲器等新興研究領域發揮著重要作用.描述憶阻器單元物理和電學特性的模型對憶阻網絡的性能仿真具有顯著影響.然而, 現有模型主要為非解析模型, 應用于憶阻網絡分析時可能存在收斂性問題.因此, 提出了一種基于同倫分析法(homotopy analysis method, HAM)的憶阻器單元解析建模策略, 該策略具有解析性和收斂性優化的特點, 可提高憶阻器單元和相應憶阻網絡的收斂性.此外, 還提出了一種面向憶阻器單元模型的驗證準則, 以驗證模型在大規模憶阻網絡中的適用性.通過憶阻器單元和憶阻矩陣網絡的長時演化實驗以及與傳統非解析(數值)方法的比較, 驗證了所提策略的解析性和收斂性優勢; 利用不同類型憶阻器單元和輸入的實驗, 驗證了該策略的擴展性.進一步地, 基于上述實驗, 揭示了憶阻網絡仿真出現收斂性問題的潛在原因.該策略可應用于基于憶阻網絡的新興研究.
憶阻器是一種由多層納米薄膜組成的新型微納器件, 具有高非線性和非易失性開關的特點, 能夠以非易失方式“記憶”流經的電荷和磁通量, 具體體現為內部狀態變量和憶阻值的動態演化[1,2].此外, 憶阻器還具有低開關功耗(約1 fJ)[3]、納米尺寸(約幾 nm)[4,5]、快速讀寫時間(約1 ps)[6]、簡單器件結構(約多層薄膜結構)[7,8]、可與現有CMOS半導體工藝兼容以實現3D集成[9]等優點.因而憶阻器在人工智能[10?13]、非易失性存儲器[14,15]、神經網絡[16?18]等新興研究領域具有廣闊的應用前景.
上述憶阻器應用的一個共同特征是大規模憶阻網絡的使用.在這些網絡中, 大量的憶阻器作為基本單元以特定的電路連接形式構建了復雜的電路拓撲架構, 從而實現了特定的電路功能, 如神經網絡中的多層感知.為了有效分析和評估憶阻網絡的性能, 憶阻器單元的模型至關重要, 其由憶阻器建模策略得來.由于描述憶阻器單元復雜工作機理的物理模型(即各物理方程, 其中以狀態方程最為重要)具有強非線性, 因而不易通過基本函數來解析求解核心變量—狀態變量的動態演化, 故傳統建模策略常采用非解析(數值)的方法以獲得近似的、適用于憶阻網絡的憶阻器單元模型.
如Mladenov[19]通過PSpice工具成功測試了5 × 5交叉桿陣列(crossbar array)和6 × 6憶阻矩陣電路的性能, 上述電路中每個基本單元均采用數值型HfO2憶阻器模型, 該模型中Sinh和Exp等高非線性函數對電路收斂性和仿真時間的影響并未被進一步研究.Biolek等[20]為了提升惠普憶阻器單元模型[1]的魯棒性, 在狀態方程中引入了磁通/電荷與狀態變量之間的非線性轉換函數.通過該函數的作用, 基于惠普模型的憶阻網絡規模可從5 × 5提升到700 × 700, 該文獻并未分析轉換函數對網絡收斂性的作用.Li等[21]驗證了現象學憶阻器單元模型在64 × 64交叉桿陣列結構中的適用性, 該模型采用基于p-Si/SiO2/n-Si憶阻器實測數據的簡化電荷空間型物理機理, 上述交叉桿陣列電路僅適用于特定材料下的多層薄膜憶阻器單元, 存在擴展性受限的問題.Yakopcic等[22]將特定條件下可“通用”的Spice模型嵌入到包含有256個憶阻器單元的憶阻電路, 通過仿真驗證了該憶阻電路的功能, 該模型需要將憶阻器單元原始物理模型與所提出的“通用”公式進行多參數擬合, 擬合方法和擬合的準確性會嚴重影響模型的收斂性.Li等[23]將氧化型憶阻器單元的簡化導電絲模型嵌入到單晶體管-單憶阻器結構中, 并以該結構為基本單元組成了從8 × 8到128 × 128的存儲陣列電路, 以仿真驗證大規模憶阻器陣列的存儲功能.該模型由于其非解析特性, 在仿真過程中易出現數據溢出的現象.
然而, 上述非解析(數值)模型由于無法得到模型的解析解表達式(也稱作顯函數或閉合函數),其求解和仿真僅限定于數值求解和數值仿真, 并具有如下限制: 1)由于解表達式的缺失, 導致數值模型無法評估、預測和優化憶阻器的物理行為; 2)數值模型的解為一系列離散的數值, 無法對所求各物理變量作定性分析, 且不具備符號計算的能力, 即無法由所得變量直接求出其他待求變量的完整表達式; 3)非解析建模所采用的核心數學方法, 如數值積分法和有限差分法, 可能導致各種離散、截斷和舍入誤差[20].盡管這些誤差值量級較小, 但對長時演化下的憶阻網絡仿真具有不可忽視的負面影響.具體而言, 采取上述方法求得的狀態變量數值近似解具有偏離理論準確解(假設該解存在)的潛在風險, 從而可導致網絡仿真的收斂性問題和數據溢出; 4)更為嚴重的是, 上述方法的仿真收斂性高度依賴于迭代步長的設置, 不同的步長會產生不同的仿真結果, 從而影響仿真結果的一致性和準確性.因此, 采用憶阻器數值近似模型的憶阻網絡在電子設計自動化工具(electronic design automation,EDA)仿真中, 經常會遇到收斂性問題, 特別是在長時演化的場景下該問題會尤其明顯(詳見本文第4節分析), 且隨著憶阻器單元數量的增加, 收斂性對仿真的影響也愈加顯著.此外, 不同的仿真平臺和仿真設置(如步長和容忍精度)下, 采用相同的數值模型可能會造成仿真結果的差異.
鑒于上述非解析建模方法的局限性, 本文提出一種新型的、適用于大規模憶阻網絡的憶阻器單元解析建模策略, 該策略可提高憶阻網絡長時演化仿真的收斂性, 從而加速憶阻網絡在神經網絡等新興領域的應用.與傳統非解析建模方法不同, 該策略首先基于廖世俊教授[24]提出的同倫分析方法(homotopy analysis method, HAM)、以解析的同倫級數的形式從狀態方程求解憶阻器單元的核心參數—狀態變量的解析近似解, 然后利用所求解得到憶阻器單元的解析近似模型, 最后將該模型嵌入到憶阻網絡中的每個憶阻器單元.該策略具有如下特點: 1)其解具有閉合表達式特征, 即解為顯函數, 具有解析性; 2)其解的近似誤差進行了優化,從而實現了模型以及基于模型的憶阻網絡的收斂性優化.通過憶阻器單元和憶阻矩陣網絡的仿真實驗, 驗證了該策略與傳統非解析方法相比所具有的解析性和收斂性優勢, 并進一步分析了該策略在不同類型憶阻單元器件和輸入下的高擴展性.
本文結構如下: 第2節提出了一種基于HAM的解析建模策略; 第3節以經典的惠普(HP)Pt/TiO2/Pt型憶阻器單元[1]為例, 具體說明了求解HAM模型的詳細步驟; 第4節提出了一種適用于憶阻網絡中憶阻器單元模型的驗證準則, 并通過憶阻器單元和大規模憶阻矩陣網絡的長時演化實驗,驗證了該策略的適用性以及解析性和收斂性優勢(與傳統非解析方法相比); 第5節就各實驗結果以及該策略的性能優勢和擴展性進行了討論; 第6節為結論部分.
HAM提出的初衷是為了求解各種非線性問題[24?27], 其基于同倫理論, 可將一個非線性方程通過同倫轉換成數個線性方程, 同時可靈活選擇適當的同倫形式以減少計算的迭代次數, 從而使復雜問題易于解析求解.結合HAM的上述特征以及憶阻器單元原始物理模型中狀態方程的非線性特性, 所提建模策略采用HAM來解析求解狀態方程, 從而得出狀態變量的解析近似解, 以構建適用于大規模憶阻網絡的憶阻器單元解析近似模型, 最終提高憶阻網絡的收斂性.
具體而言, 從憶阻器單元的定義和物理性質可知, 憶阻器單元和憶阻網絡的動態演化主要由憶阻器單元的狀態方程決定[20], 該方程描述了狀態變量的動態演化過程, 一般為非線性微分方程(nonlinear differential equation, NDE)[28].
為解析求解狀態變量, 進而通過狀態變量得到完整的解析模型, 本文所提策略包含以下步驟.
首先, 根據需要建模的憶阻單元器件類型及相應的原始物理模型(即物理方程)、器件參數, 得到狀態方程 N [x(t)]=0.其中, N 為非線性微分算子; x (t) 為狀態變量, 是憶阻器單元各物理變量中最核心的變量, 決定了憶阻值、輸出響應等其他物理變量的動態演化.
其次, 利用HAM構造高階同倫形變方程

其中, xm(t) 項表示第 m 階同倫系數( m ∈N.如果m=1 則 xm=0 , 否 則 χm=1 ); φ (t;q) 是 用 來 逼近待解狀態變量 x (t) 的近似函數; q 表示用于構造同倫轉換的嵌入參數( q ∈[0,1] ) ;是為了優化所求解的收斂性而引入的收斂控制參數; L 為輔助線性算子.
然后, 對(1)式兩邊進行逆線性算子 L?1的運算, 即

可得

其中, Dm定義為第 m 階同倫導數算子

線性算子 L 的選擇具體取決于狀態方程的表達式特征.通常情況下, 不同類型憶阻器單元的狀態方程本質為一階NDE[28], 所以取 L 為一階微分, 相應地, L?1取一階積分.
接下來, 將(3)式中 xm(t,?) 代入到原狀態方程, 通過最小化狀態方程的離散平方剩余誤差(即殘差平方和) EN(?) 以 求解最優 ? , 即

其中,

式中, N 為近似階數, NP為均勻選擇的時間點數目, 第 j 個時間點可表示為 tj=j?t , 且

其中, Γ 為原狀態方程 N [x(t)] 在具體仿真場景下的時域范圍.
為了提高解的收斂性, 最優 ? 需滿足如下要求:1)隨著同倫近似階數 N 的增加, 所得 ? 須迅速收斂到某個固定值, 避免在該值附近振蕩; 2)在給定近似階數N 的前提下, EN(?) 須盡可能小, 以減少誤差從而進一步優化 xm(t) 的收斂性.
最后, 聯立(3)式和(4)式及所得最優 ? , 可求出收斂性優化后的狀態變量 x (t) 的 N 階解析近似解 xN(t) , 其具有由解析同倫級數構成的解表達式:

其中 x0(t) 為給定的 x (t) 初始值.
以第2節所述通用求解步驟為基礎, 本節以經典惠普HP Pt/TiO2/Pt型憶阻器單元[1]為例, 詳述所提HAM建模策略的具體建模流程, 以得到HP憶阻器單元的解析近似模型.
HP Pt/TiO2/Pt型憶阻器中, Pt是兩端的金屬鉑電極, TiO2是夾在兩個Pt電極之間且摻雜了氧空位的二氧化鈦開關層薄膜.該憶阻器單元的物理機理是在輸入導致的外部電場作用下, 氧空位在二氧化鈦開關層中作非線性漂移運動, 其原始狀態方程[1]可由NDE表示為

其中, IM(t)=Asin(ωt) 為輸入電流, D 是二氧化鈦層的厚度; x (t) 是該憶阻器單元的狀態變量, 其定義為摻雜層和 D 之比; RON表示低憶阻值( x =1 );μV表示電場作用下的平均氧空位遷移率;fw[x(t)]為窗函數用于確保求解得到的 x (t) 處于合理的物理邊界, 即x∈[0,1].
本節以Joglekar窗函數[29](簡稱J窗函數)為例, 其中 fw[x(t)]=1?(2x?1)2P, 非線性控制參數 P =2 , 此時狀態方程(9)具有強非線性[29], 無法求出完全解析解.值得注意的是, 盡管J窗函數具有一定的局限性, 如存在所謂的“邊界效應”(即當 x (t)=1 或 x (t)=0 時, x (t) 將永遠固定在這兩個邊界)[30], 但由于該函數表達式簡單從而便于列出下文所求的 x (t) 表達式(由(9)式可見, 解的復雜度與窗函數的復雜度正相關), 其仍被選為例子以闡明所提HAM建模策略的具體求解流程, 對于其他重要的窗函數[30,31]而言, 求解流程是類似的。此外, 該策略也同樣適用于其他窗函數(詳見圖1和圖3—圖5仿真所采用的Prodromakis[31]和Biolek[30]窗函數以及5.3節的討論), 具體求解流程并無區別, 這是因為窗函數雖然種類繁多, 但本質上均為非線性多項式[20], 而HAM提出的初衷正是為了處理包括非線性多項式在內的各種非線性問題[24?27].
將 P =2 的J窗函數代入(9)式中, 此時HP憶阻器單元的狀態方程為

聯立(4)式和(9)式, 可得

將(11)式代入(3)式后, 再聯立(8)式, 可得到 x (t) 的三階解析近似解(注意, 此處取三階的目的僅為舉例說明HAM策略的求解流程, 也可取其他階數.其他階的解具有相似的表達式.由于篇幅所限, 本文不再羅列):

其中, px3,1—px3,3和 β 是為了化簡(12)式而形成的中間參數, 表達式為

通過最小化 E3(t,?) 求解到最優 ? 后, 可將HP憶阻器單元[1]的器件和輸入參數代入到(12)式—(15)式, 從而得出 x3(t) 的閉合表達式(本質上是解析近似解).為方便通過實例的方式更進一步闡明該策略的建模步驟, 附錄中給出了將器件和輸入參數[1]代入(12)式后得出的狀態變量 x3(t) 、相應記憶值(t) 和輸出響應(t) (等于 IM(t)×(t) )的完整表達式.上述計算過程使用Mathematica?[32]數學軟件完成.
為驗證所提HAM建模策略在大規模憶阻網絡中的適用性, 以及與傳統非解析建模方法相比具有的解析性和收斂性優勢, 利用第3節所求的HAM模型, 分別進行憶阻器單元和憶阻矩陣網絡的長時演化實驗.具體而言, 采用配置為Intel Xeon?E7-8870 CPU (2.4 GHz), 32 GB DRAM (DDR4-2666 MHz)的工作站、以及Hspice?[33]EDA工具,分別進行了長時動態演化(圖1、圖2—圖4)和演化時間統計(圖5)的仿真, 并將實驗結果與傳統非解析數值模型進行了比較(圖1和圖5)和分析.
由于憶阻器是組成憶阻網絡的最基本單元, 其收斂性對網絡整體的收斂性具有顯著影響, 因此首先在單個憶阻器(憶阻器單元)中來驗證所提策略的適用性.根據憶阻器的輸入-輸出動力學理論[34],當輸入為周期性長時信號時, 憶阻器單元的輸出響應, 即狀態變量 x (t) 的動態演化, 也必須是周期性的, 且在不同輸入周期之間須保持一致.
然而, 從圖1(a)可看出, 長時間輸入激勵下,基于傳統非解析模型[35]的狀態變量 x (t) 呈現出非周期性的演化現象, 違反了憶阻器的輸入-輸出動力學理論.該現象是由構建非解析模型過程中采用數值積分法引入的離散、截斷和舍入誤差引起的.這些誤差在每個輸入周期內不斷累積, 導致 x (t) 偏離其物理邊界 x ∈[0,1] , 最終造成仿真結果不正確.具體而言, 當輸入極性發生變化時, 非解析模型持續對輸入進行積分, 從而需要一個等值的積分量(對電流輸入而言, 該積分量為電荷), 以抵消極性改變前所累積的積分量.只有在每次抵消完成后,x(t)才會再次隨著輸入開始演變.長時場景下, 上述抵消效應對演化造成的負面影響更為明顯.
與之相反, 由于HAM模型本身的解析特性,其解不存在上述誤差.故如圖1(b)所示, 基于解析HAM模型得到的 x (t) 呈現出理想的周期性演化特征, 表明所提HAM策略比傳統非解析方法更適合于長時演化下的憶阻器單元建模.值得注意的是: 1)不同于第3節用于舉例的J窗函數, 圖1的實驗采用了優化型Prodromais窗函數[31,34], 基于該函數的HP憶阻器解析建模流程與Joglekar窗函數下的(12)式—(15)式完全一致, 只是窗函數表達式不同而已, 從而也驗證了本策略對不同窗函數的擴展性; 2)在短時仿真且x(t)沒有到達1或者0的前提下, 部分非解析模型(如HP模型[1])的演化會呈現周期性特性.但圖1的仿真前提為更加苛刻的長時演化且狀態變量演化已到達1, 故非解析模型并不適用, 這是因為非解析模型所產生的各類數值誤差會在長時間仿真下不斷累積[28], 而憶阻網絡的典型應用場景均為長時[20]; 3)同樣是基于優化型Prodromakis窗函數的HP憶阻器模型[31], 但分別采用由本策略得出的HAM解析模型和Prodromakis數值模型得到的長時演化相比, HAM模型具有明顯的周期性,從而克服了優化型Prodromakis模型由于各種數值誤差所產生的非周期性現象, 體現了該策略的價值和意義所在.

圖1 輸入電流 A sin(ωt) 激勵下, HP Pt/TiO2/Pt憶阻器單元的傳統非解析和HAM解析模型中各自狀態變量x(t)的長時演化對比 (a)非解析(數值積分)模型[35]; (b)解析的HAM模型.兩種模型采用相同的仿真配置(Hspice仿真)和器件參數: 優化型Prodromakis窗函數[31,34] (f(x)=其中非線性控制參數P=7 ), x 0=0.5,RON=1 00 Ω, δ =10 (高低憶阻值之比), A =1 mA, 以及 ω =0.1/(2π) rad/sFig.1.Dynamic evolution comparison of the state variables based on the traditional non-analytic[35] and the analytic HAM models of the HP Pt/TiO2/Pt memristor, respectively, under an input current A sin(ωt) : (a) Non-analytic (numeric integration) mode; (b) analytic 3-order HAM model.The optimized Prodromakis window[31,34] (f(x)=where nonlinearity controlling parameter P =7 ) and the device parameters:x0=0.5,RON=1 00 Ω, δ =10 (the ratio of high memristance to low resistance), A = 1 mA, and ω=0.1/(2π)rad/s, are used in the two models (Hspice simulation).
本節采用大規模憶阻矩陣網絡作為基準電路,以驗證所提HAM建模策略在憶阻網絡中的適用性, 該網絡已被廣泛用于評估憶阻網絡的收斂性[20],實驗分為直流(DC)和交流(AC)輸出響應.此外,還對基于HAM與傳統非解析方法的憶阻網絡仿真時間進行了定量對比和分析.
如圖2所示, M ×N 憶阻矩陣網絡含有(M+1)(N+1) 個 電路節點, 2 MN+M+N 個相同的憶阻器單元, M 和N 分別表示網絡中橫向和縱向單元的數目.假設每個單元模型含有 a 個內部節點, 則整個網絡總節點數為


圖2 M ×N 憶阻矩陣網絡.電流或電壓信號可作用于網絡的“in”輸入端Fig.2.A M ×N memristive matrix network.The current or voltage is inputted into the “in” terminal of the network.
4.2.1 適用于憶阻網絡中憶阻器單元模型的驗證準則
由(17)式可知, 該憶阻網絡的復雜度可通過設置 M 和 N 的取值來動態調節, 即 M 和 N 的取值越大則網絡總節點數越多, 網絡越復雜, 從而越容易產生收斂性問題.如 M =N=100 時, 該網絡具有 2 0200 個憶阻器單元和 3 0000α+20001 個節點.面對如此復雜的非線性憶阻網絡, 不易采用傳統的直接計算法對輸出響應進行分析進而驗證憶阻器單元模型的適用性, 這是因為基于憶阻器的已有理論無法直接計算出該網絡的輸出響應[28].
為解決上述問題, 受憶阻器的閉合理論[36](即任何具有地端的憶阻系統, 可等效為一個嚴格遵循憶阻器fingerprints準則的“擴展型”憶阻器)啟發,并結合Biolek等[20]的研究, 提出了一種改進的間接法驗證準則, 以對憶阻器網絡中憶阻器單元模型的適用性進行多方位驗證, 具體如下:
1)如果憶阻網絡中所有憶阻器單元都相同(包括相同的初始憶阻值), “in”輸入端的電壓或電流、以及通過零時瞬態仿真( t =0 時)得到的整體網絡憶阻值(此時網絡已等效為一“擴展型”憶阻器)可被驗證, 具體方法為: 以傳統線性電阻代替所有憶阻器單元, 采用電路分析理論(將在下文詳述)計算初始憶阻值, 然后對比理論計算和仿真實驗的結果;
2)根據憶阻器單元的定義及其物理特性, 當輸入信號為0時, “in”輸入端的電壓或電流必須為0, 即基于閉合理論, 當整個網絡等效為一個“擴展型”憶阻器后, 該憶阻器的I-V曲線須經過零點且呈現出壓縮遲滯環特征[36];
3)仿真過程中, 憶阻器單元和等效“擴展型”憶阻器的所有狀態變量必須時刻保證在合理的物理邊界內, 如HP憶阻器單元中 x (t)∈[0,1].
如仿真結果違反上述準則, 表明網絡中使用的憶阻器單元模型存在確定的收斂性問題.值得注意的是, 該準則不僅適用于如圖2所示的憶阻矩陣網絡分析, 由于其源自憶阻器基本理論[20,36,37], 理論上可被用于其他具有不同架構的憶阻網絡.
4.2.2 實驗驗證
為計算“in”輸入端與地之間的等效初始憶阻值, 采用已被驗證、可有效分析電阻網絡的Lattice Green函數[38], 作為準則1)中的電路分析理論.由該函數計算所得的憶阻值為

其中, R0是憶阻矩陣網絡中每個憶阻器單元的初始憶阻值, γ ≈0.5772 表示歐拉-馬斯克若尼常數.
圖3是將HAM解析模型嵌入到圖2憶阻矩陣網絡中的每個憶阻器單元后, 進行瞬態仿真得到的“in”輸入端電壓Vin的長時演化過程, 其中各單元均為HP Pt/TiO2/Pt憶阻器.可以看出,t=0時得到的初始電壓與依據(18)式計算出的電壓值(即250 μ A × R(50,30) )相同, 均為7.5 V, 符合驗證準則1), 從而驗證了HAM建模策略的精度, 及其在DC直流演化場景下的憶阻矩陣和大規模憶阻網絡中的適用性.由圖3還可看出, t =5 s后Vin演化趨于穩定, 符合憶阻器閉合理論[37], 即在直流輸入下, 憶阻器的輸出最終會演化到一個固定值并持續保持穩定.

圖3 250 μ A DC電流Iin輸入下, 憶阻矩陣網絡中“in”輸入端電壓Vin的瞬態長時演化 M =50 , N =30.除采用Biolek窗函數[30]( P =2 )、初始憶阻值 R 0=16.6k? ,δ=1000之外, 其余仿真參數均與圖1相同Fig.3.Transient long-time evolution of the “in” terminal voltage (Vin) of the matrix network under a 250 μ A DC current input (Iin).M =50 , N =30 , and the simulation parameters are the same as those in Fig.1 except the Biolek window[30] ( P =2 ), initial memristance R 0=1 6.6 kΩ,and δ =1000 are used.
圖4 為“in”端輸入正弦電流信號時, 采用HP Pt/TiO2/Pt憶阻器的HAM解析模型作為基本單元, 進行瞬態仿真得到的憶阻矩陣網絡輸入端電壓Vin的演化過程(圖4(a))和相應的I-V遲滯環(圖4(b)).由圖4可見, 此時I-V遲滯環呈現壓縮形態, 且當輸入電流為0時, 對應的Vin也為0, 即I-V遲滯環過零點, 上述仿真結果符合驗證方法2).另外, 依據狀態變量對輸出I-V響應的影響特性[37],以及圖4(b)中I-V響應曲線在多個周期下所呈現出的重疊(即多周期的I-V曲線重合)遲滯特性可知, 將矩陣網絡等效為單個“擴展型”憶阻器后, 該憶阻器的狀態變量在多個周期內持續地振蕩演化,且演化均在其合理物理邊界(即[0, 1])內, 滿足驗證方法3).上述圖4的分析表明, 基于HAM模型的憶阻網絡滿足驗證方法2)和3), 從而驗證了HAM建模策略對交流演化下的憶阻矩陣網絡以及大規模憶阻電路的適用性.

圖4 輸入幅度為1 mA、頻率為1 Hz的AC(正弦)電流時所對應的憶阻矩陣瞬態長時演化 (a) 憶阻矩陣網絡中“in”輸入端電壓Vin; (b) 相應的I-V曲線(重疊壓縮遲滯環).各仿真參數均與圖3相同Fig.4.Transient long-time evolutions of memristive matrix network under a 1 mA, 1 Hz AC (Sinusoidal) current input:(a) “in” voltage Vin; (b) corresponding I-V curves (compressed hysteresis loops).The simulation parameters are the same as those in Fig.3.
除適用性外, 通過對HAM和傳統非解析(數值積分)模型Hspice仿真時間的比較, 還可驗證HAM建模策略相對于傳統非解析方法的解析性和收斂性優勢.由圖5可見, 在輸入分別為DC直流(圖5(a))和AC正弦(圖5(b))輸入且同等憶阻器單元數目的情況下, HAM比傳統非解析模型所用的仿真時間更少, 因而其計算效率更高.此外, 如圖5中黑色虛線所示, 當直流和交流輸入下憶阻器單元數分別增加到5.01 × 105和5.10 × 103時, 基于非解析模型的仿真會出現崩潰和數據溢出的現象, 表明仿真失敗; 相比而言, 由于HAM策略的解析性和收斂性優化, 基于HAM模型的仿真均未出現上述現象, 表明仿真成功.需要注意的是, 除了模型的解析性和收斂性之外, 仿真硬件平臺的性能也會影響上述運行時間(例如, CPU的運行速度會影響仿真時間, DRAM的容量會影響數據溢出發生時憶阻器單元的數量), 盡管如此, 不同平臺仿真得到的演化趨勢和分析結果仍與圖5相同.
本節對圖3—圖5的仿真結果進行深入討論,進一步分析所提HAM解析建模策略相比傳統非解析建模方法所具有的收斂性優勢, 并驗證該策略對其他憶阻網絡(基于不同類型的憶阻器單元)的擴展性.

圖5 分別采用HAM和傳統非解析[37](數值積分)模型進行憶阻矩陣仿真, 隨著憶阻器單元數目的增加, 仿真時間的對比 (a)和(b)分別表示DC和AC (正弦)電流輸入下的仿真時間.此對比分析基于圖3(DC)和圖4(AC)所示Vin的動態演化Fig.5.Comparisons of running time between simulations using the HAM and the traditional non-analytic (numeric integration)[35] and models with increasing memristor cells:(a) and (b) show the time under the DC and AC (Sinusoidal) current inputs, respectively.The comparisons are adopted to analyze the dynamic evolutions Vin as shown in Figs.3(DC) and 4 (AC).
產生圖1實驗結果的內在原因為, 傳統非解析建模方法受硬件仿真平臺的機器精度限制, 存在最大有效仿真時間的問題[39], 如超過該時間, 仿真會有較大的累積誤差.這是因為, 數值積分等非解析方法的有效性建立在仿真時間步長必須有下限的基礎上, 但由于機器精度的限制, 為了將累積誤差對狀態變量 x (t) 的影響約束在合理物理邊界內, 不得不在給定的仿真時間內限制步長, 因此步長下限無法得到保證, 特別是對于長時演化仿真該問題會愈加嚴重, 故由該類方法導出的非解析模型不適用于憶阻器單元的長時演化仿真.而HAM建模策略由于求解的 x (t) 是以解析的同倫級數來表示, 其嚴格解析從而沒有步長下限的約束, 因此對 x (t) 的仿真時間不敏感, 可適用于長時演化的分析.
由4.2節分析可知, 當 M =N=100 時, 圖2所示的憶阻矩陣網絡具有 2 0200 個憶阻器單元和30000α+20001個節點.傳統的非解析憶阻器單元模型一般至少包含一個內部節點(即 α ≥1 )以及一個內部積分器(一個節點對應一個積分器)來模擬 x (t) 的動態演化[28], 故基于該模型的憶阻網絡含有大量(由(17)式可知, 具體為(1+2α)MN+(1+α)MN+1個)積分方程待求解, 而通過HAM建模策略求得的 x (t) 如(12)式所示為一閉合表達式,因此基于HAM的 x (t) 無需解積分方程且其內部節點只有一個, 即 α =1.由于 α 的不同, 導致以Pt/TiO2/Pt憶阻器為基本單元[1]的憶阻網絡, 其基于HAM和傳統非解析模型的電路節點總數分別為 3 0401 個( α =1 ) 和 5 35401 個( α =26[35]).考慮到上述差異的顯著性, 可見基于HAM模型的憶阻網絡需要求解的節點電壓方程和支路電流方程也大幅減少, 從而可提高仿真收斂性和節省仿真時間(如圖5所示為例).
圖1和圖3驗證了所提HAM建模策略對不同窗函數的擴展性, 圖3和圖4驗證了該策略對不同輸入(即直流和交流(正弦))的擴展性, 進一步地, 為驗證該策略對器件類型的擴展性以及相對于非解析方法的擴展性優勢, 將其應用于如表1所列的基于不同類型憶阻器單元所構建的矩陣網絡仿真, 這些憶阻器單元具有不同的器件結構、電極-開關層材料和物理機理.
由表1所列的仿真時長對比可看出:
1)相比傳統非解析模型, 采用HAM模型的各憶阻網絡仿真時間更少;
2)隨著原始物理模型的復雜度(即憶阻器單元的物理機理和相應物理方程)從Pt/Ta2O5/TaOx/Pt增加到Ag/TiO2/ITO, 憶阻網絡的仿真時間也相應增加, 且使用傳統非解析模型的網絡在某些情況下會遇到仿真收斂問題(即表1中所示的“不收斂”), 而采用HAM模型的所有網絡仿真均收斂;

表1 直流演化場景下, 基于不同類型憶阻器單元的憶阻網絡仿真時長比較Table 1.Running time comparisons of the DC evolution for memristive networks with different types of memristor cells.
3)由1)和2)可知, 所提HAM建模策略可適用于多種類型的憶阻器單元, 以及由這些憶阻器單元構建的大規模憶阻網絡.
值得注意的是, 憶阻器單元中狀態變量的定義取決于器件具體的器件結構、電極-開關層材料和物理機理等.除本文用于舉例的HP Pt/TiO2/Pt憶阻器單元之外, 其他憶阻器單元的狀態變量有導電絲的長度、隧穿的勢壘寬度等.雖然不同類型的憶阻器單元, 其定義具有一定差異性, 但所提HAM策略的建模步驟是類似的, 這是因為憶阻器單元的狀態方程通常是NDE, 而該策略中使用的HAM對于求解各類NDE具有極高的通用性[24?27].
通過上文實驗分析和討論可知, 非解析型模型在憶阻器單元和憶阻矩陣網絡長時演化的仿真過程中, 以及不同類型憶阻網絡的仿真時間比較中,皆存在收斂性問題, 而解析HAM模型均未出現該問題.因此, 可得出合理的推論: 憶阻器單元模型的非解析性是導致大規模憶阻網絡仿真不收斂的潛在原因之一.
綜上所述, 本文提出了一種新型的、適用于大規模憶阻網絡仿真的憶阻器單元HAM建模策略,基于其解析性和收斂性優化, 該策略具有以下優點:
1)保證了由原始物理模型導出的解析近似模型的連續收斂性;
2)能夠模擬憶阻器單元和憶阻網絡的長時動態演化;
3)對具有不同類型窗函數、輸入和憶阻器單元的大規模憶阻網絡均具有靈活的擴展性.
基于此策略和所提憶阻器單元模型的驗證準則, 通過各對比實驗, 揭示了大規模憶阻網絡仿真中不收斂現象的潛在原因之一是所用憶阻器單元模型的非解析性.該策略有望進一步用于憶阻器單元和憶阻網絡的長時性能分析和優化, 促進憶阻器單元和憶阻網絡在神經網絡等新興技術中的應用.
附 錄
利用所提HAM建模策略, 可求出Joglekar窗函數下Pt/TiO2/Pt憶阻器單元[1]的三階解析近似解為(將器件和輸入參數[1]代入(12)式后得出)

隨后, 基于(A1)式以及原始物理模型中的端口方程[1]

其中 δ 是高( x (t) =0)低( x (t) =1)憶阻值之比, 可得相應憶阻值(t) 的三階解析近似解為
