Exploiting Active Subspaces to Quantify Uncertainty in the Numerical Simulation of the HyShot II Scramjet

Paul Constantine, Colorado School of Mines, [email protected]

Michael Emory, Stanford University, [email protected]

Johan Larsson, University of Maryland, [email protected]

Gianluca Iaccarino, Stanford University, [email protected]

Ryan Howard, Colorado School of Mines, [email protected]

In this notebook, we'll be summarizing the active subspaces-based uncertainty quantification presented in [1].

A scramjet is an engine designed to intake air and combust fuel at supersonic speeds. Because physically prototyping and experimenting is expensive, we must use numerical simulations to explore the behavior of these machines. The simulation we use is based on the HyShot II vehicle and has 7 input parameters:

Parameter Lower bound Upper bound
Angle of attack 2.6 4.6
Turbulence intensity 0.1 1.9
Turbulence length scale 122.5 367.5
Stagnation pressure 1.6448e7 1.9012e7
Stagnation enthalpy 3.0551e6 3.4280e6
Cowl transition location 0.03 0.07
Ramp transition location 0.087 0.203

These parameters are used in our complex multiphysics simulation modelling the fluid flow in the scramjet. From our simulation, we compute our quantity of interest, a normalized integral of pressure over the end of the engine:

$$ \frac{1}{vol(V)}\iiint\limits_VP\ dx\ dy\ dz $$

where V is a 3D region at the end of the scramjet and vol(V) is the volume of this region. This quantity can be used as a metric for judging performance and determining whether unstart has occured (unstart is a catastrophic failure resulting in loss of thrust). For each combination of parameters used in our simulation, we compute this quantity for fuel-air equivalence ratios of .30 and .35.

We'll use the notation $f(\mathbf{x})$ to denote the quantity of interest as a function of normalized input parameters ($\mathbf x$). For the purposes of uncertainty quantification, we'd like to estimate bounds of $f$, identify sets of inputs that result in $f$ remaining below a certain threshold, and estimate a cumulative distribution function of $f$ given a density on $\mathbf x$.

The most straightforward way to achieve our objective would be to exhaustively sample the space of input parameters and use the resulting data to calculate what we want. However, because our input space is high-dimensional and the simulation is very expensive, this is infeasible. To compute what we want, we reduce the dimensionality of the problem by identifying the 'most important' direction in the input space and examining our quantity of interest along this direction. We use an active subspace approach (see [2]) based on the matrix

$$ \mathbf C = \int\nabla f(\mathbf x)\nabla f(\mathbf x)^T \rho(\mathbf x)\ d\mathbf x = \mathbf W \Lambda\mathbf W^T $$

where $\rho$ is a probability density on $\mathbf x$ and $\mathbf W \Lambda\mathbf W^T$ is the eigendecomposition of $\mathbf C$. The active subspace is the space spanned by the first several eigenvectors, and we can use it to construct an approximation of $f$:

$$ f(\mathbf x)\approx g(\mathbf W_1^T\mathbf x) $$

where $\mathbf W_1$ contains the first $n$ columns of $\mathbf W$ and $g$ is a function from $R^n$ to $R$. Since we don't have access to the gradient of $f$, we use the following least squares-based algorithm to identify the single most important direction in the parameter space:

  1. Using $M$ data points, determine the parameters $\hat u_0, \hat u_1,\dots,\hat u_m$ of the linear approximation $f(\mathbf x)\approx \hat u_0 + \hat u_1x_1 + \cdots + \hat u_mx_m$ such that $\hat{\mathbf u} = argmin_{\mathbf u}||\mathbf{Au}-\mathbf f||_2^2$ where $$ \mathbf A = \left[\begin{matrix}1 & \mathbf x_1^T\\\vdots & \vdots\\ 1 & \mathbf x_M^T\end{matrix}\right],\ \mathbf f = \left[\begin{matrix}f_1\\\vdots\\f_M\end{matrix}\right],\ \hat{\mathbf u} = \left[\begin{matrix}\hat u_0\\\vdots\\\hat u_m\end{matrix}\right] $$
  2. Compute $\mathbf w = \hat{\mathbf u}'/||\hat{\mathbf u}'||$ where $\hat{\mathbf u}' = [\hat u_1, \cdots, \hat u_m]^T$. $\mathbf w$ is the vector pointing in the most important direction in the parameter space.

