Suboptimal feedback control of PDEs by solving Hamilton-Jacobi - - PowerPoint PPT Presentation

Suboptimal feedback control of PDEs by solving Hamilton-Jacobi Bellman equations

• n sparse grids

Jochen Garcke joint work with Axel Kr¨

• ner, INRIA Saclay and CMAP, Ecole Polytechnique

Ilja Kalmykov, Universit¨ at Bonn (who is looking for a PhD position)

Department Numerical Data-Driven Prediction Institute for Numerical Simulation

Optimal control Sparse Grids Numerical Results Higher Order Methods

Outline

1 Optimal control 2 Sparse Grids 3 Numerical Results 4 Higher Order Methods

Optimal control Sparse Grids Numerical Results Higher Order Methods

Optimal Control of Low Dim. Approx. of PDEs

continuous problem described by PDE semi-discretization in space model reduction semi-discrete control problem HJB-equation sparse grids feedback operator

Optimal control Sparse Grids Numerical Results Higher Order Methods

Optimal Control of Low Dim. Approx. of PDEs

continuous problem described by PDE semi-discretization in space model reduction semi-discrete control problem HJB-equation sparse grids feedback operator

Optimal control Sparse Grids Numerical Results Higher Order Methods

Optimal Feedback Control

• ptimal feedback control of a dynamical system

         min

u∈Uad J(u) =

T l(y(t), u(t)) dt, s.t. ˙ y(t) = f (y(t), u(t)), t > 0 y(0) = y0

state y(t) ∈ ❘d, initial state y0 ∈ ❘d control u(t) ∈ Um ⊂ ❘m (often called action) Lipschitz continuous dynamics f : ❘d × ❘m → ❘d running cost with polynomial growth l : ❘d × ❘m → ❘ set of admissible controls Uad = {u ∈ L2([0, T]; Um) | Um ⊂ ❘m compact}

aim: feedback law u∗ = K(t, y∗(t))

Optimal control Sparse Grids Numerical Results Higher Order Methods

Semi-Lagrangian scheme

from DPP one derives semi-Lagrangian scheme vk(x) = min

u∈U

• vk+1 (yx(∆t)) + ∆t · l(x, u)
• (SL)

for k = N − 1, ..., 0 with vN(x) = 0, time step ∆t > 0, x ∈ ❘d yx(∆t) state obtained by time discretization scheme from x

evaluation of the right hand side in (SL), either:

(potentially expensive) comparison over a finite set Ufinite ⊂ U in linear quadratic case, i.e. f (x, u) = Ax + Bu, A ∈ ❘d×d, B ∈ ❘d×m l(x, u) = 1 2

• xTMx + uTRu
• ,

M ∈ ❘d×d, R ∈ ❘m×m the optional feedback control is given by (needing approx. for ∇v) u∗(t) = PU

• −R−1BT∇v(y ∗(t), t)

Optimal control Sparse Grids Numerical Results Higher Order Methods

Optimal Control of Low Dim. Approx. of PDEs

continuous problem described by PDE semi-discretization in space model reduction semi-discrete control problem HJB-equation sparse grids feedback operator

Optimal control Sparse Grids Numerical Results Higher Order Methods

Application to a PDE Control Problem

wave equation (following Kr¨

• ner, Kunisch, Zidani (2015))

     ˆ ytt − c∆ˆ y = Bu in (0, T) × Ω, ˆ y(0) = ˆ y0, ˆ yt(0) = ˆ y1 in Ω, ˆ y = 0

• n (0, T) × ∂Ω

initial state and velocity ˆ y0 ∈ H1(Ω) and ˆ y1 ∈ L2(Ω) control operator B := (sin(πx), . . . , sin(mπy)) (can be generalized)

formulate as first order system in time with y1 = y, y2 = ˙ y y1

t = y2,

y2

t − ∆y1 = Bu,

