Junfei Xie, Yan Wan, Yi Zhou, University of North Texas

Kevin Mills, James J. Filliben, NIST

Yu Lei, University of Texas at Arlington

Winter Simulation Conference, 2014
 Modern large-scale cyber-physical systems (CPSs) involve a large 

number of uncertain parameters.

Uncertainties

 Modulate system’s 

dynamics

 Pose significant 

challenges for real-time 

system evaluation &

decision-support

Management of physical dynamics must be 

designed in a way to achieve robust performance 

under the uncertainties

Effective 

uncertainty 

evaluation

2
 Problem Formulation

Goal: correctly 

estimate the 

mean output

Air Traffic Flow Management

Uncertain input parameters:

total traffic delay

of aircraft over a 

time span

 Example Application

3

Traffic System

Weather start times

Weather durations

Weather intensities

Output: 
 Monte Carlo Simulation Method

High Computational Cost !!

4
Original System Mapping: Low-order Mapping:

5

Reduce the number of simulations Predict the correct mean output!

from

2

m ni

m

i=1

 ni

m

i=1

 Multivariate Probabilistic Collocation Method (M-PCM) 

Y. Zhou, Y. Wan, S. Roy, C. Taylor, C. Wanke, D. Ramamurthy, J. Xie, “Multivariate Probabilistic Collocation Method for Effective 

Uncertainty Evaluation with Application to Air Traffic Management”, IEEE Transactions on Systems, Man and Cybernetics: 

System, Vol. 44, No. 10, pp.1347-1363, 2014.
 Possibility of Further Reduction

 M-PCM assumes that there exist cross-multiplication terms for all 

combinations of uncertain parameters of all degrees

 Suppose ni = 1, i = 1,2, . . , m

g x1, x2, ... , xm = a0 + a1x1 + ⋯ + amxm + am+1x1x2 + ⋯ + aNx1x2 ... xm

Computational load issue for 

real-time applications m=2

m=100

2

2

2

100

Number of 

parameters 

Number of 

simulations

Suppose ni = 2, i = 1,2, . . , m

Number of simulations: ni

m

i=1 =2

m

Still too large !!

Some of these terms may not exist (ai ≈ 0

or ai = 0 ) in realistic applications

 Limitation of M-PCM: Not scalable with the number of parameters

6
 Approach: Integration of M-PCM with the orthogonal fractional 

factorial design (OFFD)------ M-PCM-OFFD

 OFFDs meet our need to reduce the number of simulations

 Both OFFDs and our study are motivated by the same assumption

 ---high-order interactions among parameters are insignificant in real 

applications

Numerical 

Truncation

Simulation ...

points

mean output 

prediction

Complex system

Run 

simulation

This is possible !

 Challenge: existence of a practical numerical issue

 Many system simulations have constraints on the resolutions of 

input parameters

7

Incorrect
 Full Factorial Designs

8

‘−’ lower level

‘+’ higher level

x1, x2, x3 are factors 

(input parameters)

y is the output

2

3

full factorial design

 Orthogonal Fractional Factorial Designs 

(OFFDs)

 Selects a subset of experimental combinations that 

best estimate the main effects of single factors and 

low-order interaction effects.

All possible combinations of levels of all factors.

2

3−1 OFFDs
 Main Effect and Interaction

Main effect MEi of factor xi

Interaction effect MEij of xixj

Regression Model:

y = β0 + βixi

k

i=1

+ βi.jxixj

k

j=i+1

+ε

k−1

i=1

The least square estimators 

for β, denoted as β are:

β 

i =

1

2

MEi, β 

ij =

1

2

MEij

9

e.g., ME3 =

y5+y6+y7+y8

4

−

y1+y2+y3+y4

4

ME23 =

y1 + y2 + y7 + y8

4

−

y3 + y4 + y5 + y6

4
 The main effects and interactions estimated by the subset of simulations 

selected by OFFDs are aliased.

I = x1x2x3

x3 = x1x2

x1 = x2x3

x2 = x1x3

β0 + β1.2.3

β1 + β2.3

β2 + β1.3

β3 + β1.2

M E1 = ME1 + ME2.3

M E2 = ME2 + ME1.3

M E3 = ME3 + ME1.2

10

x3 = x1x2

Generator: I = x1x2x3

Regression Model:

y = β0 + βixi

k

i=1

+ βi.jxixj

k

j=i+1

+ε

k−1

i=1

These are what 

OFFD estimates
11Generate the P

m−γ

full factorial design 

for m − γ factors

Determine the levels of all 

other γ factors for each 

