Using Interval Analysis to Compute the Invariant Set of a Nonlinear Closed-Loop Control System

Abstract

In recent years, many applications, as well as theoretical properties of interval analysis have been investigated. Without any claim for completeness, such applications and methodologies range from enclosing the effect of round-off errors in highly accurate numerical computations over simulating guaranteed enclosures of all reachable states of a dynamic system model with bounded uncertainty in parameters and initial conditions, to the solution of global optimization tasks. By exploiting the fundamental enclosure properties of interval analysis, this paper aims at computing invariant sets of nonlinear closed-loop control systems. For that purpose, Lyapunov-like functions and interval analysis are combined in a novel manner. To demonstrate the proposed techniques for enclosing invariant sets, the systems examined in this paper are controlled via sliding mode techniques with subsequently enclosing the invariant sets by an interval based set inversion technique. The applied methods for the control synthesis make use of a suitably chosen Gröbner basis, which is employed to solve Bézout’s identity. Illustrating simulation results conclude this paper to visualize the novel combination of sliding mode control with an interval based computation of invariant sets.

1. Introduction

During the last few decades, interval analysis has been shown to provide efficient numerical approaches for solving various tasks in scientific computing, as well as in application-oriented disciplines in which numerical techniques with a result verification are helpful or required. In such a way, they help to solve problems in control engineering, robotics, including rigorous methods for localization and/or path planning, as well as guaranteed obstacle avoidance, computational (bio-)mechanics, civil engineering, or bio-engineering, just to mention a few of the large number of different areas [1,2,3,4,5,6,7].

In these applications, not only the possibility of interval analysis to treat the effects of round-off and truncation errors of numerical solution techniques (such as truncated temporal Taylor series expansions for the solution of initial value problems to ordinary differential equations) is exploited [8,9,10]. In addition, these applications also benefit from a further fundamental property of interval analysis, namely that it provides guaranteed outer enclosures of all possible solutions to a mathematical problem at hand in the presence of ambiguities in the solutions or in the presence of bounded uncertainty in selected parameters.

Here, a selection of useful techniques ranges from the extension of fundamental arithmetic operations to interval valued expressions [1,2,3], which were recently standardized in the IEEE standard 1788 [11], over the development of set valued counterparts for zero finding techniques for sets of algebraic equations such as the Krawczyk operator [12,13,14], interval based techniques for reachability analysis [15], the verified global optimization [16,17,18], to the solution of identification tasks by means of set inversion techniques via interval analysis (e.g., the algorithms set inversion via interval analysis (SIVIA) and problem specific generalizations) [2,19]. If control applications are concerned, interval analysis can be applied both to an offline design and verification of control procedures under consideration of feasibility and safety requirements such as input and state constraints and to an online interval evaluation in terms of real-time capable robust control strategies generalizing the ideas of variable-structure control techniques and backstepping [20,21,22,23,24,25]. Further, approaches for the control of systems with input saturations that do not rely on interval analysis were handled, for example, in [26] and the references therein. In addition to the implementation of robust controllers, also the dual task of interval based state estimation and observer design can be solved; see, for example, [27,28,29,30,31,32].

This paper addresses the following two issues: First, a method for finding a nonlinear state-feedback controller is derived for the exact tracking of predefined paths in the phase-space. These paths are assumed to be described in terms of a sliding surface for a nonlinear dynamic system. Second, this control design procedure is interfaced with techniques from the field of interval analysis that are employed to determine enclosures of the invariant sets of the resulting controlled system in a reliable and computationally efficient way. For a continuous-time system model with the state vector $\mathit{x}\in {\mathbb{R}}^n$, a subset $S\subset {\mathbb{R}}^n$ is called positively invariant if the relation $\mathit{x}\left(t_1\right)\in S$ is satisfied for all possible $\mathit{x}\left(t_0\right)\in S$ if $t_1\ge t_0$ holds [33]. Due to the fact that we are only interested in a forward-in-time propagation of state trajectories, these sets are briefly called invariant throughout this paper. For the proposed control design, an appropriate Gröbner basis is used to check if the occurring Bézout identity has a solution, such that a controller can be designed that does not contain any structural singularities.

To find the invariant set of the respective closed-loop control system, Lyapunov-like functions are employed, which transform the stability properties of dynamical systems described by a set of ordinary differential equations into inequalities. Those inequalities can be evaluated using interval computation and finally lead to the desired enclosures of the invariant sets to be determined. Note that the proposed technique does not focus on determining invariant sets by using the La Salle invariance theorem. Instead, only temporal derivatives of the aforementioned Lyapunov-like functions are required during the solution procedure.

To depict the examined problem and to distinguish it from a state-of-the-art sliding mode control design [34,35,36,37,38,39], consider a single-input system:

$$\dot{\mathit{x}}\left(t\right)=\mathit{f}(\mathit{x}\left(t\right),u\left(\mathit{x}\left(t\right)\right)),$$

(1)

where $\mathit{x}\in {\mathbb{R}}^n$ is the state vector and $\mathit{f}:{\mathbb{R}}^n\mapsto {\mathbb{R}}^n$ is the evolution function. For simplicity, assume that the input in (1) is given by the state dependent expression $u\left(\mathit{x}\right(t\left)\right)\in \mathbb{R}$. The time dependencies of the state vector $\mathit{x}:=\mathit{x}\left(t\right)\in {\mathbb{R}}^n$ and of the control vector $u:=u\left(t\right)$ are furthermore not explicitly denoted. Additionally, consider a sliding surface:

$$s\left(\mathit{x}\right)=0,\mathrm{where}s:{\mathbb{R}}^n\mapsto \mathbb{R}.$$

(2)

For the system model defined in Equation (1), the control input u should be found in such a way that the function $\frac{1}{2}{s}^2\left(\mathit{x}\right)$, which can be interpreted as a Lyapunov function candidate, decreases monotonically up to the point where perfect tracking of the predefined path (respectively, trajectory) occurs due to $s\left(\mathit{x}\right)=0$. As a consequence, the state trajectories of the controlled system will converge to an invariant set in which the sliding surface should be contained.

