diff --git a/Project.toml b/Project.toml index 7b0997823c..b0bd4381b6 100644 --- a/Project.toml +++ b/Project.toml @@ -73,14 +73,16 @@ MTKBifurcationKitExt = "BifurcationKit" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKHomotopyContinuationExt = "HomotopyContinuation" -MTKLabelledArraysExt = "LabelledArrays" MTKInfiniteOptExt = "InfiniteOpt" +MTKLabelledArraysExt = "LabelledArrays" [compat] AbstractTrees = "0.3, 0.4" ArrayInterface = "6, 7" BifurcationKit = "0.4" BlockArrays = "1.1" +BoundaryValueDiffEq = "5.12.0" +BoundaryValueDiffEqAscher = "1.1.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" @@ -139,8 +141,8 @@ SimpleNonlinearSolve = "0.1.0, 1, 2" SparseArrays = "1" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" StaticArrays = "0.10, 0.11, 0.12, 1.0" -StochasticDiffEq = "6.72.1" StochasticDelayDiffEq = "1.8.1" +StochasticDiffEq = "6.72.1" SymbolicIndexingInterface = "0.3.36" SymbolicUtils = "3.7" Symbolics = "6.19" @@ -152,6 +154,8 @@ julia = "1.9" [extras] AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d" +BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" @@ -183,4 +187,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEq", "BoundaryValueDiffEqAscher", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET", "OrdinaryDiffEqNonlinearSolve"] diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 10aba4d8a9..2d53c61401 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -150,6 +150,10 @@ include("systems/imperative_affect.jl") include("systems/callbacks.jl") include("systems/problem_utils.jl") +include("systems/optimization/constraints_system.jl") +include("systems/optimization/optimizationsystem.jl") +include("systems/optimization/modelingtoolkitize.jl") + include("systems/nonlinear/nonlinearsystem.jl") include("systems/nonlinear/homotopy_continuation.jl") include("systems/diffeqs/odesystem.jl") @@ -165,10 +169,6 @@ include("systems/discrete_system/discrete_system.jl") include("systems/jumps/jumpsystem.jl") -include("systems/optimization/constraints_system.jl") -include("systems/optimization/optimizationsystem.jl") -include("systems/optimization/modelingtoolkitize.jl") - include("systems/pde/pdesystem.jl") include("systems/sparsematrixclil.jl") diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 168260ae69..7c0fe0285d 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -983,6 +983,7 @@ for prop in [:eqs :structure :op :constraints + :constraintsystem :controls :loss :bcs diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 20293068a9..47e5c3b88c 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -827,6 +827,12 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `ODEProblem`") end + + if !isnothing(get_constraintsystem(sys)) + error("An ODESystem with constraints cannot be used to construct a regular ODEProblem. + Consider a BVProblem instead.") + end + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, u0map, parammap; t = tspan !== nothing ? tspan[1] : tspan, check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) @@ -849,6 +855,175 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = end get_callback(prob::ODEProblem) = prob.kwargs[:callback] +""" +```julia +SciMLBase.BVProblem{iip}(sys::AbstractODESystem, u0map, tspan, + parammap = DiffEqBase.NullParameters(); + constraints = nothing, guesses = nothing, + version = nothing, tgrad = false, + jac = true, sparse = true, + simplify = false, + kwargs...) where {iip} +``` + +Create a boundary value problem from the [`ODESystem`](@ref). + +`u0map` is used to specify fixed initial values for the states. Every variable +must have either an initial guess supplied using `guesses` or a fixed initial +value specified using `u0map`. + +Boundary value conditions are supplied to ODESystems +in the form of a ConstraintsSystem. These equations +should specify values that state variables should +take at specific points, as in `x(0.5) ~ 1`). More general constraints that +should hold over the entire solution, such as `x(t)^2 + y(t)^2`, should be +specified as one of the equations used to build the `ODESystem`. + +If an ODESystem without `constraints` is specified, it will be treated as an initial value problem. + +```julia + @parameters g t_c = 0.5 + @variables x(..) y(t) [state_priority = 10] λ(t) + eqs = [D(D(x(t))) ~ λ * x(t) + D(D(y)) ~ λ * y - g + x(t)^2 + y^2 ~ 1] + cstr = [x(0.5) ~ 1] + @named cstrs = ConstraintsSystem(cstr, t) + @mtkbuild pend = ODESystem(eqs, t) + + tspan = (0.0, 1.5) + u0map = [x(t) => 0.6, y => 0.8] + parammap = [g => 1] + guesses = [λ => 1] + constraints = [x(0.5) ~ 1] + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; constraints, guesses, check_length = false) +``` + +If the `ODESystem` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting +`BVProblem` must be solved using BVDAE solvers, such as Ascher. +""" +function SciMLBase.BVProblem(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem(sys::AbstractODESystem, + u0map::StaticArray, + args...; + kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, u0map, args...; kwargs...) +end + +function SciMLBase.BVProblem{true}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{false}(sys::AbstractODESystem, args...; kwargs...) + BVProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.BVProblem{iip, specialize}(sys::AbstractODESystem, u0map = [], + tspan = get_tspan(sys), + parammap = DiffEqBase.NullParameters(); + guesses = Dict(), + version = nothing, tgrad = false, + callback = nothing, + check_length = true, + warn_initialize_determined = true, + eval_expression = false, + eval_module = @__MODULE__, + kwargs...) where {iip, specialize} + + if !iscomplete(sys) + error("A completed system is required. Call `complete` or `structural_simplify` on the system before creating an `BVProblem`") + end + !isnothing(callback) && error("BVP solvers do not support callbacks.") + + has_alg_eqs(sys) && error("The BVProblem constructor currently does not support ODESystems with algebraic equations.") # Remove this when the BVDAE solvers get updated, the codegen should work when it does. + + sts = unknowns(sys) + ps = parameters(sys) + constraintsys = get_constraintsystem(sys) + + if !isnothing(constraintsys) + (length(constraints(constraintsys)) + length(u0map) > length(sts)) && + @warn "The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The BVP solvers will default to doing a nonlinear least-squares optimization." + end + + # ODESystems without algebraic equations should use both fixed values + guesses + # for initialization. + _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses)) + f, u0, p = process_SciMLProblem(ODEFunction{iip, specialize}, sys, _u0map, parammap; + t = tspan !== nothing ? tspan[1] : tspan, guesses, + check_length, warn_initialize_determined, eval_expression, eval_module, kwargs...) + + stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) + u0_idxs = has_alg_eqs(sys) ? collect(1:length(sts)) : [stidxmap[k] for (k,v) in u0map] + + bc = generate_function_bc(sys, u0, u0_idxs, tspan, iip) + return BVProblem{iip}(f, bc, u0, tspan, p; kwargs...) +end + +get_callback(prob::BVProblem) = error("BVP solvers do not support callbacks.") + +""" + generate_function_bc(sys::ODESystem, u0, u0_idxs, tspan, iip) + + Given an ODESystem with constraints, generate the boundary condition function to pass to boundary value problem solvers. +""" +function generate_function_bc(sys::ODESystem, u0, u0_idxs, tspan, iip) + iv = get_iv(sys) + sts = get_unknowns(sys) + ps = get_ps(sys) + np = length(ps) + ns = length(sts) + stidxmap = Dict([v => i for (i, v) in enumerate(sts)]) + pidxmap = Dict([v => i for (i, v) in enumerate(ps)]) + + @variables sol(..)[1:ns] p[1:np] + + conssys = get_constraintsystem(sys) + cons = Any[] + if !isnothing(conssys) + cons = [con.lhs - con.rhs for con in constraints(conssys)] + + for st in get_unknowns(conssys) + x = operation(st) + t = only(arguments(st)) + idx = stidxmap[x(iv)] + + cons = map(c -> Symbolics.substitute(c, Dict(x(t) => sol(t)[idx])), cons) + end + + for var in parameters(conssys) + if iscall(var) + x = operation(var) + t = only(arguments(var)) + idx = pidxmap[x] + + cons = map(c -> Symbolics.substitute(c, Dict(x(t) => p[idx])), cons) + else + idx = pidxmap[var] + cons = map(c -> Symbolics.substitute(c, Dict(var => p[idx])), cons) + end + end + end + + init_conds = Any[] + for i in u0_idxs + expr = sol(tspan[1])[i] - u0[i] + push!(init_conds, expr) + end + + exprs = vcat(init_conds, cons) + bcs = Symbolics.build_function(exprs, sol, p, expression = Val{false}) + if iip + return (resid, u, p, t) -> bcs[2](resid, u, p) + else + return (u, p, t) -> bcs[1](u, p) + end +end + """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, diff --git a/src/systems/diffeqs/bvpsystem.jl b/src/systems/diffeqs/bvpsystem.jl new file mode 100644 index 0000000000..ae3029fa14 --- /dev/null +++ b/src/systems/diffeqs/bvpsystem.jl @@ -0,0 +1,217 @@ +""" +$(TYPEDEF) + +A system of ordinary differential equations. + +# Fields +$(FIELDS) + +# Example + +```julia +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +@parameters σ ρ β +@variables x(t) y(t) z(t) + +eqs = [D(x) ~ σ*(y-x), + D(y) ~ x*(ρ-z)-y, + D(z) ~ x*y - β*z] + +@named de = ODESystem(eqs,t,[x,y,z],[σ,ρ,β],tspan=(0, 1000.0)) +``` +""" +struct ODESystem <: AbstractODESystem + """ + A tag for the system. If two systems have the same tag, then they are + structurally identical. + """ + tag::UInt + """The ODEs defining the system.""" + eqs::Vector{Equation} + """Independent variable.""" + iv::BasicSymbolic{Real} + """ + Dependent (unknown) variables. Must not contain the independent variable. + + N.B.: If `torn_matching !== nothing`, this includes all variables. Actual + ODE unknowns are determined by the `SelectedState()` entries in `torn_matching`. + """ + unknowns::Vector + """Parameter variables. Must not contain the independent variable.""" + ps::Vector + """Time span.""" + tspan::Union{NTuple{2, Any}, Nothing} + """Array variables.""" + var_to_name::Any + """Control parameters (some subset of `ps`).""" + ctrls::Vector + """Observed variables.""" + observed::Vector{Equation} + """ + Time-derivative matrix. Note: this field will not be defined until + [`calculate_tgrad`](@ref) is called on the system. + """ + tgrad::RefValue{Vector{Num}} + """ + Jacobian matrix. Note: this field will not be defined until + [`calculate_jacobian`](@ref) is called on the system. + """ + jac::RefValue{Any} + """ + Control Jacobian matrix. Note: this field will not be defined until + [`calculate_control_jacobian`](@ref) is called on the system. + """ + ctrl_jac::RefValue{Any} + """ + Note: this field will not be defined until + [`generate_factorized_W`](@ref) is called on the system. + """ + Wfact::RefValue{Matrix{Num}} + """ + Note: this field will not be defined until + [`generate_factorized_W`](@ref) is called on the system. + """ + Wfact_t::RefValue{Matrix{Num}} + """ + The name of the system. + """ + name::Symbol + """ + A description of the system. + """ + description::String + """ + The internal systems. These are required to have unique names. + """ + systems::Vector{ODESystem} + """ + The default values to use when initial conditions and/or + parameters are not supplied in `ODEProblem`. + """ + defaults::Dict + """ + The guesses to use as the initial conditions for the + initialization system. + """ + guesses::Dict + """ + Tearing result specifying how to solve the system. + """ + torn_matching::Union{Matching, Nothing} + """ + The system for performing the initialization. + """ + initializesystem::Union{Nothing, NonlinearSystem} + """ + Extra equations to be enforced during the initialization sequence. + """ + initialization_eqs::Vector{Equation} + """ + The schedule for the code generation process. + """ + schedule::Any + """ + Type of the system. + """ + connector_type::Any + """ + Inject assignment statements before the evaluation of the RHS function. + """ + preface::Any + """ + A `Vector{SymbolicContinuousCallback}` that model events. + The integrator will use root finding to guarantee that it steps at each zero crossing. + """ + continuous_events::Vector{SymbolicContinuousCallback} + """ + A `Vector{SymbolicDiscreteCallback}` that models events. Symbolic + analog to `SciMLBase.DiscreteCallback` that executes an affect when a given condition is + true at the end of an integration step. + """ + discrete_events::Vector{SymbolicDiscreteCallback} + """ + Topologically sorted parameter dependency equations, where all symbols are parameters and + the LHS is a single parameter. + """ + parameter_dependencies::Vector{Equation} + """ + Metadata for the system, to be used by downstream packages. + """ + metadata::Any + """ + Metadata for MTK GUI. + """ + gui_metadata::Union{Nothing, GUIMetadata} + """ + A boolean indicating if the given `ODESystem` represents a system of DDEs. + """ + is_dde::Bool + """ + A list of points to provide to the solver as tstops. Uses the same syntax as discrete + events. + """ + tstops::Vector{Any} + """ + Cache for intermediate tearing state. + """ + tearing_state::Any + """ + Substitutions generated by tearing. + """ + substitutions::Any + """ + If a model `sys` is complete, then `sys.x` no longer performs namespacing. + """ + complete::Bool + """ + Cached data for fast symbolic indexing. + """ + index_cache::Union{Nothing, IndexCache} + """ + A list of discrete subsystems. + """ + discrete_subsystems::Any + """ + A list of actual unknowns needed to be solved by solvers. + """ + solved_unknowns::Union{Nothing, Vector{Any}} + """ + A vector of vectors of indices for the split parameters. + """ + split_idxs::Union{Nothing, Vector{Vector{Int}}} + """ + The hierarchical parent system before simplification. + """ + parent::Any + + function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, + jac, ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, + torn_matching, initializesystem, initialization_eqs, schedule, + connector_type, preface, cevents, + devents, parameter_dependencies, + metadata = nothing, gui_metadata = nothing, is_dde = false, + tstops = [], tearing_state = nothing, + substitutions = nothing, complete = false, index_cache = nothing, + discrete_subsystems = nothing, solved_unknowns = nothing, + split_idxs = nothing, parent = nothing; checks::Union{Bool, Int} = true) + if checks == true || (checks & CheckComponents) > 0 + check_independent_variables([iv]) + check_variables(dvs, iv) + check_parameters(ps, iv) + check_equations(deqs, iv) + check_equations(equations(cevents), iv) + end + if checks == true || (checks & CheckUnits) > 0 + u = __get_unit_type(dvs, ps, iv) + check_units(u, deqs) + end + new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, + ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, torn_matching, + initializesystem, initialization_eqs, schedule, connector_type, preface, + cevents, devents, parameter_dependencies, metadata, + gui_metadata, is_dde, tstops, tearing_state, substitutions, complete, index_cache, + discrete_subsystems, solved_unknowns, split_idxs, parent) + end +end diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 61f16fd926..4ef9f89c6a 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -49,6 +49,8 @@ struct ODESystem <: AbstractODESystem ctrls::Vector """Observed variables.""" observed::Vector{Equation} + """System of constraints that must be satisfied by the solution to the system.""" + constraintsystem::Union{Nothing, ConstraintsSystem} """ Time-derivative matrix. Note: this field will not be defined until [`calculate_tgrad`](@ref) is called on the system. @@ -186,7 +188,7 @@ struct ODESystem <: AbstractODESystem """ parent::Any - function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, + function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, constraints, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, torn_matching, initializesystem, initialization_eqs, schedule, connector_type, preface, cevents, @@ -207,7 +209,7 @@ struct ODESystem <: AbstractODESystem u = __get_unit_type(dvs, ps, iv) check_units(u, deqs) end - new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, + new(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, constraints, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, torn_matching, initializesystem, initialization_eqs, schedule, connector_type, preface, cevents, devents, parameter_dependencies, metadata, @@ -219,6 +221,7 @@ end function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; controls = Num[], observed = Equation[], + constraints = Equation[], systems = ODESystem[], tspan = nothing, name = nothing, @@ -283,11 +286,26 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; cont_callbacks = SymbolicContinuousCallbacks(continuous_events) disc_callbacks = SymbolicDiscreteCallbacks(discrete_events) + constraintsys = nothing + if !isempty(constraints) + constraintsys = process_constraint_system(constraints, dvs′, ps′, iv, systems) + dvset = Set(dvs′) + pset = Set(ps′) + for st in get_unknowns(constraintsys) + iscall(st) ? + !in(operation(st)(iv), dvset) && push!(dvs′, st) : + !in(st, dvset) && push!(dvs′, st) + end + for p in parameters(constraintsys) + !in(p, pset) && push!(ps′, p) + end + end + if is_dde === nothing is_dde = _check_if_dde(deqs, iv′, systems) end ODESystem(Threads.atomic_add!(SYSTEM_COUNT, UInt(1)), - deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, tgrad, jac, + deqs, iv′, dvs′, ps′, tspan, var_to_name, ctrl′, observed, constraintsys, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, nothing, initializesystem, initialization_eqs, schedule, connector_type, preface, cont_callbacks, @@ -295,7 +313,7 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; metadata, gui_metadata, is_dde, tstops, checks = checks) end -function ODESystem(eqs, iv; kwargs...) +function ODESystem(eqs, iv; constraints = Equation[], kwargs...) eqs = collect(eqs) # NOTE: this assumes that the order of algebraic equations doesn't matter diffvars = OrderedSet() @@ -358,9 +376,10 @@ function ODESystem(eqs, iv; kwargs...) end end algevars = setdiff(allunknowns, diffvars) + # the orders here are very important! return ODESystem(Equation[diffeq; algeeq; compressed_eqs], iv, - collect(Iterators.flatten((diffvars, algevars))), collect(new_ps); kwargs...) + collect(Iterators.flatten((diffvars, algevars))), collect(new_ps); constraints, kwargs...) end # NOTE: equality does not check cached Jacobian @@ -770,3 +789,57 @@ function Base.show(io::IO, mime::MIME"text/plain", sys::ODESystem; hint = true, return nothing end + +# Validate that all the variables in the BVP constraints are well-formed states or parameters. +# - Any callable with multiple arguments will error. +# - Callable/delay variables (e.g. of the form x(0.6) should be unknowns of the system (and have one arg, etc.) +# - Callable/delay parameters should be parameters of the system (and have one arg, etc.) +function validate_constraint_syms(constraintsts, constraintps, sts, ps, iv) + for var in constraintsts + if !iscall(var) + occursin(iv, var) && var ∈ sts || throw(ArgumentError("Time-dependent variable $var is not an unknown of the system.")) + elseif length(arguments(var)) > 1 + throw(ArgumentError("Too many arguments for variable $var.")) + elseif length(arguments(var)) == 1 + arg = first(arguments(var)) + operation(var)(iv) ∈ sts || + throw(ArgumentError("Variable $var is not a variable of the ODESystem. Called variables must be variables of the ODESystem.")) + + isequal(arg, iv) || isparameter(arg) || arg isa Integer || arg isa AbstractFloat || + throw(ArgumentError("Invalid argument specified for variable $var. The argument of the variable should be either $iv, a parameter, or a value specifying the time that the constraint holds.")) + else + var ∈ sts && @warn "Variable $var has no argument. It will be interpreted as $var($iv), and the constraint will apply to the entire interval." + end + end + + for var in constraintps + !iscall(var) && continue + + if length(arguments(var)) > 1 + throw(ArgumentError("Too many arguments for parameter $var in equation $eq.")) + elseif length(arguments(var)) == 1 + arg = first(arguments(var)) + operation(var) ∈ ps || throw(ArgumentError("Parameter $var is not a parameter of the ODESystem. Called parameters must be parameters of the ODESystem.")) + + isequal(arg, iv) || + arg isa Integer || + arg isa AbstractFloat || + throw(ArgumentError("Invalid argument specified for callable parameter $var. The argument of the parameter should be either $iv, a parameter, or a value specifying the time that the constraint holds.")) + end + end +end + +function process_constraint_system(constraints::Vector{Equation}, sts, ps, iv, subsys::Vector{ODESystem}; name = :cons) + isempty(constraints) && return nothing + + constraintsts = OrderedSet() + constraintps = OrderedSet() + + for cons in constraints + syms = collect_vars!(constraintsts, constraintps, cons, iv) + end + validate_constraint_syms(constraintsts, constraintps, Set(sts), Set(ps), iv) + + constraint_subsys = get_constraintsystem.(subsys) + ConstraintsSystem(constraints, collect(constraintsts), collect(constraintps); systems = constraint_subsys, name) +end diff --git a/src/systems/optimization/constraints_system.jl b/src/systems/optimization/constraints_system.jl index 03225fc900..61e190e672 100644 --- a/src/systems/optimization/constraints_system.jl +++ b/src/systems/optimization/constraints_system.jl @@ -89,6 +89,10 @@ struct ConstraintsSystem <: AbstractTimeIndependentSystem tearing_state = nothing, substitutions = nothing, complete = false, index_cache = nothing; checks::Union{Bool, Int} = true) + + ##if checks == true || (checks & CheckComponents) > 0 + ## check_variables(unknowns, constraints) + ##end if checks == true || (checks & CheckUnits) > 0 u = __get_unit_type(unknowns, ps) check_units(u, constraints) @@ -253,3 +257,4 @@ function get_cmap(sys::ConstraintsSystem, exprs = nothing) cmap = map(x -> x ~ getdefault(x), cs) return cmap, cs end + diff --git a/test/bvproblem.jl b/test/bvproblem.jl new file mode 100644 index 0000000000..c6c124d09e --- /dev/null +++ b/test/bvproblem.jl @@ -0,0 +1,303 @@ +### TODO: update when BoundaryValueDiffEqAscher is updated to use the normal boundary condition conventions + +using BoundaryValueDiffEq, OrdinaryDiffEq, BoundaryValueDiffEqAscher +using BenchmarkTools +using ModelingToolkit +using SciMLBase +using ModelingToolkit: t_nounits as t, D_nounits as D + +### Test Collocation solvers on simple problems +solvers = [MIRK4] +daesolvers = [Ascher2, Ascher4, Ascher6] + +let + @parameters α=7.5 β=4.0 γ=8.0 δ=5.0 + @variables x(t)=1.0 y(t)=2.0 + + eqs = [D(x) ~ α * x - β * x * y, + D(y) ~ -γ * y + δ * x * y] + + u0map = [x => 1.0, y => 2.0] + parammap = [α => 7.5, β => 4, γ => 8.0, δ => 5.0] + tspan = (0.0, 10.0) + + @mtkbuild lotkavolterra = ODESystem(eqs, t) + op = ODEProblem(lotkavolterra, u0map, tspan, parammap) + osol = solve(op, Vern9()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) + + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] + end + + # Test out of place + bvp2 = SciMLBase.BVProblem{false, SciMLBase.AutoSpecialize}( + lotkavolterra, u0map, tspan, parammap; eval_expression = true) + + for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [1.0, 2.0] + end +end + +### Testing on pendulum +let + @parameters g=9.81 L=1.0 + @variables θ(t) = π / 2 θ_t(t) + + eqs = [D(θ) ~ θ_t + D(θ_t) ~ -(g / L) * sin(θ)] + + @mtkbuild pend = ODESystem(eqs, t) + + u0map = [θ => π / 2, θ_t => π / 2] + parammap = [:L => 1.0, :g => 9.81] + tspan = (0.0, 6.0) + + op = ODEProblem(pend, u0map, tspan, parammap) + osol = solve(op, Vern9()) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap) + for solver in solvers + sol = solve(bvp, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] + end + + # Test out-of-place + bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) + + for solver in solvers + sol = solve(bvp2, solver(), dt = 0.01) + @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) + @test sol.u[1] == [π / 2, π / 2] + end +end + +################################################################## +### ODESystem with constraint equations, DAEs with constraints ### +################################################################## + +# Test generation of boundary condition function using `generate_function_bc`. Compare solutions to manually written boundary conditions +let + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 + @variables x(..) y(..) + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), + D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] + + tspan = (0., 1.) + @mtkbuild lksys = ODESystem(eqs, t) + + function lotkavolterra!(du, u, p, t) + du[1] = p[1]*u[1] - p[2]*u[1]*u[2] + du[2] = -p[4]*u[2] + p[3]*u[1]*u[2] + end + + function lotkavolterra(u, p, t) + [p[1]*u[1] - p[2]*u[1]*u[2], -p[4]*u[2] + p[3]*u[1]*u[2]] + end + # Compare the built bc function to the actual constructed one. + function bc!(resid, u, p, t) + resid[1] = u[1][1] - 1. + resid[2] = u[1][2] - 2. + nothing + end + function bc(u, p, t) + [u[1][1] - 1., u[1][2] - 2.] + end + + u0 = [1., 2.]; p = [1.5, 1., 1., 3.] + genbc_iip = ModelingToolkit.generate_function_bc(lksys, u0, [1, 2], tspan, true) + genbc_oop = ModelingToolkit.generate_function_bc(lksys, u0, [1, 2], tspan, false) + + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, [1.,2.], tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, [1.,2.], tspan, p) + + sol1 = solve(bvpi1, MIRK4(), dt = 0.05) + sol2 = solve(bvpi2, MIRK4(), dt = 0.05) + @test sol1 ≈ sol2 + + bvpo1 = BVProblem(lotkavolterra, bc, [1,2], tspan, p) + bvpo2 = BVProblem(lotkavolterra, genbc_oop, [1,2], tspan, p) + + sol1 = solve(bvpo1, MIRK4(), dt = 0.05) + sol2 = solve(bvpo2, MIRK4(), dt = 0.05) + @test sol1 ≈ sol2 + + # Test with a constraint. + constr = [y(0.5) ~ 2.] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + + function bc!(resid, u, p, t) + resid[1] = u(0.0)[1] - 1. + resid[2] = u(0.5)[2] - 2. + end + function bc(u, p, t) + [u(0.0)[1] - 1., u(0.5)[2] - 2.] + end + + u0 = [1, 1.] + genbc_iip = ModelingToolkit.generate_function_bc(lksys, u0, [1], tspan, true) + genbc_oop = ModelingToolkit.generate_function_bc(lksys, u0, [1], tspan, false) + + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, u0, tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra!, genbc_iip, u0, tspan, p) + bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.]) + bvpi4 = SciMLBase.BVProblem{true, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.]) + + sol1 = @btime solve($bvpi1, MIRK4(), dt = 0.01) + sol2 = @btime solve($bvpi2, MIRK4(), dt = 0.01) + sol3 = @btime solve($bvpi3, MIRK4(), dt = 0.01) + sol4 = @btime solve($bvpi4, MIRK4(), dt = 0.01) + @test sol1 ≈ sol2 ≈ sol3 ≈ sol4 # don't get true equality here, not sure why + + bvpo1 = BVProblem(lotkavolterra, bc, u0, tspan, p) + bvpo2 = BVProblem(lotkavolterra, genbc_oop, u0, tspan, p) + bvpo3 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, [x(t) => 1.], tspan; guesses = [y(t) => 1.]) + + sol1 = @btime solve($bvpo1, MIRK4(), dt = 0.05) + sol2 = @btime solve($bvpo2, MIRK4(), dt = 0.05) + sol3 = @btime solve($bvpo3, MIRK4(), dt = 0.05) + @test sol1 ≈ sol2 ≈ sol3 +end + +function test_solvers(solvers, prob, u0map, constraints, equations = []; dt = 0.05, atol = 1e-3) + for solver in solvers + println("Solver: $solver") + sol = @btime solve($prob, $solver(), dt = $dt, abstol = $atol) + @test SciMLBase.successful_retcode(sol.retcode) + p = prob.p; t = sol.t; bc = prob.f.bc + ns = length(prob.u0) + if isinplace(prob.f) + resid = zeros(ns) + bc(resid, sol, p, t) + @test isapprox(zeros(ns), resid; atol) + @show resid + else + @test isapprox(zeros(ns), bc(sol, p, t); atol) + @show bc(sol, p, t) + end + + for (k, v) in u0map + @test sol[k][1] == v + end + + # for cons in constraints + # @test sol[cons.rhs - cons.lhs] ≈ 0 + # end + + for eq in equations + @test sol[eq] ≈ 0 + end + end +end + +# Simple ODESystem with BVP constraints. +let + @parameters α=1.5 β=1.0 γ=3.0 δ=1.0 + @variables x(..) y(..) + + eqs = [D(x(t)) ~ α * x(t) - β * x(t) * y(t), + D(y(t)) ~ -γ * y(t) + δ * x(t) * y(t)] + + u0map = [] + tspan = (0.0, 1.0) + guesses = [x(t) => 4.0, y(t) => 2.] + constr = [x(.6) ~ 3.5, x(.3) ~ 7.] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses) + test_solvers(solvers, bvp, u0map, constr; dt = 0.05) + + # Testing that more complicated constr give correct solutions. + constr = [y(.2) + x(.8) ~ 3., y(.3) ~ 2.] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(lksys, u0map, tspan; guesses) + test_solvers(solvers, bvp, u0map, constr; dt = 0.05) + + constr = [α * β - x(.6) ~ 0.0, y(.2) ~ 3.] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(lksys, u0map, tspan; guesses) + test_solvers(solvers, bvp, u0map, constr) +end + +# Cartesian pendulum from the docs. +# DAE IVP solved using BoundaryValueDiffEq solvers. +# let +# @parameters g +# @variables x(t) y(t) [state_priority = 10] λ(t) +# eqs = [D(D(x)) ~ λ * x +# D(D(y)) ~ λ * y - g +# x^2 + y^2 ~ 1] +# @mtkbuild pend = ODESystem(eqs, t) +# +# tspan = (0.0, 1.5) +# u0map = [x => 1, y => 0] +# pmap = [g => 1] +# guess = [λ => 1] +# +# prob = ODEProblem(pend, u0map, tspan, pmap; guesses = guess) +# osol = solve(prob, Rodas5P()) +# +# zeta = [0., 0., 0., 0., 0.] +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses = guess) +# +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.001) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test isapprox([sol[conditions][1]; sol[x][1] - 1; sol[y][1]], zeros(5), atol = 0.001) +# end +# +# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) +# for solver in solvers +# sol = solve(bvp, solver(zeta), dt = 0.01) +# @test isapprox(sol.u[end], osol.u[end]; atol = 0.01) +# conditions = getfield.(equations(pend)[3:end], :rhs) +# @test [sol[conditions][1]; sol[x][1] - 1; sol[y][1]] ≈ 0 +# end +# end + +# Adding a midpoint boundary constraint. +# Solve using BVDAE solvers. +# let +# @parameters g +# @variables x(..) y(t) [state_priority = 10] λ(t) +# eqs = [D(D(x(t))) ~ λ * x(t) +# D(D(y)) ~ λ * y - g +# x(t)^2 + y^2 ~ 1] +# constr = [x(0.5) ~ 1] +# @mtkbuild pend = ODESystem(eqs, t; constr) +# +# tspan = (0.0, 1.5) +# u0map = [x(t) => 0.6, y => 0.8] +# parammap = [g => 1] +# guesses = [λ => 1] +# +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constr) +# +# bvp2 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}(pend, u0map, tspan, parammap) +# test_solvers(daesolvers, bvp2, u0map, constr, get_alg_eqs(pend)) +# +# # More complicated constr. +# u0map = [x(t) => 0.6] +# guesses = [λ => 1, y(t) => 0.8] +# +# constr = [x(0.5) ~ 1, +# x(0.3)^3 + y(0.6)^2 ~ 0.5] +# @mtkbuild pend = ODESystem(eqs, t; constr) +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constr, get_alg_eqs(pend)) +# +# constr = [x(0.4) * g ~ y(0.2), +# y(0.7) ~ 0.3] +# @mtkbuild pend = ODESystem(eqs, t; constr) +# bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false) +# test_solvers(daesolvers, bvp, u0map, constr, get_alg_eqs(pend)) +# end diff --git a/test/odesystem.jl b/test/odesystem.jl index 85d135b338..516fa7d415 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1552,3 +1552,33 @@ end expected_tstops = unique!(sort!(vcat(0.0:0.075:10.0, 0.1, 0.2, 0.65, 0.35, 0.45))) @test all(x -> any(isapprox(x, atol = 1e-6), sol2.t), expected_tstops) end + +@testset "Constraint system construction" begin + @variables x(..) y(..) z(..) + @parameters a b c d e + eqs = [D(x(t)) ~ 3*a*y(t), D(y(t)) ~ x(t) - z(t), D(z(t)) ~ e*x(t)^2] + cons = [x(0.3) ~ c*d, y(0.7) ~ 3] + + # Test variables + parameters infer correctly. + @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + @test issetequal(parameters(sys), [a, c, d, e]) + @test issetequal(unknowns(sys), [x(t), y(t)]) + + @parameters t_c + cons = [x(t_c) ~ 3] + @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + @test_broken issetequal(parameters(sys), [a, e, t_c]) # TODO: unbreak this. + + # Test that bad constraints throw errors. + cons = [x(3, 4) ~ 3] + @test_throws ArgumentError @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + + cons = [x(y(t)) ~ 2] + @test_throws ArgumentError @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + + @variables u(t) v + cons = [x(t) * u ~ 3] + @test_throws ArgumentError @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + cons = [x(t) * v ~ 3] + @test_nowarn @mtkbuild sys = ODESystem(eqs, t; constraints = cons) +end diff --git a/test/runtests.jl b/test/runtests.jl index 4a2ec80e6a..15366dee6a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -81,6 +81,7 @@ end @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") + @safetestset "BVProblem Test" include("bvproblem.jl") @safetestset "print_tree" include("print_tree.jl") @safetestset "Constraints Test" include("constraints.jl") @safetestset "IfLifting Test" include("if_lifting.jl")