Elegant handling of possibly-state-dependent parameters in physical simulation

General Usage
,
1

For a couple years now, I have been using OrdinaryDiffEq for developing a handful of models for describing the process of pharmaceutical lyophilization. The more time I spend, the more I learn about more aspects of Julia or packages that would have made my life easier if I had known about them sooner (e.g. Parameters.jl, ComponentArrays.jl). In that light, I am wondering if my current headache already has a good solution.

These models have anywhere between 10 and 30 physical properties that need to be defined for simulation. These range from universal physical properties (e.g. universal gas constant) to things that don’t often need tinkering (e.g. density of glass) to things that are different for every simulation (input temperature/pressure profiles). As I compare these models to experiment and refine them, I periodically find that I need to add temperature-dependence (or pressure-dependence, or frequency-dependence) to a parameter that I was previously treating as constant.

Up til now, my approach has been to rewrite the code in such a way that, instead of taking a constant value, the ODE function takes a callable (or anonymous function) which it evaluates everywhere it needs the property. This works, of course. But every time it happens I inevitably end up re-running the code several times finding all the places where I am now adding an anonymous function to a number, and I feel like it clutters the input script, e.g. turning

p = 10.0u"Pa"

into

p = t->10.0u"Pa"

which, in addition to being a little uglier, also makes error messages not as nice and confuses the non-Julia-users when I need to hand off my code to someone else.

I have also experimented with creating a custom type with pretty printing, with a single-argument constructor that essentially acts as a nicer-looking version of the anonymous function.

julia> p = RampedVariable(10.0u"Pa")
RampedVariable(10.0u"Pa")
julia> p(50u"s")
10.0u"Pa"
julia> p(Inf)
10.0u"Pa"

I implemented that type as below, for any who are curious. But I find that this method still ends up being more typing, and means I need to document more interfaces myself, etc.

My question is: is there a more elegant solution to this problem than, say, implementing every single parameter every time as a callable which may not actually depend on any of its inputs?

For reference, sometimes I am thinking of time-dependent external controls, weakly temperature-dependent parameters like thermal conductivity, or things like the dielectric loss coefficient which varies with temperature and electric field frequency.

struct RampedVariable
    setpts::AbstractVector
    ramprates::AbstractVector
    holds::AbstractVector
    timestops::AbstractVector
end

function RampedVariable(initial)
    RampedVariable([initial], [], [], [0])
end

function RampedVariable(setpts, ramprate)

    if length(ramprate) == 0 || length(setpts) == 1
        @error "If no ramp necessary, construct RampedVariable with only one argument." ramprate
    end
    if length(ramprate) >= 2 || length(setpts) > 2
        @error "For multiple ramps, need at least one hold time. Construct RampedVariable with three arguments." ramprate
    end
    if length(setpts) != 2
        @error "Number of set points should be 1 more than ramps, since initial is included"
    end
    timestops = fill(0.0*setpts[1]/ramprate[1], 2)
    timestops[2] = timestops[1] + (setpts[2]-setpts[1])/ramprate
    RampedVariable(setpts, [ramprate], [], timestops)
end


function RampedVariable(setpts, ramprates, holds)

    if (length(ramprates) != length(holds) + 1 ) || (length(ramprates)==0)
        @error "Number of ramps should be zero or number of holds + 1"
    end
    if length(setpts) != length(ramprates) + 1
        @error "Number of set points should be 1 more than ramps, since initial is included"
    end
    timestops = fill(0.0*setpts[1]/ramprates[1], length(ramprates) + length(holds) + 1)
    (ramp, rest) = Iterators.peel(ramprates)
    timestops[2] = timestops[1] + (setpts[2]-setpts[1])/ramp
    for (i, ramp) in enumerate(rest)
        timestops[2i+1] = timestops[2i] + holds[i]
        timestops[2i+2] = timestops[2i+1] + (setpts[i+2]-setpts[i+1])/ramp
    end
    RampedVariable(setpts, ramprates, holds, timestops)
end

function (rv::RampedVariable)(t)
    if length(rv.timestops) == 1
        return rv.setpts[1]
    end

    im = findlast(rv.timestops .<= t)
    if im == length(rv.timestops)
        return rv.setpts[end]
    elseif iseven(im)
        return rv.setpts[im÷2+1]
    else
        ip = im+1
        return (rv.setpts[ip÷2+1] - rv.setpts[ip÷2])/(rv.timestops[ip] - rv.timestops[im])*(t - rv.timestops[im]) + rv.setpts[ip÷2]
    end
end
2

The way I would approach this would be to use a modeling tool like ModelingToolkit, and treat what you call parameters as system inputs. When you want a constant parameter, you connect a constant source, and when you want something else, you change the input signal source.

These inputs can then be propagated through the model, so that you only define the source at a single point, and if you change it it propagates through the entire model.

Take this pendulum on a cart as an example
inverted_cartpole
It’s modeled in ModelingToolkit, the main system model Cartpole includes the component

motor = TranslationalModelica.Force(use_support = true)

which is an input component. I then simulate the system with a sinusoidal input by instantiating the cartpole and an input source, Blocks.Sine and connect the two, connecting the sine to the force input

    @components begin
        world = World()
        cartpole = Cartpole()
        input = Blocks.Cosine(frequency=1, amplitude=1)
    end
    @equations begin
        connect(input.output, :u, cartpole.motor.f)
    end

Later on, I instead connect an LQG controller to the same input instead of the sine and simulate that.

This example is rather large and elaborate, the documentation of MTK includes much simpler examples. When you use MTK, MTK serves as a modeling layer, and OrdinaryDiffEq takes care of the actual integration when you are done with your model.