Chance-Constrained Yield Optimization of Photonic IC With Non-Gaussian Correlated Process Variations

Download as pdf or txt
Download as pdf or txt
You are on page 1of 12

1

Chance-Constrained Yield Optimization of Photonic


IC with Non-Gaussian Correlated Process Variations
Chunfeng Cui? , Kaikai Liu? , Zheng Zhang, Member, IEEE

Abstract—Uncertainty quantification has become an efficient optimization solvers for photonic IC. The major difficulties
tool for yield prediction, but its power in yield-aware optimization include: 1) the quantity of interest (e.g., the probability distri-
arXiv:1908.07574v1 [math.OC] 20 Aug 2019

has not been well explored from either theoretical or application bution of a bandwidth) does not admit an explicit expression.
perspectives. Yield optimization is a much more challenging task.
On one side, optimizing the generally non-convex probability Instead, we only know the simulation values at parameter sam-
measure of performance metrics is difficult. On the other side, ple points; 2) the design objectives and constraints are defined
evaluating the probability measure in each optimization iteration in a stochastic way. They are hard to compute directly and re-
requires massive simulation data, especially when the process quire massive numerical simulations to estimate their statistical
variations are non-Gaussian correlated. This paper proposes distributions; 3) practical photonic IC designs often involve
a data-efficient framework for the yield-aware optimization of
the photonic IC. This framework optimizes design performance non-Gaussian correlated process variations, which are more
with a yield guarantee, and it consists of two modules: a difficult to capture. To estimate the design yield efficiently,
modeling module that builds stochastic surrogate models for one alternative is to build a surrogate model. In [12]–[14],
design objectives and chance constraints with a few simulation posynomials were used to model statistical performance, and
samples, and a novel yield optimization module that handles geometric programming was employed to optimize the worst-
probabilistic objectives and chance constraints in an efficient
deterministic way. This deterministic treatment avoids repeatedly case performance. The reference [15] proposed a Chebyshev
evaluating probability measures at each iteration, thus it only affine arithmetic method to predict the cumulative distribution
requires a few simulations in the whole optimization flow. We function. The recent Bayesian yield optimization [10] approx-
validate the accuracy and efficiency of the whole framework by imated the probability density of the design variable under the
a synthetic example and two photonic ICs. Our optimization condition of “pass” by a kernel density estimation. The work
method can achieve more than 30× reduction of simulation cost
and better design performance on the test cases compared with [11] further approximated the yield over the design variables
a recent Bayesian yield optimization approach. directly by a Gaussian process regression. However, these
machine learning techniques may still require many simulation
Index Terms—Photonic integrated circuits, photonic design au-
tomation, uncertainty quantification, yield optimization, chance samples. Furthermore, only optimizing the yield can lead to
constraints, non-Gaussian correlations. nonoptimal (and even poor) chip performance.
Recently, uncertainty quantification methods based on gen-
eralized polynomial chaos have gained great success in mod-
I. I NTRODUCTION
eling the uncertainty caused by various process variations

T HE demand for low-power, high-speed communications


and computing have boosted the advances in photonic
integrated circuits. Based on the modern nano-fabrication
in electronic and photonic ICs [16]–[27]. A novel stochas-
tic collocation approach was further proposed in [28], [29]
to handle non-Gaussian correlated process variations, which
technology, hundreds to thousands of photonic components shows significantly better accuracy and efficiency than [30]
can be integrated on a single chip [1], [2]. However, process due to the smooth basis functions and an optimization-based
variations persist during all the fabrication processes and can quadrature rule. These techniques can construct stochastic sur-
cause a significant yield degradation in large-scale design and rogate models with a small number of simulation samples, but
manufacturing [3]–[6]. Photonic ICs are more sensitive to their power in yield optimization has not been well explored
process variations (e.g., geometric uncertainties) due to their or exploited despite recent robust optimization methods [31]
large device dimensions compared with the small wavelength. based on generalized polynomial chaos.
To achieve an acceptable yield, uncertainty-aware design op- Paper Contributions. Leveraging the chance-constrained
timization algorithms are highly desired [7]. optimization [32] and our recent uncertainty quantification
Yield optimization algorithms try to increase the success solvers [28], [29], this paper presents a data-efficient technique
ratio of a chip under random process variations, and they have to optimize photonic ICs with non-Gaussian correlated process
been studied for a long time in the electronic circuit design variations. Instead of just optimizing the yield, we optimize a
[8]–[11]. However, it is still expensive to reuse existing yield target performance metric while enforcing the probability of
?
violating some design rules to be smaller than a user-defined
C. Cui and K. Liu contributed equally to this work.
This work was partly supported by NSF Grant No. 1763699, NSF CAREER threshold. Doing so can avoid performance degradation in
Award No. 1846476, a UCSB start-up grant and a Samsung Gift Funding. yield optimization. Chance-constrained optimization [32] has
Chunfeng Cui, Kaikai Liu, and Zheng Zhang are in the Department of been widely used in system control [33], autonomous vehicles
Electrical and Computer Engineering, University of California, Santa Barbara,
CA 93106, USA (e-mail: chunfengcui@ucsb.edu, kaikailiu@ucsb.edu, and [34], and reliable power generation [35], [36], but it has
zhengzhang@ece.ucsb.edu). not been investigated for yield optimization of electronic or
2

photonic IC. Our specific contributions include: The yield optimization problem aims to find the optimal design
• A new yield optimization model that can handle non- variables x∗ such that
Gaussian correlated process variations and optimize chip x∗ = argmax Prob(S|x). (3)
performance while tolerating a certain failure probability. x∈X
Suppose that the yield requirement is represented as There are three major difficulties in solving the yield
several inequalities, we model these inequalities as a set optimization problem: 1) the indication function I(x, ξ) does
of chance constraints. not always admit an explicit formulation; 2) computing the
• A framework that approximates the stochastic objective yield Prob(S|x) involves a non-trivial numerical integration,
and constraint functions with a few simulations. In many which requires numerous simulations at each design variable
cases, both the objective function and constraints are only x; 3) Prob(S|x) is an implicit non-convex function and it is
available through an expensive black-box simulator. To difficult to get an optimal solution.
reduce the simulation time, we build a surrogate model
based on the recent uncertainty quantification solver [28], B. Chance Constraints
[29]. We propose a new three-stage process: firstly we The chance constraint is a powerful technique in
generate quadrature points for the design variables; sec- uncertainty-aware optimization [32]. In comparison with the
ondly, we compute the quadrature points for the uncer- deterministic constraints or the worst-case constraints where
tainty parameter; thirdly we re-optimize the quadrature the risk level  is zero, a chance constraint enforces the
points in their joint variable space. probability of satisfying a stochastic constraint to be above
• A deterministic reformulation. A major challenge of a certain confidence level 1 −  ( is usually not zero):
chance-constrained problems is to reformulate the
Probξ (y(x, ξ) ≤ 0) ≥ 1 −  (4)
stochastic constraints into deterministic ones [37]. We
reformulate the probabilistic objective function and con- or equivalently, the probability of violating the constraint to
straints as non-smooth deterministic functions. Afterward, be smaller than the risk level :
we transform the non-smooth deterministic formula into Probξ (y(x, ξ) ≥ 0) ≤ . (5)
an equivalent polynomial optimization problem, which
can be solved efficiently. Under strict conditions, such as the parameters being in-
• Validations on benchmarks. Finally, we validate the ef- dependent and y(x, ξ) being a linear function, (4) can be
ficiency of our proposed framework on a synthetic ex- reformulated into equivalent deterministic constraints [39].
ample, a microring band-pass filter, and a Mach-Zehnder In other words, one can reformulate the left-hand side of
filter. Preliminary numerical experiments show that our (4) by its probability density function (PDF) and substitute
proposed framework can find the optimal design variable the right-hand side by a constant related to the cumulative
efficiently. Compared with the Bayesian yield optimiza- density function (CDF). However, these conditions rarely hold
tion method [10], our proposed method can reduce the in practice. Even if the conditions hold, computing the PDF or
simulation times by 30× on the test cases and achieve CDF of an uncertain variable can be intractable [17], [37]. In
better performance while producing the same yield. these cases, we seek for deterministic reformulations that can
well approximate the chance constraints. There is a trade-off in
choosing the reformulation: if the reformulation is aggressive
II. P RELIMINARIES (the feasible domain is enlarged), it may result in an infeasible
solution; Otherwise, if the reformulation is conservative (the
A. The Yield Optimization feasible domain is decreased), the solution may be degraded.
The yield is defined as the percentage of qualified products One popular method of transforming (4) to a deterministic
overall. For a given photonic IC, denote the design variables constraint is to use the first and second moments of y(x, ξ)
by x = [x1 , x2 , ..., xd1 ]T ∈ X and the process variations [37], [39]:
by random parameters ξ = [ξ1 , ξ2 , ...ξd2 ]T ∈ Ω. Suppose x
q
Eξ [y(x, ξ)] + κ varξ [y(x, ξ)] ≥ 0. (6)
is uniformly distributed in a bound domain and ξ follows a
probability distribution ρ(ξ). Let {yi (x, ξ)}ni=1 denote a set Here Eξ [·] denotes the mean value, varξ [·] p denotes the vari-
of performance metrics of interest and I(x, ξ) denotes an ance. The constant κ is chosen as κ = (1 − )/. The
indicator function of success: detailed proof is shown in Appendix A. It is worth noting that
 (6) is a stronger condition than (4): every feasible point of (6)