experimental run

Specify γ

generators

e.g., 2

3−1 OFFD, P = 2 ,m = 3, γ = 1,

2

2

full factorial design 2

3−1 OFFD

Step 1 Step 2 Step 3

Generator 

I = −x1x2x3

x3 = −x1x2

‘−’ lower level

‘+’ higher level 11

P: the number of levels

m: the number of factors

γ: the fractionation constant
 If we view all simulation points selected by M-PCM as a full factorial 

design, the OFFDs provide systematic procedures to select a subset of 

simulation points.

 Design Procedures

Choose 2

m

M-PCM 

simulation 

points 

Produce 

the loworder 

mapping

Run simulation at 

selected points

Original System Mapping: Low-order mapping:

Check if 

1≤ τ ≤

m

2

− 1

Calculate γ and 

select simulation 

subset using 2

m−γ

OFFD

Assumption: cross-terms 

involve at most τ parameters
 Lemma 2:

If the original system mapping 

g(x1, x2, ... , xm) contains crossterms of at most τ parameters

The low-order mapping 

g

∗

(x1, x2, ... , xm) also contains crossterms of at most τ parameters

If 1≤ τ ≤

m

2

− 1

2

m−γmax OFFD can further reduce the number of 

simulations from 2

m to 2m−γmax

, 

where γmax = m − log2( 

i

m

τ

i=0

) 

Lemma 1

Lemma 2

Lemma 3

The matrix L ∈ R

loffd×l

constructed by the M-PCM-OFFD is full column 

rank, and can be represented by L = QU, where Q ∈ R

loffd×l

is an 

orthogonal matrix and U ∈ R

l×l

is an upper triangular matrix. 

Constructed by selected simulation points to 

estimate the coefficients of the low-order 

mapping. 

13

 Lemma 1:

 Lemma 3:
 Problem Formulation

 The M-PCM-OFFD involves the calculation of L

−1 or (LTL)−1LT

 L must be full column rank 

 Numerical errors may easily push L to lose rank and fail the computation

 To facilitate the calculation and minimize the impact of such numerical errorinduced disturbances, L needs to have a large margin to rank loss.

where e ∈ R

loffd×l

is a perturbation matrix

We proved that L matrix obtained using OFFD, denoted as Loffd, has the 

largest margin to rank loss, among all designs of the same size. 

Guaranteed using OFFD

 Metric: full-column-rank margin 

 The full-column rank margin for matrix L to rank loss is 

14

D L = min e F rank L + e < l}
 Original Mapping:

x1~ fX1x1 = 2e

−2x1; 

x2 ~ fX2x2 =

1

15

, 5 ≤ x2 ≤ 20; 

x3 ~ fX3x3 =

1

5

, 5 ≤ x2 ≤ 10; 

p1 = (0.2929,8.1699,6.0566),

p2 = (1.7071,8.1699,6.0566), 

p3 = (0.2929,16.8301,6.0566),

p4 = (1.7071,16.8301,6.0566), 

p5 = (0.2929, 8.1699,8.9434),

p6 =(1.7071,8.1699,8.9434), 

p7 =(0.2929,16.8301,8.9434),

p8 =(1.7071,16.8301,8.9434)

Step 1: Choose 8 M-PCM points based 

on the pdf of each parameter

Step 2: Use 2III

3−1 OFFD to select 4 

M-PCM points 

{p2, p3, p5, p8}

{p1, p4, p6, p7}

or

2III

3−1 OFFD design table

Step 3: Run simulations to evaluate 

g(x1, x2, x3)at these 4 M-PCM points 

Step 4: Estimate the coefficients of 

the low-order mapping

 Illustration of Design Procedures

15
 Illustration of Performance

Estimation of Mean Output

E g x1, x2, x3 = E g

∗ x1, x2, x3 = 3381.1

Robustness to Numerical Errors

D Loffd = 1.4142 D L = {0,0.866,1.4142}

max D(L) = D Loffd Selected by OFFD

Other possible selections

16
 An effective and scalable uncertainty evaluation method for 

large-scale complex systems

M-PCM-OFFD

reduces the number of simulations from 2

2m to at most 

2

log2(m+1)

accurately predicts the output mean under broad 

assumptions

is the most robust to numerical errors compared with 

designs of the same size

 In the future work

 Generalize the degree of uncertain input parameters by 

exploring multiple-factor OFFDs 

 Exploit parameter dependency to further reduce the number of 

simulations required.

 New interpretations of the optimality of OFFDs 

17
18

We thank National Institute of Standards and Technology for the support.