Merge pull request #1675 from SciML/myb/linearize

WIP: Add linearize function

Author: YingboMa

Project.toml

@@ -16,6 +16,7 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
@@ -54,6 +55,7 @@ DiffRules = "0.1, 1.0"
 Distributions = "0.23, 0.24, 0.25"
 DocStringExtensions = "0.7, 0.8, 0.9"
 DomainSets = "0.5"
+ForwardDiff = "0.10.3"
 Graphs = "1.5.2"
 IfElse = "0.1"
 JuliaFormatter = "1"
@@ -82,6 +84,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Ipopt_jll = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
+ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
@@ -95,4 +98,4 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]
+test = ["AmplNLWriter", "BenchmarkTools", "ForwardDiff", "Ipopt", "Ipopt_jll", "ModelingToolkitStandardLibrary", "Optimization", "OptimizationOptimJL", "OptimizationMOI", "OrdinaryDiffEq", "Random", "ReferenceTests", "SafeTestsets", "SteadyStateDiffEq", "Test", "StochasticDiffEq", "Sundials"]

src/ModelingToolkit.jl

@@ -4,7 +4,7 @@ $(DocStringExtensions.README)
 module ModelingToolkit
 using DocStringExtensions
 using AbstractTrees
-using DiffEqBase, SciMLBase, Reexport
+using DiffEqBase, SciMLBase, ForwardDiff, Reexport
 using Distributed
 using StaticArrays, LinearAlgebra, SparseArrays, LabelledArrays
 using InteractiveUtils
@@ -181,7 +181,7 @@ export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
 export independent_variables, independent_variable, states, parameters, equations, controls,
        observed, structure, full_equations
-export structural_simplify, expand_connections
+export structural_simplify, expand_connections, linearize, linear_statespace
 export DiscreteSystem, DiscreteProblem
 
 export calculate_jacobian, generate_jacobian, generate_function

src/inputoutput.jl

@@ -61,20 +61,20 @@ unbound_outputs(sys) = filter(x -> !is_bound(sys, x), outputs(sys))
 Determine whether or not input/output variable `u` is "bound" within the system, i.e., if it's to be considered internal to `sys`.
 A variable/signal is considered bound if it appears in an equation together with variables from other subsystems.
 The typical usecase for this function is to determine whether the input to an IO component is connected to another component,
-or if it remains an external input that the user has to supply before simulating the system. 
+or if it remains an external input that the user has to supply before simulating the system.
 
 See also [`bound_inputs`](@ref), [`unbound_inputs`](@ref), [`bound_outputs`](@ref), [`unbound_outputs`](@ref)
 """
 function is_bound(sys, u, stack = [])
     #=
-    For observed quantities, we check if a variable is connected to something that is bound to something further out. 
+    For observed quantities, we check if a variable is connected to something that is bound to something further out.
     In the following scenario
     julia> observed(syss)
         2-element Vector{Equation}:
         sys₊y(tv) ~ sys₊x(tv)
         y(tv) ~ sys₊x(tv)
     sys₊y(t) is bound to the outer y(t) through the variable sys₊x(t) and should thus return is_bound(sys₊y(t)) = true.
-    When asking is_bound(sys₊y(t)), we know that we are looking through observed equations and can thus ask 
+    When asking is_bound(sys₊y(t)), we know that we are looking through observed equations and can thus ask
     if var is bound, if it is, then sys₊y(t) is also bound. This can lead to an infinite recursion, so we maintain a stack of variables we have previously asked about to be able to break cycles
     =#
     u ∈ Set(stack) && return false # Cycle detected
@@ -241,7 +241,7 @@ function toparam(sys, ctrls::AbstractVector)
     ODESystem(eqs, name = nameof(sys))
 end
 
-function inputs_to_parameters!(state::TransformationState)
+function inputs_to_parameters!(state::TransformationState, check_bound = true)
     @unpack structure, fullvars, sys = state
     @unpack var_to_diff, graph, solvable_graph = structure
     @assert solvable_graph === nothing
@@ -254,7 +254,7 @@ function inputs_to_parameters!(state::TransformationState)
     input_to_parameters = Dict()
     new_fullvars = []
     for (i, v) in enumerate(fullvars)
-        if isinput(v) && !is_bound(sys, v)
+        if isinput(v) && !(check_bound && is_bound(sys, v))
             if var_to_diff[i] !== nothing
                 error("Input $(fullvars[i]) is differentiated!")
             end
@@ -270,7 +270,7 @@ function inputs_to_parameters!(state::TransformationState)
             push!(new_fullvars, v)
         end
     end
-    ninputs == 0 && return state
+    ninputs == 0 && return (state, 1:0)
 
     nvars = ndsts(graph) - ninputs
     new_graph = BipartiteGraph(nsrcs(graph), nvars, Val(false))
@@ -296,9 +296,12 @@ function inputs_to_parameters!(state::TransformationState)
 
     @set! sys.eqs = map(Base.Fix2(substitute, input_to_parameters), equations(sys))
     @set! sys.states = setdiff(states(sys), keys(input_to_parameters))
-    @set! sys.ps = [parameters(sys); new_parameters]
+    ps = parameters(sys)
+    @set! sys.ps = [ps; new_parameters]
 
     @set! state.sys = sys
     @set! state.fullvars = new_fullvars
     @set! state.structure = structure
+    base_params = length(ps)
+    return state, (base_params + 1):(base_params + length(new_parameters)) # (1:length(new_parameters)) .+ base_params
 end