Toward Feedback Control of Densities in Nonlinear Systems Abhishek - - PowerPoint PPT Presentation

Toward Feedback Control of Densities in Nonlinear Systems

Abhishek Halder

Department of Applied Mathematics University of California, Santa Cruz Santa Cruz, CA 95064

Joint work with Kenneth F. Caluya (UC Santa Cruz)

CDC 2019, Workshop: Uncertainty Synthesis, Nice

Overarching Theme

Systems-control theory for densities

What is density?

n 
 l

Probability Density Fn.

x(t) ∈ ! " x y θ # $ ∈ X ≡ R2 × S1

n 
 l

• l


 l


 l


 l


 l


 l


 l

Probability Density Fn.

x(t) ∈ ! " x y θ # $ ∈ X ≡ R2 × S1

ρ (x, t) : X × [0, ∞) $→ R≥0

!

X ρ dx = 1

for all t ∈ [0, ∞)

n 
 l

• l


 l


 l


 l


 l


 l


 l

n 
 l

• l


 l


 l


 l


 l


 l


 l

Probability Density Fn.

x(t) ∈ ! " x y θ # $ ∈ X ≡ R2 × S1

ρ (x, t) : X × [0, ∞) $→ R≥0

!

X ρ dx = 1

for all t ∈ [0, ∞)

Population Density Fn.

Why care about densities?

Prediction Problem

Process model Initial conditions Parameters Process noise State density ρ(x(t), t)

Trajectory flow:

dx(t) = f(x, t) dt + g(x, t) dw(t), dw(t) ∼ N (0, Qdt)

Density flow:

∂ρ ∂t = LFP(ρ) := −∇ · (ρf) + 1 2

n

∑

i,j=1

∂2 ∂xi∂xj !" gQg⊤#

ij ρ

$

Compute joint state PDF

ρ (x, t)

Process model Initial conditions Parameters Process noise Prior density

ρ (x(t), t)

+ x t , t

ρ+ (x(t), t)

Measurement model Posterior density Sensor noise

Filtering Problem

Trajectory flow:

dx(t) = f(x, t) dt + g(x, t) dw(t), dw(t) ∼ N (0, Qdt) dz(t) = h(x, t) dt + dv(t), dv(t) ∼ N (0, Rdt)

Density flow:

dρ+= %

LFPdt+

& h(x, t)−Eρ+{h(x, t)} '

⊤R−1&

dz(t)−Eρ+{ h(x, t)}dt '( ρ+

Compute conditional joint state PDF

Filtering Problem

ρ+ := ρ (x, t | z(s), 0 ≤ s ≤ t)

Control Problem

lem

Steer joint state PDF via feedback control over finite time horizon

minimize

u∈U

E !! 1

0 !u!2 2 dt

" subject to dx = f (x, u, t) dt + g (x, t) dw, x (t = 0) ∼ ρ0, x (t = 1) ∼ ρ1

PDFs in Mars Entry-Descent-Landing

Navigational uncertainty Heating uncertainty Chute deployment uncertainty Landing footprint uncertainty

Prediction Problem

Filtering Problem Control Problem

Predict heating rate uncertainty

Supersonic parachute

Navigational uncertainty Heating uncertainty Chute deployment uncertainty Landing footprint uncertainty

Prediction Problem

Filtering Problem Control Problem

Predict heating rate uncertainty Estimate state to deploy parachute

PDFs in Mars Entry-Descent-Landing

PDFs in Mars Entry-Descent-Landing

Gale Crater (4.49S, 137.42E)

Supersonic parachute

Navigational uncertainty Heating uncertainty Chute deployment uncertainty Landing footprint uncertainty

Prediction Problem

Filtering Problem Control Problem

Predict heating rate uncertainty Estimate state to deploy parachute Steer state PDF to achieve desired landing footprint accuracy

Solving prediction problem as Wasserstein gradient flow

What’s New?

Main idea: Solve

Infinite dimensional variational recursion:

Proximal operator: Optimal transport cost: Free energy functional:

Geometric Meaning of Gradient Flow

Gradient Flow in Rn dx dt = rϕ(x), x(0) = x0 Recursion: xk = xk1 hrϕ(xk) = arg min

