Gaussian Mixture Penalization for Trajectory Optimization Problems - - PowerPoint PPT Presentation

Gaussian Mixture Penalization for Trajectory Optimization Problems

• C. Rommel1,2, J. F. Bonnans1,
• B. Gregorutti2 and P. Martinon1

CMAP Ecole Polytechnique - INRIA1 Safety Line2

ISMP - July 2nd 2018

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 1 / 30

Motivation

20 000 airplanes — 80 000 flights per day,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Motivation

20 000 airplanes — 80 000 flights per day, Should double until 2033,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Motivation

20 000 airplanes — 80 000 flights per day, Should double until 2033, Responsible for 3% of CO2 emissions,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Motivation

20 000 airplanes — 80 000 flights per day, Should double until 2033, Responsible for 3% of CO2 emissions, Accounts for 30% of operational cost for an airline,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Motivation

20 000 airplanes — 80 000 flights per day, Should double until 2033, Responsible for 3% of CO2 emissions, Accounts for 30% of operational cost for an airline, Rectilinear climb trajectories at full thrust.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Motivation

20 000 airplanes — 80 000 flights per day, Should double until 2033, Responsible for 3% of CO2 emissions, Accounts for 30% of operational cost for an airline, Rectilinear climb trajectories at full thrust.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 2 / 30

Optimal Control Problem

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP)

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 3 / 30

Optimal Control Problem

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP)

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 3 / 30

Optimal Control Problem

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = ˆ g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP)

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 3 / 30

Dynamics are learned from QAR data

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 4 / 30

Dynamics are learned from QAR data

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 4 / 30

Dynamics are learned from QAR data

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 4 / 30

Dynamics are learned from QAR data

See e.g. [Rommel et al., 2017a] and [Rommel et al., 2017b]

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 4 / 30

Trajectory acceptability

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = ˆ g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) ③ ① ✉ ③

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 5 / 30

Trajectory acceptability

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = ˆ g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) ⇒ ˆ ③ = (ˆ ①, ˆ ✉) solution of (OCP). ③

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 5 / 30

Trajectory acceptability

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = ˆ g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) ⇒ ˆ ③ = (ˆ ①, ˆ ✉) solution of (OCP). Is ˆ ③ inside the validity region of the dynamics model ˆ g ?

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 5 / 30

Trajectory acceptability

min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt, s.t.        ˙ ① = ˆ g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) ⇒ ˆ ③ = (ˆ ①, ˆ ✉) solution of (OCP). Is ˆ ③ inside the validity region of the dynamics model ˆ g ? Does it look like a real aicraft trajectory ?

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 5 / 30

Trajectory acceptability

