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?

2

If I am understanding you correctly, it should just be:

ssu0 = [
  pk.ρ = (t) -> 0
  pk.P = P0
]
prob = SteadyStateProblem(pk, u0)
sol = solve(prob)
3

Thanks for the quick reply Chris, much appreciated!

I also thought to do something along those lines at first.
Unfortunately, there are several problems with your suggestion:

  1. pk.ρ is not a variable, so it can’t be called. How would I set this up elegantly?
  2. The SteadyStateProblem requires a guess\initial conditions for all variables (adding e.g. pk.Cₖ => [1, 2, 3, 4, 5, 6]) allows the problem to be created.
  3. When the above u0 is used (no pk.ρ and with pk.Cₖ), the retcode is Unstable and the solution is just u0. Also DifferentialEquations has to be explicitly imported for some reason for the default solver.

The P=P_0 constraint cannot simply be an initial condition, as there are solutions for any P which would most likely be the converged result.