Computing Steady State With Known Constraints

Specific Domains Modelling & Simulations
1

Hi, new to ModelingToolkit.jl, and trying to port some models in Nuclear Reactor computations. So here’s the simplest model (point reactor) I wrote:

using ModelingToolkit, Plots, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
using LinearAlgebra

"""
    dP/dt = P(ρ - β)/Λ + Cₖ⋅λₖ
    dCₖ/dt = -Cₖλₖ + P(βₖ/Λ)

Where:
    P: power
    ρ: reactivity
    β: total delayed neutron fraction (1\$)
    Λ: generation time
    Cₖ: k-group's contribution to power
    λₖ: k-group's decay rate.
"""
@mtkmodel PointKinetics begin
    begin
        β = sum(βₖ)
        n = length(βₖ)
    end
    @structural_parameters begin 
        ρ
    end
    @parameters begin
        Λ
        βₖ[1:n]
        λₖ[1:n]
    end
    @variables begin
        P(t)
        Cₖ(t)[1:n]
    end
    @equations begin
        D(P) ~ (ρ(t) - β)P/Λ + λₖ ⋅ Cₖ
        D(Cₖ) ~ -λₖ .* Cₖ .+ P * βₖ / Λ
    end
end

# ======== Example usage ========
ρin(t) = 1e-3*sin(t)
@mtkcompile pk = PointKinetics(
    βₖ = [0.00025, 0.0012, 0.00125, 0.0026, 0.00075, 0.00027],
    Λ = 1e-3,
    λₖ = [55.72, 22.72, 6.22, 2.3, 0.618, 0.23],
    ρ = ρin,
)
u0 = [
    pk.P => 10,
    pk.Cₖ => [1, 2, 3, 4, 5, 6]
]
prob = ODEProblem(pk, u0, (0, 8))
sol = solve(prob)

This works well.
Now I want to be able to compute the steady state using this model - where I know the initial power (a solution exists for any power, provided reacitivity is zero, and is easy to see).
I want to tell a solver to solve this problem under those constraints - \rho=0, P=P_0, but cannot figure out how to do it.

Of course, this example is trivial, but I would like to use MTK for much more difficult systems. Any thoughts?