Convergence Issues with Controllable Switch and Square Wave using MTK - iflifting

[junmorenodi99]

Context

I’m trying to simulate a switching circuit using an External Control Signal and a Controllable Switch. The Switch is designed as a non-ideal resistor where the resistance value varies from high to low in response to the Control Signal.

      +V  
        |  
       ┌┴┐     Control (Square Wave)
       │S│<─────●  
       └┬┘            
        |            
       [R]  ← Load resistor R  
        |            
       GND

Description

I tried a simple test circuit with two approaches for the Control Signal: a Sinusoidal Wave and a Square Wave. The Sinusoidal Wave works correctly and the system behaves as expected without any convergence issues. However, when the Square Wave is applied, there is a convergence problem stopping the execution very early (likely due to the nonlinearity induced by this signal). There isn't error messages prompted but two warnings :

┌ Warning: Rosenbrock methods on equations without differential states do not bound the error on interpolations.
└ @ OrdinaryDiffEqCore C:\Users\jnmor\.julia\packages\OrdinaryDiffEqCore\zs1s7\src\solve.jl:106

┌ Warning: At t=0.0, dt was forced below floating point epsilon 5.0e-324, and step error estimate = NaN. Aborting. There is either an error in your model specification or the true solution is unstable (or the true solution can not be represented in the precision of Float64).
└ @ SciMLBase C:\Users\jnmor\.julia\packages\SciMLBase\u2Ue2\src\integrator_interface.jl:623

As part of the suggestion to work around the issue, I tried to include the iflifting function, but when I include that into my code, I'm getting the error :

ERROR: BoundsError: attempt to access 0-element Vector{Int64} at index [1]

Expected behavior

Ideally, the model would converge, showing squared voltage on the Load resistor.
The main goal is to be able to model a Switching Circuit based on MTK.

Code

Here’s the code of the circuit model:

@mtkmodel R_Test begin
    @parameters begin
        RLoad_param = 10, [description = "Output Load (Ohm)"]
        Vin_param = 200, [description = "Vin Souce (Vdc)"]
        f = 100, [description = "Switch input Frequency (Hz)"]
    end

    @components begin
        switch_input = squareWave(frequency = f, amplitude = 1, offset = 0, dutyCycle = 0.3)
        #switch_input = Sine(frequency = f, amplitude = 0.99, offset = 1.0)
        voltage_input = Constant(k = Vin_param)
        switch = IdealSwitch(Ron = 1e-3, Roff = 1)
        source = Voltage()
        RLoad = Resistor(R = RLoad_param)
        ground = Ground()
    end

    @equations begin 
      connect(voltage_input.output, source.V)
      connect(switch_input.output, switch.conduction)
      connect(source.p, switch.p)
      connect(switch.n,RLoad.p)
      connect(source.n, RLoad.n, ground.g) 
    end
end

The Switch is defined as

@mtkmodel IdealSwitch begin
  @extend v, i = oneport = OnePort()  
  
  @parameters begin 
    Ron = 1e-3, [description = "On resistance magnitud"]
    Roff = 1e6, [description = "Off resistance magnitud"]
  end

  @components begin
    conduction = RealInput() 
  end

  @variables begin
    q(t), [description = "Control Signal "]
    R(t), [description = "Resistance "]
  end
  
  @equations begin
    q ~ conduction.u
    R ~  Ron + q * Roff
    v ~ i * R
  end
end

The Square Wave is defined as

@component function squareWave(;name, frequency,  amplitude = 1, offset = 0, dutyCycle = 0.5)
    @named output = RealOutput()
    pars = @parameters  amplitude=amplitude frequency=frequency dutyCycle=dutyCycle
    equation =  amplitude * mod1(t, 1/frequency) < dutyCycle/frequency
    eqs = output.u ~ equation
    compose(System(eqs, t, [], pars; name = name), [output])
end

and the solver configuration is

@named sys_dec = R_Test()
sys = structural_simplify(sys_dec, additional_passes = [IfLifting])
prob = ODEProblem(sys, Pair[], (0.0, 100e-3); guesses = [sys.switch.i => 5.0])
sol = solve(prob, Rodas5();
            reltol=1e-8,
            abstol=1e-10,
            dt=1e-7,        
            dtmax=1e-5,     
            saveat=1e-5)

Thanks for your help!