98 IEEE CONTROL SYSTEMS LETTERS, VOL. 1, NO. 1, JULY 2017

A Scalable Sampling Method to

High-Dimensional Uncertainties for

Optimal and Reinforcement

Learning-Based Controls

Junfei Xie, Yan Wan, Kevin Mills, James J. Filliben, and F. L. Lewis

Abstract—Modern dynamical systems often operate in

environments of high-dimensional uncertainties that modulate system dynamics in a complicated fashion. These highdimensional uncertainties, non-Gaussian in many realistic

scenarios, complicate real-time system analysis, design,

and control tasks. In this letter, we address the scalability

of computation for systems of high-dimensional uncertainties by introducing new sampling methods, the multivariate

probabilistic collocation method (M-PCM), and its extension called M-PCM-orthogonal fractional factorial design

(OFFD) which integrates M-PCM with the OFFDs to break

the curse of dimensionality. We explore the capabilities of

M-PCM and M-PCM-OFFD-based optimal control and adaptive control using the reinforcement learning approach. The

analyses and simulation studies illustrate the efficiency

and effectiveness of these two approaches.

Index Terms—Uncertain systems, optimal control, adaptive control.

I. INTRODUCTION

MODERN dynamical systems, such as complex information systems, power networks, and air traffic systems,

often operate in environments of high-dimensional uncertainties. These high-dimensional uncertainties typically have the

following features: 1) they modulate system dynamics in a

complicated fashion that cannot be captured by simple additive

white noises, and hence significantly complicate real-time system analysis, design and control tasks, and 2) their statistical

information can be obtained from environmental forecasting

Manuscript received March 6, 2017; revised May 16, 2017; accepted

May 19, 2017. Date of publication May 26, 2017; date of current version

June 8, 2017. This work was supported in part by the National Institute of

Standards and Technology under Grant 60NANB13D172, and in part by

the National Science Foundation under Grant EAGER-1522458, Grant

CAREER-1714519, and Grant CPS-1714826. Recommended by Senior

Editor F. Blanchini. (Corresponding author: Yan Wan.)

J. Xie is with the Department of Computer Science, Texas A&M

University at Corpus Christi, Corpus Christi, TX 78412 USA (e-mail:

junfei.xie@tamucc.edu).

Y. Wan and F. L. Lewis are with the Department of Electrical

Engineering, University of Texas at Arlington, Arlington, TX 76019 USA

(e-mail: yan.wan@uta.edu; lewis@uta.edu).

K. Mills and J. J. Filliben are with NIST, Gaithersburg, MD 20899 USA

(e-mail: kmills@nist.gov).

Digital Object Identifier 10.1109/LCSYS.2017.2708598

tools that are independent to system internal dynamics. For

instance, strategic air traffic flow management (ATFM) plans

traffic in the strategic time frame (with 2 − 15 hours lookahead time), during which a wide range of weather scenarios

can possibly occur [1]. The uncertain weather information is

available from probabilistic forecasting tools such as the Very

Short Range Ensemble Forecast System (VSREF). The weather

uncertainties modulate region capacities in the National

Airspace System, posing uncertain nonlinear constraints to

interweaving traffic flows [2]. Techniques are needed to

quickly design optimal and adaptive control strategies that are

robust and scalable to these high-dimensional uncertainties.

Stochastic optimal control has been widely studied in the

literature (e.g., [3]–[7]). For systems of linear dynamics, additive noise and quadratic cost, linear quadratic control can

easily find a solution [4]. However, for modern dynamical

system applications, the dynamics often may only be captured by complicated simulators rather than clean mathematics,

and uncertainties modulate system dynamics in a nonlinear

fashion. Backward-in-time dynamic programming approaches

(e.g., based on the Bellman optimality equation) and forwardin-time Monte Carlo (MC) reinforcement learning-based adaptive control [4], [6] have been used. These methods typically

involve the discretization and sampling of the uncertainty

space to estimate expected cost. As evaluating system dynamics over all discretized values of uncertain parameters is computationally expensive for large-scale dynamical systems [3],

appropriate discretization and sampling schemes are pivotal

to the accuracy and scalability of stochastic optimal control

strategies [8].

The MC simulation method and its variants such as Markov

Chain MC and Sequential MC have been used to sample the uncertainty space. However, they require a large

number of simulations to converge to the mean cost estimates, making it computationally impractical for the real-time

control of large-scale systems that typically consume considerable computation for a single simulation run. To reduce the

computational cost, quadrature schemes such as Euler [9],

Composite-Simpson [10], LGL-quadrature [11], and sparsegrid (SG) [12], have recently been applied to sample the

uncertainty space. However, their accuracy and computational

scalability are either unsatisfactory or unjustified for arbitrary

2475-1456 c 2017 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.

See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.
XIE et al.: SCALABLE SAMPLING METHOD TO HIGH-DIMENSIONAL UNCERTAINTIES 99

