Skip to content

output extra information from linearization #3429

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

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
+29 −20
Open

output extra information from linearization #3429

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
debd300
output extra information from linearization
baggepinnen Feb 27, 2025
b668014
do not output success
baggepinnen Mar 11, 2025
f24068c
handle extra arg in analysis functions
baggepinnen Mar 12, 2025
6b93b19
test op change
baggepinnen Mar 12, 2025
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
23 changes: 15 additions & 8 deletions src/linearization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,7 @@ end
"""
$(TYPEDSIGNATURES)

Linearize the wrapped system at the point given by `(u, p, t)`.
Linearize the wrapped system at the point given by `(unknowns, p, t)`.
"""
function (linfun::LinearizationFunction)(u, p, t)
if eltype(p) <: Pair
Expand All @@ -182,7 +182,7 @@ function (linfun::LinearizationFunction)(u, p, t)
linfun.prob, integ, fun, linfun.initializealg, Val(true);
linfun.initialize_kwargs...)
if !success
error("Initialization algorithm $(linfun.initializealg) failed with `u = $u` and `p = $p`.")
error("Initialization algorithm $(linfun.initializealg) failed with `unknowns = $u` and `p = $p`.")
end
uf = SciMLBase.UJacobianWrapper(fun, t, p)
fg_xz = ForwardDiff.jacobian(uf, u)
Expand Down Expand Up @@ -211,7 +211,11 @@ function (linfun::LinearizationFunction)(u, p, t)
g_u = fg_u[linfun.alge_idxs, :],
h_x = h_xz[:, linfun.diff_idxs],
h_z = h_xz[:, linfun.alge_idxs],
h_u = h_u)
h_u = h_u,
x = u,
p,
t,
success)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Isn't returning success unnecessary? If !success then we error, so at this point success is guaranteed to be true.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should this be a ReturnCode?

end

"""
Expand Down Expand Up @@ -319,7 +323,7 @@ function CommonSolve.solve(prob::LinearizationProblem; allow_input_derivatives =
p = parameter_values(prob)
t = current_time(prob)
linres = prob.f(u0, p, t)
f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u = linres
f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u, x, p, t, success = linres

nx, nu = size(f_u)
nz = size(f_z, 2)
Expand Down Expand Up @@ -356,7 +360,7 @@ function CommonSolve.solve(prob::LinearizationProblem; allow_input_derivatives =
end
end

(; A, B, C, D)
(; A, B, C, D), (; x, p, t, success)
end

"""
Expand Down Expand Up @@ -487,8 +491,8 @@ function markio!(state, orig_inputs, inputs, outputs; check = true)
end

"""
(; A, B, C, D), simplified_sys = linearize(sys, inputs, outputs; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
(; A, B, C, D) = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)
(; A, B, C, D), simplified_sys, extras = linearize(sys, inputs, outputs; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
(; A, B, C, D), extras = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)

Linearize `sys` between `inputs` and `outputs`, both vectors of variables. Return a NamedTuple with the matrices of a linear statespace representation
on the form
Expand All @@ -510,6 +514,8 @@ If `allow_input_derivatives = false`, an error will be thrown if input derivativ

`zero_dummy_der` can be set to automatically set the operating point to zero for all dummy derivatives.

The return value `extras` is a NamedTuple `(; x, p, t, success)` containing the result of the initialization problem that was solved to determine the operating point.

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

See extended help for an example.
Expand Down Expand Up @@ -616,7 +622,8 @@ function linearize(sys, inputs, outputs; op = Dict(), t = 0.0,
zero_dummy_der,
op,
kwargs...)
linearize(ssys, lin_fun; op, t, allow_input_derivatives), ssys
mats, extras = linearize(ssys, lin_fun; op, t, allow_input_derivatives)
mats, ssys, extras
end

"""
Expand Down
5 changes: 3 additions & 2 deletions test/downstream/linearize.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,14 +14,15 @@ eqs = [u ~ kp * (r - y)

@named sys = ODESystem(eqs, t)

lsys, ssys = linearize(sys, [r], [y])
lsys, ssys, extras = linearize(sys, [r], [y])
lprob = LinearizationProblem(sys, [r], [y])
lsys2 = solve(lprob)
lsys2, extras2 = solve(lprob)

@test lsys.A[] == lsys2.A[] == -2
@test lsys.B[] == lsys2.B[] == 1
@test lsys.C[] == lsys2.C[] == 1
@test lsys.D[] == lsys2.D[] == 0
@test extras == extras2

lsys, ssys = linearize(sys, [r], [r])

Expand Down
Loading