Curse-of-Complexity Attenuation Using Semidefinite Programming in Solution of HJB PDEs William M. McEneaney and Ameet Deshpande Abstract— Recently, a curse-of-dimensionality-free method was developed for solution of Hamilton-Jacobi-Bellman partial differential equations (HJB PDEs) for nonlinear control problems, using semiconvex duality and max-plus analysis. The curse-of-dimensionality-free method may be applied to HJB PDEs where the Hamiltonian is given as (or well-approximated by) a pointwise maximum of quadratic forms. Such HJB PDEs also arise in certain switched linear systems. The method constructs the correct solution of an HJB PDE from a maxplus linear combination of quadratics. The method completely avoids the curse-of-dimensionality, and is subject to cubic computational growth as a function of space dimension. However, it is subject to a curse-of-complexity. In particualr, the number of quadratics in the approximation grows exponentially with the number of iterations. Efficacy of such a method depends on the pruning of quadratics to keep the complexity growth at a reasonable level. Here we apply a pruning algorithm based on semidefinite programming. Computational speeds are exceptional.

I. I NTRODUCTION Dynamic programming is an extremely robust tool for solving nonlinear optimal control problems. In the case of deterministic optimal control, or in the case of deterministic games where one player’s feedback is prespecified, the dynamic programming equation reduces to a Hamilton-JacobiBellman (HJB) PDE. The difficulty is that one must solve the HJB PDE. Various approaches have been taken to solving the HJB PDE. The most common methods are grid-based methods (c.f., [6], [7]). Although highly refined at this point, these methods still suffer from the curse-of-dimensionality, as the number of grid points and computations grow exponentially with the space dimension. However, in recent years, entirely new classes of numerical methods for HJB PDEs have emerged (c.f., [1], [2], [8], [13], [16]). These methods exploit the max-plus linearity of the associated semigroup. In the previous work of the first author [12], [13], [15], a new method based on above semigroup linearity was proposed for certain nonlinear HJB PDEs, and this method was free from the curse-of-dimensionality. In fact, the computational growth in state-space dimension is on the order of n3 . However, there is exponential computational growth in a certain measure of complexity of the Hamiltonian. Under this measure, the minimal complexity Hamiltonian is the linear/quadratic Hamiltonian – corresponding to solution by Research partially supported by NSF grant DMS-0307229, AFOSR grant FA9550-06-1-0238 and INRIA. William McEneaney is with the Dept. of Mech. and Aero. Eng., University of California San Diego, San Diego, CA 92093-0411, USA.

[email protected] Ameet Deshpande is similarly affiliated. [email protected]

a Riccati equation. If the Hamiltonian is given as a pointwise maximum of M linear/quadratic Hamiltonians, then one could say the complexity of the Hamiltonian is M . Such PDEs also arise in switched linear systems. The algorithm transforms the original problem to its maxplus dual form, where the dual of the value function is expressed as a max-plus sum, i.e., a pointwise maximum, of certain quadratic basis functions. An infinite time-horizon problem is considered, and as such, the value function is approximated by iterating a finite-horizon semigroup until a large enough propagation horizon is reached. With the curseof-dimensionality-free method, this finite-horizon semigroup is approximated by a semigroup whose semiconvex dual is represented as a finite number of quadratic forms. The dual of the approximate value at each iteration is stored as a set of quadratic functions. Acting on this dual with the above dual semigroup leads to a new approximation, where the number of quadratics grows by a fixed factor (the number of quadratics representing the approximate finitehorizon dual semigroup) at each iteration. (The algorithm is discussed briefly below, and one may find further details in [12], [13], [15].) This is the curse-of-complexity. To contain such computational growth, a pruning method based on semidefinite programming (SDP) is considered here. II. P ROBLEM

STATEMENT AND ASSUMPTIONS

The HJB PDEs we consider arise in infinite-horizon nonlinear optimal control problems, and their Hamiltonians are given as (or well-approximated by) pointwise maxima of linear-quadratic functions. Note that pointwise maxima of quadratic forms can approximate, arbitrarily closely, any semiconvex function. More specifically, we consider e 0 = −H(x, ∇V ) = −

