# PoBO: A Polynomial Bounding Method for Chance-Constrained Yield-Aware
Optimization of Photonic ICs

Zichang He and Zheng Zhang, Member, IEEE This work was partly supported by NSF
grants #1763699 and #1846476. Zichang He and Zheng Zhang are with Department
of Electrical and Computer Engineering, University of California, Santa
Barbara, CA 93106, USA (e-mails<EMAIL_ADDRESS>zhengzhang@ece.ucsb.edu).

###### Abstract

Conventional yield optimization algorithms try to maximize the success rate of
a circuit under process variations. These methods often obtain a high yield
but reach a design performance that is far from the optimal value. This paper
investigates an alternative yield-aware optimization for photonic ICs: we will
optimize the circuit design performance while ensuring a high yield
requirement. This problem was recently formulated as a chance-constrained
optimization, and the chance constraint was converted to a stronger constraint
with statistical moments. Such a conversion reduces the feasible set and
sometimes leads to an over-conservative design. To address this fundamental
challenge, this paper proposes a carefully designed polynomial function,
called optimal polynomial kinship function, to bound the chance constraint
more accurately. We modify existing kinship functions via relaxing the
independence and convexity requirements, which fits our more general
uncertainty modeling and tightens the bounding functions. The proposed method
enables a global optimum search for the design variables via polynomial
optimization. We validate this method with a synthetic function and two
photonic IC design benchmarks, showing that our method can obtain better
design performance while meeting a pre-specified yield requirement. Many other
advanced problems of yield-aware optimization and more general safety-critical
design/control can be solved based on this work in the future.

###### Index Terms:

Photonic design automation; yield-aware optimization; chance constrained
programming; uncertainty quantification; process variation.

##  I Introduction

The increasing process variations have resulted in significant performance
degradation and yield loss in semiconductor chip design and fabrications [1,
2]. Compared with electronic ICs, photonic ICs are more sensitive to process
variations (e.g., geometric uncertainties) due to their large device
dimensions compared with the small wavelength. Therefore, yield modeling and
optimization for photonic ICs is highly desired [3, 4].

The yield optimization and yield-aware robust design have been investigated in
the electronic design automation community for a long time and have been paid
increasing attention in the photonic design automation [5, 6, 7, 8]. Typical
yield-aware design techniques include geometric approaches [9], geostatistics-
motivated performance modeling [10], yield-aware Pareto surface [11], yield-
driven iterative robust design [12], computational intelligence assisted
approaches [13], corner-based method [14], Bayesian yield optimization [15]
and so forth. Yield estimation is the key component in a yield
optimizer.Advanced yield estimators can be generally classified as Monte-
Carlo-based [16, 17] and non-Monte-Carlo-based [18, 19, 20] methods. Among the
non-Monte-Carlo ones, surrogate modeling aims to approximate some circuit
behaviors under variations to speed up the sampling and simulation process
[21, 22, 23]. Typical surrogate models include posynomial models [24], linear
quadratic models [25], support vector machine [26, 27], Gaussian process [28,
29], sparse polynomial [30], generalized polynomial chaos expansions [31, 32],
and some variants [33]. Focusing on generalized polynomial chaos expansion,
advanced techniques have also been developed to handle high-dimensional [34,
35, 36, 37, 38], mixed-integer [39] or non-Gaussian correlated [40] process
variations. The polynomial chaos expansion based modeling and optimization has
been widely used in both electronics [41, 42, 43] and photonic IC design [44,
45, 46, 47].

While most existing yield optimization approaches try to maximize the yield of
a circuit, the obtained design performance (e.g., signal gain, power
dissipation) may be far from the achievable optimal solution. Recently, an
alternative approach was proposed in [48] to achieve both excellent yield and
design performance. Instead of simply maximizing the yield, the work [48]
optimizes a design performance metric while enforcing a high yield
requirement. Specifically, the yield requirement is formulated as some chance
constraints [49], which are further transformed to tractable constraints of
the first and second statistical moments. The chance constrained modeling
itself has been widely used in many engineering fields [50, 51, 52]. The
moment bounding method offers a provably sufficient condition of the chance
constraint. However, the bounding gap may be too large in many applications
[53]. The resulting overly-reduced feasible region may lead to an over-
conservative design.

Paper contributions. This paper proposes a novel Polynomial Bounding method
for chance-constrained yield-aware Optimization (PoBO) under truncated non-
Gaussian correlated variations. Leveraging the recent uncertainty
quantification techniques [40, 48], PoBO employs and modifies the idea of
kinship functions [54] to approximate the original chance constraints with a
better polynomial bounding method. PoBO provides a less conservative design
than moment-based bounding methods [48] while ensuring a pre-specified yield
requirement. The specific contributions of this paper include:

  * •

A better bounding method of the chance constraints via optimal polynomial
kinship functions. Compared with existing work in the control community [54],
we avoid the assumption of the independence among random variables and the
convexity of kinship functions. The relaxation allows more general non-
Gaussian correlated uncertainty modeling and tightens the bounding functions.
Within a family of polynomial functions, the optimal polynomial kinship
functions can be efficiently solved via semidefinite programming. Our bounding
method preserves the polynomial formulation of the provided surrogate models.
It enables the advanced polynomial optimization solvers, which provide a
sequence of convex relaxations via semidefinite optimization and searches for
the global design.

  * •

Numerical implementation of the PoBO framework. Based on available uncertainty
quantification solvers, we implement PoBO efficiently based on some pre-
calculated optimal polynomial kinship functions and quadrature samples and
weights without requiring any additional circuit simulations.

  * •

Validations on a synthetic function and two photonic IC design examples. Our
method offers better design performance while meeting the pre-specified yield
requirements. This method requires a small number of circuit simulations due
to its compatibility with recent data-efficient uncertainty quantification
methods [36, 40].

While this work focuses on the fundamental theory, algorithms and their
validation on small-size photonic circuits, the proposed method can be
combined with sparse or low-rank surrogate modeling methods [34, 36] to handle
large-scale design cases with much more design variables and process
variations.

##  II Background

This section reviews chance-constrained yield-aware optimization and its
implementation via moment bounding [48].

###  II-A Chance-Constrained Yield-Aware Optimization

We denote the design variables by
$\mathbf{x}=[x_{1},x_{2},\ldots,x_{d_{1}}]^{T}\in\mathbf{X}$, and the process
variations by random parameters
$\boldsymbol{{\xi}}=[\xi_{1},\xi_{2},\ldots,\xi_{d_{2}}]^{T}\in\mathbf{\Xi}$.
Let $\\{y_{i}(\mathbf{x},\boldsymbol{{\xi}})\\}_{i=1}^{n}$ denote $n$
performance metrics that are considered in yield estimation,
$\\{u_{i}\\}_{i=1}^{n}$ denote their corresponding upper bounds specifying the
design requirements. An indicator function is defined as

$I(\mathbf{x},\boldsymbol{{\xi}})=\begin{cases}1,&y_{i}(\mathbf{x},\boldsymbol{{\xi}})\leq
u_{i},\forall i=[n];\\\ 0,&\text{otherwise}.\end{cases}$ (1)

Here $[n]=\\{1,2,\cdots,n\\}$. The yield at a certain design choice
$\mathbf{x}$ is defined as

$Y(\mathbf{x})=\mathbb{P}_{\boldsymbol{{\xi}}}(\mathbf{y}(\mathbf{x},\boldsymbol{{\xi}})\leq\mathbf{u})=\mathbb{E}_{\boldsymbol{{\xi}}}[I(\mathbf{x},\boldsymbol{{\xi}})].$
(2)

In conventional yield optimization, one often tries to achieve the best
possible yield. This often requires losing remarkable design performance
$f(\mathbf{x},\boldsymbol{{\xi}})$ in order to achieve a high yield.

Simply maximizing the yield may lead to an over-conservative design. As an
example shown in Fig. 1, one may lose lots of performance (from $2.2$ to
$1.4$) while just getting marginal yield improvement from $99\%$ to $100\%$.

Figure 1: The trade-off between yield and an objective design performance for
the example from V-A. Simply maximizing the yield can lead to over-
conservative performance.

