1 Introduction↩︎

Time constraints are relevant in both controller and observer performances. For control, the settling time is an important specification that enables faster transient responses and more robustness to disturbances, which can be motivated for safety or productivity reasons. For observation, faster convergence can be critical to using the estimated states in real-time applications. Therefore, time constraints are present in many engineering problems as shown in [1] and references therein.

When dealing with precisely known asymptotically stable linear systems, with no regards to state constraints, the time-optimal control problem is known to be solved by employing the “bang-bang principle” [2]. In more general settings, the Kruzkov Transformation can be employed to transform the Hamilton-Jacobi-Bellman (HJB) equations associated with the time-optimal control problem into similar HJB equations as discounted optimal control problems, which can be solved numerically [3], [4]. However, this approach can be trickier when the system under consideration does not satisfy certain regularity conditions, usually small-time controllability of the target set and Soner’s inward pointing condition, since the viscosity solution of the HJB can be discontinuous [5].

Instead of directly solving the time-optimal control problem, finite-time stability [6] defines a set of Lyapunov conditions that, if satisfied, ensures the origin of the state space is reached in a finite amount of time (in contrast to asymptotic stability which only requires the origin to be reached as \(t \rightarrow \infty\)). Even though these are useful analysis conditions, they are often hard to employ in a constructive manner for control laws. As such, it is common to employ Implicit Lyapunov Function (ILF) designs [7], or Prescribed Performance Control (PPC) [8] in order to impose finite-time stability on the closed loop, with the main disadvantage of these approaches being the fact that they require a specific structure upon the system’s dynamics.

Many approaches disregard the fact that most applications have constraints upon the system’s state and input, as well as uncertainties on the available model. Throughout the last 30 years, robust control conditions based on Linear Matrix Inequalities (LMIs) [9] have been shown to be particularly useful in this regard, specially if the uncertainties can be written as a convex set (leading to a polytopic representation). An approach that is already well-developed and employed in the literature is to make use of a saturation function in the control law, and ensuring the closed loop’s stability inside of a predetermined level set of a Lyapunov function [10].

Even though the existence of a quadratic Lyapunov function is necessary and sufficient for the stability of precisely known linear time-invariant systems, it is a well-known fact that it is a source of conservativeness when dealing with uncertain [11] or nonlinear systems [12]. An interesting approach in that regard is the use of Piecewise Quadratic Lyapunov functions [13], which, given enough partitions of the state space, is capable of approximating any nonlinear function. As shown by [14], they can be employed to approximate upper-bounds of optimal control problems, though this approach is only valid for the largest Lyapunov level set contained in the analysis region (which can be hard to find in this scenario).

Building on this idea, [15] propose a relaxed optimal control problem (for polynomial systems), and a Policy Iteration solution based on sum of squares (SOS), which is capable of finding a global suboptimal solution for the relaxed optimal control problem given an initial policy/control law. Its main disadvantage lies on the fact that it does not take into account state and input constraints.

Inspired by these works, this paper proposes a Guaranteed Time Control approach, in which the goal is to minimize an upper-bound of the time taken to reach a desired region that covers the origin of the state space, while respecting bounds on the state space and control inputs and ensuring asymptotic stability of the origin. In order to do so, we formulate a theorem with Lyapunov conditions based on the Harmonic Transformation [16], which, if satisfied, not only provides an upper-bound on the time taken to reach the desired region, but also an estimate of the origin’s domain of attraction. In order to make use of these new conditions, we propose a novel representation for piecewise quadratic Lyapunov functions, based on a simplicial partition of the state space, and a Policy Iteration approach based on Linear Matrix Inequalities relying on the structural relaxation approach [17], [18] to introduce local information into the conditions. A remark details how the proposed conditions can be easily modified to deal with piecewise polytopic systems, and nonlinear systems represented by Takagi-Sugeno fuzzy models [19]. Three examples are presented to illustrate the proposed approach and compare it against the literature.

Notation↩︎