forms of uncertainties. Another approach is to approximate

the system dynamics modulated by uncertainties as a master

Markov chain of the state space being the Cartesian product of

the system state space and the uncertainty space, and solve the

problem using Markov decision process approaches [13]. This

approach is again not scalable. It works well for systems of

simple dynamics, but is ineffective for systems of large state

and uncertainty spaces.

To address these challenges, we develop in this letter

an optimal control framework that is scalable to highdimensional uncertainties. The framework builds on two

multi-dimensional uncertainty evaluation approaches that we

developed recently, the multivariate probabilistic collocation

method (M-PCM) [14], and the scalable M-PCM-OFFD that

integrates M-PCM with orthogonal fractional factorial designs

(OFFDs) [15], [16].

In this letter, we explore the capabilities of M-PCM and

M-PCM-OFFD in facilitating optimal control for systems

of high-dimensional uncertainties, which are introduced in

Section II. Three representative scenarios are considered:

1) finite-horizon optimal control with uncertain parameters

changing independently across time (Section III-A), 2) finitehorizon optimal control with uncertain parameters evolving

according to Markov chains (Section III-B), and 3) infinitehorizon forward-in-time optimal control using reinforcement

learning (Section IV). For arbitrary system dynamics captured

by a black-box simulator, we prove that the M-PCM and MPCM-OFFD based stochastic optimal controls find accurate

solutions with very limited computational costs, under simple assumptions that hold for broad realistic systems. These

features are further validated through illustrative examples

and comparative simulation studies with existing methods in

Section V.

II. UNCERTAINTY SAMPLING METHODS:

M-PCM AND M-PCM-OFFD

M-PCM [14] and M-PCM-OFFD [15], [16] are newly developed computationally effective sampling methods to evaluate

output statistics for system mappings of high-dimensional

uncertainties.

A. M-PCM

The M-PCM smartly selects a small set of samples according to the statistics (e.g., pdfs) of uncertain parameters, and

runs simulations at these samples to produce a reduced-order

mapping that maintains precisely the output statistics of the

original mapping (see [14, Sec. II.B] for detailed design procedures). Lemma 1 shows the key result on the performance

of M-PCM in accurately estimating mapping statistics under

uncertainty.

Lemma 1 [14, Th. 2]: Consider a system mapping modulated by m independent uncertain parameters:

G(a1,..., am) =

2

n1−1

j1=0

2

n2−1

j2=0

···

2n

m−1

jm=0

j1,...,jm

m

i=1

a

ji

i (1)

where ai is an uncertain parameter with the degree up to

2ni − 1. ni are positive integers for all i ∈ {1, 2,..., m}, and

j1,...,jm ∈ R are the coefficients. Each uncertain parameter ai

follows an independent pdf fAi(ai). The M-PCM approximates

G(a1,..., am) with the following low-order mapping

G

(a1,..., am) =

n

1−1

j1=0

n

2−1

j2=0

···

n

m−1

jm=0

j1,...,jm

m

i=1

a

ji

i , (2)

with E[G(a1,..., am)] = E[G

(a1,..., am)] , where j1,...,jm ∈ R

are coefficients. The M-PCM reduces the number of simulations from 2m m

i=1 ni to mi=1 ni.

Remarks: The knowledge of uncertainty is typically available to realistic system studies in the form of probabilistic

forecasts or historical data. We note that M-PCM is not limited

by the knowledge of precise pdfs. When low-order moments

or sample data of uncertain parameters are available, samplemoment based approaches can be used to select M-PCM

simulation points (see [14, Sec. V.B]). We can also fit the sample data or low-order moments (e.g., mean and variance) with

canonical pdfs. Beyond accurate output mean estimation, MPCM also has other nice statistical properties such as accurate

cross-statistics estimation and tight connection to minimum

mean squares estimation (see [14, Sec. III.A]).

B. M-PCM-OFFD

M-PCM significantly reduces the number of simulations

to estimate output mean, however its computation cost does

not scale with the number of uncertain parameters. M-PCMOFFD [15], [16] further breaks the curse of dimensionality through leveraging the systematic procedures and nice

properties of 2-level OFFD such as balance and orthogonality [17], [18]. We show that M-PCM-OFFD has better

performance in terms of accuracy and computational scalability for estimating output statistics, compared to existing

uncertainty sampling approaches, such as stochastic response

surface method and SG [12], [19]. Lemma 2 illustrates the key

capability of M-PCM-OFFD. Please refer to [16, Sec. 5.2] for

detailed design procedures.

Lemma 2 [16, Sec. 5.2]: Consider an m-parameter system

mapping (Equation (1)) with each input parameter ai of degree

up to 3 (i.e., ni = 2, ∀i ∈ {1, 2,..., m}). Assume that its coefficients j1,...,jm = 0 if more than τ of j1,..., jm are non-zero,