x2Rn

⇢1 2kx xk1k2

2 + hϕ(x)

• =: proxk·k2

hϕ (xk1)

Convergence: xk ! x(t = kh) as h # 0 Gradient Flow in P2(Rn) ∂ρ ∂t = rW Φ(ρ), ρ(x, 0) = ρ0 Recursion: ρk = ρ(·, t = kh) = arg min

ρ2P2(Rn)

⇢1 2W 2(ρ, ρk1) + hΦ(ρ)

• =: proxW 2

hΦ (ρk1)

Convergence: ρk ! ρ(·, t = kh) as h # 0 as Lyapunov function: as Lyapunov functional:

Geometric Meaning of Gradient Flow

Gradient Flow in Rn dx dt = rϕ(x), x(0) = x0 Recursion: xk = xk1 hrϕ(xk) = arg min

x2Rn

⇢1 2kx xk1k2

2 + hϕ(x)

• =: proxk·k2

hϕ (xk1)

Convergence: xk ! x(t = kh) as h # 0 Gradient Flow in P2(Rn) ∂ρ ∂t = rW Φ(ρ), ρ(x, 0) = ρ0 Recursion: ρk = ρ(·, t = kh) = arg min

ρ2P2(Rn)

⇢1 2W 2(ρ, ρk1) + hΦ(ρ) =: proxW 2

hΦ (ρk1)

Convergence: ρk ! ρ(·, t = kh) as h # 0

D

Algorithm: Gradient Ascent on the Dual Space

Uncertainty propagation via point clouds

No spatial discretization or function approximation

Algorithm: Gradient Ascent on the Dual Space

Algorithm: Gradient Ascent on the Du

@⇢ @t = r · (r ⇢) + 1∆⇢ m Proximal Recursion ⇢k = ⇢(x, t = kh) = arg inf

⇢2P2(Rn)

⇢1 2W 2(⇢, ⇢k1) + h Φ(⇢)

Algorithm: Gradient Ascent on the Dual Space

Algorithm: Gradient Ascent on the Du

@⇢ @t = r · (r ⇢) + 1∆⇢ m Proximal Recursion ⇢k = ⇢(x, t = kh) = arg inf

⇢2P2(Rn)

⇢1 2W 2(⇢, ⇢k1) + h Φ(⇢)

• 2P

⇢

• +

Discrete Primal Formulation %k = arg min

%

⇢ min

M2Π(%k−1,%)

1 2hCk, Mi + h h k1 + 1 log %, %i

Algorithm: Gradient Ascent on the Dual Space

Algorithm: Gradient Ascent on the Du

@⇢ @t = r · (r ⇢) + 1∆⇢ m Proximal Recursion ⇢k = ⇢(x, t = kh) = arg inf

⇢2P2(Rn)

⇢1 2W 2(⇢, ⇢k1) + h Φ(⇢)

• 2P

⇢

• +

Discrete Primal Formulation %k = arg min

%

⇢ min

M2Π(%k−1,%)

1 2hCk, Mi + h h k1 + 1 log %, %i

• ⇢
• +

Entropic Regularization %k = arg min

%

⇢ min

M2Π(%k−1,%)

1 2hCk, Mi + ✏H(M) + h h k1 + 1 log %, %i

• ⇢
• m

Dualization

• pt

, opt

1

= arg max

0,10

⇢ h0, %k1i F ?(1) ✏ h ✓ exp(>

0 h/✏) exp(Ck/2✏) exp(1h/✏)

◆

Recursion on the Cone

Theorem: Consider the recursion on the cone Rn

0 ⇥ Rn

y (Γkz) = %k1, z ⇣ Γk

>y

⌘ = ⇠k1 z ✏

h ,

Then the solution (y ⇤, z⇤) gives the proximal update %k = z⇤ (Γk

>y ⇤)

Algorithmic Setup p

from the

Theorem: Block co-ordinate iteration of (y, z) recur- sion is contractive on Rn

>0 × Rn >0.

Proximal Prediction: 2D Linear Gaussian

Proximal Prediction: 2D Linear Gaussian

Proximal Prediction: Nonlinear Non-Gaussian

Computational Time: Nonlinear Non-Gaussian