1, if every yi (x, ξ) satisfies the requirements; is also a feasible point of the original chance constraint (4).
I(x, ξ) =
0, otherwise.
(1) C. Stochastic Spectral Methods
The event “pass” is defined as the events that satisfies
Assume that y(ξ) is a smooth function satisfying
I(x, ξ) = 1. Then the probability of S is the yield of the
E[y 2 (ξ)] ≤ ∞. The stochastic spectral methods can approxi-
system. Given a certain design choice x, its conditional yield
mate y(ξ) by orthonormal polynomial basis functions:
is defined as [38]
p
X
Z
y(ξ) ≈ cα Ψα (ξ), with E [Ψα (ξ)Ψβ (ξ)] = δα,β . (7)
Prob(S|x) = Eξ [I(x, ξ)] = I(x, ξ)ρ(ξ)dξ. (2)
Ω |α|=0
3

Here |α| = α1 + . . . + αd2 , Ψα (ξ) is an orthonormal basis Build the chance Input the range of
function indexed by α, and cα is its corresponding coefficient. constrained model (12) x and the PDF of ξ
If the parameters ξ are independent, ρ(ξ) equals to the
products of its one-dimensional marginal density function Reformulate (12) into Solve (28) to compute
ρi (ξi ). In this case, the basis function Ψα (ξ) is the product (14) with n constraints the quadrature rule
of multiple one-dimensional orthogonal basis functions
Reformulate (14) into Call the simulator at
Ψα (ξ) = ψ1 (ξ1 ) . . . ψd2 (ξd2 ). (8) deterministic model (17) the quadrature points

These one-dimensional basis functions ψi (ξi ) can be con-


Derive the polynomial Construct the surrogate
structed by the three term recursion [40]. Various stochastic
optimization model (31) model by (25)
spectral approaches have been proposed to compute the coef-
ficients cα , including the intrusive (i.e., non-sampling) solvers Sove (31) and output
(e.g., stochastic Galerkin [41], the stochastic testing [16]) the optimal design
and the non-intrusive (i.e., sampling) solvers (e.g., stochastic
collocation [42]). In the past few years, there has also been Fig. 1. The flowchart of our proposed framework for solving the chance
a rapid progress in handling high-dimensional parameters, constrained yield-aware optimization.
such as the tensor recovery method [19], the compressive
sensing technique [43], ANOVA (analysis of variance) or
HDMR (the high-dimensional model representation) [44], and However, the above yield maximization often contradicts
the hierarchical uncertainty quantification [18]. with our performance goals. For instance, one may have to
In practice, the random parameters may be correlated. If the reduce the clock rate of a processor significantly in order
parameters ξ are non-Gaussian correlated, the computation is to achieve a high yield. As a result, directly optimizing the
more difficult. In such cases, Ψα (ξ) can be constructed by yield may lead to an over-conservative design. Therefore, we
the Gram-Schmidt approach in [28], [29] or the Cholesky fac- instead optimize the expected value of an uncertain objective
torization in [45], [46]. The main difficulty lies in computing performance metric f (x, ξ) subject to a yield constraint:
high order moments of ξ, which can be well resolved by the min Eξ [f (x, ξ)]
functional tensor train approach proposed in [46]. x∈X
s.t. Y (x) = Probξ (y(x, ξ) ≤ u) > 1 − . (12)
III. O UR Y IELD - AWARE O PTIMIZATION M ODEL
Here  is a risk level to control the yield. The above formu-
In this section, we will show our yield optimization model lation can describe, for instance, the following design opti-
defined by a stochastic measure in the probability space, and mization problem: minimize the average power consumption
will illustrate how to convert the stochastic formulation to a of a photonic IC while ensuring a 95% yield (i.e., with 5%
deterministic one. We first present the basic assumptions in probability of violating timing and bandwidth constraints)
this paper. under process variations. Note that f (x, ξ) may also be the
Assumption 1. We made the following assumptions: function (e.g., weighted sum) of several performance metrics
that we intend to optimize simultaneously.
1) The design variable x ∈ X = [a, b]d1 follows a mutually
When the yield function Y (x) and the objective function
independent uniform distribution;
f (x, ξ) are available, we may solve the above optimization
2) The stochastic parameter ξ ∈ Ω ∈ Rd2 admits a non-
problem directly. Unfortunately, this is rarely true. Normally,
Gaussian correlated density function ρ(ξ);
one has to estimate the yield and objective at a given x by the
3) The yield is qualified by the following constraints:
Monte Carlo method [8], [9] which requires a huge number of
yi (x, ξ) ≤ ui , ∀ i ∈ [n]. (9) simulation samples at each design variable x. This is infeasible
for many simulation-expensive photonic IC design problems.
Here [n] = 1, . . . , n and E[yi (x, ξ)] ≤ ui . The quan-
Instead, we reformulate the yield as chance constraints and
tity {yi (x, ξ)}ni=1 are black-box functions: we can only
aim to propose a more data-efficient optimization framework.
obtain their function values at the given sample points.
Specifically, we set a risk threshold as i , and transform the
yield objective into the following constraints
A. The Probabilistic Yield Optimization Model
Probξ (yi (x, ξ) ≤ ui ) ≥ 1 − i , ∀ i ∈ [n]. (13)
The yield at a given design variable x can be defined as the
probability that the yield conditions (9) are satisfied, i.e., In this formulation i means the risk tolerance of violating
the i-th design specification. Any feasible point of the above
Y (x) = Probξ (y(x, ξ) ≤ u). (10)
constraint will have a high yield when i is small. Higher yield
T
Here, y(x, ξ) = [y1 (x, ξ), . . . , yn (x, ξ)] and u = can be obtained by decreasing i ’s. Consequently, we have the
[u1 , . . . , un ]T . Consequently, the yield optimization problem following chance-constrained yield-aware optimization model
can be described as:
min Eξ [f (x, ξ)]
x∈X
max Probξ (y(x, ξ) ≤ u). (11)
x∈X s.t. Probξ (yi (x, ξ) ≤ ui ) ≥ 1 − i , ∀ i ∈ [n]. (14)
4

B. From the Stochastic to Deterministic Model and p


X
The chance-constraint optimization in problem (14) is diffi- f (x, ξ) ≈ hα,β Φα (x)Ψβ (ξ). (19)
cult to solve directly. This problem is more challenging when |α|+|β|=0
yi (x, ξ) is nonlinear. In this case, it is almost impossible
Once the above surrogate models are obtained, the mean value
to formulate the chance constraints in (14) to equivalent
of yi (x, ξ) can be approximated by
deterministic formulations directly. A directly approach is to
p
replace the stochastic constraints by inequality constraints over X
Eξ [yi (x, ξ)] ≈ ciα,0 Φα (x), (20)
the expected constraints:
|α|=0
min Eξ [f (x, ξ)] and the variance is approximated by
x∈X
2
s.t. Eξ [yi (x, ξ)] ≤ ui , ∀ i ∈ [n].

