How to implement similar functionality as in Simulink/Simscape by ModelingToolkitStandardLibrary.jl

Specific Domains Modelling & Simulations
1

I try the Julia Modeling Toolkit as an alternative to Simulink/Simscape. I’m struggling a bit with the connection to the “classical” Julia code. I’m looking for something like a PS-Simulink converter block. The following picture shows my very simple example in Simulink. I’m trying to replicate a similar example in Julia, but I don’t know how to pass values between the logic part and the circuit.

Cir

The circuit is simulated using a fixed-step solver (Euler/Heun), and I assume that my logic runs in the same execution loop as the circuit. The Julia code is a slightly modified example of a RC circuit from the ModelingToolkitStandardLibrary.jl. I saw in the library RealInput, RealOutput blocks. Are they attended to this functionality? How can I use them in my example?

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant

R = 1.0
C = 1.0
V = 1.0
dt = 0.004
@variables t
@named resistor = Resistor(R = R)
@named capacitor = Capacitor(C = C)
@named source = Voltage()
@named constant = Constant(k = V)
@named ground = Ground()

rc_eqs = [connect(constant.output, source.V)
    connect(source.p, resistor.p)
    connect(resistor.n, capacitor.p)
    connect(capacitor.n, source.n, ground.g)]

@named rc_model = ODESystem(rc_eqs, t,
    systems = [resistor, capacitor, constant, source, ground])
sys = structural_simplify(rc_model)
prob = ODEProblem(sys, Pair[], (0, 10.0))
solver = init(prob, Euler(),dt =dt)
# Number of iterations
num_iterations = 1000

for i in 1:num_iterations
    step!(solver)
    u = solver.u
    #??????????????
    #??????????????
    #??????????????
    push!(data_array, u[1])
end
plot(data_array)
2

You can always create a registered function to provide numerical input to a modelingtoolkit model, like this:

sequence::Vector{Float64} = [9, 9.5, 10, 9.5, 9.0, 8.5, 8, 8.5]
"""
    stairs(t)

Generate a stair function. Width of each stair is `se.t_step` seconds.
"""
function stairs(t)
    index = floor(t/se.t_step)
    index = mod(index, length(sequence)) + 1
    sequence[Int64(index)]
end
@register_symbolic stairs(t)
3

These blocks are MTK connectors that connect signals together, i.e., they are causal, so they are indeed related to what you want to do.

The problem you describe has two parts

  • Part of the model is causal, part is acausal. The realInput / RealOutput are the correct connectors to go between these domains.
  • Part of the model is operating in discrete time, i.e., at fixed time steps.

MTK currently does not give you many options for modeling hybrid continuous/discrete time systems. The way you manually do it here, by stepping the integrator and performing computations on the state manually, is likely going to be quite fragile. MTK does not guarantee you what variables will be chosen as the state during the simulation, so you cannot know what solver.u means. Even if you figure out what this u means right now, any update to MTK can (and often do) change that.

The unfortunate answer is that MTK currently won’t allow you to model this kind of system in a nice way, but it’s being actively worked on so expect it to work within a perhaps a few months.

Depending on how complicated the “simulation logic” is, you could potentially model it using a discrete event.

4

Thank you for the answer!
The first part of your response is clear. I will try the RealInput /RealOutput

In the second part, I provide a more detailed explanation. My intention is to modify the inputs, for example, voltage, during the simulation execution. I would like to attempt something akin to “real-time” simulation. This is the reason why I stepped the integrator. I thought there might be a way to access all variables in my circuit, similar to the RC Circuit Model example, with the difference that the values are displayed after simulation execution.