In order to avoid an over-conservative design, a chance-constrained
optimization was proposed in [48]:

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}\quad$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[f(\mathbf{x},\boldsymbol{{\xi}})]$
(3a) $\displaystyle\mathop{\mathrm{s.t.}}\quad$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}(y_{i}(\mathbf{x,}\boldsymbol{{\xi}})\leq
u_{i})\geq 1-\epsilon_{i},\forall i=[n].$ (3b)

where $f(\mathbf{x},\boldsymbol{{\xi}})$ is the performance metric that we aim
to optimize, and $\epsilon_{i}\in[0,1]$ is a risk level to control the
probability of meeting each design constraint. Instead of simply maximizing
the yield, the chance-constrained optimization tries to achieve a good balance
between yield and performance: one can optimize the performance
$f(\mathbf{x},\boldsymbol{{\xi}})$ while ensuring a high yield. The circuit
yield can be controlled by $\epsilon_{i}$: reducing $\epsilon_{i}$ ensures a
lower failure rate and thus a higher yield.

The chance-constrained optimization (3) is generally hard to solve. Firstly,
the feasible set produced by the chance constraints is often non-convex and
hard to estimate. Secondly, it is also expensive to estimate the design
objective function $f(\mathbf{x},\boldsymbol{{\xi}})$ and design constraint
function $y_{i}(\mathbf{x},\boldsymbol{{\xi}})$ due to the lack of analytical
expressions. Fortunately, a moment bounding method was combined with
uncertainty quantification techniques in [48] to make the problem tractable.

###  II-B Moment Bounding Method for (3)

Employing the Chebyshev-Cantelli inequality, one can ensure the chance
constraint via a moment bounding technique [55, 48]. Specifically, with the
first and second-order statistical moments of the constraint function, one can
convert the probabilistic constraint in (3) to a deterministic one:

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[{f}(\mathbf{x},\boldsymbol{{\xi}})]$
(4) $\displaystyle\mathop{\mathrm{s.t.}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[{y}_{i}(\mathbf{x},\boldsymbol{{\xi}})]+\gamma_{\epsilon_{i}}\sqrt{\text{Var}_{\boldsymbol{{\xi}}}[{y}_{i}(\mathbf{x},\boldsymbol{{\xi}})]}\leq
u_{i},\forall i\in[n],$

where constant
$\gamma_{\epsilon_{i}}=\sqrt{\frac{1-\epsilon_{i}}{\epsilon_{i}}}$. When the
objective and constraint functions are described by certain surrogate models
such as generalized polynomial chaos [32, 40], one can easily extract their
mean and variances. This can greatly simplify the problem and reduce the
computational cost, as shown by the yield-aware optimization of photonic IC in
[48].

The moment constraint in (4) is a sufficient but unnecessary condition of the
original chance constraint in (3). Therefore, any feasible solution of (4)
should satisfy the probability constraint of (3). However, the feasible set
produced by a moment bounding can be much smaller than exact one [53]. This
usually leads to an over-conservative design solution. When the risk level is
very small, the moment bounding method may even produce an empty feasible set,
leading to an unsolvable problem (see Section V).

Figure 2: The feasible set of an original chance constraint, moment bounding
(the cross region of two orange lines), and the proposed polynomial bounding.

##  III Proposed Polynomial Bounding Method

To avoid the possible over-conservative bounding of the moment methods [48],
we propose to bound the chance constraint via a more accurate polynomial
method. Fig. 2 plots the feasible regions obtained by different bounding
methods for the synthetic example in Sec. V-A. For the same chance
constraints, the moment bounding method produces a feasible set that is much
smaller than the exact one, whereas our polynomial bounding method generates a
better approximation of the feasible set. Due to the more accurate
approximation of the feasible set, our proposed polynomial bounding method can
provide a less conservative design in yield-aware optimization.

Now we describe how to generate the polynomial bounds via kinship functions.

###  III-A Kinship Function

The kinship function was first proposed to construct a convex approximation of
an indicator function in [54]. We generalize the concepts of [54] with two
relaxations:

  * •

We do not require the convexity of a kinship function.

  * •

We do not require the random variables $\boldsymbol{{\xi}}$ to be mutually
independent. Instead, we assume that they are truncated non-Gaussian
correlated.

We slightly modify the definition of a kinship function.

######  Definition 1 (Kinship function).

A kinship function $\kappa(z):[-1,\infty)\rightarrow\mathbb{R}$ is a function
that satisfies the following constraints:

  * •

$\kappa(z)=1$ when $z=0$;

  * •

$\kappa(z)\geq 0$ for any $z\in[-1,+\infty)$;

  * •

$\kappa(z_{1})\geq\kappa(z_{2})$ for any $z_{1}\geq z_{2}$ in the range
$[-1,\infty)$.

Based on kinship functions, the following theorem offers an upper bound for
any probability of constraint violations.

######  Theorem 1 (Risk integral [54]).

Let $\kappa(\cdot)$ be a kinship function,
$g(\mathbf{x},\boldsymbol{{\xi}})\geq-1$, and $V_{\kappa}(\mathbf{x})$ be a
risk integral quantity:

$V_{\kappa}(\mathbf{x})\coloneqq\int\limits_{\mathbf{\Xi}}\kappa[g(\mathbf{x},\boldsymbol{{\xi}})]\mu(\boldsymbol{{\xi}})d\boldsymbol{{\xi}},$
(5)

then we have
$\mathbb{P}\\{\boldsymbol{{\xi}}\in\mathbf{\Xi}:g(\mathbf{x},\boldsymbol{{\xi}})>0\\}\leq
V_{\kappa}(\mathbf{x})$.

###### Proof.

According to Definition 1, $\kappa[g(\mathbf{x},\boldsymbol{{\xi}})]$ is
nonnegative in $[-1,\infty)$ and greater than 1 if
$g(\mathbf{x},\boldsymbol{{\xi}})\geq 0$. Therefore, for any probability
measure $\mu(\boldsymbol{{\xi}})$ on $\mathbf{\Xi}$, we have:

$\displaystyle V_{\kappa}(\mathbf{x})$
$\displaystyle\geq\int\limits_{\\{\boldsymbol{{\xi}}\in\mathbf{\Xi},g(\mathbf{x},\boldsymbol{{\xi}})>0\\}}\kappa[g(\mathbf{x},\boldsymbol{{\xi}})]\mu(\boldsymbol{{\xi}})d\boldsymbol{{\xi}}$
(6a)
$\displaystyle\geq\int\limits_{\\{\boldsymbol{{\xi}}\in\mathbf{\Xi}:g(\mathbf{x},\boldsymbol{{\xi}})>0\\}}\mu(\boldsymbol{{\xi}})d\boldsymbol{{\xi}}=\mathbb{P}\\{\boldsymbol{{\xi}}\in\mathbf{\Xi},g(\mathbf{x},\boldsymbol{{\xi}})>0\\}.$
(6b)

∎

There exist many possible choices of kinship functions. Next, we will show how
to pick some polynomial kinship functions. We consider the polynomial function
family because it is compatible with existing surrogate modeling techniques
[32, 36, 34] to facilitate yield-aware optimization.

Figure 3: Different kinship functions. Black curve: an indicator function.

###  III-B Optimal Polynomial Kinship Function

An optimal kinship function is defined as a kinship function that minimizes
its integral over $[-1,0]$ [54]:

$\kappa^{\star}(\cdot)\coloneqq\mathop{\mathrm{argmin}}\limits_{\kappa(\cdot)\in{\cal
K}}\int_{-1}^{0}\kappa(z)dz.$ (7)

Here ${\cal K}$ is the set of all possible kinship functions that satisfy the
constraints in Definition 1. Let $z=g(\mathbf{x},\boldsymbol{{\xi}})$, then
the above definition can be understood as minimizing the gap between the left-
and right-hand sides of (6a).

Now we consider choosing a kinship function from a family of order-$\rho$
polynomials ${\cal
K}_{\rho}\coloneqq\\{\kappa(z)|\kappa(z)=\sum_{i=0}^{\rho}\zeta_{i}z^{i}\\}$.
The optimal polynomial kinship function, denoted as $\kappa_{\rho}(\cdot)$,
can be constructed by solving the following optimization problem:

$\begin{array}[]{cc}\min\limits_{\zeta_{0},\ldots,\zeta_{\rho}}&\int\limits_{-1}^{0}\kappa_{\rho}(z)dz\\\
\mathop{\mathrm{s.t.}}&\kappa_{\rho}(z)=\zeta_{0}+\zeta_{1}z+\cdots+\zeta_{\rho}z^{\rho}\in{\cal
K}_{\rho},\\\ &\kappa_{\rho}(0)=1,\;\kappa_{\rho}(-1)=0,\\\
&\kappa_{\rho}^{\prime}(z)\geq 0,\;\forall z\geq-1.\end{array}$ (8)

The optimization problem can be reformulated as semidefinite programming and
the details are provided in Appendix A. The obtained optimal polynomial
kinship function $\kappa_{\rho}(\cdot)$ with different orders are visualized
in Fig. 3. Since we relax the convexity requirement, given the same polynomial
order, the proposed polynomial function is a tighter approximation to the
indicator function than a convex one [54] for $z\in[-1,0]$.

Based on the obtained optimal polynomial kinship function
$\kappa_{\rho}(\cdot)$, we bound the original chance constraint (3b) by
enforcing the upper bound of the failure probability below $\epsilon_{i}$:

$\mathbb{P}_{\boldsymbol{{\xi}}}(y_{i}(\mathbf{x,}\boldsymbol{{\xi}})>u_{i})\leq\\\
\int\limits_{\mathbf{\Xi}}\kappa_{\rho}(y_{i}(\mathbf{x},\boldsymbol{{\xi}})-u_{i})\mu(\boldsymbol{{\xi}})d\boldsymbol{{\xi}}\leq\epsilon_{i}.$
(9)

The order of $\kappa_{\rho}(\cdot)$ controls the upper bound of violation
probability. When $\rho$ is small, the bounding gap in (6a) is large. As
$\rho\rightarrow\infty$, the polynomial bounding leads to a worst-case robust
design optimization according to Theorem 3 in [54]. This is equivalent to
setting a risk level $\epsilon_{i}=0$, which leads to an extremely over-
conservative design. Fortunately, this is not a trouble in practice since we
do not use a very high-order polynomial due to the computational issues. In
practice, there exists an optimal order $\rho^{\star}$ for bounding the
violation probability most accurately. The optimal $\rho^{\star}$ is unknown
a-priori, but heuristically we find that setting $\rho\in[5,10]$ usually
offers an excellent bound.

The proposed bounding method can be extended to deal with joint chance
constraints by constructing a multivariate polynomial kinship function. We
plan to report the results in a future paper.

##  IV The PoBO Framework

Based on the proposed polynomial bounding for chance constraints, we further
present the novel PoBO method to achieve less conservative yield-aware
optimization.

###  IV-A Workflow of PoBO

Our PoBO framework has two weak assumptions on the design and random
variables:

  * •

The design variable $\mathbf{x}$ is box-bounded, i.e.,
$\mathbf{x}\in\mathbf{X}=[a,b]^{d_{1}}$. This is normally the case in circuit
optimization.

  * •

The process variations $\boldsymbol{{\xi}}$ are truncated and non-Gaussian
correlated with a joint probability density function
$\mu(\boldsymbol{{\xi}})$.

The 2nd assumption is not strong at all. Many practical process variations are
correlated and not guaranteed to be Gaussian. Additionally, the values of
almost all practical geometric or material parameters are bounded, although
some simplified unbounded distributions (e.g., Gaussian distributions) were
used in previous literature for the ease of implementation.

The overall flow of PoBO is summarized below.

  * •

Step 1: Surrogate modeling. We use the recent uncertainty quantification
solver [40] to construct polynomial surrogate models for the objective and
constraint functions, i.e.,
$f(\mathbf{x},\boldsymbol{{\xi}})\approx\hat{f}(\mathbf{x},\boldsymbol{{\xi}})$
and
$y_{i}(\mathbf{x},\boldsymbol{{\xi}})\approx\hat{y}_{i}(\mathbf{x},\boldsymbol{{\xi}}),\forall
i=[n]$.

  * •

Step 2: Bounding the chance constraints via the proposed optimal polynomial
kinship functions. This transforms a chance-constrained probabilistic
optimization problem to a tractable deterministic one with a high-quality
solution.

  * •

Step 3: Design optimization. We use a polynomial optimization solver, e.g.,
the semidefinite programming [56], to obtain a globally optimal solution.

The PoBO framework reformulates the original chance-constrained optimization
(3) to the following optimization:

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}\quad$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[\hat{f}(\mathbf{x},\boldsymbol{{\xi}})]$
(10a) $\displaystyle\mathop{\mathrm{s.t.}}\quad$ $\displaystyle
V_{\kappa}^{(i)}(\mathbf{x})=\int\limits_{\mathbf{\Xi}}\kappa_{\rho}(\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}}))\mu(\boldsymbol{{\xi}})d\boldsymbol{{\xi}}\leq\epsilon_{i},\;\forall
i\in[n].$ (10b)

