@@ -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
206287To 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
384465plot(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