SOSTOOLS for differential-algebraic systems?

I am trying to find a Lyapunov function for a differential-algebraic system (DAE) of the form:

,

where () represents the dynamic states and () represents the algebraic variables.

Due to certain constraints, I cannot eliminate the algebraic variables, i.e., I cannot explicitly write ().

Has anyone approached this type of problem using SOSTOOLS while keeping the algebraic variables in the formulation? Any suggestions on how to proceed, or any mini examples and references on SOS-based Lyapunov analysis for DAE systems, would be highly appreciated.

[Torsten]

The MATLAB tools to solve such systems are ode15s or ode15i (without Lyapunov function derivations).

[Sam Chak]

Hi @Husni Rois Ali

How do the algebraic variables () influence the search for a Lyapunov function, , which is primarily a function of the state variables ()?

By the way, can the in the expression be treated as coefficients so that you can continue with the regular procedure of searching an abstract, non-physical SOS Lyapunov function that yields for ?

[Husni Rois Ali]

Hi @Sam Chak

Thank you for your answer.

Regarding how (y) influences (x), this is actually the part that is still unclear to me. I am not sure whether the Lyapunov candidate should be formulated as (V(x)) only, or whether it should explicitly include the algebraic variables, i.e., (V(x,y)).

If you have any experience dealing with this type of DAE problem or SOS-based Lyapunov construction, please let me know :-)

[Sam Chak]

Hi @Husni Rois Ali

Your system is defined as dx/dt = f(x, y) subject to the constraint g(x, y) = 0.

If your Lyapunov function is a function of x only, then

dV(x)/dt = dV/dx*dx/dt.

But if your Lyapunov function is a function of x and y, then

dV(x,y)/dt = ∂V/∂x*dx/dt + ∂V/∂y*dy/dt.

However, the derivative dy/dt is unavailable. Moreover, are you trying to prove that the non-system states as via the Lyapunov function?

You can find some examples in Getting started with Sum of Squares and the SOSTOOLS user's guide.

[Torsten]

However, the derivative dy/dt is unavailable.

dg/dx*dx/dt + dg/dy*dy/dt = dg/dx*f(x,y) + dg/dy*dy/dt = 0

thus

dy/dt = -(dg/dy)^(-1) * dg/dx*f(x,y)

[Sam Chak]

Hi @Husni Rois Ali

Please follow @Torsten's insightful approach to derive dy/dt from the equality constraint, . That way, you can formulate the SOS Lyapunov function as , if you wish to also prove that as .

Just to be sure, your purpose in finding the Lyapunov function is to prove that the system eventually comes to rest at the equilibrium, , correct? Suppose you manage to find an SOS Lyapunov function and show that , which implies that and as . When this occurs, does the constraint still hold?

[Husni Rois Ali]

Hi @Sam Chak and @Torsten

Thank you for the insightful responses and helpful suggestions.

Let me clarify the problem in a bit more detail.

The system I am considering is a differential–algebraic system (DAE) of the form

where

are the differential states and

E

is the algebraic variable.

More specifically, the model is

with power equations[

where

The algebraic equation is

When I mentioned that the algebraic variables cannot be explicitly solved, that statement was not entirely accurate. In fact, the algebraic equation can be solved locally for (E). Since it is quadratic in (E), we can write

and hence

So explicit elimination is possible, at least locally on the physical branch. However, to the best of my knowledge, performing this elimination would introduce square-root terms and destroy the polynomial structure of the system, which I am trying to preserve for SOS/SOSTOOLS analysis.

Regarding the purpose, my goal is essentially to find a Lyapunov function, either

or

that can be used to show that the shifted system returns to the equilibrium (origin).

Thank you again for the helpful explanations. I will follow the approaches suggested above, particularly the manifold/equality-constraint treatment and the derivation of (dE/dt) from the algebraic constraint. I appreciate the guidance.

Answer 1

[Sam Chak]

Hi @Husni Rois Ali,

This is not a complete answer in the MATLAB or SOSTOOLS sense, but it should provide insight into approaching the sum-of-squares and stability problem.

Since the first two differential states do not explicitly depend on the third differential state, you can evaluate the stability of the first two states independently. In the following example, the system is coordinate-transformed so that the equilibrium is at the origin.

Given that

where . The system can be trigonometrically simplified to

.

Solving

to obtain the equilibrium at , where .

A Lyapunov function candidate:

Taking the derivative of V along the trajectory of the system:

You see that , V is positive definite, and is negative definite, for and over . Thus, by the Lyapunov stability theorem, you can conclude that the origin ( and ) is asymptotically stable.

The third differential equation does not contain the third differential state , but you can use the equality constraint

to infer that

.

Substituting this result back to

where

then you should get something like

.

If and is constant, then the third differential state should be asymptotically stable as well.

[Husni Rois Ali]

Thank you @Sam Chak for the illustrative example. I will go through it carefully.

Regarding the approach suggested by @Torsten, my understanding is that it leads to a division by a polynomial term. I am aware of the denominator-clearing trick, where multiplying by the denominator preserves the Lyapunov certificate for (V(x,y)) as long as the denominator remains positive. However, this positivity condition may not always hold in my case.

Do you happen to have any alternative suggestions that would avoid relying on the sign of the denominator?

[Sam Chak]

Hi @Husni Rois Ali

You will need to verify if the mathematical derivation according to @Torsten's approach is correct.

While you can eliminate the denominator of a term by multiplying both sides of an equation (dV(x,y)/dt = ∂V/∂x*dx/dt + ∂V/∂y*dy/dt) by that denominator, I don't remember there is rule stating that the trick is only valid when the denominator is positive. However, the multiplying trick is valid, provided that the denominator does not equal zero. Are you implying that the denominator dg/dE can approach zero due to electrical power oscillations? If dg/dE can approach zero, then dE/dt can approach singularity too.

What do you benefit from eliminating the denominator? Won't this making the equation monstrously long?

syms f [3 1] matrix
syms x [3 1] matrix
syms E eta X V0 R Q E0 w P0 H k Q0 D0
% w = ωb; V0 = Vg; P0 = Pref; Q0 = Qref; D0 = Dp

% vector definition
delta = x(1); 
omega = x(2); 
Phi   = x(3);

% expressions
D   = R^2 + X^2

D = 

P   = ((E^2)*R*E*V0*R*cos(delta) + E*V0*X*sin(delta))/D

P = 

Q   = ((E^2)*X*E*V0*X*cos(delta) + E*V0*R*sin(delta))/D

Q = 

% ODEs
f(1)  = w*omega;
f(2)  = (P0 - P - D0*omega)/(2*H);
f(3)  = k*(Q0 - Q)

f = 

% equality constraint
g   = eta*X*E^2 + (D - eta*V0*X*cos(delta) - eta*V0*R*sin(delta))*E - D*(eta*Q + Phi + E0)

g = 

% gradient with respect to the vector x
dgx = diff(g, x)

dgx = 

% gradient with respect to E
dgE = diff(g, E)

dgE = 

% deriving dE/dt using Torsten's approach
dE  = - dgx*f/dgE

dE =