Modeling of Memristors under Periodic Signals of Different Parameters

Abstract

In this paper, the problem of modeling memristors is studied. Two types of memristors with carbon and tungsten doping fabricated by the Knowm Inc. are tested. The memristors have been examined with either sinusoidal or triangle voltage wave periodic excitation. Some different frequencies, amplitudes and signal shapes have been applied. The collected data have been averaged and subjected to high frequency filtering. The quality of measurement data has also been discussed. The averaged measurement has been modeled using three popular memristor models: Strukov, Biolek and VTEAM. Some additional feathers to the considered models have been proposed and tested. Memristor is usually modeled by a set of algebraic-differential equations which link both electrical values (i.e., voltage and current) and the internal variable(s) responsible for the element dynamics. The interior-point with box constrains optimization method has been used to obtain the optimal parameters of the memristor model that fit best to the collected data. The results of the optimization process have been discussed and compared. The sensitivity to the different frequency range has been also examined and reviewed. Some conclusions and future work ideas have been postulated.

1. Introduction

Memristor devices are an emerging topic that have attracted many researchers to study the possibility of using them in applications including in-memory computations, neuromorphic computing and also as non-volatile memory elements. This attraction has boosted up right after the discovery of the real memristors in May 2008 by a team of scientists from HP Labs, under the leadership of R. Stanley Williams [1]. It is claimed that memristors can be used for building efficient memory structures with larger capacity and improved performance [2]. It is shown that memristor based memories may have densities two orders of magnitude larger than RAM memories. The in-memory computation concept is very interesting considering the problem well known as the von Neumann bottleneck, i.e., the restrictive limit of data bandwidth between the CPU and RAM [3]. Memristors are also considered as a good candidate for analog filters tuned by changing the resistance of the memristors [4], switches, registers, and many more. Furthermore, networks of memristor based oscillators may achieve the processing power of biological systems due to their high connectivity and capabilities of carrying out parallel computations [5]. The bio-inspired or even brain-inspired computation concept seems to be very attractive for the solution of many complex tasks with higher efficiency then contemporary computer systems [6,7].

Memristor is a two-terminal passive element that has the ability to retain its state. In general, the state is expressed by its resistance, which depends on the history of the magnetic flux or the electric charge of the element. This concept has been introduced by L. O. Chua in their famous paper published in the year 1971 [8]. In other words, when the power supply is switched off, the element retains the state (no charge or no flux in the element). This element is the complement to four known elements R, L and C. Chua has defined a number of properties expressing the element to be considered as a memristor, known as the memristor fingerprints [9]. One of the properties is that under any bipolar periodic stimulation one gets a particular hysteresis loop that passes through the origin on the $i-v$ plane. This property expresses the element to be passive (dissipative) and explains its resistance variation. The frequency of the bipolar stimulus changes the loop’s shape. Its area decreases monotonically when f grows [10]. In theory when the frequency of the excitation signal goes to infinity the loop becomes a curve passing through zero. The real memristors exhibit this phenomena already for very large frequencies.

The HP lab proposed the sandwiched structure of metal-insulator-metal (MIM). The basic idea is that when a strong electric field is applied to the element, a conductive filament is formed in a dielectric layer, which under small electric fields works as an insulator. The oxide film is less than 5 nm thick and contains two layers. First insulating ${\mathrm{TiO}}_2$ and second with the high conductivity property oxygen-poor ${\mathrm{TiO}}_{2-\mathrm{x}}$. Then under the voltage that in such small dimensions creates the high electric field in the active layer. The excited electric field is sufficient that oxygen vacancies start to drift, shifting the dividing line between ${\mathrm{TiO}}_2$ and ${\mathrm{TiO}}_{2-\mathrm{x}}$. The ions drift forms conductive filaments in the oxide film. This shift of the dividing line has a strong impact on the structure conductivity. The geometrical limitation creates two border lines of the possible memristor state. Therefore, one can distinguish two physical states: the low resistance state (LRS) and the high resistance state (HRS). Several memristor elements based on various materials have been proposed after the fabrication of the HP memristor: semiconductor spintronics devices [11], amorphous silicon based devices [12] and memristors in ferroelectric tunnel barriers [13].

In this paper, the investigated elements are the self-directed channel (SDC) memristors [14]. These elements are ion-conducting devices. The active layer of the device is made of chalcogenide material Ge${}_2$Se${}_3$/SnSe/Ag. Unlike the previously described metal–ion devices where the electric field creates conducting filaments between metal electrodes, the SDC device constrains the Ag${}^{+}$ ions movement into the chemically created channels. Permanent conductive channels in the SDC device are generated by inducing an electric field which causes the metal-catalyzed reaction. The device conductivity is determined by the amount of Ag within the channel [15]. The SDC memristor has the active layer doped to enhance and optimize the memristor properties. In this work, we consider the SDC memristors with two kinds of doping agents, i.e., tungsten or carbon. These devices have been purchased from Knowm Inc. Memristors are packaged in an encapsulated edge board with 16 elements and provided with a breakout board, which makes the measurements simpler and more reliable. In general, the behavior of the devices of a given kind is similar, but to keep the consistency of the acquired data, the authors selected the memristor element and all measurements were performed over the same two elements with the given doping agent.

2. Memristors Models

Two existing memristor models are considered in this work: the Strukov model [1] and the VTEAM model [16]. Two window functions are considered for the Strukov model: the rectangular window and the Biolek window [17]. Some improvements of the model have been proposed. It is important to add that the selected models are simple and have a compact structure, which is attractive when the nonlinear optimization algorithms are involved.

2.1. Strukov Model

In the Strukov model, also referred to as the linear ion drift model, it is assumed that the velocity of the ions drift is linearly proportional to the electric field. For this model, the relation between the current i and the voltage v is given by

$$v\left(t\right)=\left(R_{\mathrm{on}}x\left(t\right)+R_{\mathrm{off}}(1-x\left(t\right))\right)i\left(t\right)$$

(1)

where $x\left(t\right)$ represents the relative width of the low-resistance region. It also takes a role of the state variable of the element. The minimal and maximal element resistances are described by $R_{\mathrm{on}}$ and $R_{\mathrm{off}}$ accordingly. They represent the low and high resistance element states. The relative width $x\left(t\right)$ changes with the current $i\left(t\right)$ according to the formula:

$$\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}=k{R}_{\mathrm{on}}i\left(t\right),x\in [0,1]$$

(2)

The relative width $x\left(t\right)$ is the continued function on the interval $[0,1]$. The requirement that $x\in [0,1]$ is fulfilled by multiplying the RHS of the above equation by the rectangular window function