y1(0) = ˆ y0, y2(0) = ˆ y1, which we can write as yt + Ay = (0, Bu)T, y(0) = y0 with y = (y1, y2), y0 = (ˆ y0, ˆ y1)T ∈ Y 1, and A =

• Id

−c∆

Optimal control Sparse Grids Numerical Results Higher Order Methods

Semi-Discrete Formulation of Control Problem

semi-discrete formulation of control problem by method of lines for a given basis b := (ϕ1, . . . , ϕd), d ∈ N, with ϕi : Ω → ❘ we define A := ((∇ϕi(x), ∇ϕj(x))i,j=1,...,d) (stiffness matrix) M := ((ϕi(x), ϕj(x))i,j=1,...,d) (mass matrix) in our numerical examples we later choose ϕi(x) := sin(iπx), i = 1, . . . , d, and obtain A = diag((1/2(iπ)2)i=1,...,d), M = diag((1/2)i=1,...,d)

Optimal control Sparse Grids Numerical Results Higher Order Methods

Resulting Semi-Discrete PDE Control Problem

this now fits into control formulation from beginning          min

u∈Uad J(u) =

T l(y(t), u(t)) dt, s.t. ˙ y(t) = f (y(t), u(t)), t > 0 y(0) = y0 for initial value y0 ∈ ❘2d and t ∈ [0, T). dynamics f (x, u) = Ax + Bu with running cost l(x, u) := βxxT

1 Mx1 + βuuTu

if value function is differentiable in x, t optimal feedback control u∗(y∗(t), t) = PU

• − 1

βu BT∇xv(y∗(t), t)

• where PU projection on set of admissible controls

Hamilton-Jacobi Bellman equation with dimension 2d

Optimal control Sparse Grids Numerical Results Higher Order Methods

Optimal Control of Low Dim. Approx. of PDEs

continuous problem described by PDE semi-discretization in space model reduction semi-discrete problem HJB-equation sparse grids feedback operator

Optimal control Sparse Grids Numerical Results Higher Order Methods

Interpolation with Hierarchical Basis

φ3,1 φ3,5 φ3,3 φ3,2 φ3,4 φ3,6 φ3,7

nodal basis V1 ⊂ V2 ⊂ V3

φ3,7 φ3,1 φ3,3 φ2,1 φ1,1 φ3,5 φ2,3

hierarchical basis V3 = W3 W2 V1

Optimal control Sparse Grids Numerical Results Higher Order Methods

• Hier. Basis Functions in Higher Dimensions

d-dimensional piecewise d-linear functions φl,j(x) :=

d

• t=1

φlt,jt(xt) hierarchical difference space Wl (et is t-th unit vector) Wl := Vl \

d

• t=1

Vl−et,

• hier. diff. space represented by Wl = span{φl,j | j ∈ Bl} with

Bl :=

• j ∈ Nd
• jt = 1, . . . , 2lt − 1,

jt odd, t = 1, . . . , d, if lt > 1, jt = 0, 1, 2, t = 1, . . . , d, if lt = 1

• .

full grid space in hierarchical basis V s

n :=

• |l|∞≤n

Wl

Optimal control Sparse Grids Numerical Results Higher Order Methods

Hierarchical Subspaces Wl

W4,1 W4,2 W4,3 W4,4 W3,1 W3,2 W3,3 W3,4 W2,1 W2,2 W2,3 W2,4 W1,1 W1,2 W1,3 W1,4

Optimal control Sparse Grids Numerical Results Higher Order Methods

Sparse Grids

we define the sparse grid function space V s

n ⊂ Vn as

V s

n :=

• |l|1≤n+d−1

Wl every f ∈ V s

n can now be represented as

f s

n (x) =

• |l|1≤n+d−1
• j∈Bl

αl,jφl,j(x) approximation property in H2

mix

||f − f s

n ||2 = O(h2 n log(h−1 n )d−1)

sparse grid needs O(h−1

n (log(h−1 n ))d−1) points

