From e56de20aa260ffbd29bbacddcb98389c595996f7 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 25 Sep 2023 10:18:54 -0400 Subject: [PATCH 01/16] split parameters linearize bug --- test/split_parameters.jl | 73 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/test/split_parameters.jl b/test/split_parameters.jl index ef8f434dca..11608e1a7f 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -77,3 +77,76 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false) @test prob.p isa Tuple{Vector{Float64}, Vector{Int64}} sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success + + + + + + +# ------------------------- Bug +using ModelingToolkit, LinearAlgebra +using ModelingToolkitStandardLibrary.Mechanical.Rotational +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkitStandardLibrary.Blocks: t +using ModelingToolkit: connect + +"A wrapper function to make symbolic indexing easier" +function wr(sys) + ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems=[sys], name=:a_wrapper) +end +indexof(sym,syms) = findfirst(isequal(sym),syms) + +# Parameters +m1 = 1. +m2 = 1. +k = 10. # Spring stiffness +c = 3. # Damping coefficient + +@named inertia1 = Inertia(; J = m1) +@named inertia2 = Inertia(; J = m2) +@named spring = Spring(; c = k) +@named damper = Damper(; d = c) +@named torque = Torque(use_support=false) + +function SystemModel(u=nothing; name=:model) + eqs = [ + connect(torque.flange, inertia1.flange_a) + connect(inertia1.flange_b, spring.flange_a, damper.flange_a) + connect(inertia2.flange_a, spring.flange_b, damper.flange_b) + ] + if u !== nothing + push!(eqs, connect(torque.tau, u.output)) + return @named model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper, u]) + end + ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name) +end + + +model = SystemModel() # Model with load disturbance +@named d = Step(start_time=1., duration=10., offset=0., height=1.) # Disturbance +model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control +inputs = [model.torque.tau.u] +matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs) + +# Design state-feedback gain using LQR +# Define cost matrices +x_costs = [ + model.inertia1.w => 1. + model.inertia2.w => 1. + model.inertia1.phi => 1. + model.inertia2.phi => 1. +] +L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller +@named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix +@named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance + +connections = [ + [state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4] + connect(d.output, :d, add.input1) + connect(add.input2, state_feedback.output) + connect(add.output, :u, model.torque.tau) +] +closed_loop = ODESystem(connections, t, systems=[model, state_feedback, add, d], name=:closed_loop) +S = get_sensitivity(closed_loop, :u) + + From cee69ed07ad7b04dbb5977b9053b59d59bb86fed Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 25 Sep 2023 13:31:30 -0400 Subject: [PATCH 02/16] added old MatrixGain definition for comparison --- test/split_parameters.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 11608e1a7f..3ce5f7e2f0 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -137,6 +137,16 @@ x_costs = [ model.inertia2.phi => 1. ] L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller + +# This old definition of MatrixGain will work because the parameter space does not include K (an Array term) +# @component function MatrixGainAlt(K::AbstractArray; name) +# nout, nin = size(K, 1), size(K, 2) +# @named input = RealInput(; nin = nin) +# @named output = RealOutput(; nout = nout) +# eqs = [output.u[i] ~ sum(K[i, j] * input.u[j] for j in 1:nin) for i in 1:nout] +# compose(ODESystem(eqs, t, [], []; name = name), [input, output]) +# end + @named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix @named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance From b4f821eeb21b4e5ac6f92a0c522e0abd3362544d Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Mon, 25 Sep 2023 19:14:07 -0400 Subject: [PATCH 03/16] new promote_to_concrete --- src/utils.jl | 56 ++++++++++++++++++---------------------- test/split_parameters.jl | 28 ++++++++++++++++++++ 2 files changed, 53 insertions(+), 31 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index abea65ab21..452323c6ff 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -657,46 +657,40 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) end T = eltype(vs) if Base.isconcretetype(T) && (!tofloat || T === float(T)) # nothing to do - vs + return vs else sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs) isempty(sym_vs) || throw_missingvars_in_sys(sym_vs) - C = typeof(first(vs)) - I = Int8 - has_int = false - has_array = false - has_bool = false - array_T = nothing + + C = nothing for v in vs - if v isa AbstractArray - has_array = true - array_T = typeof(v) - end - E = eltype(v) - C = promote_type(C, E) - if E <: Integer - has_int = true - I = promote_type(I, E) + E = typeof(v) + if E <: Number + if tofloat + E = float(E) + end end - if E <: Bool - has_bool = true + if use_union + if C === nothing + C = E + else + C = Union{C, E} + end + else + C = E end end - if tofloat && !has_array - C = float(C) - elseif has_array || (use_union && has_int && C !== I) - if has_array - C = Union{C, array_T} - end - if has_int - C = Union{C, I} - end - if has_bool - C = Union{C, Bool} + + y = similar(vs, C) + for i in eachindex(vs) + if (vs[i] isa Number) & tofloat + y[i] = float(vs[i]) + else + y[i] = vs[i] end - return copyto!(similar(vs, C), vs) end - convert.(C, vs) + + return y end end diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 3ce5f7e2f0..1ea7161ee0 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -2,6 +2,34 @@ using ModelingToolkit, Test using ModelingToolkitStandardLibrary.Blocks using OrdinaryDiffEq + +x = [1, 2.0, false, [1,2,3], Parameter(1.0)] + +y = ModelingToolkit.promote_to_concrete(x) +@test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}} + +y = ModelingToolkit.promote_to_concrete(x; tofloat=false) +@test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}} + + +x = [1, 2.0, false, [1,2,3]] +y = ModelingToolkit.promote_to_concrete(x) +@test eltype(y) == Union{Float64, Vector{Int64}} + +x = Any[1, 2.0, false] +y = ModelingToolkit.promote_to_concrete(x; tofloat=false) +@test eltype(y) == Union{Bool, Float64, Int64} + +y = ModelingToolkit.promote_to_concrete(x; use_union=false) +@test eltype(y) == Float64 + +x = Float16[1., 2., 3.] +y = ModelingToolkit.promote_to_concrete(x) +@test eltype(y) == Float16 + + + + # ------------------------ Mixed Single Values and Vector dt = 4e-4 From e62e3ee4562414a662cb59a65d0bf90be2ad4c26 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Tue, 26 Sep 2023 05:58:01 -0400 Subject: [PATCH 04/16] improvements to protmote_to_concrete --- src/utils.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 452323c6ff..0a6c5089a3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -670,13 +670,13 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) E = float(E) end end + if C === nothing + C = E + end if use_union - if C === nothing - C = E - else - C = Union{C, E} - end + C = Union{C, E} else + @assert C == E "`promote_to_concrete` can't make type $E uniform with $C" C = E end end @@ -684,9 +684,9 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) y = similar(vs, C) for i in eachindex(vs) if (vs[i] isa Number) & tofloat - y[i] = float(vs[i]) + y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically else - y[i] = vs[i] + y[i] = vs[i] end end From ee4deb98dcc7cf83d3ccb6d4cf0edb544614af92 Mon Sep 17 00:00:00 2001 From: Brad Carman Date: Tue, 26 Sep 2023 05:58:14 -0400 Subject: [PATCH 05/16] remake issue --- test/split_parameters.jl | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 1ea7161ee0..3f231e2dd7 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -74,12 +74,28 @@ eqs = [y ~ src.output.u @named sys = ODESystem(eqs, t, vars, []; systems = [int, src]) s = complete(sys) sys = structural_simplify(sys) -prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]) +prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat=false) @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}} sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success @test sol[y][end] == x[end] +#TODO: remake becomes more complicated now, how to improve? +defs = ModelingToolkit.defaults(sys) +defs[s.src.data] = 2x +p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false) +p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake() +prob′ = remake(prob; p=p′) +sol = solve(prob′, ImplicitEuler()); +@test sol.retcode == ReturnCode.Success +@test sol[y][end] == 2x[end] + +prob′′ = remake(prob; p=[s.src.data => x]) +@test prob′′.p isa Tuple + + + + # ------------------------ Mixed Type Converted to float (default behavior) vars = @variables y(t)=1 dy(t)=0 ddy(t)=0 From a9c56b616b0a1bb57e1b8ba96fa81045eee03e4e Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 26 Sep 2023 10:10:19 -0400 Subject: [PATCH 06/16] Suppose heterogeneous parameters for linearize and remake --- Project.toml | 2 +- src/systems/abstractsystem.jl | 18 ++++- src/systems/diffeqs/abstractodesystem.jl | 10 ++- src/systems/diffeqs/odesystem.jl | 10 ++- src/utils.jl | 6 +- src/variables.jl | 26 ++++--- test/split_parameters.jl | 90 +++++++++--------------- 7 files changed, 85 insertions(+), 77 deletions(-) diff --git a/Project.toml b/Project.toml index 58d73f2921..f5035f44cb 100644 --- a/Project.toml +++ b/Project.toml @@ -85,7 +85,7 @@ PrecompileTools = "1" RecursiveArrayTools = "2.3" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "2.0.1" +SciMLBase = "1, 2.0.1" Setfield = "0.7, 0.8, 1" SimpleNonlinearSolve = "0.1.0" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 48d647f621..df0ff4773f 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -230,7 +230,8 @@ for prop in [:eqs :metadata :gui_metadata :discrete_subsystems - :unknown_states] + :unknown_states + :split_idxs] fname1 = Symbol(:get_, prop) fname2 = Symbol(:has_, prop) @eval begin @@ -1274,14 +1275,25 @@ See also [`linearize`](@ref) which provides a higher-level interface. function linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, initialize = true, + op = Dict(), + p = DiffEqBase.NullParameters(), kwargs...) sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify, kwargs...) + x0 = merge(defaults(sys), op) + u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true) + p, split_idxs = split_parameters_by_type(p) + ps = parameters(sys) + if p isa Tuple + ps = Base.Fix1(getindex, ps).(split_idxs) + ps = (ps...,) #if p is Tuple, ps should be Tuple + end + lin_fun = let diff_idxs = diff_idxs, alge_idxs = alge_idxs, input_idxs = input_idxs, sts = states(sys), - fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys), + fun = ODEFunction{true, SciMLBase.FullSpecialize}(sys, states(sys), ps; p = p), h = build_explicit_observed_function(sys, outputs), chunk = ForwardDiff.Chunk(input_idxs) @@ -1600,11 +1612,11 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0, allow_input_derivatives = false, zero_dummy_der = false, kwargs...) - lin_fun, ssys = linearization_function(sys, inputs, outputs; kwargs...) if zero_dummy_der dummyder = setdiff(states(ssys), states(sys)) op = merge(op, Dict(x => 0.0 for x in dummyder)) end + lin_fun, ssys = linearization_function(sys, inputs, outputs; op, kwargs...) linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys end diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 8528dd3c12..f11d1283c6 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -354,6 +354,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s checkbounds = false, sparsity = false, analytic = nothing, + split_idxs = nothing, kwargs...) where {iip, specialize} f_gen = generate_function(sys, dvs, ps; expression = Val{eval_expression}, expression_module = eval_module, checkbounds = checkbounds, @@ -508,6 +509,7 @@ function DiffEqBase.ODEFunction{iip, specialize}(sys::AbstractODESystem, dvs = s nothing end + @set! sys.split_idxs = split_idxs ODEFunction{iip, specialize}(f; sys = sys, jac = _jac === nothing ? nothing : _jac, @@ -765,7 +767,7 @@ Take dictionaries with initial conditions and parameters and convert them to num """ function get_u0_p(sys, u0map, - parammap; + parammap = nothing; use_union = true, tofloat = true, symbolic_u0 = false) @@ -773,7 +775,9 @@ function get_u0_p(sys, ps = parameters(sys) defs = defaults(sys) - defs = mergedefaults(defs, parammap, ps) + if parammap !== nothing + defs = mergedefaults(defs, parammap, ps) + end defs = mergedefaults(defs, u0map, dvs) if symbolic_u0 @@ -835,7 +839,7 @@ function process_DEProblem(constructor, sys::AbstractODESystem, u0map, parammap; f = constructor(sys, dvs, ps, u0; ddvs = ddvs, tgrad = tgrad, jac = jac, checkbounds = checkbounds, p = p, linenumbers = linenumbers, parallel = parallel, simplify = simplify, - sparse = sparse, eval_expression = eval_expression, + sparse = sparse, eval_expression = eval_expression, split_idxs, kwargs...) implicit_dae ? (f, du0, u0, p) : (f, u0, p) end diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index dcbecfcfb0..3b828702c2 100644 --- a/src/systems/diffeqs/odesystem.jl +++ b/src/systems/diffeqs/odesystem.jl @@ -139,6 +139,10 @@ struct ODESystem <: AbstractODESystem used for ODAEProblem. """ unknown_states::Union{Nothing, Vector{Any}} + """ + split_idxs: a vector of vectors of indices for the split parameters. + """ + split_idxs::Union{Nothing, Vector{Vector{Int}}} function ODESystem(tag, deqs, iv, dvs, ps, tspan, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, @@ -146,8 +150,8 @@ struct ODESystem <: AbstractODESystem devents, metadata = nothing, gui_metadata = nothing, tearing_state = nothing, substitutions = nothing, complete = false, - discrete_subsystems = nothing, unknown_states = nothing; - checks::Union{Bool, Int} = true) + discrete_subsystems = nothing, unknown_states = nothing, + split_idxs = nothing; checks::Union{Bool, Int} = true) if checks == true || (checks & CheckComponents) > 0 check_variables(dvs, iv) check_parameters(ps, iv) @@ -161,7 +165,7 @@ struct ODESystem <: AbstractODESystem ctrl_jac, Wfact, Wfact_t, name, systems, defaults, torn_matching, connector_type, preface, cevents, devents, metadata, gui_metadata, tearing_state, substitutions, complete, discrete_subsystems, - unknown_states) + unknown_states, split_idxs) end end diff --git a/src/utils.jl b/src/utils.jl index 0a6c5089a3..ab17bd3d8d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -661,7 +661,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) else sym_vs = filter(x -> SymbolicUtils.issym(x) || SymbolicUtils.istree(x), vs) isempty(sym_vs) || throw_missingvars_in_sys(sym_vs) - + C = nothing for v in vs E = typeof(v) @@ -676,7 +676,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) if use_union C = Union{C, E} else - @assert C == E "`promote_to_concrete` can't make type $E uniform with $C" + @assert C==E "`promote_to_concrete` can't make type $E uniform with $C" C = E end end @@ -686,7 +686,7 @@ function promote_to_concrete(vs; tofloat = true, use_union = true) if (vs[i] isa Number) & tofloat y[i] = float(vs[i]) #needed because copyto! can't convert Int to Float automatically else - y[i] = vs[i] + y[i] = vs[i] end end diff --git a/src/variables.jl b/src/variables.jl index 4cffd3fdeb..854680998e 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -145,15 +145,23 @@ function SciMLBase.process_p_u0_symbolic(prob::Union{SciMLBase.AbstractDEProblem " Please use `remake` with the `u0` keyword argument as a vector of values, paying attention to state order.")) end - # assemble defaults - defs = defaults(prob.f.sys) - defs = mergedefaults(defs, prob.p, parameters(prob.f.sys)) - defs = mergedefaults(defs, p, parameters(prob.f.sys)) - defs = mergedefaults(defs, prob.u0, states(prob.f.sys)) - defs = mergedefaults(defs, u0, states(prob.f.sys)) - - u0 = varmap_to_vars(u0, states(prob.f.sys); defaults = defs, tofloat = true) - p = varmap_to_vars(p, parameters(prob.f.sys); defaults = defs) + sys = prob.f.sys + defs = defaults(sys) + ps = parameters(sys) + if has_split_idxs(sys) && (split_idxs = get_split_idxs(sys)) !== nothing + for (i, idxs) in enumerate(split_idxs) + defs = mergedefaults(defs, prob.p[i], ps[idxs]) + end + else + # assemble defaults + defs = defaults(sys) + defs = mergedefaults(defs, prob.p, ps) + end + defs = mergedefaults(defs, p, ps) + sts = states(sys) + defs = mergedefaults(defs, prob.u0, sts) + defs = mergedefaults(defs, u0, sts) + u0, p, defs = get_u0_p(sys, defs) return p, u0 end diff --git a/test/split_parameters.jl b/test/split_parameters.jl index 3f231e2dd7..d37b9e2c13 100644 --- a/test/split_parameters.jl +++ b/test/split_parameters.jl @@ -2,34 +2,29 @@ using ModelingToolkit, Test using ModelingToolkitStandardLibrary.Blocks using OrdinaryDiffEq - -x = [1, 2.0, false, [1,2,3], Parameter(1.0)] +x = [1, 2.0, false, [1, 2, 3], Parameter(1.0)] y = ModelingToolkit.promote_to_concrete(x) @test eltype(y) == Union{Float64, Parameter{Float64}, Vector{Int64}} -y = ModelingToolkit.promote_to_concrete(x; tofloat=false) +y = ModelingToolkit.promote_to_concrete(x; tofloat = false) @test eltype(y) == Union{Bool, Float64, Int64, Parameter{Float64}, Vector{Int64}} - -x = [1, 2.0, false, [1,2,3]] +x = [1, 2.0, false, [1, 2, 3]] y = ModelingToolkit.promote_to_concrete(x) @test eltype(y) == Union{Float64, Vector{Int64}} x = Any[1, 2.0, false] -y = ModelingToolkit.promote_to_concrete(x; tofloat=false) +y = ModelingToolkit.promote_to_concrete(x; tofloat = false) @test eltype(y) == Union{Bool, Float64, Int64} -y = ModelingToolkit.promote_to_concrete(x; use_union=false) +y = ModelingToolkit.promote_to_concrete(x; use_union = false) @test eltype(y) == Float64 -x = Float16[1., 2., 3.] +x = Float16[1.0, 2.0, 3.0] y = ModelingToolkit.promote_to_concrete(x) @test eltype(y) == Float16 - - - # ------------------------ Mixed Single Values and Vector dt = 4e-4 @@ -74,7 +69,7 @@ eqs = [y ~ src.output.u @named sys = ODESystem(eqs, t, vars, []; systems = [int, src]) s = complete(sys) sys = structural_simplify(sys) -prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat=false) +prob = ODEProblem(sys, [], (0.0, t_end), [s.src.data => x]; tofloat = false) @test prob.p isa Tuple{Vector{Float64}, Vector{Int}, Vector{Vector{Float64}}} sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success @@ -83,18 +78,15 @@ sol = solve(prob, ImplicitEuler()); #TODO: remake becomes more complicated now, how to improve? defs = ModelingToolkit.defaults(sys) defs[s.src.data] = 2x -p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat=false) +p′ = ModelingToolkit.varmap_to_vars(defs, parameters(sys); tofloat = false) p′, = ModelingToolkit.split_parameters_by_type(p′) #NOTE: we need to ensure this is called now before calling remake() -prob′ = remake(prob; p=p′) +prob′ = remake(prob; p = p′) sol = solve(prob′, ImplicitEuler()); @test sol.retcode == ReturnCode.Success @test sol[y][end] == 2x[end] -prob′′ = remake(prob; p=[s.src.data => x]) -@test prob′′.p isa Tuple - - - +prob′′ = remake(prob; p = [s.src.data => x]) +@test_broken prob′′.p isa Tuple # ------------------------ Mixed Type Converted to float (default behavior) @@ -122,11 +114,6 @@ prob = ODEProblem(sys, [], tspan, []; tofloat = false) sol = solve(prob, ImplicitEuler()); @test sol.retcode == ReturnCode.Success - - - - - # ------------------------- Bug using ModelingToolkit, LinearAlgebra using ModelingToolkitStandardLibrary.Mechanical.Rotational @@ -136,51 +123,48 @@ using ModelingToolkit: connect "A wrapper function to make symbolic indexing easier" function wr(sys) - ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems=[sys], name=:a_wrapper) + ODESystem(Equation[], ModelingToolkit.get_iv(sys), systems = [sys], name = :a_wrapper) end -indexof(sym,syms) = findfirst(isequal(sym),syms) +indexof(sym, syms) = findfirst(isequal(sym), syms) # Parameters -m1 = 1. -m2 = 1. -k = 10. # Spring stiffness -c = 3. # Damping coefficient +m1 = 1.0 +m2 = 1.0 +k = 10.0 # Spring stiffness +c = 3.0 # Damping coefficient @named inertia1 = Inertia(; J = m1) @named inertia2 = Inertia(; J = m2) @named spring = Spring(; c = k) @named damper = Damper(; d = c) -@named torque = Torque(use_support=false) +@named torque = Torque(use_support = false) -function SystemModel(u=nothing; name=:model) - eqs = [ - connect(torque.flange, inertia1.flange_a) +function SystemModel(u = nothing; name = :model) + eqs = [connect(torque.flange, inertia1.flange_a) connect(inertia1.flange_b, spring.flange_a, damper.flange_a) - connect(inertia2.flange_a, spring.flange_b, damper.flange_b) - ] + connect(inertia2.flange_a, spring.flange_b, damper.flange_b)] if u !== nothing push!(eqs, connect(torque.tau, u.output)) - return @named model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper, u]) + return @named model = ODESystem(eqs, + t; + systems = [torque, inertia1, inertia2, spring, damper, u]) end ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name) end - model = SystemModel() # Model with load disturbance -@named d = Step(start_time=1., duration=10., offset=0., height=1.) # Disturbance +@named d = Step(start_time = 1.0, duration = 10.0, offset = 0.0, height = 1.0) # Disturbance model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi] # This is the state realization we want to control inputs = [model.torque.tau.u] matrices, ssys = ModelingToolkit.linearize(wr(model), inputs, model_outputs) # Design state-feedback gain using LQR # Define cost matrices -x_costs = [ - model.inertia1.w => 1. - model.inertia2.w => 1. - model.inertia1.phi => 1. - model.inertia2.phi => 1. -] -L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller +x_costs = [model.inertia1.w => 1.0 + model.inertia2.w => 1.0 + model.inertia1.phi => 1.0 + model.inertia2.phi => 1.0] +L = randn(1, 4) # Post-multiply by `C` to get the correct input to the controller # This old definition of MatrixGain will work because the parameter space does not include K (an Array term) # @component function MatrixGainAlt(K::AbstractArray; name) @@ -191,16 +175,12 @@ L = randn(1,4) # Post-multiply by `C` to get the correct input to the controller # compose(ODESystem(eqs, t, [], []; name = name), [input, output]) # end -@named state_feedback = MatrixGain(K=-L) # Build negative feedback into the feedback matrix -@named add = Add(;k1=1., k2=1.) # To add the control signal and the disturbance +@named state_feedback = MatrixGain(K = -L) # Build negative feedback into the feedback matrix +@named add = Add(; k1 = 1.0, k2 = 1.0) # To add the control signal and the disturbance -connections = [ - [state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4] +connections = [[state_feedback.input.u[i] ~ model_outputs[i] for i in 1:4] connect(d.output, :d, add.input1) connect(add.input2, state_feedback.output) - connect(add.output, :u, model.torque.tau) -] -closed_loop = ODESystem(connections, t, systems=[model, state_feedback, add, d], name=:closed_loop) + connect(add.output, :u, model.torque.tau)] +@named closed_loop = ODESystem(connections, t, systems = [model, state_feedback, add, d]) S = get_sensitivity(closed_loop, :u) - - From 1a96915f65de389430566c990e8aa308a4ae523e Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 26 Sep 2023 10:42:49 -0400 Subject: [PATCH 07/16] Fix zero_dummy_der Co-authored-by: Fredrik Bagge Carlson --- src/systems/abstractsystem.jl | 22 ++++++++++++++++------ test/input_output_handling.jl | 4 ++-- 2 files changed, 18 insertions(+), 8 deletions(-) diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index df0ff4773f..c452417409 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -1277,9 +1277,18 @@ function linearization_function(sys::AbstractSystem, inputs, initialize = true, op = Dict(), p = DiffEqBase.NullParameters(), + zero_dummy_der = false, kwargs...) - sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify, + ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; + simplify, kwargs...) + if zero_dummy_der + dummyder = setdiff(states(ssys), states(sys)) + defs = Dict(x => 0.0 for x in dummyder) + @set! ssys.defaults = merge(defs, defaults(ssys)) + op = merge(defs, op) + end + sys = ssys x0 = merge(defaults(sys), op) u0, p, _ = get_u0_p(sys, x0, p; use_union = false, tofloat = true) p, split_idxs = split_parameters_by_type(p) @@ -1612,11 +1621,12 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0, allow_input_derivatives = false, zero_dummy_der = false, kwargs...) - if zero_dummy_der - dummyder = setdiff(states(ssys), states(sys)) - op = merge(op, Dict(x => 0.0 for x in dummyder)) - end - lin_fun, ssys = linearization_function(sys, inputs, outputs; op, kwargs...) + lin_fun, ssys = linearization_function(sys, + inputs, + outputs; + zero_dummy_der, + op, + kwargs...) linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys end diff --git a/test/input_output_handling.jl b/test/input_output_handling.jl index 12f68f9f04..249a94e9d0 100644 --- a/test/input_output_handling.jl +++ b/test/input_output_handling.jl @@ -124,8 +124,8 @@ using ModelingToolkitStandardLibrary.Mechanical.Rotational t = ModelingToolkitStandardLibrary.Mechanical.Rotational.t @named inertia1 = Inertia(; J = 1) @named inertia2 = Inertia(; J = 1) -@named spring = Spring(; c = 10) -@named damper = Damper(; d = 3) +@named spring = Rotational.Spring(; c = 10) +@named damper = Rotational.Damper(; d = 3) @named torque = Torque(; use_support = false) @variables y(t) = 0 eqs = [connect(torque.flange, inertia1.flange_a) From 023e024f3ebd82991f34ea3900b9689572ac80ae Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 26 Sep 2023 14:25:27 -0400 Subject: [PATCH 08/16] Relax tests --- test/odesystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/odesystem.jl b/test/odesystem.jl index 56d436d76c..336eca11e7 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -241,7 +241,7 @@ p2 = (k₁ => 0.04, k₃ => 1e4) tspan = (0.0, 100000.0) prob1 = ODEProblem(sys, u0, tspan, p) -@test prob1.f.sys === sys +@test prob1.f.sys == sys prob12 = ODEProblem(sys, u0, tspan, [0.04, 3e7, 1e4]) prob13 = ODEProblem(sys, u0, tspan, (0.04, 3e7, 1e4)) prob14 = ODEProblem(sys, u0, tspan, p2) From 119cb27546c31ed22309240c72d14a952926f7b6 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Tue, 26 Sep 2023 16:50:06 -0400 Subject: [PATCH 09/16] Don't use SciMLBase 2 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f5035f44cb..bafa1835a3 100644 --- a/Project.toml +++ b/Project.toml @@ -85,7 +85,7 @@ PrecompileTools = "1" RecursiveArrayTools = "2.3" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "1, 2.0.1" +SciMLBase = "1" Setfield = "0.7, 0.8, 1" SimpleNonlinearSolve = "0.1.0" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" From 29aec7b3d19f6f333c863956dba07d02b070f5bf Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 2 Oct 2023 09:38:17 -0400 Subject: [PATCH 10/16] Fix eq ordering --- test/stream_connectors.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl index 9502b55be9..45d58c8552 100644 --- a/test/stream_connectors.jl +++ b/test/stream_connectors.jl @@ -157,8 +157,8 @@ eqns = [domain_connect(fluid, n1m1.port_a) pipe.port_b.P ~ sink.port.P] @test sort(equations(expand_connections(n1m1Test)), by = string) == [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow - 0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow 0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow + 0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow 0 ~ pipe.port_b.m_flow + sink.port.m_flow fluid.m_flow ~ 0 n1m1.port_a.P ~ pipe.port_a.P From addf3357af1f173a876fa735c473dd2b3ce5cba1 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 2 Oct 2023 09:53:16 -0400 Subject: [PATCH 11/16] Change bounds back --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a52124a694..470adc92b6 100644 --- a/Project.toml +++ b/Project.toml @@ -85,7 +85,7 @@ PrecompileTools = "1" RecursiveArrayTools = "2.3" Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" -SciMLBase = "1" +SciMLBase = "2.0.1" Setfield = "0.7, 0.8, 1" SimpleNonlinearSolve = "0.1.0" SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0, 2" From 2e4c1bbd6f31249a8f7045c68cde346885b04025 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 2 Oct 2023 09:55:29 -0400 Subject: [PATCH 12/16] Format --- docs/src/basics/MTKModel_Connector.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/src/basics/MTKModel_Connector.md b/docs/src/basics/MTKModel_Connector.md index 8cc106b05d..ac285d0253 100644 --- a/docs/src/basics/MTKModel_Connector.md +++ b/docs/src/basics/MTKModel_Connector.md @@ -130,6 +130,7 @@ end `@connector`s accepts begin blocks of `@components`, `@equations`, `@extend`, `@parameters`, `@structural_parameters`, `@variables`. These keywords mean the same as described above for `@mtkmodel`. !!! note + For more examples of usage, checkout [ModelingToolkitStandardLibrary.jl](https://github.com/SciML/ModelingToolkitStandardLibrary.jl/) ### What's a `structure` dictionary? From 2a553a7479da55a41d44a89b98633a2433351710 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 2 Oct 2023 14:08:07 -0400 Subject: [PATCH 13/16] Be careful with string sorting --- test/odesystem.jl | 7 +++--- test/stream_connectors.jl | 45 ++++++++++++++++++++------------------- 2 files changed, 27 insertions(+), 25 deletions(-) diff --git a/test/odesystem.jl b/test/odesystem.jl index 336eca11e7..39bcda58de 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -23,7 +23,8 @@ eqs = [D(x) ~ σ * (y - x), ModelingToolkit.toexpr.(eqs)[1] @named de = ODESystem(eqs; defaults = Dict(x => 1)) subed = substitute(de, [σ => k]) -@test isequal(sort(parameters(subed), by = string), [k, β, ρ]) +ssort(eqs) = sort(eqs, by = string) +@test isequal(ssort(parameters(subed)), [k, β, ρ]) @test isequal(equations(subed), [D(x) ~ k * (y - x) D(y) ~ (ρ - z) * x - y @@ -348,14 +349,14 @@ eqs = [0 ~ x + z D(accumulation_y) ~ y D(accumulation_z) ~ z D(x) ~ y] -@test sort(equations(asys), by = string) == eqs +@test ssort(equations(asys)) == ssort(eqs) @variables ac(t) asys = add_accumulations(sys, [ac => (x + y)^2]) eqs = [0 ~ x + z 0 ~ x - y D(ac) ~ (x + y)^2 D(x) ~ y] -@test sort(equations(asys), by = string) == eqs +@test ssort(equations(asys)) == ssort(eqs) sys2 = ode_order_lowering(sys) M = ModelingToolkit.calculate_massmatrix(sys2) diff --git a/test/stream_connectors.jl b/test/stream_connectors.jl index 45d58c8552..e4365876cf 100644 --- a/test/stream_connectors.jl +++ b/test/stream_connectors.jl @@ -136,27 +136,28 @@ eqns = [domain_connect(fluid, n1m1.port_a) @test_nowarn structural_simplify(n1m1Test) @unpack source, port_a = n1m1 -@test sort(equations(expand_connections(n1m1)), by = string) == [0 ~ port_a.m_flow +ssort(eqs) = sort(eqs, by = string) +@test ssort(equations(expand_connections(n1m1))) == ssort([0 ~ port_a.m_flow 0 ~ source.port1.m_flow - port_a.m_flow source.port1.P ~ port_a.P source.port1.P ~ source.P source.port1.h_outflow ~ port_a.h_outflow - source.port1.h_outflow ~ source.h] + source.port1.h_outflow ~ source.h]) @unpack port_a, port_b = pipe -@test sort(equations(expand_connections(pipe)), by = string) == - [0 ~ -port_a.m_flow - port_b.m_flow +@test ssort(equations(expand_connections(pipe))) == + ssort([0 ~ -port_a.m_flow - port_b.m_flow 0 ~ port_a.m_flow 0 ~ port_b.m_flow port_a.P ~ port_b.P port_a.h_outflow ~ instream(port_b.h_outflow) - port_b.h_outflow ~ instream(port_a.h_outflow)] -@test sort(equations(expand_connections(sys)), by = string) == - [0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow + port_b.h_outflow ~ instream(port_a.h_outflow)]) +@test ssort(equations(expand_connections(sys))) == + ssort([0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow 0 ~ pipe.port_b.m_flow + sink.port.m_flow n1m1.port_a.P ~ pipe.port_a.P - pipe.port_b.P ~ sink.port.P] -@test sort(equations(expand_connections(n1m1Test)), by = string) == - [0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow + pipe.port_b.P ~ sink.port.P]) +@test ssort(equations(expand_connections(n1m1Test))) == + ssort([0 ~ -pipe.port_a.m_flow - pipe.port_b.m_flow 0 ~ n1m1.source.port1.m_flow - n1m1.port_a.m_flow 0 ~ n1m1.port_a.m_flow + pipe.port_a.m_flow 0 ~ pipe.port_b.m_flow + sink.port.m_flow @@ -172,7 +173,7 @@ eqns = [domain_connect(fluid, n1m1.port_a) pipe.port_b.h_outflow ~ n1m1.port_a.h_outflow sink.port.P ~ sink.P sink.port.h_outflow ~ sink.h_in - sink.port.m_flow ~ -sink.m_flow_in] + sink.port.m_flow ~ -sink.m_flow_in]) # N1M2 model and test code. function N1M2(; name, @@ -255,12 +256,12 @@ eqns = [connect(source.port, n2m2.port_a) @named sp2 = TwoPhaseFluidPort() @named sys = ODESystem([connect(sp1, sp2)], t) sys_exp = expand_connections(compose(sys, [sp1, sp2])) -@test sort(equations(sys_exp), by = string) == [0 ~ -sp1.m_flow - sp2.m_flow +@test ssort(equations(sys_exp)) == ssort([0 ~ -sp1.m_flow - sp2.m_flow 0 ~ sp1.m_flow 0 ~ sp2.m_flow sp1.P ~ sp2.P sp1.h_outflow ~ ModelingToolkit.instream(sp2.h_outflow) - sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)] + sp2.h_outflow ~ ModelingToolkit.instream(sp1.h_outflow)]) # array var @connector function VecPin(; name) @@ -274,15 +275,15 @@ end @named simple = ODESystem([connect(vp1, vp2, vp3)], t) sys = expand_connections(compose(simple, [vp1, vp2, vp3])) -@test sort(equations(sys), by = string) == sort([0 .~ collect(vp1.i) - 0 .~ collect(vp2.i) - 0 .~ collect(vp3.i) - vp1.v[1] ~ vp2.v[1] - vp1.v[2] ~ vp2.v[2] - vp1.v[1] ~ vp3.v[1] - vp1.v[2] ~ vp3.v[2] - 0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1] - 0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]], by = string) +@test ssort(equations(sys)) == ssort([0 .~ collect(vp1.i) + 0 .~ collect(vp2.i) + 0 .~ collect(vp3.i) + vp1.v[1] ~ vp2.v[1] + vp1.v[2] ~ vp2.v[2] + vp1.v[1] ~ vp3.v[1] + vp1.v[2] ~ vp3.v[2] + 0 ~ -vp1.i[1] - vp2.i[1] - vp3.i[1] + 0 ~ -vp1.i[2] - vp2.i[2] - vp3.i[2]]) @connector function VectorHeatPort(; name, N = 100, T0 = 0.0, Q0 = 0.0) @variables (T(t))[1:N]=T0 (Q(t))[1:N]=Q0 [connect = Flow] From f9ea3a6a294cd7ddcd2ef3fda871f7109a6d9626 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Mon, 2 Oct 2023 15:21:45 -0400 Subject: [PATCH 14/16] Update LaTeXify ref tests --- test/latexify/10.tex | 6 +++--- test/latexify/20.tex | 4 ++-- test/latexify/30.tex | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/test/latexify/10.tex b/test/latexify/10.tex index 9a8335bd34..09ec74ab72 100644 --- a/test/latexify/10.tex +++ b/test/latexify/10.tex @@ -1,5 +1,5 @@ \begin{align} -\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\sigma \left( - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left( - y\left( t \right) + x\left( t \right) \right)}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\ -0 =& - y\left( t \right) + \frac{1}{10} \sigma \left( \rho - z\left( t \right) \right) x\left( t \right) \\ -\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - \beta z\left( t \right) +\frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} =& \frac{\left( - x\left( t \right) + y\left( t \right) \right) \frac{\mathrm{d}}{\mathrm{d}t} \left( x\left( t \right) - y\left( t \right) \right) \sigma}{\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t}} \\ +0 =& - y\left( t \right) + \frac{1}{10} x\left( t \right) \left( - z\left( t \right) + \rho \right) \sigma \\ +\frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} =& \left( y\left( t \right) \right)^{\frac{2}{3}} x\left( t \right) - z\left( t \right) \beta \end{align} diff --git a/test/latexify/20.tex b/test/latexify/20.tex index e32ea81eb7..9de213aafd 100644 --- a/test/latexify/20.tex +++ b/test/latexify/20.tex @@ -1,5 +1,5 @@ \begin{align} -\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left( - u(t)_1 + u(t)_2 \right) p_3 \\ -0 =& - u(t)_2 + \frac{1}{10} \left( - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\ +\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left( - u(t)_1 + u(t)_2 \right) \\ +0 =& - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\ \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3 \end{align} diff --git a/test/latexify/30.tex b/test/latexify/30.tex index bc46f41e9b..b83feeba72 100644 --- a/test/latexify/30.tex +++ b/test/latexify/30.tex @@ -1,5 +1,5 @@ \begin{align} -\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& \left( - u(t)_1 + u(t)_2 \right) p_3 \\ -\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =& - u(t)_2 + \frac{1}{10} \left( - u(t)_1 + p_1 \right) p_2 p_3 u(t)_1 \\ +\frac{\mathrm{d}}{\mathrm{d}t} u(t)_1 =& p_3 \left( - u(t)_1 + u(t)_2 \right) \\ +\frac{\mathrm{d}}{\mathrm{d}t} u(t)_2 =& - u(t)_2 + \frac{1}{10} \left( p_1 - u(t)_1 \right) p_2 p_3 u(t)_1 \\ \frac{\mathrm{d}}{\mathrm{d}t} u(t)_3 =& u(t)_2^{\frac{2}{3}} u(t)_1 - p_3 u(t)_3 \end{align} From f2d80b025965983d7d7f5fd65b6ead68ba245ee0 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 4 Oct 2023 15:27:49 -0400 Subject: [PATCH 15/16] 1.6 has a different term ordering --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index e6387e2cd9..b60d82a6ce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,4 +56,6 @@ end @safetestset "FuncAffect Test" include("funcaffect.jl") @safetestset "Constants Test" include("constants.jl") # Reference tests go Last -@safetestset "Latexify recipes Test" include("latexify.jl") +if VERSION >= v"1.9" + @safetestset "Latexify recipes Test" include("latexify.jl") +end From ae050413d35737edc0e7a404708c4159e978aac4 Mon Sep 17 00:00:00 2001 From: Yingbo Ma Date: Wed, 4 Oct 2023 17:00:00 -0400 Subject: [PATCH 16/16] Fix DiscreteProblem construction --- src/systems/jumps/jumpsystem.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/systems/jumps/jumpsystem.jl b/src/systems/jumps/jumpsystem.jl index 8e87d145c2..b37aad5d19 100644 --- a/src/systems/jumps/jumpsystem.jl +++ b/src/systems/jumps/jumpsystem.jl @@ -293,7 +293,7 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap) function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Union{Tuple, Nothing}, parammap = DiffEqBase.NullParameters(); checkbounds = false, - use_union = false, + use_union = true, kwargs...) dvs = states(sys) ps = parameters(sys)