$$f\left(x\right)=\left\{\begin{array}{ccc}\hfill & 1,\hfill & \hfill x\in [0,1],\\ \hfill & 0,\hfill & \hfill x\notin [0,1].\end{array}\right.$$

(3)

In Equation (2), the constant k depends on the properties of the material used and the geometry of the memristor element. The Strukov model depends on three real-valued parameters, i.e., $R_{\mathrm{on}}$, $R_{\mathrm{off}}$ and k.

2.2. Biolek Window Function

Proposed in the previous section, the rectangular window function does not characterize certain physical phenomena in the memristor structure. It is noticed that the dynamics of the model described in Section 2.1 slows down when x is close to zero or one. One of the possible amendments is to update the window function (3). In the literature, several different window functions have been urged [17,18,19]. Here, the Biolek window function [17] is used.

$$f(x,i)=\left\{\begin{array}{ccc}\hfill & \hfill 1-{\left|x-\mathrm{step}(-i)\right|}^p,& \hfill x\in [0,1]\\ \hfill & \hfill 0,& \hfill x\notin [0,1]\end{array}\right.$$

(4)

$\mathrm{step}\left(x\right)=(1+\mathrm{sgn}(x\left)\right)/2$, $\mathrm{sgn}(\times )$ denotes for the sign function, and p is a positive integer. The Biolek window function defined by Equation (4) is modified by introducing the absolute value instead of parentheses around the $x-\mathrm{step}(-i)$ factor [20]. This allows us to consider all integers p instead of only even numbers as in the original.

2.3. Asymmetric Strukov Model

Fabricated memristor nano-devices cannot be considered as symmetric because different switching phenomena are observed under different polarities of applied voltage or the current direction. This forces us to modify the model (2) described in Section 2.1. One of the possible options of model modification is to use different model parameters for different polarities of the input signal. Then the dynamic part of the modified Strukov model can be presented as:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}& =\left\{\begin{array}{cc}{k}_{\mathrm{on}}{R}_{\mathrm{on}}i\left(t\right)f_{\mathrm{on}}(x,i),\hfill & \hfill \mathrm{for}i\ge 0,\\ k_{\mathrm{off}}{R}_{\mathrm{on}}i\left(t\right)f_{\mathrm{off}}(x,i),\hfill & \hfill \mathrm{for}i<0,\end{array}\right.\hfill \end{array}$$

(5)

The parameters $k_{\mathrm{on}}$ and $k_{\mathrm{off}}$ define the memristor dynamics for positive and negative polarity of the current i. Additionally, the modified Strukov model is presented with two different window functions $f_{\mathrm{on}}$ and $f_{\mathrm{off}}$ which also depend on the current direction. In the case of Biolek window, it means that window functions have two additional different p parameters $p_{\mathrm{on}}$ and $p_{\mathrm{off}}$.

Six parameters define the asymmetric Strukov model with Biolek windows: $k_{\mathrm{on}}$, $k_{\mathrm{off}}$, $R_{\mathrm{on}}$, $R_{\mathrm{off}}$, $p_{\mathrm{on}}$ and $p_{\mathrm{off}}$. Four of them are real-valued and two are integer-valued numbers.

2.4. Vteam Model

In the intensive measurement of the real memristors, it has been noticed that there exists a certain current (voltage) threshold below which the device does not change its state variable, i.e., the memristor dynamic disappears. In other words, in order to change the resistance of the memristor, one has to apply the input signal exceeding a certain threshold. The Voltage ThrEshold Adaptive Memristor (VTEAM) model is based on that phenomenon [16]. In this work, author uses the modified VTEAM model proposed in [20]. The internal state variable w is confined to the interval $[0,w_{\mathrm{off}}]$. Then the VTEAM model is defined by:

$$\begin{array}{c}\hfill v\left(t\right)=\frac{R_{\mathrm{on}}\left(w_{\mathrm{off}}-w\left(t\right)\right)+R_{\mathrm{off}}w\left(t\right)}{w_{\mathrm{off}}}i\left(t\right)\end{array}$$

(6)

$$\frac{\mathrm{d}w\left(t\right)}{\mathrm{d}t}=\left\{\begin{array}{cc}{k}_{\mathrm{on}}{\left(\frac{v\left(t\right)}{v_{\mathrm{on}}}-1\right)}^{{\alpha }_{\mathrm{on}}}\hfill & \hfill 0<v_{\mathrm{on}}\le vs.\\ 0,\hfill & \hfill v_{\mathrm{off}}<v<v_{\mathrm{on}}\\ k_{\mathrm{off}}{\left(\frac{v\left(t\right)}{v_{\mathrm{off}}}-1\right)}^{{\alpha }_{\mathrm{off}}}\hfill & \hfill v\le v_{\mathrm{off}}<0\end{array}\right.$$

(7)

If the voltage v belongs to the interval of $(v_{\mathrm{off}},v_{\mathrm{on}})$, then w does not change the threshold phenomenon. Six parameters of the VTEAM model are real-valued ($R_{\mathrm{on}}$, $R_{\mathrm{off}}$, $v_{\mathrm{on}}$, $v_{\mathrm{off}}$, $k_{\mathrm{on}}$, $k_{\mathrm{off}}$) while two of them are integer-valued (${\alpha }_{\mathrm{on}}$ and ${\alpha }_{\mathrm{off}}$). Note that the variable $w_{\mathrm{off}}$ has only a scaling effect on the local resistance, so it can be fixed to a certain value, thus it does not need to be considered in the optimization process. In this work, the authors set $w_{\mathrm{off}}={10}^{-3}$.

3. Experimental

In this section, the memristor measurement setup will be explained. Moreover, the measurement results are discussed and compared.

3.1. Measurement Setup

The devices under study were SDC memristors with two kinds of doping agents, i.e., tungsten and carbon. In order to measure the characteristics of the memristor device it is connected in series with the linear resistor $R_s$. Its main role is to protect the memristor device from burning out by limiting the current. Some manufacturer recommendations for the current limits of the tested devices were applied. Two different series resistors were used. In the case of tungsten dopant, it was $R_s=5.2\mathrm{k}\mathrm{\Omega }$ which limits the current with the applied voltages amplitudes below $1\mathrm{mA}$, and for memristors with carbon dopant the series resistor was $R_s=56.1\mathrm{k}\mathrm{\Omega }$ as the recommended current needs to be lower then 50 $\mathsf{\mu }$A. The voltage $v_r$ measured over the series resistor is proportional to the current, so this measurement permits obtaining the $v-i$ relation of the element.

The setup of the measurement connection with the provided voltage indications used in this work has been presented in Figure 1. Signal $v\left(t\right)$ has been generated by the signal generator device RIGOL DG1020. The response signals $v\left(t\right)$ and $v_r\left(t\right)$ were acquired using a myDAQ oscilloscope manufactured by National Instruments. All acquisition setup was carried out with the LabVIEW software environment.

Before testing procedures, the memristor needs to be formed. Forming entails applying a gradually increasing voltage until the necessary conductive pathways are formed. After the forming procedure, the measurement circuit is excited by the periodic signal $v\left(t\right)$. In this work, the memristor was tested by the sinusoidal voltage wave signal $v\left(t\right)=V_{s}sin\left(2\pi ft\right)$ and the triangle wave signal that can be expressed as $v\left(t\right)=V_t{\int }_0^t\mathrm{sgn}(sin\left(2\pi f\tau \right))\mathrm{d}\tau$. Three different amplitudes of the sinusoidal signal were chosen, i.e, $V_s=0.75 \mathrm{V},V_s=1 \mathrm{V}$ and $V_s=1.5 \mathrm{V}$. In the case of the triangle wave signal, the amplitudes of the signal have been recalculated to keep RMS values equal for both proposed voltage shapes. This allows us to compare the results for different signal shapes but with the same energy dissipation in the memristor over the period of the signal. To fulfill this requirement, the relationship between the signal amplitudes is $V_t=\sqrt{\frac{3}{2}}{V}_s$, thus $V_t=0.92 \mathrm{V},V_t=1.22 \mathrm{V}$ and $V_t=1.83 \mathrm{V}$. Moreover, different frequencies of the exciting signal $v\left(t\right)$ were applied. In this work, the following frequencies of the signals were considered: $f\in \{1,5,10,20,50,100\}\mathrm{Hz}$. Combining all this, the author performed 36 measurements for each kind of device.

3.2. Measurement Results

Applying the setup described in Section 3.1, the data have been acquired. For each combination of excitation voltage amplitude, the shape and frequency author always acquired ${10}^5$ data points, applying the rate of acquisition to the level where 100 periods of the response signal were always measured. For example, for $f=10\mathrm{Hz}$ the rate of acquisition was $10\mathrm{kHz}$. During the measurement, some level of error has been observed. The Savitzky–Golay filter [21] has been applied to remove errors related to high frequencies components in the acquired data. To prepare the reference signal 1000 points from a single period of the output signal have been estimated. The estimation was the average signal over all 100 acquired periods of the signal. In Figure 2, one can find the $i_m=f\left(v_m\right)$ example plot of the measured data compared with the averaged data. The plot represents the acquired data excited with the sinusoidal wave voltage of the amplitude $V_s=1.5\mathrm{V}$ and frequency $f=20\mathrm{Hz}$. A typical hysteresis loop can be observed. One can also note two straight lines crossing the $v-i$ domain at the origin. These lines represent the high and low memristor resistances. Note that the transition is more stable in case of switching to LRS (the low resistance state). In general, the device with a carbon dopant has better stability then with the tungsten doping agent. This is also observed in Figure 3 where first four periods of the acquired data are shown. The averaged data are repeated over each presented period for comparison purposes. Furthermore the excitation data are limited for better visibility of the figure.

3.3. Results Discussion

For understanding the quality of acquired data, some statistical factors have been proposed. To compare the acquired data for each period j (in this work $j\in [1,100]$) with the averaged data the standard deviation normalized over the average value factor ${\delta }_j$ can be formulated as:

$$\begin{array}{c}{\delta }_j=\sqrt{{\displaystyle \frac{k_j}{\eta }}}\hfill \\ {\displaystyle k_j=\sum _{i=1}^n{\left(v_{\mathrm{av}}^i-v_r^{i,j}\right)}^2,\mathrm{and}\eta =\sum _{i=1}^n{\left(v_{\mathrm{av}}^i\right)}^2}\hfill \end{array}$$

(8)

where $n={10}^3$ is the number of data points in one period, and $v_r^{i,j}$ at the j-th period and $v_{av}^i$ are the i-th acquired and averaged voltage data points accordingly ($i\in [1,n]$). To estimate the error over all periods of the sample, another quadratic mean factor $\epsilon$ is proposed and defined as:

$$\epsilon =\sqrt{{\displaystyle \sum _{j=1}^m{\delta }_j^2}}=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{j=1}^m{k}_j}}{m\eta }}}$$