(15) p
X p−|β|
X
varξ [yi (x, ξ)] ≈  ciα,β Φα (x) . (21)
However, this treatment will lose the probability density infor-
|β|=1 |α|=0
mation and may not provide a high-quality solution, although
it can help improve the yield in practice. We will illustrate The mean value of the objective function f (x, ξ) can be
this phenomenon in numerical experiments section V-A. evaluated in the same way. Finally, the deterministic yield
Therefore, we adopt the second-order moment approach in optimization model (17) can be solved. The overall framework
[37], [39] and change (13) to is summarized in Algorithm 1. In the following, we explain
q the implementation details.
Eξ [yi (x, ξ)] + κi varξ [yi (x, ξ)] ≤ ui , ∀ i ∈ [n]. (16)
q A. Basis Functions for Design and Uncertain Variables
1−i
Here, κi = i is a scaling parameter. We present the For the uniform-distributed design variables x, their basis
detailed proof in Appendix A. We point out the following: functions Φα (x) can be decoupled into the products of 1D
• Constraint (16) is a stronger condition than (13). In other basis functions:
words, each feasible point of (16) is also a feasible Φα (x) = φ1α1 (x1 ) . . . φdα1d (xd1 ). (22)
solution of the chance constraint (13); 1

• The parameter i is a user-defined risk tolerance. When i Here, φiαi (xi ) is a Legendre polynomial [47], which can be
decreases, the feasible set will become smaller. However, constructed by the three-term recurrence relation [40].
the optimal solution may result in a higher yield; For the random vector ξ describing non-Gaussian correlated
• When the variance varξ [yi (x, ξ)] is small enough, the process variations, we construct its basis functions Ψβ (ξ) by
feasible set of (16) is close to the deterministic constraint the Gram-Schmidt approach proposed in [28], [29]. Specifi-
β
Eξ [yi (x, ξ)] ≤ ui . cally, we first reorder the monomials ξ β = ξ1β1 . . . ξd2d2 in the
Np
Consequently, the probabilistic optimization model (14) is graded lexicographic order, and denote them as {pj (ξ)}j=1 .
d2 +p

reformulated into a deterministic optimization problem: Here, Np = p is the total number of basis functions
for ξ ∈ Rd2 bounded by order p. Then we set Ψ1 (ξ) = 1
min Eξ [f (x, ξ)] Np
x∈X
q and generate the orthonormal polynomials {Ψj (ξ)}j=2 in the
s.t. Eξ [yi (x, ξ)] + κi varξ [yi (x, ξ)] ≤ ui , ∀ i ∈ [n]. correlated parameter space recursively by
j−1
(17) X
Ψ̂j (ξ) = pj (ξ) − E[pj (ξ)Ψi (ξ)]Ψi (ξ),
i=1
IV. A LGORITHM AND I MPLEMENTATION D ETAILS
Ψ̂j (ξ)
We cannot solve problem (17) directly because we do Ψj (ξ) = q , j = 2, . . . , Np . (23)
not know the mean values and variances for the black-box E[Ψ̂2j (ξ)]
functions {yi (x, ξ)}ni=1 and f (x, ξ). A direct approach is to p N
apply a Monte Carlo method to estimate the mean values and These basis functions {Ψj (ξ)}j=1 can be re-ordered into the
variances for every iterate x. However, this is not affordable graded lexicographic order {Ψβ (ξ)}p|β|=0 .
because of the large number of numerical simulations.
In this section, we build the surrogate model for f (x, ξ) B. How to Build the Surrogate Models?
and {yi (x, ξ)}ni=1 by using generalized polynomial chaos [47] By a projection approach, the coefficient ciα,β for the basis
and our recent developed uncertainty quantification solver function can be computed by
[28], [29]. Once the surrogate models are constructed, we can ciα,β = Ex,ξ [yi (x, ξ)Φα (x)Ψβ (ξ)]. (24)
perform deterministic optimization. The main task is to build
the orthogonal basis functions Φα (x), Ψβ (ξ) and compute the The above integration can be well computed given a suitable
coefficients ciα,β and hα,β such that set of quadrature points and weights {xk , ξk , wk }M
k=1 :

p M
X
ciα,β ≈ yi (xk , ξk )Φα (xk )Ψβ (ξk )wk . (25)
X
yi (x, ξ) ≈ ciα,β Φα (x)Ψβ (ξ), (18)
k=1
|α|+|β|=0
5

We need to design a proper quadrature rule. The main chal- Algorithm 1: Our Proposed Chance-Constrained Yield-
lenge here is that x is a independent vector but ξ describes aware Optimization Solver
non-Gaussian correlated uncertainties. Input: The range of the design variables x, PDF of the
In this paper, we propose a three-stage optimization method non-Gaussian correlated random parameters ρ(ξ),
to compute the quadrature points and weights: the polynomial order p, the upper bounds of
M1
• Firstly, we compute the quadrature rule {xl , vl }l=1 for performance metrics {ui }ni=1 , and the chance
the independent design variables x. constraint thresholds {i }ni=1 .
• Secondly, we employ the optimization approach pro- 1. Construct the basis functions Φα (x) and Ψβ (ξ) based on
posed in [28], [29] to calculate the quadrature points (22) and (23) independently.
and weights {ξl , ul }M 2
l=1 for the non-Gaussian correlated 2. Initialize the quadrature points for design variables
parameters ξ. {xl , vl }M 1
l=1 by (26), and quadrature points for stochastic
• Finally, we use their tensor products (M1 M2 points) parameters {ξl , ul }M 2
l=1 by the optimization problem (27),
as an initialization and recall the optimization approach respectively. Then co-optimize the quadrature rule to
proposed in [28], [29] for the coupled space of x and obtain {xk , ξk , wk }M k=1 by (28).
ξ to compute M ≤ M1 M2 joint quadrature points 3. Call the simulator to compute f (xk , ξk ), yi (xk , ξk ) for
{xk , ξk , wk }M
k=1 . all i = 1, . . . , n and k = 1, . . . , M .
The details are described below. 4. Build the coefficients hα,β and ciα,β by equation (25).
1) Initial Quadrature Points for the design variable: 5. Set up the optimization problem (31), and then solve it
One could employ the sparse grid approach [48], [49] to via the global polynomial optimization solver [50].
compute the quadrature samples and weights for the inde- Output: The optimized design variable x∗
pendent uniform-distribution variables x ∈ Rd1 . However, the
quadrature weights from a sparse grid method can be negative,
and the number of quadrature points is not small enough.
Remark: We can also solve (28) directly to obtain the
Therefore, after obtaining the sparse-grid quadrature rule, we
optimized quadrature points. However, (28) is a non-convex
propose to refine the quadrature rule by our optimization solver
optimization problem and is hard to optimize in general. The
N2p
!2
X M1
X subproblems (26) and (27) help to provide a good initial guess
min E[Φj (x)] − Φj (xl )vl . (26) for the joint optimization.
a≤xl ≤b,vl ≥0
j=1 l=1 For all optimization subproblems (26), (27), and (28), we
Specifically, we use the quadrature points and weights obtained use the block coordinate-descent optimization method de-
from a sparse-grid rule as an initialization to solve the above scribed in [29] to compute the quadrature points and weights
optimization problem. The weighted clustering method in [28], alternatively. The following theorem ensures high accuracy for
[29] is employed to reduce the number of quadrature points. our surrogate model considering the unavoidable numerical
2) Initial Quadrature Points for uncertainty parameter: optimization error and function approximation error.
For the non-Gaussian correlated parameters ξ, we adopt the
Theorem 1. [29] Assume that {xk , ξk , wk }M
k=1 are the
optimization-based quadrature rule in [28], [29]. Specifically,
numerical solution to (28).
we compute M2 quadrature points ξl and weights wl via
solving the following optimization problem 1) Suppose the objective function of (28) decays to zero. The
N2p
!2 required number of quadrature points is upper and lower
M2
X X bounded by
min E[Ψj (ξ)] − Ψj (ξl )ul . (27)
ξl ,ul ≥0
j=1 l=1
(d + p)! (d + 2p)!
3) Optimized Joint Quadrature Points: The tensor product Np = ≤ M ≤ N2p = ; (29)
p!d! 2p!d!
of the two sets of quadrature points {xl , vl }M M2
l=1 and {ξl , ul }l=1
1