We represent matrix, vectors and scalar variables by uppercase, bold text lower case, and lowercase letters respectively. Though an abuse of notation, some scalar functions will be represented by uppercase variables, either to avoid confusion, or due to that being the standard use in the literature. Given a set \(\mathcal{A}\), its boundary is represented by \(\partial\mathcal{A}\). Given matrix \(M \in \mathbb{R}^{r \times p}\), \(\boldsymbol{m}_j\) represents its \(j\)-th column, and \(m_{ij}\) represents the element in row \(i\) and column \(j\). For matrices, \(S \succeq 0\) (\(\preceq 0\)) describes that the elements of the matrix are non-negative (non-positive). For symmetric matrices, \(M > 0\) (\(< 0\)) describes that matrix \(M\) is positive (negative) definite, and \(M \geq 0\) (\(\leq 0\)) is positive (negative) semidefinite. \(M^T\) represents the transpose of a matrix, \(M^{-1}\) its inverse, and \(M^{-T} = \left(M^T\right)^{-1}=\left(M^{-1}\right)^T\) the inverse of the transpose. When describing a symmetric matrix, the symbol \(\ast\) represents block transpose terms that can be inferred from the symmetry. Throughout this paper, different convex sums will be employed in different contexts (with different convex weights), and the short notation \(S_\alpha = \sum_{i = 1}^r \alpha_i S_i\) is employed. When dealing with a set of indices up to \(s\) we will make use of the representation \(\mathcal{I}_s = \{i \in \mathbb{N} | 1 \leq i \leq s\}\). Time dependency will usually be dropped, and will only be shown explicitly when deemed necessary.

2 Methodology↩︎

Consider an uncertain polytopic system described by \[\begin{align} \dot{\boldsymbol{x}} = \sum_{i = 1}^r \alpha_i \left(A_i \boldsymbol{x} + B_i \boldsymbol{u} \right) = A_\alpha \boldsymbol{x} + B_\alpha \boldsymbol{u}, \label{eq:sys} \end{align}\tag{1}\] with \(\boldsymbol{x} \in \mathcal{X} \subset \mathbb{R}^n\) and \(\boldsymbol{u} \in \mathcal{U} \subset \mathbb{R}^m\) representing the states and the control inputs from bounded sets \(\mathcal{X}\) and \(\mathcal{U} = \{\boldsymbol{u} \in \mathbb{R}^m | \; \underline{u}_i \leq u_i \leq \bar{u}_i, \forall i = \mathcal{I}_m\}\), respectively; \(r\) the number of vertices in the polytopic representation, \(A_i \in \mathbb{R}^{n \times n}, B_i \in \mathbb{R}^{n \times m}\) the matrices representing the system’s dynamics, and \(\boldsymbol{\alpha}^T = \begin{bmatrix} \alpha_1 & \dots & \alpha_r \end{bmatrix}\) the uncertain polytopic convex weights satisfying \(\alpha_i \geq 0, \forall i\) and \(\sum_{i=1}^r \alpha_i = 1\).

Our goal is to find a control law capable of driving the largest possible subset of \(\mathcal{X}\) towards a set \(\mathcal{X}_g\), that contains the origin, in the shortest time and without violating the bounds imposed by \(\mathcal{U}\). To that end, consider the following theorem.

Theorem 1 is an analysis tool capable of estimating both a domain of attraction as well as the time to reach the set \(\mathcal{X}_g\), and any function satisfying its conditions would give us a time upper-bound. In trying to find the smallest upper-bound, we can use the conditions of the theorem as constraints and the integral of the \(\bar{V}(\boldsymbol{x})\) function over the state space (\(\int_\mathcal{X} \bar{V}(\boldsymbol{x}) \boldsymbol{dx}\)) as the cost in an optimization problem.

Similar to other approaches in the literature [15], [21], [22], employing optimization algorithms to solve for the conditions in Theorem 1 can be seen as a sort of Approximate Dynamic Programming (ADP), since the upper-bound solutions that satisfy the conditions in the theorem can be seen as supersolutions to the corresponding harmonically transformed Hamilton-Jacobi-Bellman equations for minimum-time problems [16], and solving for the minimum upper-bound can then be seen as trying to solve the Relaxed Optimal Control Problem [15].

Many of the optimization-based ADP approaches in the literature focus their attention on the global solution, considering an unrestricted control set and quadratic costs in the control and states (leading to a closed form solution for the control law/policy given a certain Lyapunov/value function). They usually consider a polynomial representation, such that the optimization problem can be cast as a sum of squares (SOS) optimization problem [15], [21], [22].

2.1 Piecewise Representation↩︎