sol = solve(prob, Tsit5())
plot(sol, idxs = [capacitor.v, resistor.i],

Is there a way to do this, or should I wait for the new functionality of MTK that you described?
The discrete events are probably not a good option for me because my simulated system contains a lot of such “subsystems.”

5

How exactly do you intend to do this? Is your input a function of the simulation state, of time t, or a function of some real-world measurement?

If you want to, e.g., implement a controller and simulate that, then this would have to either wait for the new features or use the existing discrete events. If you want to input data from a sensor you read values from in realtime, then you may find the function generate_control_function useful. This function returns a function on the form

\dot x = f(x, u, p, t)

Where x is the state and u is the input. You would then have to form a closed-loop dynamics

\dot x = g(x, p, t) = f(x, u(y), p, t)

where you compute your input u as a function of some external measurements y. The function g can finally be simulated with fixed time steps.

6

The input function is a function of some “real-world input.” I’ve attempted to modify my example. I have an RC circuit with a switch, as depicted in the picture.

image

The Switch In signal is controlled by some external input. By the external input, I mean a real-world input, such as pushing a button on the GUI. My idea is to use MTK for real-time simulation, but it is necessary to have the coupling of MTK with Julia and a third-party process (like a GUI or hardware I/O, etc.). The coupling of Julia with the third-party process is possible through a C interface, for example.

Returning to my example in Julia.

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
using ModelingToolkitStandardLibrary.Blocks: RealInput
using ModelingToolkitStandardLibrary.Blocks: Step
using BenchmarkTools

V =1
dt = 0.004
R1 = 0.536746341349443
C1 = 1.2644803335127015

@variables t


function IdealClosingSwitch(; name, Ron=0.1, Roff=0.13)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters Ron=Ron Roff=Roff
    @named s = RealInput()
    eqs = [
        v ~ ifelse(s > 0,
            i*Ron,
            i*Roff),
    ]
    extend(ODESystem(eqs, t, [], pars; name = name), [oneport,s])
end


@named constant = Constant(k = V)
@named ground = Ground()
@named source = Voltage()
@named resistor1 = Resistor(R = R1)
@named capacitor1 = Capacitor(C = C1, v = 0.0)
@named step = Step(height = 1, start_time = 0.4)
@named switch = IdealClosingSwitch()

rc_eqs = [connect(constant.output, source.V)
connect(step.output,switch.s)
connect(source.p,switch.p)
connect(switch.n,resistor1.p)
connect(resistor1.n, capacitor1.p)
connect(capacitor1.n, source.n, ground.g)]

system = [ switch, resistor1 ,capacitor1 , constant ,ground, source] 

@named rc_model = ODESystem(rc_eqs, t,
    systems = system)
sys = structural_simplify(rc_model)
prob = ODEProblem(sys, Pair[], (0, 1.0))

solver = solve(prob, ImplicitEuler(),dt=dt)

plot(solver, idxs = [switch.i, capacitor1.v, resistor1.v, capacitor1.i, switch.v],
    title = "RC Network Circuit",
       labels = ["Voltage" ])

I created a custom component called IdealClosingSwitch and then a circuit that corresponds to the picture above. I encountered the following problems:

  1. I don’t know how to create a custom component that has pin connections (OnePort) and RealInput/RealOutput. I’ve tried, but there is still something missing.

  2. In the case that my IdealClosingSwitch is working, how can I replace the Step(height = 1, start_time = 0.4) function with some variable from my script? I think sometthing like this:

.
.
.
step = 0.0
@named switch = IdealClosingSwitch()

rc_eqs = [connect(constant.output, source.V)
connect(step,switch.s)
connect(source.p,switch.p)
.
.
.
for i in 1:num_iterations
    if i>100
        step = 1.0
    else
        step = 0
    end
    step!(solver)
end 
.
.
.
  1. The last issue is with the solver. If I remove the RealInput() from IdealClosingSwitch and change the switching logic so that it depends on the simulation time only, I get the following warning for a larger difference between Ron and Roff. It causes the simulation to stop. The values should be Ron ≈ 0.01 ohms, Roff ≈ 10e6 ohms. How can I change settings of ImplicitEuler() solver that it survives such a big change in R value?
7

Are you comparing it against a fixed time step implicit euler? If so, then you might want to make it fully not adaptive

8

ImplicitEuler() with following settings is not adaptive, is it?

The warning is easily reproducible if the value of Roff is changed to 0.5

┌ Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=0.19999999999999987, and step error estimate = 144.83149945907888. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\ev\.julia\packages\SciMLBase\4XvOU\src\integrator_interface.jl:599

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: Constant
using ModelingToolkitStandardLibrary.Blocks: RealInput

V =1
dt = 0.004
@variables t

global s =0.1

function IdealClosingSwitch(; name, Ron=0.1, Roff=0.25)
    @named oneport = OnePort()
    @unpack v, i = oneport
    pars = @parameters Ron=Ron Roff=Roff
    eqs = [
        v ~ ifelse(t > 0.2,
            i*Ron,
            i*Roff),
    ]
    extend(ODESystem(eqs, t, [], pars; name = name), oneport)
end

R1 = 0.536746341349443
C1 = 1.2644803335127015
@named constant = Constant(k = V)
@named ground = Ground()
@named source = Voltage()
@named resistor1 = Resistor(R = R1)
@named capacitor1 = Capacitor(C = C1, v = 0.0)
@named switch = IdealClosingSwitch()

rc_eqs = [connect(constant.output, source.V)
connect(source.p,switch.p)
connect(switch.n,resistor1.p)
connect(resistor1.n, capacitor1.p)
connect(capacitor1.n, source.n, ground.g)]
system = [ switch, resistor1 ,capacitor1 , constant ,ground, source] 


@named rc_model = ODESystem(rc_eqs, t,
    systems = system)
sys = structural_simplify(rc_model)
prob = ODEProblem(sys, Pair[], (0, 1.0))
solver = solve(prob, ImplicitEuler(),dt=dt)

plot(solver, idxs = [switch.i, capacitor1.v, resistor1.v, capacitor1.i, switch.v],
    title = "RC Network Circuit",
       labels = ["switch.i"  "capacitor1.v" "resistor1.v" "capacitor1.i" "switch.v" ])
9

You didn’t set adaptive=false, or remove Jacobian reuse. sol = solve(prob,ImplicitEuler(nlsolve = NLNewton(always_new = true)), adaptive = false, dt = dt) gives you a fully non-adaptive form.

1 Like
10

This solver configuration works without warnings, but the results aren’t good. The change in the value of Roff causes peaks in current that are proportional to Roff.

plot_130
plot_133

11

Yes, disabling error tolerance controls is rarely the right solution, mostly showing that you can. If you need the solution to be more correct, the ifelse should be handled by a callback so that the discontinuity is exactly resolved.

12

Thank you for your response. I am currently examining the modeling toolkit, particularly the ModelingToolkitStandardLibrary.jl, for simulating electrical circuits. I am considering its use as an alternative to Simulink Simscape. However, it appears to be less straightforward than I initially anticipated. My original idea was to define a few custom components (such as switches and diodes). For my specific use case, the precise simulation of these components is not crucial. I am content with a simulation similar to that in the following Modelica library: Modelica Electrical Analog Ideal.

The primary focus lies on the complexity of the simulated system, which consists of several thousand components, such as switches. I am uncertain whether defining numerous callbacks is a good approach. Currently, my mindset is somewhat influenced by the Simscape environment, where there is no need to define callbacks for a simple switch. I expected that the definition of the simulated system would be based on components (e.g., R, C, switches) and their interconnections.

It’s possible that I am mistaken, and Simscape automatically generates this additional logic (callbacks) in the background. Alternatively, are there other solvers that could be employed for these types of problems?

13

I don’t know how simscape works, but Modelica (another acausal modeling tool) does take the approach of automatically generate callbacks (“events”) for discontinuities. See for example