From 615c7fc37ff0d76e2c6f066c5e490716a5aecaae Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Feb 2025 16:24:12 +0530 Subject: [PATCH 1/5] fix: allow passing guesses to `linearization_function` --- src/linearization.jl | 10 +++++----- test/downstream/linearize.jl | 9 ++++++--- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/src/linearization.jl b/src/linearization.jl index f83efc1642..01d8ee2f4a 100644 --- a/src/linearization.jl +++ b/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) diff --git a/test/downstream/linearize.jl b/test/downstream/linearize.jl index 3cab641d26..898040fb4f 100644 --- a/test/downstream/linearize.jl +++ b/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 From 23231d6095628ae4199537b1183f56b4b589daaa Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Feb 2025 16:26:35 +0530 Subject: [PATCH 2/5] feat: add `Base.show` methods for `LinearizationFunction` and `LinearizationProblem` --- src/linearization.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/linearization.jl b/src/linearization.jl index 01d8ee2f4a..411620e8f3 100644 --- a/src/linearization.jl +++ b/src/linearization.jl @@ -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) From bbfb34d6678b1f504ef4a5e0102b8412ec7c8766 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Feb 2025 16:26:49 +0530 Subject: [PATCH 3/5] docs: link to MTK's linear analysis doc page rather than the standard library --- docs/src/basics/Linearization.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md index 0b29beec2f..5338deda49 100644 --- a/docs/src/basics/Linearization.md +++ b/docs/src/basics/Linearization.md @@ -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" @@ -75,7 +75,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. From c17ce0a3a2bcbdb79f47c38098f8114a5af2d439 Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Feb 2025 18:03:14 +0530 Subject: [PATCH 4/5] feat: export `CommonSolve.solve` --- src/ModelingToolkit.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index b683e132a2..27d398e633 100644 --- a/src/ModelingToolkit.jl +++ b/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 From 027fe713f3e5f17003f562b24e6338e52d50484d Mon Sep 17 00:00:00 2001 From: Aayush Sabharwal Date: Mon, 24 Feb 2025 18:06:54 +0530 Subject: [PATCH 5/5] docs: add doc example for faster linearization --- docs/src/basics/Linearization.md | 71 +++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md index 5338deda49..1c06ce72d4 100644 --- a/docs/src/basics/Linearization.md +++ b/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 @@ -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. @@ -89,4 +157,5 @@ Pages = ["Linearization.md"] linearize ModelingToolkit.linearize_symbolic ModelingToolkit.linearization_function +ModelingToolkit.LinearizationProblem ```