Optimal control Sparse Grids Numerical Results Higher Order Methods

Sparse Grids in two and three dimensions

Optimal control Sparse Grids Numerical Results Higher Order Methods

Problems with Sparse Grids: Monotonicity

interpolation of peaked Gaussian fct. with sparse grid n = 2 f (x1, x2) := exp

• −100(x1 − 0.5)2

∗ exp

• −100(x2 − 0.5)2

0-level set of interpolant is pink in the right picture sparse grid interpolation does not preserve positivity

Optimal control Sparse Grids Numerical Results Higher Order Methods

Spatially Adaptive Sparse Grids

to approximate functions which either

do not fulfil smoothness condition at all or strongly vary due to finite but locally large derivatives

adaptive refinement may be used start with a regular grid of level 2 (left) to populate index set I refine one grid point by creating all children (middle) to keep grid consistent, missing parents are created (right) usually hierarchical surplus αl,j is used as refinement indicator

Optimal control Sparse Grids Numerical Results Higher Order Methods

Basic Sparse Grid SL Scheme

evaluate for x ∈ QI, ∆t > 0, K = T/∆t, k = K − 1, . . . , 0,    vk(x) = min

u∈U

• ∆tl(x, u) + vk+1(yx(∆t))
• ,

vk(x) = 0 yx(∆t) state obtained by time discretization scheme from x Algorithm 1: Adaptive SL-SG scheme Data: refinement constant ε, coarsening constant η Result: sequence of adaptive sparse grid solutions vk ∈ VI(k) initialize I(K) for k = K − 1, . . . , 0 do ⊲ iterate in time with ∆t = T/K initialize I(k − 1) with I(k) adaptively interpolate minu∈U

• vk(yx(∆t)) + ∆tl(x, u)
• ⊲ compute vk−1

coarsen vk−1 ∈ VI(k−1) see Bokanowski, G., Griebel, and Klompmaker (2013)

Optimal control Sparse Grids Numerical Results Higher Order Methods

Further Discretization Aspects

to determine minimizing control within the SL-scheme we use

minimization by comparison over finite subset Uσ ⊂ U or gradient of the value function un(x, s) = PU

• − 1

βu BT∇hv n+1(x, s)

• per finite differences
• ∇hv n+1(x, s)
• i := v n+1(x + h · ei, s) − v n+1(x − h · ei, s)

2h

for computation of yx(∆t) we use second order Heun scheme in our experiments we focus on discretization error in space, while using time resolution which is “good enough” reference solution vr in 2D computed with a higher order finite difference code on a uniform mesh by an ENO scheme compute reference trajectories yr in state space and ur in control space using a Riccati approach

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Simplified Semi-Discrete Wave Equation

dynamics f (x, u) = Ax + Bu with running cost l(x, u) := βxxT

1 Mx1 + βuuTu

example based on harmonic oscillator βx = 2, βu = 0.1, T = 1, ∆t = 0.01, A = 1 −1

• ,

B = 1

• ,

initial data x ∈ ❘2, domain Q = [−1, 1]2, U = [−3.5, 3.5].

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Convergence of value function

10−4 10−3 10−2 10−1 10−3 10−2 10−1 0.47 0.53 ε v − v∗L2 normal hat (gradient) fold out hat (gradient) normal hat (compare) fold out hat (compare)

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Convergence of value function

101 102 103 104 10−3 10−2 10−1 −0.96 −0.57 nodes (end) v − v∗L2 normal hat (adaptive) fold out hat (adaptive) normal hat (regular) fold out hat (regular)

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Value function

1.0 0.5 0.0 0.5 1.0 1.0 0.5 0.0 0.5 1.0 0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

(a) ε = 1.95−4

1.0 0.5 0.0 0.5 1.0 1.0 0.5 0.0 0.5 1.0 0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6

(b) ε = 1.95−4

Figure: adaptive sparse grid with normal and fold out hat functions

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Convergence in the Trajectory