where 1 ≤ τ ≤ m. The integrated M-PCM and 2m−γ ∗

OFFD

approximates the original system mapping with the following

low-order mapping.

G∗(a1,..., am) = 

1

j1=0



1

j2=0

··· 

1

jm=0

j1,...,jm

m

i=1

a

ji

i , (3)

such that E[G(a1,..., am)] = E[G∗(a1,..., am)], where coefficients j1,...,jm = 0 if more than τ of j1,...,jm

are non-zero, γ ∗ = max{γ | 1 ≤ γ ≤ m −

log2(

τ

i=0

 i

m



)	, and 2

m−γ

R OFFD exists, with R ≥ 2τ + 1}.

The M-PCM-OFFD reduces the number of simulations from

22m to 2m−γ ∗

in the range of 2log2(m+1)	, 2m−1

	

, making it

scalable with the number of uncertain parameters.

Remarks: The assumption on j1,...,jm reflects that the interacting effects that involve a large number of parameters are

less important to the output statistics than those that involve a

few parameters. This assumption typically holds for realistic

modern complex systems, and has been widely used in the
100 IEEE CONTROL SYSTEMS LETTERS, VOL. 1, NO. 1, JULY 2017

field of statistical experimental designs [17], [18]. Other distinguishing features of M-PCM-OFFD that have been proved

include the robustness to computational approximations in system simulators (see [16, Sec. 5.3]). We also note that in the

lemma, the term of the highest degree is of degree 3m. The

lemma can be generalized to higher degrees through replacing

the 2-level OFFD with higher-level OFFD, which we leave for

the future work.

III. OPTIMAL CONTROL FOR SYSTEMS OF HIGH

DIMENSIONAL UNCERTAINTIES

In this section, we study two optimal control problems for

systems of high dimensional uncertainties.

A. Optimal Control for Systems of Uncertain Parameters

Independent Across Time

Consider a generic dynamical system modulated by an

m-dimensional time-varying uncertain vector a[k]. Each element of a[k], ai[k], changes independently over time with pdf

fAi[k](ai[k]). The system dynamic is

x[k + 1] = hk(x, u, a), (4)

with state vector x[k] ∈ S and control input vector u[k] ∈ C,

where S and C are the known state space and control space,

respectively. hk() is the system dynamic function of x[k], u[k]

and a[k]. The total expected cost equals [3]

JN(x[0]) = Ea[0]{· · · Ea[N−1]{

N

−1

k=0

αkgk(x, u) + αNqN(x)}···}

(5)

where Ea[k]() is the expectation of function within () with

respect to the uncertain vector a[k]. α ∈ (0, 1] is a discount factor. gk() and qk() are the running and terminal cost

at time step k, respectively. Consider the problem of finding an optimal control policy π∗ = {μ∗

0, μ∗1,...,μ∗N−1},

with u∗[k] = μ∗

k (x[k]), such that the total expected cost is

minimized, i.e.,

π∗ = arg min

π {JN(x[0])}. (6)

μk is a control function that maps S into C.

This finite-horizon control problem can be solved

using backward-in-time methods, e.g., dynamic

programming. As the uncertain parameters are independent from the states, we define the value function

as Vk(x[k]) = Ea[k]{· · · Ea[N−1]{

N−1

i=k αi−kgi(x, u) +

αN−kqN(x)}···} = Ea[k]{· · · Ea[N−1]{gk(x, u) +

α[

N−1

i=k+1 αi−(k+1)gi(x, u) + αN−(k+1)qN(x)]}···} =

Ea[k]{gk(x, u) + αEa[k+1]{· · · Ea[N−1]{

N−1

i=k+1 αi−(k+1)gi(x, u) +

αN−(k+1)qN(x)} · · · }}, and VN(x[N]) = qN(x). The following

Bellman’s Equation holds.

Vk(x[k]) = Ea[k]



gk(x, u) + αVk+1(x[k + 1])

	

, (7)

The optimal cost J∗

N(x[0]) = V∗

0 (x[0]) can then be derived

by working backward in time using dynamic programming

according to the Bellman optimality equation [20]:

V∗

k (x[k]) = min

π Ea[k]




gk(x, u) + αV∗

k+1(x[k + 1])



(8)

Defining Gk(x, u, a) = gk(x, u) + αV∗

k+1(x[k + 1]), we then

notice that given an admissible state x[k] and control value

u[k], Ea[k][Gk(x, u, a)] can be approximated by the mean

output of a system mapping Gk(x, u, a) using M-PCM or

M-PCM-OFFD. When each uncertain parameter ai[k] has a

degree up to 2ni[k] − 1, Gk(x, u, a) has the following form

Gk(x, u, a) =

2n1



[k]−1

j1=0

...

2nm



[k]−1

jm=0

j1,...,jm (x[k], u[k])