The paper is structured in the following way. Firstly, the general concept of sliding mode approaches for the control design of nonlinear systems is described in Section 2 with the aim of distinguishing the investigated problem from classical techniques in this field. Then, the proposed novel method for a sliding mode-type controller design is explained in Section 3 in combination with the use of a Gröbner basis to solve Bézout’s identity. Section 4 deals with the conditions to be satisfied for an invariant set, which are evaluated numerically by using interval analysis as described in Section 5. In Section 6, several applications of the developed method are presented. These applications include cases with sliding surfaces represented by closed curves in the phase space. Moreover, straight line segments are considered as well. In addition, a concluding example is investigated where Bézout’s identity is not solvable. However, the interval based method for enclosing the resulting invariant set of the controlled system can still be employed successfully if a bounded control input for the considered nonlinear system of dimension three is analyzed. Finally, conclusions and an outlook toward future work together with an explanation of how the presented interval technique can be applied to other nonlinear control techniques are given in Section 7.

2. Classical Sliding Mode Control

The main concept in sliding mode control is to replace an nth order tracking problem by a stabilization problem of the order $(n-1)$ [40]. Looking at the dynamic system:

$$\dot{\mathit{x}}=\mathit{f}(\mathit{x},u),$$

(3)

the guaranteed stabilizing, variable structure input signal u is chosen according to:

$$u=\left\{\begin{array}{cc}{u}^{+},\hfill & \mathrm{if}s(\mathit{x},t)>0\hfill \\ u_{\mathrm{eq}},\hfill & \mathrm{if}s(\mathit{x},t)=0\hfill \\ u^{-},\hfill & \mathrm{if}s(\mathit{x},t)<0.\hfill \end{array}\right.$$

(4)

For that purpose, a switching function $s(\mathit{x},t)$ needs to be defined. In classical sliding mode procedures, which are applied to solve trajectory tracking tasks, this switching function is usually composed of the output tracking error and a weighted linear combination of the temporal derivatives of this error signal up to the order $\delta -1$, where $\delta$ denotes the system’s relative degree [24]. In the remainder of this section, it is assumed that the relative degree $\delta$ equals the system order according to $n=\delta$. Then, a stability proof of the closed-loop dynamics can be performed in terms of a square Lyapunov function candidate as was mentioned in the Introduction of this paper.

In general, the switching function $s(\mathit{x},t)$ determines under which conditions the input signal of the system changes structurally, so that u becomes discontinuous across $s(\mathit{x},t)=0$ [40,41]. In such a way, the sliding surface, which is defined by $s(\mathit{x},t)=0$, divides the state space into two mutually disjoint domains where either $s(\mathit{x},t)>0$ or $s(\mathit{x},t)<0$ holds. In those domains, different dynamics of the system occur due to different definitions of the system’s input u. The sliding surface is included in neither of these subdomains as a trajectory, whereas each point of the sliding surface $s(\mathit{x},t)=0$ has to be accessible by a trajectory from either side [41]. To make the system converge to the sliding surface, the input has to be chosen in such a way that it always makes the state trajectories approach the sliding surface, i.e., they have to be “pushed” towards $s(\mathit{x},t)=0$. This phase is called the transient or reaching phase. The possibilities for the design of controllers for this reaching phase can be found, for example, in [42,43] and the references therein. Note that the design of the corresponding control for this phase needs to account for input constraints in terms of the limited range of the admissible system inputs u. If scenarios occur, in which the system input reaches its lower or upper saturation value in such a phase, a careful stability proof is inevitable, cf. [26].

The next phase, following the completion of the reaching phase, is called sliding mode. The system’s dynamics are then represented by the switching function itself. As mentioned before, this function typically has a lower order than the open-loop system [38,39].

Thus, in this state, the system’s dynamic behavior only depends on the parameters included in the switching function, which means that this form of control is not only robust against uncertainties in plant parameters, but also against disturbances satisfying suitable matching conditions (for example, disturbances entering the state Equations (3) via the same channels as the system inputs u). From an intuitive point of view, this results from the fact that even if the system leaves the sliding surface slightly, it will immediately return to it, because of the switching-type behavior of the input u and the proven stability (i.e., attraction property) of the sliding surface.

However, the switching between structurally different inputs as defined in (4) may lead to a phenomenon called chattering, which implies high control activity and thus leads to a rapid wear of the actuator. In addition, the associated high frequencies may also excite dynamics that were not considered in the model [40]. Strategies aiming at a reduction of the chattering phenomenon by a stability preserving interval based gain adaptation strategy were investigated in [20,21,22,23,24,25].

To ensure that the system will reach the sliding surface in finite time, the condition:

$$\frac{1}{2}\frac{\mathrm{d}{s}^2\left(\mathit{x}\right)}{\mathrm{d}t}=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)\le -\eta \left|s\left(\mathit{x}\right)\right|,\mathrm{where}\eta >0,$$

(5)

is typically introduced in a classical sliding mode control synthesis. The inequality (5) means that the squared distance to the sliding surface has to decrease along all possible trajectories with a constant velocity $\eta$. This condition makes the set $s=0$ an invariant set due to the invariance theorem of La Salle. If this condition is satisfied for the complete state space, the system will always converge to the sliding surface.

To reach the sliding mode, the state dependent controller has to be designed. There are various ways to do that, where we restrict ourselves to one of the available options. To deal with possible singularities in the control law and to avoid them during perfect trajectory tracking in the sliding mode, the input u is chosen according to the concept of equivalent control. To find the equivalent control $u_{\mathrm{eq}}$ that ensures perfect tracking, $\dot{s}=0$ should be valid on the sliding surface in order for s not to change and thus the system to remain in sliding mode. For a system like (1) with a (nonlinear) feedback controller $u\left(\mathit{x}\right)$, the time derivative of $s\left(\mathit{x}\right)$ is calculated by the following Lie derivative of $s\left(\mathit{x}\right)$ in the direction of the vector field $\mathit{f}(\mathit{x},u)$, assuming that $s\left(\mathit{x}\right)$ does not explicitly depend on the time variable t:

$$\dot{s}\left(\mathit{x}\right)={\left(\frac{\partial s\left(\mathit{x}\right)}{\partial \mathit{x}}\right)}^T·\mathit{f}(\mathit{x},u).$$

(6)

To find the equivalent input $u:=u_{\mathrm{eq}}$, this expression is solved with respect to u. However, if the system states are not on the sliding surface during the reaching phase, u should make it converge to $s\left(\mathit{x}\right)=0$. To ensure this behavior, a further term depending on the sign of s resulting from the right hand side of (5) is typically added to the equivalent control $u_{\mathrm{eq}}$ to obtain the overall control signal stated in Equation (4). This can be achieved by using a sign function (often regularized by means of smooth sigmoid-type expressions) depending on the current value of $s\left(\mathit{x}\right)$ [40,44].

3. Controller Design

To find the controller for the system model (1) with a time invariant sliding surface as considered in Equation (6), the following approach is used, which is similar to the one that can be found in [45] by setting:

$$\dot{s}\left(\mathit{x}\right)=-s\left(\mathit{x}\right)·p\left(\mathit{x}\right),$$

(7)

where $p\left(\mathit{x}\right)>0$. Note that this assumption is not only useful for the classical tracking control problem mentioned in the previous section, but it can also be employed to design control procedures allowing for following paths in the phase space that satisfy the design procedures listed in the following. From a dynamics point of view, the factor $p\left(\mathit{x}\right)$ in Equation (7) can be interpreted as a state dependent decay rate for the convergence of tracking errors towards the desired sliding mode.

With the choice given in Equation (7), $s\left(\mathit{x}\right)$ will be forced towards zero, because the resulting derivative $\dot{s}\left(\mathit{x}\right)$ always has the opposite sign of $s\left(\mathit{x}\right)$, leading to the strict inequality:

$$\frac{\mathrm{d}{s}^2\left(\mathit{x}\right)}{\mathrm{d}t}=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)=-s^2\left(\mathit{x}\right)·p\left(\mathit{x}\right)<0\mathrm{for}\mathrm{all}s\left(\mathit{x}\right)\ne 0.$$

(8)

With the help of the relation (7), the control signal can be found in symbolic form. As an example, consider a two-dimensional input affine system of the form:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2\hfill \\ \hfill {\dot{x}}_2& =c\left(\mathit{x}\right)+u,\hfill \end{array}$$

(9)

for which Equation (7) becomes equal to:

$$\dot{s}\left(\mathit{x}\right)=\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{\dot{x}}_1+\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}{\dot{x}}_2=-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)$$

(10)

as well as equal to:

$$\dot{s}\left(\mathit{x}\right)=\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2+\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}(c\left(\mathit{x}\right)+u)=-s\left(\mathit{x}\right)·p\left(\mathit{x}\right).$$

(11)

The latter expression can now be solved for the control input according to:

$$u=u\left(\mathit{x}\right)=\left(-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)-\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2\right)·\frac{1}{\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}}-c\left(\mathit{x}\right).$$

(12)

To avoid singularities in the controller and in this way in the closed-loop system dynamics, the following condition has to be satisfied:

$$-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)-\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2=k\left(\mathit{x}\right)\frac{\partial s\left(\mathit{x}\right)}{\partial x_2},$$

(13)

so that all fractions in the expression for the input u cancel out analytically. After simple rearrangement, this expression becomes:

$$\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}k\left(\mathit{x}\right)+s\left(\mathit{x}\right)·p\left(\mathit{x}\right)=-\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2,$$

(14)

which is Bézout’s identity, where $p\left(\mathit{x}\right)$ and $k\left(\mathit{x}\right)$ have to be found in terms of non-singular analytic (polynomial) expressions.

Remark 1.

The design procedure summarized above is more general than the classical sliding mode design, where the goal of tracking sufficiently smooth reference trajectories $y_{\mathrm{d}}:=y_{\mathrm{d}}\left(t\right)$ for the desired temporal evolution of the state variable $x_1$ would lead to the system input:

$$u=\alpha ·\left({\dot{y}}_{\mathrm{d}}-x_2\right)+{\ddot{y}}_{\mathrm{d}}+\eta ·\mathrm{sign}\left(\tilde{s}(\mathit{x},t)\right)-c\left(\mathit{x}\right).$$

(15)

Here, the parameter α is chosen as a positive constant (generally satisfying the Hurwitz criterion [20,22]) to ensure stability on the sliding surface and $\eta >0$ ensures a stable motion during the reaching phase towards the classical (time varying) sliding surface:

$$\tilde{s}(\mathit{x},t)=\alpha ·\left(y_{\mathrm{d}}-x_1\right)+\left({\dot{y}}_{\mathrm{d}}-x_2\right)=0.$$

(16)

To solve the Bézout identity (14) in the design procedure suggested in this paper, it has to be checked firstly whether a solution for this polynomial equation exists. This can be done using a reduced Gröbner basis. If a polynomial belongs to an ideal, it can be expressed as a linear combination of the polynomials that belong to the generating set of the ideal in a polynomial ring. This means if the term:

$$-\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2$$

(17)

included in the system input defined in Equation (12) is part of the ideal:

$$I_1=〈s\left(\mathit{x}\right),\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}〉,$$

(18)

Equation (14) has a solution.

The Gröbner basis can provide this generating set of polynomials, and the reduced Gröbner basis makes sure that none of the leading coefficients of the polynomials are divided by one another. All polynomials belonging to this ideal can be expressed by linear combinations of the reduced Gröbner basis elements.

Therefore, if the reduced Gröbner bases of $I_1$ and:

$$I_2=〈s\left(\mathit{x}\right),\frac{\partial s\left(\mathit{x}\right)}{\partial x_2},-\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}{x}_2〉$$