10−4 10−3 10−2 10−1 10−1 100 101 0.51 0.49 ε y − y∗L∞ normal hat (gradient) fold out hat (gradient) normal hat (compare) fold out hat (compare)

(a) L∞ error in the trajectory vs. ε

101 102 103 104 10−1 100 101 102 −0.6 −0.62 nodes y − y∗L∞ normal hat (level) fold out hat (level) normal hat (eps) fold out hat (eps)

(b) L∞ error in the trajectory vs. nodes

Figure: Comparing the scheme on adaptive sparse grids using normal hat and fold out basis functions

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D: Convergence in the Control

10−4 10−3 10−2 10−1 10−3 10−2 10−1 100 0.52 0.56 ε u − u∗L∞ normal hat (gradient) fold out hat (gradient) normal hat (compare) fold out hat (compare)

(a) L∞ error in the control vs. ε

101 102 103 104 10−3 10−2 10−1 100 101 −0.63 −0.62 nodes u − u∗L∞ normal hat (level) fold out hat (level) normal hat (eps) fold out hat (eps)

(b) L∞ error in the control vs. nodes

Figure: Comparing the scheme using the gradient approach on adaptive and regular sparse grids using normal hat and fold out basis functions

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example: Semi-Discrete Wave Equation (4D)

reminder: wave equation as first order system in time ˙ y(t) = fw(y(t), u(t)), t > 0, y(0) = y0 with dynamics fw : ❘2d × ❘m → ❘, fw(x, u) := Ax + Bu where A :=

• Id

−cM−1A

• ,

B := b

• ,

b ∈ ❘m×d, y0 ∈ ❘2d cost lw(x, u) := βxxT

1 Mx1 + βuuTu

consider setup βx = 2, βu = 0.1, T = 4, ∆t = 0.01, c = 0.05.

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example: Convergence in the Control

10−4 10−3 10−2 10−1 10−2 10−1 100 101 0.65 0.63 ε u − u∗L∞ normal hat fold out hat

(a) L∞ error in the control vs. ε

101 102 103 104 105 10−3 10−2 10−1 100 101 102 −0.53 −1.41 nodes u − u∗L∞ normal hat (eps) fold out hat (eps) fold out hat (level)

(b) L∞ error in the control vs. nodes

Figure: Error in the control comparing the scheme based on sparse grids using normal hat and fold out basis functions, regular and adaptive sparse grids (4D).

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example: semi-discrete wave equation (4D)

time 1 2 3 4

• 0.1

0.1 0.2 0.3 0.4 SL-SG Riccati

(a) y1

time 1 2 3 4

• 0.4
• 0.3
• 0.2
• 0.1

0.1 SL-SG Riccati

(b) y3

time 1 2 3 4

• 1.5
• 1
• 0.5

0.5 SL-SG Riccati

(c) u1

time 1 2 3 4

• 0.2

0.2 0.4 0.6 SL-SG Riccati

(d) y2

time

1 2 3 4

• 0.8
• 0.6
• 0.4
• 0.2

0.2

SL-SG Riccati

(e) y4

time 1 2 3 4

• 1.5
• 1
• 0.5

0.5 1 SL-SG Riccati

(f) u2

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example: semi-discrete wave equation (6D)

10−4 10−3 10−2 10−1 10−1 100 101 0.62 ε u − u∗L∞ fold out hat (gradient)

(a) error in the control vs. ε

102 103 104 10−1 100 nodes u − u∗L∞ fold out hat (gradient)

(b) error in the control vs. nodes

Figure: The adaptive sparse grids scheme using fold out basis functions (6D). Need to decrease time step to ∆t = 0.0025 (larger entries in stiffness matrix).

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example: semi-discrete wave equation (8D)

10−4 10−3 10−2 10−1 10−1 100 101 0.57 ε u − u∗L∞ fold out hat (gradient)

(a) error in the control vs. ε

102 103 104 10−1 100 nodes u − u∗L∞ fold out hat (gradient)