m

i=1

a

ji

i [k],

(9)

Theorem 1 holds, where j1,...,jm (x[k], u[k]) ∈ R are the

coefficients determined by state x[k] and control input u[k].

Theorem 1: Consider a dynamical system described by

Equation (4), with cost and value functions given by

Equation (5) and Equations (7)–(9), respectively. By applying

dynamic programming, and sampling the uncertainty space at

each iteration using M-PCM, the true optimal control policy

can be obtained with no error.

Proof: First, we introduce a set of notions for the

optimal control obtained using M-PCM. According to

Lemma 1, for a given admissible state x[k] and control input u[k], Gk(x, u, a) can be approximated using MPCM as a low-order function of the form G

k(x, u, a) = n1[k]−1

j1=0 ···nm[k]−1

jm=0 j1,...,jm (x[k], u[k])

m

i=1 a

ji

i [k] , such that

Ea[k]



Gk(x, u, a)

	 = Ea[k]

G

k(x, u, a)

	

. Denote V∗

k (x[k]) =

min

u[k]

Ea[k][G

k(x, u, a)] and π∗ = {μ∗0 , μ∗1 ,...,μ∗N−1}. The optimal control policy can be found by evaluating V∗

k (x[k]) at

only the M-PCM samples at each iteration. Denote u∗[k] =

μ∗

k (x[k]).

In order to prove this theorem, we need to show u∗[k] =

u∗[k], i.e., for each x[k] ∈ S, arg min

u[k]

Ea[k][Gk(x, u, a)] =

arg min

u[k]

Ea[k]



G

k(x, u, a)

	

. This is equivalent to the following two statements: 1) u∗[k] = u∗[k], such that

Ea[k]



G

k(x, u∗, a)

	 < Ea[k]



Gk(x, u∗, a)

	

, and 2) u∗[k] =

u∗[k], such that Ea[k]



Gk(x, u∗, a)

	

< Ea[k]



G

k(x, u∗, a)

	

. To

prove statement 1 using contradiction, we assume that such a

u∗[k] = u∗[k] exists. Lemma 1 leads to Ea[k]



G

k(x, u∗, a)

	

=

Ea[k]



Gk(x, u∗, a)

	

< Ea[k]



Gk(x, u∗, a)

	

, which violates the

fact that u∗ = arg min

u[k]

Ea[k][Gk(x, u, a)]. Similarly, to prove

statement 2, we assume such a u∗[k] = u∗[k] exists.

Lemma 1 leads to Ea[k]



Gk(x, u∗, a)

	

= Ea[k]



G

k(x, u∗, a)

	

<

Ea[k]



G

k(x, u∗, a)

	

, which violates the fact that u∗ =

arg min

u[k]

Ea[k][Gk(x, u, a)]. The results that u∗[k] = u∗[k] and

π∗ = π∗ naturally follow.

We can also apply M-PCM-OFFD to estimate the mean cost

Ea[k][Gk(x, u, a)] using fewer samples. Following similar procedures in the above proof and through replacing G

k(x, u, a)

according to Lemma 2 instead of Lemma 1 in the proof, we

can derive Theorem 2. The proof is omitted here for the sake

of space.

Theorem 2: Consider a dynamical system described by

Equation (4), with cost and value functions given by

Equation (5) and Equations (7)–(9), respectively. Then the

control policy π∗ = {μ∗

0,...,μ∗N−1} optimal to the samples selected by M-PCM-OFFD at each time step, is also

optimal to all possible values of the uncertain parameters, if

ni[k] = 2, ∀i ∈ {1, 2,..., m} and k ∈ {0, 1,...,N − 1}, and
XIE et al.: SCALABLE SAMPLING METHOD TO HIGH-DIMENSIONAL UNCERTAINTIES 101

j1,...,jm (x[k], u[k]) = 0 when more than τ of j1,...,jm are

non-zero, where 1 ≤ τ ≤ m.

B. Optimal Control for Systems of Uncertain Parameters

Evolving According to Markov Chains

Consider a generic dynamical system modulated by an mdimensional time-varying uncertain vector a[k]. Each element

of a[k], ai[k], evolves according to a Markov chain of transition matrix Pi ∈ RMi×Mi , where Mi is the total number of

possible values (or states) of the uncertain parameter ai[k].

The (j, l)th entry of Pi, Pijl , represents the probability for

ai[k] to be at state l in the next step if the current state is

j. The system dynamics and the total expected cost are shown

in Equations (4) and (5). Similar to Equation (6), given the

probability of ai[0] at each state for all i ∈ {1, 2,..., m},

consider the problem of finding an optimal control policy

π∗ = {μ∗

0, μ∗1,...,μ∗N−1} that minimizes the total expected

cost. This problem finds broad applications where dynamical

systems operate in uncertain environments. For instance, in

realistic systems such as the air traffic systems, uncertainties

