micromachines Article Fast Fitting of the Dynamic Memdiode Model to the Conduction Characteristics of RRAM Devices Using Convolutional Neural Networks Fernando Leonel Aguirre 1, * , Eszter Piros 2 , Nico Kaiser 2 , Tobias Vogel 2 , Stephan Petzold 2 , Jonas Gehrunger 3 , Timo Oster 4 , Christian Hochberger 3 , Jordi Suñé 1 , Lambert Alff 2 and Enrique Miranda 1, * 1 Departament d’Enginyeria Electrònica, Universitat Autònoma de Barcelona, 08193 Cerdanyola del Vallès, Spain 2 Advanced Thin Film Technology Division, Institute of Materials Science, Technische Universität Darmstadt, 64289 Darmstadt, Germany 3 Computer Systems Group, Department of Electrical and Information Engineering, Technische Universität Darmstadt, 64289 Darmstadt, Germany 4 Integrated Electronic Systems, Department of Electrical and Information Engineering, Technische Universität Darmstadt, 64289 Darmstadt, Germany * Correspondence: fernando.aguirre@uab.cat (F.L.A.); enrique.miranda@uab.cat (E.M.) Abstract: In this paper, the use of Artificial Neural Networks (ANNs) in the form of Convolutional Neural Networks (AlexNET) for the fast and energy-efficient fitting of the Dynamic Memdiode Model (DMM) to the conduction characteristics of bipolar-type resistive switching (RS) devices is investigated. Despite an initial computationally intensive training phase the ANNs allow obtaining a mapping between the experimental Current-Voltage (I-V) curve and the corresponding DMM Citation: Aguirre, F.L.; Piros, E.; parameters without incurring a costly iterative process as typically considered in error minimization- Kaiser, N.; Vogel, T.; Petzold, S.; Gehrunger, J.; Oster, T.; Hochberger, based optimization algorithms. In order to demonstrate the fitting capabilities of the proposed C.; Suñé, J.; Alff, L.; et al. Fast Fitting approach, a complete set of I-Vs obtained from Y2 O3 -based RRAM devices, fabricated with different of the Dynamic Memdiode Model to oxidation conditions and measured with different current compliances, is considered. In this way, in the Conduction Characteristics of addition to the intrinsic RS variability, extrinsic variation is achieved by means of external factors RRAM Devices Using Convolutional (oxygen content and damage control during the set process). We show that the reported method Neural Networks. Micromachines provides a significant reduction of the fitting time (one order of magnitude), especially in the case of 2022, 13, 2002. https://doi.org/ large data sets. This issue is crucial when the extraction of the model parameters and their statistical 10.3390/mi13112002 characterization are required. Academic Editors: Andrey Sokolov and Haider Abbas Keywords: RRAM; neural networks; curve fitting; dynamic memdiode; memristor Received: 1 October 2022 Accepted: 8 November 2022 Published: 17 November 2022 1. Introduction Publisher’s Note: MDPI stays neutral Since the first practical description of a memristive device by HP in 2008 [1], a number with regard to jurisdictional claims in of compact models for the current-voltage (I-V) characteristic of resistive switching (RS) published maps and institutional affil- devices have been proposed [2]. However, selecting an appropriate generic model for an iations. electron device is far from being a simple and direct task. The model should not only be able to cover the basic features experimentally observed, but also the distinctive details of the device under study. The accurate representation of the electron transport mechanism Copyright: © 2022 by the authors. in the investigated device encourages the design and simulation of more complex circuits Licensee MDPI, Basel, Switzerland. and/or systems and allows the identification and organization of the elementary pieces This article is an open access article that lead to the variety of observed behaviors. For circuit simulation-oriented models, this distributed under the terms and capacity of adaptation must be achieved by means of a reduced number of parameters conditions of the Creative Commons and by well-behaved and selected equations. This is the signature of a compact behavioral Attribution (CC BY) license (https:// approach: the emphasis is on the phenomenon representation rather than the physical creativecommons.org/licenses/by/ details behind it. This does not mean that physical aspects are completely abandoned, 4.0/). but they must be conjugated with suitable ad hoc concepts and appropriate mathematical Micromachines 2022, 13, 2002. https://doi.org/10.3390/mi13112002 https://www.mdpi.com/journal/micromachines Micromachines 2022, 13, 2002 2 of 14 tools. Among this kind of behavioral approaches, the Dynamic Memdiode Model (DMM) proposed by our group in Ref. [2] has been proven capable not only of simulating isolated RRAM devices, but also simulating large RRAM-based Artificial Neural Networks (ANN) with parasitic elements used for pattern recognition tasks [3,4]. Nonetheless, the fitting of parametrized functional forms, such as those involved in the DMM, to sets of experimental data (curve fitting) is a well-known problem in data analysis. The optimal parameter values are conventionally found by minimizing an error measure, often taken to be the sum of the squares of the errors between the observed data values and those predicted by the function. If the functional form is linearly dependent on the parameters (e.g., a polynomial), then the minimization problem is linear and can be easily solved. However, in many cases, it is necessary to consider functional forms which nonlinearly depend on the unknown parameters. The error minimization procedure in such cases generally involves an iterative algorithm starting from an initial guess. Depending on the number of parameters, such iterative methods can be computationally intensive and hence, slow and for complex problems, the selection of a suitable initial guess can require human intervention to ensure convergence to the correct solution. This is the case of the Levenberg–Marquardt (LM) [5–7] and Genetic (GA) algorithms [8–11] which have been proposed to automate the model parameter calibration (extraction). While the LM method tends to be stuck in a local optimum when launched without a proper initial guess, the GA is a global optimizer for combinatorial optimization problems, but it is not recommended for tuning continuous parameters. In addition, the Pattern Search (PS) algorithm (one of the so-called direct searching algorithms) [12–14] has been considered, as it is a very simple and fast local searching algorithm since it is derivative-free and depends only on simple calculations. Although this method has the capability to step over the hillocks in the parameter space to escape the local optima, due to the independence from derivation in the pattern vector, it loses information to a certain extent in high dimensional problems since it does not utilize any information gained along the optimization trajectory. As a result, for applications involving high volumes of data or for real-time applications, there is considerable interest in techniques which can automate the curve fitting process and operate at high speed [15] without incurring too much energy consumption. In this paper, we consider the use of convolutional neural networks (CNN) as suitable and fast alternatives to determine the optimal DMM parameter values directly from the raw data [15]. ANNs, such as the CNNs, are widely used in many disciplines such as system identification, control, pattern recognition, gaming, translation, medical diagnosis, finance, etc. They are particularly important for their ability to reproduce and model nonlinear processes. Despite the computational-intensive training phase, which is ideally performed just once, ANN can infer results without the need of an iterative algorithm. Such ANN abilities make this kind of approach much faster than the iterative methods and does not require an initial guess for the solution. Furthermore, for real-time applications, it is possible to implement the network in special purpose hardware, thereby exploiting the intrinsically parallel nature of neural networks and thus, achieving very high processing speeds. To account for the coupled nature of the current and voltage signals in memristive devices, we consider an AlexNET [16] architecture and train it to estimate the DMM parameters from a graphical representation of the I-V loops. From the results reported, we show that following this approach, the DMM model accurately fits the experimental data measured under different conditions at a reduced computational cost compared with the standard error minimization method. 2. Materials and Methods 2.1. Device Fabrication and Electrical Characterization The RRAM devices being tested consist of Si/Al/TiN/Y2 O3 -x/Pt stacks as schemati- cally illustrated in Figure 1a. Oxygen-engineering [17] of the functional layer was carried out using a reactive molecular beam epitaxy setup that controls the oxygen stoichiometry of the yttria film by varying the flow of oxygen radicals and the film growth rate. A total of 2. Materials and Methods 2.1. Device Fabrication and Electrical Characterization The RRAM devices being tested consist of Si/Al/TiN/Y2O3‐x/Pt stacks as schemati‐ Micromachines 2022, 13, 2002 cally illustrated in Figure 1a. Oxygen‐engineering [17] of the functional layer was carried 3 of 14 out using a reactive molecular beam epitaxy setup that controls the oxygen stoichiometry of the yttria film by varying the flow of oxygen radicals and the film growth rate. A total of 5 different combinations of oxygen flow and growth rate were considered, including (i) 5 different combinations of oxygen flow and growth rate were considered, including (i) an an oxygen oxygen flowflow of 0.1 of 0.1 square square cubic cubic centimeter centimeter (sccm)with (sccm) witha agrowth growthraterateof of11angstrom angstrom perper second (Å s−−1), (ii) 0.2 sccm and 1 Å s−1, 0.3 sccm and 1 Å s−1, 0.5 sccm and 1 Å s−1 and 1 1 − 1 − 1 second (Å s ), (ii) 0.2 sccm and 1 Å s , 0.3 sccm and 1 Å s , 0.5 sccm and 1 Å s and − 1 1sccm sccmand and0.25 0.25ÅÅs s.−For a detailed description of the fabrication process andand the structural −1 1 . For a detailed description of the fabrication process the struc- characterization of the samples, please refer to [18]. The electrical characterization tural characterization of the samples, please refer to [18]. The electrical characterization was carried out with a Keithley 4200 semiconductor characterization system was carried out with a Keithley 4200 semiconductor characterization system (SCS) biasing (SCS) biasing the Pt top the electrode Pt top electrodeandand grounding groundingthe the TiNTiN bottom electrode. bottom TheThe electrode. internal current internal compliance current compli- (CC) of the SCS was used to prevent the hard breakdown of the ance (CC) of the SCS was used to prevent the hard breakdown of the oxide layer duringoxide layer during elec‐ troforming and SET. In addition, four different compliance levels were electroforming and SET. In addition, four different compliance levels were considered, and considered, and the I‐V characteristics were measured for over 100 consecutive cycles the I-V characteristics were measured for over 100 consecutive cycles for each fabricationfor each fabrication condition and condition and current current compliance compliance level level (Figure (Figure1b–f). 1b–f). Figure 1. Figure 1. (a) (a) Structure Structure of of the the RRAM RRAM devices devices under under test. test. The The insulating insulating layer layer varies varies from from 13 13 to to 18 18 nm nm depending on the oxidation condition. Measured I‐V loops from current compliances 50 depending on the oxidation condition. Measured I-V loops from current compliances 50 µA, 100 µA, μA, 100 μA, 1 mA, and 10 mA: (b) oxygen flow of 0.1 sccm, (c) oxygen flow of 0.2 sccm, (d) oxygen flow of 1 mA, and 10 mA: (b) oxygen flow of 0.1 sccm, (c) oxygen flow of 0.2 sccm, (d) oxygen flow of 0.3 sccm, (e) oxygen flow of 0.5 sccm, (f) oxygen flow of 1 sccm. 0.3 sccm, (e) oxygen flow of 0.5 sccm, (f) oxygen flow of 1 sccm. 2.2. Dynamic Memdiode 2.2. Memdiode Model Model (DMM) (DMM) The resistive The resistive switching switching (RS) (RS) mechanism mechanism is is the the fundamental fundamental physical physical phenomenon phenomenon behind ReRAM behind ReRAM devices. devices. In the the particular particular cases cases ofof CBRAMs CBRAMs (conducting (conducting bridge bridge random random access memories) access memories) andand OxRAMs OxRAMs (oxide-based (oxide‐based random-access random‐access memories), RS relies on on the the displacementofofmetal displacement metal ions/oxygen ions/oxygen vacancies vacancies withinwithin the dielectric the dielectric film in afilm in a Metal‐ Metal-Insulator- Metal (MIM) structure. The displacement is caused by the application of an external electrical stimulus, current or voltage [19–22]. The migration of ions originates from the alternate completion and destruction of a conductive filament (CF) spanning across the insulating film. For a ruptured CF, the device is in the high resistance state (HRS), often described by an exponential I-V relationship, while the completion of the CF leads to the low resistance state (LRS), which often exhibits a linear I-V curve [23,24]. In between these two extreme situations, the modulation of the CF transport properties renders intermediate ternal electrical stimulus, current or voltage [19–22]. The migration of ions originates from the alternate completion and destruction of a conductive filament (CF) spanning across the insulating film. For a ruptured CF, the device is in the high resistance state (HRS), often described by an exponential I‐V relationship, while the completion of the CF leads Micromachines 2022, 13, 2002 to the low resistance state (LRS), which often exhibits a linear I‐V curve [23,24]. In between 4 of 14 these two extreme situations, the modulation of the CF transport properties renders inter‐ mediate states by voltage‐controlled redox reactions. From the modelling viewpoint, the DMM is able to describe the major and minor I‐V loops and the gradual transitions in states by voltage-controlled redox reactions. From the modelling viewpoint, the DMM bipolar resistive switches for a wide variety of memristive systems. This is accomplished, is able to describe the major and minor I-V loops and the gradual transitions in bipolar as shownswitches resistive in the inset for of Figure a wide 2a, byofconsidering variety memristive asystems. nonlinear transport This equation based is accomplished, on as shown two identical opposite‐biased diodes in series with a resistor. The I‐V relationship in the inset of Figure 2a, by considering a nonlinear transport equation based on two resem‐ bles a diode identical with memorydiodes opposite-biased and that in is whywith series this device wasThe a resistor. termed memdiode. resembles I-V relationship Notice thata the anti‐parallel connected diodes allow the bi‐directional current diode with memory and that is why this device was termed memdiode. Notice thatflow through the de‐ the vice; as for both positive and negative polarities, there will always be a forward anti-parallel connected diodes allow the bi-directional current flow through the device; as biased diode. for both For the sake positive andofnegative completeness, the there polarities, DMMwillis succinctly always be areviewed forward in the next biased para‐ diode. For graphs. the sake of completeness, the DMM is succinctly reviewed in the next paragraphs. Figure Figure 2. 2. (a) (a) Memory Memorystatestate(λ, (λ, in in blue) blue) as as aa function functionof ofthe theapplied appliedvoltage, voltage,defined definedbybyEquation Equation(2). (2). The inset on the left shows the equivalent circuit model for the current equation (Equation The inset on the left shows the equivalent circuit model for the current equation (Equation (1) includ- (1) in‐ cluding the series ing the series resistance. resistance. TheThe diodes diodes are are drivendriven by bythethe memory memory state state thethe of of device device and and one one diode diode is is activated at a time. Typical I‐V characteristic for a memdiode obtained via simulation of the pro‐ activated at a time. Typical I-V characteristic for a memdiode obtained via simulation of the proposed posed model are superimposed. Current evolution is indicated by the blue arrows. (b) Experimental model are superimposed. Current evolution is indicated by the blue arrows. (b) Experimental I‐V loops of Y2O3 devices fabricated with an oxygen flow of 0.1 sccm and measured with current I-V loops of of compliances Y2 O devices 103 mA (bluefabricated squares) and with100 anμAoxygen flow of 0.1 (red circles), sccm both and fitted measured with the DMM with current model. compliances of 10 mA (blue squares) and 100 µA (red circles), both fitted with the DMM model. Physically, the memdiode is associated with a potential barrier that controls the elec‐ Physically, the memdiode is associated with a potential barrier that controls the elec- tron flow in the CF. The conduction properties of this non‐linear device change according tron flow in the CF. The conduction properties of this non-linear device change according to the variation of this barrier. Since the area of the CF is uncertain, instead of the potential to the variation of this barrier. Since the area of the CF is uncertain, instead of the potential barrier height, the diode current amplitude is used as the reference variable. Following barrier height, the diode current amplitude is used as the reference variable. Following Chua’s memristive device theory, the memdiode model comprises two equations: one for Chua’s memristive device theory, the memdiode model comprises two equations: one for the electron transport (transport equation, TE) and a second equation for the memory state the electron transport (transport equation, TE) and a second equation for the memory state of of the the device device (memory (memory equation, equation, ME). ME). The The transport transport equation equation can can be be derived derived from from the the quantum point‐contact (QPC) model [25–28], which uses the finite‐bias quantum point-contact (QPC) model [25–28], which uses the finite-bias Landauer approach Landauer ap‐ proach for the calculation for the calculation of theflow of the current current in a flow in a nanosized nanosized filamentary filamentary structure structure [29]. [29]. For a wide For a wide constriction formed by many conducting channels, neither constriction formed by many conducting channels, neither the number N of elemental the number N of elemental filamentary filamentary structures structures involved is involved known, nor is can known, theirnor can their specific specific potential potential barrier bar‐ parameters rier parameters (width 𝜃 and height 𝜑 ) be accessed individually [30]. Therefore, (width θi and height ϕi ) be accessed individually [30]. Therefore, we consider the following we consider the following heuristic approximation for TE which heuristic approximation for TE which can be derived from the QPC model: can be derived from the QPC model: I (VC ) = I0 (λ)sinh[α(λ)(VC − RS (λ) I )] (1) where I0 (λ) = Imin (1 − λ) + Imax λ is the diode current amplitude, α a fitting constant, and RS a series resistance. Imin and Imax are the minimum and maximum values of the current amplitude, respectively. While RS in (1) accounts for the contact resistance, the hyperbolic sine function expresses the barrier resistance. Making such an approximation to avoid the use of the complete Landauer formula is reasonable as Equation (1) has the same functional asymptotes for large I0 values ( RS I ≈ VC ) and low applied voltages ( RS I  VC ). As I0 Micromachines 2022, 13, 2002 5 of 14 increases in Equation (1), the I-V curve changes its shape from exponential to linear through a continuum of states as experimentally observed for these kinds of devices. For the sake of completeness, α and RS in (6) receive a similar treatment in the LTSpice script as that given to I0 (λ). Both parameters can be swept from a minimum (OFF) to a maximum (ON) if required. If not required, α and RS remain fixed. λ is a control parameter that represents the memory state of the device and runs from 0 (HRS) to 1 (LRS) and vice versa. The ME is conveniently described by the following differential equation [31]: dλ 1−λ λ = − (2) dt τS (λ, V ) τR (λ, V ) where τS,R are characteristic times associated with the SET (V > 0) and RESET (V < 0) transitions, i.e., with the ionic/defect movement within the dielectric film in one or the opposite direction; in the framework of this model, they are represented as: τS (λ, V ) = e−ηS (V −VS (λ)) (3) γ (V −V ) τR (λ, V ) = eηR λ S (4) where ηS,R and VS,R are the transition rates (ηS , ηR > 0) and the reference switching voltages (VS > 0, VR < 0), respectively. The exponential dependences of (3) and (4) on V are a consequence of the ionic/vacancy dynamics associated with the hopping mechanism [32]. The modelling of the snap back (Vs (λ) = Vs ) and snap forward (γ = 0) features of the I-V loops are de-activated for the sake of simplicity [33], mainly due to the fact that the current compliance in the SET partly masks the snap back effect. The combination of Equations (1) and (2) results in an I-V loop such as the one that is superimposed to the λ-V characteristic illustrated in Figure 2a, which starts in HRS (λ = 0) and evolves as indicated by the blue arrows. Another relevant feature of the proposed model is that it can be described by a simple SPICE script as shown in Ref. [3]. Finally, the accuracy of the model is reported in Figure 2b by fitting some of the loops reported in Figure 1. In summary, the proposed DMM not only provides a simple SPICE-compatible implementation for the resistive memory devices, but also a versatile one, as it can accurately fit the major and minor I-V loops measured in a wide variety of RRAM devices [33]. 2.3. Convolutional Neural Networks (CNNs) Besides their use in classification and prediction tasks, neural networks were also considered for parameter estimation in mathematical models. This is the case of the works by Dua [34], Morshed et al., Parikh et al. [35], Gonçalves et al. [36] and Rudi et al. [37], among others, where different kinds of neural networks were explored (Fully Connected Neural Networks—FCNN, Convolutional Neural Networks—CNN, Generative Adversar- ial Networks—GAN, etc.) for processing time series. Notably, in most of these scenarios, if not all, authors considered user-generated data with a constant sampling rate (invariant number of data points). However, experimental data from I-V measurements considering different maximal/minimal amplitudes inevitably implies waveforms with a variable num- ber of datapoints. This is also dictated by the internal electronics of the SCS. Moreover, I-V loops contain information from two different time series, that is, current and voltage as a function of time. While it would be possible to consider a CNN for multivariate time series analysis, in our case, we opted for a CNN that estimates the DMM parameters directly from a graphical representation of the I-V loop, similar to the classification of audio tracks in terms of their respective spectrogram representation [38]. This has two advantages over the use of the raw measured current and voltage signals: first, we avoid interpolating the data to fit a fixed sized array length, and second, the necessity of the normalization step required for the input signals is eliminated, which, given the exponential nature of the current signal, can be cumbersome from the practical viewpoint (logarithm of negative numbers for raw data). The input to the neural network used in this work consists of images in a constant resolution (227 × 227 pixels) containing the I-V curves generated by parameters directly from a graphical representation of the I‐V loop, similar to the classifi‐ cation of audio tracks in terms of their respective spectrogram representation [38]. This has two advantages over the use of the raw measured current and voltage signals: first, we avoid interpolating the data to fit a fixed sized array length, and second, the necessity of the normalization step required for the input signals is eliminated, which, given the Micromachines 2022, 13, 2002 exponential nature of the current signal, can be cumbersome from the practical viewpoint 6 of 14 (logarithm of negative numbers for raw data). The input to the neural network used in this work consists of images in a constant resolution (227 × 227 pixels) containing the I‐V curves generated considering by considering fixed axes for all loops.fixed Basedaxes onfor all images, these loops. Based the CNNon these images, responds withthe CNN a vector responds with a vector of floats normalized in the range (0,1) representing of floats normalized in the range (0,1) representing the 10 fitting parameters of the DMM. the 10 fitting parameters Because of theof the bigDMM. Because difference betweenof the bigand Imax difference between Ithe Imin parameters, max and Imin parameters, approach described the aboveapproach described above is recommended. is recommended. The CNN considered here The CNN considered is a variation of here is a variation the AlexNET [16] of the AlexNET network and whose[16] structure network is and whose in depicted structure Figure 3. is The depicted networkin Figure 3. Theusing was trained networkthe was trained Adam usingAlgorithm, Optimizer the Adam while Optimizer Algorithm, the also considering while also considering minimization of the the meanminimiza‐ squared tion errorofasthe themean squared for loss function error as 50 over theepochs. loss function Dropout forlayers over 50of epochs. 50% were Dropout included layers of before 50% wereconnected the fully included before thereduce layers to fully connected overfitting.layers to reduce overfitting. Figure 3. Structure of the convolutional neural network (Alex-NET) used to estimate the DMM Figure 3. Structure of the convolutional neural network (Alex‐NET) used to estimate the DMM fit‐ fitting parameters. ting parameters. It is worth mentioning that the choice of the AlexNET over other CNNs with a It is worth mentioning that the choice of the AlexNET over other CNNs with a rec‐ recognized higher accuracy, (such as the more recent ResNET) responds to the fact that ognized higher accuracy, (such as the more recent ResNET)responds to the fact that our our main target was not to find the best network in terms of accuracy but to provide a fast main target was and energy not to efficient find way tothe bestthe adjust network in terms parameters of of accuracy but a memristor to provide SPICE a fast compact and model. energy efficient way to adjust the parameters of a memristor SPICE compact For this reason, we searched for a trade-off between the accuracy and power consumption model. For this reason, we searched for a trade‐off between the accuracy and power of the neural network training/inference processes. In this connection, it has been shownconsumption of the neural network training/inference processes. In this connection, it by Yang et al. [39] that considering energy-aware pruning techniques, ResNET networkshas been shown by Yang incur et al. [39] higher that (up to considering one order ofenergy‐aware magnitude) pruning techniques, ResNET energy consumption networks than their AlexNET in‐ cur higher (up to one order of magnitude) energy consumption than their counterparts. This is due to the intrinsic, greater depth of the ResNET networks, which AlexNET coun‐ terparts. results inThis more is feature due to the mapsintrinsic, than ingreater depth of an AlexNET the that CNN ResNET needsnetworks, which to be moved results from the in more feature maps than in an AlexNET CNN that needs to be moved storage to the arithmetic unit(s) used for data computation and back to the storage. A from the storage breakdown of the energy consumption [39] reveals that most of the energy is used for such data movement and that the amount of data transferred increases with the number of feature maps. This also makes deep ResNETs much slower than the AlexNETs [40]. As such, ResNETs and other advanced CNNs, although capable of achieving higher accuracies, have proven to be not as convenient in terms of energy efficiency. Micromachines 2022, 13, 2002 7 of 14 2.4. Database Generation The database used in this work consists of a collection of 1200 I-V loop images gener- ated from pairs of experimentally measured I(t) and V(t) signals (with a variable number of points for each loop) plotted as I(t) vs. V(t) on a fixed x-axis running from −3 V to 3 V, and a fixed y-axis running from 1 nA up to 100 mA. The I-V images were stored with a fixed resolution of 227 × 227 px. in grayscale. Each image was then accompanied by a 10-element vector containing the normalized DMM parameters (indicated in Section 2.2) that fits the measured data, which acted as the image label in a supervised training procedure. The entire database is partitioned in a 5 to 1 ratio, thereby producing 1000 image-parameter pairs for the training phase and 200 image-parameter pairs for the test phase. The most challenging part of constructing such a database is obtaining the fitting parameter for each of the 1200 I-V loops. Manually fitting this number of I-V curves, with a model comprising up to 10 parameters, is impractical and prone to errors. Therefore, we developed an open-source, simulator-in-the-loop (NGSPICE) Python-approach to sys- tematize the data fitting procedure. Files containing the measured data were first parsed and stored in an SQLite database (DB). Then, each I-V loop was individually analyzed to extract the so-called “observable parameters”, such as the HRS, LRS, and HRS/LRS ratios, and fitted. The fitting phase can be split into 3 parts: (i) The I-V loop data is divided into 4 segments: HRS region (from the maximal voltage during the RESET Vmax→0, 1 in Figure 1c), LRS region (from the minimal volage applied during SET, Vmin→0, 2 ), SET region (from 0→Vmax, 3 ), and RESET region (from 0→Vmin, 4 ). (ii) Each previously mentioned segment is fitted using an approximation of the DMM mathematically derived for the region of interest (and thus neglecting or keeping constant the parameters associated with out-of-scope regions). (iii) Considering the previous fitting as an initial guess, the optimum parameter values are found by numerical optimization, consisting in an iterative simulation with the SPICE version of the DMM model. (iv) For performance comparisons, which are described later in the paper, the time re- quired to fit each loop was recorded. 2.5. ANN-Based DMM Fitting Procedure The procedure to obtain the best fitting parameters for a given experimental I-V loop is described in Figure 4. The starting point is the SQL database of measured V(t) and I(t) signals, which are first filtered based on the Mahalanobis distance [41] (with p = 0.01). This is carried out to discard anomalous curves (see Figure 5 for the case of samples fabricated with an oxygen flow of 0.1 sccm and measured with a 10 mA current compliance). V(t) and I(t) signals are then plotted together as indicated in Section 2.4. At this point, two different scenarios should be distinguished. On one hand, and with the purpose of developing/testing this methodology, the dataset described in Section 2.4 is considered. In this context, the CNN is trained following a supervised approach, and tested based on the previously fitted 1000 and 200 I-V loops, respectively. On the contrary, in a real case application scenario, only a fraction of the total I-V loops is fitted using the simulator-in-the-loop fitting procedure with the purpose of generating the dataset needed to train the CNN. The DMM parameters for the remaining measurements are then estimated using the trained CNN. Micromachines 2022, Micromachines2022, 13, 2022,13, xx 13,2002 88 of of 15 Micromachines 8 of1514 Figure Figure 4.4. Flowchart Flowchart of of the the proposed proposed fitting fitting procedure, procedure, including including the the training training of of the the fitting fitting network. network. Figure The 4. Flowchartpredicted of the proposed fitting procedure, including the training of the fittingfitting network. The The distribution distribution of of predicted values values is is analyzed analyzed in in the the following following sections sections to to discuss discuss the the fitting accu‐ accu‐ distribution racy. of predicted values is analyzed in the following sections to discuss the fitting accuracy. racy. Figure Figure Example of the 5. Example the outlierdetection detection usingthe the Mahalanobisdistance. distance. Notethat that inthis this example Figure 5. 5. Example ofof the outlier outlier detection using using the Mahalanobis Mahalanobis distance. Note Note that in in this example example we we have considered only 2 variables for the sake of a graphical representation and 4 anomalous we have considered only 2 variables for the sake of a graphical representation and 4 anomalous have considered only 2 variables for the sake of a graphical representation and 4 cases anomalous cases were were casesfound found (indicated (indicated were found in red) in red)infor (indicated for which red) which the resulting the resulting for which the resulting p‐value p-value was was above p‐value above 0.1. 0.1. 0.1. was above For For the the first first scenario scenario (the scenario (the one covered (the one one covered by covered by this by this work), this work),the work), thenetwork the networkperformance network performancefor performance for for fitting fitting the DMM DMM parameters parameters is is evaluated evaluated based based onon two two metrics: metrics: First, First, the fitting the DMM parameters is evaluated based on two metrics: First, the distribution ofthe distribution distribution of of thethe the fitted fitted fitted DMM DMMDMM parameters parameters parameters using using using the the the CNN CNN CNN is is compared is compared compared withwith with the the distribution the distribution distribution of theof of the the cor‐ cor‐ corresponding responding parameters responding parameters parameters obtained obtained obtained with with with thethe the simulator-in-the-loop simulator‐in‐the‐loop simulator‐in‐the‐loop fittingprocedure fitting fitting procedureusing procedure using using the Wasserstein distance [41]. Second, the mean squared error (MSE) between the I-V loops Micromachines 2022, 13, x 9 of 15 Micromachines 2022, 13, 2002 9 of 14 the Wasserstein distance [41]. Second, the mean squared error (MSE) between the I‐V loops generated by the predicted DMM parameters and the target I‐V loops are evaluated. generated by the predicted DMM parameters and the target I-V loops are evaluated. The The results obtained following this procedure are reported in the next section. results obtained following this procedure are reported in the next section. 3.3. Results Results and and Discussion Discussion InIn this this section, section, the the accuracy accuracy of the predicted of the predicted DMM DMM parameter parametervaluesvalueswith withthe thetwo two different metrics described in the previous section is assessed. Let us first considerthe different metrics described in the previous section is assessed. Let us first consider the superposition of the histograms corresponding to the DMM superposition of the histograms corresponding to the DMM parameters obtained from parameters obtained from the testtest the database database (used (usedas as labels or or labels target target values values ininthe thesupervised supervisedlearning learningapproach approachherehere considered) with considered) with the predicted DMM parameter values obtained by the ANNunder the predicted DMM parameter values obtained by the ANN understudy study after50 after 50training trainingepochs epochs(see(seeFigure Figure6a–j). 6a–j).Note Notethatthatthe thepredicted predictedvalues valuesfollow followveryverysimilar sim‐ ilar distributions when compared to the experimental ones; despite distributions when compared to the experimental ones; despite an occasional minor shift an occasional minor shift in theinmean, the mean, the standard the standard deviation deviation is fully is fully consistent consistent as isasthe is the shape shape of the of the distribu‐ distribution. tion. For instance, the 3 peaks observed in the histogram of the target For instance, the 3 peaks observed in the histogram of the target Imax values are nicely I max values are nicely reproducedby reproduced by the the IImax values estimated by the CNN. Similarly, the positive and negative max values estimated by the CNN. Similarly, the positive and negative skew observed in the histogram of skew observed in the histogram ofthe thetarget targetααmin and VRESET values are well captured by min and VRESET values are well captured by their CNN fitted counterpart. However, such a nicematching their CNN fitted counterpart. However, such a nice matchingisishighly highlysensitive sensitivetotothe the number of training epochs, as it can be seen by the Wasserstein distance vs. number number of training epochs, as it can be seen by the Wasserstein distance vs. numberofof epochspresented epochs presented in in the the inset inset of of each each panel panel in in Figure Figure6.6.Interestingly, Interestingly,not notevery everyparameter parameter has the same sensitivity. In opposition to the gradual reduction observed inthe has the same sensitivity. In opposition to the gradual reduction observed in therest restofofthe the parameters, the Wasserstein distance for I max drops abruptly at the very first training epoch parameters, the Wasserstein distance for Imax drops abruptly at the very first training epoch andthen and thenremains remains relatively relatively constant. constant. Figure6. Figure 6. Comparison Comparison ofof the the histograms histograms of of the the estimated estimatedDMMDMMparameters parametersusing usingthe theANN ANN(labelled (labelled as CNN fit) against the DMM parameters extracted by the simulator‐in‐the‐loop approach (labelled as CNN fit) against the DMM parameters extracted by the simulator-in-the-loop approach (labelled as LSM fit). The ANN was trained over 50 epochs. Note the good match between both histograms. as LSM fit). The ANN was trained over 50 epochs. Note the good match between both histograms. The Wasserstein distance in the inset quantifies the fitting improvement along the training. The Wasserstein distance in the inset quantifies the fitting improvement along the training. Nonetheless, it is worth pointing out that the previous metric does not necessarily demonstrate that the ANN fitted parameters properly capture the details of the experi- Micromachines 2022, 13, x 10 of 15 Micromachines 2022, 13, 2002 10 of 14 Nonetheless, it is worth pointing out that the previous metric does not necessarily demonstrate that the ANN fitted parameters properly capture the details of the experi‐ mental data, mental data, and and even even less less quantifies quantifies it. it. For For this this reason, reason, wewe evaluated evaluated thethe MSE MSE between between the experimental I‐V data and a simulated I‐V loop recreated from the the experimental I-V data and a simulated I-V loop recreated from the calculated fitting calculated fitting parameters. Figure 7 presents these results as a function of the training epoch, showing aa parameters. Figure 7 presents these results as a function of the training epoch, showing steep decrease steep decrease in the MSE in the MSE forfor the the first first 10 10 epochs. epochs. From From the 12th epoch the 12th epoch on, the MSE on, the MSE error error remains almost constant. In addition, not only does the mean MSE decrease, remains almost constant. In addition, not only does the mean MSE decrease, but so does but so does its dispersion. its dispersion. It It is is worth worth pointing pointing out out that that despite despite this this behavior, behavior, the the Wasserstein Wasserstein distance distance shows a rather constant decrease for the entire range of training epochs. This suggests that shows a rather constant decrease for the entire range of training epochs. This suggests that the MSE the MSE might might bebe less less sensitive sensitive to to certain certain parameters parameters (for (for instance instance ηηSET SET and andηηRESET). RESET ). Figure 7. Evaluation metrics metrics of of the theproposed proposedparameter parameterextraction extractiontechnique techniqueusing usingANNs. ANNs.MSE MSEasas a a functionofofthe function thetraining trainingepoch. epoch.TheTheMSE MSEisiscomputed computedbetween betweenthe themeasured measuredI-V I‐V loop loop and the one obtained by simulation using the ANN fitted fitted parameters. parameters. A final verification verificationofofthe thesuitability suitabilityofofthe the proposed proposed neural neural network network approach approach to to es‐ estimate timate thetheDMMDMM parameters parameters is illustrated is illustrated in Figure in Figure 8. Figure, 8. In this In this the Figure, the I-V I‐V loops loops resulting resulting from the from running running DMMthe DMMwith model modelthewith fittedthe fitted parameters parameters are superimposed are superimposed to the to the exper‐ experimental input data. Note that the Figure covers the 5 different oxidation imental input data. Note that the Figure covers the 5 different oxidation cases (from left cases (from left to right) to right) and and the 4the 4 different different current current compliances compliances (from(from top totop to bottom) bottom) as well. as well. In thisInway, this way, both the abrupt transition (digital) observed in the most oxidized samples both the abrupt transition (digital) observed in the most oxidized samples measured with measured with a high a high current current compliance, compliance, as well as well assoft as the the soft transition transition (quasi-analogue) (quasi‐analogue) observed observed in in the the least oxidized samples, are captured. least oxidized samples, are captured. Having verified that the DMM fitting parameters are correctly predicted by this ANN-based procedure, we proceeded to quantify the achieved time saving. To do so, we estimated the required time to fit a total of 100,000 I-V loops. For the case of the fitting time required for the simulator-in-the-loop approach, we estimated it by considering a linear regression on the fitting-time vs. fitted cycles plot (squared blue markers in Figure 9). We do so because fitting such an amount of data, considering the required time to fit each loop with this method, would take more than 106 s (or up to approximately 11 days). On the contrary, for the case of the fitting time required by the ANN-based approach, we extended our dataset by replicating the dataset roughly 100 times. This is similar to the conventional data augmentation techniques used to reduce the overfitting problems in neural networks, with the difference that the added data is not new. This is not a problem in our case since, at this point, we are not using such dataset for training or testing accuracy, but simply to quantify the time required to forward propagate the information through the network. As shown in Figure 9, the proposed method could potentially reduce the fitting time by more than 1 order of magnitude, even if we consider the required time to generate the training data and the training of the ANN. Micromachines Micromachines2022, 13,13, 2022, 2002 x 1111ofof15 14 Figure RepresentativeI-V Figure8.8.Representative I‐V loops loops for for each each oxidation-current oxidation‐current compliance compliance pair pair (totaling (totaling20 20different different combinations).Each combinations). Eachcolumn columnfrom fromleft lefttotoright rightgroups groups the the samples samples fabricated fabricated with with thethe same same oxida‐ oxidation tion condition (0.1, 0.2, 0.3, 0.5, 1 sccm) and each row summarizes the measurements condition (0.1, 0.2, 0.3, 0.5, 1 sccm) and each row summarizes the measurements carried out with carried out with the same current compliance (50 μA. 100 μA, 1 mA, and 10 mA). In each panel, one experimen‐ the same current compliance (50 µA. 100 µA, 1 mA, and 10 mA). In each panel, one experimentally tally measured I‐V loop is plotted together with a simulation of the DMM whose parameters were measured I-V loop is plotted together with a simulation of the DMM whose parameters were obtained obtained with the trained ANN. with the trained ANN. Having verified that the DMM fitting parameters are correctly predicted by this It is also worth emphasizing that further time-saving could be achieved by speeding ANN‐based procedure, we proceeded to quantify the achieved time saving. To do so, we up the training step. Training speed is determined by the learning method and the de- estimated the required time to fit a total of 100,000 I‐V loops. For the case of the fitting gree of parallelization (i.e., the hardware used to run the training phase). The learning time required for the simulator‐in‐the‐loop approach, we estimated it by considering a method employed in this work is the well-known Adam Optimizer, which has become linear regression on the fitting‐time vs. fitted cycles plot (squared blue markers in Figure the default learning method regardless of the field and has set the reference in terms of 9). We do so because fitting such an amount of data, considering the required time to fit convergence speed (although the Stochastic Gradient6 Descend method might provide a each loop with this method, would take more than 10 s (or up to approximately 11 days). better generalizing capability). As such, the remaining alternative to speed-up the training On the contrary, for the case of the fitting time required by the ANN‐based approach, we is increasing the degree of the parallelization of the hardware used for the training, i.e., extended our dataset by replicating the dataset roughly 100 times. This is similar to the using a cluster with a higher number of GPUs [42]. Another option (also exploiting the conventional data augmentation techniques used to reduce the overfitting problems in benefits of parallelization) is the use of FPGA accelerators [43–46], which can deliver a neural networks, with the difference that the added data is not new. This is not a problem reduction in the processing time per image of up to 2 orders of magnitude, as reported in in our case since, at this point, we are not using such dataset for training or testing accu‐ Ref. [43]. In that work, an end-to-end FPGA-based CNN accelerator, with all the layers racy, but simply to quantify the time required to forward propagate the information working concurrently in a pipelined structure, remarkably improves the performance. through the network. As shown in Figure 9, the proposed method could potentially Micromachines 2022, 13, x 12 of 15 Micromachines 2022, 13, 2002 12 of 14 reduce the fitting time by more than 1 order of magnitude, even if we consider the re‐ quired time to generate the training data and the training of the ANN. Figure9.9. Evaluation Figure Evaluation metrics metrics of ofthe theproposed proposedparameter parameterextraction extractiontechnique techniqueusing usingANNs. ANNs. Time Time required to fit an increasing number of I‐V loops using the simulator‐in‐the‐loop, iterative approach, required to fit an increasing number of I-V loops using the simulator-in-the-loop, iterative approach, and the ANN‐based counterpart. For a large number of loops (e.g., 100,000), the ANN approach and the ANN-based counterpart. For a large number of loops (e.g., 100,000), the ANN approach allows a significant time saving. allows a significant time saving. It is also worth emphasizing that further time‐saving could be achieved by speeding 4. Conclusions up the training step. Training speed is determined by the learning method and the degree In this work, we investigated the use of the AlexNET architecture, a particular case ofaparallelization of convolutional (i.e., neuralthenetwork, hardware forused to run the extracting the training phase). The fitting parameters learning used to model method the employed in this work is the well‐known Adam Optimizer, which I-V loops from RRAM devices with the Dynamic Memdiode Model (DMM). Experimental has become the default learning data method regardless corresponding to devicesofwith the field and has different set the oxygen reference content in terms ofwith and measured convergence different speed (although the Stochastic Gradient Descend method might provide current compliance levels were considered. In view of this scenario as a case study, a better general‐ we izing capability). As such, the remaining alternative to speed‐up the training evaluated the impact of the training epochs on the fitting accuracy and quantified it in terms is increasing the of thedegree relativeof the mean parallelization squared error. of We the proposed hardware aused for the to procedure training, train thei.e.,neural usingnetwork a cluster with labeled data—this being I-V loops with their associated DMM fitting parameters.paral‐ with a higher number of GPUs [42]. Another option (also exploiting the benefits of This lelization)was approach is the use to proven of achieve FPGA accelerators [43–46],time a reduced running which and can deliver aenergy an increased reduction in the efficiency processing when fittingtime many per image I-V loops,ofwhich up to was 2 orders moreofthan magnitude, what was astypically reportedobtained in Ref. [43]. In that using the work, an end‐to‐end FPGA‐based CNN accelerator, with all the Least-Square-Error minimization method, which needs to be performed individually on layers working concur‐ rently each in a pipelined individual structure, remarkably improves the performance. measurement. 4. Conclusions Author Contributions: Conceptualization, F.L.A. and E.M.; methodology, F.L.A. and E.M.; software, F.L.A.In and E.M.; this validation, work, F.L.A. and E.M.; we investigated formal the use analysis, of the AlexNETF.L.A.architecture, and E.M.; investigation, F.L.A., a particular case E.P., L.A. and E.M.; resources, C.H., J.S., L.A. and E.M.; data curation, E.P.; writing—original of a convolutional neural network, for extracting the fitting parameters used to model the draft preparation, F.L.A. and E.M.; writing—review and editing, F.L.A., E.P., N.K., T.V., S.P., J.G., T.O., C.H., I‐V loops from RRAM devices with the Dynamic Memdiode Model (DMM). Experimental J.S., L.A. and E.M.; visualization, F.L.A. and E.M.; supervision, L.A. and E.M.; project administration, data corresponding to devices with different oxygen content and measured with different C.H., J.S., L.A. and E.M.; funding acquisition, C.H., J.S., L.A. and E.M. All authors have read and agreed to compliance current the publishedlevels versionwere of theconsidered. manuscript. In view of this scenario as a case study, we evaluated the impact of the training epochs on the fitting accuracy and quantified it in terms ofThis Funding: the work wasmean relative conducted in the squared frameWe error. of the project StorAIge. proposed The project a procedure to trainhas thereceived neural funding from the ECSEL Joint Undertaking (JU) under grant agreement No 101007321. The JU network with labeled data—this being I‐V loops with their associated DMM fitting pa‐ receives support from the European Union’s Horizon 2020 research and innovation programme rameters. This approach was proven to achieve a reduced running time and an increased and France, Belgium, Czech Republic, Germany, Italy, Sweden, Switzerland, Turkey. Funding was energy efficiency when fitting many I‐V loops, which was more than what was typically received within the ECSEL Joint Undertaking project WAKeMeUP in collaboration with the European obtained Union’s H2020using the Least‐Square‐Error Framework minimization Program (H2020/2014-2020) andmethod, Nationalwhich needsunder Authorities, to begrant per‐ formed individually agreement on eachby No 783176. Funding individual the Federalmeasurement. Ministry of Education and Research (BMBF) under contracts 16MEE0154 and 16ESE0298 is gratefully acknowledged Also, funding from DFG under project AL 560/13-2, 384682067, and AL 560/21-1 are gratefully acknowledged. The authors also acknowledge the financial support provided by Merck KGaA, Darmstadt, Germany, in the framework of the joint research platform by the company and TU Darmstadt called “Sustainability Hub”. E.M. Micromachines 2022, 13, 2002 13 of 14 and J.S. acknowledge the support provided by the European project MEMQuD, code 20FUN06, which has received funding from the EMPIR programme co-financed by the Participating States and from the European Union’s Horizon 2020 research and innovation programme. Conflicts of Interest: The authors declare no conflict of interest. References 1. Strukov, D.B.; Snider, G.S.; Stewart, D.R.; Williams, R.S. The missing memristor found. Nature 2008, 453, 80–83. [CrossRef] [PubMed] 2. Miranda, E. Compact Model for the Major and Minor Hysteretic I-V Loops in Nonlinear Memristive Devices. IEEE Trans. Nanotechnol. 2015, 14, 787–789. [CrossRef] 3. Aguirre, F.L.; Pazos, S.M.; Palumbo, F.; Suñé, J.; Miranda, E. Application of the Quasi-Static Memdiode Model in Cross-Point Arrays for Large Dataset Pattern Recognition. IEEE Access 2020, 8, 202174–202193. [CrossRef] 4. Aguirre, F.L.; Gomez, N.M.; Pazos, S.M.; Palumbo, F.; Suñé, J.; Miranda, E. Minimization of the Line Resistance Impact on Memdiode-Based Simulations of Multilayer Perceptron Arrays Applied to Pattern Recognition. J. Low Power Electron. Appl. 2021, 11, 9. [CrossRef] 5. Levenberg, K. A Method for the Solution of Certain Non-Linear Problems in Least Squares. Q. Appl. Math. 1944, 2, 164–168. [CrossRef] 6. Marquardt, D.W. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. J. Soc. Ind. Appl. Math. 1963, 11, 431–441. [CrossRef] 7. Singh, K.; Jain, P. BSIM3v3 to EKV2.6 Model Parameter Extraction and Optimisation using LM Algorithm on 0.18 µ Technology node. Int. J. Electron. Telecommun. 2018, 64, 5–11. [CrossRef] 8. Karr, C.L.; Weck, B.; Massart, D.L.; Vankeerberghen, P. Least median squares curve fitting using a genetic algorithm. Eng. Appl. Artif. Intell. 1995, 8, 177–189. [CrossRef] 9. Fogel, D.B. Evolutionary Computation: The Fossil Record; IEEE Press: New York, NY, USA, 1998; ISBN 9780780334816. 10. Jervase, J.A.; Bourdoucen, H.; Al-Lawati, A. Solar cell parameter extraction using genetic algorithms. Meas. Sci. Technol. 2001, 12, 1922–1925. [CrossRef] 11. Spałek, T.; Pietrzyk, P.; Sojka, Z. Application of the genetic algorithm joint with the Powell method to nonlinear least-squares fitting of powder EPR spectra. J. Chem. Inf. Model. 2005, 45, 18–29. [CrossRef] 12. Hooke, R.; Jeeves, T.A. “Direct Search” Solution of Numerical and Statistical Problems. J. ACM 1961, 8, 212–229. [CrossRef] 13. Torczon, V. On the Convergence of Pattern Search Algorithms. SIAM J. Optim. 1997, 7, 1–25. Available online: http://www.siam. org/journals/ojsa.php (accessed on 25 September 2022). [CrossRef] 14. Ayala-Mató, F.; Seuret-Jiménez, D.; Escobedo-Alatorre, J.J.; Vigil-Galán, O.; Courel, M. A hybrid method for solar cell parameter estimation. J. Renew. Sustain. Energy 2017, 9, 063504. [CrossRef] 15. Bishop, C.M.; Roach, C.M. Fast curve fitting using neural networks. Rev. Sci. Instrum. 1992, 63, 4450–4456. [CrossRef] 16. Krizhevsky, A.; Sutskever, I.; Hinton, G.E. ImageNet classification with deep convolutional neural networks. Commun. ACM 2017, 60, 84–90. [CrossRef] 17. Sharath, S.U.; Vogel, S.; Molina-Luna, L.; Hildebrandt, E.; Wenger, C.; Kurian, J.; Duerrschnabel, M.; Niermann, T.; Niu, G.; Calka, P.; et al. Control of Switching Modes and Conductance Quantization in Oxygen Engineered HfOx based Memristive Devices. Adv. Funct. Mater. 2017, 27, 1700432. [CrossRef] 18. Petzold, S.; Piros, E.; Eilhardt, R.; Zintler, A.; Vogel, T.; Kaiser, N.; Radetinac, A.; Komissinskiy, P.; Jalaguier, E.; Nolot, E.; et al. Tailoring the Switching Dynamics in Yttrium Oxide-Based RRAM Devices by Oxygen Engineering: From Digital to Multi-Level Quantization toward Analog Switching. Adv. Electron. Mater. 2020, 6, 2000439. [CrossRef] 19. Lee, A.R.; Bae, Y.C.; Im, H.S.; Hong, J.P. Complementary resistive switching mechanism in Ti-based triple TiO x /TiN/TiO x and TiOx/TiOxNy/TiOx matrix. Appl. Surf. Sci. 2013, 274, 85–88. [CrossRef] 20. Duan, W.J.; Wang, J.B.; Zhong, X.L.; Song, H.J.; Li, B. Complementary resistive switching in single sandwich structure for crossbar memory arrays. J. Appl. Phys. 2016, 120, 084502. [CrossRef] 21. Yang, M.; Wang, H.; Ma, X.; Gao, H.; Hao, Y. Voltage-amplitude-controlled complementary and self-compliance bipolar resistive switching of slender filaments in Pt/HfO2/HfOx/Pt memory devices. J. Vac. Sci. Technol. B Nanotechnol. Microelectron. Mater. Process. Meas. Phenom. 2017, 35, 032203. [CrossRef] 22. Chen, C.; Gao, S.; Tang, G.; Fu, H.; Wang, G.; Song, C.; Zeng, F.; Pan, F. Effect of electrode materials on AlN-based bipolar and complementary resistive switching. ACS Appl. Mater. Interfaces 2013, 5, 1793–1799. [CrossRef] [PubMed] 23. Aguirre, F.L.; Rodriguez-Fernandez, A.; Pazos, S.M.; Sune, J.; Miranda, E.; Palumbo, F. Study on the Connection Between the Set Transient in RRAMs and the Progressive Breakdown of Thin Oxides. IEEE Trans. Electron Devices 2019, 66, 3349–3355. [CrossRef] 24. Fröhlich, K.; Kundrata, I.; Blaho, M.; Precner, M.; Tapajna, M.; Klimo, M.; Šuch, O.; Škvarek, O. Hafnium oxide and tantalum oxide based resistive switching structures for realization of minimum and maximum functions. J. Appl. Phys. 2018, 124, 152109. [CrossRef] 25. Sune, J.; Miranda, E.; Nafria, M.; Aymerich, X. Point contact conduction at the oxide breakdown of MOS devices. In Proceedings of the IEEE Technical Digest—IEDM, San Francisco, CA, USA, 6–9 December 1998; pp. 191–194. Micromachines 2022, 13, 2002 14 of 14 26. Miranda, E.; Suñé, J. Analytic modeling of leakage current through multiple breakdown paths in SiO2 films. In Proceedings of the IEEE International Reliability Physics Symposium Proceedings, Orlando, FL, USA, 30 April–3 May 2001; Institute of Electrical and Electronics Engineers Inc.: New York, NY, USA, 2001; Volume 2001, pp. 367–379. 27. Miranda, E.; Suñé, J. Electron transport through broken down ultra-thin SiO2 layers in MOS devices. Microelectron. Reliab. 2004, 44, 1–23. [CrossRef] 28. Miranda, E.; Walczyk, C.; Wenger, C.; Schroeder, T. Model for the Resistive Switching Effect in HfO2 MIM Structures Based on the Transmission Properties of Narrow Constrictions. IEEE Electron Device Lett. 2010, 31, 609–611. [CrossRef] 29. Datta, S. Electronic Transport in Mesoscopic Systems, 1st ed.; Cambridge Studies in Semiconductor Physics and Microelectronic Engineering; Cambridge University Press: Cambridge, UK, 1997; ISBN 978-0521599436. 30. Miranda, E.; Mehonic, A.; Suñé, J.; Kenyon, A.J. Multi-channel conduction in redox-based resistive switch modelled using quantum point contact theory. Appl. Phys. Lett. 2013, 103, 222904. [CrossRef] 31. Miranda, E.; Sune, J. Memristive State Equation for Bipolar Resistive Switching Devices Based on a Dynamic Balance Model and Its Equivalent Circuit Representation. IEEE Trans. Nanotechnol. 2020, 19, 837–840. [CrossRef] 32. Rodriguez-Fernandez, A.; Cagli, C.; Sune, J.; Miranda, E. Switching Voltage and Time Statistics of Filamentary Conductive Paths in HfO2-based ReRAM Devices. IEEE Electron Device Lett. 2018, 39, 656–659. [CrossRef] 33. Aguirre, F.L.; Suñé, J.; Miranda, E. SPICE Implementation of the Dynamic Memdiode Model for Bipolar Resistive Switching Devices. Micromachines 2022, 13, 330. [CrossRef] 34. Dua, V. An Artificial Neural Network approximation based decomposition approach for parameter estimation of system of ordinary differential equations. Comput. Chem. Eng. 2011, 35, 545–553. [CrossRef] 35. Parikh, J.; Kozloski, J.; Gurev, V. Integration of AI and mechanistic modeling in generative adversarial networks for stochastic inverse problems. arXiv 2020, arXiv:2009.08267. 36. Gonçalves, P.J.; Lueckmann, J.M.; Deistler, M.; Nonnenmacher, M.; Öcal, K.; Bassetto, G.; Chintaluri, C.; Podlaski, W.F.; Haddad, S.A.; Vogels, T.P.; et al. Training deep neural density estimators to identify mechanistic models of neural dynamics. Elife 2020, 9, 1–46. [CrossRef] [PubMed] 37. Rudi, J.; Bessac, J.; Lenzi, A. Parameter Estimation with Dense and Convolutional Neural Networks Applied to the FitzHugh- Nagumo ODE. arXiv 2021, arXiv:2012.06691. 38. Abdel-Hamid, O.; Mohamed, A.R.; Jiang, H.; Deng, L.; Penn, G.; Yu, D. Convolutional neural networks for speech recognition. IEEE Trans. Audio Speech Lang. Process. 2014, 22, 1533–1545. [CrossRef] 39. Yang, T.-J.; Chen, Y.-H.; Emer, J.; Sze, V. A Method to Estimate the Energy Consumption of Deep Neural Networks. In Proceedings of the 2017 51st Asilomar Conference on Signals, Systems, and Computers, Pacific Grove, CA, USA, 29 October–1 November 2017. 40. Bianco, S.; Cadene, R.; Celona, L.; Napoletano, P. Benchmark analysis of representative deep neural network architectures. IEEE Access 2018, 6, 64270–64277. [CrossRef] 41. Montgomery, D.C.; Runger, G.C. Applied Statistics and Probability for Engineers; John Wiley & Sons: Hoboken, NJ, USA, 2010; ISBN 0470053046. 42. Sun, P.; Feng, W.; Han, R.; Yan, S.; Wen, Y. Optimizing Network Performance for Distributed DNN Training on GPU Clusters: ImageNet/AlexNet Training in 1.5 min. arXiv 2019, arXiv:1902.06855. [CrossRef] 43. Li, H.; Fan, X.; Jiao, L.; Cao, W.; Zhou, X.; Wang, L. A high performance FPGA-based accelerator for large-scale convolutional neural networks. In Proceedings of the 2016 26th International Conference on Field Programmable Logic and Applications (FPL), Lausanne, Switzerland, 29 August–2 September 2016. [CrossRef] 44. Qiu, J.; Wang, J.; Yao, S.; Guo, K.; Li, B.; Zhou, E.; Yu, J.; Tang, T.; Xu, N.; Song, S.; et al. Going deeper with embedded FPGA platform for convolutional neural network. In Proceedings of the 2016 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, Monterey, CA, USA, 21–23 February 2016; pp. 26–35. [CrossRef] 45. Suda, N.; Chandra, V.; Dasika, G.; Mohanty, A.; Ma, Y.; Vrudhula, S.; Seo, J.S.; Cao, Y. Throughput-optimized openCL-based FPGA accelerator for large-scale convolutional neural networks. In Proceedings of the 2016 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, Monterey, CA, USA, 21–23 February 2016; pp. 16–25. [CrossRef] 46. Zhang, C.; Li, P.; Sun, G.; Guan, Y.; Xiao, B.; Cong, J. Optimizing FPGA-based accelerator design for deep convolutional neural networks. In Proceedings of the 2015 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays, Monterey, CA, USA, 22–24 February 2015; pp. 161–170. [CrossRef]