max

m∈{1,2,...,M}

V (0) = 0

{H m (x, ∇V )} (1) (2)

(i.e., with boundary condition V = 0 at the origin) where each of the constituent Hamiltonians has the form H m (x, p) =

1 T m 1 T m 2x D x + 2p Σ p +(l2m )T p + αm ,

+ (Am x)T p + (l1m )T x (3)

where Dm , Σm are n × n symmetric, matrices, l1m , l2m ∈ IRn and αm ∈ IR. e is associated with an optimal control problem for H . switched linear systems. Let M = {1, 2, . . . M }. The corre-

sponding value function, maximizing payoff Je is e T ; w, µ) Ve (x) = sup sup sup J(x, w∈W µ∈D∞ T <∞

= sup sup sup

w∈W µ∈D∞ T <∞

where

Z

T

Lµt (ξt ) −

0

γ2 2 2 |wt |

dt (4)

Lµt (x) = 21 xT Dµt x + (l1µt )T x + αµt , D∞ = {µ : [0, ∞) → M : measurable }, . k W = Lloc 2 ([0, ∞); IR ),

and ξ satisfies dynamics ξ˙ = Aµt ξ + l2µt + σ µt wt ,

Note that the behavior specified in (A.w) is proved in the purely quadratic case (c.f., [13]) under reasonable assumptions on the constituent-Hamiltonian matrices, but in this more general case, we need to assume it instead. Lastly, we make the following assumption. Assume there exist T , c3 ∈ (0, ∞) such that for all x ∈ IRn , all ε ∈ (0, 1], and all µε , wε which e µε , wε ) ≥ are ε–optimal for Ve (i.e., such that J(x, Ve (x) − ε ), one has (A.ξ) Z T Z T ε 2 Lµt (ξtε ) dt ≥ c3 |ξtε | dt ∀T ≥ T 0

ξ0 = x

(5)

where σ m and γ are such that Σm = γ12 σ m (σ m )T . Here µt is a switching control which appears in addition to the control w. e T ; w, µ). Let Ve (x, T ) = supw∈W supµ∈D∞ J(x, To motivate the assumptions for this rather general probe as being constructed so as to lem class, we consider H resemble some given nonlinear control problem which has a e as being chosen to (finite) solution. That is, we think of H ee resemble some given H. We suppose that problem e e 0 = −H(x, ∇V ), V (0) = 0 (6)

e to approximate has finite value, and that we are choosing H e n e H. Let QK = {φ : IR → IR | φ is semiconvex, and 0 ≤ φ(x) ≤ (K/2)|x|2 ∀x ∈ IRn }. we may take QK as the domain of the semigroup. We make the following block of assumptions. Assume there exists unique viscosity solution, e Ve , to (6) in QK for some K ∈ (0, ∞), where QK = {φ : IRn → IR | φ is semiconvex, and 0 ≤ Φ(x) ≤ (K/2)|x|2 ∀x ∈ IRn }. e Assume that H(x, p) = maxm∈M H m (x, p) ≤ e e H(x, p) for all x, p ∈ IRn . Assume H 1 (x, p) has coefficients satisfying the following: l11 = l21 = 0; α1 = 0; there exists cA ∈ (0, ∞) such that x′ A1 x ≤ −cA |x|2 ∀x ∈ IRn ; D1 (A.m) is positive definite, symmetric; and γ 2 /|σ 1 |2 > cD /c2A , where cD is such that x′ D1 x ≤ cD |x|2 ∀x ∈ IRn . Assume that system (5) is controllable in the sense that given x, y ∈ IRn and T > 0, there exist processes w ∈ W and µ measurable with range in M, such that ξT = y when ξ0 = x and one applies controls w, µ

Note that the last of these assumptions, the controllability assumption, is satisfied if there exists at least one m ∈ M such that σ m (σ m )′ (which is n × n) has n positive eigenvalues. Assume there exist c1 , c2 < ∞ such that for any e problem, one has ε–optimal pair, µε , wε for the H (A.w) kwε k2L2 [0,T ] ≤ c1 + c2 |x|2 for all ε ∈ (0, 1], all T < ∞ and all x ∈ IRn .

0

ξ˙tε

µε l2 t

ε Aµt ξtε

ε

where = + + σ µt wtε , ξ0ε = x. Note that these last two assumptions might be difficult to verify. Easily verifiable assumptions appear in [13], [15], but these generate a significantly smaller class of systems than those for which these methods apply. Now, define the operator Z T 2 e ST [φ] = sup sup Lµt (ξt ) − γ2 |wt |2 dt + φ(ξT ) w∈W µ∈DT

0

where DT = {µ : [0, T ) → M : measurable }. Under the above assumptions, the viscosity solution, Ve of (1),(2) exists, e satisfies 0 ≤ Ve ≤ Ve and is given by Ve = limT →∞ SeT [V0 ] for any V0 ∈ QK such that 0 ≤ V0 ≤ Ve , [14], [12]. In max-plus algebra, addition and multiplication are defined as a ⊕ b = max(a, b) and a ⊗ b = a + b, respectively. It is well known that SeT forms a max-plus linear semigroup. That is, if k1 , k2 ∈ IR, and T1 , T2 ≥ 0, then SeT1 [SeT2 [φ]] = SeT1 +T2 [φ] SeT [k1 ⊗ φ1 ⊕ k2 ⊗ φ2 ] = k1 ⊗ SeT [φ1 ] ⊕ k2 ⊗ SeT [φ2 ]

III. C URSE - OF - DIMENSIONALITY- FREE ALGORITHM

The key steps in the curse-of-dimensionality-free algorithm developed in [12] are given below. Since we are interested in understanding how the curse-of-complexity arises in this algorithm, we shall sidestep the theoretical foundations which are well covered in [12], [13], and focus on the algorithmic flow. A. Approximate propagation Define the consituent-Hamiltonian semigroup operators as Z τ 2 Sτm [φ] = sup Lm (ξt ) − γ2 |wt |2 dt + φ(ξτ ). w∈W

0

Importantly, propagation of a quadratic φ by an Sτm operator can be reduced to solution of a differential Riccati equation. Define the time-indexed operators M S¯τ [φ](x) = max Sτm [φ](x) = Sτm [φ](x). m∈M

m∈M

Fix any T < ∞. Under the above assumptions, we have (c.f., [14])  N lim S T /N [φ] = SeT [φ] N →∞

where the superscipt N represents repeated application of the operator, N times.

B. Duals of the semigroup operators This algorithm uses the concept of semiconvex dual. For a function φ which is uniformly semiconvex with constant c [13], the semiconvex dual, a, is given by a(z) = − max ψ(x, z) − φ(x), x∈IRN

(7)

and the dual relationship is given by

IV. P RUNING A LGORITHMS In the above curse-of-dimensionality-free algorithm, at step k, ak is represented as a max-plus sum of quadratics (the b ak ’s). Let us index the elements of this sum by i ∈ Ik (rather than by the sequences {mi }ki=1 ). That is, we have M b aki (z) ak (z) = i∈I k

φ(x) = maxn [ψ(x, z) + a(z)] z∈IR

(8)

2

where ψ(x, z) = (c/2)|x − z| . We may obtain the duals of the Sτm operators, Bτm , and these are also max-plus linear semigroup operators. In particular, they are max-plus integral operators with kernels Bτm (x, z) = − max {ψ(y, x) − Sτm [ψ(·, z)](y)} . y∈IRN

Importantly, note that as the Sτm [ψ(·, z)](y) are quadratic functions, the Bτm are quadratic functions. Each of these is obtained only once, at the outset of the algorithm. C. Dual space propagation and the curse-of-complexity Once the kernels Bτm of the dual semigroup are obtained, one can begin the iteration. One may begin with an initial quadratic (in the dual space), say b − zb) + rb. a0 (z) = b a0 (z) = (1/2)(z − zb)T Q(z Given approximation, ak at step k, one obtains the next iterate from M Z ⊕ ak+1 (z) = Bτm (z, y) ⊗ ak (y) dy m∈M

= max maxn [Bτm (z, y) + ak (y)]. m∈M y∈IR L k If a has the form ak (z) = ak{mi }k (z) k b {mi }k i=1 ∈M i=1 where each b a{mi }ki=1 is a quadratic form, then ak+1 takes the form M b ak+1 (z) ak+1 (z) = {m }k+1 k+1 {mi }k+1 i=1 ∈M

i i=1

where the b ak+1 are also quadratic. Consequently, the {mi }k+1 i=1 computations reduce to obtaining the coefficients of these quadratics at each step, and these computations are analytic (modulo matrix inverses). This is the reason that the computational growth in the space dimension is only cubic. However, note that the number of quadratics comprising the ak grows by a factor of M at each iteration — hence the curse-of-complexity. It has been noted that quite typically, most of the quadratics do not contribute to the value (as they never achieve the maximum at any point, z), and may be pruned wihtout consequence. However, over-pruning, saving only a limited number of quadratics at each iteration, has proven to be an excellent computational technique. We now proceed to discuss the use of semidefinite programming as means of pruning the constituent quadratics, the b ak{m }k+1 . i i=1 Lastly, note that once one has propagated sufficiently far (say k steps), the value function approximation is recovered from ak via (8).

where we let each b aki be given in the form

b aki (z) = b ai = z T Ai z + 2bTi z + ci

where we delete the superscript k for simplicity of notation here and in the sequel. Recall that we are reducing computational cost by pruning quadratics (b ai ) which do not contribute to the solution approximation (not achieving the maximum at any z ∈ IRn ). Further, we may intentionally over-prune, removing those that appear to contribute very little. Consequently, we want to determine whether the pth quadratic contributes to the pointwise maximum. In other words, we need to determine whether there is a region where it is greater than all other quadratics. Fix p ∈ I k . We need to check feasibility of following inequalities. That is, we want find z such that z T Ap z + 2z T bp + cp ≥ z T Ai z + 2z T bi + ci ∀i 6= p (9) Alternatively, we consider the problem: . Minimize G(z, ν) = ν subject to (10) z T (Ai − Ap )z + 2z T (bi − bp ) + (ci − cp ) ≤ ν ∀i 6= p Then, the minimum value of ν, ν¯, is the minimum amount by which the pth quadratic needs to be raised in order to contribute to the max-plus sum. If ν¯ > 0, then pth quadratic does not contribute to the max-plus sum, and hence it can be pruned without consequence. If ν¯ < 0, it implies that quadratic contributes to the max-plus sum, and moreover it is |¯ ν | above the pointwise maximum of all other quadratics at some point. In this case, |¯ ν | can serve as some measure of contribution of the pth quadratic to the value function. A. Pairwise pruning Before undertaking the pruning using semidefinite programming, pairwise pruning is used. This is a simple, fast and effective technique, which checks between all pairs of quadratic basis functions, and prunes those which are completely dominated by another. Let A = Ai1 − Ai2 , b = bi1 − bi2 , c = ci1 − ci2 , and define Q(z) = z T Az + 2bT z + c. Then Q is nonnegative everywhere if and only if the homogeneous quadratic form is nonnegative everywhere. For any t 6= 0, Q(t−1 z) ≥ 0, ∀z implies that Q(z, t) = z T Az + 2tbT z + ct2 ≥ 0, ∀z, t. Hence,  T    t c bT t ≥ 0 ∀z, t z b A z is true if and only if



c b

bT A



 0.

(11)

This pairwise pruning reduces the computational effort of the semidefinite pruning by getting rid of obviously dominated quadratics. B. Shor’s semidefinite relaxation based SDP The problem of evaluating an individual quadratic b ap (z) for pruning, (10), can be rephrased as below. Let qi (z) = b ai (z) − b ap (z) for all i 6= p. Then, b ap can be pruned if min {ν : ν − qi (z) ≥ 0 ∀i} ≥ 0. z,ν

(12)

Introducing the slack variables, ui , (and letting u = {ui }) we see that b ap can be pruned if  min ν : ν − qi (z) − ui 2 = 0 ∀i ≥ 0. (13) z,ν,u

Let fi (z, ν, u) = ν − qi (z) − ui 2 . Since fi assume value zero in the constraint set fi = 0, we can add their λ-weighted combination to the objective function without changing its k value, where λ ∈ IR#I −1 . That is,  min ν : ν − qi (z) − ui 2 = 0 ∀i z,ν,u ) ( X λi fi (z, ν, u) : fi (z, ν, u) = 0 = min ν − z,ν,u

≥ min

z,ν,u

(

i

ν−

X

)

λi fi (z, ν, u)

i

. = ζ(λ).

Note that in last equation, the domain for minimization is no longer the constraint set, but the whole space, i.e., z ∈ k IRn , ν ∈ IR, u ∈ IR#I −1 . Also, we see that ζ(λ) is a lower bound for the minimum in (13). We wish to maximize lower bound ζ(λ) by varying λ. Consequently, the dual problem is max

λ∈IR#I k −1

ζ(λ).

(14)

Next we narrow the search for P the optimal λ by considering the following arguments. If λi 6= 1, then ∀z, u  X X 
λi qi (z) + u2i = −∞. ζ(λ) ≤ min 1 − λi ν + ν

Hence, forPa finite lower bound, it is fruitful to restrict λ such that λi = 1. This makes the objective independent of ν. Now, if λj < 0, for some j, then for all z and all {ui : i 6= j},   X X λi qi (z) + λi u2i + λj u2j  = −∞ ζ(λ) ≤ min  uj

i

i6=j

Hence, we may further restrict the P domain to λi ≥ 0, ∀i. Thus, λ lies within the simplex S, λi = 1 and λi ≥ 0. Note also that this makes the objective function independent of ui , since the above minimum with respect to ui is always achieved at ui = 0. Consequently, X λi qi (z) ∀λ ∈ S. (15) ζ(λ) = min z

i

Thus, dual problem (14) can be reposed as X λi qi (z) max ζ(λ) = max ζ(λ) = max min λ

λ∈S

=

z

λ∈S

max

λ∈S,ζ∈IR

(

ζ:ζ≤

i

X

λi qi (z) ∀z

i

)

.(16)

Suppose that quadratics qi (z) are specified by parameters A¯i , ¯bi and c¯i . Then, using linear superposition and result (11), P the condition ζ ≤ i λi qi (z) ∀z from (16) can be posed as following linear matrix inequality.     P P λ c¯ − ζ Pi λi¯bTi  0 (17) max ζ : Pi i¯i ¯ λ∈S,ζ i λi Ai i λi bi

Note that if the maximal ζ value, ζ¯ satisfies ζ¯ ≥ 0, then ¯ and thus the minimal ν value in (12), ν¯, is greater than ζ, positive. This indicates that the quadratic, b ap , needs to be raised by ν¯ units before it can contribute to the max-plus sum. Hence it can be pruned. If ζ¯ < 0, then its “prunability” is not conclusive, as ζ¯ is a lower bound for ν¯, and the gap is not guaranteed to be zero. Nevertheless, it does give us a working indication of the importance of the quadratic, since ζ¯ < 0 indicates that the pth quadratic has to come down by ¯ units, before it is dominated by the convex hull at least |ζ| of the remaining quadratics. An additional way to develop intuition for result (17) is as follows. The above test evaluates the pth quadratic b ap , which can be pruned if ζ = 0 satisfies the inequality in (16). Lets now substitute qi (z) = b ai (z) − b ap (z) in (16). We see that b ap can be pruned if for some λ ∈ S 0≤

K−1 X

λi qi (z) =

i

or equivalently,  cp bp

X i6=p

bTp Ap



∀ z ∈ IRn ,

λi b ai (z) − b ap (z)



X i6=p

λi



ci bi

bTi Ai



.

(18)

Thus, if the convex hull of remaining quadratics intersects the semidefinite cone of quadratics greater than b ap , then b ap can be pruned. In our context, the substantial benefits of Shor’s relaxation are that it accommodates both concave and convex quadratics, calculates an importance/quality metric and gives a clean analytical framework for analysis of pruning. One efficient way to implement the above scheme is to prune all but the quadratics lying on the periphery of the convex hull, and then to evaluate each vertex quadratic one by one by the above semidefinite program (17). A seeming drawback of this algorithm is gives the optimum in λ space, and not z, making it difficult to assign locally weighted importance. That is, we might choose to assign more importance to the quadratics which are vital near origin than those far from it. The next method helps us do precisely the same.

C. Schur complement-based SDP In this method, we again aim to solve (10) without resorting to dual relaxation schemes. It is based on following well known Schur’s complement lemma, which can be found in [4] Lemma 4.1: Let   B CT A= C D be a symmetric matrix with k × k block B, and l × l block D. Assume that B is positive definite. Then A is positive semi-definite if and only if D − CB −1 C T is positive semidefinite. Now starting form (10), if we define A¯i = Ai − Ap , ¯bi = bi − bp and c¯i = ci − cp , and assume that A¯i ≻ 0, then (10) can be rephrased as: . Minimize G(z, ν) = ν subject to (19) T¯ T ¯ z Ai z + 2z bi + c¯i ≤ ν ∀i 6= p. By Lemma 4.1, this is equivalent to the following linear matrix inequality (LMI).   −(2z T ¯bi + c¯i ) + ν z T  0 ∀i 6= p z A¯−1 i

Note that if A¯i ≻ 0, ∀i 6= p, then (19) is equivalent to the LMI given by     −2z T ¯bi − c¯i + ν z T min ν :  0 ∀i 6= p . z,ν z A¯−1 i (20) If the minimal ν value, ν¯ satisfies ν¯ ≥ 0, then the quadratic b ap can be pruned without consequence. If it is negative, then |¯ ν | serves as a metric of importance of b ap (z), as it says that at some point z, b ap is at least |¯ ν | higher than all other quadratics. ν¯ and the z at which it occurs can be combined to create a new importance metric catered to region of interest. For example, |¯ ν |/(1 + |z|2) ensures tighter error bounds near origin than far from it. However, the assumption A¯i ≻ 0 ∀i 6= p is overly restrictive, as it only enables comparison with quadratics more convex than b ap , as we simply drop the inequalities for which this does not hold. This also ensures conservative pruning. To do a tighter pruning, it is advisable to rank the quadratics in the increasing order of lowest eigenvalues, before pruning. If A¯i  0 instead of If A¯i ≻ 0, then the above LMI can be modified to use a linear transformation of A¯i using its nonzero singular values. V. C OMPUTATIONAL COMPLEXITY Since our aim is to reduce the curse-of-complexity without letting go of our freedom from the curse-of-dimensionality, it is worthwhile to discuss the computational overhead involved in these pruning methods. They are polynomial in both dimensionality and the number of quadratic functions. In particular, we retain our freedom from the curse-ofdimensionality.

A generic semi-definite program P is given by   N   X ηj Aj ≥ 0, ||η||2 ≤ R P0 = min cT η : A0 +  η∈IRN  j=1

f diagonal where the Aj are symmetric matrices with M f blocks of size ki × ki , i = 1, . . . M . We say that η ε is an ε–optimal solution if kη ε k2 ≤ R, A0 +

N X

ηjε Aj ≥ −ǫI, cT η ε ≤ P0 + ǫ.

j=1

Reference [4] derives the computational complexity of obtaining such an η ε .  1/2 f M X C(P, ǫ) = O(1)1 + ki  (21) 

1

·N N 2 + N

f M X

ki2 +

1

f M X 1



ki3  D(P, ǫ)

where D(P, ǫ) depends on the specific problem data (as indicated by the P in the argument as well as ǫ. Using this expression, we can obtain an upper bound for arithmetic complexity of the pruning algorithms. Assuming the worst case scenario where no quadratic gets pruned, we find the complexity of testing one quadratic for pruning as follows. For the Schur complement pruning of (20), we have η = [z T ν]T . Hence N = n + 1. Further, the number of block . f = I, and k1 = k2 = . . . = LMIs is I = #Ik − 1. Thus M km = n + 1. Substituting these into (21), we obtain the complexity as C(P, ǫ) = O(1) (1 + I(n + 1))1/2 (n + 1)3 · [1 + 2(n + 1)(I − 1)] D(P, ǫ),

(22)

and we see that the complexity grows as n4.5 and I 1.5 . For the Shor semidefinite relaxation approach, we apply a similar analysis, and find that the complexity grows as n3.5 and I 3.5 . Note that these computational cost growth rates as a function of space dimension, n, are poorer than the cubic growth of the basic algorithm. However, they are still tremendously helpful for real problems, and the curse-ofdimensionality is still very far off. To find worst case complexity of one stage of the pruning, we beed to multiply above figures by I. VI. I MPORTANCE - BASED

OVER - PRUNING

Both the Shor’s relaxation and the Schur complement pruning schemes generate importance metrics for unpruned quadratics. This is very useful in controlling the complexity growth according to an error tolerance, region of interest, or computational limits. Currently, we set a bound L(k) for the number of quadratics which will remain after over-pruning at the k th step. The algorithm keeps only the L(k) quadratics, b aki with the highest values according to this measure of importance. However, there does not yet exist a theory which

allows us to map the importance measure of a quadratic at step k to an error bound in the approximation at the terminal step. Further, a bound on the errors induced by over-pruning of multiple quadratics is obviously a necessary additional step. Regardless, application of these methods has proven extremely fruitful, as can be seen in the example below.

0.15 0.1 0.05

VII. N UMERICAL E XAMPLE The example we choose to showcase here is six dimensional, and involves the maximum of two constituent Hamiltonians. It is generated by appending three 2-D systems, and then adding interaction terms coupling the dimensions. Let us define following building blocks     −1 .5 −1 .5 Aa = , Ab = , .1 −1 .3 −1   −1 .5 Ac = 0 −1     1.5 .2 0.216 −.008 Da = , Σa = , .2 1.5 −.008 .216

and let Db = 1.4Da and Σb = 0.4Σa . The parameter adjusting the interaction in the dynamics across the dimensions is γ. Consequently the dynamics for the Hamiltonians are as follows.     Aa γI γI Aa γI 0 A1 =  γI Aa 0  , A2 =  γI Ab 0  γI 0 Aa 0 0 Ac 1

0 −0.05 2 2

1 1 0 0 −1 −2

Fig. 1.

Back substitution error along axis 1 and 2

0.06 0.04 0.02 0 −0.02 −0.04 −0.06 2 1

2

2

=

l12

1

0

Σ = diag (Σa , Σa , Σa ) , Σ = diag (Σb , Σb , Σb ) l11

−1

0 −1

−1 −2

=0

−2

The parameters of the cost functions are D1 = diag (Da , Da , Da ) , D1 = diag (Db , Db , Db )

Fig. 2.

Back substitution error along axis 3 and 4

l21 = l22 = 0 For this example, with γ = 0.1 and time-discretization step-size τ = 0.2, propagation was carried out with the Schur complement pruning. The overpruning threshold was set to L(k) = 2k for k th step. That is, a maximum of 2k quadratics, b aki , were retained at the k th step. In this test, 40 iteration steps were carried out in 90 seconds, yielding a rather accurate solution in a compact domain in all six dimensions. This computation-time is for an Apple mac desktop, from roughly 2005. Slices of statistics for this value function along the 1-2, the 3-4, and the 4-5 coordinate axes are shown in the accompanying figures. The backsubstitution error depends on propagation as well as time discretization. Theoretical results [14] indicate error bounds of the form ε(1 + |x|2 ) (over the entire space) where ε ↓ 0 as the number of propagation steps goes to infinity and time-discretization go to zero, with the required relative rates being discussed in the reference. However, we find that, looking over only a compact sub-domain, the component of the error due to the number of propagation steps goes to zero after a finite number of steps, leaving the error due to time-discretization. Consequently, the error reduction on a compact sub-domain saturates after a sufficient number of iterations are performed.

0.15 0.1 0.05 0 −0.05 2 1.5 1

2

0.5 1

0 −0.5

0 −1 −1

−1.5 −2

Fig. 3.

Back substitution error along axes 5 and 6

VIII. ACKNOWLEDGMENT The authors would like to thank Dr. Stephen Gaubert who pointed out the appropriateness of this semidefinite programming approach.

2

[3] 1.5

[4]

1

0.5

[5] 0

[6]

−0.5

−1

[7] −1.5

−2 −2

0

−1

−1.5

−2 2

0

−0.5

0.5

1

1.5

2

[8] [9]

Fig. 4.

Partial w.r.t. x4 on slice along axes 2 and 3 at origin

[10] [11] [12]

0.03

0.02

[13]

0.01

0

[14]

−0.01

−0.02

[15]

−0.03 2 1

[16]

0 −1 1.5

−2 2

0

0.5

1

−0.5

−1

−1.5

−2

[17] Fig. 5.

Partial w.r.t. x6 on slice along axes 1 and 2 at origin

0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 2 1 0 −1 −2

−1.5

−1

−0.5

0

0.5

1

1.5

2

Fig. 6. Partial w.r.t. x4 on slice along axes 1 and 2 at point x4 = 0.5, x3 = x5 = x6 = 0

R EFERENCES [1] M. Akian, S. Gaubert, and A. Lakhoua, “The max-plus finite element method for optimal control problems: further approximation results”, Proc. 44th IEEE Conf. on Dec. and Control (joint with ECC’05), Sevilla (2005). [2] M. Akian, S. Gaubert and A. Lakhoua, “A max-plus finite element

method for solving finite horizon determinsitic optimal control problems”, Proc. 16th MTNS (2004). F.L. Baccelli, G. Cohen, G.J. Olsder and J.-P. Quadrat, Synchronization and Linearity, John Wiley (1992). A. Ben-Tal and A. Nemirovski, Lectures on Modern Convex Optimization: analysis, algorithms and engineering applications MPS-SIAM Series on Optimization R.A. Cuninghame-Green, Minimax Algebra, Lecture Notes in Economics and Mathematical Systems 166, Springer (1979). P. Dupuis and A. Szpiro, Second order numerical methods for first order Hamilton-Jacobi equations, SIAM J. on Num. Analysis, 40 (2002), pp. 113–118. M. Falcone and R. Ferretti, Convergence analysis for a class of highorder semi-lagrangian advection schemes, SIAM J. Num. Analysis, 35 (1998), pp. 909–940. W.H. Fleming and W.M. McEneaney, “A max-plus based algorithm for an HJB equation of nonlinear filtering”, SIAM J. Control and Optim., 38 (2000), 683–710. V.N. Kolokoltsov and V.P. Maslov, Idempotent Analysis and Its Applications, Kluwer (1997). G.L. Litvinov, V.P. Maslov and G.B. Shpiz, “Idempotent Functional Analysis: An Algebraic Approach”, Mathematical Notes, 69, No. 5 (2001), 696–729. V.P. Maslov, “On a new principle of superposition for optimization problems”, Russian Math. Surveys, 42 (1987) 43–54. W.M. McEneaney, ‘Curse-of-Dimensionality Free Method for Bellman PDEs with Semiconvex Hamiltonians,” Proc. 45th IEEE Conf. on Dec. and Control, San Diego (2006). W.M. McEneaney, Max-Plus Methods in Nonlinear Control and Estimation, Birkhauser (2006). W.M. McEneaney and J.L. Kluberg, “Convergence Rate for a Curseof-Dimensionality-Free Method for HJB PDEs Represented as Maxima fo Quadratic Forms,” Submitted to SIAM J. Control and Optim. W.M. McEneaney, “A Curse-of-Dimensionality-Free Numerical Method for Solution of Certain HJB PDEs,” SIAM J. Control and Optim., (to appear). W.M. McEneaney, “Max–Plus Eigenvector Representations for Solution of Nonlinear H∞ Problems: Basic Concepts”, IEEE Trans. Auto. Control 48, (2003) 1150–1163. R.T. Rockafellar and R.J. Wets, Variational Analysis, Springer–Verlag (1997).