2
\$\begingroup\$

I have been attempting to port a modified version of the Yakopcic memristor from Python/MATLAB to LTspice but have been running into issues with the results not being the same, as shown in the following pictures.
In Python/MATLAB the memristor is simulated by using Euler's method to solve the IVP describing the device's internal state evolution.

MATLAB simulation result
[MATLAB]

LTspice simulation result
LTspice

I think that this answer may have put me on the right track (https://electronics.stackexchange.com/a/368910/294242) by pointing out that timing in the SPICE engine could lead to issues.
In fact, if I enlarge my timestep from dt=1/10000 to dt=1 in MATLAB/Python, I get the same exact result as in LTspice.

MATLAB simulation with dt=1 (the curve in black)
MATLAB

Is there any way to solve the issue I'm seeing? I have tried .tran 0 50s 0 {1/10000} to match dt=1/10000, but the results don't change.

enter image description here

* SPICE model for memristive devices * Created by Chris Yakopcic * Last Update: 8/9/2011 * * Connections: * TE - top electrode * BE - bottom electrode * XSV - External connection to plot state variable * that is not used otherwise .subckt MEM_YAKOPCIC TE BE XSV * Fitting parameters to model different devices * gmax_p, bmax_p, gmax_n, bmax_n: Parameters for OFF state IV relationship * gmin_p, bmin_p, gmin_n, bmin_n: Parameters for OFF state IV relationship * Vp, Vn: Pos. and neg. voltage thresholds * Ap, An: Multiplier for SV motion intensity * xp, xn: Points where SV motion is reduced * alphap, alphan: Rate at which SV motion decays * xo: Initial value of SV * eta: SV direction relative to voltage .param gmax_p=9e-5 bmax_p=4.96 gmax_n=1.7e-4 bmax_n=3.23 + gmin_p=1.5e-5 bmin_p=6.91 gmin_n=4.4e-7 bmin_n=2.6 + Ap=90 An=10 + Vp=0.5 Vn=0.5 + xp=0.1 xn=0.242 + alphap=1 alphan=1 + xo=0 eta=1 * Multiplicative functions to ensure zero state * variable motion at memristor boundaries .func wp(V) = xp/(1-xp) - V/(1-xp) + 1 .func wn(V) = V/xn * Function G(V(t)) - Describes the device threshold .func G(V) = + IF(V > Vp, + Ap*(exp(V)-exp(Vp)), + IF(V < -Vn, + -An*(exp(-V)-exp(Vn)), + 0 ) ) * Function F(V(t),x(t)) - Describes the SV motion .func F(V1,V2) = + IF(eta*V1 >= 0, + IF(V2 >= xp, + exp(-alphap*(V2-xp))*wp(V2), + 1 ), + IF(V2 <= xn, + exp(alphan*(V2-xn))*wn(V2), + 1 ) ) * IV Response - Hyperbolic sine due to MIM structure .func IVRel(V1,V2) = + IF(V1 >= 0, + gmax_p*sinh(bmax_p*V1)*V2 + gmin_p*sinh(bmin_p*V1)*(1-V2), + gmax_n*sinh(bmax_n*V1)*V2 + gmin_n*sinh(bmin_n*V1)*(1-V2) + ) * Circuit to determine state variable * dx/dt = F(V(t),x(t))*G(V(t)) Cx XSV 0 {1} .ic V(XSV) = xo Gx 0 XSV value={eta*F(V(TE,BE),V(XSV,0))*G(V(TE,BE))} * Current source for memristor IV response Gm TE BE value = {IVRel(V(TE,BE),V(XSV,0))} .ends MEM_YAKOPCIC

EDIT: It turns out that the issue was not in the MATLAB to SPICE porting but rather in the original MATALB script itself. The MATLAB/Python script did not use a the same timestep as the real simulation.
I measured the experimental timestep and found it to be dt~=0.08s, so by using Euler with dt=1/10000 the MATLAB script was effectively simulating for only 60 ms, instead of the 50 s in the real experiment. Scaling the Ap and An parameters to match the updated - much larger - timestep was enough to reproduce the experimental results.

\$\endgroup\$
6
  • \$\begingroup\$ Don't use Euler, backward or forward, if you have it, use trapezoidal. That's what most SPICE engines use (with various, proprietary tweaks). \$\endgroup\$ Commented Aug 20, 2021 at 14:32
  • \$\begingroup\$ I think that's what I'm doing: in LTspice I tried trapezoidal, modified trapezoidal, and Gear integration with no discernible outcomes. I also tried switching the solver from Normal to Alternate, again to no avail. I've also tried alternate solvers in Python (for ex., LSODA, RK4) and the results are the same as plain FE. (also, Yakopcic uses FE in his MATLAB simulations) \$\endgroup\$ Commented Aug 20, 2021 at 14:41
  • \$\begingroup\$ @ThomasTiotto Are you using simulink to solve the model or just matlab code \$\endgroup\$ Commented Aug 25, 2021 at 18:40
  • \$\begingroup\$ @ThomasTiotto , Check this doc., that explains every detail, including the Spice simulations. \$\endgroup\$ Commented Aug 26, 2021 at 15:12
  • \$\begingroup\$ show us the matlab code. What is your input frequency and waveform shape. \$\endgroup\$ Commented Aug 29, 2021 at 16:49

2 Answers 2

Reset to default
1
\$\begingroup\$

Memristors are slow but hysteresis increases with frequency.

try up to 50 Hz.

enter image description here

Compare with Falstad's results and reset or choose max scale after changing sliders, here showing the VI curve of Sig Gen.

\$\endgroup\$
1
  • \$\begingroup\$ the hysteresis (more precisely the area enclosed by the pinched loop) decreases with frequency not increase. \$\endgroup\$ Commented Aug 29, 2021 at 16:38
0
\$\begingroup\$

I don't know enough about the memristor application and matlab code to comment exactly why, but I can show that experimentally, if you allow the voltage input transient to pull down further, it will display the behavior you are looking for. Also, there's no step size dependency here, either (you can just set to default). Maybe there is some factor of 2 translation to LTspice that's not being accounted for?

enter image description here

Also, just inspecting the symmetry in the matlab result, there's a range of about -.5 +/-1.5 in the x range (with -.5 being the flat center of the transition). The modified Ltspice outcome, has a similar symmetry of -1.5 +/- 2.5v centered at flat region of -1.5. Stopping ltspice sim input at -2.0 is like stopping the matlab input stimulus short at -1 to -1.5, and thus it would have a similar outcome on matlab.

You can also see that in the original time domain, the voltage negative excursion is not enough to trip the current significantly until around the -3.5v to -4.0v threshold. What are the original input stimulus values to the matlab file?

enter image description here

I'll try to look a bit, but it would help to have the matlab source and hope this helps lead you to your answer if I don't finish by two days.

*update. The more I look at it, there are a lot of parameters that are tightly fit. Even the original thesis has drastically different values and even equations. I would verify if it's a parameter/equation translation somewhere. You really need the original python or matlab source code and stimulus you are comparing to. Also, I found your github. So it seems that you were aware of the -4 v already. Would have saved some time to link there (I can see you've spent a lot of time looking at this and I mean it in a good way).

I'd also add that I don't think the time step will do much. You can lower to 1uS min step and see plenty of points in the ltspice sim across all significant regions (even though they may not be uniformly spaced). And they look fairly smooth and continuous. The link you pointed to has to do with binning samples and windowing using the dft or fft. In that case, the number of samples in your window can dramatically effect the outcome resolution if not accounted for properly. In this case, there's no special operations like that I can see.

enter image description here

Here's a small snippet of python code that seems different from ltspice sub. This python mim hyperbolic sine function (from your git model) doesn't seem to account for the gmax/gmin terms. It might be somewhere else in your code, but if you remove the bottom right gmin term in ltspice, I vs v will be a flat line like your original comment. I can't figure if you or someone ported the python from ltspice or visa versa, but the original code from yakopcic thesis has mostly ltspice solutions, so I figure it started with ltspice. Or if it's from matlab, please post, where did the wrapper and translations originate?

\$\endgroup\$
1
  • \$\begingroup\$ Thanks for all the work you put into this. It turns out that the issue was not in the MATLAB to SPICE porting but rather in the original MATALB script itself. I have updated the original answer to reflect this. \$\endgroup\$ Commented Sep 20, 2021 at 8:32

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.