diff --git a/Project.toml b/Project.toml index ba47d83eb6..8c4a58045e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ModelingToolkit" uuid = "961ee093-0014-501f-94e3-6117800e7a78" authors = ["Yingbo Ma ", "Chris Rackauckas and contributors"] -version = "9.64.1" +version = "9.64.2" [deps] AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c" @@ -65,7 +65,6 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" FMI = "14a09403-18e3-468f-ad8a-74f8dda2d9ac" -HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" InfiniteOpt = "20393b10-9daf-11e9-18c9-8db751c92c57" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" @@ -74,7 +73,6 @@ MTKBifurcationKitExt = "BifurcationKit" MTKChainRulesCoreExt = "ChainRulesCore" MTKDeepDiffsExt = "DeepDiffs" MTKFMIExt = "FMI" -MTKHomotopyContinuationExt = "HomotopyContinuation" MTKInfiniteOptExt = "InfiniteOpt" MTKLabelledArraysExt = "LabelledArrays" @@ -83,6 +81,8 @@ AbstractTrees = "0.3, 0.4" ArrayInterface = "6, 7" BifurcationKit = "0.4" BlockArrays = "1.1" +BoundaryValueDiffEqAscher = "1.1.0" +BoundaryValueDiffEqMIRK = "1.4.0" ChainRulesCore = "1" Combinatorics = "1" CommonSolve = "0.2.4" @@ -110,7 +110,6 @@ ForwardDiff = "0.10.3" FunctionWrappers = "1.1" FunctionWrappersWrappers = "0.1" Graphs = "1.5.2" -HomotopyContinuation = "2.11" InfiniteOpt = "0.5" InteractiveUtils = "1" JuliaFormatter = "1.0.47" @@ -157,6 +156,8 @@ julia = "1.9" [extras] AmplNLWriter = "7c4d4715-977e-5154-bfe0-e096adeac482" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +BoundaryValueDiffEqMIRK = "1a22d4ce-7765-49ea-b6f2-13c8438986a6" +BoundaryValueDiffEqAscher = "7227322d-7511-4e07-9247-ad6ff830280e" ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DeepDiffs = "ab62b9b5-e342-54a8-a765-a90f495de1a6" @@ -189,4 +190,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", "Logging"] +test = ["AmplNLWriter", "BenchmarkTools", "BoundaryValueDiffEqMIRK", "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", "Logging"] \ No newline at end of file diff --git a/ext/MTKChainRulesCoreExt.jl b/ext/MTKChainRulesCoreExt.jl index f153ee77de..9cf17d203d 100644 --- a/ext/MTKChainRulesCoreExt.jl +++ b/ext/MTKChainRulesCoreExt.jl @@ -103,4 +103,6 @@ function ChainRulesCore.rrule( newbuf, pullback end +ChainRulesCore.@non_differentiable Base.getproperty(sys::MTK.AbstractSystem, x::Symbol) + end diff --git a/ext/MTKHomotopyContinuationExt.jl b/ext/MTKHomotopyContinuationExt.jl deleted file mode 100644 index 8f17c05b18..0000000000 --- a/ext/MTKHomotopyContinuationExt.jl +++ /dev/null @@ -1,225 +0,0 @@ -module MTKHomotopyContinuationExt - -using ModelingToolkit -using ModelingToolkit.SciMLBase -using ModelingToolkit.Symbolics: unwrap, symtype, BasicSymbolic, simplify_fractions -using ModelingToolkit.SymbolicIndexingInterface -using ModelingToolkit.DocStringExtensions -using HomotopyContinuation -using ModelingToolkit: iscomplete, parameters, has_index_cache, get_index_cache, get_u0, - get_u0_p, check_eqs_u0, CommonSolve - -const MTK = ModelingToolkit - -""" -$(TYPEDSIGNATURES) - -Convert `expr` from a symbolics expression to one that uses `HomotopyContinuation.ModelKit`. -""" -function symbolics_to_hc(expr) - if iscall(expr) - if operation(expr) == getindex - args = arguments(expr) - return ModelKit.Variable(getname(args[1]), args[2:end]...) - else - return operation(expr)(symbolics_to_hc.(arguments(expr))...) - end - elseif symbolic_type(expr) == NotSymbolic() - return expr - else - return ModelKit.Variable(getname(expr)) - end -end - -""" -$(TYPEDEF) - -A subtype of `HomotopyContinuation.AbstractSystem` used to solve `HomotopyContinuationProblem`s. -""" -struct MTKHomotopySystem{F, P, J, V} <: HomotopyContinuation.AbstractSystem - """ - The generated function for the residual of the polynomial system. In-place. - """ - f::F - """ - The parameter object. - """ - p::P - """ - The generated function for the jacobian of the polynomial system. In-place. - """ - jac::J - """ - The `HomotopyContinuation.ModelKit.Variable` representation of the unknowns of - the system. - """ - vars::V - """ - The number of polynomials in the system. Must also be equal to `length(vars)`. - """ - nexprs::Int -end - -Base.size(sys::MTKHomotopySystem) = (sys.nexprs, length(sys.vars)) -ModelKit.variables(sys::MTKHomotopySystem) = sys.vars - -function (sys::MTKHomotopySystem)(x, p = nothing) - sys.f(x, sys.p) -end - -function ModelKit.evaluate!(u, sys::MTKHomotopySystem, x, p = nothing) - sys.f(u, x, sys.p) -end - -function ModelKit.evaluate_and_jacobian!(u, U, sys::MTKHomotopySystem, x, p = nothing) - sys.f(u, x, sys.p) - sys.jac(U, x, sys.p) -end - -SymbolicIndexingInterface.parameter_values(s::MTKHomotopySystem) = s.p - -""" - $(TYPEDSIGNATURES) - -Create a `HomotopyContinuationProblem` from a `NonlinearSystem` with polynomial equations. -The problem will be solved by HomotopyContinuation.jl. The resultant `NonlinearSolution` -will contain the polynomial root closest to the point specified by `u0map` (if real roots -exist for the system). - -Keyword arguments: -- `eval_expression`: Whether to `eval` the generated functions or use a `RuntimeGeneratedFunction`. -- `eval_module`: The module to use for `eval`/`@RuntimeGeneratedFunction` -- `warn_parametric_exponent`: Whether to warn if the system contains a parametric - exponent preventing the homotopy from being cached. - -All other keyword arguments are forwarded to `HomotopyContinuation.solver_startsystems`. -""" -function MTK.HomotopyContinuationProblem( - sys::NonlinearSystem, u0map, parammap = nothing; kwargs...) - prob = MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap; kwargs...) - prob isa MTK.HomotopyContinuationProblem || throw(prob) - return prob -end - -function MTK._safe_HomotopyContinuationProblem(sys, u0map, parammap = nothing; - fraction_cancel_fn = SymbolicUtils.simplify_fractions, kwargs...) - if !iscomplete(sys) - error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`") - end - transformation = MTK.PolynomialTransformation(sys) - if transformation isa MTK.NotPolynomialError - return transformation - end - result = MTK.transform_system(sys, transformation; fraction_cancel_fn) - if result isa MTK.NotPolynomialError - return result - end - MTK.HomotopyContinuationProblem(sys, transformation, result, u0map, parammap; kwargs...) -end - -function MTK.HomotopyContinuationProblem( - sys::MTK.NonlinearSystem, transformation::MTK.PolynomialTransformation, - result::MTK.PolynomialTransformationResult, u0map, - parammap = nothing; eval_expression = false, - eval_module = ModelingToolkit, warn_parametric_exponent = true, kwargs...) - sys2 = result.sys - denoms = result.denominators - polydata = transformation.polydata - new_dvs = transformation.new_dvs - all_solutions = transformation.all_solutions - - _, u0, p = MTK.process_SciMLProblem( - MTK.EmptySciMLFunction, sys, u0map, parammap; eval_expression, eval_module) - nlfn = NonlinearFunction{true}(sys2; jac = true, eval_expression, eval_module) - - denominator = MTK.build_explicit_observed_function(sys2, denoms) - unpack_solution = MTK.build_explicit_observed_function(sys2, all_solutions) - - hvars = symbolics_to_hc.(new_dvs) - mtkhsys = MTKHomotopySystem(nlfn.f, p, nlfn.jac, hvars, length(new_dvs)) - - obsfn = MTK.ObservedFunctionCache(sys; eval_expression, eval_module) - - has_parametric_exponents = any(d -> d.has_parametric_exponent, polydata) - if has_parametric_exponents - if warn_parametric_exponent - @warn """ - The system has parametric exponents, preventing caching of the homotopy. \ - This will cause `solve` to be slower. Pass `warn_parametric_exponent \ - = false` to turn off this warning - """ - end - solver_and_starts = nothing - else - solver_and_starts = HomotopyContinuation.solver_startsolutions(mtkhsys; kwargs...) - end - return MTK.HomotopyContinuationProblem( - u0, mtkhsys, denominator, sys, obsfn, solver_and_starts, unpack_solution) -end - -""" -$(TYPEDSIGNATURES) - -Solve a `HomotopyContinuationProblem`. Ignores the algorithm passed to it, and always -uses `HomotopyContinuation.jl`. The original solution as returned by -`HomotopyContinuation.jl` will be available in the `.original` field of the returned -`NonlinearSolution`. - -All keyword arguments except the ones listed below are forwarded to -`HomotopyContinuation.solve`. Note that the solver and start solutions are precomputed, -and only keyword arguments related to the solve process are valid. All keyword -arguments have their default values in HomotopyContinuation.jl, except `show_progress` -which defaults to `false`. - -Extra keyword arguments: -- `denominator_abstol`: In case `prob` is solving a rational function, roots which cause - the denominator to be below `denominator_abstol` will be discarded. -""" -function CommonSolve.solve(prob::MTK.HomotopyContinuationProblem, - alg = nothing; show_progress = false, denominator_abstol = 1e-7, kwargs...) - if prob.solver_and_starts === nothing - sol = HomotopyContinuation.solve( - prob.homotopy_continuation_system; show_progress, kwargs...) - else - solver, starts = prob.solver_and_starts - sol = HomotopyContinuation.solve(solver, starts; show_progress, kwargs...) - end - realsols = HomotopyContinuation.results(sol; only_real = true) - if isempty(realsols) - u = state_values(prob) - retcode = SciMLBase.ReturnCode.ConvergenceFailure - resid = prob.homotopy_continuation_system(u) - else - T = eltype(state_values(prob)) - distance = T(Inf) - u = state_values(prob) - resid = nothing - for result in realsols - if any(<=(denominator_abstol), - prob.denominator(real.(result.solution), parameter_values(prob))) - continue - end - for truesol in prob.unpack_solution(result.solution, parameter_values(prob)) - dist = norm(truesol - state_values(prob)) - if dist < distance - distance = dist - u = T.(real.(truesol)) - resid = T.(real.(prob.homotopy_continuation_system(result.solution))) - end - end - end - # all roots cause denominator to be zero - if isinf(distance) - u = state_values(prob) - resid = prob.homotopy_continuation_system(u) - retcode = SciMLBase.ReturnCode.Infeasible - else - retcode = SciMLBase.ReturnCode.Success - end - end - - return SciMLBase.build_solution( - prob, :HomotopyContinuation, u, resid; retcode, original = sol) -end - -end diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index eaa7af85d2..b683e132a2 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -155,6 +155,10 @@ include("systems/codegen_utils.jl") include("systems/problem_utils.jl") include("linearization.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") @@ -170,10 +174,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/structural_transformation/StructuralTransformations.jl b/src/structural_transformation/StructuralTransformations.jl index 1220d517cc..f0124d7f4b 100644 --- a/src/structural_transformation/StructuralTransformations.jl +++ b/src/structural_transformation/StructuralTransformations.jl @@ -11,7 +11,8 @@ using SymbolicUtils: maketerm, iscall using ModelingToolkit using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Differential, - unknowns, equations, vars, Symbolic, diff2term_with_unit, value, + unknowns, equations, vars, Symbolic, diff2term_with_unit, + shift2term_with_unit, value, operation, arguments, Sym, Term, simplify, symbolic_linear_solve, isdiffeq, isdifferential, isirreducible, empty_substitutions, get_substitutions, @@ -22,7 +23,8 @@ using ModelingToolkit: ODESystem, AbstractSystem, var_from_nested_derivative, Di get_postprocess_fbody, vars!, IncrementalCycleTracker, add_edge_checked!, topological_sort, invalidate_cache!, Substitutions, get_or_construct_tearing_state, - filter_kwargs, lower_varname_with_unit, setio, SparseMatrixCLIL, + filter_kwargs, lower_varname_with_unit, + lower_shift_varname_with_unit, setio, SparseMatrixCLIL, get_fullvars, has_equations, observed, Schedule, schedule @@ -63,6 +65,7 @@ export torn_system_jacobian_sparsity export full_equations export but_ordered_incidence, lowest_order_variable_mask, highest_order_variable_mask export computed_highest_diff_variables +export shift2term, lower_shift_varname include("utils.jl") include("pantelides.jl") diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index fe3ce45430..9818bba361 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -237,49 +237,23 @@ function check_diff_graph(var_to_diff, fullvars) end =# -function tearing_reassemble(state::TearingState, var_eq_matching, - full_var_eq_matching = nothing; simplify = false, mm = nothing, cse_hack = true, array_hack = true) - @unpack fullvars, sys, structure = state - @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure - extra_vars = Int[] - if full_var_eq_matching !== nothing - for v in 𝑑vertices(state.structure.graph) - eq = full_var_eq_matching[v] - eq isa Int && continue - push!(extra_vars, v) - end - end +""" +Replace derivatives of non-selected unknown variables by dummy derivatives. - neweqs = collect(equations(state)) - # Terminology and Definition: - # - # A general DAE is in the form of `F(u'(t), u(t), p, t) == 0`. We can - # characterize variables in `u(t)` into two classes: differential variables - # (denoted `v(t)`) and algebraic variables (denoted `z(t)`). Differential - # variables are marked as `SelectedState` and they are differentiated in the - # DAE system, i.e. `v'(t)` are all the variables in `u'(t)` that actually - # appear in the system. Algebraic variables are variables that are not - # differential variables. - # - # Dummy derivatives may determine that some differential variables are - # algebraic variables in disguise. The derivative of such variables are - # called dummy derivatives. - - # Step 1: - # Replace derivatives of non-selected unknown variables by dummy derivatives +State selection may determine that some differential variables are +algebraic variables in disguise. The derivative of such variables are +called dummy derivatives. - if ModelingToolkit.has_iv(state.sys) - iv = get_iv(state.sys) - if is_only_discrete(state.structure) - D = Shift(iv, 1) - else - D = Differential(iv) - end - else - iv = D = nothing - end +`SelectedState` information is no longer needed after this function is called. +State selection is done. All non-differentiated variables are algebraic +variables, and all variables that appear differentiated are differential variables. +""" +function substitute_derivatives_algevars!( + ts::TearingState, neweqs, var_eq_matching, dummy_sub; iv = nothing, D = nothing) + @unpack fullvars, sys, structure = ts + @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure diff_to_var = invview(var_to_diff) - dummy_sub = Dict() + for var in 1:length(fullvars) dv = var_to_diff[var] dv === nothing && continue @@ -310,294 +284,451 @@ function tearing_reassemble(state::TearingState, var_eq_matching, diff_to_var[dv] = nothing end end +end - # `SelectedState` information is no longer needed past here. State selection - # is done. All non-differentiated variables are algebraic variables, and all - # variables that appear differentiated are differential variables. - - ### extract partition information - is_solvable = let solvable_graph = solvable_graph - (eq, iv) -> eq isa Int && iv isa Int && BipartiteEdge(eq, iv) in solvable_graph - end - - # if var is like D(x) - isdervar = let diff_to_var = diff_to_var - var -> diff_to_var[var] !== nothing - end - var_order = let diff_to_var = diff_to_var - dv -> begin - order = 0 - while (dv′ = diff_to_var[dv]) !== nothing - order += 1 - dv = dv′ - end - order, dv - end - end - - #retear = BitSet() - # There are three cases where we want to generate new variables to convert - # the system into first order (semi-implicit) ODEs. - # - # 1. To first order: - # Whenever higher order differentiated variable like `D(D(D(x)))` appears, - # we introduce new variables `x_t`, `x_tt`, and `x_ttt` and new equations - # ``` - # D(x_tt) = x_ttt - # D(x_t) = x_tt - # D(x) = x_t - # ``` - # and replace `D(x)` to `x_t`, `D(D(x))` to `x_tt`, and `D(D(D(x)))` to - # `x_ttt`. - # - # 2. To implicit to semi-implicit ODEs: - # 2.1: Unsolvable derivative: - # If one derivative variable `D(x)` is unsolvable in all the equations it - # appears in, then we introduce a new variable `x_t`, a new equation - # ``` - # D(x) ~ x_t - # ``` - # and replace all other `D(x)` to `x_t`. - # - # 2.2: Solvable derivative: - # If one derivative variable `D(x)` is solvable in at least one of the - # equations it appears in, then we introduce a new variable `x_t`. One of - # the solvable equations must be in the form of `0 ~ L(D(x), u...)` and - # there exists a function `l` such that `D(x) ~ l(u...)`. We should replace - # it to - # ``` - # 0 ~ x_t - l(u...) - # D(x) ~ x_t - # ``` - # and replace all other `D(x)` to `x_t`. - # - # Observe that we don't need to actually introduce a new variable `x_t`, as - # the above equations can be lowered to - # ``` - # x_t := l(u...) - # D(x) ~ x_t - # ``` - # where `:=` denotes assignment. - # - # As a final note, in all the above cases where we need to introduce new - # variables and equations, don't add them when they already exist. +#= +There are three cases where we want to generate new variables to convert +the system into first order (semi-implicit) ODEs. + +1. To first order: +Whenever higher order differentiated variable like `D(D(D(x)))` appears, +we introduce new variables `x_t`, `x_tt`, and `x_ttt` and new equations +``` +D(x_tt) = x_ttt +D(x_t) = x_tt +D(x) = x_t +``` +and replace `D(x)` to `x_t`, `D(D(x))` to `x_tt`, and `D(D(D(x)))` to +`x_ttt`. + +2. To implicit to semi-implicit ODEs: +2.1: Unsolvable derivative: +If one derivative variable `D(x)` is unsolvable in all the equations it +appears in, then we introduce a new variable `x_t`, a new equation +``` +D(x) ~ x_t +``` +and replace all other `D(x)` to `x_t`. + +2.2: Solvable derivative: +If one derivative variable `D(x)` is solvable in at least one of the +equations it appears in, then we introduce a new variable `x_t`. One of +the solvable equations must be in the form of `0 ~ L(D(x), u...)` and +there exists a function `l` such that `D(x) ~ l(u...)`. We should replace +it to +``` +0 ~ x_t - l(u...) +D(x) ~ x_t +``` +and replace all other `D(x)` to `x_t`. + +Observe that we don't need to actually introduce a new variable `x_t`, as +the above equations can be lowered to +``` +x_t := l(u...) +D(x) ~ x_t +``` +where `:=` denotes assignment. + +As a final note, in all the above cases where we need to introduce new +variables and equations, don't add them when they already exist. + +###### DISCRETE SYSTEMS ####### + +Documenting the differences to structural simplification for discrete systems: + +In discrete systems everything gets shifted forward a timestep by `shift_discrete_system` +in order to properly generate the difference equations. + +In the system x(k) ~ x(k-1) + x(k-2), becomes Shift(t, 1)(x(t)) ~ x(t) + Shift(t, -1)(x(t)) + +The lowest-order term is Shift(t, k)(x(t)), instead of x(t). As such we actually want +dummy variables for the k-1 lowest order terms instead of the k-1 highest order terms. + +Shift(t, -1)(x(t)) -> x\_{t-1}(t) + +Since Shift(t, -1)(x) is not a derivative, it is directly substituted in `fullvars`. +No equation or variable is added for it. + +For ODESystems D(D(D(x))) in equations is recursively substituted as D(x) ~ x_t, D(x_t) ~ x_tt, etc. +The analogue for discrete systems, Shift(t, 1)(Shift(t,1)(Shift(t,1)(Shift(t, -3)(x(t))))) +does not actually appear. So `total_sub` in generate_system_equations` is directly +initialized with all of the lowered variables `Shift(t, -3)(x) -> x_t-3(t)`, etc. +=# +""" +Generate new derivative variables for the system. + +Effects on the system structure: +- fullvars: add the new derivative variables x_t +- neweqs: add the identity equations for the new variables, D(x) ~ x_t +- graph: update graph with the new equations and variables, and their connections +- solvable_graph: +- var_eq_matching: match D(x) to the added identity equation D(x) ~ x_t +""" +function generate_derivative_variables!( + ts::TearingState, neweqs, var_eq_matching; mm = nothing, iv = nothing, D = nothing) + @unpack fullvars, sys, structure = ts + @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure eq_var_matching = invview(var_eq_matching) + diff_to_var = invview(var_to_diff) + is_discrete = is_only_discrete(structure) linear_eqs = mm === nothing ? Dict{Int, Int}() : Dict(reverse(en) for en in enumerate(mm.nzrows)) + + # For variable x, make dummy derivative x_t if the + # derivative is in the system for v in 1:length(var_to_diff) dv = var_to_diff[v] dv isa Int || continue solved = var_eq_matching[dv] isa Int solved && continue - # check if there's `D(x) = x_t` already - local v_t, dummy_eq - for eq in 𝑑neighbors(solvable_graph, dv) - mi = get(linear_eqs, eq, 0) - iszero(mi) && continue - row = @view mm[mi, :] - nzs = nonzeros(row) - rvs = SparseArrays.nonzeroinds(row) - # note that `v_t` must not be differentiated - if length(nzs) == 2 && - (abs(nzs[1]) == 1 && nzs[1] == -nzs[2]) && - (v_t = rvs[1] == dv ? rvs[2] : rvs[1]; - diff_to_var[v_t] === nothing) - @assert dv in rvs - dummy_eq = eq - @goto FOUND_DUMMY_EQ - end + + # If there's `D(x) = x_t` already, update mappings and continue without + # adding new equations/variables + dd = find_duplicate_dd(dv, solvable_graph, diff_to_var, linear_eqs, mm) + if !isnothing(dd) + dummy_eq, v_t = dd + var_to_diff[v_t] = var_to_diff[dv] + var_eq_matching[dv] = unassigned + eq_var_matching[dummy_eq] = dv + continue end + dx = fullvars[dv] - # add `x_t` - order, lv = var_order(dv) - x_t = lower_varname_withshift(fullvars[lv], iv, order) - push!(fullvars, simplify_shifts(x_t)) - v_t = length(fullvars) - v_t_idx = add_vertex!(var_to_diff) - add_vertex!(graph, DST) - # TODO: do we care about solvable_graph? We don't use them after - # `dummy_derivative_graph`. - add_vertex!(solvable_graph, DST) - # var_eq_matching is a bit odd. - # length(var_eq_matching) == length(invview(var_eq_matching)) + order, lv = var_order(dv, diff_to_var) + x_t = is_discrete ? lower_shift_varname_with_unit(fullvars[dv], iv) : + lower_varname_with_unit(fullvars[lv], iv, order) + + # Add `x_t` to the graph + v_t = add_dd_variable!(structure, fullvars, x_t, dv) + # Add `D(x) - x_t ~ 0` to the graph + dummy_eq = add_dd_equation!(structure, neweqs, 0 ~ dx - x_t, dv, v_t) + + # Update matching push!(var_eq_matching, unassigned) - @assert v_t_idx == ndsts(graph) == ndsts(solvable_graph) == length(fullvars) == - length(var_eq_matching) - # add `D(x) - x_t ~ 0` - push!(neweqs, 0 ~ x_t - dx) - add_vertex!(graph, SRC) - dummy_eq = length(neweqs) - add_edge!(graph, dummy_eq, dv) - add_edge!(graph, dummy_eq, v_t) - add_vertex!(solvable_graph, SRC) - add_edge!(solvable_graph, dummy_eq, dv) - @assert nsrcs(graph) == nsrcs(solvable_graph) == dummy_eq - @label FOUND_DUMMY_EQ - var_to_diff[v_t] = var_to_diff[dv] var_eq_matching[dv] = unassigned eq_var_matching[dummy_eq] = dv end +end + +""" +Check if there's `D(x) ~ x_t` already. +""" +function find_duplicate_dd(dv, solvable_graph, diff_to_var, linear_eqs, mm) + for eq in 𝑑neighbors(solvable_graph, dv) + mi = get(linear_eqs, eq, 0) + iszero(mi) && continue + row = @view mm[mi, :] + nzs = nonzeros(row) + rvs = SparseArrays.nonzeroinds(row) + # note that `v_t` must not be differentiated + if length(nzs) == 2 && + (abs(nzs[1]) == 1 && nzs[1] == -nzs[2]) && + (v_t = rvs[1] == dv ? rvs[2] : rvs[1]; + diff_to_var[v_t] === nothing) + @assert dv in rvs + return eq, v_t + end + end + return nothing +end + +""" +Add a dummy derivative variable x_t corresponding to symbolic variable D(x) +which has index dv in `fullvars`. Return the new index of x_t. +""" +function add_dd_variable!(s::SystemStructure, fullvars, x_t, dv) + push!(fullvars, simplify_shifts(x_t)) + v_t = length(fullvars) + v_t_idx = add_vertex!(s.var_to_diff) + add_vertex!(s.graph, DST) + # TODO: do we care about solvable_graph? We don't use them after + # `dummy_derivative_graph`. + add_vertex!(s.solvable_graph, DST) + s.var_to_diff[v_t] = s.var_to_diff[dv] + v_t +end + +""" +Add the equation D(x) - x_t ~ 0 to `neweqs`. `dv` and `v_t` are the indices +of the higher-order derivative variable and the newly-introduced dummy +derivative variable. Return the index of the new equation in `neweqs`. +""" +function add_dd_equation!(s::SystemStructure, neweqs, eq, dv, v_t) + push!(neweqs, eq) + add_vertex!(s.graph, SRC) + dummy_eq = length(neweqs) + add_edge!(s.graph, dummy_eq, dv) + add_edge!(s.graph, dummy_eq, v_t) + add_vertex!(s.solvable_graph, SRC) + add_edge!(s.solvable_graph, dummy_eq, dv) + dummy_eq +end + +""" +Solve the equations in `neweqs` to obtain the final equations of the +system. + +For each equation of `neweqs`, do one of the following: + 1. If the equation is solvable for a differentiated variable D(x), + then solve for D(x), and add D(x) ~ sol as a differential equation + of the system. + 2. If the equation is solvable for an un-differentiated variable x, + solve for x and then add x ~ sol as a solved equation. These will + become observables. + 3. If the equation is not solvable, add it as an algebraic equation. + +Solved equations are added to `total_sub`. Occurrences of differential +or solved variables on the RHS of the final equations will get substituted. +The topological sort of the equations ensures that variables are solved for +before they appear in equations. + +Reorder the equations and unknowns to be: + [diffeqs; ...] + [diffvars; ...] +such that the mass matrix is: + [I 0 + 0 0]. + +Order the new equations and variables such that the differential equations +and variables come first. Return the new equations, the solved equations, +the new orderings, and the number of solved variables and equations. +""" +function generate_system_equations!(state::TearingState, neweqs, var_eq_matching; + simplify = false, iv = nothing, D = nothing) + @unpack fullvars, sys, structure = state + @unpack solvable_graph, var_to_diff, eq_to_diff, graph = structure + eq_var_matching = invview(var_eq_matching) + diff_to_var = invview(var_to_diff) + + total_sub = Dict() + if is_only_discrete(structure) + for (i, v) in enumerate(fullvars) + op = operation(v) + op isa Shift && (op.steps < 0) && + begin + lowered = lower_shift_varname_with_unit(v, iv) + total_sub[v] = lowered + fullvars[i] = lowered + end + end + end + + # if var is like D(x) or Shift(t, 1)(x) + isdervar = let diff_to_var = diff_to_var + var -> diff_to_var[var] !== nothing + end + + # Extract partition information + is_solvable = let solvable_graph = solvable_graph + (eq, iv) -> eq isa Int && iv isa Int && BipartiteEdge(eq, iv) in solvable_graph + end - # Will reorder equations and unknowns to be: - # [diffeqs; ...] - # [diffvars; ...] - # such that the mass matrix is: - # [I 0 - # 0 0]. - diffeq_idxs = Int[] - algeeq_idxs = Int[] diff_eqs = Equation[] - alge_eqs = Equation[] + diffeq_idxs = Int[] diff_vars = Int[] - subeqs = Equation[] - solved_equations = Int[] - solved_variables = Int[] - # Solve solvable equations + alge_eqs = Equation[] + algeeq_idxs = Int[] + solved_eqs = Equation[] + solvedeq_idxs = Int[] + solved_vars = Int[] + toporder = topological_sort(DiCMOBiGraph{false}(graph, var_eq_matching)) eqs = Iterators.reverse(toporder) - total_sub = Dict() idep = iv + + # Generate equations. + # Solvable equations of differential variables D(x) become differential equations + # Solvable equations of non-differential variables become observable equations + # Non-solvable equations become algebraic equations. for ieq in eqs iv = eq_var_matching[ieq] - if is_solvable(ieq, iv) - # We don't solve differential equations, but we will need to try to - # convert it into the mass matrix form. - # We cannot solve the differential variable like D(x) - if isdervar(iv) - isnothing(D) && - error("Differential found in a non-differential system. Likely this is a bug in the construction of an initialization system. Please report this issue with a reproducible example. Offending equation: $(equations(sys)[ieq])") - order, lv = var_order(iv) - dx = D(simplify_shifts(lower_varname_withshift( - fullvars[lv], idep, order - 1))) - eq = dx ~ simplify_shifts(Symbolics.fixpoint_sub( - Symbolics.symbolic_linear_solve(neweqs[ieq], - fullvars[iv]), - total_sub; operator = ModelingToolkit.Shift)) - for e in 𝑑neighbors(graph, iv) - e == ieq && continue - for v in 𝑠neighbors(graph, e) - add_edge!(graph, e, v) - end - rem_edge!(graph, e, iv) - end - push!(diff_eqs, eq) - total_sub[simplify_shifts(eq.lhs)] = eq.rhs - push!(diffeq_idxs, ieq) - push!(diff_vars, diff_to_var[iv]) - continue + eq = neweqs[ieq] + + if is_solvable(ieq, iv) && isdervar(iv) + var = fullvars[iv] + isnothing(D) && throw(UnexpectedDifferentialError(equations(sys)[ieq])) + order, lv = var_order(iv, diff_to_var) + dx = D(simplify_shifts(fullvars[lv])) + + neweq = make_differential_equation(var, dx, eq, total_sub) + for e in 𝑑neighbors(graph, iv) + e == ieq && continue + rem_edge!(graph, e, iv) end - eq = neweqs[ieq] + + total_sub[simplify_shifts(neweq.lhs)] = neweq.rhs + push!(diff_eqs, neweq) + push!(diffeq_idxs, ieq) + push!(diff_vars, diff_to_var[iv]) + elseif is_solvable(ieq, iv) var = fullvars[iv] - residual = eq.lhs - eq.rhs - a, b, islinear = linear_expansion(residual, var) - @assert islinear - # 0 ~ a * var + b - # var ~ -b/a - if ModelingToolkit._iszero(a) - @warn "Tearing: solving $eq for $var is singular!" - else - rhs = -b / a - neweq = var ~ Symbolics.fixpoint_sub( - simplify ? - Symbolics.simplify(rhs) : rhs, - total_sub; operator = ModelingToolkit.Shift) - push!(subeqs, neweq) - push!(solved_equations, ieq) - push!(solved_variables, iv) + neweq = make_solved_equation(var, eq, total_sub; simplify) + !isnothing(neweq) && begin + push!(solved_eqs, neweq) + push!(solvedeq_idxs, ieq) + push!(solved_vars, iv) end else - eq = neweqs[ieq] - rhs = eq.rhs - if !(eq.lhs isa Number && eq.lhs == 0) - rhs = eq.rhs - eq.lhs - end - push!(alge_eqs, 0 ~ Symbolics.fixpoint_sub(rhs, total_sub)) + neweq = make_algebraic_equation(eq, total_sub) + push!(alge_eqs, neweq) push!(algeeq_idxs, ieq) end end - # TODO: BLT sorting + + # Generate new equations and orderings neweqs = [diff_eqs; alge_eqs] - inveqsperm = [diffeq_idxs; algeeq_idxs] - eqsperm = zeros(Int, nsrcs(graph)) - for (i, v) in enumerate(inveqsperm) - eqsperm[v] = i - end + eq_ordering = [diffeq_idxs; algeeq_idxs] diff_vars_set = BitSet(diff_vars) if length(diff_vars_set) != length(diff_vars) error("Tearing internal error: lowering DAE into semi-implicit ODE failed!") end - solved_variables_set = BitSet(solved_variables) - invvarsperm = [diff_vars; - setdiff!(setdiff(1:ndsts(graph), diff_vars_set), - solved_variables_set)] + solved_vars_set = BitSet(solved_vars) + var_ordering = [diff_vars; + setdiff!(setdiff(1:ndsts(graph), diff_vars_set), + solved_vars_set)] + + return neweqs, solved_eqs, eq_ordering, var_ordering, length(solved_vars), + length(solved_vars_set) +end + +""" +Occurs when a variable D(x) occurs in a non-differential system. +""" +struct UnexpectedDifferentialError + eq::Equation +end + +function Base.showerror(io::IO, err::UnexpectedDifferentialError) + error("Differential found in a non-differential system. Likely this is a bug in the construction of an initialization system. Please report this issue with a reproducible example. Offending equation: $(err.eq)") +end + +""" +Generate a first-order differential equation whose LHS is `dx`. + +`var` and `dx` represent the same variable, but `var` may be a higher-order differential and `dx` is always first-order. For example, if `var` is D(D(x)), then `dx` would be `D(x_t)`. Solve `eq` for `var`, substitute previously solved variables, and return the differential equation. +""" +function make_differential_equation(var, dx, eq, total_sub) + dx ~ simplify_shifts(Symbolics.fixpoint_sub( + Symbolics.symbolic_linear_solve(eq, var), + total_sub; operator = ModelingToolkit.Shift)) +end + +""" +Generate an algebraic equation. Substitute solved variables into `eq` and return the equation. +""" +function make_algebraic_equation(eq, total_sub) + rhs = eq.rhs + if !(eq.lhs isa Number && eq.lhs == 0) + rhs = eq.rhs - eq.lhs + end + 0 ~ simplify_shifts(Symbolics.fixpoint_sub(rhs, total_sub)) +end + +""" +Solve equation `eq` for `var`, substitute previously solved variables, and return the solved equation. +""" +function make_solved_equation(var, eq, total_sub; simplify = false) + residual = eq.lhs - eq.rhs + a, b, islinear = linear_expansion(residual, var) + @assert islinear + # 0 ~ a * var + b + # var ~ -b/a + if ModelingToolkit._iszero(a) + @warn "Tearing: solving $eq for $var is singular!" + return nothing + else + rhs = -b / a + return var ~ simplify_shifts(Symbolics.fixpoint_sub( + simplify ? + Symbolics.simplify(rhs) : rhs, + total_sub; operator = ModelingToolkit.Shift)) + end +end + +""" +Given the ordering returned by `generate_system_equations!`, update the +tearing state to account for the new order. Permute the variables and equations. +Eliminate the solved variables and equations from the graph and permute the +graph's vertices to account for the new variable/equation ordering. +""" +# TODO: BLT sorting +function reorder_vars!(state::TearingState, var_eq_matching, eq_ordering, + var_ordering, nsolved_eq, nsolved_var) + @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure + + eqsperm = zeros(Int, nsrcs(graph)) + for (i, v) in enumerate(eq_ordering) + eqsperm[v] = i + end varsperm = zeros(Int, ndsts(graph)) - for (i, v) in enumerate(invvarsperm) + for (i, v) in enumerate(var_ordering) varsperm[v] = i end - deps = Vector{Int}[i == 1 ? Int[] : collect(1:(i - 1)) - for i in 1:length(solved_equations)] # Contract the vertices in the structure graph to make the structure match # the new reality of the system we've just created. - graph = contract_variables(graph, var_eq_matching, varsperm, eqsperm, - length(solved_variables), length(solved_variables_set)) + new_graph = contract_variables(graph, var_eq_matching, varsperm, eqsperm, + nsolved_eq, nsolved_var) - # Update system - new_var_to_diff = complete(DiffGraph(length(invvarsperm))) + new_var_to_diff = complete(DiffGraph(length(var_ordering))) for (v, d) in enumerate(var_to_diff) v′ = varsperm[v] (v′ > 0 && d !== nothing) || continue d′ = varsperm[d] new_var_to_diff[v′] = d′ > 0 ? d′ : nothing end - new_eq_to_diff = complete(DiffGraph(length(inveqsperm))) + new_eq_to_diff = complete(DiffGraph(length(eq_ordering))) for (v, d) in enumerate(eq_to_diff) v′ = eqsperm[v] (v′ > 0 && d !== nothing) || continue d′ = eqsperm[d] new_eq_to_diff[v′] = d′ > 0 ? d′ : nothing end + new_fullvars = state.fullvars[var_ordering] + + # Update system structure + @set! state.structure.graph = complete(new_graph) + @set! state.structure.var_to_diff = new_var_to_diff + @set! state.structure.eq_to_diff = new_eq_to_diff + @set! state.fullvars = new_fullvars + state +end - var_to_diff = new_var_to_diff - eq_to_diff = new_eq_to_diff +""" +Update the system equations, unknowns, and observables after simplification. +""" +function update_simplified_system!( + state::TearingState, neweqs, solved_eqs, dummy_sub, var_eq_matching, extra_unknowns; + cse_hack = true, array_hack = true) + @unpack solvable_graph, var_to_diff, eq_to_diff, graph = state.structure diff_to_var = invview(var_to_diff) - old_fullvars = fullvars - @set! state.structure.graph = complete(graph) - @set! state.structure.var_to_diff = var_to_diff - @set! state.structure.eq_to_diff = eq_to_diff - @set! state.fullvars = fullvars = fullvars[invvarsperm] ispresent = let var_to_diff = var_to_diff, graph = graph i -> (!isempty(𝑑neighbors(graph, i)) || (var_to_diff[i] !== nothing && !isempty(𝑑neighbors(graph, var_to_diff[i])))) end sys = state.sys - obs_sub = dummy_sub for eq in neweqs isdiffeq(eq) || continue obs_sub[eq.lhs] = eq.rhs end # TODO: compute the dependency correctly so that we don't have to do this - obs = [fast_substitute(observed(sys), obs_sub); subeqs] + obs = [fast_substitute(observed(sys), obs_sub); solved_eqs] unknowns = Any[v - for (i, v) in enumerate(fullvars) + for (i, v) in enumerate(state.fullvars) if diff_to_var[i] === nothing && ispresent(i)] - if !isempty(extra_vars) - for v in extra_vars - push!(unknowns, old_fullvars[v]) - end - end + unknowns = [unknowns; extra_unknowns] @set! sys.unknowns = unknowns obs, subeqs, deps = cse_and_array_hacks( - sys, obs, subeqs, unknowns, neweqs; cse = cse_hack, array = array_hack) + sys, obs, solved_eqs, unknowns, neweqs; cse = cse_hack, array = array_hack) @set! sys.eqs = neweqs @set! sys.observed = obs - @set! sys.substitutions = Substitutions(subeqs, deps) # Only makes sense for time-dependent @@ -606,6 +737,74 @@ function tearing_reassemble(state::TearingState, var_eq_matching, @set! sys.schedule = Schedule(var_eq_matching, dummy_sub) end sys = schedule(sys) +end + +""" +Give the order of the variable indexed by dv. +""" +function var_order(dv, diff_to_var) + order = 0 + while (dv′ = diff_to_var[dv]) !== nothing + order += 1 + dv = dv′ + end + order, dv +end + +""" +Main internal function for structural simplification for DAE systems and discrete systems. +Generate dummy derivative variables, new equations in terms of variables, return updated +system and tearing state. + +Terminology and Definition: + +A general DAE is in the form of `F(u'(t), u(t), p, t) == 0`. We can +characterize variables in `u(t)` into two classes: differential variables +(denoted `v(t)`) and algebraic variables (denoted `z(t)`). Differential +variables are marked as `SelectedState` and they are differentiated in the +DAE system, i.e. `v'(t)` are all the variables in `u'(t)` that actually +appear in the system. Algebraic variables are variables that are not +differential variables. +""" +function tearing_reassemble(state::TearingState, var_eq_matching, + full_var_eq_matching = nothing; simplify = false, mm = nothing, cse_hack = true, array_hack = true) + extra_vars = Int[] + if full_var_eq_matching !== nothing + for v in 𝑑vertices(state.structure.graph) + eq = full_var_eq_matching[v] + eq isa Int && continue + push!(extra_vars, v) + end + end + extra_unknowns = state.fullvars[extra_vars] + neweqs = collect(equations(state)) + dummy_sub = Dict() + + if ModelingToolkit.has_iv(state.sys) + iv = get_iv(state.sys) + if !is_only_discrete(state.structure) + D = Differential(iv) + else + D = Shift(iv, 1) + end + else + iv = D = nothing + end + + # Structural simplification + substitute_derivatives_algevars!(state, neweqs, var_eq_matching, dummy_sub; iv, D) + + generate_derivative_variables!(state, neweqs, var_eq_matching; mm, iv, D) + + neweqs, solved_eqs, eq_ordering, var_ordering, nelim_eq, nelim_var = generate_system_equations!( + state, neweqs, var_eq_matching; simplify, iv, D) + + state = reorder_vars!( + state, var_eq_matching, eq_ordering, var_ordering, nelim_eq, nelim_var) + + sys = update_simplified_system!(state, neweqs, solved_eqs, dummy_sub, var_eq_matching, + extra_unknowns; cse_hack, array_hack) + @set! state.sys = sys @set! sys.tearing_state = state return invalidate_cache!(sys) diff --git a/src/structural_transformation/utils.jl b/src/structural_transformation/utils.jl index bd24a1d017..f3cea9c7ba 100644 --- a/src/structural_transformation/utils.jl +++ b/src/structural_transformation/utils.jl @@ -449,13 +449,49 @@ end ### Misc ### -function lower_varname_withshift(var, iv, order) - order == 0 && return var - if ModelingToolkit.isoperator(var, ModelingToolkit.Shift) - op = operation(var) - return Shift(op.t, order)(var) +""" +Handle renaming variable names for discrete structural simplification. Three cases: +- positive shift: do nothing +- zero shift: x(t) => Shift(t, 0)(x(t)) +- negative shift: rename the variable +""" +function lower_shift_varname(var, iv) + op = operation(var) + op isa Shift || return Shift(iv, 0)(var, true) # hack to prevent simplification of x(t) - x(t) + if op.steps < 0 + return shift2term(var) + else + return var end - return lower_varname_with_unit(var, iv, order) +end + +""" +Rename a Shift variable with negative shift, Shift(t, k)(x(t)) to xₜ₋ₖ(t). +""" +function shift2term(var) + op = operation(var) + iv = op.t + arg = only(arguments(var)) + is_lowered = !isnothing(ModelingToolkit.getunshifted(arg)) + + backshift = is_lowered ? op.steps + ModelingToolkit.getshift(arg) : op.steps + + num = join(Char(0x2080 + d) for d in reverse!(digits(-backshift))) # subscripted number, e.g. ₁ + ds = join([Char(0x209c), Char(0x208b), num]) + # Char(0x209c) = ₜ + # Char(0x208b) = ₋ (subscripted minus) + + O = is_lowered ? ModelingToolkit.getunshifted(arg) : arg + oldop = operation(O) + newname = backshift != 0 ? Symbol(string(nameof(oldop)), ds) : + Symbol(string(nameof(oldop))) + + newvar = maketerm(typeof(O), Symbolics.rename(oldop, newname), + Symbolics.children(O), Symbolics.metadata(O)) + newvar = setmetadata(newvar, Symbolics.VariableSource, (:variables, newname)) + newvar = setmetadata(newvar, ModelingToolkit.VariableUnshifted, O) + newvar = setmetadata(newvar, ModelingToolkit.VariableShift, backshift) + return newvar end function isdoubleshift(var) @@ -466,6 +502,7 @@ end function simplify_shifts(var) ModelingToolkit.hasshift(var) || return var var isa Equation && return simplify_shifts(var.lhs) ~ simplify_shifts(var.rhs) + (op = operation(var)) isa Shift && op.steps == 0 && return first(arguments(var)) if isdoubleshift(var) op1 = operation(var) vv1 = arguments(var)[1] diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index 0aed6f81e7..35685d5c73 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -823,6 +823,7 @@ for prop in [:eqs :structure :op :constraints + :constraintsystem :controls :loss :bcs diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index 022a0909ed..135bf81729 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -960,3 +960,36 @@ Compute the (linearized) loop-transfer function in analysis point `ap`, from `ap See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref), [`open_loop`](@ref). """ get_looptransfer +# + +""" + generate_control_function(sys::ModelingToolkit.AbstractODESystem, input_ap_name::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}}, dist_ap_name::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}}; system_modifier = identity, kwargs) + +When called with analysis points as input arguments, we assume that all analysis points corresponds to connections that should be opened (broken). The use case for this is to get rid of input signal blocks, such as `Step` or `Sine`, since these are useful for simulation but are not needed when using the plant model in a controller or state estimator. +""" +function generate_control_function( + sys::ModelingToolkit.AbstractODESystem, input_ap_name::Union{ + Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}}, + dist_ap_name::Union{ + Nothing, Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}} = nothing; + system_modifier = identity, + kwargs...) + input_ap_name = canonicalize_ap(sys, input_ap_name) + u = [] + for input_ap in input_ap_name + sys, (du, _) = open_loop(sys, input_ap) + push!(u, du) + end + if dist_ap_name === nothing + return ModelingToolkit.generate_control_function(system_modifier(sys), u; kwargs...) + end + + dist_ap_name = canonicalize_ap(sys, dist_ap_name) + d = [] + for dist_ap in dist_ap_name + sys, (du, _) = open_loop(sys, dist_ap) + push!(d, du) + end + + ModelingToolkit.generate_control_function(system_modifier(sys), u, d; kwargs...) +end diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 9944d92a6e..205d7e4601 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -735,6 +735,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...) @@ -757,6 +763,164 @@ 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] + + fns = generate_function_bc(sys, u0, u0_idxs, tspan) + bc_oop, bc_iip = eval_or_rgf.(fns; eval_expression, eval_module) + bc(sol, p, t) = bc_oop(sol, p, t) + bc(resid, u, p, t) = bc_iip(resid, u, p, t) + + 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) + + Given an ODESystem with constraints, generate the boundary condition function to pass to boundary value problem solvers. + Expression uses the constraints and the provided initial conditions. +""" +function generate_function_bc(sys::ODESystem, u0, u0_idxs, tspan; kwargs...) + iv = get_iv(sys) + sts = unknowns(sys) + ps = parameters(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] + + 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 + 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) + _p = reorder_parameters(sys, ps) + + build_function_wrapper(sys, exprs, sol, _p..., t; output_type = Array, kwargs...) +end + """ ```julia DiffEqBase.DAEProblem{iip}(sys::AbstractODESystem, du0map, u0map, tspan, diff --git a/src/systems/diffeqs/odesystem.jl b/src/systems/diffeqs/odesystem.jl index 1590fd6ebd..33cdc59909 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. @@ -191,7 +193,8 @@ 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, @@ -212,7 +215,8 @@ 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, assertions, metadata, @@ -224,6 +228,7 @@ end function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; controls = Num[], observed = Equation[], + constraintsystem = nothing, systems = ODESystem[], tspan = nothing, name = nothing, @@ -297,9 +302,21 @@ function ODESystem(deqs::AbstractVector{<:Equation}, iv, dvs, ps; if is_dde === nothing is_dde = _check_if_dde(deqs, iv′, systems) end + + if !isempty(systems) && !isnothing(constraintsystem) + conssystems = ConstraintsSystem[] + for sys in systems + cons = get_constraintsystem(sys) + cons !== nothing && push!(conssystems, cons) + end + @show conssystems + @set! constraintsystem.systems = conssystems + end + assertions = Dict{BasicSymbolic, Any}(unwrap(k) => v for (k, v) in assertions) + 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, constraintsystem, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, description, systems, defaults, guesses, nothing, initializesystem, initialization_eqs, schedule, connector_type, preface, cont_callbacks, @@ -307,7 +324,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...) diffvars, allunknowns, ps, eqs = process_equations(eqs, iv) for eq in get(kwargs, :parameter_dependencies, Equation[]) @@ -339,8 +356,22 @@ function ODESystem(eqs, iv; kwargs...) end algevars = setdiff(allunknowns, diffvars) - return ODESystem(eqs, iv, collect(Iterators.flatten((diffvars, algevars))), - collect(new_ps); kwargs...) + consvars = OrderedSet() + constraintsystem = nothing + if !isempty(constraints) + constraintsystem = process_constraint_system(constraints, allunknowns, new_ps, iv) + for st in get_unknowns(constraintsystem) + iscall(st) ? + !in(operation(st)(iv), allunknowns) && push!(consvars, st) : + !in(st, allunknowns) && push!(consvars, st) + end + for p in parameters(constraintsystem) + !in(p, new_ps) && push!(new_ps, p) + end + end + + return ODESystem(eqs, iv, collect(Iterators.flatten((diffvars, algevars, consvars))), + collect(new_ps); constraintsystem, kwargs...) end # NOTE: equality does not check cached Jacobian @@ -436,6 +467,8 @@ an array of inputs `inputs` is given, and `param_only` is false for a time-depen """ function build_explicit_observed_function(sys, ts; inputs = nothing, + disturbance_inputs = nothing, + disturbance_argument = false, expression = false, eval_expression = false, eval_module = @__MODULE__, @@ -512,13 +545,22 @@ function build_explicit_observed_function(sys, ts; ps = setdiff(ps, inputs) # Inputs have been converted to parameters by io_preprocessing, remove those from the parameter list inputs = (inputs,) end + if disturbance_inputs !== nothing + # Disturbance inputs may or may not be included as inputs, depending on disturbance_argument + ps = setdiff(ps, disturbance_inputs) + end + if disturbance_argument + disturbance_inputs = (disturbance_inputs,) + else + disturbance_inputs = () + end ps = reorder_parameters(sys, ps) iv = if is_time_dependent(sys) (get_iv(sys),) else () end - args = (dvs..., inputs..., ps..., iv...) + args = (dvs..., inputs..., ps..., iv..., disturbance_inputs...) p_start = length(dvs) + length(inputs) + 1 p_end = length(dvs) + length(inputs) + length(ps) fns = build_function_wrapper( @@ -668,3 +710,44 @@ 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. +# - 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 process_constraint_system( + constraints::Vector{Equation}, sts, ps, iv; consname = :cons) + isempty(constraints) && return nothing + + constraintsts = OrderedSet() + constraintps = OrderedSet() + + for cons in constraints + collect_vars!(constraintsts, constraintps, cons, iv) + end + + # Validate the states. + 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 = only(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.")) + + isparameter(arg) && push!(constraintps, arg) + 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 + + ConstraintsSystem( + constraints, collect(constraintsts), collect(constraintps); name = consname) +end diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl index 5af8e24c58..8359a3e7de 100644 --- a/src/systems/discrete_system/discrete_system.jl +++ b/src/systems/discrete_system/discrete_system.jl @@ -269,15 +269,22 @@ function shift_u0map_forward(sys::DiscreteSystem, u0map, defs) for k in collect(keys(u0map)) v = u0map[k] if !((op = operation(k)) isa Shift) - error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(k)).") + isnothing(getunshifted(k)) && + error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(k)).") + + updated[Shift(iv, 1)(k)] = v + elseif op.steps > 0 + error("Initial conditions must be for the past state of the unknowns. Instead of providing the condition for $k, provide the condition for $(Shift(iv, -1)(only(arguments(k)))).") + else + updated[Shift(iv, op.steps + 1)(only(arguments(k)))] = v end - updated[Shift(iv, op.steps + 1)(arguments(k)[1])] = v end for var in unknowns(sys) op = operation(var) - op isa Shift || continue - haskey(updated, var) && continue - root = first(arguments(var)) + root = getunshifted(var) + shift = getshift(var) + isnothing(root) && continue + (haskey(updated, Shift(iv, shift)(root)) || haskey(updated, var)) && continue haskey(defs, root) || error("Initial condition for $var not provided.") updated[var] = defs[root] end diff --git a/src/systems/nonlinear/homotopy_continuation.jl b/src/systems/nonlinear/homotopy_continuation.jl index 09e11199e8..474224eacd 100644 --- a/src/systems/nonlinear/homotopy_continuation.jl +++ b/src/systems/nonlinear/homotopy_continuation.jl @@ -1,93 +1,3 @@ -""" -$(TYPEDEF) - -A type of Nonlinear problem which specializes on polynomial systems and uses -HomotopyContinuation.jl to solve the system. Requires importing HomotopyContinuation.jl to -create and solve. -""" -struct HomotopyContinuationProblem{uType, H, D, O, SS, U} <: - SciMLBase.AbstractNonlinearProblem{uType, true} - """ - The initial values of states in the system. If there are multiple real roots of - the system, the one closest to this point is returned. - """ - u0::uType - """ - A subtype of `HomotopyContinuation.AbstractSystem` to solve. Also contains the - parameter object. - """ - homotopy_continuation_system::H - """ - A function with signature `(u, p) -> resid`. In case of rational functions, this - is used to rule out roots of the system which would cause the denominator to be - zero. - """ - denominator::D - """ - The `NonlinearSystem` used to create this problem. Used for symbolic indexing. - """ - sys::NonlinearSystem - """ - A function which generates and returns observed expressions for the given system. - """ - obsfn::O - """ - The HomotopyContinuation.jl solver and start system, obtained through - `HomotopyContinuation.solver_startsystems`. - """ - solver_and_starts::SS - """ - A function which takes a solution of the transformed system, and returns a vector - of solutions for the original system. This is utilized when converting systems - to polynomials. - """ - unpack_solution::U -end - -function HomotopyContinuationProblem(::AbstractSystem, _u0, _p; kwargs...) - error("HomotopyContinuation.jl is required to create and solve `HomotopyContinuationProblem`s. Please run `Pkg.add(\"HomotopyContinuation\")` to continue.") -end - -""" - $(TYPEDSIGNATURES) - -Utility function for `safe_HomotopyContinuationProblem`, implemented in the extension. -""" -function _safe_HomotopyContinuationProblem end - -""" - $(TYPEDSIGNATURES) - -Return a `HomotopyContinuationProblem` if the extension is loaded and the system is -polynomial. If the extension is not loaded, return `nothing`. If the system is not -polynomial, return the appropriate `NotPolynomialError`. -""" -function safe_HomotopyContinuationProblem(sys::NonlinearSystem, args...; kwargs...) - if Base.get_extension(ModelingToolkit, :MTKHomotopyContinuationExt) === nothing - return nothing - end - return _safe_HomotopyContinuationProblem(sys, args...; kwargs...) -end - -SymbolicIndexingInterface.symbolic_container(p::HomotopyContinuationProblem) = p.sys -SymbolicIndexingInterface.state_values(p::HomotopyContinuationProblem) = p.u0 -function SymbolicIndexingInterface.set_state!(p::HomotopyContinuationProblem, args...) - set_state!(p.u0, args...) -end -function SymbolicIndexingInterface.parameter_values(p::HomotopyContinuationProblem) - parameter_values(p.homotopy_continuation_system) -end -function SymbolicIndexingInterface.set_parameter!(p::HomotopyContinuationProblem, args...) - set_parameter!(parameter_values(p), args...) -end -function SymbolicIndexingInterface.observed(p::HomotopyContinuationProblem, sym) - if p.obsfn !== nothing - return p.obsfn(sym) - else - return SymbolicIndexingInterface.observed(p.sys, sym) - end -end - function contains_variable(x, wrt) any(y -> occursin(y, x), wrt) end @@ -562,3 +472,92 @@ function handle_rational_polynomials(x, wrt; fraction_cancel_fn = simplify_fract end return num, den end + +function SciMLBase.HomotopyNonlinearFunction(sys::NonlinearSystem, args...; kwargs...) + ODEFunction{true}(sys, args...; kwargs...) +end + +function SciMLBase.HomotopyNonlinearFunction{true}(sys::NonlinearSystem, args...; + kwargs...) + ODEFunction{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.HomotopyNonlinearFunction{false}(sys::NonlinearSystem, args...; + kwargs...) + ODEFunction{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function SciMLBase.HomotopyNonlinearFunction{iip, specialize}( + sys::NonlinearSystem, args...; eval_expression = false, eval_module = @__MODULE__, + p = nothing, fraction_cancel_fn = SymbolicUtils.simplify_fractions, + kwargs...) where {iip, specialize} + if !iscomplete(sys) + error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationFunction`") + end + transformation = PolynomialTransformation(sys) + if transformation isa NotPolynomialError + throw(transformation) + end + result = transform_system(sys, transformation; fraction_cancel_fn) + if result isa NotPolynomialError + throw(result) + end + + sys2 = result.sys + denoms = result.denominators + polydata = transformation.polydata + new_dvs = transformation.new_dvs + all_solutions = transformation.all_solutions + + # we want to create f, jac etc. according to `sys2` since that will do the solving + # but the `sys` inside for symbolic indexing should be the non-polynomial system + fn = NonlinearFunction{iip}(sys2; eval_expression, eval_module, kwargs...) + obsfn = ObservedFunctionCache( + sys; eval_expression, eval_module, checkbounds = get(kwargs, :checkbounds, false)) + fn = remake(fn; sys = sys, observed = obsfn) + + denominator = build_explicit_observed_function(sys2, denoms) + unpolynomialize = build_explicit_observed_function(sys2, all_solutions) + + inv_mapping = Dict(v => k for (k, v) in transformation.substitution_rules) + polynomialize_terms = [get(inv_mapping, var, var) for var in unknowns(sys2)] + polynomialize = build_explicit_observed_function(sys, polynomialize_terms) + + return HomotopyNonlinearFunction{iip, specialize}( + fn; polynomialize, unpolynomialize, denominator) +end + +struct HomotopyContinuationProblem{iip, specialization} end + +function HomotopyContinuationProblem(sys::NonlinearSystem, args...; kwargs...) + HomotopyContinuationProblem{true}(sys, args...; kwargs...) +end + +function HomotopyContinuationProblem(sys::NonlinearSystem, t, + u0map::StaticArray, + args...; + kwargs...) + HomotopyContinuationProblem{false, SciMLBase.FullSpecialize}( + sys, t, u0map, args...; kwargs...) +end + +function HomotopyContinuationProblem{true}(sys::NonlinearSystem, args...; kwargs...) + HomotopyContinuationProblem{true, SciMLBase.AutoSpecialize}(sys, args...; kwargs...) +end + +function HomotopyContinuationProblem{false}(sys::NonlinearSystem, args...; kwargs...) + HomotopyContinuationProblem{false, SciMLBase.FullSpecialize}(sys, args...; kwargs...) +end + +function HomotopyContinuationProblem{iip, spec}( + sys::NonlinearSystem, u0map, pmap = SciMLBase.NullParameters(); + kwargs...) where {iip, spec} + if !iscomplete(sys) + error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `HomotopyContinuationProblem`") + end + f, u0, p = process_SciMLProblem( + HomotopyNonlinearFunction{iip, spec}, sys, u0map, pmap; kwargs...) + + kwargs = filter_kwargs(kwargs) + return NonlinearProblem{iip}(f, u0, p; kwargs...) +end diff --git a/src/systems/nonlinear/nonlinearsystem.jl b/src/systems/nonlinear/nonlinearsystem.jl index c68fefa7d1..06e843d6be 100644 --- a/src/systems/nonlinear/nonlinearsystem.jl +++ b/src/systems/nonlinear/nonlinearsystem.jl @@ -512,17 +512,10 @@ end function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem, u0map, parammap = DiffEqBase.NullParameters(); - check_length = true, use_homotopy_continuation = false, kwargs...) where {iip} + check_length = true, kwargs...) where {iip} if !iscomplete(sys) error("A completed `NonlinearSystem` is required. Call `complete` or `structural_simplify` on the system before creating a `NonlinearProblem`") end - if use_homotopy_continuation - prob = safe_HomotopyContinuationProblem( - sys, u0map, parammap; check_length, kwargs...) - if prob isa HomotopyContinuationProblem - return prob - end - end f, u0, p = process_SciMLProblem(NonlinearFunction{iip}, sys, u0map, parammap; check_length, kwargs...) pt = something(get_metadata(sys), StandardNonlinearProblem()) @@ -676,23 +669,23 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map, scc_cacheexprs = Dict{TypeT, Vector{Any}}[] scc_eqs = Vector{Equation}[] scc_obs = Vector{Equation}[] + # variables solved in previous SCCs + available_vars = Set() 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) + union!(available_vars, _dvs) + obsidxs = observed_equations_used_by(sys, _eqs; available_vars) # the ones used by previous SCCs can be precomputed into the cache setdiff!(obsidxs, prevobsidxs) _obs = obs[obsidxs] + union!(available_vars, getproperty.(_obs, (:lhs,))) # get all subexpressions in the RHS which we can precompute in the cache # precomputed subexpressions should not contain `banned_vars` banned_vars = Set{Any}(vcat(_dvs, getproperty.(_obs, (:lhs,)))) - filter!(banned_vars) do var - symbolic_type(var) != ArraySymbolic() || - all(j -> var[j] in banned_vars, eachindex(var)) - end state = Dict() for i in eachindex(_obs) _obs[i] = _obs[i].lhs ~ subexpressions_not_involving_vars!( @@ -753,9 +746,12 @@ function SciMLBase.SCCNonlinearProblem{iip}(sys::NonlinearSystem, u0map, _obs = scc_obs[i] cachevars = scc_cachevars[i] cacheexprs = scc_cacheexprs[i] + available_vars = [dvs[reduce(vcat, var_sccs[1:(i - 1)]; init = Int[])]; + getproperty.( + reduce(vcat, scc_obs[1:(i - 1)]; init = []), (:lhs,))] _prevobsidxs = vcat(_prevobsidxs, - observed_equations_used_by(sys, reduce(vcat, values(cacheexprs); init = []))) - + observed_equations_used_by( + sys, reduce(vcat, values(cacheexprs); init = []); available_vars)) if isempty(cachevars) push!(explicitfuns, Returns(nothing)) else diff --git a/src/systems/optimization/constraints_system.jl b/src/systems/optimization/constraints_system.jl index 06795916f8..094640eb37 100644 --- a/src/systems/optimization/constraints_system.jl +++ b/src/systems/optimization/constraints_system.jl @@ -123,7 +123,7 @@ function ConstraintsSystem(constraints, unknowns, ps; name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro")) - cstr = value.(Symbolics.canonical_form.(scalarize(constraints))) + cstr = value.(Symbolics.canonical_form.(vcat(scalarize(constraints)...))) unknowns′ = value.(scalarize(unknowns)) ps′ = value.(ps) diff --git a/src/systems/parameter_buffer.jl b/src/systems/parameter_buffer.jl index f2d0867a32..0954d9a40c 100644 --- a/src/systems/parameter_buffer.jl +++ b/src/systems/parameter_buffer.jl @@ -721,7 +721,7 @@ end function Base.showerror(io::IO, e::MissingParametersError) println(io, MISSING_PARAMETERS_MESSAGE) - println(io, e.vars) + println(io, join(e.vars, ", ")) end function InvalidParameterSizeException(param, val) diff --git a/src/systems/problem_utils.jl b/src/systems/problem_utils.jl index 0be41ea6f7..77f4229696 100644 --- a/src/systems/problem_utils.jl +++ b/src/systems/problem_utils.jl @@ -752,12 +752,14 @@ function process_SciMLProblem( u0Type = typeof(u0map) pType = typeof(pmap) - _u0map = u0map + u0map = to_varmap(u0map, dvs) symbols_to_symbolics!(sys, u0map) - _pmap = pmap pmap = to_varmap(pmap, parameters(sys)) symbols_to_symbolics!(sys, pmap) + + check_inputmap_keys(sys, u0map, pmap) + defs = add_toterms(recursive_unwrap(defaults(sys))) cmap, cs = get_cmap(sys) kwargs = NamedTuple(kwargs) @@ -854,6 +856,42 @@ function process_SciMLProblem( implicit_dae ? (f, du0, u0, p) : (f, u0, p) end +# Check that the keys of a u0map or pmap are valid +# (i.e. are symbolic keys, and are defined for the system.) +function check_inputmap_keys(sys, u0map, pmap) + badvarkeys = Any[] + for k in keys(u0map) + if symbolic_type(k) === NotSymbolic() + push!(badvarkeys, k) + end + end + + badparamkeys = Any[] + for k in keys(pmap) + if symbolic_type(k) === NotSymbolic() + push!(badparamkeys, k) + end + end + (isempty(badvarkeys) && isempty(badparamkeys)) || + throw(InvalidKeyError(collect(badvarkeys), collect(badparamkeys))) +end + +const BAD_KEY_MESSAGE = """ + Undefined keys found in the parameter or initial condition maps. Check if symbolic variable names have been reassigned. + The following keys are invalid: + """ + +struct InvalidKeyError <: Exception + vars::Any + params::Any +end + +function Base.showerror(io::IO, e::InvalidKeyError) + println(io, BAD_KEY_MESSAGE) + println(io, "u0map: $(join(e.vars, ", "))") + println(io, "pmap: $(join(e.params, ", "))") +end + ############## # Legacy functions for backward compatibility ############## diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 1bdc11f06a..0643f32ec4 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -140,16 +140,20 @@ get_fullvars(ts::TransformationState) = ts.fullvars has_equations(::TransformationState) = true Base.@kwdef mutable struct SystemStructure - # Maps the (index of) a variable to the (index of) the variable describing - # its derivative. + """Maps the index of variable x to the index of variable D(x).""" var_to_diff::DiffGraph + """Maps the index of an algebraic equation to the index of the equation it is differentiated into.""" eq_to_diff::DiffGraph # Can be access as # `graph` to automatically look at the bipartite graph # or as `torn` to assert that tearing has run. + """Graph that connects equations to variables that appear in them.""" graph::BipartiteGraph{Int, Nothing} + """Graph that connects equations to the variable they will be solved for during simplification.""" solvable_graph::Union{BipartiteGraph{Int, Nothing}, Nothing} + """Variable types (brownian, variable, parameter) in the system.""" var_types::Union{Vector{VariableType}, Nothing} + """Whether the system is discrete.""" only_discrete::Bool end @@ -197,7 +201,9 @@ function complete!(s::SystemStructure) end mutable struct TearingState{T <: AbstractSystem} <: AbstractTearingState{T} + """The system of equations.""" sys::T + """The set of variables of the system.""" fullvars::Vector structure::SystemStructure extra_eqs::Vector @@ -346,6 +352,8 @@ function TearingState(sys; quick_cancel = false, check = true) eqs[i] = eqs[i].lhs ~ rhs end end + + ### Handle discrete variables lowest_shift = Dict() for var in fullvars if ModelingToolkit.isoperator(var, ModelingToolkit.Shift) @@ -464,9 +472,11 @@ function shift_discrete_system(ts::TearingState) vars!(discvars, eq; op = Union{Sample, Hold}) end iv = get_iv(sys) + discmap = Dict(k => StructuralTransformations.simplify_shifts(Shift(iv, 1)(k)) for k in discvars if any(isequal(k), fullvars) && !isa(operation(k), Union{Sample, Hold})) + for i in eachindex(fullvars) fullvars[i] = StructuralTransformations.simplify_shifts(fast_substitute( fullvars[i], discmap; operator = Union{Sample, Hold})) diff --git a/src/utils.jl b/src/utils.jl index 962801622a..78c665ffca 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -222,7 +222,8 @@ function collect_ivs_from_nested_operator!(ivs, x, target_op) end function iv_from_nested_derivative(x, op = Differential) - if iscall(x) && operation(x) == getindex + if iscall(x) && + (operation(x) == getindex || operation(x) == real || operation(x) == imag) iv_from_nested_derivative(arguments(x)[1], op) elseif iscall(x) operation(x) isa op ? iv_from_nested_derivative(arguments(x)[1], op) : @@ -1028,6 +1029,8 @@ 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) +shift2term_with_unit(x, t) = _with_unit(shift2term, x, t) +lower_shift_varname_with_unit(var, iv) = _with_unit(lower_shift_varname, var, iv, iv) """ $(TYPEDSIGNATURES) @@ -1068,14 +1071,24 @@ Keyword arguments: providing this keyword is not necessary and is only useful to avoid repeatedly calling `vars(exprs)` - `obs`: the list of observed equations. +- `available_vars`: If `exprs` involves a variable `x[1]`, this function will look for + observed equations whose LHS is `x[1]` OR `x`. Sometimes, the latter is not required + since `x[1]` might already be present elsewhere in the generated code (e.g. an argument + to the function) but other elements of `x` are part of the observed equations, thus + requiring them to be obtained from the equation for `x`. Any variable present in + `available_vars` will not be searched for in the observed equations. """ function observed_equations_used_by(sys::AbstractSystem, exprs; - involved_vars = vars(exprs; op = Union{Shift, Differential}), obs = observed(sys)) + involved_vars = vars(exprs; op = Union{Shift, Differential}), obs = observed(sys), available_vars = []) obsvars = getproperty.(obs, :lhs) graph = observed_dependency_graph(obs) + if !(available_vars isa Set) + available_vars = Set(available_vars) + end obsidxs = BitSet() for sym in involved_vars + sym in available_vars && continue arrsym = iscall(sym) && operation(sym) === getindex ? arguments(sym)[1] : nothing idx = findfirst(v -> isequal(v, sym) || isequal(v, arrsym), obsvars) idx === nothing && continue @@ -1204,6 +1217,9 @@ end Find all the unknowns and parameters from the equations of a SDESystem or ODESystem. Return re-ordered equations, differential variables, all variables, and parameters. """ function process_equations(eqs, iv) + if eltype(eqs) <: AbstractVector + eqs = reduce(vcat, eqs) + end eqs = collect(eqs) diffvars = OrderedSet() @@ -1237,8 +1253,8 @@ function process_equations(eqs, iv) throw(ArgumentError("An ODESystem can only have one independent variable.")) diffvar in diffvars && throw(ArgumentError("The differential variable $diffvar is not unique in the system of equations.")) - !(symtype(diffvar) === Real || eltype(symtype(diffvar)) === Real) && - throw(ArgumentError("Differential variable $diffvar has type $(symtype(diffvar)). Differential variables should not be concretely typed.")) + !has_diffvar_type(diffvar) && + throw(ArgumentError("Differential variable $diffvar has type $(symtype(diffvar)). Differential variables should be of a continuous, non-concrete number type: Real, Complex, AbstractFloat, or Number.")) push!(diffvars, diffvar) end push!(diffeq, eq) @@ -1250,6 +1266,13 @@ function process_equations(eqs, iv) diffvars, allunknowns, ps, Equation[diffeq; algeeq; compressed_eqs] end +function has_diffvar_type(diffvar) + st = symtype(diffvar) + st === Real || eltype(st) === Real || st === Complex || eltype(st) === Complex || + st === Number || eltype(st) === Number || st === AbstractFloat || + eltype(st) === AbstractFloat +end + """ $(TYPEDSIGNATURES) diff --git a/src/variables.jl b/src/variables.jl index 83f0681a09..ddc33fa838 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -6,6 +6,9 @@ struct VariableOutput end struct VariableIrreducible end struct VariableStatePriority end struct VariableMisc end +# Metadata for renamed shift variables xₜ₋₁ +struct VariableUnshifted end +struct VariableShift end Symbolics.option_to_metadata_type(::Val{:unit}) = VariableUnit Symbolics.option_to_metadata_type(::Val{:connect}) = VariableConnectType Symbolics.option_to_metadata_type(::Val{:input}) = VariableInput @@ -13,6 +16,8 @@ Symbolics.option_to_metadata_type(::Val{:output}) = VariableOutput Symbolics.option_to_metadata_type(::Val{:irreducible}) = VariableIrreducible Symbolics.option_to_metadata_type(::Val{:state_priority}) = VariableStatePriority Symbolics.option_to_metadata_type(::Val{:misc}) = VariableMisc +Symbolics.option_to_metadata_type(::Val{:unshifted}) = VariableUnshifted +Symbolics.option_to_metadata_type(::Val{:shift}) = VariableShift """ dump_variable_metadata(var) @@ -95,7 +100,7 @@ struct Stream <: AbstractConnectType end # special stream connector Get the connect type of x. See also [`hasconnect`](@ref). """ -getconnect(x) = getconnect(unwrap(x)) +getconnect(x::Num) = getconnect(unwrap(x)) getconnect(x::Symbolic) = Symbolics.getmetadata(x, VariableConnectType, nothing) """ hasconnect(x) @@ -134,7 +139,7 @@ function default_toterm(x) if iscall(x) && (op = operation(x)) isa Operator if !(op isa Differential) if op isa Shift && op.steps < 0 - return x + return shift2term(x) end x = normalize_to_differential(op)(arguments(x)...) end @@ -263,7 +268,7 @@ end end struct IsHistory end -ishistory(x) = ishistory(unwrap(x)) +ishistory(x::Num) = ishistory(unwrap(x)) ishistory(x::Symbolic) = getmetadata(x, IsHistory, false) hist(x, t) = wrap(hist(unwrap(x), t)) function hist(x::Symbolic, t) @@ -575,7 +580,7 @@ end Fetch any miscellaneous data associated with symbolic variable `x`. See also [`hasmisc(x)`](@ref). """ -getmisc(x) = getmisc(unwrap(x)) +getmisc(x::Num) = getmisc(unwrap(x)) getmisc(x::Symbolic) = Symbolics.getmetadata(x, VariableMisc, nothing) """ hasmisc(x) @@ -594,7 +599,7 @@ setmisc(x, miscdata) = setmetadata(x, VariableMisc, miscdata) Fetch the unit associated with variable `x`. This function is a metadata getter for an individual variable, while `get_unit` is used for unit inference on more complicated sdymbolic expressions. """ -getunit(x) = getunit(unwrap(x)) +getunit(x::Num) = getunit(unwrap(x)) getunit(x::Symbolic) = Symbolics.getmetadata(x, VariableUnit, nothing) """ hasunit(x) @@ -602,3 +607,9 @@ getunit(x::Symbolic) = Symbolics.getmetadata(x, VariableUnit, nothing) Check if the variable `x` has a unit. """ hasunit(x) = getunit(x) !== nothing + +getunshifted(x::Num) = getunshifted(unwrap(x)) +getunshifted(x::Symbolic) = Symbolics.getmetadata(x, VariableUnshifted, nothing) + +getshift(x::Num) = getshift(unwrap(x)) +getshift(x::Symbolic) = Symbolics.getmetadata(x, VariableShift, 0) diff --git a/test/bvproblem.jl b/test/bvproblem.jl new file mode 100644 index 0000000000..c5451f681b --- /dev/null +++ b/test/bvproblem.jl @@ -0,0 +1,276 @@ +### TODO: update when BoundaryValueDiffEqAscher is updated to use the normal boundary condition conventions + +using OrdinaryDiffEq +using BoundaryValueDiffEqMIRK, 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) + + 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) + + 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.0, 1.0) + @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 + + # Test with a constraint. + constr = [y(0.5) ~ 2.0] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + + function bc!(resid, u, p, t) + resid[1] = u(0.0)[1] - 1.0 + resid[2] = u(0.5)[2] - 2.0 + end + function bc(u, p, t) + [u(0.0)[1] - 1.0, u(0.5)[2] - 2.0] + end + + u0 = [1.0, 1.0] + tspan = (0.0, 1.0) + p = [1.5, 1.0, 1.0, 3.0] + bvpi1 = SciMLBase.BVProblem(lotkavolterra!, bc!, u0, tspan, p) + bvpi2 = SciMLBase.BVProblem(lotkavolterra, bc, u0, tspan, p) + bvpi3 = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lksys, [x(t) => 1.0], tspan; guesses = [y(t) => 1.0]) + bvpi4 = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}( + lksys, [x(t) => 1.0], tspan; guesses = [y(t) => 1.0]) + + 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 +end + +function test_solvers( + solvers, prob, u0map, constraints, equations = []; dt = 0.05, atol = 1e-2) + 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) + guess = [x(t) => 4.0, y(t) => 2.0] + constr = [x(0.6) ~ 3.5, x(0.3) ~ 7.0] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lksys, u0map, tspan; guesses = guess) + test_solvers(solvers, bvp, u0map, constr; dt = 0.05) + + # Testing that more complicated constraints give correct solutions. + constr = [y(0.2) + x(0.8) ~ 3.0, y(0.3) ~ 2.0] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + bvp = SciMLBase.BVProblem{false, SciMLBase.FullSpecialize}( + lksys, u0map, tspan; guesses = guess) + test_solvers(solvers, bvp, u0map, constr; dt = 0.05) + + constr = [α * β - x(0.6) ~ 0.0, y(0.2) ~ 3.0] + @mtkbuild lksys = ODESystem(eqs, t; constraints = constr) + bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}( + lksys, u0map, tspan; guesses = guess) + 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/discrete_system.jl b/test/discrete_system.jl index 78afafd51d..eea0ffc36b 100644 --- a/test/discrete_system.jl +++ b/test/discrete_system.jl @@ -222,21 +222,6 @@ sol = solve(prob, FunctionMap()) @test reduce(vcat, sol.u) == 1:11 -# test that default values apply to the entire history -@variables x(t) = 1.0 -@mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t) -prob = DiscreteProblem(de, [], (0, 10)) -@test prob[x] == 2.0 -@test prob[x(k - 1)] == 1.0 - -# must provide initial conditions for history -@test_throws ErrorException DiscreteProblem(de, [x => 2.0], (0, 10)) - -# initial values only affect _that timestep_, not the entire history -prob = DiscreteProblem(de, [x(k - 1) => 2.0], (0, 10)) -@test prob[x] == 3.0 -@test prob[x(k - 1)] == 2.0 - # Issue#2585 getdata(buffer, t) = buffer[mod1(Int(t), length(buffer))] @register_symbolic getdata(buffer::Vector, t) @@ -281,3 +266,65 @@ k = ShiftIndex(t) prob = @test_nowarn DiscreteProblem(sys, nothing, (0.0, 1.0)) @test_nowarn solve(prob, FunctionMap()) end + +@testset "Initialization" begin + # test that default values apply to the entire history + @variables x(t) = 1.0 + @mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t) + prob = DiscreteProblem(de, [], (0, 10)) + @test prob[x] == 2.0 + @test prob[x(k - 1)] == 1.0 + + # must provide initial conditions for history + @test_throws ErrorException DiscreteProblem(de, [x => 2.0], (0, 10)) + @test_throws ErrorException DiscreteProblem(de, [x(k + 1) => 2.0], (0, 10)) + + # initial values only affect _that timestep_, not the entire history + prob = DiscreteProblem(de, [x(k - 1) => 2.0], (0, 10)) + @test prob[x] == 3.0 + @test prob[x(k - 1)] == 2.0 + @variables xₜ₋₁(t) + @test prob[xₜ₋₁] == 2.0 + + # Test initial assignment with lowered variable + prob = DiscreteProblem(de, [xₜ₋₁(k - 1) => 4.0], (0, 10)) + @test prob[x(k - 1)] == prob[xₜ₋₁] == 1.0 + @test prob[x] == 5.0 + + # Test missing initial throws error + @variables x(t) + @mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2) * x(k - 3)], t) + @test_throws ErrorException prob=DiscreteProblem(de, [x(k - 3) => 2.0], (0, 10)) + @test_throws ErrorException prob=DiscreteProblem( + de, [x(k - 3) => 2.0, x(k - 1) => 3.0], (0, 10)) + + # Test non-assigned initials are given default value + @variables x(t) = 2.0 + @mtkbuild de = DiscreteSystem([x ~ x(k - 1) + x(k - 2) * x(k - 3)], t) + prob = DiscreteProblem(de, [x(k - 3) => 12.0], (0, 10)) + @test prob[x] == 26.0 + @test prob[x(k - 1)] == 2.0 + @test prob[x(k - 2)] == 2.0 + + # Elaborate test + @variables xₜ₋₂(t) zₜ₋₁(t) z(t) + eqs = [x ~ x(k - 1) + z(k - 2), + z ~ x(k - 2) * x(k - 3) - z(k - 1)^2] + @mtkbuild de = DiscreteSystem(eqs, t) + u0 = [x(k - 1) => 3, + xₜ₋₂(k - 1) => 4, + x(k - 2) => 1, + z(k - 1) => 5, + zₜ₋₁(k - 1) => 12] + prob = DiscreteProblem(de, u0, (0, 10)) + @test prob[x] == 15 + @test prob[z] == -21 + + import ModelingToolkit: shift2term + # unknowns(de) = xₜ₋₁, x, zₜ₋₁, xₜ₋₂, z + vars = ModelingToolkit.value.(unknowns(de)) + @test isequal(shift2term(Shift(t, 1)(vars[1])), vars[2]) + @test isequal(shift2term(Shift(t, 1)(vars[4])), vars[1]) + @test isequal(shift2term(Shift(t, -1)(vars[5])), vars[3]) + @test isequal(shift2term(Shift(t, -2)(vars[2])), vars[4]) +end diff --git a/test/downstream/test_disturbance_model.jl b/test/downstream/test_disturbance_model.jl new file mode 100644 index 0000000000..10fcf9fc1f --- /dev/null +++ b/test/downstream/test_disturbance_model.jl @@ -0,0 +1,215 @@ +#= +This file implements and tests a typical workflow for state estimation with disturbance models +The primary subject of the tests is the analysis-point features and the +analysis-point specific method for `generate_control_function`. +=# +using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Test +using ModelingToolkitStandardLibrary.Mechanical.Rotational +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkit: connect +# using Plots + +using ModelingToolkit: t_nounits as t, D_nounits as D + +indexof(sym, syms) = findfirst(isequal(sym), syms) + +## Build the system model ====================================================== +@mtkmodel SystemModel begin + @parameters begin + m1 = 1 + m2 = 1 + k = 10 # Spring stiffness + c = 3 # Damping coefficient + end + @components begin + inertia1 = Inertia(; J = m1, phi = 0, w = 0) + inertia2 = Inertia(; J = m2, phi = 0, w = 0) + spring = Spring(; c = k) + damper = Damper(; d = c) + torque = Torque(use_support = false) + end + @equations begin + 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) + end +end + +@mtkmodel ModelWithInputs begin + @components begin + input_signal = Blocks.Sine(frequency = 1, amplitude = 1) + disturbance_signal1 = Blocks.Constant(k = 0) + disturbance_signal2 = Blocks.Constant(k = 0) + disturbance_torque1 = Torque(use_support = false) + disturbance_torque2 = Torque(use_support = false) + system_model = SystemModel() + end + @equations begin + connect(input_signal.output, :u, system_model.torque.tau) + connect(disturbance_signal1.output, :d1, disturbance_torque1.tau) + connect(disturbance_signal2.output, :d2, disturbance_torque2.tau) + connect(disturbance_torque1.flange, system_model.inertia1.flange_b) + connect(disturbance_torque2.flange, system_model.inertia2.flange_b) + end +end + +@named model = ModelWithInputs() # Model with load disturbance +ssys = structural_simplify(model) +prob = ODEProblem(ssys, [], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +# plot(sol) + +## +using ControlSystemsBase, ControlSystemsMTK +cmodel = complete(model) +P = cmodel.system_model +lsys = named_ss( + model, [:u, :d1], [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]) + +## +# If we now want to add a disturbance model, we cannot do that since we have already connected a constant to the disturbance input. We have also already used the name `d` for an analysis point, but that might not be an issue since we create an outer model and get a new namespace. + +s = tf("s") +dist(; name) = ODESystem(1 / s; name) + +@mtkmodel SystemModelWithDisturbanceModel begin + @components begin + input_signal = Blocks.Sine(frequency = 1, amplitude = 1) + disturbance_signal1 = Blocks.Constant(k = 0) + disturbance_signal2 = Blocks.Constant(k = 0) + disturbance_torque1 = Torque(use_support = false) + disturbance_torque2 = Torque(use_support = false) + disturbance_model = dist() + system_model = SystemModel() + end + @equations begin + connect(input_signal.output, :u, system_model.torque.tau) + connect(disturbance_signal1.output, :d1, disturbance_model.input) + connect(disturbance_model.output, disturbance_torque1.tau) + connect(disturbance_signal2.output, :d2, disturbance_torque2.tau) + connect(disturbance_torque1.flange, system_model.inertia1.flange_b) + connect(disturbance_torque2.flange, system_model.inertia2.flange_b) + end +end + +@named model_with_disturbance = SystemModelWithDisturbanceModel() +# ssys = structural_simplify(open_loop(model_with_disturbance, :d)) # Open loop worked, but it's a bit awkward that we have to use it here +# lsys2 = named_ss(model_with_disturbance, [:u, :d1], +# [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]) +ssys = structural_simplify(model_with_disturbance) +prob = ODEProblem(ssys, [], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +@test SciMLBase.successful_retcode(sol) +# plot(sol) + +## +# Now we only have an integrating disturbance affecting inertia1, what if we want both integrating and direct Gaussian? We'd need a "PI controller" disturbancemodel. If we add the disturbance model (s+1)/s we get the integrating and non-integrating noises being correlated which is fine, it reduces the dimensions of the sigma point by 1. + +dist3(; name) = ODESystem(ss(1 + 10 / s, balance = false); name) + +@mtkmodel SystemModelWithDisturbanceModel begin + @components begin + input_signal = Blocks.Sine(frequency = 1, amplitude = 1) + disturbance_signal1 = Blocks.Constant(k = 0) + disturbance_signal2 = Blocks.Constant(k = 0) + disturbance_torque1 = Torque(use_support = false) + disturbance_torque2 = Torque(use_support = false) + disturbance_model = dist3() + system_model = SystemModel() + + y = Blocks.Add() + angle_sensor = AngleSensor() + output_disturbance = Blocks.Constant(k = 0) + end + @equations begin + connect(input_signal.output, :u, system_model.torque.tau) + connect(disturbance_signal1.output, :d1, disturbance_model.input) + connect(disturbance_model.output, disturbance_torque1.tau) + connect(disturbance_signal2.output, :d2, disturbance_torque2.tau) + connect(disturbance_torque1.flange, system_model.inertia1.flange_b) + connect(disturbance_torque2.flange, system_model.inertia2.flange_b) + + connect(system_model.inertia1.flange_b, angle_sensor.flange) + connect(angle_sensor.phi, y.input1) + connect(output_disturbance.output, :dy, y.input2) + end +end + +@named model_with_disturbance = SystemModelWithDisturbanceModel() +# ssys = structural_simplify(open_loop(model_with_disturbance, :d)) # Open loop worked, but it's a bit awkward that we have to use it here +# lsys3 = named_ss(model_with_disturbance, [:u, :d1], +# [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]) +ssys = structural_simplify(model_with_disturbance) +prob = ODEProblem(ssys, [], (0.0, 10.0)) +sol = solve(prob, Tsit5()) +@test SciMLBase.successful_retcode(sol) +# plot(sol) + +## Generate function for an augmented Unscented Kalman Filter ===================== +# temp = open_loop(model_with_disturbance, :d) +outputs = [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w] +(f_oop1, f_ip), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function( + model_with_disturbance, [:u], [:d1, :d2, :dy], split = false) + +(f_oop2, f_ip2), x_sym, p_sym, io_sys = ModelingToolkit.generate_control_function( + model_with_disturbance, [:u], [:d1, :d2, :dy], + disturbance_argument = true, split = false) + +measurement = ModelingToolkit.build_explicit_observed_function( + io_sys, outputs, inputs = ModelingToolkit.inputs(io_sys)[1:1]) +measurement2 = ModelingToolkit.build_explicit_observed_function( + io_sys, [io_sys.y.output.u], inputs = ModelingToolkit.inputs(io_sys)[1:1], + disturbance_inputs = ModelingToolkit.inputs(io_sys)[2:end], + disturbance_argument = true) + +op = ModelingToolkit.inputs(io_sys) .=> 0 +x0, p = ModelingToolkit.get_u0_p(io_sys, op, op) +x = zeros(5) +u = zeros(1) +d = zeros(3) +@test f_oop2(x, u, p, t, d) == zeros(5) +@test measurement(x, u, p, 0.0) == [0, 0, 0, 0] +@test measurement2(x, u, p, 0.0, d) == [0] + +# Add to the integrating disturbance input +d = [1, 0, 0] +@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 1, 1] # Affects disturbance state and one velocity +@test measurement2(x, u, p, 0.0, d) == [0] + +d = [0, 1, 0] +@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 0, 1] # Affects one velocity +@test measurement(x, u, p, 0.0) == [0, 0, 0, 0] +@test measurement2(x, u, p, 0.0, d) == [0] + +d = [0, 0, 1] +@test sort(f_oop2(x, u, p, 0.0, d)) == [0, 0, 0, 0, 0] # Affects nothing +@test measurement(x, u, p, 0.0) == [0, 0, 0, 0] +@test measurement2(x, u, p, 0.0, d) == [1] # We have now disturbed the output + +## Further downstream tests that the functions generated above actually have the properties required to use for state estimation +# +# using LowLevelParticleFilters, SeeToDee +# Ts = 0.001 +# discrete_dynamics = SeeToDee.Rk4(f_oop2, Ts) +# nx = length(x_sym) +# nu = 1 +# nw = 2 +# ny = length(outputs) +# R1 = Diagonal([1e-5, 1e-5]) +# R2 = 0.1 * I(ny) +# op = ModelingToolkit.inputs(io_sys) .=> 0 +# x0, p = ModelingToolkit.get_u0_p(io_sys, op, op) +# d0 = LowLevelParticleFilters.SimpleMvNormal(x0, 10.0I(nx)) +# measurement_model = UKFMeasurementModel{Float64, false, false}(measurement, R2; nx, ny) +# kf = UnscentedKalmanFilter{false, false, true, false}( +# discrete_dynamics, measurement_model, R1, d0; nu, Ts, p) + +# tvec = 0:Ts:sol.t[end] +# u = vcat.(Array(sol(tvec, idxs = P.torque.tau.u))) +# y = collect.(eachcol(Array(sol(tvec, idxs = outputs)) .+ 1e-2 .* randn.())) + +# inds = 1:5805 +# res = forward_trajectory(kf, u, y) + +# plot(res, size = (1000, 1000)); +# plot!(sol, idxs = x_sym, sp = (1:nx)', l = :dash); diff --git a/test/extensions/Project.toml b/test/extensions/Project.toml index 5f7afe222a..5b0de73cdf 100644 --- a/test/extensions/Project.toml +++ b/test/extensions/Project.toml @@ -10,6 +10,7 @@ JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a" +NonlinearSolveHomotopyContinuation = "2ac3b008-d579-4536-8c91-a1a5998c2f8b" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" diff --git a/test/extensions/ad.jl b/test/extensions/ad.jl index 0e72b2b7b7..adaf6117c6 100644 --- a/test/extensions/ad.jl +++ b/test/extensions/ad.jl @@ -124,3 +124,13 @@ fwd, back = ChainRulesCore.rrule(remake_buffer, sys, ps, idxs, vals) nsol = solve(nprob, NewtonRaphson()) @test nsol[1] ≈ 10.0 / 1.0 + 9.81 * 1.0 / 2 # anal free fall solution is y = v0*t - g*t^2/2 -> v0 = y/t + g*t/2 end + +@testset "`sys.var` is non-differentiable" begin + @variables x(t) + @mtkbuild sys = ODESystem(D(x) ~ x, t) + prob = ODEProblem(sys, [x => 1.0], (0.0, 1.0)) + + grad = Zygote.gradient(prob) do prob + prob[sys.x] + end +end diff --git a/test/extensions/homotopy_continuation.jl b/test/extensions/homotopy_continuation.jl index 9e15ea857e..ed630ea8e3 100644 --- a/test/extensions/homotopy_continuation.jl +++ b/test/extensions/homotopy_continuation.jl @@ -1,20 +1,32 @@ -using ModelingToolkit, NonlinearSolve, SymbolicIndexingInterface +using ModelingToolkit, NonlinearSolve, NonlinearSolveHomotopyContinuation, + SymbolicIndexingInterface using SymbolicUtils import ModelingToolkit as MTK using LinearAlgebra using Test -@testset "Safe HCProblem" begin - @variables x y z - eqs = [0 ~ x^2 + y^2 + 2x * y - 0 ~ x^2 + 4x + 4 - 0 ~ y * z + 4x^2] - @mtkbuild sys = NonlinearSystem(eqs) - prob = MTK.safe_HomotopyContinuationProblem(sys, [x => 1.0, y => 1.0, z => 1.0], []) - @test prob === nothing +allrootsalg = HomotopyContinuationJL{true}(; threading = false) +singlerootalg = HomotopyContinuationJL{false}(; threading = false) + +function test_single_root(sol; atol = 1e-10) + @test SciMLBase.successful_retcode(sol) + @test norm(sol.resid)≈0.0 atol=atol end -import HomotopyContinuation +function test_all_roots(sol; atol = 1e-10) + @test sol.converged + for nlsol in sol.u + @test SciMLBase.successful_retcode(nlsol) + @test norm(nlsol.resid)≈0.0 atol=1e-10 + end +end + +function solve_allroots_closest(prob) + sol = solve(prob, allrootsalg) + return argmin(sol.u) do nlsol + return norm(nlsol.u - prob.u0) + end +end @testset "No parameters" begin @variables x y z @@ -24,19 +36,13 @@ import HomotopyContinuation @mtkbuild sys = NonlinearSystem(eqs) u0 = [x => 1.0, y => 1.0, z => 1.0] prob = HomotopyContinuationProblem(sys, u0) + @test prob isa NonlinearProblem @test prob[x] == prob[y] == prob[z] == 1.0 @test prob[x + y] == 2.0 - sol = solve(prob; threading = false) - @test SciMLBase.successful_retcode(sol) - @test norm(sol.resid)≈0.0 atol=1e-10 - - prob2 = NonlinearProblem(sys, u0; use_homotopy_continuation = true) - @test prob2 isa HomotopyContinuationProblem - sol = solve(prob2; threading = false) - @test SciMLBase.successful_retcode(sol) - @test norm(sol.resid)≈0.0 atol=1e-10 - - @test NonlinearProblem(sys, u0; use_homotopy_continuation = false) isa NonlinearProblem + sol = solve(prob, singlerootalg) + test_single_root(sol) + sol = solve(prob, allrootsalg) + test_all_roots(sol) end struct Wrapper @@ -61,9 +67,10 @@ end @test prob.ps[q] == 4 @test prob.ps[r].x == [1.0 1.0; 0.0 0.0] @test prob.ps[p * q] == 8.0 - sol = solve(prob; threading = false) - @test SciMLBase.successful_retcode(sol) - @test norm(sol.resid)≈0.0 atol=1e-10 + sol = solve(prob, singlerootalg) + test_single_root(sol) + sol = solve(prob, allrootsalg) + test_all_roots(sol) end @testset "Array variables" begin @@ -79,7 +86,7 @@ end @test prob[x] == 2ones(3) prob.ps[p] = [2, 3, 4] @test prob.ps[p] == [2, 3, 4] - sol = @test_nowarn solve(prob; threading = false) + sol = @test_nowarn solve(prob, singlerootalg) @test sol.retcode == ReturnCode.ConvergenceFailure end @@ -87,11 +94,11 @@ end @variables x = 1.0 @parameters n::Integer = 4 @mtkbuild sys = NonlinearSystem([x^n + x^2 - 1 ~ 0]) - prob = @test_warn ["parametric", "exponent"] HomotopyContinuationProblem(sys, []) - @test prob.solver_and_starts === nothing - @test_nowarn HomotopyContinuationProblem(sys, []; warn_parametric_exponent = false) - sol = solve(prob; threading = false) - @test SciMLBase.successful_retcode(sol) + prob = HomotopyContinuationProblem(sys, []) + sol = solve(prob, singlerootalg) + test_single_root(sol) + sol = solve(prob, allrootsalg) + test_all_roots(sol) end @testset "Polynomial check and warnings" begin @@ -100,45 +107,31 @@ end @test_throws ["Cannot convert", "Unable", "symbolically solve", "Exponent", "not an integer", "not a polynomial"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([x^x - x ~ 0]) @test_throws ["Cannot convert", "Unable", "symbolically solve", "Exponent", "unknowns", "not a polynomial"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([((x^2) / sin(x))^2 + x ~ 0]) @test_throws ["Cannot convert", "both polynomial", "non-polynomial", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @variables y = 2.0 @mtkbuild sys = NonlinearSystem([x^2 + y^2 + 2 ~ 0, y ~ sin(x)]) @test_throws ["Cannot convert", "recognized", "sin", "not a polynomial"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2 ~ 0, sin(x + y) ~ 0]) @test_throws ["Cannot convert", "function of multiple unknowns"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([sin(x)^2 + 1 ~ 0, cos(y) - cos(x) - 1 ~ 0]) @test_throws ["Cannot convert", "multiple non-polynomial terms", "same unknown"] HomotopyContinuationProblem( sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0]) @test_throws ["import Nemo"] HomotopyContinuationProblem(sys, []) - @test MTK.safe_HomotopyContinuationProblem(sys, []) isa MTK.NotPolynomialError - @test NonlinearProblem(sys, []) isa NonlinearProblem end import Nemo @@ -148,7 +141,8 @@ import Nemo @mtkbuild sys = NonlinearSystem([sin(x^2)^2 + sin(x^2) - 1 ~ 0]) prob = HomotopyContinuationProblem(sys, []) @test prob[1] ≈ 2.0 - sol = solve(prob; threading = false) + # singlerootalg doesn't converge + sol = solve(prob, allrootsalg).u[1] _x = sol[1] @test sin(_x^2)^2 + sin(_x^2) - 1≈0.0 atol=1e-12 end @@ -162,9 +156,7 @@ end prob = HomotopyContinuationProblem(sys, []) @test prob[x] ≈ 0.25 @test prob[y] ≈ 0.125 - sol = solve(prob; threading = false) - # can't replicate the solve failure locally, so CI logs might help - @show sol.u sol.original.path_results + sol = solve(prob, allrootsalg).u[1] @test SciMLBase.successful_retcode(sol) @test sol[a]≈0.5 atol=1e-6 @test sol[b]≈0.25 atol=1e-6 @@ -177,12 +169,12 @@ end 0 ~ (x^2 - n * x + n) * (x - 1) / (x - 2) / (x - 3) ]) prob = HomotopyContinuationProblem(sys, []) - sol = solve(prob; threading = false) + sol = solve_allroots_closest(prob) @test sol[x] ≈ 1.0 p = parameter_values(prob) for invalid in [2.0, 3.0] for err in [-9e-8, 0, 9e-8] - @test any(<=(1e-7), prob.denominator([invalid + err, 2.0], p)) + @test any(<=(1e-7), prob.f.denominator([invalid + err, 2.0], p)) end end @@ -195,7 +187,7 @@ end [n]) sys = complete(sys) prob = HomotopyContinuationProblem(sys, []) - sol = solve(prob; threading = false) + sol = solve(prob, singlerootalg) disallowed_x = [4, 5.5] disallowed_y = [7, 5, 4] @test all(!isapprox(sol[x]; atol = 1e-8), disallowed_x) @@ -205,30 +197,30 @@ end p = parameter_values(prob) for val in disallowed_x for err in [-9e-8, 0, 9e-8] - @test any(<=(1e-7), prob.denominator([val + err, 2.0], p)) + @test any(<=(1e-7), prob.f.denominator([val + err, 2.0], p)) end end for val in disallowed_y for err in [-9e-8, 0, 9e-8] - @test any(<=(1e-7), prob.denominator([2.0, val + err], p)) + @test any(<=(1e-7), prob.f.denominator([2.0, val + err], p)) end end - @test prob.denominator([2.0, 4.0], p)[1] <= 1e-8 + @test prob.f.denominator([2.0, 4.0], p)[1] <= 1e-8 @testset "Rational function in observed" begin @variables x=1 y=1 @mtkbuild sys = NonlinearSystem([x^2 + y^2 - 2x - 2 ~ 0, y ~ (x - 1) / (x - 2)]) prob = HomotopyContinuationProblem(sys, []) - @test any(prob.denominator([2.0], parameter_values(prob)) .≈ 0.0) - @test SciMLBase.successful_retcode(solve(prob; threading = false)) + @test any(prob.f.denominator([2.0], parameter_values(prob)) .≈ 0.0) + @test SciMLBase.successful_retcode(solve(prob, singlerootalg)) end @testset "Rational function forced to common denominators" begin @variables x = 1 @mtkbuild sys = NonlinearSystem([0 ~ 1 / (1 + x) - x]) prob = HomotopyContinuationProblem(sys, []) - @test any(prob.denominator([-1.0], parameter_values(prob)) .≈ 0.0) - sol = solve(prob; threading = false) + @test any(prob.f.denominator([-1.0], parameter_values(prob)) .≈ 0.0) + sol = solve(prob, singlerootalg) @test SciMLBase.successful_retcode(sol) @test 1 / (1 + sol.u[1]) - sol.u[1]≈0.0 atol=1e-10 end @@ -238,7 +230,7 @@ end @variables x=1 y @mtkbuild sys = NonlinearSystem([x^2 - 2 ~ 0, y ~ sin(x)]) prob = HomotopyContinuationProblem(sys, []) - sol = @test_nowarn solve(prob; threading = false) + sol = @test_nowarn solve(prob, singlerootalg) @test sol[x] ≈ √2.0 @test sol[y] ≈ sin(√2.0) end @@ -251,10 +243,10 @@ end @testset "`simplify_fractions`" begin prob = HomotopyContinuationProblem(sys, []) - @test prob.denominator([0.0], parameter_values(prob)) ≈ [4.0] + @test prob.f.denominator([0.0], parameter_values(prob)) ≈ [4.0] end @testset "`nothing`" begin prob = HomotopyContinuationProblem(sys, []; fraction_cancel_fn = nothing) - @test sort(prob.denominator([0.0], parameter_values(prob))) ≈ [2.0, 4.0^3] + @test sort(prob.f.denominator([0.0], parameter_values(prob))) ≈ [2.0, 4.0^3] end end diff --git a/test/odesystem.jl b/test/odesystem.jl index de166ef0a1..b9a41e1c3d 100644 --- a/test/odesystem.jl +++ b/test/odesystem.jl @@ -1589,6 +1589,10 @@ end @variables Y(t)[1:3]::String eq = D(Y) ~ [p, p, p] @test_throws ArgumentError @mtkbuild osys = ODESystem([eq], t) + + @variables X(t)::Complex + eq = D(X) ~ p - d * X + @test_nowarn @named osys = ODESystem([eq], t) end # Test `isequal` @@ -1673,3 +1677,50 @@ end prob = ODEProblem{false}(lowered_dae_sys; u0_constructor = x -> SVector(x...)) @test prob.u0 isa SVector 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), z(t)]) + + @parameters t_c + cons = [x(t_c) ~ 3] + @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + @test issetequal(parameters(sys), [a, e, t_c]) + + @parameters g(..) h i + cons = [g(h, i) * x(3) ~ c] + @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + @test issetequal(parameters(sys), [g, h, i, a, e, c]) + + # Test that bad constraints throw errors. + cons = [x(3, 4) ~ 3] # unknowns cannot have multiple args. + @test_throws ArgumentError @mtkbuild sys = ODESystem(eqs, t; constraints = cons) + + cons = [x(y(t)) ~ 2] # unknown arg must be parameter, value, or t + @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_throws ArgumentError @mtkbuild sys = ODESystem(eqs, t; constraints = cons) # Need time argument. + + # Test array variables + @variables x(..)[1:5] + mat = [1 2 0 3 2 + 0 0 3 2 0 + 0 1 3 0 4 + 2 0 0 2 1 + 0 0 2 0 5] + eqs = D(x(t)) ~ mat * x(t) + cons = [x(3) ~ [2, 3, 3, 5, 4]] + @mtkbuild ode = ODESystem(D(x(t)) ~ mat * x(t), t; constraints = cons) + @test length(constraints(ModelingToolkit.get_constraintsystem(ode))) == 5 +end diff --git a/test/problem_validation.jl b/test/problem_validation.jl new file mode 100644 index 0000000000..a0a7afaf3c --- /dev/null +++ b/test/problem_validation.jl @@ -0,0 +1,34 @@ +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D + +@testset "Input map validation" begin + import ModelingToolkit: InvalidKeyError, MissingParametersError + @variables X(t) + @parameters p d + eqs = [D(X) ~ p - d * X] + @mtkbuild osys = ODESystem(eqs, t) + + p = "I accidentally renamed p" + u0 = [X => 1.0] + ps = [p => 1.0, d => 0.5] + @test_throws MissingParametersError oprob=ODEProblem(osys, u0, (0.0, 1.0), ps) + + @parameters p d + ps = [p => 1.0, d => 0.5, "Random stuff" => 3.0] + @test_throws InvalidKeyError oprob=ODEProblem(osys, u0, (0.0, 1.0), ps) + + u0 = [:X => 1.0, "random" => 3.0] + @test_throws InvalidKeyError oprob=ODEProblem(osys, u0, (0.0, 1.0), ps) + + @variables x(t) y(t) z(t) + @parameters a b c d + eqs = [D(x) ~ x * a, D(y) ~ y * c, D(z) ~ b + d] + @mtkbuild sys = ODESystem(eqs, t) + pmap = [a => 1, b => 2, c => 3, d => 4, "b" => 2] + u0map = [x => 1, y => 2, z => 3] + @test_throws InvalidKeyError ODEProblem(sys, u0map, (0.0, 1.0), pmap) + + pmap = [a => 1, b => 2, c => 3, d => 4] + u0map = [x => 1, y => 2, z => 3, :0 => 3] + @test_throws InvalidKeyError ODEProblem(sys, u0map, (0.0, 1.0), pmap) +end diff --git a/test/runtests.jl b/test/runtests.jl index 966b02cacb..0ae66d3755 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -87,6 +87,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") @@ -120,6 +121,7 @@ end @safetestset "Linearization Dummy Derivative Tests" include("downstream/linearization_dd.jl") @safetestset "Inverse Models Test" include("downstream/inversemodel.jl") @safetestset "Analysis Points Test" include("downstream/analysis_points.jl") + @safetestset "Analysis Points Test" include("downstream/test_disturbance_model.jl") end if GROUP == "All" || GROUP == "FMI" diff --git a/test/scc_nonlinear_problem.jl b/test/scc_nonlinear_problem.jl index 57f3d72fb7..ef4638cdb0 100644 --- a/test/scc_nonlinear_problem.jl +++ b/test/scc_nonlinear_problem.jl @@ -253,3 +253,13 @@ import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC sol = solve(prob) @test SciMLBase.successful_retcode(sol) end + +@testset "Array variables split across SCCs" begin + @variables x[1:3] + @parameters (f::Function)(..) + @mtkbuild sys = NonlinearSystem([ + 0 ~ x[1]^2 - 9, x[2] ~ 2x[1], 0 ~ x[3]^2 - x[1]^2 + f(x)]) + prob = SCCNonlinearProblem(sys, [x => ones(3)], [f => sum]) + sol = solve(prob, NewtonRaphson()) + @test SciMLBase.successful_retcode(sol) +end diff --git a/test/structural_transformation/utils.jl b/test/structural_transformation/utils.jl index 863e091aad..4da3d1e924 100644 --- a/test/structural_transformation/utils.jl +++ b/test/structural_transformation/utils.jl @@ -4,6 +4,7 @@ using Graphs using SparseArrays using UnPack using ModelingToolkit: t_nounits as t, D_nounits as D +const ST = StructuralTransformations # Define some variables @parameters L g