Custom Component

In this tutorial, the creation of a custom component is demonstrated via the Chua's circuit. The circuit is a simple circuit that shows chaotic behavior. Except for a non-linear resistor, every other component already is part of ModelingToolkitStandardLibrary.Electrical.

First, we need to make some imports.

using ModelingToolkit
using ModelingToolkit: t_nounits as t
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Electrical: OnePort
using OrdinaryDiffEq
using IfElse: ifelse
using Plots

Custom Component

Now the custom component can be defined. The Modelica implementation of the NonlinearResistor looks as follows:

model NonlinearResistor "Chua's resistor"
  extends Interfaces.OnePort;

  parameter SI.Conductance Ga "conductance in inner voltage range";
  parameter SI.Conductance Gb "conductance in outer voltage range";
  parameter SI.Voltage Ve "inner voltage range limit";
equation
  i = if (v < -Ve) then Gb*(v + Ve) - Ga*Ve else if (v > Ve) then Gb*(v - Ve) + Ga*Ve else Ga*v;
end NonlinearResistor;

this can almost be directly translated to the syntax of ModelingToolkit.

function NonlinearResistor(; name, Ga, Gb, Ve)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters Ga=Ga Gb=Gb Ve=Ve
    eqs = [
        i ~ ifelse(v < -Ve,
        Gb * (v + Ve) - Ga * Ve,
        ifelse(v > Ve,
            Gb * (v - Ve) + Ga * Ve,
            Ga * v))
    ]
    extend(ODESystem(eqs, t, [], pars; name = name), oneport)
end

Explanation

All components in ModelingToolkit are created via a function that serves as the constructor and returns some form of system, in this case, an ODESystem. Since the non-linear resistor is essentially a standard electrical component with two ports, we can extend from the OnePort component of the library.

@named oneport = OnePort()

This creates a OnePort with the name = :oneport. For easier notation, we can unpack the states of the component

@unpack v, i = oneport

It might be a good idea to create parameters for the constants of the NonlinearResistor.

pars = @parameters Ga=Ga Gb=Gb Ve=Ve

The syntax looks funny but it simply creates symbolic parameters with the name Ga where its default value is set from the function's argument Ga. While this is not strictly necessary it allows the user to remake the problem easily with different parameters or allow for auto-tuning or parameter optimization without having to do all the costly steps that may be involved with building and simplifying a model. The non-linear (in this case piece-wise constant) equation for the current can be implemented using IfElse.ifelse. Finally, the created oneport component is extended with the created equations and parameters. In this case, no extra state variables are added, hence an empty vector is supplied. The independent variable t needs to be supplied as the second argument.

extend(ODESystem(eqs, t, [], pars; name = name), oneport)

Building the Model

The final model can now be created with the components from the library and the new custom component.

@named L = Inductor(L = 18)
@named Ro = Resistor(R = 12.5e-3)
@named G = Conductor(G = 0.565)
@named C1 = Capacitor(C = 10, v = 4)
@named C2 = Capacitor(C = 100)
@named Nr = NonlinearResistor(Ga = -0.757576,
    Gb = -0.409091,
    Ve = 1)
@named Gnd = Ground()

connections = [connect(L.p, G.p)
               connect(G.n, Nr.p)
               connect(Nr.n, Gnd.g)
               connect(C1.p, G.n)
               connect(L.n, Ro.p)
               connect(G.p, C2.p)
               connect(C1.n, Gnd.g)
               connect(C2.n, Gnd.g)
               connect(Ro.n, Gnd.g)]

@named model = ODESystem(connections, t, systems = [L, Ro, G, C1, C2, Nr, Gnd])

Simulating the Model

Now the model can be simulated. First, structural_simplify is called on the model and an ODEProblem is built from the result. Since the initial voltage of the first capacitor was already specified via v, no initial condition is given and an empty pair is supplied.

sys = structural_simplify(model)
prob = ODEProblem(sys, Pair[], (0, 5e4), saveat = 0.01)
sol = solve(prob, Rodas4())

Plots.plot(sol[C1.v], sol[C2.v], title = "Chaotic Attractor", label = "",
    ylabel = "C1 Voltage in V", xlabel = "C2 Voltage in V")
Plots.savefig("chua_phase_plane.png")

Plots.plot(sol; idxs = [C1.v, C2.v, L.i],
    labels = ["C1 Voltage in V" "C2 Voltage in V" "Inductor Current in A"])
Plots.savefig("chua.png")