diff --git a/Project.toml b/Project.toml index b1f5701d13..d81b3d4ef9 100644 --- a/Project.toml +++ b/Project.toml @@ -44,6 +44,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" +SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" @@ -120,13 +121,15 @@ NonlinearSolve = "3.14, 4" OffsetArrays = "1" OrderedCollections = "1" OrdinaryDiffEq = "6.82.0" -OrdinaryDiffEqCore = "1.7.0" +OrdinaryDiffEqCore = "1.13.0" +OrdinaryDiffEqNonlinearSolve = "1.3.0" PrecompileTools = "1" REPL = "1" RecursiveArrayTools = "3.26" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.64" +SCCNonlinearSolve = "1.0.0" +SciMLBase = "2.66" SciMLStructures = "1.0" Serialization = "1" Setfield = "0.7, 0.8, 1" @@ -160,6 +163,7 @@ OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" +OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" @@ -174,4 +178,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", "REPL", "Random", "ReferenceTests", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials", "StochasticDelayDiffEq", "Pkg", "JET"] +test = ["AmplNLWriter", "BenchmarkTools", "ControlSystemsBase", "DataInterpolations", "DelayDiffEq", "NonlinearSolve", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "OrdinaryDiffEqCore", "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 de8e69c41f..59be358349 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -50,6 +50,7 @@ using Distributed import JuliaFormatter using MLStyle using NonlinearSolve +import SCCNonlinearSolve using Reexport using RecursiveArrayTools import Graphs: SimpleDiGraph, add_edge!, incidence_matrix diff --git a/src/structural_transformation/bipartite_tearing/modia_tearing.jl b/src/structural_transformation/bipartite_tearing/modia_tearing.jl index cef2f5f6d7..5da873afdf 100644 --- a/src/structural_transformation/bipartite_tearing/modia_tearing.jl +++ b/src/structural_transformation/bipartite_tearing/modia_tearing.jl @@ -62,6 +62,15 @@ function tear_graph_block_modia!(var_eq_matching, ict, solvable_graph, eqs, vars return nothing end +function build_var_eq_matching(structure::SystemStructure, ::Type{U} = Unassigned; + varfilter::F2 = v -> true, eqfilter::F3 = eq -> true) where {U, F2, F3} + @unpack graph, solvable_graph = structure + var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U) + matching_len = max(length(var_eq_matching), + maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0)) + return complete(var_eq_matching, matching_len), matching_len +end + function tear_graph_modia(structure::SystemStructure, isder::F = nothing, ::Type{U} = Unassigned; varfilter::F2 = v -> true, @@ -78,10 +87,7 @@ function tear_graph_modia(structure::SystemStructure, isder::F = nothing, # find them here [TODO: It would be good to have an explicit example of this.] @unpack graph, solvable_graph = structure - var_eq_matching = maximal_matching(graph, eqfilter, varfilter, U) - matching_len = max(length(var_eq_matching), - maximum(x -> x isa Int ? x : 0, var_eq_matching, init = 0)) - var_eq_matching = complete(var_eq_matching, matching_len) + var_eq_matching, matching_len = build_var_eq_matching(structure, U; varfilter, eqfilter) full_var_eq_matching = copy(var_eq_matching) var_sccs = find_var_sccs(graph, var_eq_matching) vargraph = DiCMOBiGraph{true}(graph, 0, Matching(matching_len)) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 3a9bbd33e1..e44f250a7f 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -162,11 +162,12 @@ object. """ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys), ps = parameters(sys); wrap_code = nothing, postprocess_fbody = nothing, states = nothing, - expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, kwargs...) + expression = Val{true}, eval_expression = false, eval_module = @__MODULE__, + cachesyms::Tuple = (), kwargs...) if !iscomplete(sys) error("A completed system is required. Call `complete` or `structural_simplify` on the system.") end - p = reorder_parameters(sys, unwrap.(ps)) + p = (reorder_parameters(sys, unwrap.(ps))..., cachesyms...) isscalar = !(exprs isa AbstractArray) if wrap_code === nothing wrap_code = isscalar ? identity : (identity, identity) @@ -187,7 +188,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys postprocess_fbody, states, wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘ - wrap_array_vars(sys, exprs; dvs) .∘ + wrap_array_vars(sys, exprs; dvs, cachesyms) .∘ wrap_parameter_dependencies(sys, isscalar), expression = Val{true} ) @@ -199,7 +200,7 @@ function generate_custom_function(sys::AbstractSystem, exprs, dvs = unknowns(sys postprocess_fbody, states, wrap_code = wrap_code .∘ wrap_mtkparameters(sys, isscalar) .∘ - wrap_array_vars(sys, exprs; dvs) .∘ + wrap_array_vars(sys, exprs; dvs, cachesyms) .∘ wrap_parameter_dependencies(sys, isscalar), expression = Val{true} ) @@ -231,117 +232,47 @@ end function wrap_array_vars( sys::AbstractSystem, exprs; dvs = unknowns(sys), ps = parameters(sys), - inputs = nothing, history = false) + inputs = nothing, history = false, cachesyms::Tuple = ()) isscalar = !(exprs isa AbstractArray) - array_vars = Dict{Any, AbstractArray{Int}}() - if dvs !== nothing - for (j, x) in enumerate(dvs) - if iscall(x) && operation(x) == getindex - arg = arguments(x)[1] - inds = get!(() -> Int[], array_vars, arg) - push!(inds, j) - end - end - for (k, inds) in array_vars - if inds == (inds′ = inds[1]:inds[end]) - array_vars[k] = inds′ - end - end + var_to_arridxs = Dict() - uind = 1 - else + if dvs === nothing uind = 0 - end - # values are (indexes, index of buffer, size of parameter) - array_parameters = Dict{Any, Tuple{AbstractArray{Int}, Int, Tuple{Vararg{Int}}}}() - # If for some reason different elements of an array parameter are in different buffers - other_array_parameters = Dict{Any, Any}() - - hasinputs = inputs !== nothing - input_vars = Dict{Any, AbstractArray{Int}}() - if hasinputs - for (j, x) in enumerate(inputs) - if iscall(x) && operation(x) == getindex - arg = arguments(x)[1] - inds = get!(() -> Int[], input_vars, arg) - push!(inds, j) - end - end - for (k, inds) in input_vars - if inds == (inds′ = inds[1]:inds[end]) - input_vars[k] = inds′ - end - end - end - if has_index_cache(sys) - ic = get_index_cache(sys) else - ic = nothing - end - if ps isa Tuple && eltype(ps) <: AbstractArray - ps = Iterators.flatten(ps) - end - for p in ps - p = unwrap(p) - if iscall(p) && operation(p) == getindex - p = arguments(p)[1] - end - symtype(p) <: AbstractArray && Symbolics.shape(p) != Symbolics.Unknown() || continue - scal = collect(p) - # all scalarized variables are in `ps` - any(isequal(p), ps) || all(x -> any(isequal(x), ps), scal) || continue - (haskey(array_parameters, p) || haskey(other_array_parameters, p)) && continue - - idx = parameter_index(sys, p) - idx isa Int && continue - if idx isa ParameterIndex - if idx.portion != SciMLStructures.Tunable() - continue - end - array_parameters[p] = (vec(idx.idx), 1, size(idx.idx)) + uind = 1 + for (i, x) in enumerate(dvs) + iscall(x) && operation(x) == getindex || continue + arg = arguments(x)[1] + inds = get!(() -> [], var_to_arridxs, arg) + push!(inds, (uind, i)) + end + end + p_start = uind + 1 + history + rps = (reorder_parameters(sys, ps)..., cachesyms...) + if inputs !== nothing + rps = (inputs, rps...) + end + for sym in reduce(vcat, rps; init = []) + iscall(sym) && operation(sym) == getindex || continue + arg = arguments(sym)[1] + + bufferidx = findfirst(buf -> any(isequal(sym), buf), rps) + idxinbuffer = findfirst(isequal(sym), rps[bufferidx]) + inds = get!(() -> [], var_to_arridxs, arg) + push!(inds, (p_start + bufferidx - 1, idxinbuffer)) + end + + viewsyms = Dict() + splitsyms = Dict() + for (arrsym, idxs) in var_to_arridxs + length(idxs) == length(arrsym) || continue + # allequal(first, idxs) is a 1.11 feature + if allequal(Iterators.map(first, idxs)) + viewsyms[arrsym] = (first(first(idxs)), reshape(last.(idxs), size(arrsym))) else - # idx === nothing - idxs = map(Base.Fix1(parameter_index, sys), scal) - if first(idxs) isa ParameterIndex - buffer_idxs = map(Base.Fix1(iterated_buffer_index, ic), idxs) - if allequal(buffer_idxs) - buffer_idx = first(buffer_idxs) - if first(idxs).portion == SciMLStructures.Tunable() - idxs = map(x -> x.idx, idxs) - else - idxs = map(x -> x.idx[end], idxs) - end - else - other_array_parameters[p] = scal - continue - end - else - buffer_idx = 1 - end - - sz = size(idxs) - if vec(idxs) == idxs[begin]:idxs[end] - idxs = idxs[begin]:idxs[end] - elseif vec(idxs) == idxs[begin]:-1:idxs[end] - idxs = idxs[begin]:-1:idxs[end] - end - idxs = vec(idxs) - array_parameters[p] = (idxs, buffer_idx, sz) + splitsyms[arrsym] = reshape(idxs, size(arrsym)) end end - - inputind = if history - uind + 2 - else - uind + 1 - end - params_offset = if history && hasinputs - uind + 2 - elseif history || hasinputs - uind + 1 - else - uind - end if isscalar function (expr) Func( @@ -349,15 +280,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx].name), $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← + MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms] ), expr.body, false @@ -371,15 +298,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind].name), $v)) for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx].name), $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← + MakeArray([expr.args[bufi].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1]]) for (sym, idxs) in splitsyms] ), expr.body, false @@ -392,17 +315,11 @@ function wrap_array_vars( [], Let( vcat( - [k ← :(view($(expr.args[uind + 1].name), $v)) - for (k, v) in array_vars], - [k ← :(view($(expr.args[inputind + 1].name), $v)) - for (k, v) in input_vars], - [k ← :(reshape( - view($(expr.args[params_offset + buffer_idx + 1].name), - $idxs), - $sz)) - for (k, (idxs, buffer_idx, sz)) in array_parameters], - [k ← Code.MakeArray(v, symtype(k)) - for (k, v) in other_array_parameters] + [sym ← :(view($(expr.args[i + 1].name), $idxs)) + for (sym, (i, idxs)) in viewsyms], + [sym ← MakeArray( + [expr.args[bufi + 1].elems[vali] for (bufi, vali) in idxs], + expr.args[idxs[1][1] + 1]) for (sym, idxs) in splitsyms] ), expr.body, false diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index af96c4fcfe..f4e29346ff 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1301,6 +1301,7 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, initialization_eqs = [], fully_determined = nothing, check_units = true, + use_scc = true, 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`") @@ -1310,16 +1311,18 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, elseif isempty(u0map) && get_initializesystem(sys) === nothing isys = structural_simplify( generate_initializesystem( - sys; initialization_eqs, check_units, pmap = parammap, guesses); fully_determined) + sys; initialization_eqs, check_units, pmap = parammap, + guesses, extra_metadata = (; use_scc)); fully_determined) else isys = structural_simplify( generate_initializesystem( - sys; u0map, initialization_eqs, check_units, pmap = parammap, guesses); fully_determined) + sys; u0map, initialization_eqs, check_units, + pmap = parammap, guesses, extra_metadata = (; use_scc)); fully_determined) end ts = get_tearing_state(isys) - if warn_initialize_determined && - (unassigned_vars = StructuralTransformations.singular_check(ts); !isempty(unassigned_vars)) + unassigned_vars = StructuralTransformations.singular_check(ts) + if warn_initialize_determined && !isempty(unassigned_vars) errmsg = """ The initialization system is structurally singular. Guess values may \ significantly affect the initial values of the ODE. The problematic variables \ @@ -1341,11 +1344,17 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, neqs = length(equations(isys)) nunknown = length(unknowns(isys)) + if use_scc + scc_message = "`SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. " + else + scc_message = "" + end + 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. To make this warning into an error, pass fully_determined = true" + @warn "Initialization system is overdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)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. To make this warning into an error, pass fully_determined = true" + @warn "Initialization system is underdetermined. $neqs equations for $nunknown unknowns. Initialization will default to using least squares. $(scc_message)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) ? @@ -1381,9 +1390,20 @@ function InitializationProblem{iip, specialize}(sys::AbstractODESystem, end for (k, v) in u0map) end - if neqs == nunknown - NonlinearProblem(isys, u0map, parammap; kwargs...) + + TProb = if neqs == nunknown && isempty(unassigned_vars) + if use_scc && neqs > 0 + if is_split(isys) + SCCNonlinearProblem + else + @warn "`SCCNonlinearProblem` can only be used with `split = true` systems. Simplify your `ODESystem` with `split = true` or pass `use_scc = false` to disable this warning" + NonlinearProblem + end + else + NonlinearProblem + end else - NonlinearLeastSquaresProblem(isys, u0map, parammap; kwargs...) + NonlinearLeastSquaresProblem end + TProb(isys, u0map, parammap; kwargs...) end diff --git a/src/systems/index_cache.jl b/src/systems/index_cache.jl index ab0dd08764..aca6a0547d 100644 --- a/src/systems/index_cache.jl +++ b/src/systems/index_cache.jl @@ -594,3 +594,28 @@ function reorder_dimension_by_tunables( reorder_dimension_by_tunables!(buffer, sys, arr, syms; dim) return buffer end + +function subset_unknowns_observed( + ic::IndexCache, sys::AbstractSystem, newunknowns, newobsvars) + unknown_idx = copy(ic.unknown_idx) + empty!(unknown_idx) + for (i, sym) in enumerate(newunknowns) + ttsym = default_toterm(sym) + rsym = renamespace(sys, sym) + rttsym = renamespace(sys, ttsym) + unknown_idx[sym] = unknown_idx[ttsym] = unknown_idx[rsym] = unknown_idx[rttsym] = i + end + observed_syms_to_timeseries = copy(ic.observed_syms_to_timeseries) + empty!(observed_syms_to_timeseries) + for sym in newobsvars + ttsym = default_toterm(sym) + rsym = renamespace(sys, sym) + rttsym = renamespace(sys, ttsym) + for s in (sym, ttsym, rsym, rttsym) + observed_syms_to_timeseries[s] = ic.observed_syms_to_timeseries[sym] + end + end + ic = @set ic.unknown_idx = unknown_idx + @set! ic.observed_syms_to_timeseries = observed_syms_to_timeseries + return ic +end diff --git a/src/systems/nonlinear/initializesystem.jl b/src/systems/nonlinear/initializesystem.jl index 726d171bd0..2344727920 100644 --- a/src/systems/nonlinear/initializesystem.jl +++ b/src/systems/nonlinear/initializesystem.jl @@ -11,7 +11,7 @@ function generate_initializesystem(sys::ODESystem; default_dd_guess = 0.0, algebraic_only = false, check_units = true, check_defguess = false, - name = nameof(sys), kwargs...) + name = nameof(sys), extra_metadata = (;), kwargs...) trueobs, eqs = unhack_observed(observed(sys), equations(sys)) vars = unique([unknowns(sys); getfield.(trueobs, :lhs)]) vars_set = Set(vars) # for efficient in-lookup @@ -179,7 +179,8 @@ function generate_initializesystem(sys::ODESystem; for k in keys(defs) defs[k] = substitute(defs[k], paramsubs) end - meta = InitializationSystemMetadata(anydict(u0map), anydict(pmap), additional_guesses) + meta = InitializationSystemMetadata( + anydict(u0map), anydict(pmap), additional_guesses, extra_metadata) return NonlinearSystem(eqs_ics, vars, pars; @@ -195,6 +196,7 @@ struct InitializationSystemMetadata u0map::Dict{Any, Any} pmap::Dict{Any, Any} additional_guesses::Dict{Any, Any} + extra_metadata::NamedTuple end function is_parameter_solvable(p, pmap, defs, guesses) @@ -281,6 +283,9 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p, symbols_to_symbolics!(sys, pmap) guesses = Dict() defs = defaults(sys) + cmap, cs = get_cmap(sys) + use_scc = true + if SciMLBase.has_initializeprob(odefn) oldsys = odefn.initializeprob.f.sys meta = get_metadata(oldsys) @@ -288,6 +293,7 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p, u0map = merge(meta.u0map, u0map) pmap = merge(meta.pmap, pmap) merge!(guesses, meta.additional_guesses) + use_scc = get(meta.extra_metadata, :use_scc, true) end else # there is no initializeprob, so the original problem construction @@ -324,8 +330,11 @@ function SciMLBase.remake_initialization_data(sys::ODESystem, odefn, u0, t0, p, end filter_missing_values!(u0map) filter_missing_values!(pmap) - f, _ = process_SciMLProblem(EmptySciMLFunction, sys, u0map, pmap; guesses, t = t0) - kws = f.kwargs + + op, missing_unknowns, missing_pars = build_operating_point( + u0map, pmap, defs, cmap, dvs, ps) + kws = maybe_build_initialization_problem( + sys, op, u0map, pmap, t0, defs, guesses, missing_unknowns; use_scc) initprob = get(kws, :initializeprob, nothing) if initprob === nothing return nothing diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index 1289388197..7826e06f76 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -535,6 +535,186 @@ function DiffEqBase.NonlinearLeastSquaresProblem{iip}(sys::NonlinearSystem, u0ma NonlinearLeastSquaresProblem{iip}(f, u0, p; filter_kwargs(kwargs)...) end +struct CacheWriter{F} + fn::F +end + +function (cw::CacheWriter)(p, sols) + cw.fn(p.caches[1], sols, p...) +end + +function CacheWriter(sys::AbstractSystem, exprs, solsyms, obseqs::Vector{Equation}; + eval_expression = false, eval_module = @__MODULE__) + ps = parameters(sys) + rps = reorder_parameters(sys, ps) + obs_assigns = [eq.lhs ← eq.rhs for eq in obseqs] + cmap, cs = get_cmap(sys) + cmap_assigns = [eq.lhs ← eq.rhs for eq in cmap] + fn = Func( + [:out, DestructuredArgs(DestructuredArgs.(solsyms)), + DestructuredArgs.(rps)...], + [], + SetArray(true, :out, exprs) + ) |> wrap_assignments(false, obs_assigns)[2] |> + wrap_parameter_dependencies(sys, false)[2] |> + wrap_array_vars(sys, exprs; dvs = nothing, inputs = [])[2] |> + wrap_assignments(false, cmap_assigns)[2] |> toexpr + return CacheWriter(eval_or_rgf(fn; eval_expression, eval_module)) +end + +struct SCCNonlinearFunction{iip} end + +function SCCNonlinearFunction{iip}( + sys::NonlinearSystem, _eqs, _dvs, _obs, cachesyms; eval_expression = false, + eval_module = @__MODULE__, kwargs...) where {iip} + ps = parameters(sys) + rps = reorder_parameters(sys, ps) + + obs_assignments = [eq.lhs ← eq.rhs for eq in _obs] + + cmap, cs = get_cmap(sys) + cmap_assignments = [eq.lhs ← eq.rhs for eq in cmap] + rhss = [eq.rhs - eq.lhs for eq in _eqs] + wrap_code = wrap_assignments(false, cmap_assignments) .∘ + (wrap_array_vars(sys, rhss; dvs = _dvs, cachesyms)) .∘ + wrap_parameter_dependencies(sys, false) .∘ + wrap_assignments(false, obs_assignments) + f_gen = build_function( + rhss, _dvs, rps..., cachesyms...; wrap_code, expression = Val{true}) + f_oop, f_iip = eval_or_rgf.(f_gen; eval_expression, eval_module) + + f(u, p) = f_oop(u, p) + f(u, p::MTKParameters) = f_oop(u, p...) + f(resid, u, p) = f_iip(resid, u, p) + f(resid, u, p::MTKParameters) = f_iip(resid, u, p...) + + subsys = NonlinearSystem(_eqs, _dvs, ps; observed = _obs, + parameter_dependencies = parameter_dependencies(sys), name = nameof(sys)) + if get_index_cache(sys) !== nothing + @set! subsys.index_cache = subset_unknowns_observed( + get_index_cache(sys), sys, _dvs, getproperty.(_obs, (:lhs,))) + @set! subsys.complete = true + end + + return NonlinearFunction{iip}(f; sys = subsys) +end + +function SciMLBase.SCCNonlinearProblem(sys::NonlinearSystem, args...; kwargs...) + SCCNonlinearProblem{true}(sys, args...; kwargs...) +end + +function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map, + parammap = SciMLBase.NullParameters(); eval_expression = false, eval_module = @__MODULE__, kwargs...) where {iip} + if !iscomplete(sys) || get_tearing_state(sys) === nothing + error("A simplified `NonlinearSystem` is required. Call `structural_simplify` on the system before creating an `SCCNonlinearProblem`.") + end + + if !is_split(sys) + error("The system has been simplified with `split = false`. `SCCNonlinearProblem` is not compatible with this system. Pass `split = true` to `structural_simplify` to use `SCCNonlinearProblem`.") + end + + ts = get_tearing_state(sys) + var_eq_matching, var_sccs = StructuralTransformations.algebraic_variables_scc(ts) + + if length(var_sccs) == 1 + return NonlinearProblem{iip}( + sys, u0map, parammap; eval_expression, eval_module, kwargs...) + end + + condensed_graph = MatchedCondensationGraph( + DiCMOBiGraph{true}(complete(ts.structure.graph), + complete(var_eq_matching)), + var_sccs) + toporder = topological_sort_by_dfs(condensed_graph) + var_sccs = var_sccs[toporder] + eq_sccs = map(Base.Fix1(getindex, var_eq_matching), var_sccs) + + dvs = unknowns(sys) + ps = parameters(sys) + eqs = equations(sys) + obs = observed(sys) + + _, u0, p = process_SciMLProblem( + EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module, kwargs...) + + explicitfuns = [] + nlfuns = [] + prevobsidxs = Int[] + cachesize = 0 + for (i, (escc, vscc)) in enumerate(zip(eq_sccs, var_sccs)) + # subset unknowns and equations + _dvs = dvs[vscc] + _eqs = eqs[escc] + # get observed equations required by this SCC + obsidxs = observed_equations_used_by(sys, _eqs) + # the ones used by previous SCCs can be precomputed into the cache + setdiff!(obsidxs, prevobsidxs) + _obs = obs[obsidxs] + + # get all subexpressions in the RHS which we can precompute in the cache + banned_vars = Set{Any}(vcat(_dvs, getproperty.(_obs, (:lhs,)))) + for var in banned_vars + iscall(var) || continue + operation(var) === getindex || continue + push!(banned_vars, arguments(var)[1]) + end + state = Dict() + for i in eachindex(_obs) + _obs[i] = _obs[i].lhs ~ subexpressions_not_involving_vars!( + _obs[i].rhs, banned_vars, state) + end + for i in eachindex(_eqs) + _eqs[i] = _eqs[i].lhs ~ subexpressions_not_involving_vars!( + _eqs[i].rhs, banned_vars, state) + end + + # cached variables and their corresponding expressions + cachevars = Any[obs[i].lhs for i in prevobsidxs] + cacheexprs = Any[obs[i].lhs for i in prevobsidxs] + for (k, v) in state + push!(cachevars, unwrap(v)) + push!(cacheexprs, unwrap(k)) + end + cachesize = max(cachesize, length(cachevars)) + + if isempty(cachevars) + push!(explicitfuns, Returns(nothing)) + else + solsyms = getindex.((dvs,), view(var_sccs, 1:(i - 1))) + push!(explicitfuns, + CacheWriter(sys, cacheexprs, solsyms, obs[prevobsidxs]; + eval_expression, eval_module)) + end + f = SCCNonlinearFunction{iip}( + sys, _eqs, _dvs, _obs, (cachevars,); eval_expression, eval_module, kwargs...) + push!(nlfuns, f) + append!(cachevars, _dvs) + append!(cacheexprs, _dvs) + for i in obsidxs + push!(cachevars, obs[i].lhs) + push!(cacheexprs, obs[i].rhs) + end + append!(prevobsidxs, obsidxs) + end + + if cachesize != 0 + p = rebuild_with_caches(p, BufferTemplate(eltype(u0), cachesize)) + end + + subprobs = [] + for (f, vscc) in zip(nlfuns, var_sccs) + prob = NonlinearProblem(f, u0[vscc], p) + push!(subprobs, prob) + end + + new_dvs = dvs[reduce(vcat, var_sccs)] + new_eqs = eqs[reduce(vcat, eq_sccs)] + @set! sys.unknowns = new_dvs + @set! sys.eqs = new_eqs + sys = complete(sys) + return SCCNonlinearProblem(subprobs, explicitfuns, p, true; sys) +end + """ $(TYPEDSIGNATURES) diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl index 72af9f594d..7d64054acc 100644 --- a/src/systems/parameter_buffer.jl +++ b/src/systems/parameter_buffer.jl @@ -3,11 +3,12 @@ symconvert(::Type{T}, x) where {T} = convert(T, x) symconvert(::Type{Real}, x::Integer) = convert(Float64, x) symconvert(::Type{V}, x) where {V <: AbstractArray} = convert(V, symconvert.(eltype(V), x)) -struct MTKParameters{T, D, C, N} +struct MTKParameters{T, D, C, N, H} tunable::T discrete::D constant::C nonnumeric::N + caches::H end """ @@ -181,11 +182,18 @@ function MTKParameters( mtkps = MTKParameters{ typeof(tunable_buffer), typeof(disc_buffer), typeof(const_buffer), - typeof(nonnumeric_buffer)}(tunable_buffer, disc_buffer, const_buffer, - nonnumeric_buffer) + typeof(nonnumeric_buffer), typeof(())}(tunable_buffer, + disc_buffer, const_buffer, nonnumeric_buffer, ()) return mtkps end +function rebuild_with_caches(p::MTKParameters, cache_templates::BufferTemplate...) + buffers = map(cache_templates) do template + Vector{template.type}(undef, template.length) + end + @set p.caches = buffers +end + function narrow_buffer_type(buffer::AbstractArray) type = Union{} for x in buffer @@ -297,7 +305,8 @@ end for (Portion, field, recurse) in [(SciMLStructures.Discrete, :discrete, 1) (SciMLStructures.Constants, :constant, 1) - (Nonnumeric, :nonnumeric, 1)] + (Nonnumeric, :nonnumeric, 1) + (SciMLStructures.Caches, :caches, 1)] @eval function SciMLStructures.canonicalize(::$Portion, p::MTKParameters) as_vector = buffer_to_arraypartition(p.$field) repack = let p = p @@ -324,11 +333,13 @@ function Base.copy(p::MTKParameters) discrete = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.discrete) constant = Tuple(eltype(buf) <: Real ? copy(buf) : copy.(buf) for buf in p.constant) nonnumeric = copy.(p.nonnumeric) + caches = copy.(p.caches) return MTKParameters( tunable, discrete, constant, - nonnumeric + nonnumeric, + caches ) end @@ -640,7 +651,7 @@ end # getindex indexes the vectors, setindex! linearly indexes values # it's inconsistent, but we need it to be this way @generated function Base.getindex( - ps::MTKParameters{T, D, C, N}, idx::Int) where {T, D, C, N} + ps::MTKParameters{T, D, C, N, H}, idx::Int) where {T, D, C, N, H} paths = [] if !(T <: SizedVector{0, Float64}) push!(paths, :(ps.tunable)) @@ -654,6 +665,9 @@ end for i in 1:fieldcount(N) push!(paths, :(ps.nonnumeric[$i])) end + for i in 1:fieldcount(H) + push!(paths, :(ps.caches[$i])) + end expr = Expr(:if, :(idx == 1), :(return $(paths[1]))) curexpr = expr for i in 2:length(paths) @@ -663,12 +677,12 @@ end return Expr(:block, expr, :(throw(BoundsError(ps, idx)))) end -@generated function Base.length(ps::MTKParameters{T, D, C, N}) where {T, D, C, N} +@generated function Base.length(ps::MTKParameters{T, D, C, N, H}) where {T, D, C, N, H} len = 0 if !(T <: SizedVector{0, Float64}) len += 1 end - len += fieldcount(D) + fieldcount(C) + fieldcount(N) + len += fieldcount(D) + fieldcount(C) + fieldcount(N) + fieldcount(H) return len end @@ -691,7 +705,10 @@ end function Base.:(==)(a::MTKParameters, b::MTKParameters) return a.tunable == b.tunable && a.discrete == b.discrete && - a.constant == b.constant && a.nonnumeric == b.nonnumeric + a.constant == b.constant && a.nonnumeric == b.nonnumeric && + all(Iterators.map(a.caches, b.caches) do acache, bcache + eltype(acache) == eltype(bcache) && length(acache) == length(bcache) + end) end # to support linearize/linearization_function diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index b837eef98e..6b75fb2c6d 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -489,6 +489,89 @@ function EmptySciMLFunction(args...; kwargs...) return EmptySciMLFunction{typeof(args), typeof(kwargs)}(args, kwargs) end +""" + $(TYPEDSIGNATURES) + +Construct the operating point of the system from the user-provided `u0map` and `pmap`, system +defaults `defs`, constant equations `cmap` (from `get_cmap(sys)`), unknowns `dvs` and +parameters `ps`. Return the operating point as a dictionary, the list of unknowns for which +no values can be determined, and the list of parameters for which no values can be determined. +""" +function build_operating_point( + u0map::AbstractDict, pmap::AbstractDict, defs::AbstractDict, cmap, dvs, ps) + op = add_toterms(u0map) + missing_unknowns = add_fallbacks!(op, dvs, defs) + for (k, v) in defs + haskey(op, k) && continue + op[k] = v + end + merge!(op, pmap) + missing_pars = add_fallbacks!(op, ps, defs) + for eq in cmap + op[eq.lhs] = eq.rhs + end + return op, missing_unknowns, missing_pars +end + +""" + $(TYPEDSIGNATURES) + +Build and return the initialization problem and associated data as a `NamedTuple` to be passed +to the `SciMLFunction` constructor. Requires the system `sys`, operating point `op`, +user-provided `u0map` and `pmap`, initial time `t`, system defaults `defs`, user-provided +`guesses`, and list of unknowns which don't have a value in `op`. The keyword `implicit_dae` +denotes whether the `SciMLProblem` being constructed is in implicit DAE form (`DAEProblem`). +All other keyword arguments are forwarded to `InitializationProblem`. +""" +function maybe_build_initialization_problem( + sys::AbstractSystem, op::AbstractDict, u0map, pmap, t, defs, + guesses, missing_unknowns; implicit_dae = false, kwargs...) + guesses = merge(ModelingToolkit.guesses(sys), todict(guesses)) + has_observed_u0s = any( + k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k), + keys(op)) + solvablepars = [p + for p in parameters(sys) + if is_parameter_solvable(p, pmap, defs, guesses)] + has_dependent_unknowns = any(unknowns(sys)) do sym + val = get(op, sym, nothing) + val === nothing && return false + return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val) + end + if (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) || + !isempty(solvablepars) || has_dependent_unknowns) && + get_tearing_state(sys) !== nothing) || + !isempty(initialization_equations(sys))) && t !== nothing + initializeprob = ModelingToolkit.InitializationProblem( + sys, t, u0map, pmap; guesses, kwargs...) + initializeprobmap = getu(initializeprob, unknowns(sys)) + + punknowns = [p + for p in all_variable_symbols(initializeprob) + if is_parameter(sys, p)] + getpunknowns = getu(initializeprob, punknowns) + setpunknowns = setp(sys, punknowns) + initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) + + reqd_syms = parameter_symbols(initializeprob) + update_initializeprob! = UpdateInitializeprob( + getu(sys, reqd_syms), setu(initializeprob, reqd_syms)) + for p in punknowns + p = unwrap(p) + stype = symtype(p) + op[p] = get_temporary_value(p) + end + + for v in missing_unknowns + op[v] = zero_var(v) + end + empty!(missing_unknowns) + return (; + initializeprob, initializeprobmap, initializeprobpmap, update_initializeprob!) + end + return (;) +end + """ $(TYPEDSIGNATURES) @@ -541,6 +624,8 @@ Keyword arguments: Only applicable if `warn_cyclic_dependency == true`. - `substitution_limit`: The number times to substitute initial conditions into each other to attempt to arrive at a numeric value. +- `use_scc`: Whether to use `SCCNonlinearProblem` for initialization if the system is fully + determined. All other keyword arguments are passed as-is to `constructor`. """ @@ -554,7 +639,7 @@ function process_SciMLProblem( symbolic_u0 = false, warn_cyclic_dependency = false, circular_dependency_max_cycle_length = length(all_symbols(sys)), circular_dependency_max_cycles = 10, - substitution_limit = 100, kwargs...) + substitution_limit = 100, use_scc = true, kwargs...) dvs = unknowns(sys) ps = parameters(sys) iv = has_iv(sys) ? get_iv(sys) : nothing @@ -574,67 +659,18 @@ function process_SciMLProblem( cmap, cs = get_cmap(sys) kwargs = NamedTuple(kwargs) - op = add_toterms(u0map) - missing_unknowns = add_fallbacks!(op, dvs, defs) - for (k, v) in defs - haskey(op, k) && continue - op[k] = v - end - merge!(op, pmap) - missing_pars = add_fallbacks!(op, ps, defs) - for eq in cmap - op[eq.lhs] = eq.rhs - end - if sys isa ODESystem - guesses = merge(ModelingToolkit.guesses(sys), todict(guesses)) - has_observed_u0s = any( - k -> has_observed_with_lhs(sys, k) || has_parameter_dependency_with_lhs(sys, k), - keys(op)) - solvablepars = [p - for p in parameters(sys) - if is_parameter_solvable(p, pmap, defs, guesses)] - has_dependent_unknowns = any(unknowns(sys)) do sym - val = get(op, sym, nothing) - val === nothing && return false - return symbolic_type(val) != NotSymbolic() || is_array_of_symbolics(val) - end - if build_initializeprob && - (((implicit_dae || has_observed_u0s || !isempty(missing_unknowns) || - !isempty(solvablepars) || has_dependent_unknowns) && - get_tearing_state(sys) !== nothing) || - !isempty(initialization_equations(sys))) && t !== nothing - initializeprob = ModelingToolkit.InitializationProblem( - sys, t, u0map, pmap; guesses, warn_initialize_determined, - initialization_eqs, eval_expression, eval_module, fully_determined, - warn_cyclic_dependency, check_units = check_initialization_units, - circular_dependency_max_cycle_length, circular_dependency_max_cycles) - initializeprobmap = getu(initializeprob, unknowns(sys)) - - punknowns = [p - for p in all_variable_symbols(initializeprob) - if is_parameter(sys, p)] - getpunknowns = getu(initializeprob, punknowns) - setpunknowns = setp(sys, punknowns) - initializeprobpmap = GetUpdatedMTKParameters(getpunknowns, setpunknowns) - - reqd_syms = parameter_symbols(initializeprob) - update_initializeprob! = UpdateInitializeprob( - getu(sys, reqd_syms), setu(initializeprob, reqd_syms)) - for p in punknowns - p = unwrap(p) - stype = symtype(p) - op[p] = get_temporary_value(p) - delete!(missing_pars, p) - end + op, missing_unknowns, missing_pars = build_operating_point( + u0map, pmap, defs, cmap, dvs, ps) - for v in missing_unknowns - op[v] = zero_var(v) - end - empty!(missing_unknowns) - kwargs = merge(kwargs, - (; initializeprob, initializeprobmap, - initializeprobpmap, update_initializeprob!)) - end + if sys isa ODESystem && build_initializeprob + kws = maybe_build_initialization_problem( + sys, op, u0map, pmap, t, defs, guesses, missing_unknowns; + implicit_dae, warn_initialize_determined, initialization_eqs, + eval_expression, eval_module, fully_determined, + warn_cyclic_dependency, check_units = check_initialization_units, + circular_dependency_max_cycle_length, circular_dependency_max_cycles, use_scc) + + kwargs = merge(kwargs, kws) end if t !== nothing && !(constructor <: Union{DDEFunction, SDDEFunction}) diff --git a/src/utils.jl b/src/utils.jl index 416efd8f2c..dd7113cbe6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1001,9 +1001,143 @@ end diff2term_with_unit(x, t) = _with_unit(diff2term, x, t) lower_varname_with_unit(var, iv, order) = _with_unit(lower_varname, var, iv, iv, order) +""" + $(TYPEDSIGNATURES) + +Check if `sym` represents a symbolic floating point number or array of such numbers. +""" function is_variable_floatingpoint(sym) sym = unwrap(sym) T = symtype(sym) return T == Real || T <: AbstractFloat || T <: AbstractArray{Real} || T <: AbstractArray{<:AbstractFloat} end + +""" + $(TYPEDSIGNATURES) + +Return the `DiCMOBiGraph` denoting the dependencies between observed equations `eqs`. +""" +function observed_dependency_graph(eqs::Vector{Equation}) + for eq in eqs + if symbolic_type(eq.lhs) == NotSymbolic() + error("All equations must be observed equations of the form `var ~ expr`. Got $eq") + end + end + graph, assigns = observed2graph(eqs, getproperty.(eqs, (:lhs,))) + matching = complete(Matching(Vector{Union{Unassigned, Int}}(assigns))) + return DiCMOBiGraph{false}(graph, matching) +end + +""" + $(TYPEDSIGNATURES) + +Return the indexes of observed equations of `sys` used by expression `exprs`. +""" +function observed_equations_used_by(sys::AbstractSystem, exprs) + obs = observed(sys) + + obsvars = getproperty.(obs, :lhs) + graph = observed_dependency_graph(obs) + + syms = vars(exprs) + + obsidxs = BitSet() + for sym in syms + idx = findfirst(isequal(sym), obsvars) + idx === nothing && continue + parents = dfs_parents(graph, idx) + for i in eachindex(parents) + parents[i] == 0 && continue + push!(obsidxs, i) + end + end + + obsidxs = collect(obsidxs) + sort!(obsidxs) + return obsidxs +end + +""" + $(TYPEDSIGNATURES) + +Given an expression `expr`, return a dictionary mapping subexpressions of `expr` that do +not involve variables in `vars` to anonymous symbolic variables. Also return the modified +`expr` with the substitutions indicated by the dictionary. If `expr` is a function +of only `vars`, then all of the returned subexpressions can be precomputed. + +Note that this will only process subexpressions floating point value. Additionally, +array variables must be passed in both scalarized and non-scalarized forms in `vars`. +""" +function subexpressions_not_involving_vars(expr, vars) + expr = unwrap(expr) + vars = map(unwrap, vars) + state = Dict() + newexpr = subexpressions_not_involving_vars!(expr, vars, state) + return state, newexpr +end + +""" + $(TYPEDSIGNATURES) + +Mutating version of `subexpressions_not_involving_vars` which writes to `state`. Only +returns the modified `expr`. +""" +function subexpressions_not_involving_vars!(expr, vars, state::Dict{Any, Any}) + expr = unwrap(expr) + symbolic_type(expr) == NotSymbolic() && return expr + iscall(expr) || return expr + is_variable_floatingpoint(expr) || return expr + symtype(expr) <: Union{Real, AbstractArray{<:Real}} || return expr + Symbolics.shape(expr) == Symbolics.Unknown() && return expr + haskey(state, expr) && return state[expr] + vs = ModelingToolkit.vars(expr) + intersect!(vs, vars) + if isempty(vs) + sym = gensym(:subexpr) + stype = symtype(expr) + var = similar_variable(expr, sym) + state[expr] = var + return var + end + op = operation(expr) + args = arguments(expr) + if (op == (+) || op == (*)) && symbolic_type(expr) !== ArraySymbolic() + indep_args = [] + dep_args = [] + for arg in args + _vs = ModelingToolkit.vars(arg) + intersect!(_vs, vars) + if !isempty(_vs) + push!(dep_args, subexpressions_not_involving_vars!(arg, vars, state)) + else + push!(indep_args, arg) + end + end + indep_term = reduce(op, indep_args; init = Int(op == (*))) + indep_term = subexpressions_not_involving_vars!(indep_term, vars, state) + dep_term = reduce(op, dep_args; init = Int(op == (*))) + return op(indep_term, dep_term) + end + newargs = map(args) do arg + symbolic_type(arg) != NotSymbolic() || is_array_of_symbolics(arg) || return arg + subexpressions_not_involving_vars!(arg, vars, state) + end + return maketerm(typeof(expr), op, newargs, metadata(expr)) +end + +""" + $(TYPEDSIGNATURES) + +Create an anonymous symbolic variable of the same shape, size and symtype as `var`, with +name `gensym(name)`. Does not support unsized array symbolics. +""" +function similar_variable(var::BasicSymbolic, name = :anon) + name = gensym(name) + stype = symtype(var) + sym = Symbolics.variable(name; T = stype) + if size(var) !== () + sym = setmetadata(sym, Symbolics.ArrayShapeCtx, map(Base.OneTo, size(var))) + end + return sym +end diff --git a/test/initializationsystem.jl b/test/initializationsystem.jl index f3015f7db0..9c06dd2030 100644 --- a/test/initializationsystem.jl +++ b/test/initializationsystem.jl @@ -28,7 +28,7 @@ sol = solve(initprob) initprob = ModelingToolkit.InitializationProblem(pend, 0.0, [x => 1, y => 0], [g => 1]; guesses = ModelingToolkit.missing_variable_defaults(pend)) -@test initprob isa NonlinearProblem +@test initprob isa NonlinearLeastSquaresProblem sol = solve(initprob) @test SciMLBase.successful_retcode(sol) @test sol.u == [0.0, 0.0, 0.0, 0.0] diff --git a/test/mtkparameters.jl b/test/mtkparameters.jl index 603426aaf7..ce524fdb76 100644 --- a/test/mtkparameters.jl +++ b/test/mtkparameters.jl @@ -299,7 +299,7 @@ end # Parameter timeseries ps = MTKParameters(([1.0, 1.0],), (BlockedArray(zeros(4), [2, 2]),), - (), ()) + (), (), ()) ps2 = SciMLStructures.replace(Discrete(), ps, ones(4)) @test typeof(ps2.discrete) == typeof(ps.discrete) with_updated_parameter_timeseries_values( @@ -316,7 +316,7 @@ with_updated_parameter_timeseries_values( ps = MTKParameters( (), (BlockedArray([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], [3, 3]), BlockedArray(falses(1), [1, 0])), - (), ()) + (), (), ()) @test SciMLBase.get_saveable_values(sys, ps, 1).x isa Tuple{Vector{Float64}, Vector{Bool}} tsidx1 = 1 tsidx2 = 2 diff --git a/test/runtests.jl b/test/runtests.jl index 44846eed57..677f40c717 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,7 @@ end @safetestset "SDESystem Test" include("sdesystem.jl") @safetestset "DDESystem Test" include("dde.jl") @safetestset "NonlinearSystem Test" include("nonlinearsystem.jl") + @safetestset "SCCNonlinearProblem Test" include("scc_nonlinear_problem.jl") @safetestset "PDE Construction Test" include("pde.jl") @safetestset "JumpSystem Test" include("jumpsystem.jl") @safetestset "print_tree" include("print_tree.jl") diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl new file mode 100644 index 0000000000..fdf1646343 --- /dev/null +++ b/test/scc_nonlinear_problem.jl @@ -0,0 +1,163 @@ +using ModelingToolkit +using NonlinearSolve, SCCNonlinearSolve +using OrdinaryDiffEq +using SciMLBase, Symbolics +using LinearAlgebra, Test +using ModelingToolkit: t_nounits as t, D_nounits as D + +@testset "Trivial case" begin + function f!(du, u, p) + du[1] = cos(u[2]) - u[1] + du[2] = sin(u[1] + u[2]) + u[2] + du[3] = 2u[4] + u[3] + 1.0 + du[4] = u[5]^2 + u[4] + du[5] = u[3]^2 + u[5] + du[6] = u[1] + u[2] + u[3] + u[4] + u[5] + 2.0u[6] + 2.5u[7] + 1.5u[8] + du[7] = u[1] + u[2] + u[3] + 2.0u[4] + u[5] + 4.0u[6] - 1.5u[7] + 1.5u[8] + du[8] = u[1] + 2.0u[2] + 3.0u[3] + 5.0u[4] + 6.0u[5] + u[6] - u[7] - u[8] + end + @variables u[1:8] [irreducible = true] + eqs = Any[0 for _ in 1:8] + f!(eqs, u, nothing) + eqs = 0 .~ eqs + @named model = NonlinearSystem(eqs) + @test_throws ["simplified", "required"] SCCNonlinearProblem(model, []) + _model = structural_simplify(model; split = false) + @test_throws ["not compatible"] SCCNonlinearProblem(_model, []) + model = structural_simplify(model) + prob = NonlinearProblem(model, [u => zeros(8)]) + sccprob = SCCNonlinearProblem(model, [u => zeros(8)]) + sol1 = solve(prob, NewtonRaphson()) + sol2 = solve(sccprob, NewtonRaphson()) + @test SciMLBase.successful_retcode(sol1) + @test SciMLBase.successful_retcode(sol2) + @test sol1[u] ≈ sol2[u] +end + +@testset "With parameters" begin + function f!(du, u, (p1, p2), t) + x = (*)(p1[4], u[1]) + y = (*)(p1[4], (+)(0.1016, (*)(-1, u[1]))) + z1 = ifelse((<)(p2[1], 0), + (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))), + 0) + z2 = ifelse((>)(p2[1], 0), + (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[4]), + 0) + z3 = ifelse((>)(p2[1], 0), + (*)((*)(457896.07999999996, p1[2]), sqrt((*)(1.1686468413521012e-5, p1[3]))), + 0) + z4 = ifelse((<)(p2[1], 0), + (*)((*)((*)(0.58, p1[2]), sqrt((*)(1 // 86100, p1[3]))), u[5]), + 0) + du[1] = p2[1] + du[2] = (+)(z1, (*)(-1, z2)) + du[3] = (+)(z3, (*)(-1, z4)) + du[4] = (+)((*)(-1, u[2]), (*)((*)(1 // 86100, y), u[4])) + du[5] = (+)((*)(-1, u[3]), (*)((*)(1 // 86100, x), u[5])) + end + p = ( + [0.04864391799335977, 7.853981633974484e-5, 1.4034843205574914, + 0.018241469247509915, 300237.05, 9.226186337232914], + [0.0508]) + u0 = [0.0, 0.0, 0.0, 789476.0, 101325.0] + tspan = (0.0, 1.0) + mass_matrix = [1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0; + 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0] + dt = 1e-3 + function nlf(u1, (u0, p)) + resid = Any[0 for _ in u0] + f!(resid, u1, p, 0.0) + return mass_matrix * (u1 - u0) - dt * resid + end + + prob = NonlinearProblem(nlf, u0, (u0, p)) + @test_throws Exception solve(prob, SimpleNewtonRaphson(), abstol = 1e-9) + sol = solve(prob, TrustRegion(); abstol = 1e-9) + + @variables u[1:5] [irreducible = true] + @parameters p1[1:6] p2 + eqs = 0 .~ collect(nlf(u, (u0, (p1, p2)))) + @mtkbuild sys = NonlinearSystem(eqs, [u], [p1, p2]) + sccprob = SCCNonlinearProblem(sys, [u => u0], [p1 => p[1], p2 => p[2][]]) + sccsol = solve(sccprob, SimpleNewtonRaphson(); abstol = 1e-9) + @test SciMLBase.successful_retcode(sccsol) + @test norm(sccsol.resid) < norm(sol.resid) +end + +@testset "Transistor amplifier" begin + C = [k * 1e-6 for k in 1:5] + Ub = 6 + UF = 0.026 + α = 0.99 + β = 1e-6 + R0 = 1000 + R = 9000 + Ue(t) = 0.1 * sin(200 * π * t) + + function transamp(out, du, u, p, t) + g(x) = 1e-6 * (exp(x / 0.026) - 1) + y1, y2, y3, y4, y5, y6, y7, y8 = u + out[1] = -Ue(t) / R0 + y1 / R0 + C[1] * du[1] - C[1] * du[2] + out[2] = -Ub / R + y2 * 2 / R - (α - 1) * g(y2 - y3) - C[1] * du[1] + C[1] * du[2] + out[3] = -g(y2 - y3) + y3 / R + C[2] * du[3] + out[4] = -Ub / R + y4 / R + α * g(y2 - y3) + C[3] * du[4] - C[3] * du[5] + out[5] = -Ub / R + y5 * 2 / R - (α - 1) * g(y5 - y6) - C[3] * du[4] + C[3] * du[5] + out[6] = -g(y5 - y6) + y6 / R + C[4] * du[6] + out[7] = -Ub / R + y7 / R + α * g(y5 - y6) + C[5] * du[7] - C[5] * du[8] + out[8] = y8 / R - C[5] * du[7] + C[5] * du[8] + end + + u0 = [0, Ub / 2, Ub / 2, Ub, Ub / 2, Ub / 2, Ub, 0] + du0 = [ + 51.338775, + 51.338775, + -Ub / (2 * (C[2] * R)), + -24.9757667, + -24.9757667, + -Ub / (2 * (C[4] * R)), + -10.00564453, + -10.00564453 + ] + daeprob = DAEProblem(transamp, du0, u0, (0.0, 0.1)) + daesol = solve(daeprob, DImplicitEuler()) + + t0 = daesol.t[5] + t1 = daesol.t[6] + u0 = daesol.u[5] + u1 = daesol.u[6] + dt = t1 - t0 + + @variables y(t)[1:8] + eqs = Any[0 for _ in 1:8] + transamp(eqs, collect(D(y)), y, nothing, t) + eqs = 0 .~ eqs + subrules = Dict(Symbolics.unwrap(D(y[i])) => ((y[i] - u0[i]) / dt) for i in 1:8) + eqs = substitute.(eqs, (subrules,)) + @mtkbuild sys = NonlinearSystem(eqs) + prob = NonlinearProblem(sys, [y => u0], [t => t0]) + sol = solve(prob, NewtonRaphson(); abstol = 1e-12) + + sccprob = SCCNonlinearProblem(sys, [y => u0], [t => t0]) + sccsol = solve(sccprob, NewtonRaphson(); abstol = 1e-12) + + @test sol.u≈sccsol.u atol=1e-10 +end + +@testset "Expression caching" begin + @variables x[1:4] = rand(4) + val = Ref(0) + function func(x, y) + val[] += 1 + x + y + end + @register_symbolic func(x, y) + @mtkbuild sys = NonlinearSystem([0 ~ x[1]^3 + x[2]^3 - 5 + 0 ~ sin(x[1] - x[2]) - 0.5 + 0 ~ func(x[1], x[2]) * exp(x[3]) - x[4]^3 - 5 + 0 ~ func(x[1], x[2]) * exp(x[4]) - x[3]^3 - 4]) + sccprob = SCCNonlinearProblem(sys, []) + sccsol = solve(sccprob, NewtonRaphson()) + @test SciMLBase.successful_retcode(sccsol) + @test val[] == 1 +end diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 22c90edf7a..2f8667faa8 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -206,7 +206,7 @@ S = get_sensitivity(closed_loop, :u) BlockedArray([[1 2; 3 4], [2 4; 6 8]], [1, 1])), # (BlockedArray([[true, false], [false, true]]), BlockedArray([[[1 2; 3 4]], [[2 4; 6 8]]])), ([5, 6],), - (["hi", "bye"], [:lie, :die])) + (["hi", "bye"], [:lie, :die]), ()) @test ps[ParameterIndex(Tunable(), 1)] == 1.0 @test ps[ParameterIndex(Tunable(), 2:4)] == collect(2.0:4.0) @test ps[ParameterIndex(Tunable(), reshape(4:7, 2, 2))] == reshape(4.0:7.0, 2, 2) diff --git a/test/symbolic_events.jl b/test/symbolic_events.jl index 61690acdce..9d3aac4074 100644 --- a/test/symbolic_events.jl +++ b/test/symbolic_events.jl @@ -1063,7 +1063,7 @@ end cb = [x ~ 0.0] => [x ~ 0, y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb]) prob = ODEProblem(pend, [x => 1], (0.0, 3.0), guesses = [y => x]) - @test_throws "CheckInit specified but initialization" solve(prob, Rodas5()) + @test_throws "DAE initialization failed" solve(prob, Rodas5()) cb = [x ~ 0.0] => [y ~ 1] @mtkbuild pend = ODESystem(eqs, t; continuous_events = [cb])