1NATS UK air traffic control (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 6 / 30

Trajectory acceptability

Pilots acceptance

1NATS UK air traffic control (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 6 / 30

Trajectory acceptability

Pilots acceptance Air Traffic Control1

1NATS UK air traffic control (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 6 / 30

Trajectory acceptability

Pilots acceptance Air Traffic Control1 How can we quantify the closeness from the optimized trajectory to the set of real flights?

1NATS UK air traffic control (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 6 / 30

Likelihood

Let X be a random variable following an absolutely continuous probability distribution with density function f depending on a parameter θ. Then the function L(θ|x) = fθ(x) (1) considered as a function of θ, is the likelihood function of theta, given the

• utcome x of X.

0Picture source: wikipedia, P-Value, author: Repapetilto CC. (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 7 / 30

Likelihood

Let X be a random variable following an absolutely continuous probability distribution with density function f depending on a parameter θ. Then the function L(θ|x) = fθ(x) (1) considered as a function of θ, is the likelihood function of theta, given the

• utcome x of X.

0Picture source: wikipedia, P-Value, author: Repapetilto CC. (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 7 / 30

Likelihood

Let X be a random variable following an absolutely continuous probability distribution with density function f depending on a parameter θ. Then the function L(θ|x) = fθ(x) (1) considered as a function of θ, is the likelihood function of theta, given the

• utcome x of X.

In our case: the optimized trajectory plays the role of θ,

0Picture source: wikipedia, P-Value, author: Repapetilto CC. (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 7 / 30

Likelihood

Let X be a random variable following an absolutely continuous probability distribution with density function f depending on a parameter θ. Then the function L(θ|x) = fθ(x) (1) considered as a function of θ, is the likelihood function of theta, given the

• utcome x of X.

In our case: the optimized trajectory plays the role of θ, the set of real flights plays the role of x,

0Picture source: wikipedia, P-Value, author: Repapetilto CC. (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 7 / 30

How to apply this to functional data?

Assumption: We suppose that the real flights are observations of the same functional random variable Z = (Zt) valued in C(T, E), with E compact subset of Rd and T = [0, tf ].

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 8 / 30

How to apply this to functional data?

Assumption: We suppose that the real flights are observations of the same functional random variable Z = (Zt) valued in C(T, E), with E compact subset of Rd and T = [0, tf ]. Problem: Computation of probability densities in infinite dimensional space is untractable...

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 8 / 30

How to apply this to functional data?

Assumption: We suppose that the real flights are observations of the same functional random variable Z = (Zt) valued in C(T, E), with E compact subset of Rd and T = [0, tf ]. Problem: Computation of probability densities in infinite dimensional space is untractable... Standard approach FDA: use FPCA to decompose the data in a small number of coefficients ⇒

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 8 / 30

How to apply this to functional data?

Assumption: We suppose that the real flights are observations of the same functional random variable Z = (Zt) valued in C(T, E), with E compact subset of Rd and T = [0, tf ]. Problem: Computation of probability densities in infinite dimensional space is untractable... Standard approach FDA: use FPCA to decompose the data in a small number of coefficients Or: we can aggregate the marginal densities

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 8 / 30

Why does it make sense for this type of data?

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 9 / 30

Why does it make sense for this type of data?

Likely values of flight variables during climb are strongly dependent on the altitude

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 9 / 30

Why does it make sense for this type of data?

Likely values of flight variables during climb are strongly dependent on the altitude

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 9 / 30

How do we aggregate the marginal likelihoods?

② ② ② ②

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② ② ② ②

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② new trajectory, ② ② ②

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② new trajectory, ft(②(t)) marginal likelihood of ② at t, i.e. likelihood of observing Zt = ②(t).

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② new trajectory, ft(②(t)) marginal likelihood of ② at t, i.e. likelihood of observing Zt = ②(t). Why not average over time ?... 1 tf tf ft(②(t))dt

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② new trajectory, ft(②(t)) marginal likelihood of ② at t, i.e. likelihood of observing Zt = ②(t). Why not average over time ?... 1 tf tf ft(②(t))dt Marginal densities may have really different shapes

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

ft marginal density of Z, i.e. probability density function of Zt, ② new trajectory, ft(②(t)) marginal likelihood of ② at t, i.e. likelihood of observing Zt = ②(t).

Mean marginal likelihood [Rommel et al., 2018]

MML(Z, ②) = 1 tf tf ψ[ft, ②(t)]dt, where ψ : L1(E, R+) × R → [0; 1] is a continuous scaling map.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 10 / 30

How do we aggregate the marginal likelihoods?

Possible scalings are the normalized density ψ[ft, ②(t)] := ②(t) max

z∈E ft(z),

② ②

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 11 / 30

How do we aggregate the marginal likelihoods?

Possible scalings are the normalized density ψ[ft, ②(t)] := ②(t) max

z∈E ft(z),

• r the confidence level

ψ[ft, ②(t)] := P (ft(Zt) ≤ ft(②(t))) .

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 11 / 30

How do we deal with sampled curves?

③ ②

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 12 / 30

How do we deal with sampled curves?

In practice, the m trajectories are sampled at variable discrete times: T D := {(tr

j , zr j )} 1≤j≤n 1≤r≤m

⊂ T × E, zr

j := ③r(tr j ),

Y := {(˜ tj, yj)}˜

n j=1 ⊂ T × E,

yj := ②(˜ tj).

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 12 / 30

How do we deal with sampled curves?

In practice, the m trajectories are sampled at variable discrete times: T D := {(tr

j , zr j )} 1≤j≤n 1≤r≤m

⊂ T × E, zr

j := ③r(tr j ),

Y := {(˜ tj, yj)}˜

n j=1 ⊂ T × E,

yj := ②(˜ tj). Hence, we approximate the MML using a Riemann sum which aggregates consistent estimators ˆ f m

˜ tj of the marginal densities f˜ tj:

EMMLm(T D, Y) := 1 tf

˜ n

• j=1

ψ[ˆ f m

˜ tj , yj]∆˜

tj.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 12 / 30

How can we estimate marginal densities?

Suppose that sampling times {tr

j : j = 1, . . . , n; r = 1, . . . , m} are

i.i.d. sampled from r.v. T, indep. Z;

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 13 / 30

How can we estimate marginal densities?

Suppose that sampling times {tr

j : j = 1, . . . , n; r = 1, . . . , m} are

i.i.d. sampled from r.v. T, indep. Z; ft is the density of Zt = (ZT|T = t) = (Y |X);

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 13 / 30

How can we estimate marginal densities?

Suppose that sampling times {tr

j : j = 1, . . . , n; r = 1, . . . , m} are

i.i.d. sampled from r.v. T, indep. Z; ft is the density of Zt = (ZT|T = t) = (Y |X); Our problem can be seen as a conditional probability density learning problem with (X, Y ) = (T, ZT).

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 13 / 30

How can we estimate marginal densities?

Suppose that sampling times {tr

j : j = 1, . . . , n; r = 1, . . . , m} are

i.i.d. sampled from r.v. T, indep. Z; ft is the density of Zt = (ZT|T = t) = (Y |X); Our problem can be seen as a conditional probability density learning problem with (X, Y ) = (T, ZT). ⇒ We could apply SOA conditional density estimation techniques, such as LS-CDE [Sugiyama et al., 2010],

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 13 / 30

How can we estimate marginal densities?

Suppose that sampling times {tr

j : j = 1, . . . , n; r = 1, . . . , m} are

i.i.d. sampled from r.v. T, indep. Z; ft is the density of Zt = (ZT|T = t) = (Y |X); Our problem can be seen as a conditional probability density learning problem with (X, Y ) = (T, ZT). ⇒ We could apply SOA conditional density estimation techniques, such as LS-CDE [Sugiyama et al., 2010], ⇒ Instead, we choose to use a fine partitioning of the time domain.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 13 / 30

Partition based marginal density estimation

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 14 / 30

Partition based marginal density estimation

Idea: to average in time the marginal densities over small bins by applying classical multivariate density estimation techniques to each subset.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 14 / 30

Consistency

Assumtion 1 - Positive time density

ν ∈ L∞(E, R+) density function of T, s.t. ν+ := ess sup

t∈T

ν(t) < ∞, ν− := ess inf

t∈T ν(t) > 0.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 15 / 30

Consistency

Assumtion 1 - Positive time density

ν ∈ L∞(E, R+) density function of T, s.t. ν+ := ess sup

t∈T

ν(t) < ∞, ν− := ess inf

t∈T ν(t) > 0.

Assumtion 2 - Lipschitz in time

Function (t, z) ∈ T × E → ft(z) is continuous and |ft1(z) − ft2(z)| ≤ L|t1 − t2|, L > 0.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 15 / 30

Consistency

Assumtion 1 - Positive time density

ν ∈ L∞(E, R+) density function of T, s.t. ν+ := ess sup

t∈T

ν(t) < ∞, ν− := ess inf

t∈T ν(t) > 0.

Assumtion 2 - Lipschitz in time

Function (t, z) ∈ T × E → ft(z) is continuous and |ft1(z) − ft2(z)| ≤ L|t1 − t2|, L > 0.

Assumption 3 - Shrinking bins

The homogeneous partition {Bm

ℓ }qm ℓ=1 of [0; tf ], with binsize bm, is s.t.

lim

m→∞ bm = 0,

lim

m→∞ mbm = ∞.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 15 / 30

Consistency

Assumption 4 - i.i.d. consistency

S = {(zk)N

k=1 ∈ E N : N ∈ N∗} set of finite sequences,

Θ : S → L1(E, R+) multivariate density estimation statistic, G arbitrary family of probability density functions on E, ρ ∈ G, SN

ρ i.i.d sample of size N drawn from ρ valued in S.

The estimator obtained by applying Θ to SN

ρ , denoted by

ˆ ρN := Θ[SN

ρ ] ∈ L1(E, R+),

is a (pointwise) consistent density estimator, uniformly in ρ: For all z ∈ E, ε > 0, α1 > 0, there is Nε,α1 > 0 such that, for any ρ ∈ G, N ≥ Nε,α1 ⇒ P

• ˆ

ρN(z) − ρ(z)

• < ε
• > 1 − α1.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 16 / 30

Consistency

We denote by: ℓm(t) :=

• t

bm

• maps time to index of bin containing it;

ˆ f m

ℓm(t) := Θ[T m ℓm(t)] estimator trained using subset of data points T m ℓm(t)

whose sampling times fall in the bin containing t;

Theorem 1 - [Rommel et al., 2018]

Under assumptions 1 to 4, for any z ∈ E and t ∈ T, ˆ f m

ℓm(t)(z) consistently

approximates the marginal density ft(z) as the number of curves m grows: ∀ε > 0, lim

m→∞ P

• |ˆ

f m

ℓm(t)(z) − ft(z)| < ε

• = 1.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 17 / 30

Marginal density estimation results

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 18 / 30

Marginal density estimation results

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 19 / 30

Marginal density estimation results

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 20 / 30

How good is it compared to other methods?

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

How good is it compared to other methods?

Training set of m = 424 flights ≃ 334 531 point observations,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

How good is it compared to other methods?

Training set of m = 424 flights ≃ 334 531 point observations, Test set of 150 flights

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

How good is it compared to other methods?

Training set of m = 424 flights ≃ 334 531 point observations, Test set of 150 flights = 50 real flights (Real), 50 optimized flights with operational constraints (Opt1) and 50 optimized flights without constraints (Opt2);

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

How good is it compared to other methods?

Training set of m = 424 flights ≃ 334 531 point observations, Test set of 150 flights = 50 real flights (Real), 50 optimized flights with operational constraints (Opt1) and 50 optimized flights without constraints (Opt2); Discrimination power comparison with (gmm-)FPCA and (integrated) LS-CDE:

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

How good is it compared to other methods?

Training set of m = 424 flights ≃ 334 531 point observations, Test set of 150 flights = 50 real flights (Real), 50 optimized flights with operational constraints (Opt1) and 50 optimized flights without constraints (Opt2); Discrimination power comparison with (gmm-)FPCA and (integrated) LS-CDE:

Var. Estimated Likelihoods Real Opt1 Opt2 MML 0.63 ± 0.07 0.43 ± 0.08 0.13 ± 0.02 FPCA 0.16 ± 0.12 6.4e-03 ± 3.8e-03 3.6e-03 ± 5.4e-03 LS-CDE 0.77 ± 0.05 0.68 ± 0.04 0.49 ± 0.06

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 21 / 30

MML penalty

The MML can be used not only to assess the optimization solutions, but also to penalize the optimization itself: min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt − λ MML(Z, ①), s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP)

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 22 / 30

MML penalty

The MML can be used not only to assess the optimization solutions, but also to penalize the optimization itself: min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt − λ MML(Z, ①), s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP)

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 22 / 30

MML penalty

The MML can be used not only to assess the optimization solutions, but also to penalize the optimization itself: min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt − λ MML(Z, ①), s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) λ sets trade-off between a fuel minimization and a likelihood maximization,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 22 / 30

MML penalty

The MML can be used not only to assess the optimization solutions, but also to penalize the optimization itself: min

(①,✉)∈X×U

tf C(t, ✉(t), ①(t))dt − λ MML(Z, ①), s.t.        ˙ ① = g(t, ✉, ①), for a.e. t ∈ [0, tf ], Φ(①(0), ①(tf )) ∈ KΦ, ✉(t) ∈ Uad, ①(t) ∈ Xad, for a.e. t ∈ [0, tf ], c(✉(t), ①(t)) ≤ 0, for all t ∈ [0, tf ]. (OCP) λ sets trade-off between a fuel minimization and a likelihood maximization, If (OCP) is solved using NLP techniques, parametric estimator of MML is needed.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 22 / 30

Gaussian mixture model for marginal densities

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 23 / 30

Gaussian mixture model for marginal densities

ft(z) =

K

• k=1

wt,kφ(z, µt,k, Σt,k),

K

• k=1

wt,k = 1, wt,k ≥ 0, φ(z, µ, Σ) := 1

• (2π)d det Σ

e− 1

2 (z−µ)⊤Σ−1(z−µ). (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 23 / 30

Gaussian mixture model for marginal densities

ft(z) =

K

• k=1

wt,kφ(z, µt,k, Σt,k),

K

• k=1

wt,k = 1, wt,k ≥ 0, φ(z, µ, Σ) := 1

• (2π)d det Σ

e− 1

2 (z−µ)⊤Σ−1(z−µ).

Assuming that the number of components is known, the weights wt,k, means µt,k and covariance matrices Σt,k need to be estimated.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 23 / 30

Maximum likelihood parameters estimation

For K = 1, maximum likelihood estimates have closed form:

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 24 / 30

Maximum likelihood parameters estimation

For K = 1, maximum likelihood estimates have closed form: L(µt,1, Σt,1|z1, . . . , zN) =

N

• i=1

1

• (2π)d det Σt,1

e− 1

2 (z−µt,1)⊤Σ−1 t,1 (z−µt,1) (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 24 / 30

Maximum likelihood parameters estimation

For K = 1, maximum likelihood estimates have closed form: L(µt,1, Σt,1|z1, . . . , zN) =

N

• i=1

1

• (2π)d det Σt,1

e− 1

2 (z−µt,1)⊤Σ−1 t,1 (z−µt,1)

ˆ θ := (ˆ µt,1, ˆ Σt,1) = arg min

(µt,1,Σt,1) N

• i=1
• log det Σt,1 + (zi − µt,1)⊤Σ−1

t,1(zi − µt,1)

• (CMAP, INRIA, Safety Line)

GMM Penalization for OCP ISMP - July 2nd 2018 24 / 30

Maximum likelihood parameters estimation

For K = 1, maximum likelihood estimates have closed form: L(µt,1, Σt,1|z1, . . . , zN) =

N

• i=1

1

• (2π)d det Σt,1

e− 1

2 (z−µt,1)⊤Σ−1 t,1 (z−µt,1)

ˆ θ := (ˆ µt,1, ˆ Σt,1) = arg min

(µt,1,Σt,1) N

• i=1
• log det Σt,1 + (zi − µt,1)⊤Σ−1

t,1(zi − µt,1)

• ˆ

µt,1 = 1 N

N

• i=1

zi, ˆ Σt,1 = 1 N

N

• i=1

(zi − ˆ µt,1)(zi − ˆ µt,1)⊤.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 24 / 30

EM Algorithm

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K},

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K}, If ith observation Ji = k, then zi was drawn from the kth component,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K}, If ith observation Ji = k, then zi was drawn from the kth component, Group observations by component and compute (ˆ µt,k, ˆ Σt,k) with K = 1 maximum likelihood formulas.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K}, If ith observation Ji = k, then zi was drawn from the kth component, Group observations by component and compute (ˆ µt,k, ˆ Σt,k) with K = 1 maximum likelihood formulas.

Expectation-Maximization - [Dempster et al., 1977]

Initialization: ˆ θ = ( ˆ wt,k, ˆ µt,k, ˆ Σt,k)K

k=1 = (w0 t,k, µ0 t,k, Σ0 t,k)K k=1,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K}, If ith observation Ji = k, then zi was drawn from the kth component, Group observations by component and compute (ˆ µt,k, ˆ Σt,k) with K = 1 maximum likelihood formulas.

Expectation-Maximization - [Dempster et al., 1977]

Initialization: ˆ θ = ( ˆ wt,k, ˆ µt,k, ˆ Σt,k)K

k=1 = (w0 t,k, µ0 t,k, Σ0 t,k)K k=1,

Expectation: For k = 1, . . . , K and i = 1, . . . , N, ˆ wt,k = 1 N

N

• i=1

ˆ πk,i, ˆ πk,i := P(Ji = k|ˆ θt, Zh) = ˆ µt,kφ(zi, ˆ µt,k, ˆ Σt,k) N

j=1 ˆ

wt,kφ(zj, ˆ µt,k, ˆ Σt,k) .

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

EM Algorithm

Hidden random variable J valued on {1, . . . , K}, If ith observation Ji = k, then zi was drawn from the kth component, Group observations by component and compute (ˆ µt,k, ˆ Σt,k) with K = 1 maximum likelihood formulas.

Expectation-Maximization - [Dempster et al., 1977]

Initialization: ˆ θ = ( ˆ wt,k, ˆ µt,k, ˆ Σt,k)K

k=1 = (w0 t,k, µ0 t,k, Σ0 t,k)K k=1,

Expectation: For k = 1, . . . , K and i = 1, . . . , N, ˆ wt,k = 1 N

N

• i=1

ˆ πk,i, ˆ πk,i := P(Ji = k|ˆ θt, Zh) = ˆ µt,kφ(zi, ˆ µt,k, ˆ Σt,k) N

j=1 ˆ

wt,kφ(zj, ˆ µt,k, ˆ Σt,k) . Maximization: ˆ µt,k = N

i=1 ˆ

πk,izi N

i=1 ˆ

πk,i , ˆ Σt,k = N

i=1 ˆ

πk,i(zi − ˆ µt,k)(zi − ˆ µt,k)⊤ N

i=1 ˆ

πk,i .

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 25 / 30

Penalty effect

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 26 / 30

Consumption x Acceptability trade-off

Figure: Average over 20 flights of the fuel consumption and MML score (called acceptability here) of optimized trajectories with varying MML-penalty weight λ.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 27 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

3 Applicable to the case of aircraft climb trajectories, (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

3 Applicable to the case of aircraft climb trajectories,

Competitive with other well-established SOA approaches,

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

3 Applicable to the case of aircraft climb trajectories,

Competitive with other well-established SOA approaches,

4 Particular Gaussian mixture model implementation, (CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

3 Applicable to the case of aircraft climb trajectories,

Competitive with other well-established SOA approaches,

4 Particular Gaussian mixture model implementation,

Showed that it can be used in optimal control problems to obtain solutions close to optimal, and still realistic.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

Conclusion

1 General probabilistic criterion for quantifying the closeness between a

curve and a set random trajectories,

2 Class of consistent plug-in estimators, based on “histogram” of

multivariate density estimators,

3 Applicable to the case of aircraft climb trajectories,

Competitive with other well-established SOA approaches,

4 Particular Gaussian mixture model implementation,

Showed that it can be used in optimal control problems to obtain solutions close to optimal, and still realistic.

⇒ How could we automatically set the trade-off ?...

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 28 / 30

THANK YOU FOR YOUR ATTENTION !!

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 29 / 30

References

Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society. Series B (methodological), pages 1–38. Rommel, C., Bonnans, J. F., Gregorutti, B., and Martinon, P. (2017a). Aircraft dynamics identification for optimal control. In Proceedings of the 7th European Conference for Aeronautics and Aerospace Sciences. Rommel, C., Bonnans, J. F., Gregorutti, B., and Martinon, P. (2017b). Block sparse linear models for learning structured dynamical systems in

• aeronautics. HAL report hal-01816400.

Rommel, C., Bonnans, J. F., Gregorutti, B., and Martinon, P. (2018). Quantifying the closeness to a set of random curves via the mean marginal likelihood. HAL report hal-01816407. Sugiyama, M., Takeuchi, I., Suzuki, T., Kanamori, T., Hachiya, H., and Okanohara, D. (2010). Conditional density estimation via least-squares density ratio estimation. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, pages 781–788.

(CMAP, INRIA, Safety Line) GMM Penalization for OCP ISMP - July 2nd 2018 30 / 30