1 2 3 4

Physical time tk = kh (seconds)

10−6

Computational time (seconds)

Proximal Prediction: Satellite in Geocentric Orbit

B B B B B B @ dx dy dz dvx dvy dvz 1 C C C C C C A = B B B B B B B B B B B B B B B @ vx vy vz −µx r 3 + (fx)pert − γvx −µy r 3 + (fy)pert − γvy −µz r 3 + (fz)pert − γvz 1 C C C C C C C C C C C C C C C A dt + p 2β−1γ B B B B B B @ dw1 dw2 dw3 1 C C C C C C A , @ fx fy fz 1 A

pert

= @ sθ cφ cθ cφ −sφ sθ sφ cθ sφ cφ cθ −sθ 1 A B B B @ k 2r 4

• 3(sθ)2 − 1
• − k

r 5 sθ cθ 1 C C C A , k := 3J2R2

E, µ = constant

Computational Time: Satellite in Geocentric Orbit

0.005 0.01

Physical time tk = kh (seconds)

10−8 10−7

Computational time (seconds)

Details on Proximal Prediction

Extensions: mean-field models for nonlocal interaction, state-dependent diffusions Publications:

— K.F. Caluya, and A.H., Proximal Recursion for Solving the Fokker-Planck Equation, ACC 2019. — K.F. Caluya, and A.H., Gradient Flow Algorithms for Density Propagation in Stochastic Systems, IEEE

• Trans. Automatic Control 2020, doi: 10.1109/TAC.2019.2951348.

Git repo:

github.com/kcaluya/UncertaintyPropagation

Solving density control as Wasserstein gradient flow

Finite Horizon Feedback Density Control

lem

minimize

u∈U

E !! 1

0 !u (x, t) !2 2 dt

" subject to dx = # f (x, t)+ B(t)u (x, t) $ dt +

√

2"B(t) dw, x (t = 0) ∼ ρ0, x (t = 1) ∼ ρ1

Finite Horizon Feedback Density Control

lem

minimize

u∈U

E !! 1

0 !u (x, t) !2 2 dt

" subject to dx = # f (x, t)+ B(t)u (x, t) $ dt +

√

2"B(t) dw, x (t = 0) ∼ ρ0, x (t = 1) ∼ ρ1

Necessary conditions for optimality: coupled nonlinear PDEs (FPK + HJB)

∂ρopt ∂t

+ ∇ ·

% ρopt% f + B(t)

⊤∇ψ

&&

= %1⊤'

D(t) ⊙ Hess ' ρopt(( 1, ∂ψ ∂t + 1 2)B(t)

⊤∇ψ)2

2 + 〈∇ψ, f〉 = −%〈D(t), Hess (ψ)〉

uopt(x, t) = B(t)

⊤∇ψ

Optimal control: Boundary conditions: ∇

ρopt (x, 0) = ρ0(x), ρopt (x, 1) = ρ1(x)

Feedback Synthesis via the Schrödinger System

Schrödinger’s (until recently) forgotten papers:

Sur la théorie relativiste de l’électron

et l’interprétation de la mécanique quantique

• E. SCHRÖDINGER

• Introduction

J’ai l’intention d’exposer dans

diverses idées concer- nant

la mécanique

quantique

ralement à l’heure actuelle ; je parlerai principalement de la théorie

quantique relativiste du mouvement de l’électron. Autant

que nous pouvons nous en rendre compte aujourd’hui, il semble à peu près

sûr que la mécanique quantique de l’électron, sous

que nous ne possédons pas encore, doit former

jour la base de

toute la physique. A cet intérêt tout à fait général, s’ajoute, ici à Paris, un intérêt particulier : vous savez tous que les bases de la théorie

moderne de

l’électron ont été posées à Paris par votre célèbre compa- triote Louis de BROGLIE.

Schrödinger’s contribution: change of variable Hopf-Cole transform: (ρopt, ψ) %→ (ϕ, ˆ

ϕ) ˆ

Optimal controlled joint state PDF: Optimal control:

%→

ρopt (x, t) = ˆ ϕ(x, t)ϕ(x, t) uopt(x, t) = 2"B(t)

⊤∇ log ϕ (x, t)