(b) error in the control vs. nodes

Figure: The adaptive sparse grids scheme using fold out basis functions (8D). Need to decrease time step to ∆t = 0.00125 (larger entries in stiffness matrix).

Optimal control Sparse Grids Numerical Results Higher Order Methods

Bilinear System of Schr¨

• dinger type

we write a Schr¨

• dinger type equation as a real-valued system

fs : ❘2d × ❘ → ❘, fs(x, u) := Ax + uBy with A =

• cM−1A

−cM−1A 0,

• , B = M−1 ˆ

B, M = M M

• , ˆ

B ∈ ❘2d×2d as before discrete cost functional is given by J(u, y) := T l(y(t), u(t)) dt running cost ls(x, u) := βxxTMx + βuuTu, where x = (x1, x2) ∈ ❘2d, u ∈ U

Optimal control Sparse Grids Numerical Results Higher Order Methods

Example 2D

setup A = 0.5 −0.5

• ,

B = −0.5 0.5

• ,

M = 1 1

• with T = 1, ∆t = 1/500, βx = 2, βu = 0.1, U = [−4, 4]

in this case the control bounds are active

Optimal control Sparse Grids Numerical Results Higher Order Methods

2D-Bilinear Schr¨

• dinger Type Problem

101 102 103 104 10−2 10−1 100 −0.43 −0.59 nodes (end) error normal hat fold out hat

(a) error in the value function vs. nodes

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 1 1 2 3 4 5

(b) fold-out basis functions and ε = 7.78· 10−4

Figure: Bilinear equation: Sparse grid (2D).

Optimal control Sparse Grids Numerical Results Higher Order Methods

Higher Order Methods

we aimed to investigate (adaptive) time stepping in particular in combination with adaptive sparse grids but got side-tracked into investigating higher order time stepping which resulted into also looking at higher order (B-Splines) sparse grids we can give so far just preliminary results

Optimal control Sparse Grids Numerical Results Higher Order Methods

Higher Order Runge-Kutta Methods

look at our control formulation from beginning      min

u∈Uad J(u) =

T l(y(t), u(t)) dt, s.t. ˙ y(t) = f (y(t), u(t)), t > 0, y(0) = y0 we need time discretisation for y(t) and quadrature rule to evaluate minu∈Uad J(u) with higher order RK-scheme this could look like (misusing notation) vk(x) = min

u1,u2,u3,...

• ∆t
• i

cil(xτi, ui) + vk+1(ˆ yx(∆t))

• where ˆ

yx takes into account actions ui and RK-scheme with intermediate steps τi therefore RK4 would have O(d4

c ) complexity if complexity for

computing one control is dc see also Falcone, Ferretti (1994).

Optimal control Sparse Grids Numerical Results Higher Order Methods

Simplified Semi-Discrete Wave - Using Heun

10−3 10−2 10−1 100 10−3 10−2 10−1 ∆t v − v⋆L2 l = 4 and 2 control points l = 6 and 2 control points l = 8 and 2 control points l = 4 and 1 control point l = 6 and 1 control point l = 8 and 1 control point

(a) error in value function vs. time steps

102 103 104 105 106 10−3 10−2 10−1 Function evaluations v − v⋆L2 l = 4 and 2 control points l = 6 and 2 control points l = 8 and 2 control points l = 4 and 1 control point l = 6 and 1 control point l = 8 and 1 control point

(b) error in value function vs. evaluations

Figure: Error to the FD reference solution for Heun time integrator and different control strategies. Evaluation of control using compare approach with 40 points.

Optimal control Sparse Grids Numerical Results Higher Order Methods

Structure Preserving Runge-Kutta Methods

• bserve that for explicit or diagonally implicit RK schemes, the last

element un

s of discrete control vector affects only the last RK-stage ks

in fact, condition on RK-coefficients can be formulated, so that control

• ptimization can be applied separately to each stage

