Easily creating model variants

[cstjean]

Is your feature request related to a problem? Please describe.

Where I work, we deliver chemical plants. Each client's plant might be a bit different (eg. the pump type changes), but they're always connected in the same way.

I have my models like

@mtkmodel SchmidtPump begin
    ...
end

@mtkmodel ManualPump begin
   ...
end

and my client models

@mtkmodel ClientInSpain begin
    @components begin
        pump = SchmidtPump(eta=20)
    end
    @equations begin
        connect(...)
    end
end

@mtkmodel ClientInFrance begin
    @components begin
        pump = ManualPump(operator_strength=3)
    end
    @equations begin
        connect(...)
    end
end

Since the @equations block is always exactly the same, this is rather redundant and verbose. What I'd like is to be able to write a vanilla function like client_in_france = create_plant(pump=ManualPump(operator_strength=3), ...). This would also cover another important use case, that of figuring out which pump is best for a particular client:

best_pump, _ = findmax([simulate_cost(create_plant(pump=pump) for pump in available_pumps])

Is it possible to write create_plant in MTK? Feels like at the moment it isn't, but maybe I've missed something.

Describe the solution you’d like

I would like to be able to write ManualPump(operator_strength=3) without getting a keyword argument name not assigned error. Would it be possible to have a default name=nothing kwarg, and have MTK assign each component a name later (at building/compiling time)?

Feels like in MTK one can talk about "the plant's operating_room's pump", but one cannot talk about "a ManualPump of operator_strength=30" as a fungible object / platonic ideal.

This would also (help to) solve

@components begin
    pumps = [ManualPump(operator_strength=3), ManualPump(operator_strength=7)] 
end

which is another use case we've had.

[hersle]

This sounds somewhat familiar to what I am doing. I have cosmological models instead of chemical plants, and different "submodels" for gravity and particles instead of pump types. See e.g. the keyword arguments of ΛCDM(). The idea is that ΛCDM can be constructed with keyword arguments to replace the "default components" with other components (e.g. new theory of gravity). I am not perfectly happy with it, but it has worked for now. I don't know if it is possible with @mtkmodel. Anyway, it is one of the reasons I decided to shy away from it for now. I am not sure if this covers exactly what you want, though.

[ChrisRackauckas]

I'm not sure you can do it with the macro. The macro is defined for top-level use. Programmatic interactions should happen not on the macro but on the IR itself.

[cstjean]

How would that work? Perhaps I wasn't clear; my issue is not with @mtkmodel, which works beautifully as a DSL for communicating with engineers. My problem is that calling ManualPump() witout the name kwarg does not work. If I pirate the method to take out the error, I can achieve what I want:

@eval ModelingToolkit (m::Model)(args...; name=:unspecified, kw...) =   # ideally name=nothing, but right now that fails
      m.f(args...; name, kw...) # piracy

function create_plant(; pump=ManualPump(operator_strength=3))
      pump = Symbolics.rename(pump, :pump)
      ...   
      return System(...)
end

using The Non-DSL (non-@mtkmodel) Way of Defining an System to create my System (that's what you mean by IR?)

Even better would be to allow

@mtkmodel Plant begin
    @components begin
        pump = SchmidtPump(eta=20)
    end
    @equations begin
        connect(...)
    end
end

Plant(pump=ManualPump())    # currently not possible

but A) that's more complicated and B) it will also require a name-kwarg-free ManualPump() constructor.

The idea is that ΛCDM can be constructed with keyword arguments to replace the "default components" with other components (e.g. new theory of gravity).

That's exactly it! I see that your function massless_neutrinos(g; lmax = 6, name = :ν, kwargs...) also has a default name kwarg, as I am advocating for.

[baggepinnen]

Related discussion #2038
Modelica has a feature to mark a component as replaceable which indicates that it can be replaced with any other that has the same external interface

[AayushSabharwal]

MTK has the ability to substitute components. See https://docs.sciml.ai/ModelingToolkit/stable/API/model_building/#ModelingToolkit.substitute_component. The idea is that you have an @mtkmodel PumpInterface which defines the interface followed by all pumps, typically only the part that is used to connect to other models.

Then, the replaceable model becomes

@mtkmodel ClientModel begin
    @components begin
        pump = PumpInterface(eta=20)
    end
    @equations begin
        connect(...)
    end
end

And after instantiating the model

@named sys = ClientModel()

The pump can be swapped out

@named pump = SchmidtPump(eta=20)
spain_model = substitute_component(sys, [sys.pump => pump])

[cstjean]

I tried modifying @mtkmodel and ended up creating a separate macro @mtkbmodel. It is 95% reusing the @mtkmodel macro expansion code, but component handling is different. Demo/POC:

using MTKButter, DifferentialEquations, Test
using ModelingToolkit: t_nounits as t, D_nounits as D

@mtkbmodel LinearXProducer begin
    @variables begin
        x(t) = 0
    end
    @parameters begin
        rate = 1
    end
    @equations begin
        D(x) ~ rate
    end
end

@mtkbmodel Room begin
    @variables begin
        total_x(t)
    end
    @components begin
        producers = fill(LinearXProducer(), 4)  # @components RHS now supports arbitrary Julia code
    end
    @equations begin
        total_x ~ sum(p.x for p in producers)
    end
end

# Normal usage
@named room1 = Room()  # @mtkcompile also works, but then @set cannot be used afterwards
prob = ODEProblem(mtkcompile(room1), Dict(), (0, 100.0))

# Let's create model variations starting from `room1`, using Accessors.jl.
room2 = @set room1.producers_1.rate = 10     # use Accessors.@set to change a parameter
room2b = @set room1.producers[1].rate = 10   # equivalent; component indexing works correctly
room3 = @set room1.producers = [LinearXProducer()] # replace all the producers with a single one
room3b = Room(producers=[LinearXProducer()])  # equivalent; the constructor can specify components
room4 = @insert room1.producers[5] = LinearXProducer()  # add a fifth producer
#room5 = @set room1.producers.rate = 4   # TODO: change the `rate` of all producers at once

@mtkbmodel ExponentialXProducer begin
    @variables begin
        x(t) = 1
    end
    @equations begin
        D(x) ~ x
    end
end

room5 = @set room1.producers[2] = ExponentialXProducer() # Heterogeneous component vectors also work. 

# Confirm that the results are correct
prob = ODEProblem(mtkcompile(room5), Dict(), (0,5))
sol = solve(prob, Tsit5())
@test sol[room5.total_x][end] ≈ 5 + exp(5) + 5 + 5 rtol=0.01

Implementation-wise:

• @mtkbmodel introduces a new PreSystem type, that is convertible to System.
• Instead of storing the equations, PreSystem contains a closure-to-build-the-equations.
• MTKButter pirates the Model constructor to make some_model() work (even without the name kwarg).

Is it conceptually problematic anywhere? MTK is so huge that I have to assume that this is "nice to achieve A,B,C, but incompatible with features X,Y,Z". For that reason, I put it in a new package, MTKButter.jl, for the time being. For us it'll be useful as a clean DSL to our users.

In any case, the piracy sucks, so this is one thing I'd appreciate getting into MTK (name-less Model constructors), in one form or another.

[cstjean]

@AayushSabharwal I get the software engineering angle, but from my chemical engineer colleagues' perspective:

MTK has the ability to substitute components. See https://docs.sciml.ai/ModelingToolkit/stable/API/model_building/#ModelingToolkit.substitute_component. The idea is that you have an @mtkmodel PumpInterface which defines the interface followed by all pumps, typically only the part that is used to connect to other models.

Why do I need to create an Interface? It's a bummer that now we're losing the default pump = SchmidtPump(), and that ClientModel() now contains a meaningless PumpInterface().

@named pump = SchmidtPump(eta=20)
spain_model = substitute_component(sys, [sys.pump => pump])

Why can't I just write prob = ODEProblem(sys, [sys.pump=>SchimdtPump(eta=20)])?

Maybe the Accessors.jl stuff is a distraction. If I made a minimal PR for this, I would start with:

• name defaults to :unnamed_schmidtpump. Or :unnamed. Or nothing
• parse_components stops parsing the RHS of each equation.

That would not solve the Interface-need problem, but it would at least improve ergonomics.

[AayushSabharwal]

MTKButter is interesting. I'll have to look at it a bit closer, there might be some issues with adding it to core MTK but I do think that it is a step in the right direction.

Why do I need to create an Interface?

MTK needs to know the minimal set of features (ports/variables/parameters) the model should have to be able to swap it out cleanly. Replacing a component with a "bigger" component is easy - it has everything the old one had so stuff will still work. Replacing a component with something that adds a few things and removes a few things is much more difficult. How do we detect if the new component is missing something that would cause the model to be unsimplifiable or incorrect in some other way? I imagine there is scope for user error, and then they open an issue in MTK saying "your component substitution is broken" when in fact we just aren't able to automatically tell them that they did something wrong.

This is not to say we can't make this happen, but we can't make it happen without a concrete set of rules for when a substitution is valid and when it isn't.

Why can't I just write prob = ODEProblem(sys, [sys.pump=>SchimdtPump(eta=20)])?

Because changing a component changes the structure of the system, and it needs to be resimplified. Swapping out a component may affect the DAE index and drastically affect the simplification result. It's not just a "drop in replacement". Simplified systems also don't have the component structure - they're flattened out. We can't tell you which component each equation in the simplified system comes from.

There is actually research in this direction (being able to swap components post simplification, without resimplifying) and as before it isn't impossible but it isn't easy.

name defaults to :unnamed_schmidtpump. Or :unnamed. Or nothing

I highly doubt we will end up doing this. name used to be optional back in the olden days, and making it mandatory was a very conscious decision. If names are optional, people stop giving them and expect things to keep working when they don't.

parse_components stops parsing the RHS of each equation.

This is interesting. I'll have to look at it a little closer - could you clarify what benefit this has?

[AayushSabharwal]

I think I should mention this for clarity: improved swapping and swapping without resimplification aren't impossible even right now and might not be a whole lot of work to get running to your satisfaction. But I specifically want to avoid adding features "because we need it to work for this case". Features such as these should have well-defined semantics for when they work, and when they don't. Without this, we end up with hard-to-maintain functionality because we can't test or document it properly because if we could, those tests/docs would constitute well-defined semantics.

[cstjean]

Thank you for the answer.

Because changing a component changes the structure of the system, and it needs to be resimplified

What is the use-case for post-simplification substitution? I would be happy with pre-simplification substitution (which is the PreSystem strategy). It is a lot simpler.

The idea behind PreSystem is to minimally contain the information inside @mtkmodel, so that operations like substitutions are cleanly defined.

I highly doubt we will end up doing this. name used to be optional back in the olden days, and making it mandatory was a very conscious decision. If names are optional, people stop giving them and expect things to keep working when they don't.

Understood. How about allowing name=nothing/:unnamed in System, and erroring later if the name is not specified (at simplification time...?)

parse_components stops parsing the RHS of each equation.

This is interesting. I'll have to look at it a little closer - could you clarify what benefit this has?

The macro parsing code is a lot simpler, and the RHS becomes arbitrary Julia code, which is more general and matches user expectations better.

@components begin
    pumps = [SchmidtPump(eta=0.5), SpinPump()]  # currently fails
    pump_energy_sources = map(default_energy_source, pumps)
end