To examine the effects of variability in $\mathbf w$, we use a bootstrap algorithm:

  1. Choose the number $N$ of bootstrap replicates
  2. For $k$ from 1 to $N$:
    1. Let $\boldsymbol{\pi}^k = [\pi_1^k,\cdots,\pi_M^k]$ be integers from 1 to $M$ sampled uniformly with replacement.
    2. Let $\mathbf f^k$ be the vector whose $i^{th}$ element is $f_{\pi^k}$ and $\mathbf A^k$ be the matrix whose $i^{th}$ row is $[1,\ \mathbf x_{\pi^k}^T]$.
    3. Compute $\mathbf w^k$ using $\mathbf A = \mathbf A^k$ and $\mathbf f = \mathbf f^k$ as above.


[1] P.G. Constantine, M. Emory, J. Larsson, and G. Iaccarino. Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet. J. Comput. Phys., 302, 1-20, 2015

[2] P.G. Constantine, E. Dow, and Q. Wang. Active subspaces in theory and practice: applications to kriging surfaces. SIAM J. Sci. Comput., 36(4), A1500–A1524, 2014


This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC-0011077.

We'll now code up the analysis using the Python Active-subspaces Utility Library

In [1]:
%matplotlib inline
import numpy as np
import pandas as pn
import active_subspaces as ac
import matplotlib.pyplot as plt
np.set_printoptions(formatter={'all' : '{:<15}'.format})

First we import the data and separate inputs (X) from outputs (F)

In [2]:
df = pn.DataFrame.from_csv('HyShotII35.txt')
data = df.as_matrix()
#The last 2 rows (for the extreme values) are already normalized
X = data[:-2,0:7]
F = data[:,7:9]
labels = df.keys()
in_labels = labels[:7]
out_labels = labels[7:9]

#The data for equivalence ratio .30 is already normalized
df30 = pn.DataFrame.from_csv('HyShotII30.txt', header=None, index_col=None)
data30 = df30.as_matrix()
XX30 = data30[:,0:7]
F30 = data30[:,7]

Now we normalize each of the input parameters to the interval [-1, 1] and separate our outputs

In [3]:
#normalize the ratio = .35 data
xl = np.array([2.6, 0.1, 122.5, 1.6448e7, 3.0551e6, 0.6*0.05, 0.6*0.145])
xu = np.array([4.6, 1.9, 367.5, 1.9012e7, 3.4280e6, 1.4*0.05, 1.4*0.145])
XX = ac.utils.misc.BoundedNormalizer(xl, xu).normalize(X)
XX = np.vstack((XX, data[-2:,:7]))

f30 = F30[:]
out_label30 = out_labels[0]
f35 = F[:,1]
out_label35 = out_labels[1]

Set up our subspaces. The last two rows of each data set were from simulations done after initial subspace analysis was conducted, for the purpose of verifying predicted extrema; we exclude those data points initially to mirror the process of discovering the subspace.

In [4]:
ss30 = ac.subspaces.Subspaces()
ss35 = ac.subspaces.Subspaces()

#Compute OLS subspaces
ss30.compute(X=XX30[:-2,:], f=f30[:-2,np.newaxis], sstype='OLS', nboot=500)
ss35.compute(X=XX[:-2,:], f=f35[:-2,np.newaxis], sstype='OLS', nboot=500)

Plot the quantity of interest against $\mathbf w^T\mathbf x$ for each equivalence ratio

In [5]:
#active variable values associated with XX, XX30
y30 = np.dot(XX30, ss30.W1)
y35 = np.dot(XX, ss35.W1)

#plots of response (pressure) vs. active variable
ac.utils.plotters.sufficient_summary(y30[:-2], f30[:-2], out_label30)
ac.utils.plotters.sufficient_summary(y35[:-2], f35[:-2], out_label35)

We can see a very tight univariate trend in each case (see also Figure 7 of [1]); we can use this low-dimensional structure to perform our uncertainty quantification.

