Tropical Kraus maps for optimal control of switched systems

TL;DR

This paper introduces tropical Kraus maps, a novel mathematical tool for approximating value functions in switched optimal control problems, offering improved scalability over existing methods.

Contribution

The paper develops tropical analogues of Kraus maps and demonstrates their application in approximating value functions and solving switched control problems more efficiently.

Findings

Major scalability improvements over previous schemes

Effective approximation of value functions in switched control

Successful application to joint spectral radius and Hamilton-Jacobi PDEs

Abstract

Kraus maps (completely positive trace preserving maps) arise classically in quantum information, as they describe the evolution of noncommutative probability measures. We introduce tropical analogues of Kraus maps, obtained by replacing the addition of positive semidefinite matrices by a multivalued supremum with respect to the L\"owner order. We show that non-linear eigenvectors of tropical Kraus maps determine piecewise quadratic approximations of the value functions of switched optimal control problems. This leads to a new approximation method, which we illustrate by two applications: 1) approximating the joint spectral radius, 2) computing approximate solutions of Hamilton-Jacobi PDE arising from a class of switched linear quadratic problems studied previously by McEneaney. We report numerical experiments, indicating a major improvement in terms of scalability by comparison with earlier numerical schemes, owing to the "LMI-free" nature of our method.

Tables

Table 1. Table 1. Comparison of the methods with respect to the size of the dimension of the matrices.

Dimension $n$	CPU time (tropical)	CPU time (LMI)	Upper bound on $\rho \left(𝒜\right)$ (tropical)	Upper bound on $\rho \left(𝒜\right)$ (LMI)
$5$	$0.9$ s	$3.1$ s	$2.767$	$2.7627$
$10$	$1.5$ s	$4.2$ s	$3.797$	$3.7426$
$20$	$3.5$ s	$31$ s	$5.4093$	$5.3891$
$30$	$7.9$ s	$3$min	$6.2038$	$6.1942$
$40$	$13.7$ s	$18$min	$7.3402$	$7.3363$
$45$	$18.1$ s	$-$	$7.687$	$-$
$50$	$25.2$ s	$-$	$8.1591$	$-$
$100$	$1$min	$-$	$11.487$	$-$
$500$	$8$min	$-$	$25.44$	$-$

Table 2. Table 2. Comparison of the methods w.r.t. the size of the graph D d subscript 𝐷 𝑑 D_{d} .

Order $d$	$2$	$4$	$6$	$8$	$10$
Size of $𝒲$	$8$	$32$	$128$	$512$	$2048$
CPU time (tropical)	$0.03$s	$0.07$s	$0.4$s	$2.0$s	$9.0$s
CPU time (LMI)	$1.9$s	$4.0$s	$24$s	$1$min	$10$min
Upper bound on $\rho \left(𝒜\right)$ (tropical)	$1.842$	$1.821$	$1.804$	$1.800$	$1.801$
Upper bound on $\rho \left(𝒜\right)$ (LMI)	$1.8216$	$1.7974$	$1.7957$	$1.7922$	$1.7905$

Table 3. Table 3. Numerical benchmarks of the tropical Kraus method applied to McEneaney’s switched linear quadratic problem

Example	$1$	$2$	$2$	$3$	$4$
Dimension	$2$	$6$	$6$	$20$	$20$
Size of $\Sigma$	$3$	$6$	$6$	$2$	$4$
$\tau$	$0.05$ s	$0.2$	$0.1$	$0.1$	$0.1$
Size of $𝒲$	$81$	$216$	$1296$	$128$	$256$
Initial error	$0.78$	$1.12$	$1.12$	$4.2$	$4.79$
Final error	$0.047$	$0.071$	$0.090$	$0.0006$	$0.17$
Iterations	$194$	$115$	$200$	$55$	$288$
CPU time	$8$ s	$41$ s	$5$ min	$5$ s	$2.5$ min

Equations

$$x_{k+1}=A_{\sigma (k)}{x}_k,\sigma (k)\in \{1,\dots ,m\}$$

$$\rho (\mathcal{A})$$

$$\lambda v(x)=\underset{1⩽i⩽m}{\max }v(A_{i}x),\forall x\in {\mathbb{R}}^n,$$

$$\underset{1⩽i⩽m}{\max }v(A_{i}x)⩽\mu v(x),\forall x\in {\mathbb{R}}^n,$$

$$A_i^{\top }Q{A}_i\preccurlyeq {\mu }^{2}Q\forall i,$$

$$V(x)=\underset{u\in \mathcal{U}}{\sup }\underset{\mu \in \mathcal{D}}{\sup }\underset{t>0}{\sup }{\int }_0^t\frac{1}{2}\xi (s{)}^{\top }{D}^{\mu (s)}\xi (s)-\frac{{\gamma }^2}{2}|u(s){|}^{2}ds.$$

$$\dot{\xi }(s)=A^{\sigma }\xi (s)+B^{\sigma }u(s),\xi (0)=x,$$

$$H(x,\nabla V)=0,x\in {\mathbb{R}}^n.$$

$$H^{\sigma }(x,p)=(A^{\sigma }x{)}^{\top }p+\frac{1}{2}{x}^{\top }{D}^{\sigma }x+\frac{1}{2}{p}^{\top }{Q}^{\sigma }p,$$

$$\displaystyle\begin{split}S_{t}[V^{0}](x)=\sup_{u\in\mathcal{U}}\sup_{\mu\in\mathcal{D}}&\int_{0}^{t}\frac{1}{2}\xi(s)^{\top}D^{\mu(s)}\xi(s)\\ &-\frac{\gamma^{2}}{2}|u(s)|^{2}\,ds+V^{0}\big{(}\xi(t)\big{)}\,.\end{split}$$

$$\displaystyle\begin{split}S_{t}^{\sigma}[V^{0}](x)=\sup_{u\in\mathcal{U}}&\int_{0}^{t}\frac{1}{2}\xi(s)^{\top}D^{\sigma}\xi(s)\\ &-\frac{\gamma^{2}}{2}|u(s)|^{2}\,ds+V^{0}\big{(}\xi(t)\big{)}\,.\end{split}$$

$$\dot{P}=(A^{\sigma }{)}^{\top }P+P{A}^{\sigma }+P{Q}^{\sigma }P+D^{\sigma },P(0)=P_0,$$

$$P\preccurlyeq Q⟺Q-P\in {\mathcal{S}}_n^{+}⟺Q-P\succcurlyeq 0.$$