(9)

where m is the number of acquired periods of the measured signal (in this work $m=100$), and $k_j$ and $\eta$ are defined in Equation (8). For example, in Figure 4 the distribution of the deviation factor ${\delta }_j$ defined in (8) over each j-th period is presented. The graphs correspond to the samples presented in Figure 2. For better visibility the $\delta$ axis is in logarithmic scale. The maximum deviation is bigger for tungsten than for the memristor with carbon dopant. The graphs also show the quite big probabilistic behavior of the memristive elements. For the tungsten doped element, the maximum and minimum values of the deviation factor ${\delta }_j$ is 0.106491 and 0.019543 accordingly for 31st and 66th period. In case of carbon doping, these values are 0.070173 and 0.011106 for 61st and 41st period. It can be also noted that carbon doped element is more stable than doped with tungsten which can be measured using the error estimation defined in (9): ${\epsilon }_{\mathrm{T}}=4.3638\times {10}^{-2}$ and ${\epsilon }_{\mathrm{C}}=2.8217\times {10}^{-2}$.

In

Table 1. The values of the error $\epsilon$ defined in Equation (9) of the acquired signals for two tested kinds of memristors and different combinations of voltage amplitude, shapes, and frequencies of the excitation signals. The highest and lowest numbers are indicated with bold font type.

$\mathit{f} \left[\mathbf{Hz}\right]$	Sine			Triangle
	${\mathbf{V}}_{\mathbf{s}}=\mathbf{0}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{0}.\mathbf{92}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{22}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{83}\mathrm{V}$
	tungsten
1	$9.97\times {10}^{-2}$	$7.78\times {10}^{-2}$	$8.59\times {10}^{-2}$	$6.93\times {10}^{-2}$	$\mathbf{1}.\mathbf{69}\times {\mathbf{10}}^{-\mathbf{1}}$	$8.93\times {10}^{-2}$
5	$1.07\times {10}^{-1}$	$5.45\times {10}^{-2}$	$4.97\times {10}^{-2}$	$7.36\times {10}^{-2}$	$4.46\times {10}^{-2}$	$8.35\times {10}^{-2}$
10	$6.11\times {10}^{-2}$	$4.27\times {10}^{-2}$	$3.81\times {10}^{-2}$	$1.10\times {10}^{-1}$	$5.24\times {10}^{-2}$	$1.22\times {10}^{-1}$
20	$6.98\times {10}^{-2}$	$3.85\times {10}^{-2}$	$4.36\times {10}^{-2}$	$5.89\times {10}^{-2}$	$3.93\times {10}^{-2}$	$4.83\times {10}^{-2}$
50	$5.20\times {10}^{-2}$	$2.94\times {10}^{-2}$	$3.51\times {10}^{-2}$	$4.83\times {10}^{-2}$	$1.07\times {10}^{-1}$	$5.64\times {10}^{-2}$
100	$9.08\times {10}^{-2}$	$7.12\times {10}^{-2}$	$8.21\times {10}^{-2}$	$7.03\times {10}^{-2}$	$1.05\times {10}^{-1}$	$1.05\times {10}^{-1}$
	carbon
1	$6.86\times {10}^{-2}$	$4.67\times {10}^{-2}$	$6.71\times {10}^{-2}$	$7.16\times {10}^{-2}$	$1.05\times {10}^{-1}$	$9.34\times {10}^{-2}$
5	$1.01\times {10}^{-1}$	$6.52\times {10}^{-2}$	$6.86\times {10}^{-2}$	$1.00\times {10}^{-1}$	$1.12\times {10}^{-1}$	$6.97\times {10}^{-2}$
10	$6.04\times {10}^{-2}$	$4.54\times {10}^{-2}$	$4.79\times {10}^{-2}$	$1.04\times {10}^{-1}$	$7.89\times {10}^{-2}$	$8.58\times {10}^{-2}$
20	$8.26\times {10}^{-2}$	$\mathbf{2}.\mathbf{42}\times {\mathbf{10}}^{-\mathbf{2}}$	$2.82\times {10}^{-2}$	$8.43\times {10}^{-2}$	$3.81\times {10}^{-2}$	$6.98\times {10}^{-2}$
50	$8.58\times {10}^{-2}$	$6.49\times {10}^{-2}$	$5.50\times {10}^{-2}$	$6.33\times {10}^{-2}$	$5.21\times {10}^{-2}$	$4.37\times {10}^{-2}$
100	$8.31\times {10}^{-2}$	$4.14\times {10}^{-2}$	$4.80\times {10}^{-2}$	$5.99\times {10}^{-2}$	$4.75\times {10}^{-2}$	$4.92\times {10}^{-2}$

, one can find all the errors $\epsilon$ defined in Equation (9) for all possible combinations of the excitation signal. The worst sample is for the memristor with tungsten dopant and for a triangle signal of frequency $f=10\mathrm{Hz}$ and the amplitude $V_t=1.83$ where $\epsilon =0.122$. The best of the acquired signal with $\epsilon =0.0245$ is for the signal measured on memristor with carbon dopant and for sinusoidal excitation of the frequency $f=20\mathrm{Hz}$ and amplitude $V_s=1\mathrm{V}$. Note that in the case of memristor with carbon dopant the higher frequency is, the device responds with better stability.