like convective weather have correlations across time, which

can be captured by Markov chains.

Consider the i-th uncertain parameter, ai[k], in a system of

m uncertain parameters. Suppose initially ai[0] is at state j, i.e.,

ai[0] = j, then the conditional probability of ai[k] = l can be

obtained using p(ai[k] = l | ai[0] = j) = Pk

ijl . Suppose ri[k] ∈

R1×Mi is a vector of the possibilities of ai[k] at each state,

with its j-th element, denoted as rij[k], equal to the possibility

of ai[k] at state j, i.e., rij[k] = p(ai[k] = j). Then we have

ri[k] = ri[0]Pk

i . Since ri[0] and the transition matrix Pi are

given, we can calculate ri[k] at each time step k. Let us now

apply dynamic programming to solve this problem. Denoting

the set {1, 2,..., Mi} as Mi, the optimal value function is

V∗

k (x[k]) = min

u[k]

{Ea1[k]{· · · Eam[k]{gk(x, u)

+ αV∗

k+1(x[k + 1])} · · · }}

= min

u[k]

{Ea1[k]{· · · Eam−1[k]{



j∈Mm

p(am[k] = j)

{gk(x, u) + αV∗

k+1(x[k + 1])}} · · · }}

= min

u[k]

{Ea1[k]{· · · Eam−1[k]{



j∈Mm

rmj[k]

{gk(x, u) + αV∗

k+1[hk(x, u, a1,..., am = j)]}} · · · }}

= min

u[k]

{



l∈M1

r1j[k]{· · · { 

j∈Mm

rmj[k]{gk(x, u)

+ αV∗

k+1[hk(x, u, a1 = l,..., am = j)]}} · · · }} (10)

This problem then transforms to the first problem discussed

in Section III-A, with the pdfs of the uncertain parameters replaced by the probability mass functions (pmf) that

can then be solved by M-PCM or M-PCM-OFFD based

stochastic optimal controls. The sample points can be selected

according to the pdfs approximated by pmfs, or directly

using the sample-moment based approach [14]. The proof of

Theorem 3 is omitted as it can be easily derived following the

proof procedure of Theorem 1 with V∗

k (x[k]) represented by

Equation (10).

Theorem 3: Consider a dynamical system shown in

Equation (4), with each element of a[k], ai[k], evolves according to a Markov chain and the total expected cost shown

in Equation (5). The optimal control policy can be found

by sampling the uncertainty space using M-PCM under the

assumptions in Theorem 1, or using M-PCM-OFFD under the

assumptions in Theorem 2.

IV. REINFORCEMENT LEARNING-BASED

INFINITE HORIZON CONTROL

The uncertainty sampling methods can also be integrated

with reinforcement learning to bring the backward-in-time

optimal control to forward-in-time adaptive control for systems of high-dimensional uncertainties [5]–[7]. For the infinite

horizon case, the adaptive solution is also the optimal solution.

Consider a generic dynamical system described by

Equation (4), but with time horizon N → ∞ and ai[k] following a time-invariant pdf fAi(ai[k]). Consider the problem

of finding the optimal control policy π∗ = {μ, μ, . . .} that

minimizes the total expected cost:

J(x[k]) = Ea[k]{Ea[k+1]{· · · Ea[∞]{

∞

i=k

αi−kgi(x, u)} · · · }} (11)

The value function V(x[k]) = J(x[k]), and the Bellman’s

equation becomes

V(x[k]) = Ea[k]



gk(x, u) + αV(x[k + 1])




, (12)

where the same value function V() appears on both sides.

Therefore, given a control policy π, we can calculate the total

expected cost by solving Equation (12). The optimal control

policy

π∗ = arg min

π V(x[k]) (13)

can be found by repeatedly iterating two procedures policy

evaluation and policy improvement. In the policy evaluation

step, the value function V(x[k]) is solved using Equation (12)

for a set of admissible states x[k] ∈ S ⊆ S, given a current

control policy π. A best control policy is then derived based

on the value function V(x[k]) determined in the previous step

using Equation (13) in the policy improvement step. Various

approaches can be used to implement these two procedures,

and here we use value iteration. In particular, starting from an

arbitrary initial control policy π0, the following two steps are

iterated until convergence [5], [6]:

Value Update: ∀x[k] ∈ S

j ⊆ S

Vj+1(x[k]) = Ea[k]




gk(x, μj(x)) + αVj(x[k + 1])



(14)

Policy Improvement: ∀x[k] ∈ S

j ⊆ S

μj+1(x[k]) = arg min

μ(·)

Ea[k]




gk(x, μ(x)) + αVj+1(x[k + 1])



(15)

where j is the iteration step index. For systems of infinite

state and control spaces, approximators (e.g., neural network and polynomial functions) can be used to estimate

the value function by V(x) = WT(x), where (x) =

