Skip to content

Commit c3487a5

Browse files
docs: document parameter initialization in the tutorial
1 parent 8d41eb8 commit c3487a5

File tree

1 file changed

+134
-0
lines changed

1 file changed

+134
-0
lines changed

docs/src/tutorials/initialization.md

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -201,6 +201,87 @@ long enough you will see that `λ = 0` is required for this equation, but since
201201
problem constructor. Additionally, any warning about not being fully determined can
202202
be suppressed via passing `warn_initialize_determined = false`.
203203

204+
## Initialization of parameters
205+
206+
Parameters may also be treated as unknowns in the initialization system. Doing so works
207+
almost identically to the standard case. For a parameter to be an initialization unknown
208+
(henceforth referred to as "solved parameter") it must represent a floating point number
209+
(have a `symtype` of `Real` or `<:AbstractFloat`) or an array of such numbers. Additionally,
210+
it must have a guess and one of the following conditions must be satisfied:
211+
212+
1. The value of the parameter as passed to `ODEProblem` is an expression involving other
213+
variables/parameters. For example, if `[p => 2q + x]` is passed to `ODEProblem`. In
214+
this case, `p ~ 2q + x` is used as an equation during initialization.
215+
2. The parameter has a default (and no value for it is given to `ODEProblem`, since
216+
that is condition 1). The default will be used as an equation during initialization.
217+
3. The parameter has a default of `missing`. If `ODEProblem` is given a value for this
218+
parameter, it is used as an equation during initialization (whether the value is an
219+
expression or not).
220+
4. `ODEProblem` is given a value of `missing` for the parameter. If the parameter has a
221+
default, it will be used as an equation during initialization.
222+
223+
All parameter dependencies (where the dependent parameter is a floating point number or
224+
array thereof) also become equations during initialization, and the dependent parameters
225+
become unknowns.
226+
227+
`remake` will reconstruct the initialization system and problem, given the new
228+
constraints provided to it. The new values will be combined with the original
229+
variable-value mapping provided to `ODEProblem` and used to construct the initialization
230+
problem.
231+
232+
### Parameter initialization by example
233+
234+
Consider the following system, where the sum of two unknowns is a constant parameter
235+
`total`.
236+
237+
```@example paraminit
238+
using ModelingToolkit, OrdinaryDiffEq # hidden
239+
using ModelingToolkit: t_nounits as t, D_nounits as D # hidden
240+
241+
@variables x(t) y(t)
242+
@parameters total
243+
@mtkbuild sys = ODESystem([D(x) ~ -x, total ~ x + y], t;
244+
defaults = [total => missing], guesses = [total => 1.0])
245+
```
246+
247+
Given any two of `x`, `y` and `total` we can determine the remaining variable.
248+
249+
```@example paraminit
250+
prob = ODEProblem(sys, [x => 1.0, y => 2.0], (0.0, 1.0))
251+
integ = init(prob, Tsit5())
252+
@assert integ.ps[total] ≈ 3.0 # hide
253+
integ.ps[total]
254+
```
255+
256+
Suppose we want to re-create this problem, but now solve for `x` given `total` and `y`:
257+
258+
```@example paraminit
259+
prob2 = remake(prob; u0 = [y => 1.0], p = [total => 4.0])
260+
initsys = prob2.f.initializeprob.f.sys
261+
```
262+
263+
The system is now overdetermined. In fact:
264+
265+
```@example paraminit
266+
[equations(initsys); observed(initsys)]
267+
```
268+
269+
The system can never be satisfied and will always lead to an `InitialFailure`. This is
270+
due to the aforementioned behavior of retaining the original variable-value mapping
271+
provided to `ODEProblem`. To fix this, we pass `x => nothing` to `remake` to remove its
272+
retained value.
273+
274+
```@example paraminit
275+
prob2 = remake(prob; u0 = [y => 1.0, x => nothing], p = [total => 4.0])
276+
initsys = prob2.f.initializeprob.f.sys
277+
```
278+
279+
The system is fully determined, and the equations are solvable.
280+
281+
```@example
282+
[equations(initsys); observed(initsys)]
283+
```
284+
204285
## Diving Deeper: Constructing the Initialization System
205286

206287
To get a better sense of the initialization system and to help debug it, you can construct
@@ -383,3 +464,56 @@ sol[α * x - β * x * y]
383464
```@example init
384465
plot(sol)
385466
```
467+
468+
## Solving for parameters during initialization
469+
470+
Sometimes, it is necessary to solve for a parameter during initialization. For example,
471+
given a spring-mass system we want to find the un-stretched length of the spring given
472+
that the initial condition of the system is its steady state.
473+
474+
```@example init
475+
using ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica: Fixed, Mass, Spring,
476+
Force, Damper
477+
using ModelingToolkitStandardLibrary.Blocks: Constant
478+
479+
@named mass = Mass(; m = 1.0, s = 1.0, v = 0.0, a = 0.0)
480+
@named fixed = Fixed(; s0 = 0.0)
481+
@named spring = Spring(; c = 2.0, s_rel0 = missing)
482+
@named gravity = Force()
483+
@named constant = Constant(; k = 9.81)
484+
@named damper = Damper(; d = 0.1)
485+
@mtkbuild sys = ODESystem(
486+
[connect(fixed.flange, spring.flange_a), connect(spring.flange_b, mass.flange_a),
487+
connect(mass.flange_a, gravity.flange), connect(constant.output, gravity.f),
488+
connect(fixed.flange, damper.flange_a), connect(damper.flange_b, mass.flange_a)],
489+
t;
490+
systems = [fixed, spring, mass, gravity, constant, damper],
491+
guesses = [spring.s_rel0 => 1.0])
492+
```
493+
494+
Note that we explicitly provide `s_rel0 = missing` to the spring. Parameters are only
495+
solved for during initialization if their value (either default, or explicitly passed
496+
to the `ODEProblem` constructor) is `missing`. We also need to provide a guess for the
497+
parameter.
498+
499+
If a parameter is not given a value of `missing`, and does not have a default or initial
500+
value, the `ODEProblem` constructor will throw an error. If the parameter _does_ have a
501+
value of `missing`, it must be given a guess.
502+
503+
```@example init
504+
prob = ODEProblem(sys, [], (0.0, 1.0))
505+
prob.ps[spring.s_rel0]
506+
```
507+
508+
Note that the value of the parameter in the problem is zero, similar to unknowns that
509+
are solved for during initialization.
510+
511+
```@example init
512+
integ = init(prob)
513+
integ.ps[spring.s_rel0]
514+
```
515+
516+
The un-stretched length of the spring is now correctly calculated. The same result can be
517+
achieved if `s_rel0 = missing` is omitted when constructing `spring`, and instead
518+
`spring.s_rel0 => missing` is passed to the `ODEProblem` constructor along with values
519+
of other parameters.

0 commit comments

Comments
 (0)