$$Z\mapsto ⟨C,Z⟩$$

$$\displaystyle X_{C}=\frac{P+Q}{2}+\frac{1}{2}C^{-1/2}\,\big{|}C^{1/2}(P-Q)C^{1/2}\big{|}\,C^{-1/2}\,.$$

$$\displaystyle T_{j}(X)\coloneqq\bigvee\Big{\{}A_{\sigma}^{\top}X_{i}A_{\sigma}\colon(i,\sigma)\in\mathcal{W}\times\Sigma\,,i\cdot\sigma=j\Big{\}}\,.$$

$$\displaystyle(M_{\tau})_{j}(X)\coloneqq\bigvee\Big{\{}\operatorname{ricc}_{\tau,\sigma}X_{i}\colon(i,a)\in\mathcal{W}\times\Sigma\,,i\cdot\sigma=j\Big{\}}\,.$$

$$\displaystyle\begin{split}\bigvee\Big{\{}f_{\alpha}(X)&\colon\alpha=(i,\sigma,j)\;,i\cdot\sigma=j\Big{\}}\\ &\text{with}\;f_{\alpha}(X)=A_{\alpha}^{\top}XA_{\alpha}\,.\end{split}$$

$$v(z):=\underset{j\in \mathcal{W}}{\sup }(z^{\top }{X}_{j}z{)}^{1/2}$$

$$\underset{\sigma }{\max }{S}_t^{\sigma }[V](z)⩽V(z)\text{for all}z.$$

$$\displaystyle T^{\text{sel}}_{j}(X)\coloneqq\bigsqcup\Big{\{}A_{\sigma}^{\top}X_{i}A_{\sigma}\colon i\cdot\sigma=j\Big{\}}\,.$$

$$\displaystyle\widehat{T^{\text{sel}}}(X)\coloneqq\frac{1}{2}\Big{[}\frac{1}{\langle I_{n}^{p},T^{\text{sel}}(X)\rangle}T^{\text{sel}}(X)+X\Big{]}\,.$$

$$\displaystyle X^{(k+1)}=\frac{1}{2}\Big{[}\frac{1}{\langle I_{n}^{p},T^{\text{sel}}(X^{(k)})\rangle}T^{\text{sel}}(X^{(k)})+X^{(k)}\Big{]}\,.$$

$$\displaystyle X^{(k+1)}=\Big{[}\frac{T^{\text{sel}}(X^{(k)})}{\langle I_{n}^{p},T^{\text{sel}}(X^{(k)})\rangle}\Big{]}\#X^{(k)}\,,$$

$$X^{(k+1)}=M_{\tau }^{\text{sel}}(X^{(k)}).$$

$$\displaystyle T^{\text{sel}}_{j}(X)\coloneqq\bigsqcup\Big{\{}A_{\sigma}^{\top}X_{i}A_{\sigma}+\varepsilon I_{n}\colon i\cdot\sigma=j\Big{\}}\,.$$

$$(i,\sigma ,j)\in E⟺A_{\sigma }^{\top }{X}_i{A}_{\sigma }\preccurlyeq {\rho }^2{X}_j.$$

$$V\approx \underset{{\sigma }_1,\dots ,{\sigma }_N\in \Sigma }{\sup }{S}_{\tau }^{{\sigma }_1}\cdots S_{\tau }^{{\sigma }_N}[V^0],$$

$$\rho (\mathcal{A})⩽\rho (\mathcal{A},\mathcal{G})⩽\hat{\rho }(\mathcal{A},\mathcal{G}).$$

