You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Do you know a good references/textbook about nonlinear kalman filter and DAE systems @baggepinnen ?
The funny part is the only straightforward one to implement will be the MHE, since it's only implementing a transcription that support DAE. Come to think of it, I may start to work on #293, since it will pave the way to OrthogonalCollocation, thus at least one StateEstimator that supports DAE.
I do not, I haven't found any good references when looking into this myself. I generally approach this by explicitly modeling the disturbance inputs rather than a blanket $w$ like $$dx = Ax + Bu + w$$
since such a $w$ may bring the state out of the manifold described by the algebraic equations. An UKF with noise inputs works reasonably well on the formulation $$\dot x = f(x,u,p,t,w)$$
as long as $w$ is included in the integration. It sometimes requires some creativity to handle semidefinite covariance matrices though
I almost finished the support of OrthogonalCollocation for the MovingHorizonEstimator. The next step will be to support DAE systems with this specific estimator + a NonLinMPC with the same transcription method.
I was looking at the documentation of SimpleColloc. Is there a specific reason why you "merged" the dynamics equation $f$ and the algebraic equation $g$ into one single dyn function (apart being consistent with ModelingToolkit convention)?
I'm asking the question since I would implement this API ($\mathbf{x}$ is the state variable vector and $\mathbf{z}$ is the algebraic variable vector):
since it seems more intuitive and natural to me. I will not name the algebraic equation $\mathbf{g}$ to avoid any confusions with the nonlinear inequality constraint equation. About this, the algebraic vector $\mathbf{z}$ could be potentially confused with the decision vector $\mathbf{Z}$, but at least the first one is always lower case, and the second one always upper case in the documentation. Does it seems distinct enough to you? I could name the algebraic vector $\mathbf{a}$ but its not a common notation at all.
The user would provide 3 functions:
f (or in-place version f!)
q (or in-place version q!)
h or (or in-place version h!)
Are you aware of any limitations with this standard DAE form?
Having it in a single function is very advantageous since there may be time consuming computations that are shared between $f$ and $g$ in your formulation, and you would have to perform those twice if you split into two functions.
About the API for DAES, I'm thinking right now of having a single function, but still keeping the state and algebraic variable cleanly separated, since it is clearer for the state estimation problem (the MHE). Something like f_q!(ẋ, RHS, x, z, u, d, p) -> nothing or f_q(x, z, u, d, p) -> ẋ, RHS.
About initialization, similarly to what MTK.jl does, the algebraic variables will still need to be consistently initialized before starting the control loop (for a given set of differential state). I will likely implement this in the initstate! function, maybe by calling Ipopt as a root solver, to avoid new dependencies? What do you do in Dyad Control for intializing DAEs before starting the control loop?
Ok I think initialization of the algebraic variable through initstate! is not strictly necessary. The initial $\mathbf{z}(t=0)$ value can be included in the decision variables at $t=0$. This is another reason to keep the state cleanly separated from the algebraic variables: the $\mathbf{z}(t=0)$ will be included in the decision variable but not $\mathbf{\hat{x}}(t=0)$.
Ok, before anything else, I'm implementing a simple DAE solver based on the two collocation methods currently available, for simulating DAE model in open-loop in the package. 1. It will improve my understanding of DAEs. 2. This is useful to the user for sanity-checking the plant model outside of a MPC/MHE framework, and also to do simple closed-loop simulations.
I don't want to add an additional dependency just for this simple functionality (a root solver package). So I will use JuMP and Ipopt as a root solver, since everything is already there. It also provides an identical "solving environment" for open-loop simulation and for MPC/MHE solving. This is what gekko does BTW, not computationally optimal, but there are other advantages. So the following strictly refers to open-loop simulations of DAE systems.
As I said above, I just understood that I don't even need an initial value for the algebraic variable $\mathbf{z}$ to start the simulations. Let's take the simplest case: a TrapezoidalCollocation. If I know the current differential state $\mathbf{x}(k=0)$, we solve an optimization problem (root solving) at $k=0$ with $\mathbf{z}(k=0)$, $\mathbf{z}(k=1)$ and $\mathbf{x}(k=1)$ as the decision variables. We solve and extract the optimal $\mathbf{x}(k+1)$ to update the state. At time $k=1$ we warm-start Ipopt with $\mathbf{z}(k=1)$ value from the last time step.
Do you agree with my understanding of the problem @baggepinnen ? Thanks!
Do you agree with my understanding of the problem @baggepinnen ?
Yes. Having a start value or initial guess for the algebraic variables may still matter in some cases since the algebraic equations may have multiple solutions, but only one of them may be desirable.
Ok no, I made a mistake. We need the initial $\mathbf{z}$ value to compute the output $\mathbf{y}$. I want to support the current API with the new NonLinModelDAE type, that is:
model =NonLinModelDAE(...)
u, d =# starting valuesinitstate!(model, u, d)
for k=0:10
d =# some simulated values
y =evaloutput(model, d) # we need the current algebraic value `z` to compute the output `y`, so it needs to be initialized in `initstate!`updatestate!(model, u, d) # compute next state `xnext` and algebraic var. `znext` and store it internally in `model`end
I cannot compute $\mathbf{z}(k+1)$ in initstate! or updatestate! and use this precomputed value afterward in evalouput. The measured disturbance d may change between this computation and the evaloutput call, and this change may violate the algebraic equations $\mathbf{q}$. The only valid way is to solve another root solving problem directly inside evaloutput method, that is, solve the algebraic equation to extract the current algebraic variable $\mathbf{z}(k)$, and then extract the current model output $\mathbf{y}(k)$.
In terms of API, it will look like this:
model =NonLinModelDAE(...)
u, d =# starting valuesinitstate!(model, u, d) # this will do nothing at all, similarly to `NonLinModel` right now (for ODEs)for k=0:10
d =# some simulated values
y =evaloutput(model, d) # will solve 0 = q(x, z, u, d, p) to compute the current `z`, and then it will call h(x, z, d, p)updatestate!(model, u, d) # solve the collocation problem to compute next state `xnext`end
An important precision: I will need to use the previous u value to solve the root problem in evaloutput. But it is not an issue here since the strictly proper systems w.r.t. u restriction will still apply for NonLinModelDAE. So, by definition, if u impacts some specific algebraic variables, they are "orthogonal" to the output equation and will not corrupt y value. I will likely validate that the NonLinModelDAE is strictly proper by using SparseConnectivityTracer.jl at construction and output a clear error message otherwise (it's impossible to impose strictly proper DAEs just by forcing specific function signatures without being too restrictive).
Activity
franckgaga commentedon Mar 22, 2026
Do you know a good references/textbook about nonlinear kalman filter and DAE systems @baggepinnen ?
The funny part is the only straightforward one to implement will be the MHE, since it's only implementing a transcription that support DAE. Come to think of it, I may start to work on #293, since it will pave the way to
OrthogonalCollocation, thus at least oneStateEstimatorthat supports DAE.baggepinnen commentedon Mar 23, 2026
I do not, I haven't found any good references when looking into this myself. I generally approach this by explicitly modeling the disturbance inputs rather than a blanket$w$ like
$$dx = Ax + Bu + w$$ $w$ may bring the state out of the manifold described by the algebraic equations. An UKF with noise inputs works reasonably well on the formulation
$$\dot x = f(x,u,p,t,w)$$ $w$ is included in the integration. It sometimes requires some creativity to handle semidefinite covariance matrices though
since such a
as long as
franckgaga commentedon Apr 11, 2026
I just found this for EKF and UKF:
https://www.sciencedirect.com/science/article/pii/S1474667015303724?via%3Dihub
franckgaga commentedon Aug 12, 2026
Hello @baggepinnen!
I almost finished the support of
OrthogonalCollocationfor theMovingHorizonEstimator. The next step will be to support DAE systems with this specific estimator + aNonLinMPCwith the same transcription method.I was looking at the documentation of$f$ and the algebraic equation $g$ into one single
SimpleColloc. Is there a specific reason why you "merged" the dynamics equationdynfunction (apart being consistent withModelingToolkitconvention)?I'm asking the question since I would implement this API ($\mathbf{x}$ is the state variable vector and $\mathbf{z}$ is the algebraic variable vector):
since it seems more intuitive and natural to me. I will not name the algebraic equation$\mathbf{g}$ to avoid any confusions with the nonlinear inequality constraint equation. About this, the algebraic vector $\mathbf{z}$ could be potentially confused with the decision vector $\mathbf{Z}$ , but at least the first one is always lower case, and the second one always upper case in the documentation. Does it seems distinct enough to you? I could name the algebraic vector $\mathbf{a}$ but its not a common notation at all.
The user would provide 3 functions:
f(or in-place versionf!)q(or in-place versionq!)hor (or in-place versionh!)Are you aware of any limitations with this standard DAE form?
baggepinnen commentedon Aug 24, 2026
Having it in a single function is very advantageous since there may be time consuming computations that are shared between$f$ and $g$ in your formulation, and you would have to perform those twice if you split into two functions.
franckgaga commentedon Aug 24, 2026
Hum very good point. I will change the API on my WIP branch.
franckgaga commentedon Aug 28, 2026
Ok DAEs are more complex than I first though...
About the API for DAES, I'm thinking right now of having a single function, but still keeping the state and algebraic variable cleanly separated, since it is clearer for the state estimation problem (the MHE). Something like
f_q!(ẋ, RHS, x, z, u, d, p) -> nothingorf_q(x, z, u, d, p) -> ẋ, RHS.About initialization, similarly to what MTK.jl does, the algebraic variables will still need to be consistently initialized before starting the control loop (for a given set of differential state). I will likely implement this in the
initstate!function, maybe by callingIpoptas a root solver, to avoid new dependencies? What do you do in Dyad Control for intializing DAEs before starting the control loop?franckgaga commentedon Aug 31, 2026
Ok I think initialization of the algebraic variable through$\mathbf{z}(t=0)$ value can be included in the decision variables at $t=0$ . This is another reason to keep the state cleanly separated from the algebraic variables: the $\mathbf{z}(t=0)$ will be included in the decision variable but not $\mathbf{\hat{x}}(t=0)$ .
initstate!is not strictly necessary. The initialfranckgaga commentedon Sep 2, 2026
Ok, before anything else, I'm implementing a simple DAE solver based on the two collocation methods currently available, for simulating DAE model in open-loop in the package. 1. It will improve my understanding of DAEs. 2. This is useful to the user for sanity-checking the plant model outside of a MPC/MHE framework, and also to do simple closed-loop simulations.
I don't want to add an additional dependency just for this simple functionality (a root solver package). So I will use JuMP and Ipopt as a root solver, since everything is already there. It also provides an identical "solving environment" for open-loop simulation and for MPC/MHE solving. This is what gekko does BTW, not computationally optimal, but there are other advantages. So the following strictly refers to open-loop simulations of DAE systems.
As I said above, I just understood that I don't even need an initial value for the algebraic variable$\mathbf{z}$ to start the simulations. Let's take the simplest case: a $\mathbf{x}(k=0)$ , we solve an optimization problem (root solving) at $k=0$ with $\mathbf{z}(k=0)$ , $\mathbf{z}(k=1)$ and $\mathbf{x}(k=1)$ as the decision variables. We solve and extract the optimal $\mathbf{x}(k+1)$ to update the state. At time $k=1$ we warm-start Ipopt with $\mathbf{z}(k=1)$ value from the last time step.
TrapezoidalCollocation. If I know the current differential stateDo you agree with my understanding of the problem @baggepinnen ? Thanks!
baggepinnen commentedon Sep 3, 2026
Yes. Having a start value or initial guess for the algebraic variables may still matter in some cases since the algebraic equations may have multiple solutions, but only one of them may be desirable.
franckgaga commentedon Sep 4, 2026
Ok no, I made a mistake. We need the initial$\mathbf{z}$ value to compute the output $\mathbf{y}$ . I want to support the current API with the new
NonLinModelDAEtype, that is:franckgaga commentedon Sep 5, 2026
Still wrong! Wow ODEs are so much simpler! 🤯
I cannot compute$\mathbf{z}(k+1)$ in $\mathbf{q}$ . The only valid way is to solve another root solving problem directly inside $\mathbf{z}(k)$ , and then extract the current model output $\mathbf{y}(k)$ .
initstate!orupdatestate!and use this precomputed value afterward inevalouput. The measured disturbancedmay change between this computation and theevaloutputcall, and this change may violate the algebraic equationsevaloutputmethod, that is, solve the algebraic equation to extract the current algebraic variableIn terms of API, it will look like this:
doc: correct `Z` for DAEs
franckgaga commentedon Sep 7, 2026
An important precision: I will need to use the previous
uvalue to solve the root problem inevaloutput. But it is not an issue here since the strictly proper systems w.r.t.urestriction will still apply forNonLinModelDAE. So, by definition, ifuimpacts some specific algebraic variables, they are "orthogonal" to the output equation and will not corruptyvalue. I will likely validate that theNonLinModelDAEis strictly proper by usingSparseConnectivityTracer.jlat construction and output a clear error message otherwise (it's impossible to impose strictly proper DAEs just by forcing specific function signatures without being too restrictive).