(19)

are the same, Equation (14) has a solution. As soon as a first solution for $p\left(\mathit{x}\right)$ and $k\left(\mathit{x}\right)$ has been found, multiple further solutions can be obtained by exploiting the fact that if:

$$p\left(\mathit{x}\right)·G_0+k\left(\mathit{x}\right)·G_1=B\left(\mathit{x}\right)$$

(20)

holds true, it directly implies that also the equality:

$$\begin{array}{cc}\hfill (p\left(\mathit{x}\right)+Y\left(\mathit{x}\right)·G_1)·G_0+(k\left(\mathit{x}\right)-Y\left(\mathit{x}\right)·G_0)·G_1& =\hfill \\ \hfill p\left(\mathit{x}\right)·G_0+Y\left(\mathit{x}\right)·G_0·G_1+k\left(\mathit{x}\right)·G_1-Y\left(\mathit{x}\right)·G_0·G_1& =B\left(\mathit{x}\right)\hfill \end{array}$$

(21)

is satisfied, where $G_0$ and $G_1$ are the elements of the Gröbner basis [46,47].

4. Computation of the Invariant Set of the Closed-Loop Control System

Once the controller $u\left(\mathit{x}\right)$ in Equation (12) is found, it should be shown that the system that converges to the sliding surface $s\left(\mathit{x}\right)=0$ that is characterized by (7) stays on the surface for all times. This is the definition of an invariant set due to La Salle. If a trajectory enters an invariant set, it will stay inside this set forever [33,40]. To find this set, the switching function $s\left(\mathit{x}\right)$ is used as a Lyapunov-like function $\tilde{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)$, which does not necessarily comply with all requirements of a Lyapunov function, as its derivative may also be zero for states $\mathit{x}$, where $V=\frac{1}{2}{s}^2\left(\mathit{x}\right)\ne 0$. Moreover, the function $\tilde{V}\left(\mathit{x}\right)$ may exhibit a change of sign.

Under this preliminary consideration of the characteristics of invariant sets, states are not part of such sets if they are related to an inflection point in the form of a saddle point of the function $V\left(\mathit{x}\right)$ (see Figure 1).

In fact, it suffices if any finite time derivative of $V\left(\mathit{x}\right)$ is not equal to zero for a point not belonging to the invariant set. This can also be confirmed when thinking of the Barbalat lemma, which states that, if a function has a finite limit for $t\to \infty$, its uniformly continuous derivative will tend to zero [40]. This can also be applied for the first derivative being the function and the second derivative being its first derivative, and so on:

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}V\left(\mathit{x}\right)=c\Rightarrow \underset{t\to \infty }{lim}\dot{V}\left(\mathit{x}\right)=0\Rightarrow \\ \hfill \underset{t\to \infty }{lim}\ddot{V}\left(\mathit{x}\right)=0\Rightarrow \underset{t\to \infty }{lim}\stackrel{⃛}{V}\left(\mathit{x}\right)=0\dots ,\\ \hfill \mathrm{where}c\mathrm{is}\mathrm{a}\mathrm{constant}.\end{array}$$

(22)

Using the fact that the desired sliding surface $s\left(\mathit{x}\right)$ itself and the first three derivatives of $s\left(\mathit{x}\right)$ have to be zero as a necessary condition if $V\left(\mathit{x}\right)=0$ holds constantly, the invariant limit set of the system can be found.

This statement can be verified by means of the following relations:

$$\begin{array}{cc}\hfill V\left(\mathit{x}\right)=\frac{1}{2}{s}^2\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill \dot{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\hfill \end{array}$$

(24)

satisfied due to (23)

$$\begin{array}{cc}\hfill \ddot{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\ddot{s}\left(\mathit{x}\right)+{\dot{s}}^2\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\hfill \end{array}$$

(25)

satisfied due to (23) ∧$\dot{s}\left(\mathit{x}\right)$ = 0

$$\begin{array}{cc}\hfill \stackrel{⃛}{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\stackrel{⃛}{s}\left(\mathit{x}\right)+3\dot{s}\left(\mathit{x}\right)·\ddot{s}\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\wedge \dot{s}\left(\mathit{x}\right)=0.\hfill \end{array}$$

(26)

satisfied due to (23) and (25).

Requiring furthermore that $\ddot{s}\left(\mathit{x}\right)=0$ and $\stackrel{⃛}{s}\left(\mathit{x}\right)=0$ hold, all further higher order derivatives of $V\left(\mathit{x}\right)$ up to order seven are directly proven to be equal to zero.

Moreover, accounting for the fact that $\dot{s}\left(\mathit{x}\right)=-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)$ holds by construction (see Equation (7)), the relations (23)–(25) above turn into:

$$\begin{array}{cc}\hfill V\left(\mathit{x}\right)=\frac{1}{2}{s}^2\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \dot{V}\left(\mathit{x}\right)=-s^2\left(\mathit{x}\right)·p\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0\hfill \end{array}$$

(28)

satisfied due to (27)

$$\begin{array}{cc}\hfill \ddot{V}\left(\mathit{x}\right)=2s^2\left(\mathit{x}\right)·p^2\left(\mathit{x}\right)-s^2\left(\mathit{x}\right)·\dot{p}\left(\mathit{x}\right)=0& ⟹s\left(\mathit{x}\right)=0,\hfill \end{array}$$

(29)

satisfied due to (27), where $\ddot{V}\left(\mathit{x}\right)$ in (29) can be rewritten into:

$$\ddot{V}\left(\mathit{x}\right)=s^2\left(\mathit{x}\right)·\left(2p^2\left(\mathit{x}\right)-\dot{p}\left(\mathit{x}\right)\right).$$

(30)

Hence, the equivalent control designed in this paper guarantees that all temporal derivatives of $V\left(\mathit{x}\right)$ are perfectly equal to zero as soon as $s\left(\mathit{x}\right)=0$ has been reached. Enforcing furthermore that $\dot{s}\left(\mathit{x}\right)=\ddot{s}\left(\mathit{x}\right)=\stackrel{⃛}{s}\left(\mathit{x}\right)=0$ (involved in (23)–(26)) holds after expressing these terms as functions of the controlled system dynamics on the invariant set during the interval based solution procedure in Section 5 helps to improve the computational efficiency by means of defining so-called contractors that quickly shrink the search space for the invariant set towards the actual solution domain.

Remark 2.

If an invariant set is computed according to the definition:

$$S_1=\left\{\mathit{x}\right|\dot{V}\left(\mathit{x}\right)=\ddot{V}\left(\mathit{x}\right)=\stackrel{⃛}{V}\left(\mathit{x}\right)=0\},$$

(31)

it is not guaranteed that the system converges under all circumstances to the set:

$$S_0=\left\{\mathit{x}\right|V\left(\mathit{x}\right)=\dot{V}\left(\mathit{x}\right)=\ddot{V}\left(\mathit{x}\right)=\stackrel{⃛}{V}\left(\mathit{x}\right)=0\},$$

(32)

where the system is on the sliding surface. If the condition $s\left(\mathit{x}\right)=0$ is not enforced explicitly in the interval based procedure described in the following section, the set $S_1$ may not only consist of the invariant set that is characterized by the sliding mode, but may also contain further fixed points.

For example, for systems of dimension two, in which the examined sliding surface describes a closed trajectory in the state space in terms of a limit cycle, the Poincaré index theorem can be considered to visualize this fact. This theorem states that every limit cycle encloses at least one equilibrium point of the system, which will therefore also belong to the invariant set defined by $S_1$ [40].

5. Interval Computation

The outer interval enclosure of the invariant set $S_1$ of the controlled system in $\mathit{x}$, denoted exemplarily by $(x_1,x_2)$, is calculated using interval computation. This is possible due to the fact that the Lyapunov-like functions defined before provide inequalities. Therefore, it is checked for subintervals (i.e, boxes) covering the operating domain of interest if the conditions:

$$\begin{array}{cc}\hfill s\left(\mathit{x}\right)& =0(\mathrm{which}\mathrm{is}\mathrm{not}\mathrm{included}\mathrm{in}\mathrm{the}\mathrm{definition}\mathrm{of}{S}_1,\mathrm{but}\mathrm{only}\mathrm{in}{S}_0),\hfill \\ \hfill \dot{s}\left(\mathit{x}\right)& =0,\hfill \\ \hfill \ddot{s}\left(\mathit{x}\right)& =0,\mathrm{and}\hfill \\ \hfill \stackrel{⃛}{s}\left(\mathit{x}\right)& =0\hfill \end{array}$$

(33)

are satisfied in a domain that is a priori restricted to the set $s^2\left(\mathit{x}\right)\le \mu ,\mu >0$. This last condition restricts the search space to a maximum admissible squared distance from the sliding surface.

To solve a problem according to the specification (33), a general approach for interval techniques aiming at set inversion should be briefly described. The following description is based on Figure 2. In the examined case, a set is known in the domain of $(y_1,y_2)$, and the set that should be found by the inversion technique is located in the $(x_1,x_2)$-plane. The function $\mathit{f}$ performing the mapping, $\mathit{y}=\mathit{f}\left(\mathit{x}\right),\mathit{x}\in {\mathbb{R}}^2,\mathit{y}\in {\mathbb{R}}^2$, from the $(x_1,x_2)$-plane to the $(y_1,y_2)$-plane is assumed to be known together with a suitable inclusion function $\left[\mathit{f}\right]$ that finds a rectangular box that encloses the set (gray). This function will not necessarily lead to the tightest enclosing box, but it is convergent with respect to interval bisections. The inverse mapping of $\mathit{f}$ does not necessarily have to be known explicitly, so $(x_1,x_2)$ either does not need to be or even cannot be determined analytically.

Starting with a sufficiently large box $\left[{\mathit{x}}_0\right]$ in the $(x_1,x_2)$-plane that encloses the whole solution domain of interest, the corresponding box in $(y_1,y_2)$ is calculated. It is then checked if it is located either outside the set (blue), inside the set (red), or if it is overlapping with both, the outer and interior domains, i.e., if the location is undecided (yellow). Whenever it is not sure if the box is inside or outside, the box in the $(x_1,x_2)$-plane will be split at the midpoint of each edge, and the procedure will be repeated for each of the new boxes. This is continued until a certain accuracy is reached. The red boxes will represent the corresponding set in the $(x_1,x_2)$-plane. This algorithm is called set inversion via interval analysis (SIVIA) [2].

Note that interval computations are guaranteed in the following sense: If a box in the $(x_1,x_2)$-plane has been found that is mapped into the interior of $\mathcal{Y}$, it is guaranteed to be part of the true solution set $\mathcal{Y}$, and vice versa, that if a box lies outside the solution set $\mathcal{Y}$ in the $(y_1,y_2)$-plane, its counterpart in $(x_1,x_2)$ certainly contains no admissible solutions. If both cases are included in one box, it will, as mentioned before, either be split or, if the volume of the respective hypercuboid falls below some threshold value, be labeled as undecided.

As this successive bisectioning scheme leads to lots of boxes that have to be checked ($2^n$ boxes for a simultaneous subdivision of an n-dimensional box in each of its edges) the concept of contractors can effectively be applied to shrink a box as far as possible to the solution set. Here, the interval contraction is performed by means of each component of the vector valued function $\mathit{f}$.

A contractor can be built by looking at the syntactic tree of the function $\mathit{f}$, which means a decomposition into primitive equations by defining auxiliary variables [2]. In this way, the intervals can be contracted in a forward-backward manner. The forward part is where the definitions of the auxiliary variables are executed, and the backward part is where the original intervals are intersected with the newly calculated ones through the auxiliary variables’ definition.

This type of algorithm is described in more detail in [2] and has been employed for an interval based parameter identification scheme in [19]. Using this contractor, the first box will be blue, and inside this box, there will be a yellow one that is the smallest one to contain the whole solution set included in the first blue box, if the contractor is minimal. Then, the yellow box containing the solution set is split, and the procedure is repeated for each of the boxes. In this way, the boxes become smaller more quickly than in the case where every box is repeatedly divided. Figure 3 illustrates this fact [2].

To transfer this contractor approach to the problem defined in this paper according to Equation (33), the following analogues are exploited: (i) the first condition $\dot{\tilde{V}}=0$ corresponds to $\mathit{f}$; (ii) $(x_1,x_2)$ is generalized in terms of the n-dimensional state vector; (iii) the value of zero is considered as the known set in the image space $\mathit{y}$. Analogously, contractors are also built for the further derivatives in (33). All of these contractors are intersected with each other because all conditions should be satisfied simultaneously to provide an (outer) interval enclosure of the system’s invariant set.

6. Application Scenarios

6.1. Pendulum Motion on a Circle in the Phase Plane

A pendulum with velocity proportional linear friction, as well as a nonlinear restoring sine-type force shall be described by the following state equations:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =x_2\hfill \\ \hfill {\dot{x}}_2& =-sin\left(x_1\right)-x_2+u,\hfill \end{array}$$

(34)

where $x_1$ is the angle, $x_2$ the angular velocity, and u is the control torque. A controller $u\left(\mathit{x}\right)$ should be found to make the pendulum move on a circle in the phase plane. This means, according to the sliding mode control design presented in Section 3, the switching function $s\left(\mathit{x}\right)$ is defined as:

$$s\left(\mathit{x}\right)=x_1^2+x_2^2-1$$

(35)

describing a circle with a radius of one around the origin.

6.1.1. Controller Design

From Equations (12) and (35), the control input can be found as:

$$u=\left(-\left(x_1^2+x_2^2-1\right)·p\left(\mathit{x}\right)-2x_1{x}_2\right)·\frac{1}{2x_2}+sin\left(x_1\right)+x_2.$$

(36)

The Bézout identity to solve is hence:

$$\left(x_1^2+x_2^2-1\right)·p\left(\mathit{x}\right)+2x_{2}k\left(\mathit{x}\right)=-2x_1{x}_2.$$

(37)

The reduced Gröbner bases for both ideals:

$$I_1=〈x_1^2+x_2^2-1,2x_2〉\mathrm{and}{I}_2=〈x_1^2+x_2^2-1,2x_2,-2x_1{x}_2〉$$

(38)

are the same as:

$$G_1=G_2=\left\{x_1^2-1,x_2\right\}.$$

(39)

This means there is a solution for Equation (37). Choose, e.g., $p\left(\mathit{x}\right)=4x_2^2={\left({\displaystyle \frac{\partial s}{\partial x_2}}\right)}^2$, which makes:

$$k\left(\mathit{x}\right)=-2x_1^2{x}_2-x_1-2x_2^3+2x_2$$

(40)

and, finally, leads to the nonlinear equivalent control:

$$u=-2x_2·\left(x_1^2+x_2^2-1\right)-x_1+sin\left(x_1\right)+x_2.$$

(41)

6.1.2. Verification of Stability

The stability of the closed-loop control structure can be verified by computing the total time derivative of the sliding surface $s\left(\mathit{x}\right)$ defined in (35) according to:

$$\begin{array}{cc}\hfill \dot{s}\left(\mathit{x}\right)& =\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}·{\dot{x}}_1+\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}·{\dot{x}}_2\hfill \\ \hfill & =2x_1{\dot{x}}_1+2x_2{\dot{x}}_2\hfill \\ \hfill & =2x_1{x}_2+2x_2·\left(-sin\left(x_1\right)-x_2+\left(-\left(x_1^2+x_2^2-1\right)·p\left(\mathit{x}\right)-2x_1{x}_2\right)·\frac{1}{2x_2}+sin\left(x_1\right)+x_2\right)\hfill \\ \hfill & =-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)\mathrm{with}p\left(\mathit{x}\right)=4x_2^2>0\mathrm{for}s\left(\mathit{x}\right)\ne 0\wedge x_2\ne 0\hfill \end{array}$$