4. Optimization Procedure

In this section, the optimization procedure is going to be described. The goal of this procedure is to find the optimal fitting parameters of the described in Section 2 memristor models that minimize some cost function. Cost function compares the results obtained from the model simulation with the reference signal. The models considered have parameters either real or integer valued. Let $\mathbf{x}$ and $\mathbf{y}$ be, respectively, the vector of real- and integer-valued parameters of the selected model. The cost function is defined as:

$$F(\mathbf{x},\mathbf{y})=\frac{{\displaystyle \sum _{k=1}^n{\left(i_k(\mathbf{x},\mathbf{y})-i_{\mathrm{ref}}^k\right)}^2}}{{\displaystyle \sum _{k=1}^n{\left(i_{\mathrm{ref}}^k\right)}^2}},$$

(10)

where n is the number of points in the reference signal (here one signal period with $n=1000$ data points), $i_{\mathrm{ref}}^k$ for $k\in \{1,\dots ,n\}$ is the current at $k-\mathrm{th}$ point in the reference signal, and $i_k(\mathbf{x},\mathbf{y})$ is the current at $k-\mathrm{th}$ point computed using the model for the given optimal parameters $\mathbf{x}$ and $\mathbf{y}$. The reference current $i_{\mathrm{ref}}$ is simply computed from the averaged voltage data $v_{\mathrm{av}}$ divided by the linear series resistor $R_s$.

During the optimization process, one needs to find the the optimal parameters $(\mathbf{x},\mathbf{y})$ that minimize the cost function $F(\mathbf{x},\mathbf{y})$. Those optimal parameters represent the best model fit to the acquired real memristor data. Parameters represented as reals and integers are treated differently. Optimal real-valued parameter represented by the vector $\mathbf{x}$ is determined for fixed integer-valued parameter $\mathbf{y}$. This process is carried out for all integer-valued parameters from the proposed set from the fixed range. The vector $\mathbf{x}$ and $\mathbf{y}$ minimizing the cost function Equation (10) is chosen as the optimal solution. The integer-valued parameters range from 1 to 10. Therefore, in the case of the asymmetric Strukov model and the VTEAM model, this means that the optimization process has to be repeated 100 times, which significantly extends the computation time. Similarly to the integer-valued parameters, the real-valued parameters are also subject to some box constrains. As the excitation signal is described by the continuous and periodic function, then the clearly modeled problem needs to be periodic in general. The internal state variables of the model also need to be expressed by the continuous function. To ensure that the continuity condition holds, the process of integration the dynamical part of the model is repeated until the condition $w\left(0\right)=w\left(t_{\mathrm{end}}\right)$ is met, i.e., the transient time is finished and the system works in steady state. This condition requires to estimate an additional real-valued parameter $w_{\mathrm{init}}$ that tells when the integration should start to get the memristor in steady state. All applied constraints to the real-valued parameters are indicated in

Table 2. Box constrains for the real-valued parameters of vector $\mathbf{x}$.

Model	${\mathit{R}}_{\mathbf{on}}$	${\mathit{R}}_{\mathbf{off}}$	${\mathit{k}}_{\mathbf{on}}$, ${\mathit{k}}_{\mathbf{off}}$	${\mathit{v}}_{\mathbf{on}}$, ${\mathit{v}}_{\mathbf{off}}$	x, w
Strukov	$10÷{10}^4$	${10}^3÷5\times {10}^6$	$0.1÷100$	—	$0÷1$
VTEAM	${10}^2÷{10}^4$	${10}^3÷{10}^6$	$0.1÷{10}^3$	$0÷1$	$0÷{10}^{-3}$

.

The interior-point algorithm [22] has been used in the optimization process. This algorithm is implemented as the fmincon function in the MATLAB optimization toolbox. This Matlab function is recommended for the optimization problems dealt in this work, it can handle the box constrains and the situations when the constraints are not valid or cannot be determined. Here the author noticed some problems with the internal variables continuity condition. In some cases, due to the lack of continuity condition, the fulfilling the cost function could not be determined, so this second Matlab function’s property was very useful. In [23], it was shown that the standard explicit Runge–Kutta (4,5) formula, which is implemented in Matlab as ode45 demonstrates some numerical errors while integrating the examined models with some discontinuities. The first-order Euler method was used to simulate memristor models. The Euler method offers reasonably fast integration with controlled step size and can deal with the problems mentioned above.

5. Results

In this section, the results of the optimization process are presented and discussed. The accumulated results are presented in

Table 3. The results of the optimization procedure over all models considered. The maximum and minimum values of the cost function for each model are indicated with bold font type.

Model	$\mathit{f} \left[\mathbf{Hz}\right]$	Tungsten						Carbon
		Sine			Triangle			Sine			Triangle
		${\mathbf{V}}_{\mathbf{s}}=\mathbf{0}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{0}.\mathbf{92}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{22}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{83}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{0}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}\mathrm{V}$	${\mathbf{V}}_{\mathbf{s}}=\mathbf{1}.\mathbf{5}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{0}.\mathbf{92}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{22}\mathrm{V}$	${\mathbf{V}}_{\mathbf{t}}=\mathbf{1}.\mathbf{83}\mathrm{V}$
Strukov	1	$3.61\times {10}^{-2}$	$3.34\times {10}^{-2}$	$2.33\times {10}^{-2}$	$8.46\times {10}^{-3}$	$\mathbf{2}.\mathbf{14}\times {\mathbf{10}}^{-\mathbf{2}}$	$8.90\times {10}^{-3}$	$5.86\times {10}^{-2}$	$5.03\times {10}^{-2}$	$3.47\times {10}^{-2}$	$3.57\times {10}^{-1}$	$3.96\times {10}^{-2}$	$6.83\times {10}^{-2}$
	5	$1.33\times {10}^{-1}$	$1.14\times {10}^{-1}$	$9.78\times {10}^{-2}$	$4.56\times {10}^{-2}$	$5.16\times {10}^{-2}$	$4.31\times {10}^{-2}$	$3.88\times {10}^{-1}$	$3.23\times {10}^{-1}$	$3.66\times {10}^{-1}$	$3.31\times {10}^{-1}$	$3.20\times {10}^{-1}$	$3.92\times {10}^{-1}$
	10	$4.82\times {10}^{-1}$	$4.81\times {10}^{-1}$	$1.51\times {10}^{-1}$	$2.05\times {10}^{-1}$	$2.96\times {10}^{-1}$	$5.87\times {10}^{-2}$	$3.27\times {10}^{-1}$	$2.92\times {10}^{-1}$	$3.49\times {10}^{-1}$	$3.93\times {10}^{-1}$	$3.22\times {10}^{-1}$	$2.82\times {10}^{-1}$
	20	$4.95\times {10}^{-1}$	$\mathbf{5}.\mathbf{38}\times {\mathbf{10}}^{-\mathbf{1}}$	$4.47\times {10}^{-1}$	$3.82\times {10}^{-1}$	$3.26\times {10}^{-1}$	$1.11\times {10}^{-1}$	$4.22\times {10}^{-1}$	$4.36\times {10}^{-1}$	$4.14\times {10}^{-1}$	$3.33\times {10}^{-1}$	$3.11\times {10}^{-1}$	$4.00\times {10}^{-1}$
	50	$4.71\times {10}^{-1}$	$4.02\times {10}^{-1}$	$3.85\times {10}^{-1}$	$3.24\times {10}^{-1}$	$3.62\times {10}^{-1}$	$3.08\times {10}^{-1}$	$2.98\times {10}^{-1}$	$3.59\times {10}^{-1}$	$3.31\times {10}^{-1}$	$1.65\times {10}^{-1}$	$1.64\times {10}^{-1}$	$1.59\times {10}^{-1}$
	100	$3.87\times {10}^{-1}$	$3.33\times {10}^{-1}$	$3.18\times {10}^{-1}$	$3.70\times {10}^{-1}$	$3.40\times {10}^{-1}$	$2.92\times {10}^{-1}$	$3.90\times {10}^{-1}$	$3.75\times {10}^{-1}$	$3.20\times {10}^{-1}$	$1.70\times {10}^{-1}$	$1.24\times {10}^{-1}$	$1.38\times {10}^{-1}$
