ModelingToolkit for nonlinear circuits and interface with fixed time-step simulators

Specific Domains Modelling & Simulations
1

mtk_for_nonlin_circuit
mtk_for_nonlin_circuit1044×502 8.8 KB

Hi all,
I am relatively new to Julia and believe MTK offers an efficient way to model nonlinear circuits/communication channels for an existing link simulation framework I am building.

As shown in the picture, most sub-blocks in the system is fixed time-step based and creates oversampled time-domain waveforms. It operates on a block basis, and every for-loop iteration processes some signal in batch.

My questions regarding MTK involve the following

  1. I believe it can fit well into the fixed time-step framework by using saveat, but I am not sure if it will significantly impact the performance of the solver. Can anyone comment on this?

  2. I couldn’t find examples showing how to provide input vectors to an acasually modelled circuit (for example, a simple RC circuit with arbitrary input from somewhere else). Most examples I found involve just constant voltages and getting the step response of RC circuits. Is it done through something like RealInput?

  3. To fit within the block-processing framework, the states of the ODE from the previous block needs to be used as the initial condition of the current block. Is there a way to solve the same ODE problem with continuously updated initial conditions without redefining the problem?

  4. Lastly, can one define a nonlinear component with any arbitrary function as shown in the picture? The example I found in the documentation uses piecewise linear functions, but still numbers as parameters. I guess I could extract parameters and define the function internally to the component, but was just thinking perhaps it’s the most general to just pass in a function.

I apologize if this is a bloated question, but all these seem related within the simulation framework I am working with. I really appreciate your help! Truly having a blast with Julia so far.

2

saveat doesn’t change the stepping, just the saving. You can use adaptive=false on solvers to force it though. Generally not recommended, instead you can handle the control inputs in callbacks.

Yes, you just make a component that connects to it.

Callbacks

Yes, you can either define it symbolically or register a function for it.

3

Thanks Chris!

I tried using the TimeVaryingFunction component and Linear_Interpolation to pass in external data vectors as inputs in the following example, and it seems to work. My goal now is to wrap this in a for-loop so that I can supply new data vectors in each iteration.

using ModelingToolkit, OrdinaryDiffEq, Plots, DataInterpolations
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks: TimeVaryingFunction
using ModelingToolkit: t_nounits as t

dt = 1/10e9/32
tend = (32*100)*dt
tt = 0:dt:tend
vin = zeros(length(tt))
f_itp_vin = LinearInterpolation(vin, tt)

R = 1e3
C = 1.0e-12
V = 1.0
@named resistor = Resistor(R = R)
@named capacitor = Capacitor(C = C, v = 0.0)
@named source = Voltage()
@named driver = TimeVaryingFunction(f_itp_vin)
@named ground = Ground()

rc_eqs = [connect(driver.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, driver, source, ground])
sys = structural_simplify(rc_model)
prob = ODEProblem(sys, Pair[], (0, tend), saveat=dt)

lastval = vin[end]