Here
$\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})\coloneqq\hat{y}_{i}(\mathbf{x},\boldsymbol{{\xi}})-u_{i}$,
$\hat{f}(\mathbf{x},\boldsymbol{{\xi}})$ and
$\hat{y}_{i}(\mathbf{x},\boldsymbol{{\xi}})$ are the polynomial surrogate
models of $f(\mathbf{x},\boldsymbol{{\xi}})$ and
$y_{i}(\mathbf{x},\boldsymbol{{\xi}})$, respectively.

In the following subsections, we will describe the implementation details of
this PoBO framework.

###  IV-B Building Surrogate Models

We employ the advanced stochastic collocation method with non-Gaussian
correlated uncertainty [40, 48] in Step 1. This method approximates a smooth
stochastic function as the linear combination of some orthogonal and
normalized polynomial basis functions:

$f(\mathbf{x},\boldsymbol{{\xi}})\approx\hat{f}(\mathbf{x},\boldsymbol{{\xi}})=\sum_{\lvert\boldsymbol{\alpha}\rvert+\lvert\boldsymbol{\beta}\rvert=0}^{p}c_{\boldsymbol{\alpha},\boldsymbol{\beta}}\mathbf{\Phi}_{\boldsymbol{\alpha}}(\mathbf{x})\mathbf{\Psi}_{\boldsymbol{\beta}}(\boldsymbol{{\xi}}).$
(11)

Here $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are two index vectors,
$\mathbf{\Phi}_{\boldsymbol{\alpha}}(\mathbf{x})$ and
$\mathbf{\Psi}_{\boldsymbol{\beta}}(\boldsymbol{{\xi}})$ are two series of
orthogonal polynomial basis functions, and $p$ upper bounds the total order of
the product of two basis functions. The corresponding coefficients
$c_{\boldsymbol{\alpha},\boldsymbol{\beta}}$ are calculated via a projection
method using some optimization-based quadrature samples and weights of
$\mathbf{x}$ and $\boldsymbol{{\xi}}$ [40]. When the parameter dimensionality
is not high, this method only needs a small number of simulation samples to
produce a highly accurate surrogate model with a provable error bound.

###  IV-C Scaling the Yield Metrics
${\upsilon}_{i}(\mathbf{x},\boldsymbol{{\xi}})$

To bound the failure probability
$\mathbb{P}\\{\boldsymbol{{\xi}}\in\mathbf{\Xi}:\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})>0\\}$
via the optimal kinship function,
$\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})$ must be in the range
$[-1,\infty)$ according to Definition 1. Once
$\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})$ is lower bounded, we can always
scale it to meet this requirement.

Since $\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})$ is a polynomial function
in our problem setting, and both $\mathbf{x}$ and $\boldsymbol{{\xi}}$ are
assumed bounded, we can compute the minimum value
${\upsilon_{i}^{-}}\coloneqq{\min\limits_{\mathbf{x}\in\mathbf{X},\boldsymbol{{\xi}}\in\mathbf{\Xi}}\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})}$.
Then we change the lower bound of
$\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})$ to -1 as follows:

$\begin{matrix}\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})\longleftarrow-\frac{1}{\upsilon_{i}^{-}}\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}}),\;\forall
i\in[n].\end{matrix}$ (12)

