output extra information from linearization

Author: baggepinnen

src/linearization.jl

@@ -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
@@ -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)
@@ -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)
 end
 
 """
@@ -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)
@@ -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
 
 """
@@ -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
@@ -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.
@@ -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
 
 """

test/downstream/linearize.jl

@@ -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])