ϕ (x, t) = exp %ψ (x, t) 2" & , ˆ ϕ (x, t) = ρopt (x, t) exp %

−ψ (x, t)

2" & ,

Feedback Synthesis via the Schrödinger System

2 coupled nonlinear PDEs boundary-coupled linear PDEs!! Wasserstein proximal algorithm fixed point recursion over (Contractive in Hilbert metric)

( ˆ

ϕ0, ϕ1)

— Y. Chen, T.T. Georgiou, and M. Pavon, Entropic and displacement interpolation: a computational approach using the Hilbert metric, SIAM J. Applied Mathematics, 2016.

∂ ˆ ϕ ∂t = −∇ · ( ˆ ϕ f) + %1⊤(D(t) ⊙ Hess ( ˆ ϕ))1 ) *+ ,

forward Kolmogorov PDE

, ϕ0 ˆ ϕ0 = ρ0, ∂ϕ ∂t = −〈∇ϕ, f〉 − %〈D(t), Hess (ϕ)〉 ) *+ ,

backward Kolmogorov PDE

, ϕ1 ˆ ϕ1 = ρ1.

Feedback Density Control: Zero Prior Dynamics

Feedback Density Control: LTI Prior Dynamics

Feedback Density Control: Nonlinear Prior Dyn.

How to solve the Schrödinger System with nonlinear drift? — No analytical handle on the transition kernel — The backward Kolmogorov PDE cannot be written as Wasserstein gradient flow

Feedback Density Control: Nonlinear Prior Dyn.

How to solve the Schrödinger System with nonlinear drift? — No analytical handle on the transition kernel — The backward Kolmogorov PDE cannot be written as Wasserstein gradient flow Can we exploit some structural nonlinearities in practice?

dx = {−∇V(x) + u(x, t)} dt +

√

2% dw Assume: x ∈ Rn, V ∈ C2 (Rn)

• dξ

dη .

=

• η

−∇ξV(ξ) − κη + u(x, t)

. dt+

√

2%κ

• 0m×m

Im×m . dw Assume: ξ, η ∈ Rm, x := (ξ, η)⊤∈ Rn, n = 2m, V ∈ C2 (Rm), inf V > −∞, Hess (V) unif. bounded

Gradient drift: Mixed conservative

• dissipative drift:

Feedback Density Control: Gradient Drift

Theorem For t ∈ [0, 1], let s := 1 − t. Define the change-of-variables ϕ $→ q $→ p as q(x, s) := ϕ(x, s) = ϕ(x, 1 − t), p(x, s) := q(x, s) exp (−V(x)/%). Then the pair ( ˆ ϕ, p) solves ∂ ˆ ϕ ∂t = ∇ · ( ˆ ϕ∇V) + %∆ ˆ ϕ, ˆ ϕ (x, 0) = ˆ ϕ0(x), ∂p ∂s = ∇ · (p∇V) + %∆p, p (x, 0) = ϕ1(x) exp (−V(x)/%).

Feedback Density Control: Mixed Conservative- Dissipative Drift

Theorem For t ∈ [0, 1], let s := 1 − t. Also, let ϑ := −η. Define the change-of-variables ϕ $→ q $→ / p $→ p as q(ξ, η, s) := ϕ(ξ, η, s) = ϕ(ξ, η, 1 − t), / p(ξ, −η, s) := q(ξ, η, s) exp '− 1

!

'1

2)η)2 2 + V(ξ)