result in M1 M2 simulation points in total, which may be still 2) For any smooth and square integral function y(ξ), the
unaffordable for large-scale photonic design problems. In order approximation error of its p-th order stochastic approxi-
to further reduce the simulation cost of building surrogate mation ỹ(ξ) satisfies
models, we propose an optimization model to compute the
joint quadrature rule for both the design variables x and the ky(x, ξ) − ỹ(x, ξ)k2 ≤ α1 δ1 + α2 δ2 . (30)
uncertain parameters ξ:
!2 Pp
N2p N2p −j1
X X M
X Here, ỹ(x, ξ) = |α|+|β|=0 cα,β Φα (x)Ψβ (ξ), δ1
min δ0j1 δ0j2 − Φj1 (xk )Ψj2 (ξk )wk . is the `1 -norm of the objective function of (28)
a≤xk ≤b
ξk ,wk ≥0 j1 =1 j2 =1 k=1 evaluated at its final numerical solution, δ2 is the
(28) distance of y(x, ξ) to the p-th order polynomial space,
α1 = Np LT , α2 = 1 + Np W , L = max ky(x, ξ)k2 , T =
Here δ0j1 δ0j2 = 1 if j1 = j2 = 0 and zero otherwise. Our maxj1 +j2 ,l1 +l2 =1,...,N2p kΦj1 (x)Ψj2 (ξ)Ψl1 (x)Ψl2 (ξ)k2
numerical experiments show that the total number of optimized are constants,.
quadrature points is M is significantly smaller than M1 M2 .
6

p
X
min hα,0 Φα (x)
x∈X
|α|=0
 2  2
p p−|β| p p
X X X X
s.t. κ2i  ciα,β Φα (x) ≤ ui − ciα,0 Φα (x) , ciα,0 Φα (x) ≤ ui , ∀ i ∈ [n]. (31)
|β|=1 |α|=0 |α|=0 |α|=0

(a) (b) (c)


1 1 1

e e
se as as
0.5 rea 0.5 cre 0.5 cr
e
dec de de

0 0 0

-0.5 -0.5 -0.5

-1 -1
-1 -1 -0.5 0 0.5 1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1

Fig. 2. The feasible set of the synthetic example with risk tolerance levels  ∈ [10−2 , 10−0.1 ] under different uncertainty distributions. (a): a positive-
correlated non-Gaussian distribution; (b): a Gaussian independent distribution; (c): a negative correlated non-Gaussian distribution. The domain between the
red lines are the deterministic feasible set x21 ± x2 ≤ 1, and the blue lines demonstrate the effects of chance constraints.

C. The Proposed Polynomial Optimization nomial optimization sub-problem (31) is solved by the global
optimization solver GloptiPoly 3 [50]. The yield is defined as
With the formula for the mean value (20) and the variance
(21), we obtain the following deterministic formula for the the number of ξj s.t. yi (x, ξj ) ≤ ui , ∀ i ∈ [n]
yield(x) = .
chance-constrained optimization: the total number of random parameters ξj
(34)
p
X We set all risk thresholds to , i.e.,  = i , ∀ i ∈ [n]. For
min hα,0 Φα (x)
x∈X the synthetic example, we will compare our method with the
|α|=0
v deterministic formulation (15). For the photonic IC examples,
2
we will compare our method with the state-of-the-art Bayesian
u 
uXp p−|β|
X
ciα,β Φα (x) yield optimization method [10]. We summarize the key idea
u
s.t. κi t 
|β|=1 |α|=0 of the Bayesian yield optimization in Appendix B.
p
X
+ ciα,0 Φα (x) ≤ ui , ∀ i ∈ [n]. (32) A. Synthetic Example
|α|=0
Firstly, we consider a synthetic example with two design
However, the constraints are non-smooth because of the variables and two non-Gaussian correlated random parameters.
square-root terms. This non-smooth optimization is hard to The design variables x admits a uniform distribution U[−1, 1]2
solve because it may not admit a gradient at some points [51]. and the uncertain parameter ξ follows a Gaussian mixture
Instead, we use the equivalent smooth polynomial formula: distribution. We define the yield criterion as (x1 +ξ1 )2 ±(x2 +
ξ2 ) ≤ 1 and our goal is to maximize Eξ [3(x1 +ξ1 )+(x2 +ξ2 )].
κ2i varξ [yi (x, ξ)] ≤ (ui − Eξ [yi (x, ξ)])2 . (33) We formulate the yield into chance constraints and derive the
following problem
Consequently, (17) can be reduced to a deterministic and
smooth optimization problem of x in (31). max Eξ [3(x1 + ξ1 ) − (x2 + ξ2 )]
x
Noting that both the objective function and the constraints s.t. Probξ (x1 + ξ1 )2 − (x2 + ξ2 ) ≤ 1 ≥ 1 − ,

of (31) are polynomials, we can obtain the global optimal
Probξ (x1 + ξ1 )2 + (x2 + ξ2 ) ≤ 1 ≥ 1 − .

solution by using any polynomial solvers, such as those based (35)
on semi-definite relaxation [52], [53]. To illustrate the effects of different parameter distributions,
we study three probability density functions: the independent
V. N UMERICAL E XPERIMENTS distribution N (0, I), the non-Gaussian positive correlations

1 1 1 0.75
In this section, we verify our proposed approach by a 2 N (1, Σ) + 2 N (−1, Σ) with Σ = 0.75 1
, and
1 T
synthetic example and two photonic IC examples. The poly- the non-Gaussian negative correlations 2 N ([1, −1] , Σ) +
7

(a) (b) (c) (d)


1 0.15 1 0.15
0.3 0.1

0.5 0.1 0.25 0.5 0.1 0.08

0.2 0.06
x2

2
x2
0 0.05

2
0.05 0
0.15
0.04
0.1
-0.5 0 -0.5 0
0.02
0.05
-1 -0.05 -1 -0.05
-1 0 1 -0.05 0 0.05 0.1 -1 0 1 -0.05 0 0.05 0.1
x1 x1 1
1

Fig. 3. The quadrature points and weights in the synthetic experiment. (a) and (b): The initial 2-D quadrature points for the design variables x and uncertain
parameters ξ by solving (26) and (27), respectively. (c) and (d): The optimized quadrature points for the joint 4-D space of x and ξ by solving (28). Here
we project the optimized 4-D quadrature points to the 2-D sub-space of x and ξ, respectively. The quadrature weights are shown in colors.

TABLE I
T HE OPTIMAL SOLUTION FOR THE SYNTHETIC EXPERIMENT UNDER
DIFFERENCE RISK THRESHOLD .

Algorithm x∗ Objective Yield (%)


Proposed ( = 0.01) 0.8630 -0.1172 2.4717 100
Proposed ( = 0.05) 0.9379 -0.0522 2.7616 100
Proposed ( = 0.10) 0.9587 -0.0402 2.8360 99.42
Proposed ( = 0.15) 0.9689 -0.0351 2.8717 93.84
Proposed ( =0.20) 0.9751 -0.0293 2.8959 87.49
(15) 1.0000 0 3.0000 41.34

 
1 T 1 −0.75
2 N ([−1, 1] , Σ) with Σ =
−0.75 1
. The feasi- Fig. 4. An optical band-pass filter with three microrings coupled in series.
ble sets under three probability density distributions are shown
in Fig. 2. The comparison clearly shows that the effects of
different uncertainties. For all three density functions, the bandwidth and extinction ratio [54]. A broad and flat passband
feasible regions are reduced when the risk level  decreases. with a high extinction ratio can be achieved by optimizing
Next we take the non-Gaussian positive correlated distribu- the coupling strengths between the microrings [54]. In this
tion as an example to compute the optimal solution of (35). example, we employ silicon as the waveguide material and
We first build the surrogate models for both the objective and assume the effective refractive index to be neff = 2.44 and the
constraints by the second-order polynomial basis functions. effective group index to be ng = 4.19 near the wavelength of
The optimized quadrature points {xl , vl }6l=1 for the design 1.55 µm. The design variables are the coupling coefficients
variables by (26) and {ξl , ul }6l=1 for the random parameter x = [K1 , K2 , K3 , K4 ] that are to be optimized within the
by (27) are shown in Fig. 3 (a) and (b), respectively. Directly interval of [0.3, 0.6]. The random variables are the small
tensorizing the two sets of quadrature points generates 36 deviations of the coupling coefficients, ξ = δx, and we assume
samples. We further solve (28) to reduce them to M = 19 that ξ is non-Gaussian correlated and has a variance of 0.006.
optimized samples and weights. According to Theorem 1, We mainly focus on three metrics of the microring filter:
the number of quadrature samples for d = 4, p = 2 should the 3dB bandwidth (BW, in GHz), the extinction ratio (RE,
be in the range [15, 70]. Our optimization algorithm obtains in dB) of the transmission at the drop port, and the roughness
M = 19, which is close to the theoretical lower bound. (σpass , in dB) of the passband that takes a standard deviation
We further show the results for different risk tolerance levels of the passband and measures the roughness of the passband.
 in Table I. A smaller  results in a smaller feasible domain (as The yield-aware optimization problem of the microring filter
shown in Fig. 2), and generates a higher yield with a smaller design can be formulated as:
objective value. Compared with the solution x̃ = [1, 0]T from
solving (15), our method can achieve a significantly higher
yield: our optimized yield is above 87% while solving (15) max Eξ [BW(x, ξ)]
x∈X
only leads to a yield of 41.34%.
s.t. Probξ (RE(x, ξ) ≤ RE0 ) ≥ 1 − ,
Probξ (σpass (x, ξ) ≤ σ0 ) ≥ 1 − , (36)
B. Microring Band-Pass Filter
We continue to consider the design of an optical band-pass where the yield is defined through the extinction ratio and the
filter consisting of three identical silicon microrings coupled roughness of the passband and is specified to be above the
in series, as shown in Fig. 4. In designing such a broadband yield level 1 − . In our simulation, the threshold extinction
optical filter, the coupling coefficients play an important role ratio (RE0 ) and the roughness of the passband (σ0 ) are -25dB
in determining some key performance metrics, such as the and 0.5dB, respectively.
8

(a) (b) (c)


0.2 1 20
MC
0.15 15 Proposed
PDF

PDF

PDF
0.1 0.5 10

0.05 5

0 0 0
110 115 120 -28 -27 -26 -25 0.2 0.3 0.4 0.5 0.6
Bandwidth (GHz) Extinction ratio (dB) Roughness (dB)

Fig. 5. The probability density functions (PDF) of the bandwidth, extinction ratio and roughness for the ring filter example at the optimal solution x∗ =
[0.5582, 0.4208, 0.3000, 0.6000] by our proposed optimization method with  = 0.05. Our surrogate model uses only 64 simulations, and Monte Carlo
(MC) uses 103 simulation points.

Fig. 6. The transmission curves of the microring filter at different design choices. The grey lines show the uncertainties caused by the process variations.
The orange and blue curves show the mean transmission rates at the drop port and the through port, respectively. Here RE, BW and σpass denote the mean
values of extinction ratio, bandwidth and roughness, respectively. (a) The transmission at x0 = [0.45, 0.45, 0.45, 0.45] without any optimization. It doesn’t
have a clear passband because σpass is too large. (b) The results after the Bayesian yield optimization; (c) The results obtained from our chance-constrained
optimization with  = 0.05.

TABLE II 0.25
O PTIMIZATION RESULTS FOR THE PHOTONIC RING FILTER . BYO
proposed
0.2
Algorithm Simulations Eξ [BW] (GHz) Yield(%)
BYO [10] 2020 112.3 99.8 0.15
PDF

Proposed ( = 0.03) 64 113.4 100


Proposed ( = 0.05) 64 115.6 99.8 0.1
Proposed ( = 0.07) 64 117.2 99.5
Proposed ( = 0.10) 64 118.4 98.1 0.05

0
108 110 112 114 116 118 120
Bandwidth (GHz)
We first build the second-order polynomial surrogate model
by our proposed Algorithm 1. We only need 17 initial quadra- Fig. 7. The optimized bandwidth probability density distribution of the
ture points for the variable x by solving (26), 16 quadrature microring filter. Our chance-constrained optimization obtain an expected
points for the parameters ξ by solving (27), and 64 quadrature value of 115.6 GHz while the Bayesian yield optimization only produces
an expected value of 112.3 GHz.
points for the joint optimization of x and ξ by solving (28).
Fig. 5 shows that our surrogate model can well approximate
the probabilistic distributions of the performance metrics with
the comparison of 103 Monte Carlo simulations, although our
method only needs 64 simulation samples for this example.
We summarize the results of our proposed method with
different choices of  and the results obtained by Bayesian Fig. 8. The scheme of the third-order MZ filter.
yield optimization (BYO) in Table II. It is shown that with
lower risk tolerance level  our proposed method can achieve
higher yield while the expected value may decrease, because a optimal solution from the polynomial optimization solver, with
lower risk level  demands higher yield and shrinks the feasible  = 0.05 our proposed method can achieve the same level
region. Because our proposed method can return the global of yield at 99.8% but a higher value of Eξ [BW] than BYO.
9

(a) (b) (c)


0.4 1.5 6
MC
Proposed
0.3
1 4

PDF

PDF
PDF

0.2
0.5 2
0.1

0 0 0
182 184 186 188 190 -10 -9 -8 -7 0.2 0.4 0.6 0.8 1
Bandwidth (GHz) Crosstalk (dB) Attenuation (dB)

Fig. 9. The probability density functions (PDF) for the bandwidth, crosstalk, and attenuation of the MZ filter at our optimized design variable x∗ =
[0.300, 0.5036, 0.300]. Here our surrogate model uses only 36 simulations and Monte Carlo (MC) uses 103 simulations.

Fig. 10. The transmission curves of the bandwidth for the MZ filter at different design variables. The grey lines show the uncertainties caused by the process
variations. The orange curve and blue curve show the transmission rates at the drop port and through port, respectively. The mean values of the bandwidth,
crosstalk and attenuation are denoted as BW, CT and α, respectively. (a) The result of the initial design x0 = [0.45, 0.45, 0.45]; (b) The result after Bayesian
yield optimization; (c) The design obtained from the proposed chance-constrained yield optimization.

Our proposed method requires only 64 simulation points and and the attenuation. The yield-aware optimization problem of
achieves 32× reduction in terms of simulation cost compared the MZ filter design can be formulated as:
to 2020 simulation points required by BYO. In Fig. 6, we
compare the transmission lines before and after the yield-aware max Eξ [BW(x, ξ)]
x
optimization. After the optimization, we can achieve very high s.t. Probξ (CT(K, ξ) ≤ CT0 ) ≥ 1 − ,
bandwidth with a very smooth passband and a low extinction
Probξ (α(x, ξ) ≤ α0 ) ≥ 1 − , (37)
ratio. Moreover, our proposed method can have a higher band-
width and a smoother passband compared to BYO. In Fig. 7, where the yield risk level is . In our simulation, the threshold
we further plot the probability density of the bandwidth at the crosstalk (CT0 ) and attenuation (α0 ) are -4 dB and 1 dB,
optimal design by our chance-constrained optimization with respectively.
 = 0.05 and by the Bayesian yield optimization, respectively. We first build three second-order polynomial surrogate
It clearly shows that our proposed method can increase the models for the bandwidth, crosstalk, and the attenuation of
bandwidth while achieving the same yield. ctionMach-Zehnder the transmission at the drop port of the MZ filter by our
Bandpass Filter proposed Algorithm 1. We used 11 quadrature points in the
We apply the same framework to optimize a third-order design variable x, 10 quadrature points for the uncertainty
Mach-Zehnder (MZ) filter which consists of three port cou- parameter ξ, and 36 quadrature points for the joint space after
pling and two MZ arms, as shown in Fig. 8, where the the co-optimization. Fig. 9 shows that our surrogate models
coupling coefficients between the MZ arms play the most constructed with 36 quadrature points can well approximate
important role in the designing of the MZ bandpass filter. the density functions of all three performance metrics com-
The design variables x = [K1 , K2 , K3 ] are to be optimized pared with Monte Carlo with 103 sample points.
within the interval of [0.3, 0.6]. The non-Gaussian correlated We also compare our proposed method and BYO in Ta-
random variables ξ = δx are the small deviations of the design ble III. Similar to the result in Table II, lower risk tolerance
variable x and has a variance of 0.006. We consider three results in higher yield and lowers down the expected value
metrics of the MZ filter: the 3dB bandwidth (BW, in GHz), of bandwidth. Our method requires 56× fewer simulation
the crosstalk (CT, in dB) and the attenuation (α, in dB) of the points than BYO, which is a great advantage for design
peak transmission. The yield is defined through the crosstalk cases with the time-consuming simulations. For  = 0.05, the
10

0.4 process variations. How to handle non-smoothness in this


BYO
proposed optimization framework is a critical issue.
0.3 • High Dimensionality. Large-scale photonic ICs may have
a huge number of design variables and process variation
PDF

0.2 parameters. This brings new challenges to the surrogate


modeling and the resulting polynomial optimization in
0.1
our framework.
0
175 180 185 190
Bandwidth (GHz)
A PPENDIX A
D ETAILED D ERIVATION OF E QUATION (6)
Fig. 11. The optimized bandwidth of the MZ filter by the Bayesian yield We show that for u > Eξ [y(x, ξ)] the following determin-
optimization and our proposed method, respectively. The expectation bandwith
of the Bayesian yield optimization is 175.4 GHz while our proposed method istic constraint
with  = 0.05 can get 186.4 GHz. p q
Eξ [y(x, ξ)] + (1 − )/ varξ [y(x, ξ)] ≤ u
TABLE III is a sufficient but not necessary condition for the following
O PTIMIZATION RESULT FOR THE MZ FILTER .
probability constraint:
Algorithm Simulations Eξ [BW] (GHz) Yield(%)
BYO [10] 2020 175.4 100 Probξ (y(x, ξ) ≤ u) ≥ 1 − .
Proposed ( = 0.03) 36 183.8 100
Proposed ( = 0.05) 36 186.4 100 In other words, we want to show that each feasible point of
Proposed ( = 0.07) 36 187.9 100 (16) is a feasible point of the chance constraint (13).
Proposed ( = 0.10) 36 189.6 97.3
Denote the random variable as X = y(x, ξ). Cantelli’s
inequality [55] states that for any random variable X with
a mean value E[X] = Eξ [y(x, ξ)] and variance σ 2 =
optimized nominal design is x∗ = [0.3000, 0.5036, 0.3000] varξ [y(x, ξ)], it holds that the probability of a single tail can
and its expected bandwidth is 186.4GHz. In Fig. 10, we be bounded as follows:
compare the transmission lines before and after the yield-
aware optimization. Our proposed method can have a higher σ2
Prob(X − E[X] ≤ λ) ≥ 1 − if λ > 0. (38)
bandwidth and a smoother passband compared to Bayesian σ 2 + λ2
yield optimization and the initial design. Fig. 11 further Therefore, for any constant u ≥ E[X] we have
shows the probability density of the optimized bandwidth by
our chance-constrained optimization and the Bayesian yield Prob(X ≤ u) = Prob(X − E[X] ≤ u − E[X])
optimization, respectively. It clearly shows that our proposed σ2
≥1− 2 .
method produces higher bandwidth. σ + (u − E[X])2
For any , a sufficient condition for Prob(X ≤ u) ≥ 1 −  is
σ2
VI. C ONCLUSIONS AND R EMARKS 1 − σ2 +(u−E[X]) 2 ≥ 1 − , i.e.,

This paper has presented a data-efficient framework for


p
E[X] + (1 − )/σ ≤ u. (39)
the yield-aware optimization of the photonic ICs under non-
Gaussian correlated process variations. We have proposed Substituting X = y(x, ξ) in the above equation will arrive (6).
to reformulate the stochastic chance-constrained optimization The proof is completed.
into a deterministic polynomial optimization problem. Our
framework only requires simulation at a small number of A PPENDIX B
important points and admits a surrogate model for yield- BAYESIAN Y IELD O PTIMIZATION (BYO)
aware optimization. In the experiments by the microring filter Bayesian yield optimization (BYO) is a state-of-the-art
and the Mach Zehnder filter, we have demonstrated that our tool for the yield optimization of electronic devices and
optimization scheme can give high yield and high bandwidth. circuits [10]. This method approximates and optimizes the
Compared with Bayesian yield optimization, our method has posterior distribution of design variable under the condition
consumed much fewer simulation samples and produced better of “pass” events S. Here, S stands for the samples that
design performance while achieving the same yield. pass the yield criterion constraints in (1). Specifically, with
This work should be regarded as a presentation of some the Bayes’ theorem, Prob(S|x) = Prob(S)
Prob(x) Prob(x|S). In our
preliminary results in this direction. Many problems are worth problem setting, Prob(x) is a constant because we assume
further investigation in the future, for instance: that x follows a uniform distribution and Prob(S) should
• Non-Smoothness. Similar to generalized polynomial also be a constant without the dependence on the variable x.
chaos [47], the surrogate modeling techniques in [28], Therefore, we have Prob(S|x) ∝ Prob(x|S) and the original
[29] require the stochastic functions to be smooth. How- yield optimization problem (3) is equivalent to
ever, some performance metrics of a photonic IC may
be non-smooth with repsect to the design variables and xBY O = argmax Prob(x|S). (40)
x∈X
11

The paper [10] proposed an expectation-maximization frame- [13] X. Li, P. Gopalakrishnan, Y. Xu, and L. T. Pileggi, “Robust ana-
work to solve it. At the t-th iteration, the expectation step log/RF circuit design with projection-based performance modeling,”
IEEE Trans. CAD of Integr. Circ. Syst., vol. 26, no. 1, pp. 2–15, 2006.
approximates the probability by the kernel density estimation. [14] X. Li, J. Le, L. T. Pileggi et al., “Statistical performance modeling
Specifically, we generate N = 100 samples randomly and and optimization,” Foundations and Trends R in Electronic Design
call the simulator to compute the quantity of interests at those Automation, vol. 1, no. 4, pp. 331–480, 2007.
[15] X. Li, J. Sun, F. Xiao, and J.-S. Tian, “An efficient Bi-objective
samples. Then we choose M ≤ N “pass” samples to perform optimization framework for statistical chip-level yield analysis under
the kernel density estimation parameter variations,” Frontiers of Information Technology & Electronic
Engineering, vol. 17, no. 2, pp. 160–172, 2016.
M
1 X 1 1 [16] Z. Zhang, T. A. El-Moselhy, I. A. M. Elfadel, and L. Daniel, “Stochastic
Prob(x|S) ≈ √ exp (− (x − µi )T (x − µi )), testing method for transistor-level uncertainty quantification based on
M i=1 2πh 2h generalized polynomial chaos,” IEEE Trans. Computer-Aided Design
Integr. Circuits Syst., vol. 32, no. 10, pp. 1533–1545, Oct. 2013.
where {µi }Mi=1 ∈ S are design variable samples can pass [17] Z. Zhang, T. A. El-Moselhy, I. M. Elfadel, and L. Daniel, “Calculation
the yield constraints and h = 0.3 is a bandwidth parameter. of generalized polynomial-chaos basis functions and Gauss quadrature
rules in hierarchical uncertainty quantification,” IEEE Trans. Computer-
Afterward, the maximization step returns an updated design Aided Design of Integrated Circuits and Systems, vol. 33, no. 5, pp.
variable xBY O,t . We will call the simulator again at this 728–740, 2014.
design variable to record its objective value and “pass” status. [18] Z. Zhang, I. Osledets, X. Yang, G. E. Karniadakis, and L. Daniel,
“Enabling high-dimensional hierarchical uncertainty quantification by
We terminate the algorithm if the maximal iteration number ANOVA and tensor-train decomposition,” IEEE Trans. CAD of Inte-
20 is reached, or the residue of two consecutive iterations is grated Circuits and Systems, vol. 34, no. 1, pp. 63 – 76, Jan 2015.
below 10−6 . After the whole optimization process, we return [19] Z. Zhang, T.-W. Weng, and L. Daniel, “Big-data tensor recovery for
high-dimensional uncertainty quantification of process variations,” IEEE
the design variable that can pass the yield constraint with the Trans. Components, Packaging and Manufacturing Technology, vol. 7,
highest bandwidth no. 5, pp. 687–697, 2017.
[20] P. Manfredi, D. V. Ginste, D. De Zutter, and F. G. Canavero, “Un-
xBY O = arg max BW(x) s.t. pass(xBY O,t ) = 1. certainty assessment of lossy and dispersive lines in SPICE-type en-
x∈xBY O,t
vironments,” IEEE Trans. Components, Packaging and Manufacturing
Technology, vol. 3, no. 7, pp. 1252–1258, 2013.
R EFERENCES [21] S. Vrudhula, J. M. Wang, and P. Ghanta, “Hermite polynomial based
interconnect analysis in the presence of process variations,” IEEE Trans.
[1] S. C. Nicholes, M. L. Masanovic, B. Jevremovic, E. Lively, L. A. Computer-Aided Design of Integrated circuits and systems, vol. 25,
Coldren, and D. J. Blumenthal, “The world’s first InP 8× 8 monolithic no. 10, pp. 2001–2011, 2006.
tunable optical router (MOTOR) operating at 40 Gbps line rate per port,”
[22] K. Strunz and Q. Su, “Stochastic formulation of SPICE-type electronic
in Proc. Optical Fiber Communication, 2009, pp. 1–3.
circuit simulation with polynomial chaos,” ACM Trans. Modeling and
[2] M. Kato, R. Nagarajan, J. Pleumeekers, P. Evans, A. Chen, A. Mathur,
Computer Simulation, vol. 18, no. 4, p. 15, 2008.
A. Dentai, S. Hurtt, D. Lambert, P. Chavarkar et al., “40-channel
[23] M. R. Rufuie, E. Gad, M. Nakhla, and R. Achar, “Generalized Hermite
transmitter and receiver photonic integrated circuits operating at a per
polynomial chaos for variability analysis of macromodels embeddedin
channel data rate 12.5 Gbit/s,” in National Fiber Optic Engineers
nonlinear circuits,” IEEE Trans. Components, Packaging and Manufac-
Conference. Optical Society of America, 2007, p. JThA89.
turing Technology, vol. 4, no. 4, pp. 673–684, 2013.
[3] X. Chen, M. Mohamed, Z. Li, L. Shang, and A. R. Mickelson, “Process
variation in silicon photonic devices,” Applied optics, vol. 52, no. 31, [24] J. Tao, X. Zeng, W. Cai, Y. Su, D. Zhou, and C. Chiang, “Stochas-
pp. 7638–7647, 2013. tic sparse-grid collocation algorithm (SSCA) for periodic steady-state
[4] T. Lipka, J. Müller, and H. K. Trieu, “Systematic nonuniformity anal- analysis of nonlinear system with process variations,” in Proc. Asia and
ysis of amorphous silicon-on-insulator photonic microring resonators,” South Pacific Design Automation Conference, 2007, pp. 474–479.
Journal of Lightwave Technology, vol. 34, no. 13, pp. 3163–3170, 2016. [25] R. Shen, S. X.-D. Tan, J. Cui, W. Yu, Y. Cai, and G.-S. Chen,
[5] Z. Lu, J. Jhoja, J. Klein, X. Wang, A. Liu, J. Flueckiger, J. Pond, and “Variational capacitance extraction and modeling based on orthogonal
L. Chrostowski, “Performance prediction for silicon photonics integrated polynomial method,” IEEE transactions on very large scale integration
circuits with layout-dependent correlated manufacturing variability,” (VLSI) systems, vol. 18, no. 11, pp. 1556–1566, 2009.
Optics express, vol. 25, no. 9, pp. 9712–9733, 2017. [26] A. Waqas, D. Melati, Z. Mushtaq, and A. Melloni, “Uncertainty quantifi-
[6] J. Pond, J. Klein, J. Flückiger, X. Wang, Z. Lu, J. Jhoja, and L. Chros- cation and stochastic modelling of photonic device from experimental
towski, “Predicting the yield of photonic integrated circuits using statis- data through polynomial chaos expansion,” in Integrated Optics: De-
tical compact modeling,” in Integrated Optics: Physics and Simulations vices, Materials, and Technologies XXII, vol. 10535. International
III, vol. 10242, 2017, p. 102420S. Society for Optics and Photonics, 2018, p. 105351A.
[7] T. W. Weng, D. Melati, A. I. Melloni, L. Daniel et al., “Stochastic [27] A. Waqas, D. Melati, P. Manfredi, and A. Melloni, “Stochastic process
simulation and robust design optimization of integrated photonic filters,” design kits for photonic circuits based on polynomial chaos augmented
Nanophotonics, vol. 6, no. 1, pp. 299–308, 2017. macro-modelling,” Optics express, vol. 26, no. 5, pp. 5894–5907, 2018.
[8] T.-K. Yu, S.-M. Kang, J. Sacks, and W. J. Welch, “An efficient method [28] C. Cui, M. Gershman, and Z. Zhang, “Stochastic collocation with non-
for parametric yield optimization of MOS integrated circuits,” in Proc. Gaussian correlated parameters via a new quadrature rule,” in Proc. 27th
Intl. Conf. on Computer-Aided Design, 1989, pp. 190–193. IEEE Conference on Electrical Performance of Electronic Packaging
[9] Y. Li, H. Schneider, F. Schnabel, R. Thewes, and D. Schmitt-Landsiedel, and Systems, 2018, pp. 57–59.
“DRAM yield analysis and optimization by a statistical design ap- [29] C. Cui and Z. Zhang, “Stochastic collocation with non-Gaussian cor-
proach,” IEEE Transactions on Circuits and Systems I: Regular Papers, related process variations: Theory, algorithms and applications,” IEEE
vol. 58, no. 12, pp. 2906–2918, 2011. Trans. Components, Packaging and Manufacturing Technology, vol. 9,
[10] M. Wang, F. Yang, C. Yan, X. Zeng, and X. Hu, “Efficient Bayesian no. 7, pp. 1362–1375, July 2019.
yield optimization approach for analog and SRAM circuits,” in Proc. [30] T.-W. Weng, Z. Zhang, Z. Su, Y. Marzouk, A. Melloni, and L. Daniel,
Design Automation Conference, 2017, pp. 1–6. “Uncertainty quantification of silicon photonic devices with correlated
[11] M. Wang, W. Lv, F. Yang, C. Yan, W. Cai, D. Zhou, and X. Zeng, and non-Gaussian random parameters,” Optics Express, vol. 23, no. 4,
“Efficient yield optimization for analog and SRAM circuits via Gaus- pp. 4242–4254, 2015.
sian process regression and adaptive yield estimation,” IEEE Trans. [31] T.-W. Weng, D. Melati, A. Melloni, and L. Daniel, “Stochastic sim-
Computer-Aided Design of Integrated Circuits and Systems, vol. 37, ulation and robust design optimization of integrated photonic filters,”
no. 10, pp. 1929–1942, 2018. Nanophotonics, vol. 6, no. 1, pp. 299–308, 2017.
[12] Y. Xu, K.-L. Hsiung, X. Li, I. Nausieda, S. Boyd, and L. Pileggi, [32] P. Li, H. Arellano-Garcia, and G. Wozny, “Chance constrained program-
“OPERA: optimization with ellipsoidal uncertainty for robust analog IC ming approach to process optimization under uncertainty,” Computers
design,” in Proc. Design Automation Conference, 2005, pp. 632–637. & chemical engineering, vol. 32, no. 1-2, pp. 25–45, 2008.
12

[33] A. Mesbah, S. Streif, R. Findeisen, and R. D. Braatz, “Stochastic Chunfeng Cui received the Ph.D. degree in com-
nonlinear model predictive control with probabilistic constraints,” in putational mathematics from Chinese Academy of
American Control Conference, 2014, pp. 2413–2419. Sciences, Beijing, China, in 2016 with a specializa-
[34] L. Blackmore, M. Ono, A. Bektassov, and B. C. Williams, “A proba- tion in numerical optimization. From 2016 to 2017,
bilistic particle-control approximation of chance-constrained stochastic she was a Postdoctoral Fellow at City University of
predictive control,” IEEE Trans. Robotics, vol. 26, no. 3, pp. 502–517, Hong Kong, Hong Kong. In 2017 She joined the
2010. Department of Electrical and Computer Engineering
[35] H. Akhavan-Hejazi and H. Mohsenian-Rad, “Energy storage planning in at University of California Santa Barbara as a Post-
active distribution grids: A chance-constrained optimization with non- doctoral Scholar.
parametric probability functions,” IEEE Transactions on Smart Grid, Dr. Cui’s research activities are mainly focused
vol. 9, no. 3, pp. 1972–1985, 2018. on the areas of tensor computing, uncertainty quan-
[36] Z. Wang, C. Shen, F. Liu, X. Wu, C.-C. Liu, and F. Gao, “Chance- tification, machine learning, and their interface. She has been working on
constrained economic dispatch with non-Gaussian correlated wind power tensor data analysis by convex and non-convex optimization, high-dimensional
uncertainty,” IEEE Trans. Power Systems, vol. 32, no. 6, pp. 4880–4893, uncertainty quantification with non-Gaussian correlations for electronics and
2017. photonics IC, and theoretical structural analysis of deep learning. She is the
[37] A. Nemirovski and A. Shapiro, “Convex approximations of chance recipient of the 2019 Rising Stars in Computational and Data Sciences, 2019
constrained programs,” SIAM Journal on Optimization, vol. 17, no. 4, Rising Stars in Electrical Engineering and Computational Sciences, the 2018
pp. 969–996, 2006. Best Paper Award of IEEE Electrical Performance of Electronic Packaging
[38] P. Gupta and E. Papadopoulou, “Yield analysis and optimization,” in and Systems (EPEPS), and the Best Journal Paper Award of Scientia Sinica
Handbook of Algorithms for Physical Design Automation. Auerbach Mathematica.
Publications, 2008, pp. 788–807.
[39] G. Calafiore, L. El Ghaoui et al., “Distributionally robust chance-
constrained linear programs with applications,” Techical Report, DAUIN,
Politecnico di Torino, Torino, Italy, 2005.
[40] W. Gautschi, “On generating orthogonal polynomials,” SIAM J. Sci. Stat.
Comput., vol. 3, no. 3, pp. 289–317, Sept. 1982.
[41] R. Ghanem and P. Spanos, Stochastic finite elements: a spectral ap-
proach. Springer-Verlag, 1991.
[42] D. Xiu and J. S. Hesthaven, “High-order collocation methods for
differential equations with random inputs,” SIAM J. Sci. Comp., vol. 27, Kaikai Liu received the B.S. degree in physics
no. 3, pp. 1118–1139, Mar 2005. in 2018 from Huazhong University of Science and
[43] X. Li, “Finding deterministic solution from under-determined equation: Technology, Wuhan, China. In 2018 he joined the
large-scale performance variability modeling of analog/RF circuits,” Department of Electrical and Computer Engineering
IEEE Transactions on Computer-Aided Design of Integrated Circuits at University of California Santa Barbara as a Ph.D.
and Systems, vol. 29, no. 11, pp. 1661–1668, 2010. student.
[44] Z. Zhang, X. Yang, G. Marucci, P. Maffezzoni, I. M. Elfadel, G. Karni- Kaikai’s research focuses on developing the novel
adakis, and L. Daniel, “Stochastic testing simulator for integrated circuits data-driving algorithms for photonic integrated cir-
and MEMS: Hierarchical and sparse techniques,” in Proc. IEEE Custom cuits design. He has been working on the tensorized
Integrated Circuits Conf. San Jose, CA, Sept. 2014, pp. 1–8. Bayesian optimization and the chance-constraint op-
[45] C. Cui and Z. Zhang, “High-dimensional uncertainty quantification timization method.
of electronic and photonic IC with non-Gaussian correlated process
variations,” IEEE Trans. Computer Aided Design of Integrated Circuits
and Systems, June 2019, doi: 10.1109/TCAD.2019.2925340.
[46] ——, “Uncertainty quantification of electronic and photonic ICs
with non-Gaussian correlated process variations,” in Proc. Intl. Conf.
Computer-Aided Design, 2018, pp. 1–8.
[47] D. Xiu and G. E. Karniadakis, “The Wiener–Askey polynomial chaos for
stochastic differential ezquations,” SIAM journal on scientific computing,
vol. 24, no. 2, pp. 619–644, 2002.
[48] V. Barthelmann, E. Novak, and K. Ritter, “High dimensional polynomial Zheng Zhang (M’15) received his Ph.D degree in
interpolation on sparse grids,” Adv. Comput. Math., vol. 12, no. 4, pp. Electrical Engineering and Computer Science from
273–288, Mar. 2000. the Massachusetts Institute of Technology (MIT),
[49] T. Gerstner and M. Griebel, “Numerical integration using sparse grids,” Cambridge, MA, in 2015. He is an Assistant Profes-
Numer. Algor., vol. 18, pp. 209–232, Mar. 1998. sor of Electrical and Computer Engineering with the
[50] D. Henrion, J.-B. Lasserre, and J. Löfberg, “GloptiPoly 3: moments, University of California at Santa Barbara (UCSB),
optimization and semidefinite programming,” Optimization Methods & CA. His research interests include uncertainty quan-
Software, vol. 24, no. 4-5, pp. 761–779, 2009. tification with applications to the design automa-
[51] F. H. Clarke, Optimization and nonsmooth analysis. Siam, 1990, vol. 5. tion of multi-domain systems (e.g., nano-scale elec-
[52] H. Waki, S. Kim, M. Kojima, and M. Muramatsu, “Sums of squares and tronics, integrated photonics, and autonomous sys-
semidefinite program relaxations for polynomial optimization problems tems), and tensor computational methods for high-
with structured sparsity,” SIAM Journal on Optimization, vol. 17, no. 1, dimensional data analytics.
pp. 218–242, 2006. Dr. Zhang received the Best Paper Award of IEEE Transactions on
[53] J. Nie, “An exact Jacobian SDP relaxation for polynomial optimization,” Computer-Aided Design of Integrated Circuits and Systems in 2014, the
Mathematical Programming, vol. 137, no. 1-2, pp. 225–255, 2013. Best Paper Award of IEEE Transactions on Components, Packaging and
[54] R. Orta, P. Savi, R. Tascone, and D. Trinchero, “Synthesis of multiple- Manufacturing Technology in 2018, two Best Paper Awards (IEEE EPEPS
ring-resonator filters for optical systems,” IEEE Photonics Technology 2018 and IEEE SPI 2016) and three additional Best Paper Nominations (CICC
Letters, vol. 7, no. 12, pp. 1447–1449, 1995. 2014, ICCAD 2011 and ASP-DAC 2011) at international conferences. His
[55] F. P. Cantelli, “Sui confini della probabilita,” in Atti del Congresso Ph.D. dissertation was recognized by the ACM SIGDA Outstanding Ph.D.
Internazionale dei Matematici: Bologna del 3 al 10 de settembre di Dissertation Award in Electronic Design Automation in 2016, and by the
1928, 1929, pp. 47–60. Doctoral Dissertation Seminar Award (i.e., Best Thesis Award) from the
Microsystems Technology Laboratory of MIT in 2015. He was a recipient
of the Li Ka-Shing Prize from the University of Hong Kong in 2011, and the
NSF CAREER Award in 2019.

You might also like