A Study of the Applicability of Existing Compact Models to the Simulation of Memristive Structures Characteristics on Low-Dimensional Materials

The use of low-dimensional materials is a promising approach to improve the key characteristics of memristors. The development process includes modeling, but the question of the most common compact model applicability to the modeling of device characteristics with the inclusion of low-dimensional materials remains open. In this paper, a comparative analysis of linear and nonlinear drift as well as threshold models was conducted. For this purpose, the assumption of the relationship between the results of the optimization of the volt–ampere characteristic loop and the descriptive ability of the model was used. A global random search algorithm was used to solve the optimization problem, and an error function with the inclusion of a regularizer was developed to estimate the loop features. Based on the characteristic features derived through meta-analysis, synthetic volt–ampere characteristic contours were built and the results of their approximation by different models were compared. For every model, the quality of the threshold voltage estimation was evaluated, the forms of the memristor potential functions and dynamic attractors associated with experimental contours on graphene oxide were calculated.


Introduction
One of the most promising solutions in the field of building neuromorphic systems is associated with a memristor [1], a nonlinear element of electric circuits, which resistance can reversibly change depending on the electrical signal entering its input. The resistive switching effect was first implemented in 2008 [2] through the atomic structure rearrangement in local regions of nanometer size in a metal-oxide-metal structure, which demonstrated the potential complementary metal-oxide-semiconductor (CMOS) compatibility, unique scalability and high switching speed of memristive devices [3]. To date, memristor research has led to the construction of complex high-performance devices capable of simulating the behavior of complex neural networks [4], building multi-level logic elements [5], and other devices for high-speed real-time operations [6].
Depending on the material of the functional layer and the design features of the memristive element, their main drawbacks vary from structure to structure, but most of them are determined by the random nature of the processes occurring inside the device [7][8][9][10][11]. One important direction in the development of memristors is the use of promising lowdimensional materials to increase the stability of device characteristics and form its new properties [12][13][14].
However, to the best of our knowledge, the current research in the field of memristive elements based on low-dimension materials does not include the study of existing memristor models' applicability to them. In particular, the study of compact models is important because such models are used as a basis for circuit models in the development of microelectronic devices, and related to investigations of the properties of new experimental structures, for example, for investigations of dynamic attractors.
Due to the presence of a number of experimentally determined peculiarities of compact models, it is necessary to perform extraction of model parameters for each structure in order to use them. The main method used for parameter extraction is the solution of the inverse problem based on experimental measurements of the device current-voltage (I-V) curves. In general, due to the type of evolution equation, an analytical solution of the inverse problem cannot be obtained, and parameter extraction is an optimization problem.
Thus, the quality of the solution of the optimization problem in the study of memristors on low-dimensional materials determines the compact model applicability to a particular device. Therefore, introducing a function corresponding to the optimization error, it is possible to perform a comparison of the models with each other. A synthetic sample containing a set of key features of different structures, described below, was used as a simulation set to compensate for the contribution of random effects of each particular structure.
The task of evaluation of the approximation quality on a set of experimental results is important because, according to the conclusions of the "No Free Lunch" theorem [15,16], the effectiveness of models and optimization algorithms developed for them is inextricably linked to the task for which they are designed, and using them on other tasks may yield reduced effectiveness.
Examples of benchmark studies conducted to solve the optimization problem in different subject areas are presented in [17][18][19].
In this paper we performed a feature meta-analysis of memristor based on lowdimensional materials. Based on the results obtained, we conducted a study comparing linear drift models [2], nonlinear drift models [20][21][22] and threshold adaptive models [10,23,24]. In addition to the synthetic sample of the volt-ampere characteristics built on the basis of the features of the devices highlighted in various studies, the results of experiments in the study of the graphene oxide memristor characteristics, first presented in [25], are also used to investigate the descriptive power of the models. Finally, an error function is chosen and a general optimization method is proposed for the numerical experiment. Obtained results include the quality of the approximation of the volt-ampere characteristic loop, threshold voltages and the suitability of the model for the study of dynamic attractors.
The paper is organized as follows. Section 2 contains review of memristors based on low-dimensional materials and main criteria for model comparison, description of analyzed models, optimization problem statement for volt-ampere characteristic approximation, sampling of synthetic curves for model analysis, difference appearance of threshold voltages in various models, method of dynamic attractor calculation, and proposed optimization methods. Section 3 is devoted to the results of model comparison related to the approximation of synthetic and experimental volt-ampere characteristic (VAC) approximations, assessment of threshold voltages, and dynamic attractors. In Section 4 we make some concluding remarks.
At the end of the introduction, the authors would like to emphasize the novelty of the paper. To the best of our knowledge, we are the first to compare the applicability of popular compact models to the low-dimensional memristor characteristics. The article includes the requirements for a comparison based on set of test volt-ampere characteristics, proposed a general method for solving the optimization problem, and modification of the error functional to solve the optimization problem of a particular lobe of the volt-ampere characteristic. In this work we obtained the solution of the optimization problem for ten lobes of the volt-ampere characteristic of a unique shape, for the first time the study of the optimization problem for a set of close loops actually corresponding to the change of device parameters from switching to switching [11] was carried out, the threshold voltages were estimated, and for each model, for the first time for the new experimental result, an estimate of potential functions was obtained.

