一种适用于大规模忆阻网络的忆阻器单元解析建模策略*

胡炜 廖建彬 杜永乾

1) (福州大学物理与信息工程学院, 福州 350116)

2) (集美大学轮机工程学院, 福建省船舶与海洋工程重点实验室, 厦门 361021)

3) (西北工业大学深圳研究院, 深圳 518057)

4) (西北工业大学电子信息学院, 西安 710072)

忆阻网络是一种基于忆阻器单元的大规模非线性电路, 在下一代人工智能、生物电子、高性能存储器等新兴研究领域发挥着重要作用.描述忆阻器单元物理和电学特性的模型对忆阻网络的性能仿真具有显著影响.然而, 现有模型主要为非解析模型, 应用于忆阻网络分析时可能存在收敛性问题.因此, 提出了一种基于同伦分析法(homotopy analysis method, HAM)的忆阻器单元解析建模策略, 该策略具有解析性和收敛性优化的特点, 可提高忆阻器单元和相应忆阻网络的收敛性.此外, 还提出了一种面向忆阻器单元模型的验证准则, 以验证模型在大规模忆阻网络中的适用性.通过忆阻器单元和忆阻矩阵网络的长时演化实验以及与传统非解析(数值)方法的比较, 验证了所提策略的解析性和收敛性优势; 利用不同类型忆阻器单元和输入的实验, 验证了该策略的扩展性.进一步地, 基于上述实验, 揭示了忆阻网络仿真出现收敛性问题的潜在原因.该策略可应用于基于忆阻网络的新兴研究.

1 引 言

忆阻器是一种由多层纳米薄膜组成的新型微纳器件, 具有高非线性和非易失性开关的特点, 能够以非易失方式“记忆”流经的电荷和磁通量, 具体体现为内部状态变量和忆阻值的动态演化[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节为结论部分.

2 基于HAM的解析建模策略

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) 初始值.

3 求解HAM解析模型

以第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]数学软件完成.

4 验 证

为验证所提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)和分析.

4.1 忆阻器单元中状态变量的长时演化实验

由于忆阻器是组成忆阻网络的最基本单元, 其收敛性对网络整体的收敛性具有显著影响, 因此首先在单个忆阻器(忆阻器单元)中来验证所提策略的适用性.根据忆阻器的输入-输出动力学理论[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).

4.2 忆阻矩阵网络输出响应的长时演化实验

本节采用大规模忆阻矩阵网络作为基准电路,以验证所提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相同.

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).

5.1 忆阻器单元的长时演化

产生图1实验结果的内在原因为, 传统非解析建模方法受硬件仿真平台的机器精度限制, 存在最大有效仿真时间的问题[39], 如超过该时间, 仿真会有较大的累积误差.这是因为, 数值积分等非解析方法的有效性建立在仿真时间步长必须有下限的基础上, 但由于机器精度的限制, 为了将累积误差对状态变量 x (t) 的影响约束在合理物理边界内, 不得不在给定的仿真时间内限制步长, 因此步长下限无法得到保证, 特别是对于长时演化仿真该问题会愈加严重, 故由该类方法导出的非解析模型不适用于忆阻器单元的长时演化仿真.而HAM建模策略由于求解的 x (t) 是以解析的同伦级数来表示, 其严格解析从而没有步长下限的约束, 因此对 x (t) 的仿真时间不敏感, 可适用于长时演化的分析.

5.2 忆阻矩阵网络的长时演化

由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所示为例).

5.3 扩展性

图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].

5.4 导致忆阻网络仿真不收敛的潜在原因

通过上文实验分析和讨论可知, 非解析型模型在忆阻器单元和忆阻矩阵网络长时演化的仿真过程中, 以及不同类型忆阻网络的仿真时间比较中,皆存在收敛性问题, 而解析HAM模型均未出现该问题.因此, 可得出合理的推论: 忆阻器单元模型的非解析性是导致大规模忆阻网络仿真不收敛的潜在原因之一.

6 结 论

综上所述, 本文提出了一种新型的、适用于大规模忆阻网络仿真的忆阻器单元HAM建模策略,基于其解析性和收敛性优化, 该策略具有以下优点:

1)保证了由原始物理模型导出的解析近似模型的连续收敛性;

2)能够模拟忆阻器单元和忆阻网络的长时动态演化;

3)对具有不同类型窗函数、输入和忆阻器单元的大规模忆阻网络均具有灵活的扩展性.

基于此策略和所提忆阻器单元模型的验证准则, 通过各对比实验, 揭示了大规模忆阻网络仿真中不收敛现象的潜在原因之一是所用忆阻器单元模型的非解析性.该策略有望进一步用于忆阻器单元和忆阻网络的长时性能分析和优化, 促进忆阻器单元和忆阻网络在神经网络等新兴技术中的应用.

附 录

利用所提HAM建模策略, 可求出Joglekar窗函数下Pt/TiO2/Pt忆阻器单元[1]的三阶解析近似解为(将器件和输入参数[1]代入(12)式后得出)

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

其中 δ 是高( x (t) =0)低( x (t) =1)忆阻值之比, 可得相应忆阻值(t) 的三阶解析近似解为