# Computing Steady State With Known Constraints
**URL:** https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489
**Category:** Modelling & Simulations
**Tags:** modelingtoolkit
**Created:** [September 18, 2025, 12:32pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489 "2025-09-18T12:32:57Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Aviv\_Barnea](https://discourse.julialang.org/u/Aviv_Barnea)
#### Post date: [September 18, 2025, 12:32pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/1 "2025-09-18T12:32:57Z")
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:
```julia
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?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 18, 2025, 11:39pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/2 "2025-09-18T23:39:38Z")
If I am understanding you correctly, it should just be:
```julia-auto
ssu0 = [
pk.ρ = (t) -> 0
pk.P = P0
]
prob = SteadyStateProblem(pk, u0)
sol = solve(prob)
```
---
### Author:  [@Aviv\_Barnea](https://discourse.julialang.org/u/Aviv_Barnea)
#### Post date: [September 19, 2025, 9:10am UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/3 "2025-09-19T09:10:00Z")
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.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 17, 2025, 12:29pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/4 "2025-10-17T12:29:19Z")
Did you ever open up an issue with this? I just found this in an old tab and wondered if it was solved.
---
### Author:  [@Aviv\_Barnea](https://discourse.julialang.org/u/Aviv_Barnea)
#### Post date: [October 18, 2025, 4:32pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/5 "2025-10-18T16:32:57Z")
Thanks for following up! I have not opened an issue with this (is the right repository ModelingToolkit.jl?)
I will open an issue, but does that mean this is not the expected way to work with MTK?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 19, 2025, 6:51pm UTC](https://discourse.julialang.org/t/computing-steady-state-with-known-constraints/132489/6 "2025-10-19T18:51:10Z")
It’s hard to tell what the issue is but if you open up the issue in a way where you resummarize where you are at then I might have a better answer, but tbh I am not completely following what the ask is right now and it would likely be easier if given a second try and added to my response queue