Unlike previous approaches in the literature, in this work we focus our attention to a local constrained (both in states and control) setting. To accomplish this, we consider a piecewise quadratic representation for \(\bar{V}(\boldsymbol{x})\) [13], [14], which is more directly amenable to deriving LMI conditions, and recasting our optimization as a Semi-Definite Programming Problem. However, similarly to [23], we do not make use of the parametrization defined by the constraint matrices used in [13], and present a new parametrization for piecewise quadratic functions.

Consider that we have a simplicial partition of the set \(\mathcal{X}\), \[\mathcal{X} = \bigcup_{q = 1}^{r_s} \boldsymbol{\Delta}_q, \label{eq:union95simplex}\tag{2}\] representing our local analysis region in the state space with \(\boldsymbol{\Delta}_q\) representing each simplex in our partition, and \(r_s\) the total number of simplices. This partition can always be found from a (possibly nonregular) grid of points \(\Gamma = \{\boldsymbol{x}_1, \dots , \boldsymbol{x}_{r_p}\}\) in the state space, with \(r_p\) the number of grid points, by using a Delaunay Triangulation [24]. We consider that the origin is always part of this generating grid, and that the set \(\mathcal{X}_g\) is composed by the simplices which contain the origin as a vertex. This description is illustrated in Figure 1 for a 2-dimensional state space.

In addition to this, the triangulation of the state space defines a mapping \(\eta(q,j):\mathcal{I}_{r_s} \times \mathcal{I}_{n+1} \rightarrow \mathcal{I}_{r_p}\), which given indices \(q\) and \(i\) representing a simplex \(\boldsymbol{\Delta}_q\) and one of the vertices of said simplex respectively, relates this vertex to a point on the grid \(\Gamma\). Given this mapping, each simplex is defined as the convex hull of a set of \(n+1\) vertex points, such that \[\begin{align} \boldsymbol{\Delta}_q &= \operatorname{co}\left(\boldsymbol{x}_{\eta(q,1)}, \dots, \boldsymbol{x}_{\eta(q,n+1)}\right), \end{align}\] which implies that there exists \(\boldsymbol{\beta} \in \mathbb{R}^{n+1}\) unique convex weights that can be used to represent any point in the simplex from its vertices as \(\boldsymbol{x} = \sum_{j = 1}^{n+1} \beta_j \boldsymbol{x}_{\eta(q,j)}\), with \(\boldsymbol{\beta}^T = \begin{bmatrix} \beta_1 & \dots & \beta_{n+1}\end{bmatrix}\), \(\beta_{j} \geq 0 \; \forall j, \; \sum_{j=1}^{n+1} \beta_j = 1\).

For each point \(\boldsymbol{x}_i\) in our grid, we define an affine function \[\bar{v}_i(\boldsymbol{x}) = \boldsymbol{p}_{i}^T \begin{bmatrix} \boldsymbol{x} \\ 1 \end{bmatrix} = \boldsymbol{p}_i^T \bar{\boldsymbol{x}},\] with \(\bar{\boldsymbol{x}}\) representing an augmented state vector, \(\boldsymbol{p}_i \in \mathbb{R}^{n+1}\) the parameters defining the function, and \(\boldsymbol{p}_{i(n+1)} = 0\) if \(\boldsymbol{x}_i = \boldsymbol{0}\), such that \(\bar{v}_i(\boldsymbol{0}) = 0\). Our Piecewise Quadratic function is given by the finite-element linear interpolation of these points, using our simplices as the finite-elements. Note that, for any point belonging to simplex \(\boldsymbol{\Delta}_q\), we can write that \[\begin{align} \begin{bmatrix} \boldsymbol{x} \\ 1 \end{bmatrix} &= \begin{bmatrix} \boldsymbol{x}_{\eta(q,1)} & \dots & \boldsymbol{x}_{\eta(q,n+1)} \\ 1 & \dots & 1 \end{bmatrix} \begin{bmatrix} \beta_1 \\ \vdots \\ \beta_{n+1}\end{bmatrix},\\ \bar{\boldsymbol{x}} &= \bar{X}_q \boldsymbol{\beta}, \end{align}\] which allow us to calculate the convex weights by \[\begin{align} \boldsymbol{\beta} &= \bar{X}_q^{-1} \bar{\boldsymbol{x}}. \end{align}\] Note that, so long as there are no degenerate simplices in the triangulation (simplices whose hypervolume is zero), the inverse is always guaranteed to exist (since the hypervolume of simplex \(\boldsymbol{\Delta}_q\) is proportional to \(\det(\bar{X}_q)\)).