class of RK methods, which fulfils this are diagonally implicit symplectic Runge-Kutta (DISRK) schemes “re-use” implicit collocation RK-scheme for quadrature

• riginally developed for long time integration of Hamiltonian systems

are equivalent to composition of implicit midpoint schemes Ψh = Φbsh ◦ · · · ◦ Φb2h ◦ Φb1h we can construct a SL scheme, which has O(dc · s) complexity for the minimization problem problem for s > 1 DISRK goes backwards in time for some steps !?

Optimal control Sparse Grids Numerical Results Higher Order Methods

Simplified Semi-Discrete Wave - Revisited

10−3 10−2 10−1 100 10−3 10−2 10−1 ∆t v − v⋆L2 SG level 2 SG level 4 SG level 6 SG level 8 SG level 10 SG level 12

(a) Heun, d2

c cost

10−3 10−2 10−1 100 10−3 10−2 10−1 ∆t v − v⋆L2 SG level 2 SG level 4 SG level 6 SG level 8 SG level 10 SG level 12

(b) implicit midpoint, 2 · dc cost

Figure: Error against FD reference, T = 1.0

Optimal control Sparse Grids Numerical Results Higher Order Methods

Simplified Semi-Discrete Wave - Revisited

10−4 10−3 10−2 10−1 10−7 10−6 10−5 10−4 ∆t v − v⋆L2 l = 6 l = 8 l = 10 l = 12 l = 14

(a) implicit midpoint / DISRK1

10−4 10−3 10−2 10−1 10−8 10−7 10−6 10−5 10−4 2.07 ∆t v − v⋆L2 l = 10 eval. on [−1.0, 1.0]2 l = 10 eval. on [−0.5, 0.5]2 l = 12 eval. on [−1.0, 1.0]2 l = 12 eval. on [−0.5, 0.5]2

(b) implicit midpoint / DISRK1

Error against Riccati-Solution, T = 0.1 measured on full domain [−1, 1]2 and smaller [−0.5, 0.5]2 similarly we observed for DISRK3 an order 3 with error to 10−6

Optimal control Sparse Grids Numerical Results Higher Order Methods

Simplified Semi-Discrete Wave Higher Order Discretization

10−4 10−3 10−2 10−1 10−11 10−10 10−9 10−8 10−7 10−6 10−5 ∆t v − v⋆L2 l = 10 modified linear l = 4 mod. BSpline deg. 3 l = 5 mod. BSpline deg. 3 l = 6 mod. BSpline deg. 3 l = 7 mod. BSpline deg. 3

(c) linear SG vs. B-Splines -SG

results using Sparse Grids with B-Splines using DISRK1, measure error on smaller domain

Optimal control Sparse Grids Numerical Results Higher Order Methods

Conclusion

• ptimal feedback control PDE problems lead to HJB equations

use model reduction to reduce complexity SL-scheme on sparse grids for HJB equations value function in coefficients of spectral basis can be ”smooth“ higher order time stepping schemes and higher order discretisation can be effective for smooth problems sparse grids with B-splines converge ”nicer“, but have higher computing costs due to need to solve a linear equation system are there are other composite methods, which are better suitable ?

Optimal control Sparse Grids Numerical Results Higher Order Methods

Key References

• O. Bokanowski, J. G., M. Griebel, and I. Klompmaker An adaptive

sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations, J. of Sci. Comp., (2013)

• M. Falcone, R. Ferretti, Discrete time high-order schemes for viscosity

solutions of Hamilton-Jacobi-Bellman equations, Numerische

• Mathematik. (1994).
• J. G., A. Kr¨
• ner, Suboptimal feedback control of PDEs by solving

HJB equations on adaptive sparse grids, J. of Sci. Comp., (2016).

• A. Kr¨
• ner, K. Kunisch, H. Zidani, Optimal feedback control for the

undamped wave equation equation, ESAIM: Control Optim. Calc. Var.

(2015).