fix: allow passing guesses to linearize, improve linearization docs

docs/src/basics/Linearization.md

@@ -22,7 +22,7 @@ The `linearize` function expects the user to specify the inputs ``u`` and the ou
 ```@example LINEARIZE
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
-@variables x(t)=0 y(t)=0 u(t)=0 r(t)=0
+@variables x(t)=0 y(t) u(t) r(t)=0
 @parameters kp = 1
 
 eqs = [u ~ kp * (r - y) # P controller
@@ -43,7 +43,7 @@ using ModelingToolkit: inputs, outputs
 
 !!! note "Inputs must be unconnected"
 
-    The model above has 4 variables but only three equations, there is no equation specifying the value of `r` since `r` is an input. This means that only unbalanced models can be linearized, or in other words, models that are balanced and can be simulated _cannot_ be linearized. To learn more about this, see [How to linearize a ModelingToolkit model (YouTube)](https://www.youtube.com/watch?v=-XOux-2XDGI&t=395s). Also see [ModelingToolkitStandardLibrary: Linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/) for utilities that make linearization of completed models easier.
+    The model above has 4 variables but only three equations, there is no equation specifying the value of `r` since `r` is an input. This means that only unbalanced models can be linearized, or in other words, models that are balanced and can be simulated _cannot_ be linearized. To learn more about this, see [How to linearize a ModelingToolkit model (YouTube)](https://www.youtube.com/watch?v=-XOux-2XDGI&t=395s). Also see [ModelingToolkitStandardLibrary: Linear analysis](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/linear_analysis/) for utilities that make linearization of completed models easier.
 
 !!! note "Un-simplified system"
 
@@ -57,6 +57,74 @@ The operating point to linearize around can be specified with the keyword argume
 
 If linearization is to be performed around multiple operating points, the simplification of the system has to be carried out a single time only. To facilitate this, the lower-level function [`ModelingToolkit.linearization_function`](@ref) is available. This function further allows you to obtain separate Jacobians for the differential and algebraic parts of the model. For ODE models without algebraic equations, the statespace representation above is available from the output of `linearization_function` as `A, B, C, D = f_x, f_u, h_x, h_u`.
 
+All variables that will be fixed by an operating point _must_ be provided in the operating point to `linearization_function`. For example, if the operating points fix the value of
+`x`, `y` and `z` then an operating point with constant values for these variables (e.g. `Dict(x => 1.0, y => 1.0, z => 1.0)`) must be provided. The constant values themselves
+do not matter and can be changed by subsequent operating points.
+
+One approach to batch linearization would be to call `linearize` in a loop, providing a different operating point each time. For example:
+
+```@example LINEARIZE
+using ModelingToolkitStandardLibrary
+using ModelingToolkitStandardLibrary.Blocks
+
+@parameters k=10 k3=2 c=1
+@variables x(t)=0 [bounds = (-0.5, 1.5)]
+@variables v(t) = 0
+
+@named y = Blocks.RealOutput()
+@named u = Blocks.RealInput()
+
+eqs = [D(x) ~ v
+       D(v) ~ -k * x - k3 * x^3 - c * v + 10u.u
+       y.u ~ x]
+
+@named duffing = ODESystem(eqs, t, systems = [y, u], defaults = [u.u => 0])
+
+# pass a constant value for `x`, since it is the variable we will change in operating points
+linfun, simplified_sys = linearization_function(duffing, [u.u], [y.u]; op = Dict(x => NaN));
+
+println(linearize(simplified_sys, linfun; op = Dict(x => 1.0)))
+println(linearize(simplified_sys, linfun; op = Dict(x => 0.0)))
+
+@time linearize(simplified_sys, linfun; op = Dict(x => 0.0))
+
+nothing # hide
+```
+
+However, this route is still expensive since it has to repeatedly process the symbolic map provided to `op`. `linearize` is simply a wrapper for creating and solving a
+[`ModelingToolkit.LinearizationProblem`](@ref). This object is symbolically indexable, and can thus integrate with SymbolicIndexingInterface.jl for fast updates.
+
+```@example LINEARIZE
+using SymbolicIndexingInterface
+
+# The second argument is the value of the independent variable `t`.
+linprob = LinearizationProblem(linfun, 1.0)
+# It can be mutated
+linprob.t = 0.0
+# create a setter function to update `x` efficiently
+setter! = setu(linprob, x)
+
+function fast_linearize!(problem, setter!, value)
+    setter!(problem, value)
+    solve(problem)
+end
+
+println(fast_linearize!(linprob, setter!, 1.0))
+println(fast_linearize!(linprob, setter!, 0.0))
+
+@time fast_linearize!(linprob, setter!, 1.0)
+
+nothing # hide
+```
+
+Note that `linprob` above can be interacted with similar to a normal `ODEProblem`.
+
+```@repl LINEARIZE
+prob[x]
+prob[x] = 1.5
+prob[x]
+```
+
 ## Symbolic linearization
 
 The function [`ModelingToolkit.linearize_symbolic`](@ref) works similar to [`ModelingToolkit.linearize`](@ref) but returns symbolic rather than numeric Jacobians. Symbolic linearization have several limitations and no all systems that can be linearized numerically can be linearized symbolically.
@@ -75,7 +143,7 @@ If the modeled system is actually proper (but MTK failed to find a proper realiz
 
 ## Tools for linear analysis
 
-[ModelingToolkitStandardLibrary](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/) contains a set of [tools for more advanced linear analysis](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/API/linear_analysis/). These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems.
+ModelingToolkit contains a set of [tools for more advanced linear analysis](https://docs.sciml.ai/ModelingToolkit/stable/tutorials/linear_analysis/). These can be used to make it easier to work with and analyze causal models, such as control and signal-processing systems.
 
 Also see [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK.jl/dev/) for an interface to [ControlSystems.jl](https://github.com/JuliaControl/ControlSystems.jl) that contains tools for linear analysis and frequency-domain analysis.
 
@@ -89,4 +157,5 @@ Pages = ["Linearization.md"]
 linearize
 ModelingToolkit.linearize_symbolic
 ModelingToolkit.linearization_function
+ModelingToolkit.LinearizationProblem
 ```