Inside of simplex \(\boldsymbol{\Delta}_q\), the finite-element linear interpolation gives us \[\begin{align} V_q(\boldsymbol{x}) &= \sum_{j=1}^{n+1} \beta_j \boldsymbol{p}_j^T \bar{\boldsymbol{x}} = \begin{bmatrix} \beta_1 & \dots & \beta_{n+1} \end{bmatrix} \begin{bmatrix} \boldsymbol{p}_1^T \\ \vdots \\ \boldsymbol{p}_{n+1}^T \end{bmatrix} \bar{\boldsymbol{x}}, \\ &= \boldsymbol{\beta}^T P_q^T \bar{\boldsymbol{x}} = \bar{\boldsymbol{x}}^T \bar{X}_q^{-T} P_q \bar{\boldsymbol{x}}, \\&= \bar{\boldsymbol{x}}^T \left(\dfrac{1}{2}\left(P_q^T \bar{X}_q^{-1} + \bar{X}_q^{-T} P_q\right)\right)\bar{\boldsymbol{x}}, \\ V_q(\boldsymbol{x}) &= \bar{\boldsymbol{x}}^T S_q \bar{\boldsymbol{x}}, \end{align}\] with \(P_q \in \mathbb{R}^{(n+1) \times (n+1)}\) a matrix whose columns correspond to the \(\boldsymbol{p}_i\) vectors of each vertex composing simplex \(\boldsymbol{\Delta}_q\), and \[S_q = \dfrac{1}{2}\left(P_q^T \bar{X}_q^{-1} + \bar{X}_q^{-T} P_q\right) \label{eq:Sq}\tag{3}\] the local quadratic parametrization of the Piecewise Quadratic function.

With this description of the finite-element linear interpolation as a Piecewise Quadratic function, we can write that \[\bar{V}(\boldsymbol{x}) = V_q(\boldsymbol{x}), \; \textrm{if } \boldsymbol{x} \in \boldsymbol{\Delta}_q. \label{eq:switch95V}\tag{4}\] Since this parametrization arises from a linear interpolation of the linear functions defined for every grid point, it is naturally continuous over the state space.

In this work, we consider that a control law, satisfying \(\boldsymbol{u}(\boldsymbol{x}) \in \mathcal{U}, \forall \boldsymbol{x} \in \mathcal{X}\), is given by a (possibly discontinuous) Piecewise Linear function of the form \[\boldsymbol{u}(\boldsymbol{x}) = \begin{bmatrix} K_q & \boldsymbol{k}_q \end{bmatrix} \begin{bmatrix} \boldsymbol{x} \\ 1 \end{bmatrix} = \bar{K}_q \bar{\boldsymbol{x}}, \; \textrm{if } \boldsymbol{x} \in \boldsymbol{\Delta}_q,\label{eq:K}\tag{5}\] with \(\bar{K}_q \in \mathbb{R}^{m \times (n+1)}\) attributing a gain for each simplex in the triangulation and \(\boldsymbol{k}_{q} = 0\) if \(\boldsymbol{\Delta}_q \in \mathcal{X}_g\), to ensure that \(\boldsymbol{u}(\boldsymbol{0}) = \boldsymbol{0}\). With this control law, we can write that the closed loop dynamics are given by \[\begin{align} \dot{\bar{\boldsymbol{x}}} = \sum_{k=1}^r \alpha_i \bar{A}_{kq} \bar{\boldsymbol{x}}, \end{align}\] with \[\begin{align} \bar{A}_{kq} = \begin{bmatrix} A_k + B_k K_q & B_k \boldsymbol{k}_q \\ \boldsymbol{0} & \boldsymbol{0} \end{bmatrix}. \label{eq:x95Akq} \end{align}\tag{6}\]

Since we do not impose continuity upon the control law, the closed loop vector field may be discontinuous on the boundaries of the simplices. Ergo, the generalized solution of the closed loop dynamics will lie in the convex hull of both possible switching dynamics, and additional constraints need to be added to the optimization problem to ensure they will also satisfy our conditions. Furthermore, in order to impose that \(\bar{V}(\boldsymbol{x}) \geq 1\) on the boundaries of \(\mathcal{X}\), we also require a way to describe these boundary regions.