In [6]:
#Display the components of w
print "For ratio .30, the components of w are:"
print np.hstack((np.round(ss30.W1, 4), in_labels[:,np.newaxis]))
print "\nFor ratio .35, the components of w are:"
print np.hstack((np.round(ss35.W1, 4), in_labels[:,np.newaxis]))
For ratio .30, the components of w are:
[[0.6506          AoA            ]
 [0.5565          Turb int       ]
 [-0.0002         Turb len       ]
 [0.3685          Stag pres      ]
 [-0.3566         Stag enth      ]
 [-0.0196         Cowl trans     ]
 [0.0607          Ramp trans     ]]

For ratio .35, the components of w are:
[[0.6996          AoA            ]
 [0.4823          Turb int       ]
 [-0.027          Turb len       ]
 [0.1997          Stag pres      ]
 [-0.4738         Stag enth      ]
 [-0.0602         Cowl trans     ]
 [0.0957          Ramp trans     ]]

The components of $\mathbf w$ tell us how sensitive the quantity of interest is to each parameter. From the above output, we can see that Angle of Attack, Turbulence Intensity, and Stagnation conditions (pressure and enthalpy) dominate the active subspace, and (perhaps by coincidence) the components are roughly the same between the fuel-air equivalence ratios, except for an exchange of weight between stagnation pressure and stagnation enthalpy.

We'll now demonstrate the bootstrap for equivalence ratio .30.

In [7]:
N = int(1e5) #Number of bootstrap replicates
M,m = XX30.shape
w_boot = np.empty(shape=(m, N)) #matrix containing replicates
for k in range(N):
    mask = np.random.randint(0, M, size=(M,))
    f = np.matrix(f30[mask].reshape(M,1))
    A = np.matrix(np.hstack((np.ones((M,1)), XX30[mask])))
    u = (A.T.dot(A)).I.dot(A.T).dot(f)
    w_boot[:,k] = np.squeeze(u[1:])/np.linalg.norm(np.squeeze(u[1:]))