(( , p (ξ, ϑ, s) := / p(ξ, −η, s). Then the pair ( ˆ ϕ, p) solves ∂ ˆ ϕ ∂t = −〈η, ∇ξ ˆ ϕ〉 + ∇η · ' ˆ ϕ '∇ξV (ξ) + κη (( + %κ∆η ˆ ϕ, ∂p ∂s = −〈ϑ, ∇ξp〉 + ∇ϑ · ' p '∇ξV (ξ) + κϑ (( + %κ ∆ϑp, ˆ ϕ (ξ, η, 0) = ˆ ϕ0(ξ, η), p(ξ, ϑ, 0) = ϕ1(ξ, −ϑ) exp '− 1

!

' 1

2)ϑ)2 2 + V(ξ)

(( .

Feedback Density Control via Wasserstein prox.

Details:

— K.F. Caluya, and A.H., Wasserstein Proximal Algorithms for the Schrödinger Bridge Problem: Density Control with Nonlinear Drift, arXiv 1912.01244.

Design proximal recursions over discrete time pair:

(tk−1, sk−1) := ((k − 1)τ, (k − 1)σ), k ∈ N, and

1)τ,

1)σ)are step-sizes. The recursions are of the form:

• ˆ

φtk−1 ϖsk−1 .

$→

• ˆ

φtk ϖsk . := ! " arg inf

ˆ φ∈P2(Rn) 1 2d2 ' ˆ

φtk−1, ˆ φ ( + τF( ˆ φ) arg inf

ϖ∈P2(Rn) 1 2d2 '

ϖsk−1, ϖ ( + σF(ϖ) # 1 $ Consistency guarantees:

ˆ φtk−1(x) → ˆ ϕ(x, t = (k − 1)τ) in L1(Rn) as τ ↓ 0, ϖsk−1(x) → p(x, s = (k − 1)σ) in L1(Rn) as σ ↓ 0.

Feedback Density Control: Gradient Drift

Uncontrolled joint PDF evolution: Optimal controlled joint PDF evolution:

Feedback Density Control: Mixed Conservative- Dissipative Drift

Density Control with Det. Path Constraints

Reflecting Schrödinger Bridge Contraction in the Hilbert metric * Ongoing work

Density Control with Feedback Linearizable Dyn.

Setting:

For x ∈ Rn, u ∈ Rm, and given ρ0, ρ1, consider inf

u∈U

E 2! 1 1 2)u(x, t))2

2 dt

3 , subject to ˙ x = f(x) + G(x)u, x(t = 0) ∼ ρ0(x) x(t = 1) ∼ ρ1(x), with (f(x), G(x)) feedback linearizable, i.e., there exists a triple (δ(x), Γ(x), τ(x)) such that

(∇τ (f(x) + G(x)δ(x)))x=τ−1(z) = Az, (∇τ (G(x)Γ(x)))x=τ−1(z) = B,

where (A, B) is controllable. So, (x, u) $→ (z, v) with ˙ z = Az + Bv, u = δ(x) + Γ(x)v.

Density Control with Feedback Linearizable Dyn.

Main idea: Push-forward the endpoint PDFs via diffeomorphism τ : X !→ Z

σi(z) := τ!ρi = ρi(τ−1(z))

|det(∇xτx=τ−1(z))|,

i ∈ {0, 1}.

Define maps

|

δτ := δ ◦ τ−1, Γτ := Γ ◦ τ−1

Rewrite the problem in feedback linearized coordinates as

minimize

σ,v

! 1 !

Z

1 2L(z, v)σ(z, t) dzdt, subject to ∂σ ∂t + ∇z · ((Az + Bv)σ) = 0 σ(z, t = 0) = σ0, σ(z, t = 1) = σ1, where L(z, v) := &δτ(z) + Γτ(z)v&2

2.

Density Control with Feedback Linearizable Dyn.

— K.F. Caluya, and A.H., Finite Horizon Density Control for Static State Feedback Linearizable Systems, arXiv 1904.02272. — K.F. Caluya, and A.H., Finite Horizon Density Steering for Multi-input State Feedback Linearizable Systems, arXiv 1909.12511.

Details: Optimality: Optimal control: vopt(z, t) = (Γ⊤

τ Γτ(z))−1B⊤∇zψ − Γ−1 τ (z)δτ(z)

HJB:

∂ψ ∂t + 〈∇zψ, Az〉 − 〈∇zψ, BΓτ −1(z)δτ(z)〉 + 1 2 〈∇zψ, B ! Γ⊤

τ (z)Γτ(z)

"−1 B⊤∇zψ〉 = 0.

Solve by dynamic stochastic regularization SBP fixed point recursion

⇝ ⇝

Take Home Message

Emerging system-control theory for densities Wasserstein gradient flow: one unifying framework for the prediction, estimation, and feedback control Feedback density control theory: many recent progress, much remains to be done Several applications: controlling biological and robotic swarm, process control

Thank You

Support: A