Note that, in both cases (either boundaries with neighboring simplices or with the boundary of the analysis region), these regions will be described by the \(n\) dimensional faces of the simplices, \(\boldsymbol{\delta}_\ell\). Considering that there are \(r_f\) faces, representing all of the \(n\) dimensional faces on the simplices of our triangulation, similarly to the simplex case, we define a mapping \(\eta_f(\ell,j): \mathcal{I}_{r_f} \times \mathcal{I}_{n} \rightarrow \mathcal{I}_{r_p}\), which given an index \(\ell\) representing a face \(\boldsymbol{\delta}_\ell\), and an index \(j\) representing one vertex of said face, returns an index corresponding to a point on the grid \(\Gamma\). Given this mapping, each face can be defined as the convex hull of a set of \(n\) vertex points, such that \[\begin{align} \boldsymbol{\delta}_\ell &= \operatorname{co}\left(\boldsymbol{x}_{\eta_f(\ell,1)}, \dots, \boldsymbol{x}_{\eta_f(\ell,n)}\right), \end{align}\] which implies that there exists \(\bar{\boldsymbol{\beta}} \in \mathbb{R}^{n}\) unique convex weights that can be used to represent any point in the face from its vertices as \(\boldsymbol{x} = \sum_{j = 1}^{n} \bar{\beta}_j \boldsymbol{x}_{\eta_f(q,j)}\), with \(\bar{\boldsymbol{\beta}}^T = \begin{bmatrix} \bar{\beta}_1 & \dots & \bar{\beta}_{n}\end{bmatrix}\), \(\bar{\beta}_{j} \geq 0 \; \forall j, \; \sum_{j=1}^{n} \bar{\beta}_j = 1\). We also consider that there exists a function \(\nu(\ell) : \mathcal{I}_{r_f} \rightarrow \mathcal{I}_{r_s} \times \mathcal{I}_{r_s}\), which, given an index representing a face, returns the indices of the two neighboring simplices that share said face, or the same simplex twice (if the face is on the border of \(\mathcal{X}\), it is only contained by a single simplex, which gets returned twice).

2.2 Policy Iteration↩︎

Like many other LMI based approaches, if we simply tried to find a set of solutions from Theorem 1 that directly tried to solve for both the Lyapunov function \(\bar{V}(\boldsymbol{x})\) and the control law, we would end up with a set of Bilinear Matrix Inequalities (BMIs). In order to avoid that, as in many other works in the literature, we tackle the problem using Policy Iteration, which amounts to the Policy Evaluation and Policy Improvement steps being repeated iteratively.

Policy Evaluation is a step in which, given a control law/policy, \(\boldsymbol{u}(\boldsymbol{x})\) in 5 , for the closed loop system, we try to find the smallest upper-bound Lyapunov function, \(\bar{V}(\boldsymbol{x})\) in 4 . In contrast to this, Policy Improvement is a step in which, given a Lyapunov function, \(\bar{V}(\boldsymbol{x})\) in 4 , finds the control law that tries to maximize its decay. Whereas most optimization-based ADP approaches have a closed-form solution for the Policy Improvement step [14], [15], [21], [22], due to the fact that the control set \(\mathcal{U}\) is bounded, and that our control law is parametrized by 5 , we instead also propose a set of LMI conditions to solve this step.

The conditions for both steps are presented in Theorems 2 and 3 in the following.

3 Examples↩︎

In this section we present 3 examples to illustrate the proposed approach. Throughout them, the grid of points was chosen so that the region around the origin would be reasonably small and the proposed approach would steer the system towards the origin fast. In addition to this, \(\gamma\) was chosen as a small value to ensure the positivity of the Lyapunov function inside of \(\mathcal{X}_g\) while not affecting much the results outside of it.

3.1 Linear System - Double Tank System↩︎