$$\delta :=(\hat{\rho }(\mathcal{A},{\mathcal{G}}^{\prime })-\rho (\mathcal{A},\mathcal{G})/\rho (\mathcal{A},\mathcal{G}).$$

Full text

The authors were partially supported by the ANR projects CAFEIN, ANR-12-INSE-0007 and MALTHY, ANR-13-INSE-0003, by ICODE and by the academic research chair “Complex Systems Engineering” of École polytechnique - THALES - FX - DGA - DASSAULT AVIATION - DCNS Research - ENSTA ParisTech - Télécom ParisTech - Fondation ParisTech - FDO ENSTA and by the PGMO programme of EDF and FMJH.

Stéphane Gaubert and Nikolas Stott are with INRIA and CMAP, École polytechnique, CNRS, France, firstname.lastname@inria.fr

1. INTRODUCTION

1.1. Curse of dimensionality attenuation methods

Dynamic programming is one of the main methods to solve optimal control problems. It characterizes the value function as the solution of a functional equation or of a Hamilton-Jacobi partial differential equation. It provides a feedback law that is guaranteed to be globally optimal. However, it is subject to the “curse of dimensionality”. Indeed, the main numerical methods, including monotone finite difference or semi-Lagrangean schemes [10, 8, 11, 9], and the anti-diffusive schemes [6], are grid-based. It follows that the time needed to obtain an approximate solution with a given accuracy is exponential in the dimension of the state space.

Recently, two numerical methods have been shown, theoretically or practically, to attenuate the curse of dimensionality, for specific classes of optimal control problems with switches.

McEneaney considered hybrid optimal control problems in which a discrete control allows one to switch between different linear quadratic models. The method he developed [24] approximates the value function by a supremum of elementary functions like quadratic forms, hence it belongs to the family of “max-plus basis methods” [13, 2]. The method of [24] has a remarkable feature: it attenuates the curse of dimensionality, as shown by the complexity estimates of Kluberg and McEneaney [27] and of Qu [28]. McEaneney’s method [24] has been studied and extended in a series of works [31, 14, 19, 26].

A different problem consists in computing the joint spectral radius of a finite set of matrices [17]. This can be formulated as an ergodic optimal control problem for a switched system. In this context, the ergodic value function is known as the Barabanov norm. Specific numerical methods have been developed, which approximate the Barabanov ball by a polytope [16], or are of semi-Lagrangean type [20]. Ahmadi et al. [1] developed a new method, based on a path complete automaton. It approximates the Barabanov norm by a supremum of quadratic norms. Whereas the worst case complexity estimates in [1] are still subject to a curse of dimensionality, in practice, the efficiency of the method is determined by the complexity of the optimal switching law rather than by the dimension itself. It allows one to solve instances of dimension inaccessible by grid-based method.

The method of McEneaney [24], like the one of Ahmadi et al. [1], provide max-plus basis expansions of approximate value functions. In both methods, semidefinite programming (solution of LMI, linear matrix inequalities) is the bottleneck. Indeed, LMI arise in the “pruning step”, i.e., the elimination of redundant quadratic forms, an essential ingredient of McEaneney’s method. It was observed in [14] that 99% of the CPU time was spent in the solution of semidefinite programs. The method of Ahmadi et al. involves a truncation method, considering switching sequences of a fixed length, and it requires the solution of a semidefinite program whose size is exponential in this length.

1.2. Contribution

In this paper, we introduce a new approximation method for optimal control of switched systems. This method still relies on the approximation of the value function by a supremum of quadratic forms. However, it avoids the recourse to LMI. We exploit the geometric properties of the space of positive semidefinite matrices equipped with the Löwner order. A key ingredient is the introduction of the tropical analogues of the Kraus maps arising in quantum information and control [30]. The latter are quantum Markov operators, describing the evolution of density matrices (the quantum analogues of probability measures). They act on the space of positive semidefinite matrices, and can be written as $T(X)=\sum_{i}A_{i}XA_{i}^{\dagger}$, where $(\cdot)^{\dagger}$ denotes the adjoint of a matrix. Tropical Kraus maps are defined by replacing the sum in the definition of $T(X)$ by a multivalued operator, providing the set of minimal upper bounds with respect to the Löwner order. Tropical Kraus maps may be thought of as “1-player” non-commutative dynamic programming operators (classical Kraus maps correspond to the 0-player case). We show that every non-linear eigenvector of the tropical Kraus map yields an approximation of the value function. Moreover, the non-linear eigenvalue yields an upper bound for the joint spectral radius. We show that a non-linear eigenvector does exist. We compute non-linear eigenvectors through an iterative scheme, in which at each stage, a specific minimal upper bound of a collection of matrices is selected. The latter operation is implemented in an algebraic way, leading to a fast algorithm (not relying on semidefinite programming). We report numerical experiments, showing a major speedup, allowing us to treat instances of dimension inaccessible by earlier dynamic programming methods.

The paper is organized as follows. In Section 2, we recall the definitions of the joint spectral radius and switched linear quadratic control problems, which will serve as benchmarks. In Section 3, we recall some basic properties of the Löwner order on the space of positive semidefinite matrices. In Section 4, we introduce tropical Kraus maps, establish the existence of non-linear eigenvectors, and present the iterative scheme. In Section 5, we present numerical experiments.

2. CLASSES OF SWITCHED SYSTEMS

In this section, we describe the two optimal control problems to which we will apply our method.

2.1. Joint spectral radius

We consider here the stability under arbitrary switching of discrete-time linear switched systems as studied in e.g. [23, 33, 7].

Let $\mathcal{A}=\{A_{1},\dots,A_{m}\}$ be a set of real $n\times n$ matrices. A discrete-time switched linear system is described by:

[TABLE]

where $x_{k}\in\mathbb{R}^{n}$ denotes the trajectory of the system, and $\sigma$ is the switching mechanism, which selects one of the matrices in $\mathcal{A}$ at each instant.

We are interested in the approximation of the joint spectral radius [17] associated to this set. The latter is defined by

[TABLE]

A fundamental result of Barabanov [4] shows that if $\mathcal{A}$ is irreducible, meaning that there is no nontrivial subspace of $\mathbb{R}^{n}$ that is left invariant by every matrix in $\mathcal{A}$, then there is a norm $v$ on $\mathbb{R}^{n}$ such that

[TABLE]

for some $\lambda>0$. The scalar $\lambda$ is unique and it coincides with the joint spectral radius $\rho(\mathcal{A})$. This shows that, when $\mathcal{A}$ is irreducible, all the trajectories of the switched linear system converge to zero if and only if $\rho<1$.

The norm $v$ is known as a Barabanov norm. A norm which satisfies the inequality $\max_{i}v(A_{i}x)\leqslant\rho(\mathcal{A})v(x)$ for all $x\in\mathbb{R}^{n}$ is called an extremal norm.

Extremal norms and Barabanov norms are generally non unique and cannot be computed exactly, except in special examples. Hence, we shall be content with an approximate extremal norm $v$, i.e., a solution of

[TABLE]

where $\mu>0$. Then, it is readily seen that $\mu\geqslant\rho(\mathcal{A})$, so that an approximate extremal norm yields a guaranteed approximation of the joint spectral radius.

For instance, an approximate extremal norm can be obtained via an ellipsoidal norm, as described in [3], meaning that there is some positive definite matrix $Q$ and a real $\mu$ such that

[TABLE]

where $\preccurlyeq$ denotes the Löwner order (see Section 3) and $(\cdot)^{\top}$ denotes the transpose. The approximation is usually coarse, since $1\leqslant\mu/\rho(\mathcal{A})\leqslant\sqrt{n}$ is a tight estimate. This approach can be refined by lifting the set of matrices $\mathcal{A}$ into another set $\mathcal{A}^{\prime}$ which has the same joint spectral radius, but where the quadratic form $Q$ gives a better approximation of an extremal norm [1]. We adopt such an approach in this paper, see Remark 2 in Section 4.1.

2.2. Linear Quadratic Optimal Control Problems with Switches

We also consider the following problem of optimal switching between linear quadratic models, studied by McEneaney [24], namely approximating the value function $V$ of an optimal control problem having both a control $u$ taking values in $\mathbb{R}^{p}$ and a discrete control (switches between different modes) $\mu$ taking values in $\Sigma\coloneqq\{1,\dots,m\}$:

[TABLE]

Here, $\mathcal{D}$ denotes the set of measurable functions from $[0,+\infty)$ to $\Sigma$ (i.e. switching functions), $\mathcal{U}\coloneqq L^{2}([0,+\infty),\mathbb{R}^{p})$ is the space of $\mathbb{R}^{p}$-valued control functions, and the state $\xi$ is subject to

[TABLE]

where $\sigma=\mu(s)$ denotes the mode that is selected at time $s$.

It is known [24] that, under some assumptions on the parameters, the value function $V$ takes finite values and is the unique viscosity solution of the stationary Hamilton-Jacobi-Bellman PDE:

[TABLE]

The Hamiltonian $H(x,p)$ in the latter equation is the point-wise maximum of simpler Hamiltonians $H^{\sigma}(x,p)$ given for $\sigma\in\Sigma$ by

[TABLE]

and $Q^{\sigma}=\gamma^{-2}B^{\sigma}(B^{\sigma})^{\top}$.

We associate with this problem the Lax-Oleinik semi-group $\{S_{t}\}_{t\geqslant 0}$ defined by

[TABLE]

McEneaney showed in [24] that $V(x)$ coincides with $\lim_{t\to+\infty}S_{t}[V^{0}](x)$ and that the latter limit is uniform on compact sets if $V^{0}$ satisfies a quadratic growth condition (one requires that $\epsilon|x|^{2}\leqslant V^{0}(x)\leqslant\lambda|x|^{2}$ for some positive constants $\epsilon,\lambda$ that are determined from the parameters).

We also associate to every value $\sigma\in\Sigma$ the semi-group $\{S_{t}^{\sigma}\}_{t\geqslant 0}$ corresponding to the unswitched control problem obtained by setting $\mu(s)\equiv\sigma$, i.e.,

[TABLE]

Computing $S_{\tau}^{\sigma}[V^{0}]$ when $V^{0}(x)=x^{\top}P_{0}x$, reduces to solving the following indefinite Riccati differential equation,

[TABLE]

with $P(s)\in\mathcal{S}_{n}$. Indeed, we have $S_{\tau}^{\sigma}[V^{0}](x)=x^{\top}P(\tau)x$. We denote by $\operatorname{ricc}_{\tau,\sigma}$ the flow of this equation, so that $\operatorname{ricc}_{\tau,\sigma}[P_{0}]\coloneqq P(\tau)$.

3. Minimal upper bounds in the Loewner order

We begin by recalling some standard notation and definitions. We denote by $\wp(X)$ the powerset of a set $X$. We denote by $\mathcal{S}_{n}$ the space of symmetric matrices, which is equipped with the Frobenius scalar product defined by $\langle P,Q\rangle=\operatorname{trace}(PQ)$. The product space $(\mathcal{S}_{n})^{p}$ is equipped with the scalar product $\langle(P_{k})_{k},(Q_{k})_{k}\rangle=\sum\langle P_{k},Q_{k}\rangle$. The $n\times n$ identity matrix is denoted by $I_{n}$ and we use the shorthand $I_{n}^{p}$ to mean the $p$-tuple $(I_{n},\dots,I_{n})$.

A symmetric matrix $P$ is positive semidefinite when the quadratic form $x^{\top}Px=\sum_{i,j}P_{i,j}x_{i}x_{j}$ takes non-negative values for all vectors $x\in\mathbb{R}^{n}$, or equivalently when all the eigenvalues of $P$ are non-negative. Then, we write $P\succcurlyeq 0$. The set of positive semidefinite matrices is denoted by $\mathcal{S}_{n}^{+}$. When $x^{\top}Px$ is positive for all nonzero vectors $x\in\mathbb{R}^{n}$, we say that the matrix $P$ is positive definite. The set of positive semidefinite matrices constitutes a convex cone in $\mathcal{S}_{n}$, meaning that $\lambda P+\mu Q\in\mathcal{S}_{n}^{+}$ for all $P,Q\in\mathcal{S}_{n}^{+}$ and non-negative $\lambda,\mu$. It is also closed and pointed ($\mathcal{S}_{n}^{+}\cap-\mathcal{S}_{n}^{+}=\{0\}$), thus it defines an order relation on $\mathcal{S}_{n}$ by

[TABLE]

The partial order $\preccurlyeq$ is called the Löwner order.

A classical result by Kadison [18] shows that the set $\mathcal{S}_{n}$ equipped with this order constitutes an antilattice, meaning that two matrices $P,Q\in\mathcal{S}_{n}$ have a supremum (least upper bound) if and only if they are comparable, meaning that $P\succcurlyeq Q$ or $Q\succcurlyeq P$. If the matrices $P,Q$ are not comparable, then they possess a continuum of minimal upper bounds, i.e. upper bounds $S$ such that $P,Q\preccurlyeq X\preccurlyeq S$ implies $X=S$, see [32] for more information.

Given a finite set of symmetric matrices $\mathcal{Q}$, we denote by $\mathcal{Q}^{\uparrow}\coloneqq\{X\in\mathcal{S}_{n}\colon X\succcurlyeq Q_{i}\,,Q_{i}\in\mathcal{Q}\}$ the set of upper bounds of the matrices in $\mathcal{Q}$. This set is convex as an intersection of convex sets. We also denote by $\bigvee{\mathcal{Q}}\in\wp(\mathcal{S}_{n})$ the subset of $\mathcal{Q}^{\uparrow}$ consisting of all minimal upper bounds of $\mathcal{Q}$. We use the symbol $\bigvee$ to denote a “supremum” operation which is multivalued owing to the antilattice character of $\mathcal{S}_{n}$. The set $\bigvee{\mathcal{Q}}$ coincides with the set of positively exposed points $\mathcal{Q}^{\uparrow}$:

Theorem 1**.**

The matrix $X$ is a minimal upper bound of a finite set of matrices $\mathcal{Q}$ if and only if there is a positive definite matrix $C$ such that $X$ minimizes the map

[TABLE]

*over the set $\mathcal{Q}^{\uparrow}$. The minimizer, denoted by $X_{C}$, is unique, and, when the set $\mathcal{Q}$ consists of two matrices $P,Q$, *

[TABLE]

Here, $C^{1/2}$ denotes the unique positive definite solution to the equation $X^{2}=C$ and $|X|=(XX^{\top})^{1/2}$.

*Sketch of proof: * This is deduced from the optimality condition of the associated semidefinite program and from a generalization of the characterization of minimal upper bounds in [32, Theorem 3.1]. When $C=I_{n}$, it can be checked that the matrix $X_{C}=P+|Q-P|$ satisfies the optimality conditions. The formula in the general case is obtained by a change of variable $X\mapsto C^{1/2}XC^{1/2}$ and a symetrization in $P,Q$.

 ∎

We say that the minimal upper bound $X_{C}$ is selected by the matrix $C$. Note that the expression of $X_{C}$ is similar to that of the maximum of two scalars: $\max(a,b)=(a+b+|a-b|)/2$. When more than two matrices are involved, minimal upper bounds can be computed by solving a semidefinite program. By choosing $C=I_{n}$ in Theorem 1, we obtain a special minimal upper bound of $\mathcal{Q}$, denoted by $\sqcup_{\operatorname{tr}}\mathcal{Q}$.

Finally, we point out a remarkable selection of a minimal upper bound of two positive semidefinite matrices. If all matrices in $Q_{i}\in\mathcal{Q}$ are positive definite, a minimal upper bound is given by $(X^{*})^{-1}$, where $X^{*}$ denotes the unique matrix that maximizes $\log\det X$ over all positive definite matrices $X$ such that $X\preccurlyeq Q_{i}^{-1}$ for all $i$. We denote this selection by $\sqcup_{\det}\mathcal{Q}$. This minimal upper bound corresponds to the minimum volume ellipsoid enclosing the ellipsoids $\{x\in\mathbb{R}^{n}\colon x^{\top}Q_{i}^{-1}x\leqslant 1\}$ and has thus received much attention across several fields, see [3] and references therein. In particular, when the set $\mathcal{Q}$ consists of two matrices $P,Q$, it has been shown in [3] that $\sqcup_{\det}\mathcal{Q}$ can be obtained by selecting $C=P^{-1}$ (or $C=Q^{-1}$) in (1).

4. Tropical Kraus Maps

4.1. Definitions

In the sequel, we assume that we are given an index set $\Sigma=\{1,\dots,m\}$, $m$ matrices $\mathcal{A}=\{A_{\sigma}\}_{\sigma\in\Sigma}$. We also assume that we are given a finite set $\mathcal{W}$ with $p$ elements and a map sending $\mathcal{W}\times\Sigma$ to $\mathcal{W}$, denoted by $(i,\sigma)\mapsto i\cdot\sigma$, for $i\in\mathcal{W}$ and $\sigma\in\Sigma$. We say that $(i,\sigma,j)\in\mathcal{W}\times\Sigma\times\mathcal{W}$ is an admissible transition when $i\cdot\sigma=j$.

We now introduce the tropical Kraus map $T$, from $(\mathcal{S}_{n}^{+})^{p}$ to $(\wp(\mathcal{S}_{n}^{+}))^{p}$, whose $j$-th coordinate maps $X=(X_{1},\dots,X_{p})\in(\mathcal{S}_{n}^{+})^{p}$ to the subset of $\mathcal{S}_{n}^{+}$:

[TABLE]

Remark 1**.**

In the setting of quantum information [21], a Kraus map is given by $X\mapsto\sum_{i}A_{i}XA_{i}^{\dagger}$ and acts on the set of density matrices (positive semidefinite matrices of trace $1$). We say that the map $T$ is a tropical Kraus map since it ressembles the latter, except the sum has been replaced with the “supremum” operation $\bigvee$, and the matrices have been transposed. The transposition is not surprising: classical Kraus maps provide a forward propagation of density matrices, whereas we are interested in Lyapunov functions, whose propagation follows a backward scheme. I.e., the present tropical Kraus maps are analogues to the adjoints of classical Kraus maps. Finally, as we work with real quadratic forms, instead of hermitian forms, the hermitian conjugate ${\dagger}$ is replaced by transposition ⊤.

We also consider a variant of the tropical Kraus map, adapted to the optimal control problem in Section 2.2, defining the map $M_{\tau}$ from $(\mathcal{S}_{n}^{+})^{p}$ to $(\wp(\mathcal{S}_{n}^{+}))^{p}$ by:

[TABLE]

Remark 2**.**

A tropical Kraus map from $(\mathcal{S}_{n}^{+})^{p}$ to $(\wp(\mathcal{S}_{n}^{+}))^{p}$ can be represented by a new tropical Kraus map from $\mathcal{S}_{np}^{+}$ to $\wp(\mathcal{S}_{np}^{+})$, preserving the set of block-diagonal matrices. Indeed, for $X=\operatorname{diag}(X_{1},\dots,X_{m})$, these maps can be written as

[TABLE]

where $A_{\alpha}=E_{ij}\otimes A_{\sigma}$ denotes the lifted set of matrices. Here, $E_{ij}$ denotes the matrix with $1$ in the $(i,j)$-th entry and [math] everywhere else and $\otimes$ is the Kronecker product. Indeed, whenever the selection matrix $C$ in Theorem 1 is block-diagonal $C=(C_{1},\dots,C_{m})$, the associated minimal upper bound is also block diagonal and the value of the $j$-th block is exactly the minimal upper bound in $T_{j}(X)$ that is selected by $C_{j}$. The case where the matrix $C$ is not block-diagonal does not appear in our analysis. The same remark applies to the variant $M_{\tau}$ of the tropical Kraus map.

4.2. Non linear eigenvalue and fixed points problems associated to Tropical Kraus Maps

The tropical Kraus map $T$ is positively homogeneous, meaning that $T(\alpha X)=\alpha T(X)$ for all $X\in\mathcal{S}_{n}^{+}$ and $\alpha\geqslant 0$. This suggests to consider a multivalued eigenproblem. A (non-linear) eigenvector of $T$, associated to the eigenvalue $\lambda$ is a nonzero matrix $X\in\mathcal{S}_{n}^{+}$ such that $\lambda X_{j}\in T_{j}(X)$ holds for all $j\in\mathcal{W}$. We write $\lambda X\in T(X)$ for brevity. This notation is licit since we can identify $T(X)$ which is an element of $(\wp(\mathcal{S}_{n}))^{p}$ to an element of \wp\big{(}(\mathcal{S}_{n})^{p}\big{)}.

The following result shows that a non-linear eigenvalue of the tropical Kraus map provides an upper bound for the joint spectral radius.

Theorem 2**.**

If the multivalued eigenvector problem $\lambda X\in T(X)$ has a solution such that the matrix $\sum_{j\in\mathcal{W}}X_{j}$ is positive definite, then, the map

[TABLE]

is a norm, and $v(A_{\sigma}z)\leqslant\sqrt{\lambda}v(z)$ holds for all $z\in\mathbb{R}^{n}$ and $\sigma\in\Sigma$. In particular, the joint spectral radius of $\mathcal{A}$ does not exceed $\sqrt{\lambda}$.

Proof.

The proof of this result is similar to the proof Theorem 2.4 by Ahmadi et al. [1]. Indeed, a non-linear eigenvector $X$ of the tropical Kraus map $T$ provides a feasible point of the semidefinite program considered in [1, Theorem 2.4]. ∎

We have the following analogous result for the switched linear quadratic control problem.

Theorem 3**.**

If the multivalued fixed point problem $X\in M_{\tau}(X)$ has a solution, then the map $V\coloneqq z\mapsto\sup_{i}z^{\top}X_{i}z$ determines a sub-invariant function of the Lax-Oleinik semi-group $S_{t}$, meaning that:

[TABLE]

4.3. Existence of non-linear eigenvectors of tropical Kraus maps and computation by a Krasnoselskii-Mann iteration

For a completely positive map, $X\mapsto\sum_{i}A_{i}XA_{i}^{\top}$, the existence of a positive semidefinite eigenvector follows from the Perron-Frobenius theorem [22]. Moreover, such an eigenvector is necessarily positive definite as soon as the map is irreducible in the Perron-Frobenius sense, meaning that the map does not leave invariant a non-trivial face of the closed cone $\mathcal{S}_{n}^{+}$. As shown in [12], the latter condition holds if and only if the set of matrices $\{A_{i}\}$ is irreducible in the algebraic sense, meaning that there is no non-trivial subspace invariant by each matrix in this set.

In order to show that tropical Kraus maps have eigenvectors, we specialize the multivalued map $T$ defined in Section 4.1 by fixing a selection of minimal upper bound $\sqcup$. We obtain the map $T^{\text{sel}}$ defined on $(\mathcal{S}_{n}^{+})^{p}$ by

[TABLE]

We will prove that the map $T^{\text{sel}}$ has a non-linear eigenvector if the selection $\sqcup$ is the ”minimum volume” selection $\sqcup_{\det}$. We introduce the “non-commutative simplex” $\Delta_{p}\coloneqq\{X\in(\mathcal{S}_{n}^{+})^{p}\colon\langle I_{n}^{p},X\rangle=1\}$ and the map $\widehat{T^{\text{sel}}}$ sending $\Delta_{p}$ to itself:

[TABLE]

Observe that, independently of the selection $\sqcup$, a fixed point $X\in\Delta_{p}$ of the map $\widehat{T^{\text{sel}}}$ yields an eigenvector for the map $T^{\text{sel}}$ associated with the eigenvalue $\langle I_{n}^{p},T^{\text{sel}}(X)\rangle$. We can now state the theorem.

Theorem 4**.**

If the set of matrices $\{E_{ij}\otimes A_{\sigma}\colon i\cdot\sigma=j\}$ is irreducible, then the map $\widehat{T^{\text{sel}}}$ has a positive definite fixed point.

*Sketch of proof: * The proof relies on the inequality $p^{-1}\sum_{1\leqslant k\leqslant p}Q_{k}\preccurlyeq\sqcup_{\det}Q_{k}\preccurlyeq\sum_{1\leqslant k\leqslant p}Q_{k}$ which can be deduced from [3, Theorem 4.1] when the set $\mathcal{A}$ is reduced to two matrices. The irreducibility of $\mathcal{A}$ implies that there is an integer $q$ such that the map $\widehat{T^{\text{sel}}}$ iterated $q$ times sends $\Delta_{p}$ into its interior. We then show the existence of a convex compact set $K$ included in the interior of $\Delta$ which is invariant by $\widehat{T^{\text{sel}}}$. The operator $(Q_{1},\dots,Q_{p})\mapsto\sqcup_{\det}\{Q_{1},\dots,Q_{p}\}$ is continuous on the interior of $\Delta$, hence, $\widehat{T^{\text{sel}}}$ is continuous on $K$. We conclude by applying Brouwer’s fixed-point theorem to $\widehat{T^{\text{sel}}}$.

 ∎

We obtain as an immediate corollary:

Corollary 1**.**

If the set of matrices $\{E_{ij}\otimes A_{\sigma}\colon i\cdot\sigma=j\}$ is irreducible, then, the tropical Kraus map $T$ has a positive definite eigenvector.

Remark 3**.**

*Several basic methods allow one to prove non-linear extensions of the Perron-Frobenius theorem. These involve contraction properties with respect to Hilbert’s projective metric, Brouwer fixed point theorem, or monotonicity properties, see [22]. A direct application of all these methods fails in the case of the tropical Kraus maps $T^{\text{sel}}$, which are not contracting, not monotone, and which do not have continuous extensions to the closure of the cone on which they act. This is why the proof of Corollary 1 relies on the detour through $\widehat{T^{\text{sel}}}$ in Theorem 4. *

In order to compute a fixed point of the map $\widehat{T^{\text{sel}}}$, we compute successive iterates starting from a positive definite matrix $X^{(0)}$. This yields a Krasnoselskii-Mann-type scheme:

[TABLE]

This is a power-type iteration, involving a renormalization and a “damping term” (addition of $X^{(k)}$) to avoid oscillations. This should be compared with the classical Krasnoselskii-Mann iteration, which applies to non-expansive mappings $T$, and takes the form $X^{(k+1)}=(T(X^{(k)})+X^{(k)})/2$, see [29].

Remark 4**.**

There is a multiplicative variant of the iteration, defined by

[TABLE]

where P\#Q\coloneqq P^{1/2}\big{(}P^{-1/2}QP^{-1/2}\big{)}^{1/2}P^{1/2} denotes the Riemannian barycenter of the positive definite matrices $P,Q$, see [5, Chapter 2] for more information. We can show that this multiplicative version does converge in the “commutative case”, i.e., when $n=1$. Then, the map $X\mapsto{\langle I_{n}^{p},T^{\text{sel}}(X)\rangle}^{-1}{T^{\text{sel}}(X)}$ is nonexpansive in the Hilbert metric [22], and then, the general result of [29] can be applied. The additive version can also be shown to be converging when $n=1$, by a reduction to the same result, but the proof is more involved.

We use a different iteration scheme to compute fixed points of the variant $M_{\tau}$. First, we specialize again the multivalued map $M_{\tau}$ with a minimal upper bound selection $\sqcup$ to obtain the map $M_{\tau}^{\text{sel}}$. Then, we compute iteratively

[TABLE]

We can show that this iteration converges on some ”good” instances of the problem, when we choose $\sqcup\coloneqq\sqcup_{\det}$. Indeed, the (indefinite) Riccati flow is a contraction in the Thompson metric, with contraction rate $\alpha>0$ determined by the parameters of the flow [15, Corollary 4.7]. Moreover, as stated by Allamigeon et al. [3], the selection $\sqcup_{\det}$ has a Lipschitz constant in the Thompson metric which is not larger than $1+(4/\pi)\log n$ (this bound is conservative). Combining these results, the iteration in Equation 2 is guaranteed to converge locally when $\exp(\alpha\tau)>1+(4/\pi)\log n$. The contraction rate $\alpha$ depends on an interval $\{X\in\mathcal{S}_{n}\colon\lambda_{1}I_{n}\preccurlyeq X\preccurlyeq\lambda_{2}I_{n}\}$ and this interval is not preserved by $\sqcup_{\det}$, hence we do not have global convergence. Determining whether there is a minimal upper bound selection that preserves this interval and that has a finite Lipschitz constant in Thompson’s metric, to be used instead of $\sqcup_{\det}$ remains an open problem.

4.4. Implementation issues

We describe in this section the resolution to several issues that arise in the implementation of the iterative scheme.

First, in the iterative scheme to approximate the joint spectral radius, we introduce a small positive perturbation $\varepsilon$ in the computation:

[TABLE]

In practice, we use values for $\varepsilon$ in the range $10^{-4}-10^{-2}$. This additional parameter allows us to obtain, in a finite number of iterations, a solution $(X,\rho)$ that satisfies $\rho^{2}X_{j}\succcurlyeq A_{\sigma}^{\top}X_{i}A_{\sigma}$ for all admissible $(i,\sigma,j)$. Moreover, this parameter absorbs numerical imprecisions that may appear during the computation and ensures that the matrices $X_{j}$ are positive definite, so the assumptions of Theorem 2 and Theorem 3 are satisfied.

We choose the selection $\sqcup\coloneqq\sqcup_{\text{tr}}$. Then, when the set $\Sigma$ contains only two elements, $T^{\text{sel}}_{j}(X)$ can be computed analytically thanks to Theorem 1. When $\Sigma$ has more than two elements, instead of computing the true minimal upper bound $\sqcup_{\operatorname{tr}}\mathcal{Q}$, we compute an approximation by sequential evaluation: $Q_{1}\sqcup_{\operatorname{tr}}(Q_{2}\sqcup_{\operatorname{tr}}(\dots\sqcup_{\operatorname{tr}}Q_{p}))$.

Finally, as pointed out in [14], the propagation of the Riccati operator $\operatorname{ricc}_{\tau,\sigma}$ on a single quadratic form $P_{0}$ is computed analytically by $\operatorname{ricc}_{\tau,\sigma}=Y(\tau)X(\tau)^{-1}$, with $\mathcal{M}^{\sigma}=\begin{pmatrix}-A^{\sigma}&-Q^{\sigma}\\ D^{\sigma}&A^{\sigma}\end{pmatrix}$ and \big{(}X(\tau)\,;\,Y(\tau)\big{)}^{\top}=\exp\big{(}\mathcal{M}^{\sigma}\tau)\big{(}I_{n}\,;\,P_{0}\big{)}^{\top}.

5. Experimental results

5.1. Path-complete graph Lyapunov functions

In [1], Ahmadi and al. developed a method to compute an overapproximation of the joint spectral radius of a finite set of matrices, to which we shall compare our method.

Given a set of states $\mathcal{W}$ and an alphabet $\Sigma$, an edge of a labeled graph is a triple $(i,\sigma,j)\in\mathcal{W}\times\Sigma\times\mathcal{W}$. The set of edges is denoted $E$. Such a graph is called path-complete if for every state $i$ and letter $\sigma$, there is some state $j$ such that $(i,\sigma,j)$ is an edge.

Let $\mathcal{A}=\{A_{\sigma}\}_{\sigma\in\Sigma}$ denote a finite set of $n\times n$ matrices and $\rho$ a non-negative real number. In [1], the authors examine graphs, denoted $\mathcal{G}(X,\rho)$, whose states are positive definite matrices $\{X_{i}\}_{i}$ and whose edges are determined by

[TABLE]

The main theorem in [1] shows that the construction of a path-complete graph $\mathcal{G}(X,\rho)$ gives an upper bound of the joint spectral radius:

Theorem 5** (Theorem 2.4 [1]).**

If the graph $\mathcal{G}(X,\rho)$ is path-complete for some set of positive definite matrices $\{X_{i}\}_{i}$, then $\rho(\mathcal{A})\leqslant\rho$. Moreover, the map $V\colon z\mapsto\max_{i}z^{\top}X_{i}z$ is a Lyapunov-type function: it satisfies $V(A_{\sigma}x)\leqslant\rho^{2}V(x)$ for all $\sigma\in\Sigma$ and $x\in\mathbb{R}^{n}$.

In practice, for a fixed value of $\rho$ and a given path-complete graph $\mathcal{G}$, checking the existence of a path-complete graph $\mathcal{G}(X,\rho)$ whose edges coincide with $\mathcal{G}$ amounts to checking the feasibility of an LMI. A bisection scheme is then implemented to refine $\rho$. For brevity, we shall refer to this method as the LMI method.

A class of graphs which provides good theoretical and experimental approximations is the class of De Bruijn graphs. The set of states of the De Bruijn graph of order $d$ is the set $\Sigma^{d}$ of words built on $\Sigma$ which have length $d$. There is an edge $(i,\sigma,j)$ between states $i$ and $j$ if and only if $i=\sigma_{1}\dots\sigma_{d}$ and $j=\sigma_{2}\dots\sigma_{d}\sigma$. This graph, denoted by $D_{d}$, is path-complete by construction.

5.2. McEneaney’s curse of dimensionality attenuation scheme

We assume that $V^{0}$ is a quadratic function $V^{0}(x)=x^{\top}P^{0}x$. The method of [24] that solves the linear quadratic optimal control problem described in Section 2.2 approximates the value function $V$ by a finite supremum of quadratic forms

[TABLE]

where $\tau$ is a (small) time discretization step and $N$ is a maximal number of switches. The latter supremum represents the value of a modified optimal control problem, in horizon $\tau N$, in which switches occur only at times multiple of $\tau$. We have $S_{\tau}^{\sigma_{1}}\cdots S_{\tau}^{\sigma_{N}}[V^{0}](x)=x^{\top}Qx$, where $Q=\operatorname{ricc}_{\tau,\sigma_{1}}\circ\dots\circ\operatorname{ricc}_{\tau,\sigma_{N}}(P_{0})$, can be computed by integrating successive Riccati equations, which allows us to evaluate the expression in Equation 3.

The propagation of a quadratic form by the Lax-Oleinik semi-group has only a cubic cost in terms of the dimension $n$, contrary to classical grid-based methods whose cost is exponential in the dimension. In this sense, the curse of dimensionality has been reduced. However, the memory footprint of this method is exponential in the number of switches, since $m^{N}$ quadratic forms are computed after $N$ iterations. Several pruning schemes have been proposed in [14] to limit this growth. This is a costly operation, indeed, $99\%$ of the computation time is spent solving LMIs inside the pruning procedure [14].

5.3. Application to the joint spectral radius

Given that the approximation of the joint spectral radius $\rho(\mathcal{A})$ depends on the graph $\mathcal{G}$ that underlies the analysis, we denote by $\widehat{\rho}(\mathcal{A},\mathcal{G})$ the approximation obtained as (the square root of) an eigenvalue of a tropical Kraus map and by $\rho(\mathcal{A},\mathcal{G})$ the one obtained by solving LMIs.

The map $\cdot$ sending $\mathcal{W}\times\Sigma$ to $\mathcal{W}$ defined in Section 4.1 can be interpreted as a path-complete graph. For this reason, our method, when applied to the joint spectral radius, is a relaxation of the path-complete Lyapunov function framework, and thus we always have

[TABLE]

However, we shall see that the tropical method is much more tractable, so we may use a bigger graph and sometimes get a better approximation than by solving LMIs.

We compare the performance of our algorithm with the path-complete graph Lyapunov method, in terms of computation time and accuracy of the approximation of the joint spectral radius, measured by

[TABLE]

All the experiments were implemented in Matlab, running on one core of a $2.2$ GHz Intel Core i7 with $8$ GB RAM. The semidefinite programs were solved using YALMIP (R$20160930$), calling SeDuMi $1.3$.

5.3.1. Accuracy of the approximation

We generate $600$ pairs $\mathcal{A}=\{A_{1},A_{2}\}$ of random $6\times 6$ matrices. For each of these pairs, we compare the approximation of the joint spectral radius obtained by the LMI method on the graph $D_{3}$ (involving $8$ positive semidefinite matrices) and by the tropical Kraus method on the graph $D_{6}$ (involving $64$ positive semidefinite matrices). On these examples, we report that the tropical method obtains a similar approximation of the joint spectral radius, within a margin of $2.5\%$, and outperforms the LMI-method on $25\%$ of these examples. Moreover, whereas the LMI-method requires between $3$s and $5$s to obtain this approximation, the tropical method consistently returns an approximation in $1$s.

5.3.2. Scalability - dimension

We generate random pairs of $n\times n$ matrices, for $n$ ranging from $5$ to $500$. We use again the De Bruijn graph $D_{3}$ in the LMI-method and the graph $D_{6}$ in the tropical Kraus method. We show in Table 1 the mean computation time required to obtain an overapproximation and the mean relative accuracy of the tropical method with respect to the LMI-method, when it applies. First, one can observe the major speedup provided by the tropical method, from $4$ times faster when $n=5$ to $80$ times faster for $n=40$.

Also note that the tropical method is using $8$ times more quadratic forms in its analysis and remains much faster than the LMI-method. Thus, given a fixed time budget, the tropical method enjoys more flexibility regarding the size of the graph that is used in the analysis.

Moreover, observe that the LMI-method cannot provide estimates on the joint spectral radius for values of $n$ greater than $45$, whereas the tropical method easily reaches values of $n$ greater than $100$.

Finally, the accuracy of the tropical approximation remains within a $3\%$ margin of the one obtained by the LMI-method.

5.3.3. Scalability - graphs

We now analyze the influence of the order of the De Bruijn graph $D_{d}$ used in the analysis on the computation of the upper bound on the joint spectral radius obtained by both methods. We use the matrices $A_{1}=\begin{pmatrix}[r]-1&1&-1\\ -1&-1&1\\ 0&1&1\end{pmatrix}$ and $A_{2}=\begin{pmatrix}[r]-1&1&-1\\ -1&-1&0\\ 1&1&1\end{pmatrix}$, introduced in [16]. Their joint spectral radius is $\rho(\mathcal{A})=1.78893$.

We show in Table 2 the upper bound on the joint spectral radius and the computation time with respect to the length order $d$ of the De Bruijn graphs.

5.4. A faster curse of dimensionality attenuation scheme

We now apply the iteration scheme described in Section 4.3 to the approximation of the value function $V$. In all examples, we measure the quality of the approximation of the value function as in [25, 14] with the H-infinity back-substitution error $\max_{x^{\top}x\leqslant 1}|H(x,\nabla V(x))|$ on the subspace spanned by the canonical vectors $e_{1}$ and $e_{2}$.

The first example is Example $1$ in [24] and we use the instance of [14] in the second example. Examples $3$ and $4$ are randomly generated examples that satisfy the technical assumptions in [24].

Table 3 depicts the results of the computations. In particular, we give the backsubstitution error at the beginning of the computation, when the value function is approximated by a single quadratic form ($Q(x)=0.1|x|^{2}$ in all cases) and the final backsubstitution error when the scheme has converged.

6. CONCLUDING REMARKS

We introduced a new method to approximate the value function of optimal control problems for switched systems. This method applies to situations in which the evolution semi-group of the unswitched problem preserves the space of quadratic forms. This includes the computation of the joint spectral radius and a class of linear quadratic control problems with switches considered by McEneaney. Our scheme belongs to the family of max-plus methods as it approximates the value function by a supremum of quadratic forms. It avoids the recourse to semidefinite programming (which was the bottleneck of earlier max-plus methods) by a reduction to a non-linear eigenproblem, exploiting the geometry of the Löwner order. This leads to a major speedup, allowing us to obtain approximate solutions of instances in dimension up to $100$ in the case of the joint spectral radius, and 20 for McEneaney’s problem, hardly accessible by other methods.

Let us now point out the limitations of the present approach, together with possible ways to overcome them.

A key ingredient in our method is the replacement of LMI formulations by a selection of minimal upper bounds in the Löwner order. This induces a “relaxation gap”, which is difficult to estimate as it depends on the specific selection which is used. Another difficulty is that these selections may be expansive, resulting in a potential instability or lack of convergence of the iterative scheme. In the case of the joint spectral radius, experiments indicate that the scheme does converge (although a proof of convergence in the general case is missing). In the case of switched linear quadratic control problems, we do have a proof of convergence in the discrete time case. This proof requires the contraction rates of the Riccati flows arising in our problem to be sufficiently small to absorb the expansiveness of the selection of a joint. In practice, the scheme converges in more general circumstances. However, in the continuous time case, the precision of the scheme becomes limited as it can be unstable for small values of the time discretization step. In other words, the scheme currently allows one to compute quickly a coarse approximation of the solution of a Hamilton-Jacobi PDE. The most promising improvement of the scheme may be to adapt dynamically the selection of a minimal upper bound, which will reduce the relaxation gap, and might also improve the convergence. This is left for further work.

Acknowledgements

We thank the reviewers for their helpful comments.