The scaling factor $-\frac{1}{\upsilon_{i}^{-}}$ is positive as long as the
problem (10) is solvable. This is because the existence of
$\mathbf{x}\in\mathbf{X}$ and $\boldsymbol{{\xi}}\in\mathbf{\Xi}$ such that
$\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}})<0$ is the necessary condition to
satisfy the yield constraint. We can easily avoid $\upsilon_{i}^{-}=0$ by
adding a sufficiently small perturbation.

###  IV-D Calculating Risk Integral $V_{\kappa}^{(i)}(\mathbf{x})$

Figure 4: An example of optimized quadrature rule for a two-dimensional
truncated Gaussian mixture model with two components. The color bar of the
right figure represents the weights of all samples. The shown quadrature rule
satisfies the exact integration up to order 6.

In order to upper bound the probability of violating a design constraint, we
need to calculate the integration in (10b). This can be a challenging task for
a truncated non-Gaussian correlated random vector $\boldsymbol{{\xi}}$ since
classical numerical quadrature rules [57, 58] do not work for non-Gaussian
correlated variables. Fortunately, we can reuse the quadrature rule of
$\boldsymbol{{\xi}}$ as the by-product of building the surrogate models in
Sec. IV-B. Specifically, the quadrature points and weights
$\\{\boldsymbol{{\xi}}_{l},w_{l}\\}_{l=1}^{M}$ compute the exact integration
up to order $2q$, obtained by solving the following optimization problem [40]:

$\min_{\boldsymbol{{\xi}}_{l},w_{l}\geq
0}\;\sum_{\lvert\boldsymbol{\beta}\rvert=0}^{2q}{\left(\mathbb{E}_{\boldsymbol{{\xi}}}\left[\mathbf{\Psi}_{\boldsymbol{\beta}}(\boldsymbol{{\xi}})\right]-\sum_{l=1}^{M}\mathbf{\Psi}_{\boldsymbol{\beta}}(\boldsymbol{{\xi}}_{l})w_{l}\right)}^{2}.$
(13)

An example of the solved quadrature is shown in Fig. 4. The detailed accuracy
analysis and the bound of $M$ are provided in [40], which is omitted here.

1:Box-bounded design variable $\mathbf{x}\in\mathbf{X}$, truncated non-
Gaussian correlated variations $\boldsymbol{{\xi}}\in\mathbf{\Xi}$, risk
levels $\mathbf{\epsilon}$

2:Optimized design $\mathbf{x}^{\star}$

3:Formulate the chance-constrained problem (3).

4:Obtain surrogate models for the design objective function
$f(\mathbf{x},\boldsymbol{{\xi}})\approx\hat{f}(\mathbf{x},\boldsymbol{{\xi}})$
and constraint functions
$y_{i}(\mathbf{x},\boldsymbol{{\xi}})\approx\hat{y}_{i}(\mathbf{x},\boldsymbol{{\xi}}),\forall
i\in[n]$.

5:Scale the yield metrics ${\upsilon_{i}}(\mathbf{x},\boldsymbol{{\xi}})$ via
Eq. (12).

6:Compute the optimal polynomial kinship function $\kappa_{\rho}(\cdot)$ via
(8).

7:Calculate risk integral $V_{\kappa}^{(i)}(\mathbf{x})$ via the quadrature
rule obtained from (13).

8:Seek the optimal design of problem (14) via a global polynomial optimization
solver.

Algorithm 1 Flow of the proposed PoBO.

Theoretically, we need a quadrature rule to exactly calculate the integration
up to order $p\rho$ in (10b). The exact quadrature rule can be computed
offline via solving (13) with $q=\lceil\frac{\rho p}{2}\rceil$. In practice, a
low-order quadrature rule, like $q=p$, often offers sufficient numerical
accuracy. Therefore, we can directly use the quadrature rule used when
building surrogate model (11) to calculate the risk integral
$V_{\kappa}^{(i)}(\mathbf{x})$. Since the quadrature points will not be
simulated, it does not introduce any additional computational burdens. Based
on the quadrature rule, problem (10) can be converted to the following
deterministic constrained polynomial optimization (14):

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[\hat{f}(\mathbf{x},\boldsymbol{{\xi}})]$
(14) $\displaystyle\mathop{\mathrm{s.t.}}$
$\displaystyle\sum_{l=1}^{M}w_{l}\kappa_{\rho}(\upsilon_{i}(\mathbf{x},\boldsymbol{{\xi}}_{l}))\leq\epsilon_{i},\quad\forall
i\in[n].$

Note that the expectation value in the objective function can be easily
obtained since $\hat{f}(\mathbf{x},\boldsymbol{{\xi}})$ is a generalized
polynomial-chaos expansion [32, 40].

###  IV-E Algorithm Summary

We summarize PoBO in Alg. 1. Below are some remarks:

  * •

In line 4 of Alg. 1, the optimal polynomial kinship functions (8) can be
computed offline.

  * •

In line 5 of Alg. 1, the optimization-based quadrature rule in [40] can be
used to calculate the risk integral without any additional simulations. The
quadrature rule can be computed offline as well.

  * •

This method enables a global polynomial optimization solver to obtain the
optimal design of (14).

Since a high polynomial order can increase the computational burden, the
current PoBO framework becomes inefficient for high-dimensional design
problems. This challenge may be overcome in the future by using other
nonlinear optimization solvers, or a better polynomial optimization solver
(e.g. a sparse polynomial solver) that can exploit the sparse structure of the
polynomial surrogate.

##  V Numerical Results

In this section, we validate the proposed PoBO framework via the synthetic
example and two realistic photonic IC examples from [48]. The polynomial
optimization is solved via GloptiPoly 3 [56], which is a global optimization
solver based on hierarchical semidefinite programming. Our codes are
implemented in MATLAB and run on a computer with a 2.3 GHz CPU and 16 GB
memory.

TABLE I: Optimization Results for Synthetic Function

Risk level $\epsilon$ Method Objective $\Delta_{1}$ (%) $\Delta_{2}$ (%) Yield
(%) 0.01 Moment [48] N/A* N/A* N/A* N/A* Proposed 1.14 1.01 0.99 99.98 0.05
Moment [48] 1.88 5.25 5.26 99.98 Proposed 2.11 5.21 3.97 98.76 0.1 Moment [48]
2.19 10.47 10.98 99.36 Proposed 2.26 10.67 7.98 97.08 * The algorithm fails
with no feasible solution.

Figure 5: The schematic of a third-order Mach-Zehnder interferometer. Figure
6: The transmission curves of the MZI. The grey lines show the performance
uncertainties. The orange and blue curves show the transmission rates at the
drop and through ports, respectively. The mean values of the bandwidth,
crosstalk, and attenuation are denoted as BW, XT, and $\alpha$, respectively.
(a) The initial design: $\mathbf{x}$=[150, 150, 150]; (b) Design after
Bayesian yield optimization [15]: $\mathbf{x}$=[286.63, 170.59, 299.3]; (c)
Design with the moment-bounding yield-aware optimization [48]:
$\mathbf{x}$=[300, 149.67, 300]; (d) Design with the proposed PoBO method:
$\mathbf{x}$=[300, 112.15, 300]. TABLE II: Optimization Results for MZI
Benchmark

Risk level $\epsilon$ | Method | Simulation # |  $\mathbb{E}_{\boldsymbol{{\xi}}}[\text{BW}]$ (GHz) |  $\Delta_{1}$ (%) |  $\Delta_{2}$ (%) | Yield (%)
---|---|---|---|---|---|---
0.05 | Moment [48] | 35 | 184.53 | 5.26 | 5.26 | 100
Proposed | 35 | 190.99 | 5.26 | 5.26 | 100
0.07 | Moment [48] | 35 | 187.02 | 7.53 | 7.53 | 100
Proposed | 35 | 192.10 | 7.53 | 7.2 | 99.7
0.1 | Moment [48] | 35 | 189.35 | 11.11 | 11.11 | 100
Proposed | 35 | 193.28 | 11.11 | 4.22 | 93.8
N/A* | BYO [15] | 2020 | 182.82 | N/A* | N/A* | 100

  * *