Consider the double tank system illustrated in Figure 2. A simplified version of its dynamics can be described by \[\begin{align} \dot{h}_1 &= \frac{-1}{R_{12}C_1}h_1 + \frac{1}{R_{12}C_1}h_2 + \frac{\overline{q}_i}{C_1}\delta, \\ \dot{h}_2 &= \frac{1}{R_{12}C_2}h_1 - \left(\frac{R_{12} + R_o}{R_{12}R_o C_2}\right)h_2, \end{align}\] with \(h_1\) and \(h_2\) representing the level of both tanks, \(C_1\) and \(C_2\) are the area of the section of each tank, \(\bar{q}_i\) represents the maximum input flow, \(\delta \in [0,1]\) represents the opening of the input valve, and \(R_{12}\) and \(R_o\) are hydraulic resistances. Without loss of generality, we will consider that the problem is driving the level of the second tank to a desired equilibrium value \(\bar{h}_2\), with a corresponding equilibrium level \(\bar{h}_1\) for the first tank and equilibrium value for the input valve \(\bar{\delta}\). With a small change of coordinates, such that this desired equilibrium point is the origin of the state space, the dynamics of the system are given by \[\begin{align} \dot{\boldsymbol{x}} = \begin{bmatrix} \dfrac{-1}{R_{12}C_1} & \dfrac{1}{R_{12}C_1} \\ \dfrac{1}{R_{12}C_2} & - \left(\dfrac{R_{12} + R_o}{R_{12}R_o C_2}\right) \end{bmatrix} \boldsymbol{x} + \begin{bmatrix} \dfrac{\overline{q}_i}{C_1} \\ 0 \end{bmatrix} u, \end{align}\] with \(\boldsymbol{x} = \begin{bmatrix} h_1 - \bar{h}_1 & h_2 - \bar{h}_2 \end{bmatrix}^T\) and \(u = \delta - \bar{\delta}\). The parameters employed for this system were \(C_1 = 2 m^2\), \(C_2 = 4 m^2\), \(\bar{q}_i = 0.4 m^3/s\), \(R_{12} = 1 s/m^2\), \(R_o = 10 s/m^2\), \(\bar{h}_1 = 2.2 m\), \(\bar{h}_2 = 2 m\), and \(\bar{\delta} = 0.5\). In these new coordinates, we have that \(x_1 \in [-2.2, 2.2], x_2 \in [-2, 2]\) and \(u \in [-0.5, 0.5]\).

In order to solve this problem, we employ the Guaranteed Time Control approach proposed in Theorems 2 and 3, by making use of a regular grid of \(15 \times 15\) points in the state space, leading to 225 points and 392 simplices. Since the system is open-loop stable, the initial controller (used to start the policy iteration) is chosen as not using any control (\(K_q = 0 \; \forall q)\), and 7 iterations were performed (alternating between Theorems 2 and 3), considering \(\gamma = 10^{-3}\).

The results were initially compared against a linear control law that imposes the maximum decay rate it can, without ever saturating inside of the analysis region and a saturated linear control law [10] that imposes the maximum decay rate while guaranteeing the largest estimated domain of attraction. This comparison is shown in Figure 3, for an initial condition that corresponds to both tanks completely empty, as well as in Figure 4, which depicts the comparison of the estimates of the domains of attraction from the proposed conditions and the saturated linear control law [10].

As shown in Figure 3, the proposed control law is capable of saturating the control signal for a significantly larger amount of time than the saturated linear control law, leading to a significantly faster time to reach the origin than both the linear and saturated linear control laws. In addition to this, as shown in Figure 4, the estimate of the domain of attraction of the origin in closed loop is considerably larger in our approach (consisting of almost all of the state space) compared against the estimate from the conditions in [10] (which is the largest area a quadratic Lyapunov function could estimate in this example).

In a different set of comparisons, regarding the behavior of the closed loop system starting from the empty tanks initial condition, the proposed approach was also compared against two different finite-time control approaches, one based on an implicit Lyapunov function [7], and another based on the idea of Prescribed Performance Control [8]. The results are shown in Figure 5, with the performance functions used for the prescribed performance control tuned slightly slower in order for us to be able to set the responses apart.

As shown in the figure, the control law based on the proposed approach is capable of driving the system towards the desired goal with a similar closed loop behavior as the finite-time control approaches, however with a considerably smoother control action.

3.2 Unstable Uncertain Polytopic System↩︎

