I’m working with a DAE system in ModelingToolkit.jl, and I’d like to run only the algebraic part of the integrator-that is, solve the algebraic constraints for given differential states, without advancing the differential equations part.
My main motivation is to use FiniteDiff.jl to compute the Jacobian with respect to the differential states only, and not with respect to the algebraic variables.
Is this possible, and how should I do this?
Hi Bart,
Would this be the consistent initialization of the algebraic states for different initial conditions of the differential states? If so, you could build the initialization problem using MTK for some ODESystem in semi-explicit mass matrix form. Then you grab the system of equations from the initialization problem and specify the finite_difference_jacobian on only the differential states. Let me know if there is any lack of clarity on grabbing the system of equations for initialization.
Edit: I see now that you have specified that you are using a DAE system. It may be beneficial to use the ODE System and then convert the problem to the implicit form after working with structural simplify (if you are keen on still using a solver like IDA).
Would you mind sharing what you need this for?
A linearized kite model where the quasi-static tether points (tether_point_acc ~ 0) are not part of this linearized model. So pulling on the tether moves the kite instantly in the linearized model.
This function returns the individual jacobians for differential and algebraic variables
The autodiff kwarg allows you to select FiniteDiff as differentiation backend.
You can also construct an ODEProblem and then linearize prob.f manually
So I can just use the f_x jacobian from that function and ignore the g_x jacobian? Will the sensitivity from winch inputs to kite movement still be correct, even though there are quasi-static tether point masses in between the winch and the kite?
If you call linearize, you should get the statespace model you’re after directly, with the algebraic equations simplified away, does this not work?
You can see the linear algebra used to perform this simplification here
I am getting an error:
julia> [s.full_sys.set_values...]
3-element Vector{Num}:
sys₊set_values[1]
sys₊set_values[2]
sys₊set_values[3]
julia> [s.full_sys.ω_b...]
3-element Vector{Num}:
(sys₊ω_b(t))[1]
(sys₊ω_b(t))[2]
(sys₊ω_b(t))[3]
julia> lin_fun, simplified_sys = linearization_function(s.full_sys, [s.full_sys.set_values...], [s.full_sys.ω_b...])
┌ Warning: An empty operating point was passed to `linearization_function`. An operating point containing the variables that will be changed in `linearize` should be provided. Disable this warning by passing `warn_empty_op = false`.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/linearization.jl:51
ERROR: Some specified inputs were not found in system. The following variables were not found Any[sys₊set_values[3], sys₊set_values[2], sys₊set_values[1]]
Stacktrace:
[1] error
@ ./error.jl:44
[2] #markio!#824
@ ~/.julia/packages/ModelingToolkit/oBkpS/src/linearization.jl:577
[3] markio!
@ ~/.julia/packages/ModelingToolkit/oBkpS/src/linearization.jl:548 [inlined]
[4] _structural_simplify!(state::TearingState{…}, io::Tuple{…}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, dummy_derivative::Bool, kwargs::@Kwargs{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systemstructure.jl:735
[5] _structural_simplify!
@ ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systemstructure.jl:723 [inlined]
[6] structural_simplify!(state::TearingState{…}, io::Tuple{…}; simplify::Bool, check_consistency::Bool, fully_determined::Bool, warn_initialize_determined::Bool, kwargs::@Kwargs{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systemstructure.jl:681
[7] __structural_simplify(sys::ODESystem, io::Tuple{…}; simplify::Bool, sort_eqs::Bool, kwargs::@Kwargs{…})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systems.jl:93
[8] structural_simplify(sys::ODESystem, io::Tuple{…}; additional_passes::Vector{…}, simplify::Bool, split::Bool, allow_symbolic::Bool, allow_parameter::Bool, conservative::Bool, fully_determined::Bool, kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systems.jl:34
[9] structural_simplify
@ ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/systems.jl:29 [inlined]
[10] io_preprocessing(sys::ODESystem, inputs::Vector{Any}, outputs::Vector{Any}; simplify::Bool, kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/abstractsystem.jl:2565
[11] io_preprocessing
@ ~/.julia/packages/ModelingToolkit/oBkpS/src/systems/abstractsystem.jl:2563 [inlined]
[12] linearization_function(sys::ODESystem, inputs::Vector{…}, outputs::Vector{…}; simplify::Bool, initialize::Bool, initializealg::Nothing, initialization_abstol::Float64, initialization_reltol::Float64, op::Dict{…}, p::SciMLBase.NullParameters, zero_dummy_der::Bool, initialization_solver_alg::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{…}, autodiff::AutoForwardDiff{…}, eval_expression::Bool, eval_module::Module, warn_initialize_determined::Bool, guesses::Dict{…}, warn_empty_op::Bool, kwargs::@Kwargs{})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/linearization.jl:61
[13] linearization_function(sys::ODESystem, inputs::Vector{Num}, outputs::Vector{Num})
@ ModelingToolkit ~/.julia/packages/ModelingToolkit/oBkpS/src/linearization.jl:33
[14] top-level scope
@ REPL[32]:1
Some type information was truncated. Use `show(err)` to see complete types.
I reported a bug for this here.
Edit: this was just a mistake on my end.
When running linearize on my ODESystem the algebraic variables are not simplified away.
The simplified system has 121 unknowns, of which around 20 are differential (the rest is algebraic), but I am getting a matrix A of shape 121x121.
Some of the equations:
julia> equations(simplified_sys)
121-element Vector{Equation}:
Differential(t)((Q_b_w(t))[1]) ~ (Q_vel(t))[1]
Differential(t)((Q_b_w(t))[2]) ~ (Q_vel(t))[2]
Differential(t)((Q_b_w(t))[3]) ~ (Q_vel(t))[3]
Differential(t)((Q_b_w(t))[4]) ~ (Q_vel(t))[4]
Differential(t)((free_twist_angle(t))[1]) ~ (twist_ω(t))[1]
⋮
0 ~ 0.5(drag_force(t))[3, 2] + 0.5(drag_force(t))[3, 5] + 0.5(drag_force(t))[3, 9] - (point_force(t))[3, 9] - (spring_force_vec(t))[3, 2] - (spring_force_vec(t))[3, 5] + (spring_force_vec(t))[3, 9]
0 ~ 0.5(drag_force(t))[3, 31] + 0.5(drag_force(t))[3, 32] - (point_force(t))[3, 31] - (spring_force_vec(t))[3, 31] + (spring_force_vec(t))[3, 32]
0 ~ 0.5(drag_force(t))[3, 33] + 0.5(drag_force(t))[3, 34] - (point_force(t))[3, 33] - (spring_force_vec(t))[3, 33] + (spring_force_vec(t))[3, 34]
0 ~ 0.5(drag_force(t))[3, 35] + 0.5(drag_force(t))[3, 36] - (point_force(t))[3, 35] - (spring_force_vec(t))[3, 35] + (spring_force_vec(t))[3, 36]
0 ~ 0.5(drag_force(t))[3, 37] + 0.5(drag_force(t))[3, 38] - (point_force(t))[3, 37] - (spring_force_vec(t))[3, 37] + (spring_force_vec(t))[3, 38]
But you should get 121 linear differential equations, no algebraic variables left. You can then perform Model-order reduction on that model to reduce it to your desired size