No risk level is defined for BYO method. Correspondingly, no gap $\Delta$ is
defined.

Figure 7: The probability density function of the optimized bandwidth of the
MZI by Bayesian yield optimization [15], moment bounding [48] and the proposed
PoBO (with $\epsilon=0.1$).

Baseline Methods. We choose the moment-bounding chance-constraint optimization
[48] as the baseline for comparison. On the photonic IC benchmarks, we further
compare our with the Bayesian yield optimization (BYO) method [15], a recent
state-of-the-art yield optimization approach.

Gap of chance constraints. We modify the indicator function (1) to define an
indicator function ${I_{i}}(\mathbf{x},\boldsymbol{{\xi}})$ for each
individual design constraint in yield definition:

${{I_{i}}(\mathbf{x},\boldsymbol{{\xi}})=\begin{cases}1,&\hat{y}_{i}(\mathbf{x},\boldsymbol{{\xi}})\leq
u_{i};\\\ 0,&\text{otherwise}\end{cases},\forall i=[n].}$

With $N$ random samples, the individual success rate for each design
constraint is evaluated as
${Y_{i}}(\mathbf{x})=\begin{matrix}\sum_{j=1}^{N}{I_{i}}(\mathbf{x},\boldsymbol{{\xi}}_{j})/N.\end{matrix}$
The gap for the $i$-th chance constraint is the relative difference between
$Y_{i}$ and the pre-specified success rate $1-\epsilon_{i}$:

$\Delta_{i}=\frac{{Y_{i}}(\mathbf{x})-(1-\epsilon_{i})}{1-\epsilon_{i}},\;\forall
i=[n].$ (15)

The chance-constrained optimization can always provide a solution to ensure
the yield requirement controlled by $\epsilon_{i}$’s if a feasible solution
exists. Therefore, $\Delta_{i}$ is always non-negative. Notice that we do not
attempt to achieve the highest yield. Instead, our goal is to avoid over-
conservative design while ensuring the pre-specified yield requirement.
Therefore, given a certain risk level, we prefer a smaller gap $\Delta_{i}$
and a less conservative design solution with better objective performance.

###  V-A Synthetic Function

We first consider a synthetic function with design variables
$\mathbf{x}\in\mathbf{X}={[-1,1]}^{2}$ and random parameters
$\boldsymbol{{\xi}}$ following a truncated Gaussian mixture model.
Specifically, we assume
$\mu(\boldsymbol{{\xi}})=\frac{1}{2}\mathcal{TN}(\bar{\boldsymbol{\mu}}_{1},\mathbf{\Sigma}_{1},\mathbf{a}_{1},\mathbf{b}_{1})+\frac{1}{2}\mathcal{TN}(\bar{\boldsymbol{\mu}}_{2},\mathbf{\Sigma}_{2},\mathbf{a}_{2},\mathbf{b}_{2})$,
with
$\bar{\boldsymbol{\mu}}_{1}=-\bar{\boldsymbol{\mu}}_{2}={[0.1,-0.1]}^{T}$,
$\mathbf{\Sigma}_{1}=\mathbf{\Sigma}_{2}=10^{-2}\begin{bmatrix}1&-0.75\\\
-0.75&1\end{bmatrix}$, $\mathbf{a}_{1}=-{[0.2,0.4]}^{T}$,
$\mathbf{a}_{2}=-{[0.4,0.2]}^{T}$, $\mathbf{b}_{1}={[0.4,0.2]}^{T}$ and
$\mathbf{b}_{2}={[0.2,0.4]}^{T}$. Here
$\mathcal{TN}(\bar{\boldsymbol{\mu}},\mathbf{\Sigma},\mathbf{a},\mathbf{b})$
denotes a distribution that is a normal distribution with mean
$\bar{\boldsymbol{\mu}}$ and variance $\mathbf{\Sigma}$ in the box
$[\mathbf{a},\mathbf{b}]$.

We consider the following chance-constrained optimization:

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[3(x_{1}+\xi_{1})-(x_{2}+\xi_{2}))]$
(16) $\displaystyle\mathop{\mathrm{s.t.}}$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}({(x_{1}+\xi_{1})}^{2}+(x_{2}+\xi_{2})\leq
1)\geq 1-\epsilon_{1},$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}({(x_{1}+\xi_{1})}^{2}-(x_{2}+\xi_{2})\leq
1)\geq 1-\epsilon_{2},$

where the two risk levels are set to be equal
$\epsilon_{1}=\epsilon_{2}=\epsilon$.

We use 2nd-order polynomials to approximate the three analytical functions and
bound the chance constraint with an order-10 optimal polynomial kinship
function. As shown in Table I, compared with the moment method [48], the
proposed PoBO method produces a better objective value and smaller gaps for
chance constraints while meeting the pre-specified yield requirement. Clearly,
a smaller $\epsilon_{i}$ produces a higher yield. The moment bounding method
fails to work when $\epsilon=0.01$ while our PoBO can still solve this
problem.

Figure 8: The schematic of a microring add-drop filter. Figure 9: The
probability density function of the optimized bandwidth of the microring
filter by Bayesian yield optimization [15], moment bounding [48] and the
proposed PoBO (with $\epsilon=0.1$).

###  V-B Mach-Zehnder Interferometer

We consider a third-order Mach-Zehnder interferometer (MZI) which consists of
three port coupling and two arms, as shown in Fig. 5. The coupling
coefficients $\tau$ between the MZ arms play an important role in the design,
whose relationship with the the gap $g$ (nm) is $\tau=\exp(-\frac{g}{260})$.
The gap variables $\mathbf{x}=[g_{1},g_{2},g_{3}]$ have the design space of
$\mathbf{X}={[100,300]}^{3}$. Their corresponding process variations
$\boldsymbol{{\xi}}$ follows a truncated Gaussian mixture distribution (see
Appendix B). We aim to maximize the expected 3-dB bandwidth (BW, in GHz) with
probability constraints on the crosstalk (XT, in dB) and the attenuation
($\alpha$, in dB) of the peak transmission. Therefore, the yield-aware chance-
constrained design is formulated as

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[\text{BW}(\mathbf{x},\boldsymbol{{\xi}})]$
(17) $\displaystyle\mathop{\mathrm{s.t.}}$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}(\text{XT}(\mathbf{x},\boldsymbol{{\xi}})\leq\text{XT}_{0})\geq
1-\epsilon_{1},$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}(\alpha(\mathbf{x},\boldsymbol{{\xi}})\leq\alpha_{0})\geq
1-\epsilon_{2}.$

The two risk levels are set to be equal $\epsilon_{1}=\epsilon_{2}=\epsilon$.
The thresholds of the crosstalk ($\text{XT}_{0}$) and the attenuation
($\alpha_{0}$) are -4 dB and 1.6 dB respectively.

We build three 2nd-order polynomial surrogate models for BW, XT, and $\alpha$,
respectively. We further bound the probabilistic yield constraints via an
order-5 optimal polynomial kinship function. The optimized results and
comparisons are listed in Table II. It shows that at the same risk level, the
proposed PoBO method can achieve larger bandwidth while meeting the yield
requirements and having smaller gaps for the chance constraints. The
simulation samples are the ones used for building surrogate models. The
proposed PoBO requires the same number of simulations as the moment bounding
method [48], and both of them require much fewer simulation samples than the
Bayesian yield optimization due to the efficient surrogate modeling. Fig. 6
compares the frequency response before and after the yield-aware optimization
with $\epsilon=0.1$. Our PoBO method has a higher expected bandwidth compared
with the Bayesian yield optimization, the moment bounding method, and the
initial design. Fig. 7 further shows the probability density of the optimized
bandwidth by different models. It clearly shows that our proposed method
produces the highest bandwidth while meeting the yield requirement.

###  V-C Microring Add-Drop Filter

We further consider the design of an optical add-drop filter consisting of
three identical silicon microrings coupled in series, as shown in Fig. 8. The
design variables are the coupling coefficients
$\mathbf{x}=[K_{1},K_{2},K_{3},K_{4}]$ that are to be optimized within the
interval of $\mathbf{X}={[0.3,0.6]}^{4}$. The process variations
$\boldsymbol{{\xi}}$ are described by a truncated Gaussian mixture model (see
Appendix B). The design problem is to maximize the expected 3-dB bandwidth
(BW, in GHz) with constraints on the extinction ratio (RE, in dB) of the
transmission at the drop port and the roughness ($\sigma_{\text{pass}}$, in
dB) of the passband that takes its standard deviation, formulated as:

$\displaystyle\max_{\mathbf{x}\in\mathbf{X}}$
$\displaystyle\mathbb{E}_{\boldsymbol{{\xi}}}[\text{BW}(\mathbf{x},\boldsymbol{{\xi}})]$
(18) $\displaystyle\mathop{\mathrm{s.t.}}$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}(\text{RE}(\mathbf{x},\boldsymbol{{\xi}})\geq\text{RE}_{0})\geq
1-\epsilon_{1},$
$\displaystyle\mathbb{P}_{\boldsymbol{{\xi}}}(\sigma_{\text{pass}}(\mathbf{x},\boldsymbol{{\xi}})\leq\sigma_{0})\geq
1-\epsilon_{2}.$

The two risk levels are set to be equal $\epsilon_{1}=\epsilon_{2}=\epsilon$.
The thresholds of the extinction ratio ($\text{RE}_{0}$) and the roughness of
the passband ($\sigma_{0}$) are 20 dB and 0.65 dB respectively.

Figure 10: The transmission curves of the microring add-drop filter. The grey
lines show the performance uncertainties. The orange and blue curves show the
transmission rates at the drop and through ports, respectively. The mean
values of the extinction ratio, bandwidth, and roughness are denoted as RE,
BW, and $\sigma_{\text{pass}}$, respectively. (a) The initial design:
$\mathbf{x}$=[0.45, 0.45, 0.45, 0.45]; (b) Design after Bayesian yield
optimization [15]: $\mathbf{x}$=[0.5758, 0.3718, 0.3720, 0.5746]; (c) Design
with the moment-bounding yield-aware optimization [48]: $\mathbf{x}$=[0.6,
0.3751, 0.3642 0.6]; (d) Design with the proposed PoBO optimization:
$\mathbf{x}$=[0.6, 0.3971, 0.3642, 0.6]. TABLE III: Optimization Results for
Microring Add-drop Filter benchmark

Risk level $\epsilon$ | Method | Simulation # |  $\mathbb{E}_{\boldsymbol{{\xi}}}[\text{BW}]$ (GHz) |  $\Delta_{1}$ (%) |  $\Delta_{2}$ (%) | Yield (%)
---|---|---|---|---|---|---
0.05 | Moment [48] | 65 | N/A* | N/A* | N/A* | N/A*
Proposed | 65 | 116.85 | 5.26 | 4.84 | 99.6
0.07 | Moment [48] | 65 | 112.64 | 7.53 | 7.42 | 99.9
Proposed | 65 | 120.05 | 7.53 | 6.67 | 99.2
0.1 | Moment [48] | 65 | 118.47 | 11.11 | 10.78 | 99.7
Proposed | 65 | 123.05 | 11.11 | 8.33 | 97.5
N/A | BYO [15] | 2020 | 117.42 | N/A | N/A | 95.1

  * *

The algorithm fails with no feasible solution.

Similarly, we build three 2nd-order polynomial surrogate models for BW, RE,
and $\sigma_{\text{pass}}$ and bound the chance constraints via an order-5
optimal polynomial kinship function. The optimized results and comparisons are
shown in Table III. The moment bounding [48] fails when $\epsilon=0.05$ since
no feasible solution is found under its over-conservative bounding. At all
risk levels, the proposed method can achieve larger bandwidth while meeting
the yield requirements and having smaller gaps for the chance constraints. As
shown in Fig. 9, the proposed PoBO has a higher expected bandwidth compared to
Bayesian yield optimization, existing yield-aware chance-constraint
optimization via moment bounding [48]. Fig. 10 shows the frequency response
before and after the yield-aware optimization with $\epsilon=0.1$. All the
reported design problems can be solved within a few minutes on a regular
laptop computer.

##  VI Conclusion and Future Work

This paper has proposed a novel Polynomial Bounding method for chance-
constrained yield-aware Optimization (PoBO) of photonic ICs with truncated
non-Gaussian correlated uncertainties. In PoBO, we first construct surrogate
models with a few simulation samples for the interested quantities based on
available uncertainty quantification solvers. To avoid over-conservative
design, we have proposed an optimal polynomial kinship function to tightly
bound the chance constraints. This bounding method can be efficiently
implemented without additional simulations. It also preserves the polynomial
form and enables seeking a globally optimal design. The proposed PoBO is
verified with a synthetic function, a Mach-Zehnder interferometer, and a
microring add-drop filter. In all experiments, the proposed PoBO has achieved
the yield requirements, produced tighter bounds on the chance constraints than
the state-of-the-art moment bounding method, and led to better design
objective performances with a few simulation samples. On the two photonic IC
examples, the proposed method has also reduced the simulation samples by
$58\times$ compared with Bayesian yield optimization.

The theoretical and numerical results of this work have laid the foundation of
many future topics. Possible extensions of this work include, but are not
limited to: (1) improved algorithms to handle many design parameters and
process variation, (2) formulations and algorithms to handle joint chance
constraint for yield descriptions, (3) PoBO with non-polynomial surrogates.
The proposed framework is very generic, and it can also be employed in other
applications beyond EDA, including probabilistic control of energy systems,
safety-critical control of autonomous systems, and so forth.

##  Appendix A Solution to Optimal Polynomial Kinship (8)

Given a $\rho$-order optimal kinship function $\kappa_{\rho}(\cdot)$, we
introduce two positive semidefinite matrices
$\mathbf{Y}_{1}\in\mathbb{R}^{(n_{1}+1)\times(n_{1}+1)}$ and
$\mathbf{Y}_{2}\in\mathbb{R}^{(n_{2}+1)\times(n_{2}+1)}$ with
$n_{1}=\lfloor(\rho-1)/2\rfloor$ and $n_{2}=\lfloor(\rho-2)/2\rfloor$. We
further define two series of Hankel matrices
$\mathbf{H}_{1,m}\in\mathbb{R}^{(n_{1}+1)\times(n_{1}+1)}$ and
$\mathbf{H}_{2,m}\in\mathbb{R}^{(n_{2}+1)\times(n_{2}+1)}$ as

$\mathbf{H}_{k,m}\left(i,j\right)=\begin{cases}1,&i+j=m+1\\\
0,&\text{otherwise}\end{cases},k=1,2.$ (19)

Based on the sum-of-square representation of a nonnegative univariate
polynomial, we can reformulate (8) as a finite dimensional semidefinite
programming (20), which can be handled by many efficient solvers and toolboxes
[59]. The detailed proof will be similar to the Corollary 1 of [54], where the
difference is the order of nonnegative polynomial.

$\centering\begin{array}[]{c}\min_{{\zeta_{0}},\ldots,{\zeta_{\rho}},{\mathbf{Y}_{1}},{\mathbf{Y}_{2}}}\;\sum_{i=0}^{\rho}\frac{(-1)^{i}}{i+1}\zeta_{i}\\\
\mathop{\mathrm{s.t.}}\;\zeta_{0}=1,\\\
\sum_{i=0}^{\rho}(-1)^{i}\zeta_{i}=0,\\\
\text{Tr}(\mathbf{Y}_{1}\mathbf{H}_{1,m})+\text{Tr}(\mathbf{Y}_{2}\mathbf{H}_{2,m})=\\\
\sum_{i=m+1}^{\rho}\frac{i!(-1)^{i-m-1}}{k!(i-m-1)!}\zeta_{i},m=0,1,\ldots,\rho-2,\\\
\mathbf{Y}_{k}\succeq 0,\;k=1,2.\end{array}\@add@centering$ (20)

##  Appendix B Details about benchmark setup

In the MZI benchmark (Sec. V-B), the process variations on the coupling
coefficients are described by a truncated Gaussian mixture model with two
components:

$\mu(\mathbf{\xi})=\frac{1}{2}\mathcal{TN}_{1}(\bar{\boldsymbol{\mu}}_{1},\mathbf{\Sigma}_{1},\mathbf{a}_{1},\mathbf{b}_{1})+\frac{1}{2}\mathcal{TN}_{2}(\bar{\boldsymbol{\mu}}_{2},\mathbf{\Sigma}_{2},\mathbf{a}_{2},\mathbf{b}_{2}),$
(21)

where $\bar{\boldsymbol{\mu}}_{1}=-\bar{\boldsymbol{\mu}}_{2}={[3,3,3]}^{T}$,
$\mathbf{\Sigma}_{1}=\mathbf{\Sigma}_{2}=3^{2}\begin{bmatrix}1&0.4&0.1\\\
0.4&1&0.4\\\
0.1&0.4&1\end{bmatrix},\mathbf{a}_{1}=-\mathbf{b}_{2}={[-6,-6,-6]}^{T},\mathbf{a}_{2}=-\mathbf{b}_{1}={[-12,-12,-12]}^{T}.$

In the microring benchmark (Sec. V-C), the process variations on the coupling
coefficients are described the same as Eq. (21) with different parameters
$\bar{\boldsymbol{\mu}}_{1}=-\bar{\boldsymbol{\mu}}_{2}={[0.03,0.03,0.03,0.03]}^{T}$,
$\mathbf{\Sigma}_{1}=\mathbf{\Sigma}_{2}=0.03^{2}\begin{bmatrix}1&0.4&0.1&0.4\\\
0.4&1&0.4&0.1\\\ 0.1&0.4&1&0.4\\\
0.4&0.1&0.4&1\end{bmatrix},\mathbf{a}_{1}=-\mathbf{b}_{2}={[-0.06,-0.06,-0.06,-0.06]}^{T},\mathbf{a}_{2}=-\mathbf{b}_{1}={[-0.12,-0.12,-0.12,-0.12]}^{T}.$