In this example, we consider a numerical problem, in which we have an open loop unstable uncertain polytopic system composed by 2 vertices, described by 1 with matrices \[\begin{align} &A_1 = \begin{bmatrix} -0.2868 & -3.5353 \\ -3.5353 & -8.7132 \end{bmatrix}, && A_2 = \begin{bmatrix} -0.2868 & 3.5353 \\ 3.5353 & -8.7132 \end{bmatrix}, \end{align}\] \(B_1 = B_2 = \begin{bmatrix} 1 & -1 \end{bmatrix}^T\). We consider that \(x_1 \in [-5, 5]\), \(x_2 \in [-5, 5]\) and \(u \in [-10, 10]\).

Similarly to the previous example, the Guaranteed Time Control approach proposed in Theorems 2 and 3 is employed, making use of a regular grid of \(15 \times 15\) points in the state space, leading to 225 points and 392 simplices.

Since the system is open-loop unstable, we make use of an initial controller that imposes a minimum decay rate of 1 to a quadratic Lyapunov function, while minimizing the controller gain (but ignoring the control bounds). This controller can be found by solving the optimization problem \[\begin{align} &\min_{X,Y,s} s, \\ s.t.\\ &X\geq I, \\ &A_iX+XA_i^T + BY+Y^TB^T + 2X \leq 0, \\ &\begin{bmatrix} -s & Y \\ \ast &-X \end{bmatrix} \leq 0, \end{align}\] which leads to the initial gain \(K = YX^{-1} = \begin{bmatrix} -3.2668 & -1.0985\end{bmatrix}\). After this initialization, 10 iterations were performed (alternating between Theorems 2 and 3), considering \(\gamma = 10^{-3}\). The results were compared against an adaptation of the Semi-Lagragian approximation of time-optimal control problems in [16] to a zero-sum dynamic game setting in which one agent, controlling \(\boldsymbol{\alpha}\), tries to maximize the time taken to reach the origin, and the other, controlling \(u\), tries to minimize the time.

This comparison is shown in Figures 6 and 7, which show respectively the time evolution of the closed loop system, and the \(\bar{V}(\boldsymbol{x})\) Lyapunov function as well as the behavior from 4 different initial conditions. As can be seen in the figures, the results obtained are comparable to the ones found from the optimal control approach (though slightly slower) with the added benefit that the control law found is already guaranteed to work in closed loop (since unlike the optimal control approach no approximation of the problem are done to solve it).

3.3 Nonlinear System - Chua’s circuit↩︎

