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
A clarification, the simplified_sys object you printed hasn’t yet been linearized, but the matrices A,B,C,D you obtain when that system is linearized represent a linearized system with only differential equations.
I understand that the matrices represent the state space system, but these A, B, C, D matrices clearly contain non-differential variables as well. Maybe I understand “differential equations” wrong, but I would not see this equation for example as a differential equation, but as an algebraic equation:
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]
So I would not expect the variable related to this equation to show up in the A, B, C, D matrices. But they clearly do show up, as the A matrix has size 121x121, so the linear model contains all variables, also the algebraic ones.
What are algebraic variables in the nonlinear model may not be algebraic in the linearized model, you can differentiate the index-1 algebraic equations until you arrive at differential equations (index 0), which is what the linearized model represents.
MTK performs index reduction down to index 1 for algebraic equations hat are originally of index greater than 1, but you can do this down to index 0 as well.
Thanks for the explanation. I will try to use ModelOrderReduction.jl to reduce the order of the model.
Since you have an input-output system, I’m not sure the package you link is the most suitable. See Home · RobustAndOptimalControl Documentation for model-order reduction of linear statespace systems with inputs, in particular, try
using ControlSystemsBase
sys = ss(matrices...)
baltrunc(sys, residual=true)
baltrunc(sys, n=desired_state_dimension, residual=true)
or
using RobustAndOptimalControl
using DSP
w1 = 1e-4 # Lower frequency limit
w2 = 1e1 # upper
fc = DSP.analogfilter(DSP.Bandpass(w1, w2), DSP.Butterworth(2))
tfc = DSP.PolynomialRatio(fc)
W = tf(DSP.coefb(tfc), DSP.coefa(tfc))
rsys, hs = frequency_weighted_reduction(sys, W, 1)
for a reduction that focuses on a particular frequency range
I am not getting good results with using the continuous time linearized system, I think because the tether model is very stiff. Even without reducing the model order, the linear model is a lot more unstable than the nonlinear model. In the image below I re-linearize every 0.5 seconds. Reducing the model order gives similar results.
I want to remove the tether states directly by running the algebraic part of the solver inside the linearization process in the hope that the system will get less nonlinear and stiff.
(; A, B, C, D) = KiteModels.linearize(s)
csys = ss(A, B, C, D)
dsys = c2d(csys, dt)
I will try linearizing the discrete nonlinear control function instead of linearizing the continuous control function. This will probably give better results, but is also way slower, so I might need to implement gain scheduling for it to be real-time capable.
It’s unclear to me what the plot is showing, how is the linear simulation performed? Is it open or closed loop? Why does the state restart at zero every 0.5s?
If you linearize around the trajectory of the nonlinear simulation, do the linearizations look very different?
The simulation is open loop. The linear model state is reset to the nonlinear state at the same time the model is linearized every 0.5s, because the error will get very large otherwise.
Thanks for the link! I will try the ControlSystemsMTK.jl functions, also to make sure I didn’t do anything wrong with the linearization.
If the model is unstable, you will always get divergence of the linear system in open-loop simulation, always. See, e.g., Identification of unstable systems for an example and some motivation. A more reasonable way to judge the linearization in such a case may be to use closed-loop simulation, is this an option?
Thanks, I did not know that… I can try to use the linear system in MPC, and hope for the best. To close the loop I could try to add a simple PID controller to the system?