w_br = np.empty(shape=(m, 2)) #min-to-max range of bootstrap values
w_br[:,0] = np.squeeze(np.apply_over_axes(np.amin, w_boot, axes=(1)))
w_br[:,1] = np.squeeze(np.apply_over_axes(np.amax, w_boot, axes=(1)))
In [8]:
#plot eigenvector components with bootstrap ranges
ac.utils.plotters.eigenvectors(ss30.W1, w_br, in_labels, out_label30)
plt.figure(figsize=(14, 5*4))
for k in range(m): #plot histograms
    plt.subplot(4, 2, k+1)
    plt.xlim([-1, 1])
    plt.vlines(ss30.W1[k], 0, .5, 'r', lw=2)
    plt.hist(w_boot[k,:], bins=10,\

We can see the histograms are tightly concentrated around our estimated values (the red vertical lines), indicating that our estimates are close to their true values. Because we have only a few data points for the equivalence ratio of .35, the bootstrap ranges would be much wider in that case, but our confirmation of the results in the .30 case and the similarities between the two cases lead us to be confident for both ratios (see also Figures 8 and 9 of [1]).

Now that we are confident in our subspace, we can begin the uncertainty quantification (estimating ranges on $f$, identifying safe operating conditions, and estimating a cdf on $f$). Range estimation is sraightforward given the monotonic nature of the subspace: we can maximize $f$ by maximizing the value of our active variable, $\mathbf w^T\mathbf x$, over our input space, and running our simulation at these values of $\mathbf x$ (four runs: 1 max and 1 min for each equivalence ratio). We plot the initial data (blue dots) and the acquired extreme values (red squares) below (see also Figure 7 of [1]).

In [9]:
#Plot initial data (blue dots) and extrema (red squares)
plt.figure(figsize=(7, 7))
plt.plot(y30[:-2], f30[:-2], 'bo', y30[-2:], f30[-2:], 'rs', ms=12); plt.grid(True)
plt.xlabel('Active Variable'); plt.ylabel(out_label30)
plt.figure(figsize=(7, 7))
plt.plot(y35[:-2], f35[:-2], 'bo', y35[-2:], f35[-2:], 'rs', ms=12); plt.grid(True)
plt.xlabel('Active Variable'); plt.ylabel(out_label35)
<matplotlib.text.Text at 0x7f3763740f90>

Suppose now that the engine remained safe so long as exit pressure is below 2.8 and becomes unsafe when pressure rises above this. We might reasonably be interested in what regions of our parameter space result in safe operating conditions and which don't. To constrain the exit pressure, we can construct response surfaces from the data points $\{(y_i,f_i)\}$, where $y_i$ are the active variable values and $f_i$ are the values of the quantity of interest, and use them to place a bound on how high the active variable can become before the engine becomes unsafe, and map the safe active-variable region to the original parameter space.

Below, we plot our response surfaces (obtained from OLS on a quadratic model) and an upper prediction limit that is 2.33 standard errors above our mean response (see also Figure 10 of [1]). We also find the active variable value that marks the boundary between safe and unsafe operation. The shaded regions denote safe operation.

In [10]:
#Quadratic response surfaces
surface30 = ac.utils.response_surfaces.PolynomialApproximation(2)
surface35 = ac.utils.response_surfaces.PolynomialApproximation(2)

#Training the surfaces with our data
surface30.train(y30, f30[:,None])
surface35.train(y35, f35[:,None])

M30 = len(y30); M35 = len(y35)#Number of data points for each ratio

#sums of squares required for computing standard error
sf30 = np.sqrt(sum((ac.utils.misc.atleast_2d_col(f30) -\
                    surface30.predict(y30)[0])**2)/(M30 - 3))
sy30 = sum((y30 - y30.mean())**2)
#predicted response using mean response + 2.33 std errors
predlim30 = lambda x: surface30.predict(x)[0] + 2.33*sf30*np.sqrt(1 +\
                      1/M30 + (x - y30.mean())**2/sy30)

sf35 = np.sqrt(sum((ac.utils.misc.atleast_2d_col(f35) -\
                    surface35.predict(y35)[0])**2)/(M35 - 3))
sy35 = sum((y35 - y35.mean())**2)
predlim35 = lambda x: surface35.predict(x)[0] + 2.33*sf35*np.sqrt(1 +\
                      1/M35 + (x - y35.mean())**2/sy35)

x = ac.utils.misc.atleast_2d_col(np.linspace(-2, 3, 800))
#maximum Active variable value s.t. f(x) <= 2.8 according to prediction interval
avmaxval30 = x[np.argmin(abs(predlim30(x) - 2.8))]
plt.figure(figsize=(7, 7))
plt.plot(x, surface30.predict(x)[0], 'b-', x, predlim30(x), 'r:', lw=2)
plt.plot(x, 2.8*np.ones_like(x), 'k-')
plt.vlines(avmaxval30, 2, 3, 'k')
plt.fill_between(x[x <= avmaxval30], 2, 2.8, alpha=.3)
plt.ylim([2, 3])
plt.xlabel('Active Variable'); plt.ylabel(out_label30)
plt.legend(['mean response', 'prediction limit'], loc=2)

avmaxval35 = x[np.argmin(abs(predlim35(x) - 2.8))]
plt.figure(figsize=(7, 7))
plt.plot(x, surface35.predict(x)[0], 'b-', x, predlim35(x), 'r:', lw=2)
plt.plot(x, 2.8*np.ones_like(x), 'k-')
plt.vlines(avmaxval35, 2, 3.6, 'k')
plt.fill_between(x[x <= avmaxval35], 2, 2.8, alpha=.3)
plt.ylim([2, 3.6]); plt.xlim([-2, 3])
plt.xlabel('Active Variable'); plt.ylabel(out_label35)
plt.legend(['mean response', 'prediction limit'], loc=2)
<matplotlib.legend.Legend at 0x7f376350c450>

The set of normalized parameters resulting in safe operation is $S = \{\mathbf x : \mathbf w^T\mathbf x \leq y_{max},\ -\mathbf 1 \leq \mathbf x\leq \mathbf 1\}$, where $y_{max}$ is the maximum safe active variable value. We will find ranges of input parameters that guarantee safety by finding the largest 'box' that can fit within the set $S$ by solving the optimization problem: maximize$_{\mathbf x}\prod_{i=1}^m|x_i-x_{i,min}|,\ s.t.\ \mathbf x\in S$, where the $x_{i,min}$ are the components of xmin above (the minimizer of the active variable, $\mathbf w^T\mathbf x$).

In [11]:
import scipy.optimize as op

#Functions to minimize (-log(prod(x_i-x_{i,min})))
minfun30 = lambda x: -1.0*np.sum(np.log(abs(x - XX30[-2,:].reshape(x.shape))))
minfun35 = lambda x: -1.0*np.sum(np.log(abs(x - XX[-2,:].reshape(x.shape))))
#constraint functions of the form c(x) >= 0
cfun30 = lambda x: avmaxval30 - x.dot(ss30.W1)
cfun35 = lambda x: avmaxval35 - x.dot(ss35.W1)
#constraint dictionaries passed to minimizer
cdict30 = {'type': 'ineq', 'fun': cfun30}
cdict35 = {'type': 'ineq', 'fun': cfun35}

#Normalized parameter values defining the edge of the largest box assuring safe operation
xop30 = op.minimize(minfun30, XX30[-1,:],\
                    bounds=zip(-1.0*np.ones_like(xl), 1.0*np.ones_like(xl)),\
xop35 = op.minimize(minfun35, XX[-1,:],\
                    bounds=zip(-1.0*np.ones_like(xl), 1.0*np.ones_like(xl)),\

#un-normalizing inputs to original parameter space
xop30 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(xop30.x[None,:]).reshape(m, 1)
xop35 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(xop35.x[None,:]).reshape(m, 1)
min30 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(XX30[-2,:][None,:]).reshape(m, 1)
min35 = ac.utils.misc.BoundedNormalizer(xl, xu).unnormalize(XX[-2,:][None,:]).reshape(m, 1)

#sorting bounds for display convenience
A30 = np.round(np.sort(np.hstack((min30, xop30)), axis=1), 4)
A35 = np.round(np.sort(np.hstack((min35, xop35)), axis=1), 4)
In [12]:
#Display the safe operating ranges
print "For ratio .30, the ranges providing safe operation are:"
for i in range(7):
    print "{:0.2e} \t {:0.2e}\t {}".format(A30[i,0],A30[i,1],in_labels[i])
print "\nFor ratio .35, the ranges providing safe operation are:"
for i in range(7):
    print "{:0.2e} \t {:0.2e}\t {}".format(A35[i,0],A35[i,1],in_labels[i])
For ratio .30, the ranges providing safe operation are:
2.60e+00 	 4.60e+00	 AoA
1.00e-01 	 1.90e+00	 Turb int
1.22e+02 	 3.67e+02	 Turb len
1.64e+07 	 1.90e+07	 Stag pres
3.06e+06 	 3.43e+06	 Stag enth
3.00e-02 	 7.00e-02	 Cowl trans
8.70e-02 	 2.03e-01	 Ramp trans

For ratio .35, the ranges providing safe operation are:
2.60e+00 	 3.25e+00	 AoA
1.00e-01 	 9.53e-01	 Turb int
1.22e+02 	 3.67e+02	 Turb len
1.64e+07 	 1.90e+07	 Stag pres
3.25e+06 	 3.43e+06	 Stag enth
3.00e-02 	 7.00e-02	 Cowl trans
8.70e-02 	 2.03e-01	 Ramp trans

Lastly, we can construct a cumulative distribution of $f$ by randomly sampling the (normalized) parameter space, computing the active variable, and computing the response surface value at each point (see also Figure 11 of [1]).

In [13]:
N = int(1e5) #Number of random samples
x = np.random.uniform(-1, 1, (N, m))
a30 = x.dot(ss30.W1); a35 = x.dot(ss35.W1) #active variable values
s30 = surface30.predict(a30)[0]; s35 = surface35.predict(a35)[0] #response surface values
x = np.linspace(1.5, 3.5, 300)
c30 = np.vectorize(lambda x: sum(s30 <= x)/(1.0*N)) #CDF functions
c35 = np.vectorize(lambda x: sum(s35 <= x)/(1.0*N))

plt.plot(x, c30(x), 'k-')
plt.ylim([0, 1.1])
plt.vlines(f30[-2:], 0, 1.1, 'k', 'dashed')

plt.plot(x, c35(x), 'k-')
plt.ylim([0, 1.1])
plt.vlines(f35[-2:], 0, 1.1, 'k', 'dashed')
<matplotlib.collections.LineCollection at 0x7f37637852d0>