Support DAE

[franckgaga]

The hardest part will be how I'll be able to support all the state estimators. Maybe a subset of it, as a start.

[franckgaga]

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.

[baggepinnen]

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

[franckgaga]

I just found this for EKF and UKF:

https://www.sciencedirect.com/science/article/pii/S1474667015303724?via%3Dihub

[franckgaga]

Hello @baggepinnen!

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):

$$\begin{aligned} \mathbf{\dot{x}} &= \mathbf{f}(\mathbf{x, z, u, d, p}) \\\ \mathbf{0} &= \mathbf{q}(\mathbf{x, z, u, d, p}) \\\ \mathbf{y} &= \mathbf{h}(\mathbf{x, z, d, p}) \\\ \end{aligned}$$

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:

1. f (or in-place version f!)
2. q (or in-place version q!)
3. h or (or in-place version h!)

Are you aware of any limitations with this standard DAE form?

[baggepinnen]

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]

Hum very good point. I will change the API on my WIP branch.

[franckgaga]

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) -> 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?

[franckgaga]

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)$.

[franckgaga]

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!

[baggepinnen]

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.

[franckgaga]

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 values
initstate!(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

[franckgaga]

Still wrong! Wow ODEs are so much simpler! 🤯

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 values
initstate!(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