From e61770376e260b5f717aa940a708e15e13ab7217 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 19 Jul 2024 19:16:44 -0400 Subject: [PATCH 1/4] Allow for setting fully_determined=true in initialization --- src/systems/diffeqs/abstractodesystem.jl | 9 ++++++--- test/initializationsystem.jl | 18 ++++++++++++++++++ 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index ab0565eb31..f4268c506a 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -762,6 +762,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; warn_initialize_determined = true, build_initializeprob = true, initialization_eqs = [], + fully_determined = false, kwargs...) eqs = equations(sys) dvs = unknowns(sys) @@ -835,7 +836,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; end initializeprob = ModelingToolkit.InitializationProblem( sys, t, u0map, parammap; guesses, warn_initialize_determined, - initialization_eqs, eval_expression, eval_module) + initialization_eqs, eval_expression, eval_module, fully_determined) initializeprobmap = getu(initializeprob, unknowns(sys)) zerovars = Dict(setdiff(unknowns(sys), keys(defaults(sys))) .=> 0.0) @@ -1478,6 +1479,7 @@ InitializationProblem{iip}(sys::AbstractODESystem, u0map, tspan, simplify = false, linenumbers = true, parallel = SerialForm(), initialization_eqs = [], + fully_determined = false, kwargs...) where {iip} ``` @@ -1527,6 +1529,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, check_length = true, warn_initialize_determined = true, initialization_eqs = [], + fully_determined = false, kwargs...) where {iip, specialize} if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") @@ -1535,10 +1538,10 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, isys = get_initializesystem(sys; initialization_eqs) elseif isempty(u0map) && get_initializesystem(sys) === nothing isys = structural_simplify( - generate_initializesystem(sys; initialization_eqs); fully_determined = false) + generate_initializesystem(sys; initialization_eqs); fully_determined) else isys = structural_simplify( - generate_initializesystem(sys; u0map, initialization_eqs); fully_determined = false) + generate_initializesystem(sys; u0map, initialization_eqs); fully_determined) end uninit = setdiff(unknowns(sys), [unknowns(isys); getfield.(observed(isys), :lhs)]) diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index f4097ecd97..5111381121 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -17,6 +17,10 @@ sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test maximum(abs.(sol[conditions])) < 1e-14 +@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem(pend, 0.0, [], [g => 1]; + guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2], + fully_determined = true) + initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) @test initprob isa NonlinearProblem @@ -31,6 +35,10 @@ initprob = ModelingToolkit.InitializationProblem( sol = solve(initprob) @test !SciMLBase.successful_retcode(sol) +@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem( + pend, 0.0, [], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend), + fully_determined = true) + prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend)) prob.f.initializeprob isa NonlinearProblem @@ -47,6 +55,10 @@ sol = solve(prob.f.initializeprob) sol = solve(prob, Rodas5P()) @test maximum(abs.(sol[conditions][1])) < 1e-14 +@test_throws ModelingToolkit.ExtraVariablesSystemException ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], + guesses = ModelingToolkit.missing_variable_defaults(pend), + fully_determined = true) + @connector Port begin p(t) dm(t) = 0, [connect = Flow] @@ -225,6 +237,8 @@ initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) @test maximum(abs.(initsol[conditions])) < 1e-14 +@test_throws ModelingToolkit.ExtraEquationsSystemException ModelingToolkit.InitializationProblem(sys, 0.0, fully_determined = true) + allinit = unknowns(sys) .=> initsol[unknowns(sys)] prob = ODEProblem(sys, allinit, (0, 0.1)) sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit()) @@ -232,7 +246,11 @@ sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit()) @test sol.retcode == SciMLBase.ReturnCode.Unstable @test maximum(abs.(initsol[conditions][1])) < 1e-14 +prob = ODEProblem(sys, allinit, (0, 0.1)) prob = ODEProblem(sys, [], (0, 0.1), check = false) + +@test_throws ModelingToolkit.ExtraEquationsSystemException ODEProblem(sys, [], (0, 0.1), fully_determined = true) + sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure @test sol.retcode == SciMLBase.ReturnCode.Unstable From 5fe46e69b7363d0b897e1c91d15c1a2d207e7e7a Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 19 Jul 2024 19:23:04 -0400 Subject: [PATCH 2/4] document the fully_determined behavior --- docs/src/tutorials/initialization.md | 9 ++++++++- src/systems/diffeqs/abstractodesystem.jl | 4 ++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index 8af6a512e2..4886a47a73 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -194,6 +194,13 @@ may not be analytically satisfiable!**. In our case here, if you sit down with a long enough you will see that `λ = 0` is required for this equation, but since we chose `λ = 1` we end up with a set of equations that are impossible to satisfy. +!!! note + + If you would prefer to have an error instead of a warning in the context of non-fully + determined systems, pass the keyword argument `fully_determined = true` into the + problem constructor. Additionally, any warning about not being fully determined can + be supressed via passing `warn_initialize_determined = false`. + ## Diving Deeper: Constructing the Initialization System To get a better sense of the initialization system and to help debug it, you can construct @@ -271,7 +278,7 @@ sol = solve(iprob) ``` !!! note - + For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems, check out the [NonlinearSolve.jl tutorials!](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/). diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index f4268c506a..11abadad5f 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1556,10 +1556,10 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, nunknown = length(unknowns(isys)) if warn_initialize_determined && neqs > nunknown - @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" end if warn_initialize_determined && neqs < nunknown - @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false." + @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true" end parammap = parammap isa DiffEqBase.NullParameters || isempty(parammap) ? From 7627ebbb34194bece2e57a0d0c4d06bac7385c1c Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 19 Jul 2024 19:23:33 -0400 Subject: [PATCH 3/4] format --- docs/src/tutorials/initialization.md | 4 ++-- test/initializationsystem.jl | 12 ++++++++---- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index 4886a47a73..136f671281 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -195,7 +195,7 @@ long enough you will see that `λ = 0` is required for this equation, but since `λ = 1` we end up with a set of equations that are impossible to satisfy. !!! note - + If you would prefer to have an error instead of a warning in the context of non-fully determined systems, pass the keyword argument `fully_determined = true` into the problem constructor. Additionally, any warning about not being fully determined can @@ -278,7 +278,7 @@ sol = solve(iprob) ``` !!! note - + For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems, check out the [NonlinearSolve.jl tutorials!](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/). diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index 5111381121..1d1ae9bf77 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -17,7 +17,8 @@ sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test maximum(abs.(sol[conditions])) < 1e-14 -@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem(pend, 0.0, [], [g => 1]; +@test_throws ModelingToolkit.ExtraVariablesSystemException ModelingToolkit.InitializationProblem( + pend, 0.0, [], [g => 1]; guesses = [ModelingToolkit.missing_variable_defaults(pend); x => 1; y => 0.2], fully_determined = true) @@ -55,7 +56,8 @@ sol = solve(prob.f.initializeprob) sol = solve(prob, Rodas5P()) @test maximum(abs.(sol[conditions][1])) < 1e-14 -@test_throws ModelingToolkit.ExtraVariablesSystemException ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], +@test_throws ModelingToolkit.ExtraVariablesSystemException ODEProblem( + pend, [x => 1], (0.0, 1.5), [g => 1], guesses = ModelingToolkit.missing_variable_defaults(pend), fully_determined = true) @@ -237,7 +239,8 @@ initsol = solve(initprob, reltol = 1e-12, abstol = 1e-12) @test SciMLBase.successful_retcode(initsol) @test maximum(abs.(initsol[conditions])) < 1e-14 -@test_throws ModelingToolkit.ExtraEquationsSystemException ModelingToolkit.InitializationProblem(sys, 0.0, fully_determined = true) +@test_throws ModelingToolkit.ExtraEquationsSystemException ModelingToolkit.InitializationProblem( + sys, 0.0, fully_determined = true) allinit = unknowns(sys) .=> initsol[unknowns(sys)] prob = ODEProblem(sys, allinit, (0, 0.1)) @@ -249,7 +252,8 @@ sol = solve(prob, Rodas5P(), initializealg = BrownFullBasicInit()) prob = ODEProblem(sys, allinit, (0, 0.1)) prob = ODEProblem(sys, [], (0, 0.1), check = false) -@test_throws ModelingToolkit.ExtraEquationsSystemException ODEProblem(sys, [], (0, 0.1), fully_determined = true) +@test_throws ModelingToolkit.ExtraEquationsSystemException ODEProblem( + sys, [], (0, 0.1), fully_determined = true) sol = solve(prob, Rodas5P()) # If initialized incorrectly, then it would be InitialFailure From 704f08531c06bac2e4a701fa58105c599fc8e486 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Fri, 19 Jul 2024 19:59:19 -0400 Subject: [PATCH 4/4] spell check --- docs/src/tutorials/initialization.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index 136f671281..8852fc7ac0 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -195,11 +195,11 @@ long enough you will see that `λ = 0` is required for this equation, but since `λ = 1` we end up with a set of equations that are impossible to satisfy. !!! note - + If you would prefer to have an error instead of a warning in the context of non-fully determined systems, pass the keyword argument `fully_determined = true` into the problem constructor. Additionally, any warning about not being fully determined can - be supressed via passing `warn_initialize_determined = false`. + be suppressed via passing `warn_initialize_determined = false`. ## Diving Deeper: Constructing the Initialization System @@ -278,7 +278,7 @@ sol = solve(iprob) ``` !!! note - + For more information on solving NonlinearProblems and NonlinearLeastSquaresProblems, check out the [NonlinearSolve.jl tutorials!](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/getting_started/).