Strukov asymmetric	1	$4.94\times {10}^{-1}$	$4.78\times {10}^{-1}$	$4.82\times {10}^{-1}$	$2.19\times {10}^{-1}$	$\mathbf{5}.\mathbf{86}\times {\mathbf{10}}^{-\mathbf{3}}$	$4.06\times {10}^{-2}$	$2.76\times {10}^{-1}$	$3.46\times {10}^{-1}$	$4.09\times {10}^{-1}$	$2.12\times {10}^{-1}$	$1.33\times {10}^{-1}$	$1.15\times {10}^{-1}$
	5	$4.77\times {10}^{-1}$	$4.82\times {10}^{-1}$	$4.82\times {10}^{-1}$	$3.04\times {10}^{-1}$	$3.17\times {10}^{-1}$	$3.01\times {10}^{-1}$	$3.96\times {10}^{-1}$	$3.25\times {10}^{-1}$	$3.68\times {10}^{-1}$	$3.64\times {10}^{-1}$	$3.43\times {10}^{-1}$	$3.95\times {10}^{-1}$
	10	$4.84\times {10}^{-1}$	$4.68\times {10}^{-1}$	$4.71\times {10}^{-1}$	$2.06\times {10}^{-1}$	$2.96\times {10}^{-1}$	$2.38\times {10}^{-1}$	$3.23\times {10}^{-1}$	$2.95\times {10}^{-1}$	$3.50\times {10}^{-1}$	$4.06\times {10}^{-1}$	$3.24\times {10}^{-1}$	$3.53\times {10}^{-1}$
	20	$4.95\times {10}^{-1}$	$\mathbf{5}.\mathbf{37}\times {\mathbf{10}}^{-\mathbf{1}}$	$4.47\times {10}^{-1}$	$3.83\times {10}^{-1}$	$3.28\times {10}^{-1}$	$2.52\times {10}^{-1}$	$4.25\times {10}^{-1}$	$4.38\times {10}^{-1}$	$4.13\times {10}^{-1}$	$3.42\times {10}^{-1}$	$3.14\times {10}^{-1}$	$3.97\times {10}^{-1}$
	50	$4.72\times {10}^{-1}$	$4.02\times {10}^{-1}$	$3.85\times {10}^{-1}$	$3.25\times {10}^{-1}$	$3.64\times {10}^{-1}$	$3.05\times {10}^{-1}$	$2.91\times {10}^{-1}$	$3.60\times {10}^{-1}$	$3.35\times {10}^{-1}$	$1.60\times {10}^{-1}$	$1.68\times {10}^{-1}$	$1.66\times {10}^{-1}$
	100	$3.90\times {10}^{-1}$	$3.34\times {10}^{-1}$	$3.19\times {10}^{-1}$	$3.71\times {10}^{-1}$	$3.42\times {10}^{-1}$	$2.94\times {10}^{-1}$	$3.85\times {10}^{-1}$	$3.76\times {10}^{-1}$	$3.20\times {10}^{-1}$	$1.68\times {10}^{-1}$	$1.27\times {10}^{-1}$	$1.44\times {10}^{-1}$
Biolek window	1	$3.65\times {10}^{-2}$	$3.36\times {10}^{-2}$	$2.35\times {10}^{-2}$	$\mathbf{7}.\mathbf{08}\times {\mathbf{10}}^{-\mathbf{3}}$	$2.11\times {10}^{-2}$	$8.69\times {10}^{-3}$	$6.68\times {10}^{-2}$	$4.95\times {10}^{-2}$	$3.50\times {10}^{-2}$	$8.01\times {10}^{-2}$	$3.92\times {10}^{-2}$	$6.75\times {10}^{-2}$
	5	$1.34\times {10}^{-1}$	$1.16\times {10}^{-1}$	$9.90\times {10}^{-2}$	$4.63\times {10}^{-2}$	$5.23\times {10}^{-2}$	$4.37\times {10}^{-2}$	$3.86\times {10}^{-1}$	$3.23\times {10}^{-1}$	$8.50\times {10}^{-1}$	$3.33\times {10}^{-1}$	$3.38\times {10}^{-1}$	$1.03\times {10}^{-1}$
	10	$1.76\times {10}^{-1}$	$1.61\times {10}^{-1}$	$1.52\times {10}^{-1}$	$6.37\times {10}^{-2}$	$9.24\times {10}^{-2}$	$5.99\times {10}^{-2}$	$3.26\times {10}^{-1}$	$2.92\times {10}^{-1}$	$3.49\times {10}^{-1}$	$3.72\times {10}^{-1}$	$2.40\times {10}^{-1}$	$5.88\times {10}^{-1}$
	20	$4.95\times {10}^{-1}$	$5.38\times {10}^{-1}$	$2.10\times {10}^{-1}$	$2.53\times {10}^{-1}$	$2.09\times {10}^{-1}$	$1.18\times {10}^{-1}$	$\mathbf{6}.\mathbf{67}\times {\mathbf{10}}^{-\mathbf{1}}$	$4.36\times {10}^{-1}$	$4.15\times {10}^{-1}$	$3.35\times {10}^{-1}$	$3.05\times {10}^{-1}$	$3.96\times {10}^{-1}$
	50	$4.71\times {10}^{-1}$	$4.05\times {10}^{-1}$	$3.15\times {10}^{-1}$	$3.28\times {10}^{-1}$	$3.82\times {10}^{-1}$	$3.06\times {10}^{-1}$	$6.06\times {10}^{-1}$	$6.06\times {10}^{-1}$	$3.34\times {10}^{-1}$	$1.65\times {10}^{-1}$	$1.69\times {10}^{-1}$	$1.61\times {10}^{-1}$
	100	$3.87\times {10}^{-1}$	$3.40\times {10}^{-1}$	$3.63\times {10}^{-1}$	$3.71\times {10}^{-1}$	$3.51\times {10}^{-1}$	$3.55\times {10}^{-1}$	$3.88\times {10}^{-1}$	$3.74\times {10}^{-1}$	$3.20\times {10}^{-1}$	$1.70\times {10}^{-1}$	$1.27\times {10}^{-1}$	$1.42\times {10}^{-1}$
Biolek asymmetric	1	$1.42\times {10}^{-1}$	$1.22\times {10}^{-1}$	$9.06\times {10}^{-2}$	$1.96\times {10}^{-2}$	$\mathbf{9}.\mathbf{83}\times {\mathbf{10}}^{-\mathbf{3}}$	$4.12\times {10}^{-2}$	$5.79\times {10}^{-2}$	$6.17\times {10}^{-2}$	$6.44\times {10}^{-2}$	$8.82\times {10}^{-2}$	$3.84\times {10}^{-2}$	$4.93\times {10}^{-2}$
	5	$2.10\times {10}^{-1}$	$2.72\times {10}^{-1}$	$1.49\times {10}^{-1}$	$7.74\times {10}^{-2}$	$9.13\times {10}^{-2}$	$6.35\times {10}^{-2}$	$2.18\times {10}^{-1}$	$1.49\times {10}^{-1}$	$1.20\times {10}^{-1}$	$1.87\times {10}^{-1}$	$1.34\times {10}^{-1}$	$1.18\times {10}^{-1}$
	10	$4.47\times {10}^{-1}$	$3.73\times {10}^{-1}$	$2.29\times {10}^{-1}$	$9.19\times {10}^{-2}$	$1.77\times {10}^{-1}$	$8.29\times {10}^{-2}$	$2.66\times {10}^{-1}$	$2.77\times {10}^{-1}$	$2.54\times {10}^{-1}$	$2.93\times {10}^{-1}$	$2.09\times {10}^{-1}$	$1.28\times {10}^{-1}$
	20	$4.94\times {10}^{-1}$	$\mathbf{5}.\mathbf{37}\times {\mathbf{10}}^{-\mathbf{1}}$	$2.87\times {10}^{-1}$	$3.67\times {10}^{-1}$	$2.21\times {10}^{-1}$	$1.25\times {10}^{-1}$	$4.01\times {10}^{-1}$	$4.16\times {10}^{-1}$	$4.10\times {10}^{-1}$	$3.32\times {10}^{-1}$	$2.85\times {10}^{-1}$	$3.28\times {10}^{-1}$
	50	$4.72\times {10}^{-1}$	$4.02\times {10}^{-1}$	$3.85\times {10}^{-1}$	$3.25\times {10}^{-1}$	$3.63\times {10}^{-1}$	$3.05\times {10}^{-1}$	$2.90\times {10}^{-1}$	$3.53\times {10}^{-1}$	$3.06\times {10}^{-1}$	$1.57\times {10}^{-1}$	$1.55\times {10}^{-1}$	$1.40\times {10}^{-1}$
	100	$3.89\times {10}^{-1}$	$3.34\times {10}^{-1}$	$3.19\times {10}^{-1}$	$3.70\times {10}^{-1}$	$3.41\times {10}^{-1}$	$2.94\times {10}^{-1}$	$3.82\times {10}^{-1}$	$3.70\times {10}^{-1}$	$3.12\times {10}^{-1}$	$1.63\times {10}^{-1}$	$1.21\times {10}^{-1}$	$1.37\times {10}^{-1}$
VTEAM	1	$2.02\times {10}^{-3}$	$1.56\times {10}^{-3}$	$1.44\times {10}^{-3}$	$3.28\times {10}^{-3}$	$2.18\times {10}^{-2}$	$5.12\times {10}^{-3}$	$1.62\times {10}^{-3}$	$3.88\times {10}^{-3}$	$3.20\times {10}^{-3}$	$1.16\times {10}^{-2}$	$3.00\times {10}^{-3}$	$4.62\times {10}^{-3}$
	5	$2.81\times {10}^{-3}$	$8.93\times {10}^{-4}$	$1.13\times {10}^{-3}$	$3.97\times {10}^{-3}$	$4.45\times {10}^{-3}$	$1.86\times {10}^{-2}$	$1.36\times {10}^{-3}$	$3.31\times {10}^{-3}$	$3.24\times {10}^{-3}$	$1.38\times {10}^{-2}$	$6.17\times {10}^{-3}$	$2.55\times {10}^{-3}$
	10	$3.02\times {10}^{-3}$	$1.49\times {10}^{-3}$	$\mathbf{7}.\mathbf{96}\times {\mathbf{10}}^{-\mathbf{4}}$	$4.29\times {10}^{-3}$	$4.88\times {10}^{-3}$	$1.45\times {10}^{-2}$	$3.12\times {10}^{-3}$	$2.57\times {10}^{-3}$	$3.83\times {10}^{-3}$	$3.61\times {10}^{-3}$	$5.31\times {10}^{-3}$	$6.75\times {10}^{-3}$
	20	$4.20\times {10}^{-2}$	$1.29\times {10}^{-1}$	$1.10\times {10}^{-3}$	$1.53\times {10}^{-2}$	$4.44\times {10}^{-3}$	$7.33\times {10}^{-3}$	$5.23\times {10}^{-3}$	$4.33\times {10}^{-3}$	$3.36\times {10}^{-3}$	$8.00\times {10}^{-3}$	$1.12\times {10}^{-2}$	$4.53\times {10}^{-3}$
	50	$1.69\times {10}^{-3}$	$5.51\times {10}^{-3}$	$8.19\times {10}^{-3}$	$1.02\times {10}^{-2}$	$6.82\times {10}^{-3}$	$6.45\times {10}^{-3}$	$8.63\times {10}^{-3}$	$4.12\times {10}^{-3}$	$3.27\times {10}^{-3}$	$5.77\times {10}^{-3}$	$6.78\times {10}^{-3}$	$9.19\times {10}^{-3}$
	100	$7.25\times {10}^{-3}$	$3.41\times {10}^{-3}$	$\mathbf{2}.\mathbf{34}\times {\mathbf{10}}^{-\mathbf{2}}$	$4.17\times {10}^{-3}$	$1.41\times {10}^{-2}$	$2.04\times {10}^{-2}$	$1.25\times {10}^{-3}$	$4.00\times {10}^{-3}$	$3.36\times {10}^{-3}$	$5.57\times {10}^{-3}$	$6.14\times {10}^{-3}$	$7.78\times {10}^{-3}$

. The best and worst modeled memristors for each model are indicated with the pink background. The memristor models are nonlinear, so the cost function (10) is also strongly nonlinear. Obtained results usually are rather local and it would be very hard to obtain the global minimum. Thus, the performance of the optimization procedure strongly depends on the initial point ${\mathbf{x}}_{\mathrm{init}}$ in the possible solution space provided to the algorithm. The initial points are prepared by hand performing some optimization tests before and then fixed for one model with all combinations of measurement data.

5.1. Strukov Model, Rectangular Window

Comparison of the signal obtained from the optimal Strukov model with the rectangular window and the reference signal is presented in Figure 5 and Figure 6. The asymmetric Strukov model with rectangular window outperforms the symmetric Strukov model. The cost function value is $3.65$ times lower in case of the asymmetric model. The fitting parameters for the symmetric models are: $R_{\mathrm{on}}=929.5\mathrm{\Omega }$, $R_{\mathrm{off}}=17.99\mathrm{k}\mathrm{\Omega }$, $k_{\mathrm{on}}=499.97$ and $x_{\mathrm{init}}=0.00016$. For the asymmetric model, the parameters are: $R_{\mathrm{on}}=818.5\mathrm{\Omega }$, $R_{\mathrm{off}}=477.86\mathrm{k}\mathrm{\Omega }$, $k_{\mathrm{on}}=97.23$, $k_{\mathrm{off}}=4.736$ and $x_{\mathrm{init}}=0.938754$. The presented data for the best model are obtained using the Strukov model. From Table 3, one can notice that while the frequency increases, the model performance decreases and the proposed model is not sufficiently reflecting the real data.

5.2. Strukov Model, Biolek Window