[φ1(x), φ2(x), . . . , φL(x))] is the basis vector, W is the weight
102 IEEE CONTROL SYSTEMS LETTERS, VOL. 1, NO. 1, JULY 2017

vector, and the estimation error approaches 0 when the number

of terms in (x), L → ∞. The value update then becomes

WT

j+1(x[k]) = Ea[k]




gk(x, u) + αWT

j (x[k + 1])



(16)

which can be solved using the least squares estimation.

Similarly, the control function can be estimated by μ(x) =

UT σ (x), where σ (x) and U are the basis and weight vectors,

respectively. For systems of x[k + 1] = A(x[k]) + B(x[k])u[k]

and gk(x, u) = Q1(x[k]) + uT [k]Q2u[k], where Q1(x[k]) > 0

and Q2 > 0, the policy improvement step can then be performed using the gradient descent method according to the

following equation [5], [6]:

Ui+1

j+1 = Ui

j+1 − βσ(x[k])



2Q2(Ui

j+1)

T σ (x[k])

+ Ea[k]



αB(x[k])

TT (x[k + 1])Wj+1

	
T

, (17)

where (x) = ∂(x)/∂x, and β is a tuning parameter and

i is the tuning index [5], [6]. For systems of finite state and

control spaces, value iteration can be achieved by storing and

updating lookup tables. Please refer to [5] and [6] for detailed

descriptions of different approaches.

We note that given an admissible state x[k] and control

function μj(x), the value function Vj+1(x[k]) in Equation (14)

can be approximated by the mean output of a system mapping Gk(x, a) = gk(x, μ(x)) + αVj(x[k + 1]), using M-PCM

or M-PCM-OFFD. Similarly, to derive μj+1(x) in the policy

improvement step, the expected cost of a candidate control

function can also be approximated by M-PCM or M-PCMOFFD. Theorem 4 describes the utilization of M-PCM and

M-PCM-OFFD to find the optimal control policy, the proof

of which can be easily derived based on above analysis and

proof of Theorem 1 and is thus omitted.

Theorem 4: Consider a dynamical system shown in

Equation (4), with time horizon N → ∞ and ai[k] following a

time-invariant pdf fAi(ai[k]). The optimal control policy can be

found by applying value iteration of reinforcement learning,

and approximating the value function using M-PCM under the

assumptions in Theorem 1, or using M-PCM-OFFD under the

assumptions in Theorem 2.

V. ILLUSTRATIVE EXAMPLES

In this section, we use three simple examples to illustrate

and validate the M-PCM and M-PCM-OFFD based stochastic

optimal controls. As the validation requires extensive Monte

Carlo simulations, we limit the dimension of uncertainties

in the first two examples. The third example illustrates the

capability of the M-PCM-OFFD based approach in addressing

high-dimensional uncertainties.

A. Dynamic Programming Based Finite-Horizon Control

Consider a system of the dynamics x[k + 1] = a1[k]x[k] +

a2[k]u[k] + a3[k], where ai[k] are uncertain parameters that

follow independent non-Gaussian distributions: fA1 (a1[k]) = 5

3

with −0.1 ≤ a1[k] ≤ 0.5, fA2 (a2[k]) = 5

3 with −0.2 ≤ a2[k] ≤

0.4 and fA3 (a3[k]) = 10

7 , 0 ≤ a3[k] ≤ 0.7. The optimal control

policy π∗ = {μ∗

0,...,μ∗N−1} is sought to minimize the total

expected cost given by Equation (5) with qN(x) = (x[N]−10)2,

gk(x, u) = (x[k] − 10)2 + u[k]

2, and α = 1. The state space

Fig. 1. The optimal control policies found by applying the four

approaches to sample the uncertainty space.

and control space are S = {−2, −1.9,..., 1.9, 2} and C =

{−1, −0.9,..., 0.9, 1}, respectively.

We apply M-PCM and M-PCM-OFFD (with τ = 1) based

stochastic optimal controls to find the optimal solution. Linear

interpolation is used to approximate costs at intermediate

points [7]. As a comparative study, we also apply the SG

with Gauss-Legendre quadrature rules to sample the uncertainty space [21]. For a fair comparison, the accuracy level1

of the SG is set to 2 such that a similar number of samples

is produced. The MC method is also applied to estimate the

true optimal solution for the validation purpose.

The optimal control policies found by the four methods are

shown in Figure 1, with N = 2. The number of samples n

selected by MPCM, MPCM-OFFD, SG and MC at each iteration are 8, 4, 7 and 100000, respectively. M-PCM-OFFD (or

M-PCM) based method finds accurate control solutions, with

significantly reduced number of simulations (4 or 8). The SG

method is relatively efficient but is less accurate.

B. Reinforcement Learning Based Infinite-Horizon

Control

Consider a system of the dynamic x[k + 1] = A[k]x[k] +

