WIP: Add linearize function #1675
Conversation
Co-authored-by: Fredrik Bagge Carlson <baggepinnen@gmail.com>
…nto myb/linearize
|
The test failure appears to be Coveralls having a bad day :/ |
| @test lsys.A == [0 0; 0 -10] | ||
| @test lsys.B == [2 -2; 10 -10] | ||
| @test lsys.C == [400 -4000] | ||
| @test lsys.D == [4400 -4400] |
There was a problem hiding this comment.
Is it possible to get a symbolic version of it?
There was a problem hiding this comment.
WSMLinearize is able to do this.
There was a problem hiding this comment.
You can get a function that calculates the matrices based on an operating point, but the differentiation is done using ForwardDiff so you can not get Symbolic Jacobians. In some cases, it would be possible to provide symbolic jacobians, but right now there is no support for it since symbolic differentiation is less applicable than numeric.
There was a problem hiding this comment.
What do you mean by less applicable?
There was a problem hiding this comment.
There are functions that are not differentiable symbolically (at least not by Symbolics.jl), many of which are differentiable with ForwardDiff. Examples are user-registered functions, stuff with ifelse etc.
There was a problem hiding this comment.
WSM also defaults to numeric diff, probably for the same reason
There was a problem hiding this comment.
There is nothing preventing symbolic jacobians being added in a separate PR, and just error if there is something that can not be differentiated by Symbolics
There was a problem hiding this comment.
By the way, there is calculate_jacobian already that can be used for symbolic linearization of simple systems.
src/systems/abstractsystem.jl
Outdated
| function linearize(sys, lin_fun; op = Dict(), allow_input_derivatives = false) | ||
| x0 = merge(defaults(sys), op) | ||
| prob = ODEProblem(sys, x0, (0.0, 1.0)) | ||
| linres = lin_fun(prob.u0, prob.p, 0.0) |
There was a problem hiding this comment.
Note that in this case, we only use the initial condition and parameters from the ODEProblem. You can use the internal function to compute that instead of creating an ODEProblem each time which would be pretty wasteful.
There was a problem hiding this comment.
hopefully solved by 23c81c5
|
I'm quite happy with the interface by now, and tests pass for simple systems. The Cartesian pendulum below does, however, not work to linearize ## Model with high index
@parameters t L g
@variables x(t)=0 y(t)=0 w(t)=0 z(t)=0 T(t) u(t)=0 [input=true]
D = Differential(t)
# Simple pendulum in cartesian coordinates
eqs = [D(x) ~ w,
D(y) ~ z,
D(w) ~ T * x + u,
D(z) ~ T * y - g,
0 ~ x^2 + y^2 - L^2]
pendulum = ODESystem(eqs, t, name = :pendulum)
lsys, ssys = linearize(pendulum, [u], [x,y])julia> lsys, ssys = linearize(pendulum, [u], [x,y])
ERROR: The LHS operator must be unique. Please run `structural_simplify` or use the DAE form. Differential(t)(y(t)) appears in LHS more than once.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] check_operator_variables(eqs::Vector{Equation}, op::Type{Differential})
@ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/utils.jl:296
[3] linearization_function(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1015
[4] linearization_function
@ ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:994 [inlined]
[5] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num}; op::Dict{Any, Any}, allow_input_derivatives::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1228
[6] linearize(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
@ ModelingToolkit ~/.julia/dev/ModelingToolkit/src/systems/abstractsystem.jl:1226
[7] top-level scope
@ REPL[12]:1 |
baggepinnen
force-pushed
the
myb/linearize
branch
from
23c81c5 to
ddd09dc
Compare
|
Yeah, because it's a DAE. It will be handled later. |
|
Then this PR should be good to go for now |
| function linearize(sys, lin_fun; t = 0.0, op = Dict(), allow_input_derivatives = false, | ||
| p = DiffEqBase.NullParameters()) | ||
| x0 = merge(defaults(sys), op) | ||
| f, u0, p = process_DEProblem(ODEFunction{true}, sys, x0, p) |
There was a problem hiding this comment.
The actual wasteful part is the function construction, since f is exactly lin_fun. Also, beware that ODEProblem/process_DEProblem does NOT necessarily give consistent initial condition.
|
We can merge this for now and add improvements later. |
No description provided.