Comparison of the signal obtained from the optimal Strukov model with the Biolek window and the reference signal is presented is Figure 7. The reference data are obtained for a carbon doped memristor with the sine wave of the applied excitation voltage with the frequency $f=1 \mathrm{Hz}$ and the amplitude $V_s=0.75 \mathrm{V}$. The obtained fitting parameters are: $R_{\mathrm{on}}=9.99\mathrm{k}\mathrm{\Omega }$, $R_{\mathrm{off}}=258.24\mathrm{k}\mathrm{\Omega }$, $k_{\mathrm{on}}=499.99$, $p_{\mathrm{on}}=5$ and $x_{\mathrm{init}}=0.00017$. This result can be compared with the asymmetric model which is presented in Figure 8 where the obtained fitting parameters are: $R_{\mathrm{on}}=9.99\mathrm{k}\mathrm{\Omega }$, $R_{\mathrm{off}}=4.99\mathrm{M}\mathrm{\Omega }$, $k_{\mathrm{on}}=99.99$, $k_{\mathrm{off}}=7.64$, $p_{\mathrm{on}}=10$, $p_{\mathrm{off}}=8$ and $x_{\mathrm{init}}=0.00023$. The cost function value for the asymmetric model is 15% lower then for the symmetric model. Unfortunately, as for Strukov with a rectangular window, it can be observed from Table 3 that when the frequency of the signal increases, the model performance decreases. Additionally, it has to be pointed out that in general, the asymmetric model is not always better than the symmetric one and the fitting cost function values are on similar levels.

5.3. VTEAM

Figure 9 presents the comparison plot of the solution using the VTEAM model and the reference data. The fitting parameters of the model are: $R_{\mathrm{on}}=1.3578\mathrm{k}\mathrm{\Omega }$, $R_{\mathrm{off}}=352.196\mathrm{k}\mathrm{\Omega }$, $k_{\mathrm{on}}=-997.642$, $k_{\mathrm{off}}=0.2969$, ${\alpha }_{\mathrm{on}}=2$, ${\alpha }_{\mathrm{off}}=1$, $v_{\mathrm{on}}=-0.340661\mathrm{V}$, $v_{\mathrm{off}}=0.044669\mathrm{V}$ and $w_{\mathrm{init}}=0.00099943$. In this case the value of the cost function $F(\mathbf{x},\mathbf{y})=3.3647\times {10}^{-3}$ is significantly lower than the result from the rest of the considered models. In general, the VTEAM model fits much better then the other models considered in this work. One needs always to remember that the solution might be better if different initial points are provided to the optimization algorithm.

5.4. Model Performance with Varying Frequency

In this subsection, the behavior of the VTEAM model is studied when frequencies of the input signal are varied. The authors consider the VTEAM model only as it outperforms the other considered models. The optimal model obtained for a given frequency is applied to other frequencies. The results are presented in

Table 4. Cost function values if the optimal model parameters obtained for certain frequency are applied for the data of different frequency. The VTEAM model is considered.

Wave	${\mathit{V}}_{\mathit{s},\mathit{t}}\left[\mathbf{V}\right]$	$\mathit{f} \left[\mathbf{Hz}\right]$	Tungsten						Carbon
			1	5	10	20	50	100	1	5	10	20	50	100
sine	$0.5$	1	$0.00202$	$0.00528$	$0.08234$	$0.19195$	$0.05030$	$0.07946$	$0.00162$	$0.05078$	$0.12275$	$0.32129$	$0.36847$	$0.72049$
		5	$0.00410$	$0.00281$	$0.06955$	$0.14980$	$0.02493$	$0.04155$	$0.02232$	$0.00136$	$0.02233$	$0.07843$	$0.07436$	$0.27856$
		10	$0.04493$	$0.04796$	$0.00302$	$0.10312$	$0.01155$	$0.01302$	$0.02721$	$0.01354$	$0.00312$	$0.02527$	$0.00702$	$0.01995$
		20	$0.06857$	$0.05963$	$0.02514$	$0.04199$	$0.03364$	$0.03270$	$0.04361$	$0.01868$	$0.01732$	$0.00523$	$0.02024$	$0.01433$
		50	$0.00898$	$0.01828$	$0.01407$	$0.11139$	$0.00169$	$0.00959$	$0.01651$	$0.01581$	$0.01206$	$0.04643$	$0.00863$	$0.01941$
		100	$0.02860$	$0.02529$	$0.00882$	$0.08274$	$0.00382$	$0.00725$	$0.02185$	$0.01301$	$0.01235$	$0.01588$	$0.01868$	$0.00125$
	1	1	$0.00156$	$0.00685$	$0.03065$	$0.23832$	$0.05171$	$0.07386$	$0.00388$	$0.02676$	$0.07722$	$0.41649$	$0.56214$	$0.73261$
		5	$0.00362$	$0.00089$	$0.01033$	$0.21392$	$0.01445$	$0.02901$	$0.00594$	$0.00331$	$0.01184$	$0.09024$	$0.15026$	$0.34278$
		10	$0.00560$	$0.00657$	$0.00149$	$0.21238$	$0.00597$	$0.02043$	$0.00966$	$0.00777$	$0.00257$	$0.04661$	$0.02910$	$0.07968$
		20	$0.03739$	$0.02961$	$0.01599$	$0.12908$	$0.02403$	$0.04214$	$0.01728$	$0.02449$	$0.02891$	$0.00433$	$0.01129$	$0.00841$
		50	$0.01125$	$0.00479$	$0.00204$	$0.19447$	$0.00551$	$0.01977$	$0.01803$	$0.01185$	$0.01215$	$0.01192$	$0.00412$	$0.01229$
		100	$0.02749$	$0.02485$	$0.02543$	$0.24276$	$0.00925$	$0.00341$	$0.12617$	$0.01378$	$0.01544$	$0.00913$	$0.00739$	$0.00400$
	$1.5$	1	$0.00144$	$0.00263$	$0.00386$	$0.01026$	$0.01107$	$0.02769$	$0.00320$	$0.03052$	$0.09668$	$0.29213$	$0.42642$	$0.62507$
		5	$0.00260$	$0.00113$	$0.00170$	$0.00532$	$0.01276$	$0.03036$	$0.00770$	$0.00324$	$0.02066$	$0.05778$	$0.07646$	$0.20810$
		10	$0.00408$	$0.00160$	$0.00080$	$0.00305$	$0.00854$	$0.02597$	$0.01458$	$0.01297$	$0.00383$	$0.01361$	$0.00763$	$0.04172$
		20	$0.00848$	$0.00377$	$0.00215$	$0.00110$	$0.01102$	$0.02595$	$0.01340$	$0.01931$	$0.01030$	$0.00336$	$0.01710$	$0.01607$
		50	$0.00365$	$0.00145$	$0.00111$	$0.00320$	$0.00819$	$0.02508$	$1.30218$	$0.00837$	$0.00669$	$0.02131$	$0.00327$	$0.02399$
		100	$0.00491$	$0.00163$	$0.00107$	$0.00266$	$0.01365$	$0.02344$	$0.01842$	$0.01884$	$0.00744$	$0.01319$	$0.01605$	$0.00336$