(42)

along the trajectories of the controlled system. Hence, $\dot{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)<0$ is guaranteed for all states that do not belong to the desired invariant set, except for the fixed point at the origin of the phase space.

Remark 3.

The control law defined in (41) ensures perfect trajectory tracking as soon as the system has reached the sliding surface $s\left(\mathit{x}\right)=0$. However, it does not possess the variable structure characteristics defined in (5). Therefore, typically, it only leads to an asymptotic convergence towards the desired invariant set. However, here, and analogously in all other examples shown in the remainder of the paper, finite time convergence towards $s\left(\mathit{x}\right)=0$ can be achieved if u is replaced by:

$$\tilde{u}=u-\frac{\eta }{2x_2}·\mathrm{sign}\left(s\left(\mathit{x}\right)\right)\mathrm{with}\eta >0$$

(43)

for states that are not contained in the invariant set. Then, the control $\tilde{u}$ turns into the equivalent control u as soon as $s\left(\mathit{x}\right)=0$ starts to hold, and hence, the computation procedure for invariant sets derived in Section 4 becomes again applicable.

6.1.3. Computation of the Invariant Set

The invariant set $S_1$ of the closed-loop system, where (33) is true, can be calculated via interval computation (see Figure 4 in which the same color code as in Figure 3 is employed).

The vector field in Figure 5 supports this result. As is shown in these figures, the pendulum converges to the circle or the equilibrium point at the origin, that is inside the circle due to the Poincaré fixed point theorem.

6.2. Pendulum Motion on an Ellipse in the Phase Plane

The pendulum described by Equation (34) will now move on an ellipse in the phase plane. Thus, the switching function (35) is replaced by:

$$s\left(\mathit{x}\right)=x_1^2+x_2^2+x_1{x}_2-1.$$

(44)

Controller Design

The Bézout identity to be solved for this case is:

$$\left(2x_2+x_1\right)·k\left(\mathit{x}\right)+\left(x_1^2+x_2^2+x_1{x}_2-1\right)·p\left(\mathit{x}\right)=-\left(2x_1+x_2\right)·x_2.$$

(45)

As the reduced Gröbner bases:

$$G_1=\left\{x_1+2x_2,3x_2^2-1\right\}\mathrm{and}{G}_2=\left\{1\right\}$$

(46)

are not the same, there is no solution for an input without singularity, and thus, the pendulum cannot be made to move on this specific ellipse.

6.3. Pendulum with an Iso-Energy Constraint

Now, the pendulum should keep a constant energy. The switching function is then defined as a combination of kinetic and potential energy:

$$E\left(\mathit{x}\right)=x_2^2-cos\left(x_1\right)+1=2$$

(47)

leading to:

$$s\left(\mathit{x}\right)=x_2^2-cos\left(x_1\right)-1.$$

(48)

6.3.1. Controller Design

Due to Equation (12), the controller is now parameterized by:

$$u=\frac{1}{2x_2}\left(-\left(x_2^2-cos\left(x_1\right)-1\right)·p\left(\mathit{x}\right)-sin\left(x_1\right)·x_2\right)+sin\left(x_1\right)+x_2.$$

(49)

As this is not a polynomial equation, the Bézout identity is not of any help. However, luckily, due to the surface definition, no singularity occurs, if $p\left(\mathit{x}\right)$ is chosen as $p\left(\mathit{x}\right)=2x_2^2$.

6.3.2. Verification of Stability

The stability of the closed-loop control structure can be verified by computing the total time derivative of the sliding surface $s\left(\mathit{x}\right)$ defined in (48) according to:

$$\begin{array}{cc}\hfill \dot{s}\left(\mathit{x}\right)& =\frac{\partial s\left(\mathit{x}\right)}{\partial x_1}·{\dot{x}}_1+\frac{\partial s\left(\mathit{x}\right)}{\partial x_2}·{\dot{x}}_2\hfill \\ \hfill & =sin\left(x_1\right)·{\dot{x}}_1+2x_2{\dot{x}}_2\hfill \\ \hfill & =sin\left(x_1\right)·x_2+2x_2·(-sin\left(x_1\right)-x_2\hfill \\ \hfill & \phantom{sin\left(x_1\right)·x_2+2x_2·\left(\right)}+\frac{1}{2x_2}\left(-\left(x_2^2-cos\left(x_1\right)-1\right)·p\left(\mathit{x}\right)-sin\left(x_1\right)·x_2\right)+sin\left(x_1\right)+x_2)\hfill \\ \hfill & =-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)\mathrm{with}p\left(\mathit{x}\right)=2x_2^2>0\mathrm{for}s\left(\mathit{x}\right)\ne 0\wedge x_2\ne 0\hfill \end{array}$$

(50)

along the trajectories of the controlled system. Hence, $\dot{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)<0$ is guaranteed for all states that do not belong to the desired invariant set, except for the fixed point at the origin of the phase space.

6.3.3. Computation of the Invariant Set

The invariant set $S_1$ for the iso-energy constraint of the controlled pendulum is depicted in Figure 6.

This corresponds to the pendulum turning all around, and from the top $x_1=\pm \pi$, it is possible to turn to either of the sides. Similarly, as in Figure 4, a fixed point at the origin of the phase space was verified by means of the interval evaluation in this scenario.

6.4. Stabilization of the Pendulum Motion towards Zero

In this case, the pendulum should be brought to its naturally stable equilibrium, where:

$$x_1=0\wedge x_2=0.$$

(51)

Therefore, the switching function is given as:

$$s\left(\mathit{x}\right)=x_1+x_2,$$

(52)

which corresponds to a straight line passing through the origin. This is a typical choice of s to bring the system to the origin.

6.4.1. Controller Design

Due to Equation (12), the control input has to be chosen as:

$$u=\left(-(x_1+x_2)·p\left(\mathit{x}\right)-x_2\right)+sin\left(x_1\right)+x_2.$$

(53)

Then, the Bézout identity trivially yields:

$$(x_1+x_2)·p\left(\mathit{x}\right)+k\left(\mathit{x}\right)=-x_2.$$

(54)

Choose $p\left(\mathit{x}\right)=1$, such that the input:

$$u=sin\left(x_1\right)-x_2-x_1$$

(55)

is obtained.

6.4.2. Verification of Stability

The stability of the closed-loop control structure can be verified by computing the total time derivative of the sliding surface $s\left(\mathit{x}\right)$ defined in (52) according to:

$$\begin{array}{cc}\hfill \dot{s}\left(\mathit{x}\right)& ={\dot{x}}_1+{\dot{x}}_2\hfill \\ \hfill & =x_2-sin\left(x_1\right)-x_2+sin\left(x_1\right)-x_2-x_1=-s\left(\mathit{x}\right)\hfill \end{array}$$

(56)

along the trajectories of the controlled system. Hence, $\dot{V}\left(\mathit{x}\right)=-s^2\left(\mathit{x}\right)<0$ is guaranteed for all states that do not belong to the desired invariant set.

6.4.3. Computation of the Invariant Set

The invariant set for this case is depicted in Figure 7. As can be seen, the whole line is part of the invariant set, implying that once on the surface, the system will stay there for ever. To see that the pendulum will actually converge to its equilibrium point, the vector field in Figure 8 is helpful. It can be seen that once on the surface, the trajectories tend to the origin.

6.5. Line Following for a Dubins Car

As a final example, consider the Dubins car [45], which is described by the state equations:

$$\begin{array}{cc}\hfill {\dot{x}}_1& =cos\left(x_3\right)\hfill \\ \hfill {\dot{x}}_2& =sin\left(x_3\right)\hfill \\ \hfill {\dot{x}}_3& =u,\hfill \end{array}$$

(57)

where $x_3$ is the heading of the car and $(x_1,x_2)$ are the position coordinates of its center. It should follow a line $x_2=0$, which is defined by the switching function:

