Skip to content

add symbolic linearization function #2119

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Mar 15, 2023
Merged

add symbolic linearization function #2119

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
233a743
add symbolic linearization function
baggepinnen Mar 14, 2023
b27da2e
add docstring
baggepinnen Mar 14, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 37 additions & 3 deletions src/systems/abstractsystem.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1287,6 +1287,35 @@ function linearization_function(sys::AbstractSystem, inputs,
return lin_fun, sys
end

"""
(; A, B, C, D), simplified_sys = linearize_symbolic(sys::AbstractSystem, inputs, outputs; simplify = false, kwargs)

Similar to [`linearize`](@ref), but returns symbolic matrices `A,B,C,D` rather than numeric. While `linearize` uses ForwardDiff to perform the linearization, this function uses `Symbolics.jacobian`.

See [`linearize`](@ref) for a description of the arguments.
"""
function linearize_symbolic(sys::AbstractSystem, inputs,
outputs; simplify = false,
kwargs...)
sys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs; simplify,
kwargs...)
sts = states(sys)
t = get_iv(sys)
p = parameters(sys)

fun = generate_function(sys, sts, p; expression = Val{false})[1]
dx = fun(sts, p, t)

A = Symbolics.jacobian(dx, sts)
B = Symbolics.jacobian(dx, inputs)

h = build_explicit_observed_function(sys, outputs)
y = h(sts, p, t)
C = Symbolics.jacobian(y, sts)
D = Symbolics.jacobian(y, inputs)
(; A, B, C, D), sys
end

function markio!(state, orig_inputs, inputs, outputs; check = true)
fullvars = state.fullvars
inputset = Dict{Any, Bool}(i => false for i in inputs)
Expand Down Expand Up @@ -1343,7 +1372,7 @@ the default values of `sys` are used.

If `allow_input_derivatives = false`, an error will be thrown if input derivatives (``u̇``) appear as inputs in the linearized equations. If input derivatives are allowed, the returned `B` matrix will be of double width, corresponding to the input `[u; u̇]`.

See also [`linearization_function`](@ref) which provides a lower-level interface, and [`ModelingToolkit.reorder_states`](@ref).
See also [`linearization_function`](@ref) which provides a lower-level interface, [`linearize_symbolic`](@ref) and [`ModelingToolkit.reorder_states`](@ref).

See extended help for an example.

Expand Down Expand Up @@ -1405,14 +1434,19 @@ connections = [f.y ~ c.r # filtered reference to controller reference

@named cl = ODESystem(connections, t, systems = [f, c, p])

lsys, ssys = linearize(cl, [f.u], [p.x])
lsys0, ssys = linearize(cl, [f.u], [p.x])
desired_order = [f.x, p.x]
lsys = ModelingToolkit.reorder_states(lsys, states(ssys), desired_order)
lsys = ModelingToolkit.reorder_states(lsys0, states(ssys), desired_order)

@assert lsys.A == [-2 0; 1 -2]
@assert lsys.B == [1; 0;;]
@assert lsys.C == [0 1]
@assert lsys.D[] == 0

## Symbolic linearization
lsys_sym, _ = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])

@assert substitute(lsys_sym.A, ModelingToolkit.defaults(cl)) == lsys.A
```
"""
function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives = false,
Expand Down
29 changes: 25 additions & 4 deletions test/linearize.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,15 +68,22 @@ connections = [f.y ~ c.r # filtered reference to controller reference

@named cl = ODESystem(connections, t, systems = [f, c, p])

lsys, ssys = linearize(cl, [f.u], [p.x])
lsys0, ssys = linearize(cl, [f.u], [p.x])
desired_order = [f.x, p.x]
lsys = ModelingToolkit.reorder_states(lsys, states(ssys), desired_order)
lsys = ModelingToolkit.reorder_states(lsys0, states(ssys), desired_order)

@test lsys.A == [-2 0; 1 -2]
@test lsys.B == reshape([1, 0], 2, 1)
@test lsys.C == [0 1]
@test lsys.D[] == 0

## Symbolic linearization
lsyss, _ = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])

@test substitute(lsyss.A, ModelingToolkit.defaults(cl)) == lsys.A
@test substitute(lsyss.B, ModelingToolkit.defaults(cl)) == lsys.B
@test substitute(lsyss.C, ModelingToolkit.defaults(cl)) == lsys.C
@test substitute(lsyss.D, ModelingToolkit.defaults(cl)) == lsys.D
##
using ModelingToolkitStandardLibrary.Blocks: LimPID
k = 400
Expand All @@ -86,16 +93,24 @@ Nd = 10
@named pid = LimPID(; k, Ti, Td, Nd)

@unpack reference, measurement, ctr_output = pid
lsys, ssys = linearize(pid, [reference.u, measurement.u], [ctr_output.u])
lsys0, ssys = linearize(pid, [reference.u, measurement.u], [ctr_output.u])
@unpack int, der = pid
desired_order = [int.x, der.x]
lsys = ModelingToolkit.reorder_states(lsys, states(ssys), desired_order)
lsys = ModelingToolkit.reorder_states(lsys0, states(ssys), desired_order)

@test lsys.A == [0 0; 0 -10]
@test lsys.B == [2 -2; 10 -10]
@test lsys.C == [400 -4000]
@test lsys.D == [4400 -4400]

lsyss, _ = ModelingToolkit.linearize_symbolic(pid, [reference.u, measurement.u],
[ctr_output.u])

@test substitute(lsyss.A, ModelingToolkit.defaults(pid)) == lsys.A
@test substitute(lsyss.B, ModelingToolkit.defaults(pid)) == lsys.B
@test substitute(lsyss.C, ModelingToolkit.defaults(pid)) == lsys.C
@test substitute(lsyss.D, ModelingToolkit.defaults(pid)) == lsys.D

# Test with the reverse desired state order as well to verify that similarity transform and reoreder_states really works
lsys = ModelingToolkit.reorder_states(lsys, states(ssys), reverse(desired_order))

Expand Down Expand Up @@ -151,6 +166,12 @@ lsys, ssys = linearize(sat, [u], [y])
@test isempty(lsys.C)
@test lsys.D[] == 1

@test_skip lsyss, _ = ModelingToolkit.linearize_symbolic(sat, [u], [y]) # Code gen replaces ifelse with if statements causing symbolic evaluation to fail
# @test substitute(lsyss.A, ModelingToolkit.defaults(sat)) == lsys.A
# @test substitute(lsyss.B, ModelingToolkit.defaults(sat)) == lsys.B
# @test substitute(lsyss.C, ModelingToolkit.defaults(sat)) == lsys.C
# @test substitute(lsyss.D, ModelingToolkit.defaults(sat)) == lsys.D

# outside the linear region the derivative is 0
lsys, ssys = linearize(sat, [u], [y]; op = Dict(u => 2))
@test isempty(lsys.A) # there are no differential variables in this system
Expand Down