triangle	$0.92$	1	$0.00328$	$0.02756$	$0.01345$	$0.06863$	$0.03902$	$0.05916$	$0.01162$	$0.05320$	$0.15727$	$0.36263$	$0.24198$	$0.33194$
		5	$0.02250$	$0.00397$	$0.01784$	$0.01500$	$0.01453$	$0.03309$	$0.02694$	$0.01375$	$0.07155$	$0.12952$	$0.09625$	$0.17665$
		10	$0.01124$	$0.01798$	$0.00429$	$0.05015$	$0.03001$	$0.05441$	$0.02624$	$0.08068$	$0.00361$	$0.07752$	$0.09537$	$0.10014$
		20	$0.07781$	$0.03306$	$0.08420$	$0.01526$	$0.02954$	$0.02089$	$0.02373$	$0.06826$	$0.05062$	$0.00800$	$0.04952$	$0.03280$
		50	$0.02490$	$0.00812$	$0.02863$	$0.01645$	$0.01025$	$0.01696$	$0.08447$	$0.13602$	$0.10410$	$0.07456$	$0.00577$	$0.00836$
		100	$0.04533$	$0.01032$	$0.03942$	$0.00420$	$0.00969$	$0.00417$	$0.07098$	$0.12083$	$0.09650$	$0.06306$	$0.00652$	$0.00557$
	$1.22$	1	$0.02379$	$0.29691$	$0.51993$	$0.78530$	$0.92088$	$0.82673$	$0.00300$	$0.12867$	$0.14644$	$0.45204$	$0.31136$	$0.29323$
		5	$0.03804$	$0.00445$	$0.00537$	$0.00629$	$0.00685$	$0.00663$	$0.05600$	$0.00617$	$0.05412$	$0.13171$	$0.10372$	$0.11275$
		10	$0.03668$	$0.00505$	$0.00488$	$0.00695$	$0.00983$	$0.00996$	$0.01944$	$0.06724$	$0.00531$	$0.08045$	$0.07953$	$0.08451$
		20	$0.04579$	$0.00585$	$0.00698$	$0.00444$	$0.00540$	$0.00904$	$0.03641$	$0.03867$	$0.03776$	$0.01121$	$0.05271$	$0.06869$
		50	$0.04583$	$0.00935$	$0.01256$	$0.00853$	$0.00682$	$0.01288$	$0.06496$	$0.12495$	$0.05243$	$0.06924$	$0.00678$	$0.01374$
		100	$0.07679$	$0.02742$	$0.03541$	$0.02772$	$0.01491$	$0.01411$	$0.10404$	$0.17346$	$0.08817$	$0.10400$	$0.01124$	$0.00614$
	$1.83$	1	$0.00512$	$0.00533$	$0.01501$	$0.00838$	$0.00710$	$0.00831$	$0.00462$	$0.22123$	$0.22634$	$0.81042$	$0.32819$	$0.30914$
		5	$0.02298$	$0.01858$	$0.03450$	$0.03431$	$0.05711$	$0.11863$	$0.06020$	$0.00255$	$0.02584$	$0.04995$	$0.03050$	$0.04747$
		10	$0.03422$	$0.02305$	$0.01449$	$0.06873$	$0.29556$	$0.41337$	$0.04828$	$0.03286$	$0.00675$	$0.10665$	$0.04991$	$0.06241$
		20	$0.00763$	$0.00946$	$0.01384$	$0.00733$	$0.01075$	$0.01193$	$0.06889$	$0.01262$	$0.06221$	$0.00453$	$0.09330$	$0.08848$
		50	$0.00633$	$0.00536$	$0.01834$	$0.00960$	$0.00645$	$0.00886$	$0.09882$	$0.06801$	$0.03454$	$0.10145$	$0.00919$	$0.01805$
		100	$0.03956$	$0.03337$	$0.06507$	$0.05181$	$0.03147$	$0.02044$	$0.13061$	$0.10931$	$0.05451$	$0.13825$	$0.01409$	$0.00778$

. The model is optimized for the fixed frequency f reported in the third column of Table 4. The values of the cost function obtained for other frequencies are presented in the following columns. All this is presented for all combinations of the excitation voltage and the type of memristor dopant. In general, it can be noticed that the cost function values for the frequency combinations laying outside the diagonals are higher. It is interesting that there are some exceptions, e.g., for the tungsten doped memristor and with triangle excitation voltage of amplitude $V_t=1.83\mathrm{V}$ the performance of the model with the parameters obtained for frequency $100\mathrm{Hz}$ and applied to the data of the frequencies $1\mathrm{Hz}$, $20\mathrm{Hz}$ and $50\mathrm{Hz}$ performs even better then for original one. There are also some worse scenarios. For example, the carbon doped memristor with sine voltage function and the parameters obtained for frequency $100\mathrm{Hz}$ applied to $1\mathrm{Hz}$ produce almost a 600 times worse fit compared to the reference data.

6. Conclusions

The main goal of this work was to check the performance of the memristor models. As possible memristor applications are very wide, scientists and engineers expect a model that can provide some level of certainty independent of the excitation conditions. In the modeled devices, one can observe the bistability phenomenon. The problem is the transient state when the memristor switches. From the collected data, one can observe that switching from the HRS to LRS is more stable, but when the state changes in the opposite direction, one can find a much bigger stochastic behavior. Memristors with a carbon doping agent seem to be more stable than the device doped with tungsten. However, it is an interesting fact that the stability is on the higher level in comparison with the Pt/HfO${}_2$/Ti/Pt memristor [24]. This uncertain behavior is modeled using some stochastic models [25,26,27,28]. There are some interesting approaches for the stochastic models adding the additional noise source [29,30,31]. It was also noticed that the considered noise and periodic signals can generate much more interesting and complex phenomena [32]. This manifests as an additional memristor property of stochastic multistable nonlinear systems. In the author’s opinion, this is interesting, but on the other hand this approach can model the element only with the certain fixed data measurement. Actually, the electronic engineers dispose of the appropriate models for standard electronic devices like diodes, transistors, operational amplifiers etc. The main goal of the model is to be as general as possible and it should fit to the modeled device under a big variety of excitation signals or even physical conditions. That is why in this work some memristor models have been tested. The response of the model parameters obtained for one frequency to the other frequency has also been analyzed. The result looks to be promising. In some cases, the performance is weak but there are also some contrary positive examples.

In this work, it has been shown that among all memristor models considered, the VTEAM model significantly outperforms the other ones. This conclusion can be seen analyzing the data presented in

Table 5. Average cost function values for all tested models grouped by the doping agent and the wave shape of the excitation voltage.

Model	Tungsten		Carbon
	Sine	Triangle	Sine	Triangle
Strukov	$0.29595$	$0.19754$	$0.30751$	$0.24822$
Strukov asymmetric	$0.44997$	$0.27185$	$0.35735$	$0.25717$
Biolek window	$0.24764$	$0.17035$	$0.48473$	$0.22962$
Biolek window asymmetric	$0.31406$	$0.18703$	$0.26208$	$0.17020$
VTEAM	$0.01316$	$0.00956$	$0.00354$	$0.00680$

where the average cost function values of the models are presented grouped by the excitation voltage, wave shape, and the doping agent. The threshold idea was the right choice to model the memristor behavior; on the other hand, the proposed asymmetric models do not increase significantly the model performance. It can be noticed that in the case of Strukov model with Biolek window, the symmetric model is even a little better then the asymmetric one for the case of memristor with tungsten dopant. Moreover, it has to be mentioned that in general, the Strukov model performance is limited to the frequency in the low range.

In the future work, authors can introduce some amendments to the tested models and apply the multi-objective optimization approach to improve their performance over a wider frequency range. The improvement of the measurement set-up is going to be introduced to avoid some data acquisition errors, especially for low currents. The author is going to implement the achieved results to test the memory devices as multi-bit memory cells.

The total results of the models with their fitting parameters are not presented in this work for paper clarity, but they can be provided by the author on any request.