for _ in 0:nblks
  #new data will be generated by some other function
  vin = randn(length(tt))
  f_itp_vin = LinearInterpolation(vin, tt)

  #######
  # How to update the interpolation function in the TimeVaryingFunction 
  # component in the existing ODE problem?
  ######

  #remake the problem with updated initial condition from previous iteration and new interpolated data
  prob = remake(prob, u0=[capacitor.v=>lastval) 
  sol = solve(prob, Tsit5())
  lastval = sol[end][1]
end

Is it possible to update f_itp_vin? I am not sure what the best way to achieve this is. The data vector (driver’s output in this example) will be generated by some other functions in each for-loop iteration. Does callback help in this case, knowing that we are inside a for-loop?

On a side note, seems the f_itp_vin function is immutable. For my use case, it might be more efficient to simply update the interpolator with new data (something like f_ipt_vin.u = vin) every time new vin is generated. I have a feeling this capability might be handy when the same interpolation method is used just for different data.

4

Check out

5

Thanks Fredrik. This is very helpful.

After experimenting, I could now update the force data vector continuously with the Ref approach, but seems it needs to reference a global variable rdata and dt. Is there anyway to still use the TimeVaryingFunction component and somehow pass a get_sampled_data(t; dt, rdata) function with dt and rdata set?

I also tried SampledData component but got stuck at the remaking ODE problem part.

My update_prob function looks like below

function update_prob(prob, u_last, data, dt; sys, defs)
    defs[sys.driver.buffer] = Parameter(data, dt)
    # ensure p is a uniform type of Vector{Parameter{Float64}} (converting from Vector{Any})
    u_new  = ModelingToolkit.varmap_to_vars([capacitor.v=>u_last], unknowns(sys))
    p = ModelingToolkit.MTKParameters(sys, defs, u_new; tofloat = false)
    remake(prob; p)
end

My for loop now looks like

sys = structural_simplify(rc_model)
defs = ModelingToolkit.defaults(sys)
prob = ODEProblem(sys_sim, [capacitor.v=>0.0], (0, tend), saveat=dt)

vc_last = 0.0
@time for n = 1:10

#dummy clock signal
vin = [vin_last; kron(kron(ones(Int(nsym/2)), [-1, 1]), ones(osr))] 

prob = update_prob(prob, vc_last, vin, dt; sys, defs)
sol = solve(prob, ImplicitEuler(), adaptive=false, dt=dt)
vc_last = sol[capacitor.v][end]
vin_last = vin[end]

#plotting to see how the boundary looks like between each iteration
lines!(ax, n*tend .+ tt[2:end], sol[capacitor.v][2:end])

end

I get the following error

caused by: MethodError: no method matching Parameter(::Parameter{Float64}, ::Float64, ::Bool)

Closest candidates are:
  Parameter(::SymbolicUtils.Symbolic{<:Vector}, ::Real, ::Bool)
   @ ModelingToolkitStandardLibrary C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\wrapper-types.jl:140
  Parameter(::Vector{T}, ::T, ::Bool) where T<:Real
   @ ModelingToolkitStandardLibrary C:\Users\Kevin Z\.julia\packages\ModelingToolkitStandardLibrary\OIKXb\src\Blocks\sources.jl:449
  Parameter(::Symbolics.Arr{<:Any, 1}, ::Real, ::Bool)
   @ ModelingToolkitStandardLibrary C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\wrapper-types.jl:140
  ...

Stacktrace:
  [1] fast_substitute(expr::SymbolicUtils.BasicSymbolic{Real}, subs::Dict{Any, Any}; operator::Type)
    @ Symbolics C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\variable.jl:489
  [2] fast_substitute
    @ C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\variable.jl:473 [inlined]
  [3] fixpoint_sub(x::SymbolicUtils.BasicSymbolic{Real}, dict::Dict{Any, Any}; operator::Type)
    @ Symbolics C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\variable.jl:429
  [4] fixpoint_sub(x::SymbolicUtils.BasicSymbolic{Real}, dict::Dict{Any, Any})
    @ Symbolics C:\Users\Kevin Z\.julia\packages\Symbolics\PAFGz\src\variable.jl:428
  [5] (::ModelingToolkit.var"#115#126"{Dict{Any, Any}})(::Pair{Any, Any})
    @ ModelingToolkit .\none:0
  [6] iterate
    @ .\generator.jl:47 [inlined]
  [7] grow_to!(dest::Dict{…}, itr::Base.Generator{…}, st::Int64)
    @ Base .\dict.jl:131
  [8] grow_to!(dest::Dict{Any, Any}, itr::Base.Generator{Dict{Any, Any}, ModelingToolkit.var"#115#126"{Dict{Any, Any}}})
    @ Base .\dict.jl:125
  [9] dict_with_eltype
    @ .\abstractdict.jl:590 [inlined]
 [10] Dict(kv::Base.Generator{Dict{Any, Any}, ModelingToolkit.var"#115#126"{Dict{Any, Any}}})
    @ Base .\dict.jl:109
 [11] ModelingToolkit.MTKParameters(sys::ODESystem, p::Dict{Any, Any}, u0::Vector{Float64}; tofloat::Bool, use_union::Bool)
    @ ModelingToolkit C:\Users\Kevin Z\.julia\packages\ModelingToolkit\IVr5X\src\systems\parameter_buffer.jl:94
 [12] update_prob(prob::ODEProblem{…}, u_last::Float64, data::Vector{…}, dt::Float64; sys::ODESystem, defs::Dict{…})
    @ Main c:\JuliaWS\mtk_test\mtk_test.jl:37
 [13] macro expansion
    @ c:\JuliaWS\mtk_test\mtk_test.jl:93 [inlined]
 [14] macro expansion
    @ .\timing.jl:279 [inlined]
 [15] top-level scope
    @ c:\JuliaWS\mtk_test\mtk_test.jl:90
Some type information was truncated. Use `show(err)` to see complete types.

Not sure I am regenerating the parameter maps correctly