BMS Current Control using Callbacks in ModelingToolKit

Specific Domains Modelling & Simulations
, , , ,
1

I want to implement a battery management system that simulates the internal behaviour of a cell (ECM), and controls the incoming current based on a voltage limit (here 3.8V). I’m trying to do this using ModelingToolkit.jl and callbacks, but i always get the same error: ERROR: KeyError: key cell1₊i(t) not found at prob = ODEProblem(sys, Pair[], tspan). Here’s the code:

using OrdinaryDiffEq
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using Plots: plot, plot!
using DataFrames
using CSV
using DataInterpolations
using LinearAlgebra
using CairoMakie
using ModelingToolkitStandardLibrary.Electrical
using ModelingToolkitStandardLibrary.Blocks
using DiffEqCallbacks

soc_vec = range(0, 1.0; length=1000)
ocv_vec = range(2.5, 4.2; length=1000)
focv = LinearInterpolation(ocv_vec, soc_vec, extrapolation=ExtrapolationType.Constant)

t_vec = range(0.0, 24*3600; length=1000)
i_vec = rand(-10.0:10.0, 1000)
fi = ConstantInterpolation(i_vec, t_vec)

@mtkmodel ECM begin
    @extend OnePort()
    @parameters begin
        Q
        R0
        R1
        τ1
        R2
        τ2
    end
    @structural_parameters begin
        focv
    end
    @variables begin
        vr(t)
        v1(t) = 0
        v2(t) = 0
        ocv(t)
        soc(t) = 0.5
        test(t) 
    end
    @equations begin
        D(soc) ~ i / (Q * 3600.0)
        D(v1) ~ -v1 / τ1 + i * (R1 / τ1)
        D(v2) ~ -v2 / τ2 + i * (R2 / τ2)
        vr ~ i * R0
        ocv ~ focv(soc)
        v ~ ocv + vr + v1 + v2
        test ~ ifelse(v > 3.8, 0, i)
    end
    @discrete_events begin #BMS
        (v >= 3.8) => [i ~ 0] 
        (v < 3.8) => [i ~ i]
    end
end

@mtkmodel Series_Connection begin
    @components begin
        cell1 = ECM(
            Q=4.8,
            R0=15e-3,
            R1=15e-3,
            τ1=60.0,
            R2=15e-3,
            τ2=600.0,
            focv=focv,
        )
        source = Current()
        ground = Ground()
    end
    @structural_parameters begin
        fi
    end
    @equations begin
        connect(cell1.p, source.p)
        connect(source.n, cell1.n, ground.g)
        source.i ~ -fi(t)
    end
  #= @discrete_events begin
        (cell1.v >= 3.8) => [source.i ~ 0] 
        (cell1.v < 3.8) => [source.i ~ -fi(t)]
    end=#
end

@mtkbuild sys = Series_Connection(; fi)
tspan = (0, 24 * 3600) 
prob = ODEProblem(sys, Pair[], tspan)
sol = solve(prob; saveat=1.0)

I also tried implementing the callback in the Series_Connection block but I get the same error.
I’d rather use @continuous_events but I can’t find a way to set the conditional v>=3.8V other than with @discrete_events.
In general, I need the current of the ECM model to behave as the test variable. Any ideas on how can I make it work?

2

There’s multiple parts to this. We should do better on this front (and MTKv10 solves exactly this problem) but the ERROR: KeyError: key cell1₊i(t) not found you’re getting is because MTK prior to v10 can only update variables that are part of unknowns(sys) in callbacks. Since cell.i isn’t part of unknowns(sys), code generation errors. I agree, the error message is arcane and should have been handled better.

The way to handle discontinuities like ifelse in your model is by turning the condition of the ifelse into a callback. Something like this:

@mtkmodel ECM begin
    @extend OnePort()
    @parameters begin
        Q
        R0
        R1
        τ1
        R2
        τ2
        cond(t)
    end
    @structural_parameters begin
        focv
    end
    @variables begin
        vr(t)
        v1(t) = 0
        v2(t) = 0
        ocv(t)
        soc(t) = 0.5
        test(t) 
    end
    @equations begin
        D(soc) ~ i / (Q * 3600.0)
        D(v1) ~ -v1 / τ1 + i * (R1 / τ1)
        D(v2) ~ -v2 / τ2 + i * (R2 / τ2)
        vr ~ i * R0
        ocv ~ focv(soc)
        v ~ ocv + vr + v1 + v2
        test ~ ifelse(cond, 0, i)
    end
    @continuous_events begin
        [v ~ 3.8] => # ...
    end
end

Where #... is replaced with a ModelingToolkit.FunctionalAffect (or, even better, ModelingToolkit.ImperativeAffect) that sets cond to true on an upcrossing and false on a downcrossing using affect! and affect_neg!.

ModelingToolkit has a simplification pass called IfLifting that can do this automatically. It is not enabled by default because it is still experimental with some known bugs. Most notably, it doesn’t always work if the discontinuity (the point where the ifelse changes) is at the initial time (tspan[1]). You can try out this pass using

@mtkbuild sys = Series_Connection(; fi) additional_passes = [ModelingToolkit.IfLifting]

Hopefully this helps. Let me know if you have further questions.

1 Like
3

Thank you for your fast reply! I tried implementing the ImperativeAffect parameter:

@mtkmodel ECM begin
    @parameters begin
       # ...
        condit(t)
    end
    @equations begin
       # ...
        i ~ fi(t) * condit
    end
    @continuous_events begin
        [v ~ 3.8] => ModelingToolkit.ImperativeAffect(modified = (; 
            i)) do x, o, i, c
            @set! condit(t) = 0
        end
    end
end
#...

But I get the following:
Warning: Expression i(t) cannot be assigned to; currently only unknowns and parameters may be updated by an affect. and then ERROR: Modified equation i(t) in affect refers to missing variable(s) Any[i(t)]; the variables may not have been added (e.g. if a component is missing) or they may have been reduced away.

I also tried taking out the @mtkmodel and modeling as in this Heater example:

@parameters Q = 4.8 R0 = 15e-3 R1 = 15e-3 τ1 = 60.0 R2 = 15e-3 τ2 = 600.0 
@parameters focv::Function fi::Function bms_on(t)::Bool=true
@variables soc(t) v1(t) = 0.0 v2(t) = 0.0 v(t) vr(t) ocv(t) i(t)

#...

using Setfield
bms_disable = ModelingToolkit.SymbolicContinuousCallback(
    [v ~ 3.8],
    ModelingToolkit.ImperativeAffect(modified = (; bms_on)) do x, o, c, i
        @set! x.bms_on = false
    end)
bms_enable = ModelingToolkit.SymbolicContinuousCallback(
    [v ~ 2.5],
    ModelingToolkit.ImperativeAffect(modified = (; bms_on)) do x, o, c, i
        @set! x.bms_on = true
    end)

eqs = [
    D(soc) ~ i / (Q * 3600.0),
    D(v1) ~ -v1 / τ1 + i * (R1 / τ1),
    D(v2) ~ -v2 / τ2 + i * (R2 / τ2),
    vr ~ i * R0,
    ocv ~ focv(soc),
    i ~ fi(t) * bms_on,
    v ~ ocv + vr + v1 + v2,
]
@named ecm = ODESystem(
    eqs, t, [soc, v1, v2, v, vr, ocv, i], [Q, R0, R1, τ1, R2, τ2, bms_on]; continuous_events = [bms_disable, bms_enable])
ss = structural_simplify(ecm)

prob = ODEProblem(ss, [ss.soc => 0.5, bms_on => true], (0.0, 10.0))
sol = solve(ode, Tsit5())
plot(sol)
plot(sol.t / 3600, sol[ss.v], axis=(; xlabel="Time in h", ylabel="Voltage in V"))

And it manages to solve the problem, but when trying to plot ss.v following comes:
ERROR: Invalid indexing operation: tried to access object of type ODESolution{...} (which is a parameter timeseries object) with variables having mixed timeseries indexes Any[SymbolicIndexingInterface.ContinuousTimeseries(), 1].