Bu[k] + Ca3[k], where x[k] = [x1[k], x2[k]]T , A[k] =

[a1[k], 0; 0, a2[k]], B = [1, 0.5]T , and C = [1, 1]T . a1[k], a2[k]

and a3[k] are uncertain parameters following fA1 (a1[k]) = 5

2

with 0.1 ≤ a1[k] ≤ 0.5, fA2 (a2[k]) = 2 with −0.5 ≤ a2[k] ≤ 0

and fA3 (a3[k]) = 10

3 , 0.2 ≤ a3[k] ≤ 0.5, respectively.

The optimal control policy π∗ is sought to minimize the

total expected cost given by Equation (11), with α = 0.8,

gk(x, u) = 8x2

1[k] + 2x22[k] + u[k]

2. Both state and control

spaces are continuous and infinite.

We use value iteration, with cost and control functions

approximated as polynomials, i.e., J(x) = WT(x) and

μ(x) = UT σ (x), where (x) = [1, x1, x2, x2

1, x22, x1x2]

T ,

1The accuracy level reflects the accuracy of the underlying quadrature rule [21]. Higher accuracy level requires more sample runs to reach

convergence.
XIE et al.: SCALABLE SAMPLING METHOD TO HIGH-DIMENSIONAL UNCERTAINTIES 103

Fig. 2. Illustration of the a) convergence of U and b) trajectory of total

expected cost J(x[0]) with x[0] = [0.1, 0.1]T .

σ (x) = [1, x1, x2]

T , W ∈ R6×1 and U ∈ R3×1 are weight

vectors. We iteratively update W and U vectors according

to Equations (16)–(17) (β = 0.1) until convergence, with

mean values Ea[k]() in Equations (16) and (17) estimated

using M-PCM, M-PCM-OFFD (τ = 1), and MC. The number of samples n used by M-PCM, M-PCM-OFFD, and MC

to estimate each mean value are 8, 4 and 10000, respectively.

Figure 2 shows the convergence of U along with the associated total expected cost J(x[0]) using M-PCM-OFFD. The

weight vectors derived by M-PCM (or M-PCM-OFFD) based

method are W = [1.076, 0.118, −0.133, 8.192, 2.136, 0.104]T

and U = [−0.329, −0.246, 0.028]T , which are close to

W = [1.076, 0.118, −0.130, 8.203, 2.130, 0.103]T and U =

[−0.327, −0.249, 0.026]T obtained by the MC based method.

When x[0] = [0.1, 0.1]T , the optimal total expected cost found

by M-PCM (or M-PCM-OFFD) and MC based methods are

1.175 and 1.180, respectively.

C. Control Under High-Dimensional Uncertainties

Consider a system modulated by 50 time-invariant

uncertain parameters and 2 time-invariant control

inputs. The system dynamics are described by x[k] =

0.2

10

i=1 ai −20i=11 ai

0 25

i=1 a2

i



x[k] + 0.1



0 0

0 50

i=26 ai



u,

where x[k] = [x1[k], x2[k]]T and u = [u1, u2]

T are state

and control vectors, respectively. ai, i ∈ {1, 2,..., 50},

is an uncertain variable that follows a uniform distribution fAi(ai) = 1, 0 ≤ ai ≤ 1. The state space

is x[k] ∈ R2 and the control space is of size 121:

u1 ∈ {1, 1.1,..., 1.9, 2} and u2 ∈ {−0.5, −0.4,..., 0.4, 0.5}.

The optimal control inputs u∗ are sought to minimize the following total expected cost: JN(x[0]) =

Ea[0]{· · · Ea[N−1][

N−1

k=0 αk(xT [k]Q1x[k] + uTQ2u)] ···}

where α = 0.8, N = 5, x[0] = [0, 1]T , Q1 =



2 0

0 3

and

Q2 =



0.5 0

0 1

.

With the high number of uncertain parameters, MC simulations are impractical to use. M-PCM based control also

requires 250 simulations to estimate the total expected cost for

each admissible control vector u. Theorem 2 proves the feasibility of M-PCM-OFFD in handling this high dimensionality.

Integrating M-PCM with 250−44

III OFFD (τ = 1), the M-PCMOFFD can further reduce the number of simulations from 250

to 26 = 64. The optimal control inputs and corresponding

total expected cost obtained by M-PCM-OFFD based control

are u∗ = [1, −0.5]T and J∗

N(x[0]) = 167.351, respectively.

VI. CONCLUSION

This letter develops two multi-dimensional uncertainty evaluation based approaches to address the scalability issue of

stochastic optimal control for systems of high-dimensional

uncertainties. For three cases covering varying scenarios,

we prove that the control solution optimal to the sampled

uncertainty space produced by M-PCM or M-PCM-OFFD is

also optimal to the original uncertainty space under simple

assumptions on the forms of the cost functions and orders