## References

  * [1] G. Gielen, P. De Wit, E. Maricau, J. Loeckx, J. Martin-Martinez, B. Kaczer, G. Groeseneken, R. Rodriguez, and M. Nafria, “Emerging yield and reliability challenges in nanometer CMOS technologies,” in _Proc. Design, Autom. Test Eur. Conf. Exhibit._ , 2008, pp. 1322–1327.
  * [2] X. Chen, M. Mohamed, Z. Li, L. Shang, and A. R. Mickelson, “Process variation in silicon photonic devices,” _Appl. Opt._ , vol. 52, no. 31, pp. 7638–7647, 2013.
  * [3] T. Lipka, J. Müller, and H. K. Trieu, “Systematic nonuniformity analysis of amorphous silicon-on-insulator photonic microring resonators,” _J. Light. Technol._ , vol. 34, no. 13, pp. 3163–3170, 2016.
  * [4] W. Bogaerts, Y. Xing, and U. Khan, “Layout-aware variability analysis, yield prediction, and optimization in photonic integrated circuits,” _IEEE J. Sel. Top. Quantum Electron._ , vol. 25, no. 5, pp. 1–13, 2019.
  * [5] K. J. Antreich, H. E. Graeb, and C. U. Wieser, “Circuit analysis and optimization driven by worst-case distances,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 13, no. 1, pp. 57–71, 1994\.
  * [6] F. Gong, Y. Shi, H. Yu, and L. He, “Variability-aware parametric yield estimation for analog/mixed-signal circuits: Concepts, algorithms, and challenges,” _IEEE Des. Test_ , vol. 31, no. 4, pp. 6–15, 2014.
  * [7] Z. Lu, J. Jhoja, J. Klein, X. Wang, A. Liu, J. Flueckiger, J. Pond, and L. Chrostowski, “Performance prediction for silicon photonics integrated circuits with layout-dependent correlated manufacturing variability,” _Opt. Express_ , vol. 25, no. 9, pp. 9712–9733, 2017.
  * [8] J. Pond, J. Klein, J. Flückiger, X. Wang, Z. Lu, J. Jhoja, and L. Chrostowski, “Predicting the yield of photonic integrated circuits using statistical compact modeling,” in _Integrated Optics: Physics and Simulations III_ , vol. 10242, 2017, p. 102420S.
  * [9] Y. Xu, K.-L. Hsiung, X. Li, L. T. Pileggi, and S. P. Boyd, “Regular analog/RF integrated circuits design using optimization with recourse including ellipsoidal uncertainty,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 28, no. 5, pp. 623–637, 2009.
  * [10] G. Yu and P. Li, “Yield-aware analog integrated circuit optimization using geostatistics motivated performance modeling,” in _Proc. Intl. Conf. Computer Aided Design_ , 2007, pp. 464–469.
  * [11] S. K. Tiwary, P. K. Tiwary, and R. A. Rutenbar, “Generation of yield-aware Pareto surfaces for hierarchical circuit design space exploration,” in _Proc. Design Autom. Conf_ , 2006, pp. 31–36.
  * [12] Y. Li and V. Stojanovic, “Yield-driven iterative robust circuit optimization algorithm,” in _Proc. Design Autom. Conf_ , 2009, pp. 599–604.
  * [13] B. Liu, F. V. Fernández, and G. G. Gielen, “Efficient and accurate statistical analog yield optimization and variation-aware circuit sizing based on computational intelligence techniques,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 30, no. 6, pp. 793–805, 2011\.
  * [14] M. Barros, J. Guilherme, and N. Horta, “Analog circuits optimization based on evolutionary computation techniques,” _Integration_ , vol. 43, no. 1, pp. 136–155, 2010.
  * [15] M. Wang, F. Yang, C. Yan, X. Zeng, and X. Hu, “Efficient Bayesian yield optimization approach for analog and SRAM circuits,” in _Proc. Design Autom. Conf_ , 2017, pp. 1–6.
  * [16] A. Singhee, S. Singhal, and R. A. Rutenbar, “Practical, fast Monte Carlo statistical static timing analysis: Why and how,” in _Proc. Intl. Conf. Computer Aided Design_ , 2008, pp. 190–195.
  * [17] A. Papoulis and H. Saunders, _Probability, random variables and stochastic processes_. McGraw-Hill, 2001.
  * [18] C. Gu and J. Roychowdhury, “An efficient, fully nonlinear, variability-aware non-Monte-Carlo yield estimation procedure with applications to SRAM cells and ring oscillators,” in _Proc. Asia South Pac. Design Autom. Conf._ , 2008, pp. 754–761.
  * [19] F. Gong, X. Liu, H. Yu, S. X. Tan, J. Ren, and L. He, “A fast non-Monte-Carlo yield analysis and optimization by stochastic orthogonal polynomials,” _ACM Trans. Des. Autom. Electron. Syst._ , vol. 17, no. 1, pp. 1–23, 2012.
  * [20] Z. Gao and R. Rohrer, “Efficient non-Monte-Carlo yield estimation,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , pp. 1–1, 2021\.
  * [21] X. Shi, H. Yan, Q. Huang, J. Zhang, L. Shi, and L. He, “Meta-model based high-dimensional yield analysis using low-rank tensor approximation,” in _Proc. Design Autom. Conf_ , 2019, pp. 1–6.
  * [22] J. Yao, Z. Ye, and Y. Wang, “An efficient SRAM yield analysis and optimization method with adaptive online surrogate modeling,” _IEEE Trans. Very Large Scale Integr. (VLSI) Syst._ , vol. 23, no. 7, pp. 1245–1253, 2014.
  * [23] X. Li, J. Le, P. Gopalakrishnan, and L. T. Pileggi, “Asymptotic probability extraction for nonnormal performance distributions,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 26, no. 1, pp. 16–37, 2006\.
  * [24] X. Li, P. Gopalakrishnan, Y. Xu, and T. Pileggi, “Robust analog/RF circuit design with projection-based posynomial modeling,” in _Proc. Intl. Conf. Computer Aided Design_ , 2004, pp. 855–862.
  * [25] X. Li, Y. Zhan, and L. T. Pileggi, “Quadratic statistical $max$ approximation for parametric yield estimation of analog/rf integrated circuits,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 27, no. 5, pp. 831–843, 2008.
  * [26] A. Ciccazzo, G. Di Pillo, and V. Latorre, “A SVM surrogate model-based method for parametric yield optimization,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 35, no. 7, pp. 1224–1228, 2015.
  * [27] H. Ma, E.-P. Li, A. C. Cangellaris, and X. Chen, “Support vector regression-based active subspace (SVR-AS) modeling of high-speed links for fast and accurate sensitivity analysis,” _IEEE Access_ , vol. 8, pp. 74 339–74 348, 2020.
  * [28] A. C. Sanabria-Borbón, S. Soto-Aguilar, J. J. Estrada-López, D. Allaire, and E. Sánchez-Sinencio, “Gaussian-process-based surrogate for optimization-aided and process-variations-aware analog circuit design,” _Electronics_ , vol. 9, no. 4, p. 685, 2020.
  * [29] M. Wang, W. Lv, F. Yang, C. Yan, W. Cai, D. Zhou, and X. Zeng, “Efficient yield optimization for analog and SRAM circuits via Gaussian process regression and adaptive yield estimation,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 37, no. 10, pp. 1929–1942, 2017.
  * [30] Y. Wang, M. Orshansky, and C. Caramanis, “Enabling efficient analog synthesis by coupling sparse regression and polynomial optimization,” in _Proc. Design Autom. Conf_ , 2014, pp. 1–6.
  * [31] D. Xiu and G. E. Karniadakis, “Modeling uncertainty in flow simulations via generalized polynomial chaos,” _J. Comput. Phys._ , vol. 187, no. 1, pp. 137–167, 2003.
  * [32] Z. Zhang, T. A. El-Moselhy, I. M. Elfadel, and L. Daniel, “Stochastic testing method for transistor-level uncertainty quantification based on generalized polynomial chaos,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 32, no. 10, pp. 1533–1545, 2013.
  * [33] R. Trinchero and F. G. Canavero, “Combining LS-SVM and GP regression for the uncertainty quantification of the EMI of power converters affected by several uncertain parameters,” _IEEE Trans. Electromagn. Compat._ , vol. 62, no. 5, pp. 1755–1762, 2020.
  * [34] X. Li, “Finding deterministic solution from underdetermined equation: large-scale performance variability modeling of analog/RF circuits,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 29, no. 11, pp. 1661–1668, 2010.
  * [35] Z. Zhang, X. Yang, I. V. Oseledets, G. E. Karniadakis, and L. Daniel, “Enabling high-dimensional hierarchical uncertainty quantification by ANOVA and tensor-train decomposition,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 34, no. 1, pp. 63–76, 2014.
  * [36] Z. Zhang, T.-W. Weng, and L. Daniel, “Big-data tensor recovery for high-dimensional uncertainty quantification of process variations,” _IEEE Trans. Compon. Packag. Manuf. Technol._ , vol. 7, no. 5, pp. 687–697, 2016.
  * [37] Z. He and Z. Zhang, “High-dimensional uncertainty quantification via active and rank-adaptive tensor regression,” in _Proc. Electr. Perform. Electron. Packag. Syst._ , 2020, pp. 1–3.
  * [38] Z. He and Z. Zhang, “High-dimensional uncertainty quantification via tensor regression with rank determination and adaptive sampling,” _accepted by IEEE Trans. Compon. Packag. Manuf. Technol., doi: 10.1109/TCPMT.2021.3093432_ , 2021.
  * [39] Z. He, W. Cui, C. Cui, T. Sherwood, and Z. Zhang, “Efficient uncertainty modeling for system design via mixed integer programming,” in _Proc. Intl. Conf. Computer Aided Design_ , 2019, pp. 1–8.
  * [40] C. Cui and Z. Zhang, “Stochastic collocation with non-Gaussian correlated process variations: Theory, algorithms, and applications,” _IEEE Trans. Compon. Packag. Manuf. Technol._ , vol. 9, no. 7, pp. 1362–1375, 2018.
  * [41] P. Manfredi, D. V. Ginste, D. De Zutter, and F. G. Canavero, “Stochastic modeling of nonlinear circuits via SPICE-compatible spectral equivalents,” _IEEE Trans. Circuits Syst. I, Reg. Papers_ , vol. 61, no. 7, pp. 2057–2065, 2014.
  * [42] M. Ahadi and S. Roy, “Sparse linear regression (SPLINER) approach for efficient multidimensional uncertainty quantification of high-speed circuits,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , vol. 35, no. 10, pp. 1640–1652, 2016.
  * [43] A. Kaintura, T. Dhaene, and D. Spina, “Review of polynomial chaos-based methods for uncertainty quantification in modern integrated circuits,” _Electronics_ , vol. 7, no. 3, p. 30, 2018.
  * [44] A. Waqas, P. Manfredi, and D. Melati, “Performance variability analysis of photonic circuits with many correlated parameters,” _J. Light. Technol._ , 2021.
  * [45] A. Waqas, D. Melati, P. Manfredi, and A. Melloni, “Stochastic process design kits for photonic circuits based on polynomial chaos augmented macro-modelling,” _Opt. Express_ , vol. 26, no. 5, pp. 5894–5907, 2018.
  * [46] T.-W. Weng, Z. Zhang, Z. Su, Y. Marzouk, A. Melloni, and L. Daniel, “Uncertainty quantification of silicon photonic devices with correlated and non-Gaussian random parameters,” _Opt. Express_ , vol. 23, no. 4, pp. 4242–4254, 2015.
  * [47] T.-W. Weng, D. Melati, A. Melloni, and L. Daniel, “Stochastic simulation and robust design optimization of integrated photonic filters,” _Nanophotonics_ , vol. 6, no. 1, pp. 299–308, 2017.
  * [48] C. Cui, K. Liu, and Z. Zhang, “Chance-constrained and yield-aware optimization of photonic ICs with non-Gaussian correlated process variations,” _IEEE Trans. Comput.-Aided Design Integr. Circuits Syst._ , 2020.
  * [49] A. Shapiro, D. Dentcheva, and A. Ruszczyński, _Lectures on stochastic programming: modeling and theory_. SIAM, 2014.
  * [50] A. Mesbah, S. Streif, R. Findeisen, and R. D. Braatz, “Stochastic nonlinear model predictive control with probabilistic constraints,” in _Proc. Am. Control Conf._ , 2014, pp. 2413–2419.
  * [51] M. P. Vitus, Z. Zhou, and C. J. Tomlin, “Stochastic control with uncertain parameters via chance constrained control,” _IEEE Trans. Autom. Control_ , vol. 61, no. 10, pp. 2892–2905, 2015.
  * [52] Z. Wang, C. Shen, F. Liu, X. Wu, C.-C. Liu, and F. Gao, “Chance-constrained economic dispatch with non-Gaussian correlated wind power uncertainty,” _IEEE Trans. Power Syst._ , vol. 32, no. 6, pp. 4880–4893, 2017.
  * [53] B. P. Van Parys, P. J. Goulart, and D. Kuhn, “Generalized Gauss inequalities via semidefinite programming,” _Math. Program._ , vol. 156, no. 1-2, pp. 271–302, 2016.
  * [54] C. Feng, F. Dabbene, and C. M. Lagoa, “A kinship function approach to robust and probabilistic optimization under polynomial uncertainty,” _IEEE Trans. Autom. Control_ , vol. 56, no. 7, pp. 1509–1523, 2010.
  * [55] G. C. Calafiore and L. El Ghaoui, “On distributionally robust chance-constrained linear programs,” _J. Optim. Theory Appl._ , vol. 130, no. 1, pp. 1–22, 2006.
  * [56] D. Henrion, J.-B. Lasserre, and J. Löfberg, “Gloptipoly 3: moments, optimization and semidefinite programming,” _Optim. Methods Softw._ , vol. 24, no. 4-5, pp. 761–779, 2009.
  * [57] G. H. Golub and J. H. Welsch, “Calculation of Gauss quadrature rules,” _Math. Comp._ , vol. 23, no. 106, pp. 221–230, 1969.
  * [58] T. Gerstner and M. Griebel, “Numerical integration using sparse grids,” _Numer. Algorithms_ , vol. 18, no. 3, pp. 209–232, 1998.
  * [59] J. Lofberg, “YALMIP : a toolbox for modeling and optimization in MATLAB,” in _Proc. Intl. Conf. Robot. Autom._ , 2004, pp. 284–289.