Review of Low-Dimensional Structures
The works on the study of vertical and lateral memristive structures [12][13][14][26][27][28][29] based on combinations of molybdenum disulphide (MoS 2 ) with different materials demonstrate a significant variation of key characteristics. For example, the R off to R on ratio varies from 10 2 to 10 8 depending on the composition and thickness of the active layer, threshold voltages range from 0.1 V to 5 V, endurance varies from 200 cycles to more than 10 4 , retention up to 10 4 s [13], and the presence of dynamic attractors [14].
Based on the results obtained from the literature analysis [12][13][14][27][28][29] and the conducted experiments, the main criteria by which it is possible to compare models are: 1.
Ability to simulate a wide range of R off R on ∼ 10 2 − 10 8 ; 2.
Ability to accurately simulate threshold voltages both for the transition fromR on to R off and vice versa; 3.
Ability to simulate a dynamic change of volt-ampere characteristic curvature near threshold voltages; 4.
Ability to simulate the dynamic attractors; 5.
Stability of key values of memristor simulation result on 2D structure relatively insignificant changes of the contour. In particular, maintaining a stable R off R on 6.
The result of the simulation should make sense in the original treatment of the parameters by the authors of the model, including the associated physically justified constraints.

Models Analysed
To simulate a memristor using compact models, it is necessary to describe the relationship between the parameters of the equations and the real parameters of the device. The focus of compact modeling is the simulation of the I-V curves of the system, the general approach involves solving the following system of equations: where X is the control signal, Y is the response, F is the function linking the input signal and the response (resistance or conductivity), g is the evolution function, f is the window function, x is the state variable of the memristor. In many cases [30] these equations are simplified to: The first and most popular memristor model was proposed in [2]. The linear drift model is still relevant today and is used as the basis for other models. It is described by a system of equations: The above equations are simple, and their influence on the appearance of the VAC can be deduced trivially from general considerations.
However, in addition to the advantages of good interpretability, computational simplicity, and stability, the model has a number of key drawbacks that make it difficult to use.
The simplicity of the model is due to strong simplifying assumptions about the processes in the active layer of the memristor, which affect the accuracy of the model. Also, the model explicitly includes assumptions of switching physics that severely limit the range of approximated VACs.
The emergence of nonlinear drift models was designed to compensate for these drawbacks. An essential part of nonlinear drift models is based on [2]. Some of the first models were proposed by Joglekar and Wolf, Prodromakis, Biolek, Zha and their colleagues [31][32][33][34]. Each model was designed to correct for state evolution in the linear drift model using a window function with internal parameters. At present, there is a tendency to extend the original window functions and build generalized window functions on their basis [20][21][22], and it is the generalized nonlinear drift models that we will analyze.
Shi's [21] nonlinear drift model is a fractional-order memristor model with the introduction of a new window function and a fractional time derivative of the state variable into the evolution equation: where the parameter a lies in (0, 1), p is a positive real number, C 0 α t is the partial derivative of Caputo, defined by the formula Here n is an integer, n − 1 < α ≤ n, and Γ(z) is the Euler gamma function. In the analysis of this model, α is assumed to be equal to one for simplicity.
According to the authors, models with fractional derivatives better describe physical phenomena related to nonlocality, time memory, power law and weak singularity. The proposed window function avoids problems related to stick effect, and by varying the parameter a one can achieve a change in the maximum value of the window function.
Zha's generalized window function [20] is given by the formula: where the parameter b is a positive integer, p is a positive real number, and the parameter a belongs to the interval (0, 1). This window function solves the problem of the stick effect and scaling. The authors of the original paper propose a fine-tuning of this window function to obtain adjustable resolution in memristors using a sequence of input pulses. A nonlinear drift model with the Li's window function [22] was developed to overcome the stick effect, the inflexible parameter and the distorted pinched hysteresis loop problems. The window function has the following form: As the authors note, the general form of the window function allows one to adapt it to spike neural network simulations.
Another group of nonlinear drift models includes the Yang's [35] and Simmons' barrier models [36]. The latter model is considered to be the most accurate memristor model from Hewlett Packard (HP) Labs at the moment. Models from this group were omitted in the analysis, because in the Yang's model the nonlinear evolution and boundary effects are introduced with window functions, and the Simmons' barrier model is inefficient in terms of Simulation Program with Integrated Circuit Emphasis (SPICE) simulations [37].
There are also adaptive models with varying degrees of adaptivity. Such models include Yakopcic's model [24], defined by equations: Here x is the state variable, a 1,2 and b are constants, A p and A n are constants that determine the change rate of the state variable after overcoming threshold voltages, V p and V n are the absolute values of the upper and lower threshold voltages, respectively. Parameters x p and x n are restricted only to the range (0, 1).
The motivation for creating this model was the experimental confirmation of threshold voltages as well as the inability of the models [35,[38][39][40] to simulate a decrease in the lobe width of the current-voltage characteristic as the average conductance per cycle increases.
The VTEAM model [23] was created by the authors as an adaptive threshold model of a voltage controlled memristor. The current-voltage relationship in this model is defined differently for each case. We use the following combination of equations for this model: where x on < x < x off . The VTEAM threshold model is a modification of the earlier TEAM [41] model used to simulate current-controlled memristors. The reason for creating the VTEAM model was to avoid performance and reliability problems in crossbars of threshold voltage memristors, since the voltage across all memristors is the same, and the variable resistance does not affect the procedure of switching to the high impedance state.
We previously developed a mobility modification model (MMM). This model is based on the Yakopcic's model, which demonstrates competitive simulation accuracy and speed of the experimental VAC approximation. The model was obtained by multiplying the evolution equation by functions to account for inhomogeneities that modify the mobility of charge carriers: Here i Yakopcic (v) is the current-voltage relation in the Equation (15), U i (x) is the accounting function of the i-th inhomogeneity, is the number of inhomogeneities, is the effective position of the i-th inhomogeneity in the state space of a memristor, is the effective width of a given inhomogeneity. The motivation for creating this model was the presence of inhomogeneities in the active layer of a memristor that change the mobility of charge carriers [42].
The above equations were solved using the 4th order predictor-corrector method. The idea behind the predictor-corrector method is to use a suitable combination of an explicit and an implicit technique to obtain a method with better convergence characteristics. The explicit and implicit methods are third order Adams Methods [43], called Adams-Bashforth and Adams-Moulton methods respectively.

Optimization Problem
To compare the models with each other, it is necessary to solve the inverse problem for each model using the method of optimization of some functional. In general terms, with relation to an arbitrary criterion the functional can be written as: where P is the vector of model parameters, Aim(P) is the error for the criterion, λ is regularization parameter, S 0 is area of fitted i-v curve, S(P) is area of fitting i-v curve, G is the regularizer represented in Figure 1 as a function of areas.
We propose to use a regularizer related to the physical property of the curve to obtain a single solution in the parameter space of the model. This is due to the fact that the simulation result of some key features can be directly determined only by a part of the model parameters. Due to uniqueness of characteristics, regularizer based on the value of approximated contour area allows simultaneously to take into account uncharacteristic for usual contours features and to optimize values of other parameters, implementing physically adequate constraint independent from model.

Synthetic Sampling
Analyzing the results of both the experiment performed and those obtained in the works of other authors, we can conclude that the lowest accuracy of the obtained results refers to the region of switching between states and some of its vicinity. Moreover, the different character of switching can be observed not only between devices on different materials, but also within one series of switching. This phenomenon is connected with a set of effects: with peculiarities of experiment, with the equipment signal of registration, different nature of the functional domain effects, and so on. Thus, a significant part of the features that are possible to distinguish in the switching domain can refer not just to a structure of a certain type, but even to a specific device or even experiment. Consequently, due to the ambiguous nature of the switching region features (e.g., set voltage [26]), when comparing models, priority should be given to the pre-switching regions.
Thus, parts one and four of Figure 2 are considered indirectly, through the value of the regularizer, as a function of the area. The choice of a particular area is conditioned, among other things, by the estimation of the characteristic values of the readout voltage.
Therefore, the parts two and three of the I-V contour are used to construct the basic error function. The points over which the approximation functional was constructed were selected proportionally to the maximum voltage amplitude in the structure during the switching process.
For 100 voltage points proportional to the maximum voltage, the square of the current value difference (MSE) was calculated and averaged for the approximation result and the testing sample using the formula: where I k is the model current at point k, i k is the loop current at point k, n is the total number of points for which the error is estimated. Based on experimental contours, including [12,14,[44][45][46], we designed asymmetric elements of the VAC contours (Figure 3), significantly different for right and left parts.

Threshold Voltages
Threshold voltages in memristor models correspond to the values of the electric field strengths, above which the memristor resistance dynamics change significantly. This effect physically corresponds to changes within the conducting regions.
Depending on the type of model, thresholds can be included explicitly or implicitly. In the VTEAM and Yakopcic's threshold models, threshold voltages appear in the evolution equations explicitly: In the exponential drift model [47], instead of threshold voltage, a characteristic value of the electric field is accounted in the vacancy drift rate as follows: The threshold voltage is not specified in the linear drift model. In the nonlinear drift models based on it, the threshold voltage is not explicitly specified either. However, depending on the values of the window function parameters, one can observe a nonlinear evolution simulating the presence of threshold characteristics. In this case the evolution of the state variable occurs along the entire switching cycle.
The presence of the threshold voltages is confirmed experimentally, so it is necessary to account them in modeling. Taking into account the threshold values generate additional computational complexity, because due to the asymmetry of the I-V curve, we have to additionally determine a group of related parameters separately for the left and right lobe of the I-V curve.

Dynamic Attractors
Dynamic attractors are defined as dynamic equilibrium points in the state space of the memristor at fixed model parameters [48].
To find the attractors, minima of the so-called potential function of the memristor should be calculated as minima of the following function: This formula implies that the memristor is voltage controlled and f is the right part of the evolution Equation (2).
In particular, we have previously explicitly demonstrated the possibility to obtain attractors for the Yakopcic's and MMM models [10]. Also, the original paper by Slipko and Pershin [49] considers dynamic attractors for some nonlinear drift models with different window functions.

Optimization Method
As the most popular method for solving the optimization problem is a combination of multistart method for some parameters and gradient descent for the rest. Application of this approach is justified for a small number of considered switching cycles for each particular model.
The search method for an optimal solution for a set of models is determined by the complexity of each individual model, which is associated with unique parameter sets and equations complexity.
Thus, it is important to use a general global optimization algorithm. One of the well-studied global optimization methods is random search [50,51].
By random search we will imply an iterative algorithm for finding the minimum of the target function f (x), x ∈ X ⊂ R n , with its' step being defined as follows: (31) where {ζ k } +∞ k=1 is a sequence of random vectors.
A common choice of ζ k is: where {ξ k } +∞ k=1 are independent equally distributed random vectors, A k and b k are parameters (matrix and vector respectively) depending only on the results of previous iterations, e.g., in the most general case: However, a simpler dependence is usually chosen. In the following, we will assume that: Usually, ξ k is selected as normally distributed or uniformly distributed on a unit sphere.
There are theoretical results guaranteeing the convergence of this algorithm to the probability optimum under rather general assumptions, as well as estimates of the expectation of the number of steps to convergence [50,51]. For example, if we assume that functions f , A and b are continuous, and ξ k have a distribution with nonzero density on all space R n , then there is a probability convergence to a global minimum: (37) Earlier studies have been conducted [52] on estimation of expectation of number of steps τ of algorithm before achievement of accuracy > 0, i.e., f ( be the density of distribution of random variables ξ k , where µ is Lebesgue measure in R n .
We used the adaptive random search algorithm implemented in the Python programming language. The algorithm step depends on the current argument value. The step has normal distribution with zero expected value and standard deviation proportional to the absolute value of the argument.
As an experimental curve for the calculation, we used the results of experiments on the study of the characteristics of the graphene oxide memristor, first presented in [25]. The I-V contour of the memristor is presented in Figure 4.

Approximation
The feature approximation was conducted for the synthetic sample contours characteristics ( Figure 3). To separate out the features of the right and left branches of the characteristics, which approximation quality we want to obtain, the calculation of the MSE on the branches R on and R off of the corresponding lobe was used as the error function on the criterion from Equation (26). From the results of the approximation (Figures 5-8, Table 1) of the right and left parts of the VAC, it can be seen that the adaptive models show significantly better results to the nonlinear models. It is important to note that, despite this, for some features the nonlinear drift models show superior results, in particular Shi's [21] model, which justifies their application to some modeling problems.     However, in order to evaluate the optimization quality, it is necessary to make an assessment of the model stability to small changes of the test contours. To assess how sensitive the solution of the optimization problem is under such conditions, the curvature of the synthetic contours was artificially increased to 110% with a step of 2% and an insignificant change in amplitude (Figure 9).
Adaptive models remain their approximation accuracy given in the Table 1 with the increase of synthetic curves curvature. In turn, both the linear model and the nonlinear models built on its basis showed (Table 2) either no response to the change of contour curvature (which corresponds to the smooth growth of the error with curvature increase), or revealed other, significantly different from the previous state, shapes of the resulting contours (the change of metric occurred by leaps) with corresponding sharp change of model parameters.

Threshold Voltages
The threshold voltage was estimated within the framework of approximation by models with an explicit threshold voltage parameter. The results of the approximation are presented in the Table 3. The best result was obtained by the mobility modification model for the positive threshold voltage on the second loop: 0.87 V on the model loop and ∼1 V on the experiment. In other cases, the estimates of the threshold voltages on the loops are given, as a rule, inaccurately and differ from the experimental ones by an order of magnitude: ∼0.1 V instead of 1 V.

Attractors
To construct the potential function of the memristor for the experimentally obtained curve, it was approximated by the basic models. The optimal values by which the potential function is constructed are given in the corresponding tables (Tables 4-10). x start 0.20  x start 0.14 For the most accurate approximating models a similar result was obtained: the attractor is located at the right end of the region of acceptable values of the state variable ( Figure 10). Less accurate models show divergent results ( Figure 11): both different positions of attractors and their complete absence. The Yakopcic's model is omitted as its potential function has the same form as that of MMM.

Discussion
In this paper we performed a feature meta-analysis of memristor based on lowdimensional materials. Based on the results obtained, we conducted a study comparing linear drift models [2], nonlinear drift models [20][21][22] and threshold adaptive models [10,23,24]. A global random search algorithm was used to solve the optimization problem, and an error function with the inclusion of a regularizer was developed to estimate the loop features. Based on the characteristic features derived through meta-analysis, synthetic volt-ampere characteristic contours were built and the results of their approximation by different models were compared. For every model, the quality of the threshold voltage estimation was evaluated, the forms of the memristor potential functions and dynamic attractors associated with experimental contours on graphene oxide were calculated.
As a result of the study, we identified several main requirements to the model. The basic ones are the possibility of simulation of a wide range of R on /R off , simulation of curve conductivity dynamics, stability of simulation result with respect to small perturbations of a curve. To evaluate these parameters, we performed optimization by the mean-square error of the simulation result from the approximated contour with a regularizer. To evaluate the approximation quality of different features, the optimization was performed separately for the left and right parts of the volt-ampere characteristic. The use of the regularizer was required to bound the parameters weakly related to the features of a particular area of the VAC, as well as to keep the shape of the VAC contour physically reasonable. As a result, it was demonstrated that modeling of memristor features should be based on the application of adaptive models that can, in addition to high approximation accuracy, take into account the small change in device characteristics from switching to switching.
The experimental contour of the volt-ampere characteristic was used to evaluate the accuracy of threshold voltage calculations as a result of solving the approximation problem for threshold models and to evaluate the existence of dynamic attractors as a result of the approximation by all models considered. For considered adaptive models a similar result was obtained: the attractor is located at the right end of the region of admissible values of the state variable. On the contrary, linear and nonlinear models show divergent results: both different positions of attractors and their complete absence. It is shown that for the tests given in the article the best results with respect to the chosen comparison parameters were shown by the models presented in the articles [10,23].