src/ModelingToolkit.jl

@@ -262,6 +262,7 @@ export independent_variable, equations, controls, observed, full_equations
 export initialization_equations, guesses, defaults, parameter_dependencies, hierarchy
 export structural_simplify, expand_connections, linearize, linearization_function,
        LinearizationProblem
+export solve
 
 export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function
 export calculate_control_jacobian, generate_control_jacobian

src/linearization.jl

@@ -41,6 +41,7 @@ function linearization_function(sys::AbstractSystem, inputs,
         initialization_solver_alg = TrustRegion(),
         eval_expression = false, eval_module = @__MODULE__,
         warn_initialize_determined = true,
+        guesses = Dict(),
         kwargs...)
     op = Dict(op)
     inputs isa AbstractVector || (inputs = [inputs])
@@ -66,11 +67,10 @@ function linearization_function(sys::AbstractSystem, inputs,
         initializealg = initialize ? OverrideInit() : NoInit()
     end
 
-    fun, u0, p = process_SciMLProblem(
-        ODEFunction{true, SciMLBase.FullSpecialize}, sys, op, p;
-        t = 0.0, build_initializeprob = initializealg isa OverrideInit,
-        allow_incomplete = true, algebraic_only = true)
-    prob = ODEProblem(fun, u0, (nothing, nothing), p)
+    prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+        sys, op, (nothing, nothing), p; allow_incomplete = true,
+        algebraic_only = true, guesses)
+    u0 = state_values(prob)
 
     ps = parameters(sys)
     h = build_explicit_observed_function(sys, outputs; eval_expression, eval_module)
@@ -148,6 +148,12 @@ function SymbolicIndexingInterface.parameter_values(f::LinearizationFunction)
 end
 SymbolicIndexingInterface.current_time(f::LinearizationFunction) = current_time(f.prob)
 
+function Base.show(io::IO, mime::MIME"text/plain", lf::LinearizationFunction)
+    printstyled(io, "LinearizationFunction"; bold = true, color = :blue)
+    println(io, " which wraps:")
+    show(io, mime, lf.prob)
+end
+
 """
     $(TYPEDSIGNATURES)
 
@@ -265,6 +271,12 @@ mutable struct LinearizationProblem{F <: LinearizationFunction, T}
     t::T
 end
 
+function Base.show(io::IO, mime::MIME"text/plain", prob::LinearizationProblem)
+    printstyled(io, "LinearizationProblem"; bold = true, color = :blue)
+    println(io, " at time ", prob.t, " which wraps:")
+    show(io, mime, prob.f.prob)
+end
+
 """
     $(TYPEDSIGNATURES)
 

test/downstream/linearize.jl

@@ -317,12 +317,15 @@ matrices = linfun([1.0], Dict(p => 3.0), 1.0)
     @test matrices.f_u[] == 3.0
 end
 
-@testset "Issue #2941" begin
-    @variables x(t) y(t) [guess = 1.0]
+@testset "Issue #2941 and #3400" begin
+    @variables x(t) y(t)
     @parameters p
     eqs = [0 ~ x * log(y) - p]
     @named sys = ODESystem(eqs, t; defaults = [p => 1.0])
     sys = complete(sys)
-    @test_nowarn linearize(
+    @test_throws ModelingToolkit.MissingVariablesError linearize(
         sys, [x], []; op = Dict(x => 1.0), allow_input_derivatives = true)
+    @test_nowarn linearize(
+        sys, [x], []; op = Dict(x => 1.0), guesses = Dict(y => 1.0),
+        allow_input_derivatives = true)
 end