Consider the chaotic oscillator, known as Chua’s circuit [27], with a controlled current source added illustrated in Figure 8. If we consider that \(V_{C_1}\) and \(V_{C_2}\) are expressed in \(V\), whereas \(i_L\) and \(u\) are expressed in \(mA\), it follows that the system can be described by \[\begin{align} \dot{V}_{C_1} &= \frac{1}{C_1} \left(\frac{V_{C_2} - V_{C_1}}{R} - G(V_{C_1}) V_{C_1} + \frac{u}{1000}\right), \\ \dot{V}_{C_1} &= \frac{1}{C_2} \left(\frac{V_{C_1} - V_{C_2}}{R} + \frac{i_L}{1000}\right), \\ \frac{d}{dt} i_L &= \frac{-1000 V_{C_2}}{L}. \end{align}\] with \(G(V_{C_1})\) the conductance of the Chua’s diode, given by \[\begin{align} G(V_{C_1}) = \left\{\begin{array}{ll} G_a, & \textrm{if } |V_{C_1}| < E \\ G_b + (G_a-G_b)\dfrac{E}{|V_{C_1}|}, & \textrm{otherwise}. \end{array}\right. \end{align}\] The parameters employed for this system were \(C_1 = 30.14 \mu F\), \(C_2 = 185.66 \mu F\), \(L = 52.28 H\), \(R = 1673 \Omega\), \(G_a = -0.801 mS\), \(G_b = -0.365 mS\), \(E = 1.74 V\). We consider that \(V_{C_1} \in [-6, 6]V\), \(V_{C_2} \in [-3, 3]V\), \(i_{L} \in [-3, 3]mA\) and \(u \in [-200, 200]mA\).

Similarly to the previous examples, the Guaranteed Time Control approach proposed in Theorems 2 and 3 is employed, making use of a regular grid of \(11 \times 5 \times 5\) points in the state space, leading to 275 points and 960 simplices. Since we are dealing with a nonlinear system, for each simplex we employ a local Takagi-Sugeno model (as discussed in Remark 2) by considering the local maximum and minimum values of \(G(V_{C_1})\). Since the system is open-loop unstable, we make use of an initial linear stabilizing controller (that does not saturate in the region \(\mathcal{X}_g\) around the origin) given by \(K = \begin{bmatrix}-38.8017 & -31.1041 & 25.3298\end{bmatrix}\) and 18 iterations were performed (alternating between Theorems 2 and 3), considering \(\gamma = 10^{-3}\).

The level sets for the Lyapunov function found for this system are illustrated in Figure 9, whereas the 0.99 sublevel-set (an estimate of the closed loop’s domain of attraction) is illustrated in Figure 10 along with several trajectories with initial conditions on the border of this set. In addition to this, Figure 11 presents the closed loop behavior of the system, starting from initial condition \(\begin{bmatrix}-2.5 & 2.5 & 2.5 \end{bmatrix}^T\), and Figure 12 compares the estimate of the domain of attraction found using the Guaranteed Time Control approach against the estimate found from using the conditions in [28] (modified to maximize the volume of an ellipsoidal volume inside of the domain of attraction, using \(\max \operatorname{logdet}\) instead of maximizing the volume of a ball contained inside of the domain of attraction), with \(|\dot{h}_i| \leq 40\), and the one found using the conditions in [29].

Note that, similar to the other examples, the Guaranteed Time Control approach proposed in this paper was capable of finding a control law whose estimated domain of attraction covers most of the desired analysis region, and completely encompasses the estimate found from the competing approaches. Figure 10 shows that the estimate found is indeed a valid approximation of the domain of attraction as the system states converge towards the origin from the initial conditions picked on the border of the estimate found. Finally, Figure 11 show that the control bounds are respected, and that the extra conditions, imposed over the faces of the simplices, are indeed needed in this example since it exhibits a switching behavior in its control signal (likely due to discontinuities on the control law) but is still ensured to stabilize the system.

4 Conclusion and Future Directions↩︎

This paper proposed a Guaranteed Time Control approach for dynamical systems. In order to do so, Theorem 1 proposed a set of Lyapunov conditions, which, if satisfied, ensures asymptotic stability, finds an upper-bound for the time taken to reach the desired set \(\mathcal{X}_g\) that contains the origin, and ensures that the strict 1-sublevel set of the Lyapunov function is an estimate of the origin’s domain of attraction. Theorems 2 and 3 make use of these conditions, as well as the novel piecewise representation presented in Section 2.1, to propose a Policy Iteration approach for the problem (circumventing the bilinear terms that appear if we try to solve for both Lyapunov function and control law simultaneously).

Three examples were presented to illustrate the proposed approach and compare it against the literature: a linear system, an unstable uncertain polytopic system, and a nonlinear system. The first example showed that the approach is capable of achieving considerably better results than simple linear conditions (both in terms of time taken to reach the origin, and on the guaranteed domain of attraction) and that it is capable of finding comparable results with finite time control laws (both based on ILF and PPC projects) while finding a considerably smoother control law. The second example illustrates that, given a dense enough discretization of the state space, the proposed approach is capable of finding results comparable to the solution of a time-optimal control law. The third example shows that the proposed approach is capable of handling highly nonlinear systems (given that it is a chaotic system in open-loop) and that the guaranteed domain of attraction is superior to the current state-of-the-art conditions.

Even though the proposed approach was capable of handling the examples presented in this paper, like many other approaches based on discretizing the state space in some manner, it suffers from the curse of dimensionality with the number of simplices increasing exponentially as the dimension of the problem grows. In that regard, the adaptation of Domain Decomposition Methods [30] and possible adaptive grids [31] are very interesting research directions as they would allow for a decomposition/parallelization of the solution, as well as a partition with a reduced number of simplices.

Other possible future directions include the use of control laws that incorporate further information about the system, either via Parallel Distributed Compensation [26] for the nonlinear systems, or adaptive polytopic controllers [32] for the uncertain systems.

Acknowledgments↩︎

This paper was partially funded by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) through grants 308597/2023-0, 408419/2023-7, Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG) through grants APQ-04246-25, Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), and Université Polytechnique Hauts-de-France (UPHF).

References↩︎