of uncertain parameters. Simulation and comparison studies demonstrate the accuracy and computational efficiency of

these two approaches.

REFERENCES

[1] C. Taylor, C. Wanke, Y. Wan, and S. Roy, “A framework for flow contingency management,” in Proc. AIAA Aviation Technol. Integr. Oper.

Conf. (ATIO), 2011, pp. 6904–6926.

[2] Y. Wan and S. Roy, “A scalable methodology for evaluating and designing coordinated air-traffic flow management strategies under uncertainty,” IEEE Trans. Intell. Transp. Syst., vol. 9, no. 4, pp. 644–656,

Dec. 2008.

[3] D. P. Bertsekas and S. E. Shreve, Stochastic Optimal Control: The

Discrete Time Case, vol. 23. New York, NY, USA: Academic Press,

1978.

[4] H. J. Kappen, “An introduction to stochastic control theory, path integrals

and reinforcement learning,” in Proc. Cooperative Behav. Neural Syst.,

vol. 887. 2007, pp. 149–181.

[5] F. L. Lewis and D. Vrabie, “Reinforcement learning and adaptive

dynamic programming for feedback control,” IEEE Circuits Syst. Mag.,

vol. 9, no. 3, pp. 32–50, 3rd Quart., 2009.

[6] F. L. Lewis, D. Vrabie, and K. G. Vamvoudakis, “Reinforcement learning and feedback control: Using natural decision methods to design

optimal adaptive controllers,” IEEE Control Syst., vol. 32, no. 6,

pp. 76–105, Dec. 2012.

[7] F. L. Lewis, L. Xie, and D. Popa, Optimal and Robust Estimation: With

an Introduction to Stochastic Control Theory, vol. 29. Hoboken, NJ,

USA: CRC Press, 2007.

[8] C. D. Phelps, “Computational optimal control of nonlinear systems with

parameter uncertainty,” Ph.D. dissertation, Dept. Appl. Math. Stat., Univ.

California, Santa Cruz, CA, USA, 2014.

[9] H. Chung, E. Polak, J. O. Royset, and S. Sastry, “On the optimal detection of an underwater intruder in a channel using unmanned underwater

vehicles,” Naval Res. Logistics, vol. 58, no. 8, pp. 804–820, 2011.

[10] J. C. Foraker, “Optimal search for moving targets in continuous time

and space using consistent approximations,” Ph.D. dissertation, Naval

Postgraduate School, Monterey, CA, USA, 2011.

[11] J. Ruths and S. Li, Jr., “Optimal control of inhomogeneous ensembles,”

IEEE Trans. Autom. Control, vol. 57, no. 8, pp. 2021–2032, Aug. 2012.

[12] Y. Matsuno, T. Tsuchiya, J. Wei, I. Hwang, and N. Matayoshi,

“Stochastic optimal control for aircraft conflict resolution under wind

uncertainty,” Aerosp. Sci. Technol., vol. 43, pp. 77–88, Jun. 2015.

[13] W. Liu and I. Hwang, “Probabilistic aircraft midair conflict resolution using stochastic optimal control,” IEEE Trans. Intell. Transp. Syst.,

vol. 15, no. 1, pp. 37–46, Feb. 2014.

[14] Y. Zhou et al., “Multivariate probabilistic collocation method for

effective uncertainty evaluation with application to air traffic flow management,” IEEE Trans. Syst., Man, Cybern., Syst., vol. 44, no. 10,

pp. 1347–1363, Oct. 2014.

[15] J. Xie et al., “Effective and scalable uncertainty evaluation for largescale complex system applications,” in Proc. IEEE Winter Simulat. Conf.

(WSC), Savannah, GA, USA, 2014, pp. 733–744.

[16] J. Xie et al., “Effective uncertainty evaluation in large-scale systems,”

in Principles of Cyber Physical Systems. Cambridge, U.K.: Cambridge

Univ. Press, 2017.

[17] G. E. Box, J. S. Hunter, and W. G. Hunter, Statistics for Experimenters:

Design, Innovation, and Discovery, vol. 2. New York, NY, USA: Wiley,

2005.

[18] Fractional Factorial Design Specifications and Design

Resolution. Accessed on May 15, 2017. [Online]. Available:

http://www.itl.nist.gov/div898/handbook/pri/section3/pri3344.htm

[19] S. S. Isukapalli, “Uncertainty analysis of transport-transformation models,” Ph.D. dissertation, Dept. Chem. Biochem. Eng., Rutgers Univ.-New

Brunswick, New Brunswick, NJ, USA, 1999.

[20] R. E. Bellman, Dynamic Programming. Princeton, NJ, USA: Princeton

Univ. Press, 1957.

[21] F. Heiss and V. Winschel, “Likelihood approximation by numerical integration on sparse grids,” J. Econometrics, vol. 144, no. 1, pp. 62–80,

2008.