# Simulating MOSFET equations with ModelingToolkit
**URL:** https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408
**Category:** New to Julia
**Tags:** modelingtoolkit
**Created:** [December 26, 2020, 4:00pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408 "2020-12-26T16:00:24Z")
**Posts on this page:** 19
**Page:** 1
### Author:  [@pepijndevos](https://discourse.julialang.org/u/pepijndevos)
#### Post date: [December 26, 2020, 4:00pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/1 "2020-12-26T16:00:24Z")
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:
```julia
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](https://mtk.sciml.ai/dev/tutorials/nonlinear/)
```julia
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:
```julia
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.
---
### Author:  [@klaff](https://discourse.julialang.org/u/klaff)
#### Post date: [December 27, 2020, 5:58pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/2 "2020-12-27T17:58:22Z")
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.
---
### Author:  [@klaff](https://discourse.julialang.org/u/klaff)
#### Post date: [December 27, 2020, 7:25pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/3 "2020-12-27T19:25:44Z")
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`.
```julia
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.
---
### Author:  [@contradict](https://discourse.julialang.org/u/contradict)
#### Post date: [December 27, 2020, 8:46pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/4 "2020-12-27T20:46:31Z")
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.
```julia
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
```
---
### Author:  [@MA\_Laforge](https://discourse.julialang.org/u/MA_Laforge)
#### Post date: [December 27, 2020, 9:33pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/5 "2020-12-27T21:33:35Z")
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:
→ [https://youtu.be/gfHMPIxEyv4](https://youtu.be/gfHMPIxEyv4)
I personally was leaning towards Modia.jl as a solution to this kind of thing:
[https://github.com/ModiaSim/Modia.jl](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 🙂 .
---
### Author:  [@klaff](https://discourse.julialang.org/u/klaff)
#### Post date: [December 27, 2020, 10:08pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/6 "2020-12-27T22:08:02Z")
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.
---
### Author:  [@MA\_Laforge](https://discourse.julialang.org/u/MA_Laforge)
#### Post date: [December 27, 2020, 10:26pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/7 "2020-12-27T22:26:21Z")
> [@klaff](#):
>
> 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.
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()?
---
### Author:  [@klaff](https://discourse.julialang.org/u/klaff)
#### Post date: [December 27, 2020, 10:31pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/8 "2020-12-27T22:31:12Z")
I don’t know enough to answer that.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 28, 2020, 1:25am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/9 "2020-12-28T01:25:19Z")
> [@MA\_Laforge](#):
>
> ### 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()?
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:
[ ](https://www.youtube.com/watch?v=UNkXNZZ3hSw)
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](https://en.wikipedia.org/wiki/Pantelides_algorithm) 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](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/](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.](https://github.com/SciML/Catalyst.jl). 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:
[ ](https://www.youtube.com/watch?v=BGrJchagGLg)
based on [a new SciML method](https://arxiv.org/abs/2010.04004) 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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 28, 2020, 1:30am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/10 "2020-12-28T01:30:09Z")
> [@MA\_Laforge](#):
>
> I thought everything under the SciML umbrella was for machine learning. I thought it was mainly for creating models from empirical data.
> **[GitHub - SciML/DifferentialEquations.jl: Multi-language suite for...](https://github.com/SciML/DifferentialEquations.jl)**
>
> Multi-language suite for high-performance solvers of differential equations and scientific machine learning (SciML) components. Ordinary differential equations (ODEs), stochastic differential equat...
? 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](https://github.com/SciML/NBodySimulator.jl/pull/37), [PDE solvers](https://github.com/SciML/NeuralPDE.jl), [quadrature methods](https://github.com/SciML/Quadrature.jl), 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/](https://sciml.ai/).
---
### Author:  [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)
#### Post date: [December 28, 2020, 10:07am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/11 "2020-12-28T10:07:25Z")
> [@MA\_Laforge](#):
>
> I thought everything under the SciML umbrella was for machine learning.
> [@ChrisRackauckas](#):
>
> Everything is about scientific simulators, like [n-body simulators](https://github.com/SciML/NBodySimulator.jl/pull/37), [PDE solvers](https://github.com/SciML/NeuralPDE.jl), [quadrature methods](https://github.com/SciML/Quadrature.jl), and more which all allow for incorporating differentiable programming and machine learning into and around the core methods.
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 😃). 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](https://github.com/SciML/DifferentialEquations.jl) 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:](https://github.com/SciML/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:](https://github.com/SciML/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.
---
### Author:  [@pepijndevos](https://discourse.julialang.org/u/pepijndevos)
#### Post date: [December 28, 2020, 10:59am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/12 "2020-12-28T10:59:47Z")
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?
My code so far:
```julia
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 28, 2020, 3:06pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/13 "2020-12-28T15:06:27Z")
This is off-topic, but I’m willing to discuss it. It should probably be kicked into a new thread though.
> [@carstenbauer](#):
>
> 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](https://github.com/SciML/DifferentialEquations.jl) 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
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.
> [@carstenbauer](#):
>
> [RecursiveArrayTools.jl:](https://github.com/SciML/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?
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?
> [@carstenbauer](#):
>
> [ExponentialUtilities.jl:](https://github.com/SciML/ExponentialUtilities.jl) “Utility functions for exponential integrators for the SciML **scientific machine learning** ecosystem”.
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.
> [@carstenbauer](#):
>
> 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).
See the roadmap:
> **[SciML: Open Source Software for Scientific Machine Learning](https://sciml.ai/roadmap/)**
>
> Open Source Software for Scientific Machine Learning
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!
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 28, 2020, 3:10pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/14 "2020-12-28T15:10:28Z")
> [@pepijndevos](#):
>
> ```julia
> # 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)
>
> ```
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?
> [@pepijndevos](#):
>
> 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.
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.
> [@pepijndevos](#):
>
> One solver even complained it could only be used on a DAE and not on an ODE.
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)](https://diffeq.sciml.ai/stable/solvers/dae_solve/#OrdinaryDiffEq.jl-(Mass-Matrix))
But that should change soon.
---
### Author:  [@pepijndevos](https://discourse.julialang.org/u/pepijndevos)
#### Post date: [December 28, 2020, 4:02pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/15 "2020-12-28T16:02:58Z")
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.
---
### Author:  [@MA\_Laforge](https://discourse.julialang.org/u/MA_Laforge)
#### Post date: [December 28, 2020, 9:58pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/16 "2020-12-28T21:58:16Z")
### Off topic
> [@ChrisRackauckas](#):
>
> This is off-topic, but I’m willing to discuss it. It should probably be kicked into a new thread though.
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.
> [@ChrisRackauckas](#):
>
> 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.
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
> [@ChrisRackauckas](#):
>
> The extra tools in: [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.
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
> [@ChrisRackauckas](#):
>
> so a DSL over ModelingToolkit would give a high performance SPICE (note: one is in development).
> [@ChrisRackauckas](#):
>
> 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.
Very excited to hear this!
### My pre-conceptions
> [@ChrisRackauckas](#):
>
> ? 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](https://github.com/SciML/NBodySimulator.jl/pull/37), [PDE solvers](https://github.com/SciML/NeuralPDE.jl), [quadrature methods](https://github.com/SciML/Quadrature.jl), 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/](https://sciml.ai/).
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 🙂 .
### 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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 29, 2020, 1:13am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/17 "2020-12-29T01:13:51Z")
> [@MA\_Laforge](#):
>
> 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 🙂 .
“MLDifferentialEquations.jl”, or helper functions for using DifferentialEquations.jl with ML libraries like Flux.jl is:
> **[GitHub - SciML/DiffEqFlux.jl: Universal neural differential equations with...](https://github.com/SciML/DiffEqFlux.jl)**
>
> Universal neural differential equations with O(1) backprop, GPUs, and stiff+non-stiff DE solvers, demonstrating scientific machine learning (SciML) and physics-informed machine learning methods - G...
> [@MA\_Laforge](#):
>
> ### 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?
Scientific machine learning as a term comes from the DOE which put out a special report essentially defining the descipline:
> **[Workshop Report on Basic Research Needs for Scientific Machine Learning: Core...](https://www.osti.gov/biblio/1478744)**
>
> The U.S. Department of Energy's Office of Scientific and Technical Information
At around the same time, we had been starting to put out blog posts and talks around the same topic:
> **[DiffEqFlux.jl – A Julia Library for Neural Differential Equations](https://julialang.org/blog/2019/01/fluxdiffeq/)**
>
> DiffEqFlux.jl – A Julia Library for Neural Differential Equations...
[ ](https://www.youtube.com/watch?v=FGfx8CQHdQA)
And actually we had been doing some of this since the start of “JuliaDiffEq”:
> **[NeuralNetDiffEq.jl: A Neural Network solver for ODEs](https://julialang.org/blog/2017/10/gsoc-NeuralNetDiffEq/)**
>
> NeuralNetDiffEq.jl: A Neural Network solver for ODEs | My GSoC 2017 project (https://summerofcode.withgoogle.com/projects/#5850956641075200) was to implement a package for Julia to solve Ordinary Differential Equations using Neural Networks....
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](https://cran.r-project.org/web/packages/diffeqr/index.html) and [in Python](https://pypi.org/project/diffeqpy/0.1/). It had [optimization libraries](https://galacticoptim.sciml.ai/dev/), [quadrature libraries](https://github.com/SciML/Quadrature.jl), etc. mixed with [libraries for reservoir computing](https://github.com/SciML/ReservoirComputing.jl), [physics-informed neural networks](https://neuralpde.sciml.ai/dev/), 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.
> [@MA\_Laforge](#):
>
> Was it that otherwise, a layered solution became hacky and unmanageable?
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:
[ ](https://www.youtube.com/watch?v=AIX8ardhAxA)
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](https://github.com/SciML/DifferentialEquations.jl/issues/47), 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](https://github.com/SciML/SciMLBenchmarks.jl) 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](https://arxiv.org/abs/2010.04004), 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.
---
### Author:  [@pepijndevos](https://discourse.julialang.org/u/pepijndevos)
#### Post date: [December 29, 2020, 8:59pm UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/18 "2020-12-29T20:59:05Z")
Success??
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.
```julia
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")
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 30, 2020, 12:10am UTC](https://discourse.julialang.org/t/simulating-mosfet-equations-with-modelingtoolkit/52408/19 "2020-12-30T00:10:04Z")
> [@pepijndevos](#):
>
> 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.
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.