Simulating MOSFET equations with ModelingToolkit

New to Julia
1

Hello, I’m trying to model a MOSFET with ModelingToolkit using the most basic quadratic model, but having a lot of issues. The equations I came up with so far are as follows:

using ModelingToolkit, OrdinaryDiffEq


function NMOS(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    return ((Vgs > Vth) * (Vds < Vgs-Vth))*uCox*W/L*((Vgs-Vth)*Vds-Vds^2/2) +
           ((Vgs > Vth) * (Vds > Vgs-Vth))*uCox/2*W/L*(Vgs-Vth)^2
end

@parameters t Vth W L KP C R Vcc f
@variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
@derivatives D'~t

eqs = [
    Vg ~ 1+sin(2*pi*f*t),
    Id ~ NMOS(Vg, Vc, 0, Vth, W, L, KP),
    Ic ~ C*D(Vc),
    Ir ~ (Vcc-Vc)/R,
    0 ~ Ir - Id - Ic,
]

Basically I tried to generate a sine input to the mosfet gate, which is connected to a resistor to vcc, and a capacitor at the output. The MOSFET model is the most basic quadratic model, where I implemented the different operating regions by multiplying with booleans.

So far, baring any silly mistakes, no problems. But then I try to make a system of equations from this and am not having much luck.

Despite knowing that my MOSFET is nonlinear, I started by just putting it in ODESystem to see what happens. Apparently ERROR: LoadError: type Sym has no field op. Okay, time to move on to Modeling Nonlinear Systems · ModelingToolkit.jl

ns = NonlinearSystem(eqs, [Vg, Vc, Id, Ic , Ir], [t, Vth, W, L, KP, C, R, Vcc, f])
nlsys_func = generate_function(ns)[2] # second is the inplace version
fn = eval(nlsys_func)
du = zeros(5); u = ones(5)
params = (0, 0.7, 10e-6, 1e-6, 2.0718e-5, 1e-6, 10e3, 5, 1e3)
fn(du,u,params)
du

Everything seems fine until I call the eval’d function. At that point I get all sorts of type errors that I don’t understand:

ERROR: LoadError: MethodError: objects of type Term{Real} are not callable
ERROR: LoadError: MethodError: Cannot `convert` an object of type Term{Float64} to an object of type Float64

They seem to come from the MOSFET and capacitor line.

2 Likes
2

Hi and welcome!

I’ve been wanting to learn ModelingToolkit for some time but haven’t even kicked the tires yet.

I suspect your issue stems from the fact that your circuit really has only one differential equation due to the single reactive element, C.

I think I could solve your problem using Differential equations.jl. It would require manually writing an equation for dVC/dt.

3

I rearranged your code into the style I’ve been using for ODEs lately. I use Parameters to avoid retyping parameter names any more than necessary and StaticArrays for performance. It seems to be running okay. I made a guess about tspan.

using DifferentialEquations
using Plots
using Parameters
using StaticArrays

@with_kw struct Params
	Vth:: Float64 = 0.7
	W:: Float64   = 10e-6
	L:: Float64   = 1e-6
	KP:: Float64  = 2.0718e-5
	C:: Float64   = 1e-6
	R:: Float64   = 10e3
	Vcc:: Float64 = 5
	f:: Float64   = 1e3
end

p = Params()

function NMOS(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    return ((Vgs > Vth) * (Vds < Vgs-Vth))*uCox*W/L*((Vgs-Vth)*Vds-Vds^2/2) +
           ((Vgs > Vth) * (Vds > Vgs-Vth))*uCox/2*W/L*(Vgs-Vth)^2
end

function ode(u,p,t)
	Vc, = u
	@unpack_Params p # unpacks Params fields by name into this function's scope
	Vg = 1+sin(2*pi*f*t)
	Id = NMOS(Vg, Vc, 0, Vth, W, L, KP)
	Ir = (Vcc-Vc)/R
	Ic = Ir - Id
	dVc = 1/C*Ic
	@SVector[dVc]
end

function sim()
	u0 = @SVector[0.0]
	tspan = (0.0, 100e-3)
	prob = ODEProblem(ode, u0, tspan, p)
	solve(prob,abstol=1e-8, reltol=1e-8))
end

sol=sim()

plot(sol)

EDIT: t should not have been in Params, set tolerances to better values.

4

Here are some tricks to get ModelingToolkit to almost work. I was able to get the NonlinearSystem to work, but it seems this example doesn’t work with the ODE solvers yet. I have mostly stopped using ModelingToolkit in favor of the style @klaff used here.

using ModelingToolkit
using OrdinaryDiffEq
# Core.ifelse cannot have new methods added, so you need this package
using IfElse: ifelse
using Plots


function NMOS(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Voth = Vgs-Vth # over threshold
    # ModelingToolkit doesn't handle Boolean expressions, this sometimes works
    ifelse(signbit(Voth),
           zero(uCox/2*W/L*Voth^2), # below threshold is 0
           ifelse(signbit(Vds-Voth),
                  uCox*W/L*(Voth*Vds-Vds^2/2),
                  uCox/2*W/L*Voth^2))
    # return ((Vgs > Vth) * (Vds < Voth))*uCox*W/L*(Voth*Vds-Vds^2/2) +
    #        ((Vgs > Vth) * (Vds > Voth))*uCox/2*W/L*Voth^2
end

# This works and prints a value for t==0
function solvepoint()
    @parameters t Vth W L KP C R Vcc f
    @variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
    @derivatives D'~t

    eqs = [
        Vg ~ 1+sin(2π*f*t),
        Id ~ NMOS(Vg, Vc, 0, Vth, W, L, KP),
        # Ic ~ C*D(Vc), # I'm not sure why, but differentials must appear at
        # least once on the LHS.
        D(Vc) ~ Ic/C,
        Ir ~ (Vcc-Vc)/R,
        0 ~ Ir - Id - Ic,
    ]

    ns = NonlinearSystem(eqs, [Vg, Vc, Id, Ic, Ir], [t, Vth, W, L, KP, C, R, Vcc, f])
    # expression=Val{false} to return function instead of expresion
    fn = generate_function(ns; expression=Val{false})[2] # second is the inplace version
    du = zeros(5); u = ones(5)
    params = (0, 0.7, 10e-6, 1e-6, 2.0718e-5, 1e-6, 10e3, 5, 1e3)
    fn(du,u,params)
    du
end

function plotsolution()
    @parameters t Vth W L KP C R Vcc f
    @variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
    @derivatives D'~t

    eqs = [
        Vg ~ 1+sin(2*pi*f*t),
        Id ~ NMOS(Vg, Vc, 0, Vth, W, L, KP),
        D(Vc) ~ Ic/C,
        Ir ~ (Vcc-Vc)/R,
        0 ~ Ir - Id - Ic,
    ]

    sys = ODESystem(eqs, t, [Vg, Vc, Id, Ic , Ir], [Vth, W, L, KP, C, R, Vcc, f])
    u0 = ones(5)
    tspan = (0, 0.010)
    params = [
              Vth => 0.7,
              W => 10e-6,
              L => 1e-6,
              KP => 2.0718e-5,
              C => 1e-6,
              R => 10e3,
              Vcc => 5,
              f => 1e3
             ]
    # fails here with
    # ERROR: TypeError: non-boolean (Term{Bool}) used in boolean context
    # Looks like a bug. Maybe worth a report?
    prob = ODEProblem(sys, u0, tspan, params, jac=true)
    soln = solve(prob, Tsit5())
    plot(soln)
end
5

Very interesting. I looked a bit at ModelingToolkit.jl but I did not think it was useful to perform circuit simulations.

I thought everything under the SciML umbrella was for machine learning. I thought it was mainly for creating models from empirical data. I must admit I have recently heard discussions on using these tools with “physics-based modelling” (Karen Willcox) to get improvements - but I don’t really have a grasp on how things fit together:

→ JuliaCon 2020 Live Session Wednesday - YouTube

I personally was leaning towards Modia.jl as a solution to this kind of thing:

https://github.com/ModiaSim/Modia.jl

Modia appears to follow in the footsteps of the Modelica toolset. I find the solutions look alot more like what you write in VerilogA/VerilogAMS - so I had a better idea how to approach the problem. I suppose I should look more into this ModelingToolkit package, then.

RFI

But if someone can point me to some useful information comparing these tools (Modia, ModelingToolkit, SciML) in the context of creating circuit simulations/simulators, I would very much appreciate it :slight_smile: .

1 Like
6

DifferentialEquations.jl can definitely be used to do circuit simulations but you can’t (today) import a netlist and instantiate models of ICs by part number and be off and running like you can with a SPICE program. Instead you need to be able to write your equations and do a bit of algebra up front and in back to get all the answers.

On the other hand, you can describe a modest amount of circuitry that way, along with equations for some physical things, and simulate control algorithms using callbacks (which can even call C code if you wish), letting you do things that are difficult or impossible with SPICE.

7

That’s fair. I don’t expect to be able to read in a netlist in a package that is not designed specifically for circuit simulations. I’ll look into that.

ModelingToolkit

But what about ModelingToolkit? It sounds like it is not ready for this kind of problem. Is there a reason why this is? Why would I want to use ModelingToolkit instead of just DifferentialEquations+solve()?

8

I don’t know enough to answer that.

9

I’ll be putting a lot more out on this in the near future, but the last JuliaCon video starts to describe one aspect of the story:

The problem is that the best way to express a model is not necessarily the most efficient way to solve a model. For example, a lot of circuit models are index-2 DAEs (due to Kirchhoff’s law), and so while you could in theory solve it sometimes, in many cases there are numerical issues if you solve the model as the user wishes to write it. So there are tricks that a domain specific language will use in order to make it solvable. And there are the tricks that you know: write the model in log-space to keep variables positive for example.

All of these tricks form a set of transformations on a model which allow for solving the model in a numerically better form. Automating these transformations would essentially be a compiler on models for performing symbolic changes to numerical codes. That is ModelingToolkit.

There are many flavors of this. Automatically transforming code into equivalent faster code (the last JuliaCon video). Automatically parallelizing code. Automatically performing index reduction, removing variables which have explicit algebraic relationships, and solving the smallest version of the model while deducing the equations for how to transform back into the original variables. This is the symbolic-numeric system of ModelingToolkit as it applies to differential equations, nonlinear solving, optimization, and more.

The full ecosystem is going to be unveiled over the next few months, and I hope to detail it at the next JuliaCon. But essentially ModelingToolkit is a composable symbolic-numeric modeling language. The extra tools in:

https://github.com/JuliaComputing/StructuralTransformations.jl

make it a superset of the Modelica language (and thus Modia), where with a specific choice of transformations you get the standard Modelica behavior. Similarly, standard high performance SPICE implementations perform subsets of these tricks to reach their speed, so a DSL over ModelingToolkit would give a high performance SPICE (note: one is in development). And for this reason, there are a ton of both open and proprietary modeling projects built on ModelingToolkit, from https://pumas.ai/ and a new thermodynamics library to GitHub - SciML/Catalyst.jl: Chemical reaction network and systems biology interface for scientific machine learning (SciML). High performance, GPU-parallelized, and O(1) solvers in open source software.. Launching a high performance DSL on ModelingToolkit takes away the mathematical burden and makes it into a representation burden: how do I represent the model symbolically, since the symbolic system will generate the most efficient numerical code with the hooks into DifferentialEquations.jl etc. to then simulate fast.

But there’s so much more coming, with automated generation of surrogates of models:

based on a new SciML method forming the backbone of approximate transformations.

So it’s a huge project, a lot is documented but what’s seen is only a small part of its scope. Early 2021 will be all about sharing the demonstrations that are now beginning to be completed.

9 Likes
10

? One of the core tools for incorporating models with empirical data is the ability to simulate said model. Everything is about scientific simulators, like n-body simulators, PDE solvers, quadrature methods, and more which all allow for incorporating differentiable programming and machine learning into and around the core methods. So the starting point is always methods for simulating models. How this all comes together is documented at https://sciml.ai/.

2 Likes
11

I appreciate your explanation and, more importantly, your great work on all of this. But it seems to me that SciML might have a marketing/perception problem. After all it is named “Scientific Machine Learning” so it shouldn’t come as a surprise that someone not familiar with the story you try to tell will assume it is all about ML.

As a longish-time Julia user I loosely followed the emergence of the SciML org in Julia as a super cool project. However, I got really confused when I started to read “scientific machine learning” everywhere (in the @ChrisRackauckas universe :smiley:). Above you link DifferentialEquations.jl as an example of SciML. It hasn’t been part of the org a few years back. I seem to have missed the point where the good old diff eq solver package became “DifferentialEquations.jl: Scientific Machine Learning (SciML) Enabled Simulation and Estimation” (title taken from the docs). The very first indication of what DiffEq is about mentions scientific machine learning. With a bit of exaggeration, a new user might think: Huh, machine learning? I thought I was checking out a package for solving ODEs? (Note that on GitHub the subtitle is missing and the first mention of SciML comes under “Additionally, DifferentialEquations.jl comes with built-in analysis features, including:” which, IMHO, seems much more appropriate.

Other examples:
RecursiveArrayTools.jl: “Tools for easily handling objects like arrays of arrays and deeper nestings in scientific machine learning (SciML) and other applications”. It is really necessary to mention, and potentially confuse, a reader with a SciML reference here?
ExponentialUtilities.jl: “Utility functions for exponential integrators for the SciML scientific machine learning ecosystem”.

Personally, I basically take SciML as a branding effort. And of course, this is both understandable and perhaps necessary. But, personally, I hope that the grouping of efforts under the SciML org doesn’t come at the expense of making packages with a rather well-defined purpose less well-defined (or at least make it seem this way).

Anyways, these are just my two cents and maybe I’m exaggerating a bit. Hope that it is still a valuable perspective from a non-SciML scientist.

All the best.

5 Likes
12

Thanks for the improvements and explanations. I’m very new to all the Julia stuff so understanding how the libraries work together is very useful. I’m very excited about the SPICE on top of Julia idea, which is in fact how I got here. (It was mentioned on the Sky130 Slack)

Based on @contradict’s improvements I was able to get a little bit further.
Not sure if it’s a bug or a feature, but ModelingToolkit kinda seems sensitive to the form of my equations.
It was already mentioned that it wants the D(Vc) on the LHS, but looking into the error I was getting a bit more, it seems to also expect the LHS of an algebraic equation to be 0. Is it just expecting you to follow certain conventions?

After that it complained that it could not make a mass matrix out of my system. After a bit of trial and error with different solvers, I got it to kinda run: Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.

From the partial plot it seemed to go a bit crazy below the threshold voltage. Increasing Vgs so that it stays above Vth made the simulation run! Maybe it ends up trying to relate Id to Vg, which has no relation at all below Vth. Maybe the best workaround is actually to use a better transistor model that models weak inversion?

model

My code so far:

using ModelingToolkit
using StructuralTransformations
using OrdinaryDiffEq
using IfElse: ifelse
using Plots


function NMOS(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Vo = Vgs-Vth # overdrive
    # ModelingToolkit doesn't handle Boolean expressions, this sometimes works
    ifelse(signbit(Vo),
           zero(Vo), # below threshold is 0
           ifelse(signbit(Vds-Vo),
                  uCox*W/L*(Vo*Vds-Vds^2/2),
                  uCox/2*W/L*Vo^2))
end

@parameters t Vth W L KP C R Vcc f
@variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
@derivatives D'~t

# It does not like it AT ALL if you put functions on LHS or derivatives on RHS
eqs = [
    0 ~ -Vg+2+sin(2*pi*f*t),
    0 ~ -Id+NMOS(Vg, Vc, 0, Vth, W, L, KP),
    D(Vc) ~ Ic/C,
    0 ~ -Ir+(Vcc-Vc)/R,
    0 ~ Ir - Id - Ic,
]

sys = ODESystem(eqs, t, [Vg, Vc, Id, Ic , Ir], [Vth, W, L, KP, C, R, Vcc, f])
u0 = zeros(5)
tspan = (0.0, 10e-3)
params = [
    Vth => 0.7,
    W => 100e-6,
    L => 1e-6,
    KP => 2.0718e-5,
    C => 1e-6,
    R => 10e3,
    Vcc => 5,
    f => 1e3
]
# Turn into a first order differential equation system
# sys = ModelingToolkit.ode_order_lowering(sys)

# Perform index reduction to get an index 1 DAE
# sys = StructuralTransformations.dae_index_lowering(sys)
prob = ODEProblem(sys, u0, tspan, params, jac=true)
# soln = solve(prob,abstol=1e-8, reltol=1e-8)
soln = solve(prob, RadauIIA5(), abstol=1e-8, reltol=1e-8)
plot(soln)
png("model")

Any feedback or improvements appreciated.

I’m realizing I’m kind of rusty on all the algebra definitions. Like… with the nonlinear transistor model and the algebraic Kirchhoff laws, I was sorta expecting to have a nonlinear Differential Algebraic Equation, but in the end it’s handled just fine as an Ordinary Differential Equations. One solver even complained it could only be used on a DAE and not on an ODE.

So seems like I’m using “algebraic” and “nonlinear” in a general sense, while they have a very specific meaning in this context? Can imagine for example that the only non-linearity that matters is that of a differential, while my MOS is purely algebraic and the capacitor perfectly linear. Or that an DAE refers to something particular like algebraic loops, rather than just a system of equations that has purely algebraic parts. I’m sure this was in some lectures I attended, but seems I need to refresh those definitions a bit.

1 Like
13

This is off-topic, but I’m willing to discuss it. It should probably be kicked into a new thread though.

People looking for a differential equation solver don’t need to be told it’s a differential equation solver. That’s clear and obvious from the name and now thousands of materials out there. What we need to do early and often is emphasize “what else?”. Model discovery, Bayesian estimation, optimization with respect to parameter uncertainty, etc. Not only are these things on the outside, but there’s also neural networks “on the inside”, with NNODE solving ODEs via neural networks and the new Qu&Co methods using differentiable quantum circuits to generate ODE solvers for quantum computers. That it has ODE solvers is obvious and doesn’t need too much writing, but that these things exist needs a bit of note.

If you look at the uses, it’s a core library of the SciML ecosystem since almost all outputs are some form of AbstractVectorOfArray. For people looking for where these arrays of arrays and deeper nestings in the SciML ecosystem are defined, yes this is the place to look. So yes, referencing the organization that it was made it as a base tool and where more than half of the direct dependents are? Seems like it makes sense?

They are the utility functions for exponential integrators in SciML organization. There is an effort to make them all differentiable for these model discovery purposes as well. That is a pretty good description of what it is, is it not? Maybe we could reference more on expmv, but most people don’t know what those tools are and are looking for exponential integrators.

See the roadmap:

The packages aren’t “less well-defined”, the issue is the true scope of the work and the organization. Most of the repositories in the organization are not differential equation tools. Some of the most used packages in the organization are not written in Julia or for Julia users. So as a general principle, JuliaDiffEq was outgrown. What pulled together this whole organization was this thread of general differentiability throughout scientific computing tools and integration with machine learning stacks, both in and around the methods. If you look at the tools as a whole and the roadmap, this as the central theme of the organization is absolutely clear and well-defined. If you think that our work is just in differential equation solving, then you are exactly the person who needs to be shown that it’s not, because that means there’s so much more that you’re missing and that you need help finding!

14

Is there a reason why you commented out the index lowering? @shashi @YingboMa you might want to take a look at making this into the tutorial we were discussing. The key might be missing the simple alias elimination?

A differential algebraic equation (DAE) is a differential equation which mixes in algebraic equations. A purely algebraic system is a NonlinearSystem, and a purely differential equation is an index-0 DAE. So indeed you can think of a DAE as an ODE with algebraic loops.

Right now MTK always targets ODEProblem and thus needs to use the methods for mass matrices:

https://diffeq.sciml.ai/stable/solvers/dae_solve/#OrdinaryDiffEq.jl-(Mass-Matrix)

But that should change soon.

15

My understanding is that these lowering passes mainly server to make things numerically faster and more stable. For this small problem they did not seem to make any difference, so eliminated them as a source of bugs while testing.

Since I have only one first order differential equation, I don’t think the order_lowering actually does anything? And for the index_lowering I’m not completely sure I thoroughly understand what it does. One of them actually got stuck in an infinite recursion, so I just got rid of them for now.

Your description of DAE makes sense… So basically… I’m confused by what type of system I’m having vs what kinds of systems the various functions in ModelingToolkit.jl support. Not really sure where the root of my confusion lies. Maybe it’s just that my understanding of an ODE is more narrow than the actual definition, which seems to be just that it has one independent variable, and not necessarily of the basic y=ax+bx’+cx’'… form. Maybe I was just thrown off by the special DAE and NonlinearSystem functions, thinking that therefore ODE only handles purely linear differential equations.

So yea… I think it’s going in the right direction. Once I find some time I’ll try implementing the most basic actual spice mosfet model, and see if that solves the convergence issue with the subthreshold region.

1 Like
16

Off topic

Agreed. I find this discussion really useful (unless, of course me & @carstenbauer were the only two people who struggle in understanding the SciML universe).

Thanks for clarifying

Thank you very much @ChrisRackauckas. I cannot say I fully understand everything you said, but it significantly helped me to understand how SciML is structred and how fits in with everything else. Also: Thank you for linking in relevant videos.

Sadly, I have very limited understanding of this side of the modeling world. I mostly did circuit design, so I typically only use SPICE to run simulations. I have spent little-to-no time writing actual simulator code. At best, I have written VerilogA/AMS models for ideal circuit components that could dramatically speed up simulation times. In the end, I have not had the time to invest in core of simulator algorithms. My work was needed elsewhere.

Excellent. Thanks for that! If this was already explained in Modeling toolkit, I’m sorry I never reached that point. Again: once I noticed ModelingToolkit pas part of the SciML umbrella, I started to believe it was part of a machine learning layer that found model parameters that best fit the input data. …So I decided not to waste my time learning another new concept I didn’t need at that time.

Comparison to Modelica

That is an extremely useful description to me! I hope it will be easy to locate in the SciML documentation. I’m guessing other non-experts would be interested in knowing these correspondence points.

Future development

Very excited to hear this!

My pre-conceptions

Indeed this all makes sense. However, I assumed these activities would have been split into two separate layers:

  • A layer for scientific computation/simulation
  • A layer for machine learning/model fitting

…So when I saw that DifferentialEquations.jl was under the SciML organization, I (incorrectly assumed) this was the “machine-learning layer” that managed model fitting algorithms.

I should have known better: In Julia, package names are global (not relative to an organization/problem space). So, it makes sense that DifferentialEquations.jl was not limited to my imagined “machine-learning layer”. In that case, I guess it would have been called something like MLDifferentialEquations.jl :slight_smile: .

SciML origin question

So just out of curiosity, is that why the term “SciML” came to exist? People found out that the it was easier to implement these algorithms if the base libraries (like DifferentialEquations.jl) included hooks for machine learning algorithms? Was it that otherwise, a layered solution became hacky and unmanageable?

17

“MLDifferentialEquations.jl”, or helper functions for using DifferentialEquations.jl with ML libraries like Flux.jl is:

Scientific machine learning as a term comes from the DOE which put out a special report essentially defining the descipline:

At around the same time, we had been starting to put out blog posts and talks around the same topic:

And actually we had been doing some of this since the start of “JuliaDiffEq”:

That blog post on neural network differential equation solvers is from 2017, going to show how early this was embedded into the organization.

So in 2020 when we were going to apply for NumFOCUS support, we had to rethink the name a bit because it no longer worked. “JuliaDiffEq” had packages in R and in Python. It had optimization libraries, quadrature libraries, etc. mixed with libraries for reservoir computing, physics-informed neural networks, and more. But there was driving principle behind all of the work: building a high performance scientific computing stack that allows for this future where there is no difference between scientific computing and machine learning, where the two can be freely mixed, and where physical information overcomes the shortcomings of incomplete data. This of course directly aligned with what the DOE and others had just declared as the scientific goal for the next decade, so it made since to directly say, this is the SciML open source software organization. This is where you will see high performance differential equation solvers, and neural networks discovering differential equations from data, and symbolic-numeric modeling languages which augment a Modelica-like experience with neural surrogates. That is SciML.

And the final part.

What if with every new project you had to write the differential equation solvers from scratch? While this might be tenable for non-stiff equations (with some efficiency loss), when you start having to handle stiff equations there’s a reason why there are only a few production-quality general purpose stiff ODE solvers that have ever seen significant use! So when going down this route, we knew we had to create a foundation that would allow reuse of methods. On top of that, just doing the simplest connections does not work. This talk goes into detail as to how some methods that are used in PyTorch and Tensorflow libraries will straight up fail if used on stiff systems:

So this problem of mixing inference with the solvers is a very difficult problem that will take a lot of work to perfect. We started working on inverse problems in 2016, more than half of our development has always been towards handling of inverse problems, and here we are entering 2021 saying we’ve only started handling the hardest generalized forms of inference. From the SciMLBenchmarks you can see we’ve cracked the code with forward solvers, and while we will continue to improve things like stiff ODE solvers indefinitely, we can show on most equations that the most efficient method is something pure Julia at this point.

But even with all of that, there are still some inverse problems that are out of our reach. There are ways to get there, but a lot of my ideas for how to do it require deep integration between the ML/optimization areas and the differential equation areas. That’s a bit hard to describe right now because not all of the research is completed, so I can only point to a few scattered examples like the continuous-time echo state network examples, but I can assure you that throughout the 2020’s what you’ll see is a lot more here. Differential equation solvers in the optimization routines and the neural network layers, neural networks integrated into the stepping behaviors of some of the differential equation solvers, and symbolic techniques transforming models to and from mechanistic and data-driven forms. So basically, there’s a lot more where that came from and I think down the road about 50% of the organizations activities will continue to be on the basic traditional solver methods, but the other 50% will be pushing the front line of this interplay between the two fields.

7 Likes
18

Success??

model

I implemented the level 6 model, which is supposedly faster than level 3 but more accurate than level 1. It still does not have weak inversion, but with some hacks I got it to work. (level 3 seems a bit much for a toy)

I had to define my own max function to clamp sub-threshold voltages. And then I had to allow a bit of “convergence leakage current” to make the simulator not minstep and crash. And then I switched to a simple Trapezoidal solver, which did not freak out as much as the Runge-Kutta I had before. Probably less efficient, but fine for my toy.

It would be interesting to learn more about how to solve these kind of issues more properly, but that goes a bit deep for now maybe.

using ModelingToolkit
using StructuralTransformations
using OrdinaryDiffEq
using IfElse: ifelse
using Plots


function NMOS_L1(Vg, Vd, Vs, Vth, W, L, uCox)
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Vo = Vgs-Vth # overdrive
    # ModelingToolkit doesn't handle Boolean expressions, this sometimes works
    ifelse(signbit(Vo),
           zero(Vds), # below threshold is 0
           ifelse(signbit(Vds-Vo),
                  uCox*W/L*(Vo*Vds-Vds^2/2),
                  uCox/2*W/L*Vo^2))
end

mymax(x, y) = ifelse(signbit(x-y), y, x)

# https://www2.eecs.berkeley.edu/Pubs/TechRpts/1990/ERL-90-19.pdf
function NMOS_L6(Vg, Vd, Vs, Vb,
                 B, n, K, m, λ₀, λ₁, Vt0, γ, ϕf, W, Leff)
    Vbs = Vb-Vs
    Vgs = Vg-Vs
    Vds = Vd-Vs
    Vth=Vt0+γ*(√(2ϕf-Vbs)-√(2ϕf))
    Vdsat=K*mymax(1e-9, Vgs-Vth)^m
    Idsat=W/Leff*B*mymax(1e-9, Vgs-Vth)^n
    λ=λ₀-λ₁*Vbs
    Id5=Idsat*(1+λ*Vds) # Vds > Vdsat (saturated region)
    Id3=Id5*(2-(Vds/Vdsat))*Vds/Vdsat
    sat = Vds-Vdsat
    ifelse(signbit(sat), Id3, Id5)
end
# x = LinRange(0, 2, 100)
# res = NMOS_L6.(1, x, 0, 0,
#     4.9721e-5, 1.0484, 0.83496, 0.6193, 0.066265,
#     0.0038573, 0.85502, 0.29648, 0.20556/2, 100e-6, 1e-6)
# plot(x, hcat(res...)')

@parameters t W L C R Vcc f
@variables Vg(t) Vc(t) Id(t) Ic(t) Ir(t)
@derivatives D'~t

# It does not like it AT ALL if you put functions on LHS or derivatives on RHS
eqs = [
    0 ~ -Vg+1+sin(2*pi*f*t),
    0 ~ -Id+NMOS_L6(Vg, Vc, 0, 0,
                    4.9721e-5, 1.0484, 0.83496, 0.6193, 0.066265,
                    0.0038573, 0.85502, 0.29648, 0.20556/2, W, L),
    D(Vc) ~ Ic/C,
    0 ~ -Ir+(Vcc-Vc)/R,
    0 ~ Ir - Id - Ic,
]

sys = ODESystem(eqs, t, [Vg, Vc, Id, Ic , Ir], [W, L, C, R, Vcc, f])
u0 = zeros(5)
tspan = (0.0, 10e-3)
params = [
    W => 100e-6,
    L => 1e-6,
    C => 1e-6,
    R => 10e3,
    Vcc => 5,
    f => 1e3
]
# Turn into a first order differential equation system
# sys = ModelingToolkit.ode_order_lowering(sys)

# Perform index reduction to get an index 1 DAE
# sys = StructuralTransformations.dae_index_lowering(sys)
prob = ODEProblem(sys, u0, tspan, params, jac=true)
# soln = solve(prob,abstol=1e-8, reltol=1e-8)
soln = solve(prob, Trapezoid(), abstol=1e-8, reltol=1e-8)
plot(soln)
png("model")
3 Likes
19

Indeed Trapezoidal will be more conservative, so it’s a common method used in SPICE simulators. You might want to look into using a Gauss-Legendre method and see how those do in this context.