$$s\left(\mathit{x}\right)=x_3+arctan\left(x_2\right),$$

(58)

with the heading being zero, if the car is on $x_2=0$ and $arctan\left(x_2\right)=0$ for $x_2=0$. The saturation function arctan is used to exclude the (unbounded) case $x_2=-x_3\ne 0$ from the surface.

6.5.1. Controller Design

Due to Equation (7), the system’s input in terms of a nonlinear feedback controller becomes:

$$u=-\left(x_3+arctan\left(x_2\right)\right)·p\left(\mathit{x}\right)-\frac{sin\left(x_3\right)}{1+x_2^2}.$$

(59)

6.5.2. Verification of Stability

Stability of the closed-loop control structure can be verified by computing the total time derivative of the sliding surface $s\left(\mathit{x}\right)$ defined in (58) according to:

$$\begin{array}{cc}\hfill \dot{s}\left(\mathit{x}\right)& =\frac{{\dot{x}}_2}{1+x_2^2}+{\dot{x}}_3\hfill \\ \hfill & =\frac{sin\left(x_3\right)}{1+x_2^2}-\left(x_3+arctan\left(x_2\right)\right)·p\left(\mathit{x}\right)-\frac{sin\left(x_3\right)}{1+x_2^2}\hfill \\ \hfill & =-s\left(\mathit{x}\right)·p\left(\mathit{x}\right)\mathrm{with}p\left(\mathit{x}\right)>0\mathrm{for}s\left(\mathit{x}\right)\ne 0\hfill \end{array}$$

(60)

along the trajectories of the controlled system. Hence, $\dot{V}\left(\mathit{x}\right)=s\left(\mathit{x}\right)·\dot{s}\left(\mathit{x}\right)<0$ is guaranteed for all states that do not belong to the desired invariant set.

6.5.3. Consideration of Input Range Constraints

As a new aspect, range constraints for this input should be examined, as it is not a polynomial equation and thus not Bézout’s identity. Let $u\in [-1;1]$. This makes:

$$p\left(\mathit{x}\right)\in \frac{sin\left(x_3\right)}{-(x_3+arctan\left(x_2\right))·(1+x_2^2)}+\frac{1}{x_3+arctan\left(x_2\right)}·[-1;1].$$

(61)

The extremal choice for $p\left(\mathit{x}\right)$ is:

$$p\left(\mathit{x}\right)=\frac{sin\left(x_3\right)}{-(x_3+arctan\left(x_2\right))·(1+x_2^2)}+\frac{1}{|x_3+arctan\left(x_2\right)|}>0.$$

(62)

The extremal choice of $p\left(\mathit{x}\right)$ leads to the largest absolute value of u. If this is inserted in u, the result is:

$$\begin{array}{cc}\hfill u& =-(x_3+arctan\left(x_2\right))·\left(\frac{sin\left(x_3\right)}{-(x_3+arctan\left(x_2\right))·(1+x_2^2)}+\frac{1}{|x_3+arctan\left(x_2\right)|}\right)-\frac{sin\left(x_3\right)}{1+x_2^2}\hfill \\ \hfill & =\frac{-(x_3+arctan\left(x_2\right))}{|x_3+arctan\left(x_2\right)|}\hfill \\ \hfill & =-\mathrm{sign}(x_3+arctan\left(x_2\right)).\hfill \end{array}$$

(63)

To make the input differentiable, the hard switching function is replaced by:

$$u=-tanh(x_3+arctan\left(x_2\right)).$$

(64)

6.5.4. Computation of the Invariant Set

As this example is three dimensional, the result is projected into the $(x_1,x_2)$-plane. As can be seen in Figure 9, the invariant set is the line, where $x_2=0$ as desired.

The control input (64) stays within its required bounds for a simulation with randomly chosen initial conditions $x\left(0\right)$, as can be seen in Figure 10.

7. Conclusions and Outlook on Future Work

In this paper, a method to compute the invariant set of a nonlinear system controlled by sliding mode approaches, which were designed in such a way that structural singularities were guaranteed to be prevented, was suggested. Numerically, the invariant set could be obtained using interval computation, if the controller was designed according to the equation:

$$\dot{s}\left(\mathit{x}\right)=-s\left(\mathit{x}\right)·p\left(\mathit{x}\right),p\left(\mathit{x}\right)>0.$$

(65)

Avoiding singularities in u was ensured using a suitable Gröbner basis representation. The various investigated examples showed a surface represented by a closed curve, a case where Bézout’s identity had no solution and thus no controller can be found, and a surface that was a line. Additionally, a three-dimensional example was given and bounds, on the control input were considered.

The presented interval technique for the identification of invariant sets can easily be generalized towards the computation of the region of attraction of asymptotically stable equilibria of nonlinear feedback control systems by verifying those domains in the state space in which the time derivative $\dot{V}\left(\mathit{x}\right)$ of a Lyapunov function candidate was strictly non-positive within a level set characterized by $V\left(\mathit{x}\right)=c>0$, where this curve needed to be proven by the suggested approach not to coincide with an invariant set. Straightforwardly, also generalizations to the identification of domains in the state space, for which the stability of stochastic differential equations cannot be verified, became possible. Such domains were derived in [48].

In future work, the suggested procedure for the computation of invariant sets can be extended by the consideration of set valued uncertainty in the state variables included in the control signal $u\left(\mathit{x}\right)$, which result either from an imperfect state measurement or from uncertainty if the state information included in the feedback control law is reconstructed by means of an interval observer [49]. Moreover, procedures for enhanced control parameterizations can be developed on the basis of the presented identification of invariant sets so that the volume of the intervals containing these sets is minimized. Hence, the aim will be to make the guaranteed provable stability domains as large as possible, as well as to maximize the regions of attraction of asymptotically stable equilibria by control re-parameterizations, to minimize the non-provable stability domains in the aforementioned stochastic settings, and to tune Lyapunov functions so that the provable stability domains show a minimum amount of pessimism.