Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Initialization systems don't support parametric guesses #2716

Closed
MarcBerliner opened this issue May 15, 2024 · 1 comment · Fixed by #2775
Closed

Initialization systems don't support parametric guesses #2716

MarcBerliner opened this issue May 15, 2024 · 1 comment · Fixed by #2775
Assignees
Labels
bug Something isn't working

Comments

@MarcBerliner
Copy link
Contributor

MarcBerliner commented May 15, 2024

For example,

using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t
D = Differential(t)

@parameters a = -1.0
@variables x(t) [guess = a]

eqs = [D(x) ~ a]

initialization_eqs = [1 ~ exp(1 + x)]

@named sys = ODESystem(eqs, t; initialization_eqs)
sys = complete(structural_simplify(sys))

tspan = (0.0, 0.2)
prob = ODEProblem(sys, [], tspan, [])

@show prob.f.initializeprob.u0

sol = solve(prob.f.initializeprob; show_trace=Val(true))

the exact solution to the initialization_eqs is x(0) = -1. If we show the trace, you can see that it begins with an initial value of x(0) = 0

prob.f.initializeprob.u0 = [0.0]

Algorithm: NewtonRaphson(
   descent = NewtonDescent()
)

----     -------------        -----------         
Iter     f(u) inf-norm        Step 2-norm         
----     -------------        -----------         
0        1.71828183e+00       6.46802053e-314     
...
6        0.00000000e+00       1.22346577e-12      
Final    0.00000000e+00      
----------------------      
retcode: Success
u: 1-element Vector{Float64}:
 -1.0

When we change the guess for x(t) to a Float64, i.e., @variables x(t) [guess = -1.0], it respects the guess:

prob.f.initializeprob.u0 = [-1.0]

Algorithm: NewtonRaphson(
   descent = NewtonDescent()
)

----     -------------        -----------         
Iter     f(u) inf-norm        Step 2-norm         
----     -------------        -----------         
0        0.00000000e+00       4.94065646e-324     
1        0.00000000e+00       0.00000000e+00      
Final    0.00000000e+00      
----------------------      
retcode: Success
u: 1-element Vector{Float64}:
 -1.0

cc: @AayushSabharwal

@MarcBerliner MarcBerliner added the bug Something isn't working label May 15, 2024
@AayushSabharwal
Copy link
Member

Interestingly, a defaults to 0 in the initialization system:

julia> prob.f.initializeprob.f.sys
Model sys with 1 equations
Unknowns (1):
  x(t) [defaults to a]
Parameters (2):
  a [defaults to 0.0]
  t

This is because schedule.dummy_sub contains D(x) => a and as a result it sets the default value for a to be the default dummy derivative guess. I'm confused as to why this happens.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants