\section{Introduction}
\label{sec:introduction}

In deterministic computer experiments, we use computer codes to study the input/output relationship of some complex physical, economical, or engineering models, e.g. large eddy simulations for rocket engine injector design \parencite{mak2018surrogate}. However, computer simulations are often time-consuming. Thus, the first step is to build a computationally cheap surrogate model that approximates the expensive computer code using some offline simulation runs \parencite{santner2018cdoe}.

Space-filling designs are commonly used for constructing the experimental designs where we run the computer simulations. Since we do not have a priori information about the input/output relationship, it is important to have the design points well spread out across the entire design region $\mathcal{X} \subseteq\mathbb{R}^p$. Minimax and maximin are the two popular space-filling measures proposed by \textcite{johnson1990design}. A minimax design aims to minimize the maximum distance from any point in $\mathcal{X}$ to the closest design point, whereas a maximin design maximizes the minimum distance between any two design points. Due to computational tractability, the maximin measure is more commonly used in the literature, which is followed in this paper as well. However, maximin designs are often collapsing, that is, some design points share the same value in one-dimensional projections. Latin hypercube designs \parencite[LHDs;][]{mckay1979lhd} are developed for having good projection of each factor, which can be further improved by integrating it with other space-filling criteria such as maximin \parencite{morris1995MmLHD}. However, maximin LHDs can only ensure good one-dimensional projection and full-dimensional space-fillingness. The maximum projection (MaxPro) designs \parencite{joseph2015maxpro}, on the other hand, are able to achieve good space-filling properties on projections to all subsets of factors.

Most literature on space-filling designs focus on  bounded rectangular region $\mathcal{X} = \prod_{d=1}^{p}[a_d,b_d]\subseteq\mathbb{R}^{p}$. However, in real world applications such as the welded beam design problem \parencite{dong2018cgo} and the NASA speed reducer design problem \parencite{liu2017cgo}, we frequently need to deal with non-rectangular bounded design space:
\begin{equation}
  \label{eq:space}
  \mathcal{X} = \bigg\{x \in \prod_{d=1}^{p}[a_d,b_d]: g_k(x) \leq 0\; \forall k = 1,\ldots,K \bigg\} \; ,
\end{equation}
where the rectangular shape is jeopardized by the $K$ inequality constraints $\{g_k(x)\leq 0\}_{k=1}^{K}$. For simplicity, let us consider the bounded space of a unit hypercube, that is $a_d = 0,\; b_d = 1\; \forall d = 1,\ldots,p$. This is possible since we can always re-scale the factors. Figure~\ref{fig:mot-space} shows a two-dimensional design space $\mathcal{X}$ obtained by three nonlinear inequality constraints \eqref{eq:MOT}. The non-convex, non-rectangular shape with extremely small feasibility ratio makes it challenging to construct space-filling designs. \par

\begin{figure}[t!]
    \centering
    \includegraphics[width=0.4\textwidth]{figures/motivation/feasible.png}
    \caption{A two-dimensional non-convex, non-rectangular bounded design space $\mathcal{X}$ (in red) with feasibility ratio of only 0.53\% (Table \ref{tab:benchmark_info}) due to the three nonlinear inequality constraints \eqref{eq:MOT}.}
    \label{fig:mot-space}
\end{figure}

Two primary approaches are proposed in the literature for constructing designs in non-rectangular design space. One approach is to directly employ general purpose constrained optimization techniques  \parencite{trosset1999cMm, stinstra2003cMm, kang2019}. However, this approach can be computationally very expensive and can be limited by the type of constraints and design properties (such as projections) it can handle. The alternative approach instead relies on a two-step process:
\begin{itemize}
  \item \textbf{Candidate Generation:} generate a large set of uniformly distributed candidates in $\mathcal{X}$.
  \item \textbf{Design Construction:} choose points from the set of candidates by a desired criterion.
\end{itemize}
The flexibility of choosing any design criterion aforementioned in the construction step easily allows for both space-filling and noncollapsing properties in the resulting designs, but how to efficiently generate good quality candidate points remains the \textit{key difficulty} of this approach.  The main objective of this paper is to propose an efficient method to generate good quality candidate points that are suitable for constructing maximin designs.

Several candidates generation methods have been discussed in the literature. If the desired space $\mathcal{X}$ is regularly-shaped, e.g. simplex and circle, where the closed-form inverse Rosenblatt transform exists, we can obtain uniform samples in $\mathcal{X}$ by applying the inverse transform on a set of low-discrepancy sequence in $[0,1]^{p}$ \parencite{Fang1994}. However, we generally cannot compute the inverse Rosenblatt transform for arbitrary irregularly-shaped space $\mathcal{X}$. An alternative solution is to perform acceptance/rejection sampling on a large set of uniformly distributed points in $[0,1]^{p}$, such as grid points \parencite{pratola2017cMm}, Latin hypercube samples \parencite{wu2019cMm}, and quasi-random points \parencite{joseph2016maxpro}. For the design space $\mathcal{X}$ considered in Figure~\ref{fig:mot-space}, given its small feasibility ratio of 0.53\%, on average only 5 out of 1{,}000 samples in the unit hypercube would land in the design space, indicating that the one-step acceptance/rejection approach can be highly inefficient. One remedy is to iterate between acceptance/rejection sampling and candidate augmentation \parencite{draguljic2012cMm}. To benefit from simulated annealing \parencite{kirkpatrick1983sa}, the multi-step acceptance/rejection can be performed on a sequence of shrinking regions \parencite[Subset Simulation;][]{bect2017ss}, and this idea is further improved using the probabilistic constraint, leading to the Sequentially Constrained Monte Carlo \parencite[SCMC;][]{golchi2015scmc,golchi2016scmc}. However, SCMC suffers the same issue of Monte Carlo sampling: many samples are repeated or are very close to each other, which add minimal value for the ultimate goal of constructing a maximin design. Moreover, by having the samples well spread out, fewer proposed samples are required to cover the entire design space, and thus fewer evaluations of the constraints, which is beneficial when the constraints are expensive to evaluate. Minimum Energy Design (MinED) is a state-of-the-art deterministic sampling method for simulating well-spaced samples for any distribution \parencite{joseph2015mined,joseph2019mined}. When the target distribution is uniform, the MinED is equivalent to the maximin design, showing its strong connection to the problem considered in this paper. Thus, by incorporating the probabilistic constraints from SCMC in MinED, we propose the Constrained Minimum Energy Design (CoMinED) as a more efficient approach for generating good quality design candidate samples in arbitrarily constrained space. \par

The paper is organized as follows. Section~\ref{sec:existing_algorithms} reviews the existing candidates generation and designs construction algorithms. Section~\ref{sec:constrained_minimum_energy_design} discusses minimum energy design and proposes the constrained minimum energy design (CoMinED). Section~\ref{sec:simulation_results} demonstrates the improvement of the proposed CoMinED with extensive simulation studies. We conclude the article with some remarks and future research directions in Section~\ref{sec:conclusion}.

\section{Existing Algorithms}
\label{sec:existing_algorithms}

\subsection{Candidate Generation}
\label{subsec:candidates_generation}

\subsubsection{Acceptance/Rejection Sampling}
\label{subsubsec:acceptance_rejection_sampling}

The simplest approach to generate candidates in any constrained space $\mathcal{X}$ is to first simulate large set of uniformly distributed samples in a rectangular region that contains $\mathcal{X}$, and then apply acceptance/rejection sampling based on the constraints to keep only the feasible samples. Figure~\ref{fig:mot-lhd} shows that only 14 out of the 2{,}385 randomized Latin hypercube samples in $[0,1]^2$ are feasible for the motivation problem in Figure~\ref{fig:mot-space}. The Latin hypercube samples are generated using R package \texttt{lhs} \parencite{carnell2019lhs}. Apart from the low percentage of feasible samples, we can see that these 14 points do not cover the feasible space uniformly well, showing that the one-step acceptance/rejection sampling performs poorly on constrained design problems with very small feasibility ratio. In fact, this issue becomes more serious in higher dimensional problems.

\begin{figure}[t!]
  \centering
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/lhd_all.png}
      \caption{2{,}385 candidate samples}
    \end{subfigure}%
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/lhd_feasible.png}
      \caption{14 feasible samples}
    \end{subfigure}%
  \caption{Left panel shows 2{,}385 randomized Latin hypercube samples from $[0,1]^2$. Right panel shows the 14 feasible candidate samples in $\mathcal{X}$ after applying acceptance/rejection sampling.}
  \label{fig:mot-lhd}
\end{figure}

The inefficiency of the one-step acceptance/rejection sampling results from wasting majority of resources exploring the unit hypercube rather than the target space $\mathcal{X}$. Thus, a $T$-step acceptance/rejection sampling on a sequence of shrinking regions $[0,1]^{p} = \mathcal{X}_0 \supset \mathcal{X}_1 \supset \cdots \supset \mathcal{X}_T = \mathcal{X}$ would allow for exploitation of the important region that is likely feasible. This is known as the subset simulation for estimating the probability of failure in reliability analysis \parencite{bect2017ss}. Consider a non-rectangular bounded design space that is defined by one inequality constraint, $\mathcal{X} = \{x\in[0,1]^{p}:g(x) \leq 0\}$. The subset simulation defines $\{\mathcal{X}_{t}\}_{t=0}^{T}$ by introducing a decreasing sequence of thresholds, $\infty = u_0 > u_1 > \cdots > u_T = 0$, such that $\mathcal{X}_{t} = \{x\in[0,1]^{p}:g(x) \leq u_t\}$. The acceptance/rejection sampling is equivalent to sampling from a indicator function, so we can also view the subset simulation as sampling from the sequence of distributions $\{\mathbbm{1}_{\mathcal{X}_{t}}(x) = \mathbbm{1}(g(x) \leq u_t)\}_{t=0}^{T}$, and Sequential Monte Carlo (SMC) sampling can be applied. The sampling step in SMC is usually by Markov Chain Monte Carlo \parencite[MCMC;][]{robert2013mc}, but poor performances of MCMC on indicator function are observed in practice. One solution is to replace the hard constraint $g(x)\leq 0$ with a probabilistic constraint, leading to the Sequentially Constrained Monte Carlo, which is discussed next.

\subsubsection{Sequentially Constrained Monte Carlo}
\label{subsubsec:sequentially_constrained_monte_carlo}

\begin{algorithm}[t!]
\SetAlgoLined
  \textbf{Design Space:} $\mathcal{X} = \{x\in[0,1]^{p}:g_k(x)\leq 0 \; \forall k = 1,\ldots,K\}$. \\
  \textbf{Initialization:}
  \begin{itemize}
    \setlength{\itemsep}{0pt}
    \setlength{\parskip}{0pt}
    \setlength{\parsep}{0pt}
    \item set the increasing sequence of rigidity parameters $0 = \tau_{0} < \tau_{1} < \cdots < \tau_{T}$.
    \item simulate the initial $M$ samples $\{x_m^{(0)}\}_{m=1}^{M}$ from $[0,1]^{p}$.
  \end{itemize}
  \vspace{-\topsep}
  \vspace{2mm}
  \For{$t = 1,\ldots,T$}{
    $\bullet$ \textbf{Weighting:} compute the importance weight,
    \begin{equation}
      \label{eq:scmc1}
      w_{m}^{(t)} = \rho_{\tau_{t}}(x_{m}^{(t-1)}) / \rho_{\tau_{t-1}}(x_{m}^{(t-1)}) \; ,
    \end{equation}
    where $\rho_{\tau}(\cdot)$ is defined in \eqref{eq:pc3}. Normalize the weight by $\bar{w}_{m}^{(t)} = w_{m}^{(t)} / \sum_{i=1}^{M}w_{i}^{(t)}$. \\
    $\bullet$ \textbf{Resample:} draw $M$ samples $\{y_m^{(t)}\}_{m=1}^{M}$ from $\sum_{m=1}^{M}\bar{w}_{m}^{(t)}\delta(x - x_m^{(t-1)})$.\\
    $\bullet$ \textbf{Sampling:} for $m = 1,\ldots,N$, draw $x_m^{(t)} \sim K_{\sigma^{(t)}}(y_m^{(t)}, \cdot)$ where $K_{\sigma^{(t)}}(y_m^{(t)}, \cdot)$ is a Markov kernel with target distribution $\rho_{\tau_{t}}$ and adaptive scale $\sigma^{(t)}$ where
    \begin{equation}
      \label{eq:scmc2}
      \sigma^{(t)} = \bigg(\mbox{75\% quantile of } \{\min_{j\neq m}\lVert x_{m}^{(t-1)} - x_{j}^{(t-1)}\rVert\}_{m=1}^{M}\bigg)/\sqrt{p} \;.
    \end{equation} \\

  }
  \vspace{2mm}
  \textbf{Return:} all particles $\{x_{m}^{(t)}\}_{m=1}^{M}{}_{t=0}^{T}$ that are in $\mathcal{X}$. \\
 \caption{Adaptive Sequentially Constrained Monte Carlo.}
 \label{algo:scmc}
\end{algorithm}

For any inequality constraint $g(x) \leq 0$, \textcite{golchi2015scmc} proposed the following probabilistic relaxation using the probit function,
\begin{equation}
  \label{eq:pc1}
  \rho_{\tau}(x) = \Phi(-\tau g(x))\; ,
\end{equation}
where $\Phi$ is the standard normal cumulative distribution function and $\tau$ is the parameter that controls the \textit{rigidity} of the constraint. We can see that the function $\rho_{\tau}$ assigns value close to 1 for $x$ that meets the constraint and value close to 0 otherwise. Moreover, in the limit,
\begin{equation}
  \label{eq:pc2}
  \lim_{\tau\to\infty}\rho_{\tau}(x) = \lim_{\tau\to\infty}\Phi(-\tau g(x)) = \mathbbm{1}(g(x) \leq 0)\; .
\end{equation}
The above can be generalized to multiple inequality constraints $\{g_{k}(x)\leq 0\}_{k=1}^{K}$ by
\begin{equation}
  \label{eq:pc3}
  \rho_{\tau}(x) = \prod_{k=1}^{K}\Phi(-\tau g_{k}(x))\; .
\end{equation}
This leads to the Sequentially Constrained Monte Carlo (SCMC) that replaces the sequence of indicator functions $\{\mathbbm{1}_{\mathcal{X}_{t}}\}_{t=0}^{T}$ in the subset simulation by the sequence of probabilistic constraint functions $\{\rho_{\tau_{t}}\}_{t=0}^{T}$ defined in \eqref{eq:pc3} with an increasing sequence $0 = \tau_{0} < \tau_{1} < \cdots < \tau_{T}$ where $\tau_{T}$ is a large constant, e.g. $10^{6}$. However, in the SCMC algorithm of \textcite{golchi2015scmc}, a pre-fixed normal distribution proposal is used in the Markov kernel of the MCMC step, but how to pick the scale of the normal proposal is difficult for high dimensional problem with small feasible region. Thus, for a more robust comparison to our proposed approach, we improve the SCMC method by allowing adaptation of the Markov kernel. At each iteration, the scale (standard deviation) of the normal proposal is adapted to be the 75\% quantile of the prior step samples' distances to their closest neighbors divided by the square root of the problem dimension. This adaptive kernel shows robust performance for majority of the benchmark problems considered in this paper. Algorithm~\ref{algo:scmc} details the \textit{adaptive SCMC} algorithm for generating large number of uniformly distributed samples from any design space $\mathcal{X}$ with arbitrary number of constraints.

Figure~\ref{fig:mot-scmc} shows the performance of the adaptive SCMC on the motivation problem in Figure~\ref{fig:mot-space} with $M = 265$, $T = 8$, and $\{\tau_{t}\}_{t=0}^{8} = [0,e^{1},e^{2},e^{3},e^{4},e^{5},e^{6},e^{7},10^{6}]$. The 265 Sobol' points in $[0,1]^2$ simulated by the R package \texttt{randtoolbox} \parencite{christophe2019randtoolbox} are used as the initial candidate set. We can see that now with 2{,}385 samples, the adaptive SCMC yields 1{,}205 feasible samples, and they cover the feasible space much better than the the one-step acceptance/rejection sampling on the Latin hypercube samples (Figure~\ref{fig:mot-lhd}), but we can still see large gaps left unexplored in the feasible space.

\begin{figure}[t!]
  \centering
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/scmc_all.png}
      \caption{2{,}385 candidate samples}
    \end{subfigure}%
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/scmc_feasible.png}
      \caption{1{,}205 feasible samples}
    \end{subfigure}%
  \caption{Left panel shows 2{,}385 candidate samples from applying adaptive SCMC (Algorithm~\ref{algo:scmc}) on the design space $\mathcal{X}$ defined by \eqref{eq:MOT} with $M = 265$, $T = 8$, and $\{\tau_{t}\}_{t=0}^{8} = [0,e^{1},e^{2},e^{3},e^{4},e^{5},e^{6},e^{7},10^{6}]$. Red circles indicate the initial candidate set of 265 Sobol' points. Right panel shows the 1{,}205 feasible candidate samples.}
  \label{fig:mot-scmc}
\end{figure}

\subsection{Design Construction}
\label{subsec:designs_construction}
The ultimate goal is to construct an $n$-point design $\mathcal{D}_{n} = \{x_{i}\in\mathcal{X}\}_{i=1}^{n}$ in $\mathcal{X}$ that achieves some good design properties. From the candidate generation step, we obtain a finite set of $N$ ($N \geq n$) candidate points $\mathcal{C}_{N} = \{y_j\in\mathcal{X}\}_{j=1}^{N}$ that are approximately uniformly distributed in $\mathcal{X}$. The next step is to find the $n$ samples from the candidate set that maximize a desired design criterion $\psi$, that is to solve
\begin{equation}
  \label{eq:dc1}
  \arg\max_{\mathcal{D}_{n}\subseteq \mathcal{C}_{N}} \psi(\mathcal{D}_{n}) \; .
\end{equation}
In the case of the maximin design, $\psi(\mathcal{D}_{n}) = \min_{x_i,x_j\in \mathcal{D}_{n}; i\neq j}\lVert x_i - x_j\rVert_{2}$ where $\lVert \cdot \rVert_{2}$ is the Euclidean distance. Many stochastic optimization algorithms have been developed to efficiently address the combinatorial optimization problem in \eqref{eq:dc1}, including local search, threshold accepting, simulated annealing \parencite{morris1995MmLHD}, enhanced stochastic evolutionary \parencite{jin2003MmLHD,wu2019cMm}. See \textcite{fang2005design} Chapter 4 for a detailed review of the above methods. On the other hand, \textcite{kennard1969cdoe} proposed a one-point-at-a-time greedy algorithm for solving \eqref{eq:dc1}. The idea is that by having a $m$-point design $\mathcal{D}_{m}$ ($m < n$), we generate the $(m+1)$-th point by
\begin{equation}
  \label{eq:dc2}
  x_{m+1} = \arg\max_{x\in \mathcal{C}_{N}\backslash \mathcal{D}_{m}}\psi(\mathcal{D}_{m}\cup\{x\}) \; .
\end{equation}
The one-point-at-a-time greedy procedure is also employed in the R package \texttt{mined} \parencite{wang2019mined} for generating minimum energy design and the R package \texttt{MaxPro} \parencite{ba2018maxpro} for design augmentation. Although the one-point-at-a-time greedy algorithm results in a local optimum, it is efficient and shows good empirical performance in practice. We also employ this greedy procedure in the constrained minimum energy design discussed in the next section.

\section{Constrained Minimum Energy Design}
\label{sec:constrained_minimum_energy_design}

\subsection{Minimum Energy Design}
\label{subsec:minimum_energy_design}

We begin by formally defining the minimum energy design (MinED).

\begin{definition}{\textbf{\parencite[Minimum Energy Design;][]{joseph2015mined}}}
\label{df:minimum_energy_designs}
Let $\pi$ be the target probability density function. An n-point minimum energy design of $\pi$ is the optimal solution of
\begin{equation}
  \label{eq:mined1}
  \arg\min_{\mathcal{D}_n\in\mathbb{D}_n}\sum_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{q(x_i)q(x_j)}{\lVert x_i - x_j\rVert_{2}}\; ,
\end{equation}
where $\mathbb{D}_{n} = \{\{x_i\}_{i=1}^{n}:x_i\in\mathbb{R}^{p}\}$ is the set of all unordered n-tuple in $\mathbb{R}^{p}$. $q(\cdot) = 1/\pi^{1/(2p)}(\cdot)$ is the charge function and $\lVert \cdot \rVert_{2}$ is the Euclidean distance. Under the proposed charge function, the limiting distribution of the design points converges to $\pi$.
\end{definition}
\noindent However, the optimization problem in \eqref{eq:mined1} is difficult to solve and numerically unstable. To circumvent this issue, \textcite{joseph2019mined} recognize that \eqref{eq:mined1} is closely related to
\begin{equation}
  \label{eq:mined2}
  \arg\min_{\mathcal{D}_n\in\mathbb{D}_n}\bigg[\sum_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\bigg(\frac{q(x_i)q(x_j)}{\lVert x_i - x_j\rVert_{2}}\bigg)^{k}\bigg]^{1/k}
\end{equation}
for $k>0$. As $k\to\infty$, the optimization problem becomes
\begin{equation}
  \label{eq:mined3}
  \arg\min_{\mathcal{D}_n\in\mathbb{D}_n}\max_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{q(x_i)q(x_j)}{\lVert x_i - x_j\rVert_{2}}\; .
\end{equation}
By substituting $q(\cdot) = 1/\pi^{1/(2p)}(\cdot)$ into \eqref{eq:mined3}, we have
\begin{equation}
  \begin{aligned}
    \label{eq:mined4}
    & \arg\min_{\mathcal{D}_n\in\mathbb{D}_n}\max_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{1}{\pi^{1/(2p)}(x_i)\pi^{1/(2p)}(x_j)\lVert x_i - x_j\rVert_{2}} \\
    = & \arg\max_{\mathcal{D}_n\in\mathbb{D}_n}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\pi^{1/(2p)}(x_i)\pi^{1/(2p)}(x_j)\lVert x_i - x_j\rVert_{2} \\
    = & \arg\max_{\mathcal{D}_n\in\mathbb{D}_n}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\gamma^{1/(2p)}(x_i)\gamma^{1/(2p)}(x_j)\lVert x_i - x_j\rVert_{2} \\
    = & \arg\max_{\mathcal{D}_n\in\mathbb{D}_n}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{1}{2p}\log\gamma(x_i) + \frac{1}{2p}\log\gamma(x_j) + \log \lVert x_i - x_j\rVert_{2} \; , \\
  \end{aligned}
\end{equation}
where $\gamma \propto \pi$ is the unnormalized probability density function. Now we only need to work with the log-unnormalized density and therefore, the objective function is numerically more stable. Intuitively, \eqref{eq:mined4} wants the design points to be as far apart as possible while are still placed in the high density regions. If we take $\pi = \mbox{Uniform}[0,1]^{p}$, then \eqref{eq:mined4} reduces to
\begin{equation}
  \label{eq:mined5}
  \arg\max_{\mathcal{D}_n\in\mathbb{D}_n^{u}}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\lVert x_i - x_j\rVert_{2} \; ,
\end{equation}
where $\mathbb{D}_n^{u} = \{\{x_i\}_{i=1}^{n}:x_i\in[0,1]^{p}\}$ is the set of all unordered n-tuple in $[0,1]^{p}$. \eqref{eq:mined5} is the same optimization problem for the maximin design in the unit hypercube \parencite{johnson1990design}. \textcite{joseph2019mined} further propose a generalized distance,
\begin{equation}
  \label{eq:mined6}
  \lVert u\rVert_{s} = \bigg(\frac{1}{p}\sum_{l=1}^{p}|u_{l}|^{s}\bigg)^{1/s},\;  s > 0.
\end{equation}
Under the distance measure defined in \eqref{eq:mined6}, the limiting distribution of the MinED points converge to $\pi$ for all $s > 0$. When $s \to 0$, the distance measure converge to $\lVert u\rVert_{0} = \prod_{l=1}^{p}|u_{l}|^{1/p}$. If $\pi$ is the uniform distribution, MinED with $\lVert \cdot\rVert_{0}$ is the limiting case of the MaxPro design \parencite{joseph2019mined}, showing that noncollapsing property can also be easily achieved by carefully choosing the distance measure.

\subsection{Constrained Minimum Energy Design}
\label{subsec:constrained_minimum_energy_design}
Now consider the case that we need to generate MinED for $\gamma \propto \pi$, an unnormalized probability density function, in some non-rectangular bounded space $\mathcal{X}=\{x\in[0,1]^{p}:g_{k}(x)\leq 0\; \forall k = 1,\ldots,K\}$, then the optimization problem becomes
\begin{equation}
  \label{eq:comined1}
  \arg\max_{\mathcal{D}_n\in\mathbb{D}_n^{\mathcal{X}}}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{1}{2p}\log\gamma(x_i) + \frac{1}{2p}\log\gamma(x_j) + \log \lVert x_i - x_j\rVert_{s} \; ,
\end{equation}
where $\mathbb{D}_n^{\mathcal{X}} = \{\{x_i\}_{i=1}^{n}:x_i\in\mathcal{X}\}$ is the set of all unordered n-tuple in $\mathcal{X}$. However, constraint optimization is generally hard to solve, especially when some of the constraint functions are nonlinear. Similar to Sequentially Constrained Monte Carlo, we can simplify the optimization problem \eqref{eq:comined1} by introducing the probabilistic relaxation $\rho_{\tau}$, \eqref{eq:pc3}, for the inequality constraints $\{g_{k}\}_{k=1}^{K}$, leading to the constrained minimum energy design (CoMinED) defined below.
\begin{definition}{\textbf{(Constrained Minimum Energy Design)}}
\label{df:constrained_minimum_energy_design}
Let $\gamma \propto \pi$ be the target unnormalized probability density function. An n-point minimum energy design of $\pi$ in any non-rectangular bounded space $\mathcal{X}=\{x\in[0,1]^{p}:g_{k}(x)\leq 0\; \forall k = 1,\ldots,K\}$ is the optimal solution of
\begin{equation}
  \label{eq:comined2}
  \arg\max_{\mathcal{D}_n\in\mathbb{D}_n^{u}}\min_{\substack{x_i,x_j\in\mathcal{D}_n\\i\neq j}}\frac{1}{2p}\log\tilde{\gamma}_{\tau}(x_i) + \frac{1}{2p}\log\tilde{\gamma}_{\tau}(x_j) + \log \lVert x_i - x_j\rVert_{s} \; ,
\end{equation}
where $\mathbb{D}^{u}_{n} = \{\{x_i\}_{i=1}^{n}:x_i\in [0,1]^{p}\}$ is the set of all unordered n-tuple in unit hypercube, $\lVert \cdot \rVert_{s}$ is the distance measure function defined in \eqref{eq:mined6}, and
\begin{equation}
  \label{eq:comined3}
  \tilde{\gamma}_{\tau}(\cdot) = \gamma(\cdot) \times \rho_{\tau}(\cdot) = \gamma(\cdot) \prod_{k=1}^{K}\Phi(-\tau g_{k}(\cdot)) \; ,
\end{equation}
where $\tau$ controls the rigidity of the constraints. As $\tau\to\infty$, \eqref{eq:comined2} is equivalent to \eqref{eq:comined1} in the limit, and $\tau = 10^{6}$ is sufficient to achieve the limit numerically in practice, provided that the constraints are properly scaled, a point that we will discuss in detail in Section~\ref{sec:simulation_results}.
\end{definition}

Although the CoMinED is applicable for any distribution $\pi$, in this paper, we mainly focus on $\pi = \mbox{Uniform}[0,1]^{p}$ for a direct comparison to the existing methods for generating space-filling design in non-rectangular bounded regions. As pointed out in MinED \parencite{joseph2015mined,joseph2019mined}, solving the optimization directly using nonlinear programming solver is difficult and computationally expensive. The proposed remedy is to (i) generate the design from a set of candidate samples, and (ii) apply simulated annealing on $\tau$ by starting with an ``easier" problem and slowly increasing the rigidity of the constraints, which is also employed in the Sequentially Constrained Monte Carlo. Suppose we do a $T$ step-simulated annealing, we need to first define the increasing sequence of rigidity parameters $0 = \tau_0 < \tau_1 < \cdots < \tau_T = 10^{6}$. At each step, we generate the $n$-point intermediate CoMinED as follows. Let $\tau_{t}$ be the rigidity parameter, $\mathcal{C}^{t} = \{y_{j}^{t}\}_{j=1}^{N_{t}}$ be the candidate samples, and $\mathcal{D}^{t} = \{x_{i}^{t}\}_{i=1}^{n} \subseteq \mathcal{C}^{t}$ be the CoMinED at the $t$-th step. To construct $\mathcal{D}^{t+1}$, we first augment the candidate samples to $\mathcal{C}^{t+1}$ by including the linear combinations of nearby points in $\mathcal{D}^{t}$ which we call adaptive lattice grid refinement, and then apply the one-point-at-a-time greedy algorithm \eqref{eq:dc2} to solve \eqref{eq:comined2} with $\tau_{t+1}$ as the rigidity parameter. Algorithm~\ref{algo:comined} presents the detail procedures of generating CoMinED.

\begin{algorithm}[t!]
\SetAlgoLined
  \textbf{Design Space:} $\mathcal{X} = \{x\in[0,1]^{p}:g_k(x)\leq 0 \; \forall k = 1,\ldots,K\}$. \\
  \textbf{Initialization:}
  \begin{itemize}
    \setlength{\itemsep}{0pt}
    \setlength{\parskip}{0pt}
    \setlength{\parsep}{0pt}
    \item set the increasing sequence of rigidity parameters $0 = \tau_{0} < \tau_{1} < \cdots < \tau_{T} = 10^{6}$.
    \item set the number of nearest neighbors $Q$ to consider for the candidate augmentation.
    \item simulate $N_{1} > n$ (prime number) lattice points $\{y_j^{(1)}\}_{j=1}^{N_{1}}$ from $[0,1]^{p}$ as the initial set of candidate samples $\mathcal{C}^{1}$.
  \end{itemize}
  \vspace{-\topsep}
  \vspace{2mm}
  \For{$t = 1,\ldots,T$}{
    $\bullet$ \textbf{Construction:} solve \eqref{eq:comined2} with $\tau = \tau_{t}$ by one-point-at-a-time greedy algorithm \eqref{eq:dc2} to obtain the CoMinED $\mathcal{D}^{t} = \{x_{i}^{t}\}_{i=1}^{n}$, i.e., with $\{x_{1}^{t},\ldots,x_{m}^{t}\}$, $x_{m+1}^{t}$ is given by
    \begin{equation}
      \label{eq:comined4}
      \begin{aligned}
      x_{m+1}^{t} = \arg\max_{x\in \mathcal{C}^{t}\backslash\{x_{l}^{t}\}_{l=1}^{m}}\min_{i=1:m}\frac{1}{2p}\sum_{k=1}^{K}\log\Phi(-\tau_{t}g_{k}(x)) & + \\
      \frac{1}{2p}\sum_{k=1}^{K}\log\Phi(-\tau_{t}g_{k}(x_i)) &+ \log \lVert x_i - x_j\rVert_{s} \; .
      \end{aligned}
    \end{equation}
    \If{$t < T$}{
      $\bullet$ \textbf{Augmentation:} augment the set of candidate samples $\mathcal{C}^{t+1} = \mathcal{C}^{t} \cup \tilde{\mathcal{C}}^{t}$ where $\tilde{\mathcal{C}}^{t}$ is the set of linear combinations of nearby points in $\mathcal{D}^{t}$. We construct $\tilde{\mathcal{C}}^{t}$ as follows.
      \For{$i = 1,\ldots,n$}{
        $\bullet$ find the $Q$ nearest neighbors of $x_{i}^{t}$ in $\mathcal{D}^{t}$. \\
        $\bullet$ for each nearest neighbor $\tilde{x}_{i,q}$ ($q = 1,\ldots,Q$), compute the mid-point
        \begin{equation}
          \label{eq:comined5}
            \tilde{y}^{(m)}_{i,q} = x_{i} + \frac{1}{2}(\tilde{x}_{i,q} - x_{i}) = \frac{x_{i}+\tilde{x}_{i,q}}{2} \; ,
        \end{equation}
        and the reflection mid-point
        \begin{equation}
          \label{eq:comined6}
          \tilde{y}^{(r)}_{i,q} = x_{i} - \frac{1}{2}(\tilde{x}_{i,q} - x_{i}) = \frac{3x_{i} - \tilde{x}_{i,q}}{2}\; .
        \end{equation}
        $\bullet$ Update $\tilde{\mathcal{C}}^{t} = \tilde{\mathcal{C}}^{t} \cup \{\tilde{y}^{(m)}_{i,q},\tilde{y}^{(r)}_{i,q}\}$.
      }
      Remove repeated points in $\tilde{\mathcal{C}}^{t}$, and only keep points in $\tilde{\mathcal{C}}^{t}$ that are not in $\mathcal{C}^{t}$.
    }
  }
  \vspace{2mm}
  \textbf{Return:} 1. feasible candidate samples $\{y\in \mathcal{C}^{T}:y\in\mathcal{X}\}$ and 2. the CoMinED $\mathcal{D}^{T}$. \\
 \caption{Algorithm for Generating $n$-point CoMinED.}
 \label{algo:comined}
\end{algorithm}

Now let us discuss the advantage of the proposed adaptive lattice grid refinement over the local maximin LHDs for candidate augmentation. In the MinED algorithm, \textcite{joseph2019mined} augment the candidate points by maximin LHDs in the local region of each MinED point, where the local region is defined as the hypercube inscribed in the ball with center being the MinED point and radius being the distance to its nearest neighbor in the design. The left panel of Figure~\ref{fig:mmlhd_vs_algr} shows the candidate augmentation by the local maximin LHDs. The initial set of points (red circles) are generated using ``Lattice" function in R package \texttt{mined} \parencite{wang2019mined} and maximin LHDs are generated using ``maximinLHS" function in R package \texttt{lhs} \parencite{carnell2019lhs}. We see that (i) some of the augmented candidate samples are arbitrarily close because local regions are overlapping, and (ii) some regions are left unexplored since the proposed local hypercubes cannot fill the space fully. Given the aforementioned shortcoming of the local maximin LHDs, we propose the adaptive lattice grid refinement (ALGR) for a better ``space-filling" candidate augmentation. The ALGR is based on the good rank-1 lattice rule \parencite{nuyens2006lattice,nuyens2007lattice}, one popular type of quasi-Monte Carlo (QMC) methods. Different from other QMC samples such as Halton' and Sobol' points \parencite{niederreiter1992qmc}, lattice points possess a grid structure (see the red circles of Figure~\ref{fig:mmlhd_vs_algr}) that make them advantageous for candidate augmentation. Under the grid structure, all lattice points share the same distance $\delta$ to their nearest neighbors. Augmenting the candidate samples by the mid-point \eqref{eq:comined5} and reflection mid-point \eqref{eq:comined6}, we can ensure that the minimal interpoint spacing of the new candidate samples is $\delta/2$. If we do a $T$-step simulated annealing, then the minimum interpoint distance for the final set of candidates would be $\delta/2^{T}$, which agrees with the minimum interpoint distance constraints used in bridge design \parencite{jones2015bridge}.

The right panel of Figure~\ref{fig:mmlhd_vs_algr} shows the candidate augmentation of the ALGR, which is better space-filling than the local maximin LHDs augmentation. Furthermore, the ALGR exhibits a good trade-off between exploration and exploitation. It starts with a space-filling but sparse grid as the candidate set for good exploration of the whole hypercube. As the rigidity parameter $\tau$ increases, the intermediate CoMinED would only occupy the key regions, leading to refinement of the lattice grid in those regions exclusively for exploitation, and that is why we call it the adaptive lattice grid refinement.

\begin{figure}[t!]
  \centering
  \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/mmlhd.png}
      \caption{Maximin LHDs}
  \end{subfigure}%
  \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/algr.png}
      \caption{ALGR}
  \end{subfigure}%
  \caption{One step candidate augmentation (in green diamonds) on 53 lattice points (in red circles) by 11 maximin LHDs in local regions (left panel) and adaptive lattice grid refinement (ALGR) considering 11 nearest neighbors (right panel).}
  \label{fig:mmlhd_vs_algr}
\end{figure}

\begin{figure}[t!]
  \centering
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/comined_all.png}
      \caption{2{,}155 candidate samples}
    \end{subfigure}%
    \begin{subfigure}{0.33\textwidth}
      \centering
      \includegraphics[width=0.9\textwidth]{figures/motivation/comined_feasible.png}
      \caption{915 feasible samples}
    \end{subfigure}%
  \caption{Left panel shows 2{,}155 candidate samples from applying Algorithm~\ref{algo:comined} to generate $n=53$ points CoMinED on the design space $\mathcal{X}$ defined by \eqref{eq:MOT} with $Q = 5$, $N_{1} = 263$, $T = 8$, and $\{\tau_{t}\}_{t=0}^{8} = [0,e^{1},e^{2},e^{3},e^{4},e^{5},e^{6},e^{7},10^{6}]$. Red circles indicate the initial candidate set of 263 lattice points. Right panel shows the 915 feasible candidate samples.}
  \label{fig:mot-comined}
\end{figure}

Figure~\ref{fig:mot-comined} shows the candidate samples of applying Algorithm~\ref{algo:comined} to generate $n=53$ points CoMinED with $Q=5$ nearest neighbors and $N_{1} = 263$ (the largest prime number that is less than $Qn = 265$) lattice points as the initial candidate set in $T=8$ steps. With 2{,}155 evaluations of the constraints, it yields 915 feasible samples that cover the design space uniformly well: almost no gap spotted visually in the right panel of Figure~\ref{fig:mot-comined}, showing its significant improvement over one-step acceptance/rejection sampling on Latin hypercube samples (Figure~\ref{fig:mot-lhd}) and the adaptive Sequentially Constrained Monte Carlo approach (Figure~\ref{fig:mot-scmc}), which both conduct 2{,}385 evaluations of the constraints to generate the feasible samples. See Figure~\ref{fig:mot-comined_evolution} in Appendix~\ref{appendix:simulation} for the 8-step evolution of the CoMinED.

\section{Simulation Results}
\label{sec:simulation_results}

\begin{table}[t!]
  \centering
  \begin{tabular}{|c|cccc|}
  \hline
  Problem & Dimension ($p$) & No. of LIC & No. of NIC & Feasibility Ratio \\
  \hline
  MOT \eqref{eq:MOT} & 2 & 0 & 3 & 0.0053 \\
  G01 \eqref{eq:G01} & 13 & 9 & 0 & 0.0000 \\
  G04 \eqref{eq:G04} & 5 & 0 & 6 & 0.2696 \\
  G06 \eqref{eq:G06} & 2 & 0 & 2 & 0.0001 \\
  G07 \eqref{eq:G07} & 10 & 3 & 5 & 0.0000 \\
  G08 \eqref{eq:G08} & 2 & 0 & 2 & 0.0086 \\
  G09 \eqref{eq:G09} & 7 & 0 & 4 & 0.0053 \\
  G10 \eqref{eq:G10} & 8 & 3 & 3 & 0.0000 \\
  IBD \eqref{eq:IBD} & 4 & 0 & 3 & 0.0015 \\
  PVD \eqref{eq:PVD} & 4 & 3 & 1 & 0.4032 \\
  SRD \eqref{eq:SRD} & 7 & 0 & 11 & 0.0019 \\
  TSD \eqref{eq:TSD} & 3 & 1 & 3 & 0.0075 \\
  TTD \eqref{eq:TTD} & 2 & 0 & 3 & 0.2179 \\
  WBD \eqref{eq:WBD} & 4 & 1 & 5 & 0.0010 \\
  SCBD \eqref{eq:SCBD} & 10 & 0 & 11 & 0.0005 \\
  \hline
  \end{tabular}
  \caption{Basic information for the constraints of the benchmark problems. LIC stands for linear inequality constraints and NIC stands for nonlinear inequality constraints. Feasibility ratio is estimated using $10^{7}$ Sobol' points in $[0,1]^{p}$ by R package \texttt{randtoolbox} \parencite{christophe2019randtoolbox}.}
  \label{tab:benchmark_info}
\end{table}

In this section, we report the simulation results of applying CoMinED (Algorithm~\ref{algo:comined}) and adaptive SCMC (Algorithm~\ref{algo:scmc}) to 15 benchmark problems with dimensions ranging from 2 to 13 that are popular in the constrained Bayesian optimization literature \parencite{liu2017cgo,dong2018cgo,chaiyotha2020cgo,tao2020cgo}. Since the one-step acceptance/rejection sampling on
Latin hypercube samples is at a clear disadvantage compared to both adaptive SCMC and CoMinED, we do not include the acceptance/rejection sampling simulation results in this section, but are provided in Appendix~\ref{appendix:simulation} for the interested readers. Table~\ref{tab:benchmark_info} provides some basic information of the 15 problems, and their formulations are provided in Appendix~\ref{appendix:benchmark_problems}. All problems are re-scaled to be in the unit hypercube, and the design measures are also compared under the scaling to $[0,1]^{p}$. \par

Throughout the simulations, CoMinED is ran with $s=2$ for the distance measure \eqref{eq:mined6}, the euclidean distance. Also, we take $N_1$, the number of the initial candidate samples, to be the largest prime number that is less than the product of the number of CoMinED points $n$ and the number of neighbors to be considered for candidate augmentation $Q$. For ease of comparison, we fix $T=8$ and the set of rigidity parameters $\{\tau_{t}\}_{t=0}^{8} = [0,e^{1},e^{2},e^{3},e^{4},e^{5},e^{6},e^{7},10^{6}]$, which shows stable performance on both CoMinED and adaptive SCMC. Ways of choosing $Q$ and the rigidity parameters $\tau_{t}$'s are discussed in Appendix~\ref{appendix:comined_details}. Thus, other than the number of design points $n$, $Q$ is the only free parameter that we alter during the simulation for CoMinED. \par

Let us now consider the setting for the adaptive SCMC comparison. Given that we cannot control the total number of constraint evaluations in CoMinED as the repeated samples from augmentation are discarded, we choose the number of samples per iteration in adaptive SCMC to be
\begin{equation}
  \label{eq:scmc_M}
  M = \max\bigg\{nQ, \bigg\lceil\frac{N_{T}}{T+1}\bigg\rceil\bigg\}\;,
\end{equation}
such that (i) $M$ is larger than the number of initial samples of CoMinED ($N_1 < nQ$), and (ii) the total adaptive SCMC samples $M(T+1)$ is larger than the total CoMinED samples $N_{T}$. This puts CoMinED in a slightly \textit{disadvantageous} position in the comparison for assuredly demonstrating its effectiveness over the adaptive SCMC.

Since CoMinED is a deterministic algorithm by the initial candidate of lattice points and the one-point-at-a-time greedy algorithm for designs construction, only one simulation run is performed. For the adaptive SCMC, 50 runs are used for the comparison. Source codes and tutorials can be found at \url{https://github.com/BillHuang01/CoMinED}.

\subsection{Motivation Example}
\label{subsec:motivation}

Let us first consider the two-dimensional motivation example presented in Figure~\ref{fig:mot-space}.
\begin{equation}
  \label{eq:MOT}
    \begin{aligned}
      &g_{1}(x) = x_{1} - \sqrt{50(x_{2}-0.52)^2+2} + 1 \leq 0 \\
      &g_{2}(x) = \sqrt{120(x_{2}-0.48)^2+1) - 0.75 - x_{1}} \leq 0 \\
      &g_{3}(x) = 0.65^2 - x_{1}^2 - x_{2}^2 \leq 0 \\
      \mbox{where} \quad & 0 \leq x_{i} \leq 1\; (i=1,2).
    \end{aligned}
\end{equation}
We compare the performance of CoMinED and adaptive SCMC on generating a 53-point design from their feasible candidate set under the following settings: CoMinED is ran with $Q = 5, 11, 17$, and the corresponding adaptive SCMC comparison is ran with $M$ computed by \eqref{eq:scmc_M}. \par

\begin{figure}[t!]
  \centering
  \begin{subfigure}{0.24\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/moto_feasible_ratio.png}
    \caption{Feasible Ratio}
  \end{subfigure}%
  \begin{subfigure}{0.24\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/moto_fill_dist.png}
    \caption{Fill Distance}
  \end{subfigure}%
  \begin{subfigure}{0.24\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/moto_maximin.png}
    \caption{Maximin}
  \end{subfigure}%
  \begin{subfigure}{0.24\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/moto_maxpro.png}
    \caption{MaxPro}
  \end{subfigure}%
  \caption{Comparisons of the candidates quality and the resulted 53-point design from applying CoMinED (red squares) and adaptive SCMC (violin plots over 50 runs) on the motivation problem \eqref{eq:MOT}. CoMinED is ran with $Q=5,11,17$, and the corresponding adaptive SCMC comparison is ran with $M$ computed by \eqref{eq:scmc_M}. For feasible ratio and maximin measure, the larger the better, and for the fill distance and MaxPro measure, the smaller the better.}
  \label{fig:moto}
\end{figure}

First, let us look at the quality of the feasible samples. Two metrics are considered for the evaluation. One is the \textit{feasible ratio}, the percentage of total samples that are feasible. The larger the feasible ratio, the better the algorithm identifying the feasible region for exploitation. However, the feasible ratio alone could be misleading, as we can always restrict the sampling in an arbitrarily small ball around each feasible point, yielding many samples in the feasible region but also leaving unexplored gaps. For example, Figure~\ref{fig:mot-scmc} shows that 50.5\% ($1{,}205/2{,}385$) of the total adaptive SCMC samples are feasible, but it does not cover the feasible region well as the CoMinED does (Figure~\ref{fig:mot-comined}), which has the feasible ratio of only 42.0\%. Thus, the \textit{fill distance}, the largest distance of any point in $\mathcal{X}$ to the closest feasible samples, is proposed as the other metric to assess how good the algorithm explores the feasible region completely. The smaller the fill distance the better. In simulation, the fill distance is approximated via $10^{4}$ feasible samples from acceptance/rejection sampling on a very large set of Sobol' points in unit hypercube. From Figure~\ref{fig:moto}, we can see that CoMinED exhibits significant improvement in the fill distance over the adaptive SCMC, though CoMinED has smaller feasible ratio. One intuitive explanation is that by the way of candidate augmentation, CoMinED only refines the lattice grid up to certain granularity such that further refinement would not yield candidate samples that add value to the final space-filling design construction, even though those samples are likely feasible. This shows that the CoMinED naturally comes with the heuristic for adaptive resource allocation between exploration and exploitation during candidate augmentation.

After getting the feasible candidate samples, the next step is to construct the desired $n$-point design using the candidate set. In this paper, we consider constructing both maximin and MaxPro designs by the one-point-at-a-time greedy algorithm. Given that the greedy approach likely results in local optimum, we allow for 10 restarts to obtain the design with best measure. For maximin measure, the larger the better \parencite{johnson1990design}; and for the Maxpro measure, the smaller the better \parencite{joseph2015maxpro}. From Figure~\ref{fig:moto}, we can see that CoMinED outperforms adaptive SCMC on both maximin and MaxPro design construction using the candidate set, and the improvement is more significant when the total sample size is small. From both the candidates quality metrics (feasible ratio and fill distance) and the resulted design measures (maximin and MaxPro), additional samples bring minimal benefit for CoMinED as it already cover the feasible space well with $Q = 5$ presented in Figure~\ref{fig:mot-comined}. The robust performance of CoMinED under small sample size makes it the favorable option especially when the constraints are expensive to evaluate. \par

\begin{figure}[t!]
  \centering

  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/mot_feasible_ratio.png}
    \caption{Feasible Ratio (Larger is Better)}
  \end{subfigure}%
  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/mot_fill_dist.png}
    \caption{Fill Distance (Smaller is Better)}
  \end{subfigure}%

  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/mot_maximin.png}
    \caption{Maximin (Larger is Better)}
  \end{subfigure}%
  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/mot_maxpro.png}
    \caption{MaxPro (Smaller is Better)}
  \end{subfigure}%

  \caption{Comparisons of the candidates quality and the resulted 53-point design from applying CoMinED (squares) and adaptive SCMC (violin plots over 50 runs) on the motivation problem (MOT-O, \ref{eq:MOT}) and scaled motivation problem (MOT-S, \ref{eq:MOT-S}). CoMinED is ran with $Q=5,11,17$, and the corresponding adaptive SCMC comparison is ran with $M$ computed by \eqref{eq:scmc_M}. Results without (in blue) and with (in green) constraint value normalization are both presented.}
  \label{fig:mot_all}
\end{figure}

In many real world problems such as the welded beam design \parencite[WBD;][]{dong2018cgo} problem \eqref{eq:WBD}, it is common that the constraints yield values in very different scales, as $g_{3}$ of WBD is in the scale of 10s, but $g_{2}$ is in the scale of 10{,}000s. It is natural to test out how CoMinED and adaptive SCMC perform under the aforementioned circumstance. Consider the scaled version of the motivation example (MOT-S) presented below:
\begin{equation}
  \label{eq:MOT-S}
    \begin{aligned}
      &g_{1}(x) = 10^{-3}(x_{1} - \sqrt{50(x_{2}-0.52)^2+2} + 1) \leq 0 \\
      &g_{2}(x) = \sqrt{120(x_{2}-0.48)^2+1) - 0.75 - x_{1}} \leq 0 \\
      &g_{3}(x) = 10^{3}(0.65^2 - x_{1}^2 - x_{2}^2) \leq 0 \\
      \mbox{where} \quad & 0 \leq x_{i} \leq 1\; (i=1,2).
    \end{aligned}
\end{equation}
The blue squares and violin plots in Figure~\ref{fig:mot_all} are the results from applying CoMinED and adaptive SCMC directly on the motivation problem (MOT-O, \ref{eq:MOT}) and the scaled motivation problem (MOT-S, \ref{eq:MOT-S}). Comparing the MOT-S and the MOT-O facet, we can see that both CoMinED (blue squares) and adaptive SCMC (blue violin plots) perform substantially worse on all evaluation metrics for the scaled problem. Thus, we propose the \textit{constraint value normalization} using the median absolute deviation (MAD) centered at zero for remedy. Suppose that we have evaluated the constraints on some samples $\{x_{i}\}_{i=1}^{N^{\prime}}$, we replace the constraints $g_{k}$ in CoMinED and adaptive SCMC by
\begin{equation}
  \label{eq:cvn}
    \tilde{g}_{k}(\cdot) = g_{k}(\cdot) / \sigma_{k}, \mbox{ where } \sigma_{k} = \mbox{Median}_{i=1:N^{\prime}}(|g_{k}(x_{i})-0|).
\end{equation}
The use of MAD instead of standard deviation is for robustness against large absolute constraint values. The green squares and violin plots in Figure~\ref{fig:mot_all} shows the performance of CoMinED and adaptive SCMC with the constraint value normalization. Compared to the results without the constraint value normalization (blue squares and violin plots), substantial improvements are observed, especially on the scaled motivation example. Hence, for the rest of the simulations, we only compare CoMinED to adaptive SCMC after incorporating constraint value normalization.

\subsection{More Benchmark Problems}
\label{subsec:benchmark_problems}

\begin{figure}[t!]
  \centering
  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/benchmark_feasible_ratio.png}
    \caption{Feasible Ratio (Larger is Better)}
  \end{subfigure}%
  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/benchmark_fill_dist.png}
    \caption{Fill Distance (Smaller is Better)}
  \end{subfigure}%

  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/benchmark_maximin.png}
    \caption{Maximin (Larger is Better)}
  \end{subfigure}%
  \begin{subfigure}{0.45\textwidth}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/experiment/benchmark_maxpro.png}
    \caption{MaxPro (Smaller is Better)}
  \end{subfigure}%

  \caption{Comparisons of the candidates quality and the resulted 109-point design from applying CoMinED (squares) and adaptive SCMC (violin plots over 50 runs) on 14 benchmark problems. The problems are in ascending order by number of dimensions and descending order by the feasibility ratio. CoMinED is ran with $Q = 19$ for problems with dimension smaller than 10 and $Q = 27$ for problems with dimensions at least 10 (SCBD,G07,G01), and the corresponding adaptive SCMC comparison is ran with $M$ computed by \eqref{eq:scmc_M}. For each problem, the evaluation metrics are shifted such that the median of the adaptive SCMC results is 0 and re-scaled such that the IQR of adaptive SCMC results is 1. Both metrics are truncated at $\pm 5$ for visualization purpose.}

  \label{fig:benchmark_all}
\end{figure}

To further demonstrate the improvement of CoMinED, let us look at the simulation results on the 14 benchmark problems with dimensions ranging from 2 to 13, including 7 real world engineering problems. Top panels of Figure~\ref{fig:benchmark_all} compare the quality of the candidate samples. Except for the two-dimensional problems (TTD, G08, G06), CoMinED shows strong improvement in the feasible sample ratio, especially on the high dimensional problems. For example, 23.12\% of the total CoMinED samples are feasible for problem G01, while adaptive SCMC only yields in average 0.78\% feasible samples (See Table~\ref{tab:candidate_quality} in Appendix~\ref{appendix:simulation} for the actual feasible ratio value). On the other hand, looking at the fill distance\footnote{Fill distance is again approximated by $10^{4}$ feasible samples by acceptance/rejection sampling on a larger set of Sobol' points in unit hypercube. However, for problem G07, only $10^{3}$ samples are used since it is too expensive to generate more feasible samples by acceptance/rejection given its extremely small feasibility ratio (less than 1e-6).}, CoMinED outperforms the adaptive SCMC for all benchmark problems except G04, which is a 5 dimensional problem with feasibility ratio of 27.0\% (Table~\ref{tab:benchmark_info}). Given that G04 is an ``easy" problem, with more than 15{,}000 evaluations of the constraints (Table~\ref{tab:candidate_quality} in Appendix~\ref{appendix:simulation}), both CoMinED and adaptive SCMC should yield samples that cover the feasible region well. Next, let us look at the performance of the 109-point designs constructed from the feasible candidate samples. Similar to the motivation examples, we allow for 10 restarts of the one-point-at-a-time greedy algorithm to obtain the best design. From the bottom panels of Figure~\ref{fig:benchmark_all}, we can see that both maximin and MaxPro designs generated using the CoMinED candidates significantly outmatch the corresponding designs by the adaptive SCMC samples. In summary, from the extensive simulation results on the 14 benchmark problems, CoMinED is more robust than the adaptive SCMC for generating good space-filling design candidate points in arbitrary non-rectangular bounded design space, especially when the space is high dimensional and highly constrained.

\section{Conclusion}
\label{sec:conclusion}
This paper proposes the Constrained Minimum Energy Design (CoMinED; Algorithm~\ref{algo:comined}) for constructing space-filling designs in any non-rectangular bounded space defined by inequality constraints. The key idea of the CoMinED is to employ the state-of-the-art deterministic sampling algorithm, Minimum Energy Design (MinED), on the target distribution using the the probabilistic constraints proposed in Sequentially Constrained Monte Carlo (SCMC). Different from the use of local maximin LHDs for candidate augmentation in the MinED algorithm, we propose the adaptive lattice grid refinement that would impose restriction on the minimal interpoint spacing for the candidate samples, making them more favorable as the candidate set for space-filling designs construction. The extensive simulations on the 15 benchmark problems with dimensions ranging from 2 to 13 demonstrate the significant improvement of CoMinED over adaptive SCMC, the best candidate generation approach we can find from the existing literature. CoMinED also enjoys from fewer number of constraint evaluations by avoiding the sampling of the arbitrarily close points that add minimal value for the space-filling designs construction. However, many real world applications involve discrete variables, but CoMinED and adaptive SCMC can only handle continuous variables. One future research direction is to investigate how to construct constrained space-filling design on the set of mixed discrete and continuous variables. On the other hand, as some recent papers have proposed that ``sampling can be faster than optimization" \parencite{ma2019sampling}, another future research direction is to investigate how CoMinED might be useful for solving expensive constrained optimization problem.

\bigskip
\section*{\Large{Acknowledgments}}
This research is supported by a U.S. Army Research Office grant W911NF-17-1-0007.

\bigskip
\printbibliography

\bigskip

\section*{\LARGE{Appendix}}


\section{Introduction}

Let $G$ be a group with some finite generating set $\mathcal G$. We define the metric $d_{\mathcal G}$ on $G$ by taking $d_{\mathcal G} (g_1, g_2)$ to be the infimum over all
$k\ge 0$ such that there exist $f_1,\ldots, f_k\in\mathcal G$ and $\epsilon_1,\ldots,\epsilon_k\in\{-1, 1\}$ satisfying $g_2=f_1^{\epsilon_1}\cdots f_k^{\epsilon_k}g_1$.

Now let $H$ be an arbitrary group. An element $f\in H$ is called {\it distorted} in $H$  if there exists
a finitely generated subgroup $G\subset H$ containing $f$ such that
$$
\lim_{n\to\infty}\frac{d_{\mathcal G}(f^n, {\rm id})}{n}=0
$$
for some (and hence for every) generating set $\mathcal G$. Since the limit always exists, it is enough to verify it for some subsequence.
The notion of distortion comes from geometric group theory and it was introduced by M. Gromov in \cite{Gromov}.
\vskip2mm

The problem of the existence of distorted elements in some groups of homeomorphisms has been intensively studied for many years (see \cite{C_F, Franks, F_H, G_L, Navas1, Pol}). Substantial progress has been achieved for groups of diffeomorphisms of manifolds. In particular, A. Avila \cite{Avila} proved that rotations with irrational rotation number are distorted in the group of smooth diffeomorphisms of the circle.  In this note we give a constructive proof that all irrational rotations are distorted both in the group of piecewise-affine circle homeomorphisms, ${\rm PAff}_+({\bf S}^1)$, and in the group of smooth circle diffeomorphisms, ${\Diff}^{\infty}(\mathbb R/\mathbb Z)$.
The result gives an answer to {\bf Question 11} in \cite{Navas} (see also {\bf Question 2.5} in \cite{Franks}). So far it has not even been known whether there exist distorted elements in ${\rm PAff}_+({\bf S}^1)$. Now from \cite{G_L} it follows that distorted elements, if they exist, are conjugate to rotations.

\vskip2mm
From now on let $G$ be either the group ${\rm PAff}_+({\bf S}^1)$ or
 ${\Diff}^{\infty}(\mathbb R/\mathbb Z)$. We say that $g\in G$ is {\it trivial on some set} if there exists a non-empty open set $I\subset {\bf S}^1$  such that $g(x)=x$ for $x\in I$. The set of all homemomorphisms in $G$ which are trivial on some set will be denoted by $G_{{\rm {triv}}}$. By ${\rm T}$ we denote the set of all rotations.

 \vskip2mm
This note is devoted to the proof of

\begin{thm}
All irrational rotations are distorted in $G$.
\end{thm}

\section{Proofs}

We first formulate two lemmas and deduce the theorem.
The proofs of the lemmas will be given at the end of the paper.

\begin{lem}\label{L1}
For any irrational  rotation $T_{\alpha}$ and $g\in G_{{\rm {triv}}}\cup {\rm T}$  there exist a finite generating set $\mathcal G_g\subset G$ and a constant $C>0$ such that
\[
d_{\mathcal G}(T_{\alpha}^n g T_{\alpha}^{-n}, {\rm id})\le C\log n\quad\text{for all $n\ge 1$}.
\]
\end{lem}

\vskip2mm

\begin{lem}\label{Prop1_18.04.21}  In $G$ there exist
$g_1,\ldots, g_l\in G_{{\rm {triv}}}\cup{\rm T}$ and $k, k_1,\ldots, k_l\in {\mathbb Z}$
with $k\neq k_1+\cdots+k_l$, such that for each sufficiently small $\beta>0$ the element
$x=T_{\beta}$ satisfies
\begin{equation}\label{e1_17.04.21}
x^{k_1}g_1x^{k_2}g_2\cdots x^{k_l}g_l=x^k.
\end{equation}
\end{lem}

{\it Proof of the theorem.}  Fix an irrational rotation $T_{\alpha}$. From Lemma \ref{Prop1_18.04.21} it follows that in $G$ there exists an equation of the form (\ref{e1_17.04.21}) such that $x=T_{\beta}$, for all sufficiently small $\beta$, is its solution. Let $\mathcal G=\mathcal G_{g_1}\cup\cdots\cup \mathcal G_{g_l}$, where $\mathcal G_{g_i}$, $i=1,\ldots, l$, are finite generating sets derived from Lemma \ref{L1} for $T_{\alpha}$.
We may rewrite equation (\ref{e1_17.04.21}) in the form
\begin{equation}\label{e2_18.04.21}
x^{k_1}g_1x^{-k_1} x^{k_2+k_1}g_2x^{-k_2-k_1}\cdots x^{k_1+\cdots+k_l}g_l x^{-k_1-\cdots-k_l}=x^{k-k_1-\cdots-k_l}.
\end{equation}
Let $\beta_0$ be a positive constant such that $x=T_{\beta}$ for $\beta\in (0, \beta_0)$ satisfies (\ref{e2_18.04.21}). Set $m:=k-k_1-\cdots-k_l$, and let $(n_i)$ be an increasing sequence of  integers such that $n_i\alpha \in (0, \beta_0)\, ({\rm mod}\, 1)$. From Lemma \ref{L1} it follows that
$$
d_{\mathcal G}(T_{\alpha}^{n_i(k_1+\cdots+k_j)} g_j T_{\alpha}^{-n_i(k_1+\cdots+k_j)}, {\rm id})\le C_j\log n_i\quad\text{for all $i\ge 1$ and $j=1,\ldots, l$}.
$$
Since $x=T_{n_i\alpha}$ satisfies (\ref{e2_18.04.21}), we obtain
$$
d_{\mathcal G}(T_{\alpha}^{n_i m}, {\rm id})\le \sum_{j=1}^l C_j\log n_i:= C \log n_i\quad\text{for all $i\ge 1$.}
$$
Hence
$$
\lim_{n\to\infty}\frac{d_{\mathcal G}(T_{\alpha}^n, {\rm id})}{n}=
\lim_{i\to\infty}\frac{d_{\mathcal G}(T_{\alpha}^{n_i m}, {\rm id})}{m n_i}\le \lim_{i\to\infty}\frac{C \log n_i}{m n_i}=0
$$
and the proof is complete. $\square$
\vskip4mm

{\it Proof of Lemma 1.} The proof relies on the observation that for a given interval $I\subset (0, 1)$  there exists a finite generating set $\mathcal G\subset G$ such that for any $n\ge 1$ there exists a
homemorphism $h_n$ with $d_{\mathcal G}(h_n, {\rm id})\le C\log n$ for some constant $C>0$ independent of $n$, and $h_n(x)=T_{\alpha}^n(x)$ for $x\notin I$.
Without loss of generality we may assume that $I=(a, 1).$ Let $m\ge 1$ be an integer such that $a+2/m<1$. Let $h\in G $ be any homeomorphism such that $h(x)=x/2$ for $x\in [0, a+2/m)$, and let $r(x)=x+1/m$.

We shall define $h_n$ by induction. Set $h_0=h$. If $n$ is odd we put $h_n=T_{\alpha} h_{n-1}$. If $n$ is even, we take $s_n:=h_{n/2} h$
and observe that $s_n((0, a))=(n\alpha/2, a/2+n\alpha/2)$. Let $k\in\{1,\ldots, m\}$ be such that $n\alpha/2+k/m\in [0, 1/m)$\,(mod $1$). Then $r^k s_n((0, a))\subset (0, a/2+1/m)$. Therefore
$$
h^{-1} r^k h_{n/2} h(x)=2(x/2+n\alpha/2+k/m)=x+n\alpha+2k/m=T_{\alpha}^n(x)+2k/m
$$
for $x\in (0, a)$. Put $h_n:=r^{-2k}h^{-1} r^k h_{n/2} h$, and let $\mathcal G:=\{T_{\alpha}, h, r\}$. Note that
$$
d_{\mathcal G}(h_n, {\rm id})\le 3m+3+d_{\mathcal G}(h_{\lfloor n/2\rfloor}, {\rm id}).
$$
Thus we obtain $d_{\mathcal G}(h_n, {\rm id})\le C\log n$. Finally, observe that
for any $g\in G_{{\rm {triv}}}$ such that $g(x)={\rm id}$ on $I$ we have
$$
T_{\alpha}^n g T_{\alpha}^{-n}=h_n g h_n^{-1}.
$$
Therefore, we obtain
$$
d_{\mathcal G}(T_{\alpha}^n g T_{\alpha}^{-n}, {\rm id})\le C\log n.
$$
In the case when $g$ is a rotation the conclusion of the lemma is obvious. $\square$
\vskip4mm
{\it Proof of Lemma 2.} Let $\beta\in (0, 10^{-3})$, and let $f_1\in G_{{\rm {triv}}}$ be arbitrary such that
$$
f_1(x)=0.4+2(x-0.4)\,\,\text{for $x\in [0.4, 0.6]$ and } f_1(x)=x\,\,\text{for $x\in [0.9, 1.1]$}.
$$
Set
$$
H_1=T_{2\beta}^{-1} f_1 T_{2\beta} f_1^{-1}.
$$
It is obvious that
$$
H_1(x)=x+2\beta\,\,\text{for $x\in [0.41, 0.79]$  and }  H_1(x)=x\,\,\text{for $x\in [0.91, 1.09]$}.
$$
Define
$$
H_2=T_{1/2} H_1^{-1} T_{1/2} H_1,
$$
and observe that
$$
H_2(x)=x-2\beta\quad\text{for $x\in [0.95, 1]$}.
$$
Simple computation gives
$$
T_{1/2} H_2T_{1/2} H_2={\rm id}.
$$
Set
$$
H_3 =T_{2\beta} H_2.
$$
Then we have
$$
H_3(x)=x\quad\text{for $x\in [0.95, 1]$}
$$
and
\begin{equation}\label{e1_22.04.21}
T_{2\beta+1/2}H_3 T_{-2\beta-1/2} H_3=T_{4\beta}.
\end{equation}
Take an arbitrary $f_2\in G_{{\rm {triv}}}$ satisfying
$$
f_2(x)=2x\quad\text{for $x\in [0, 0.49]$},
$$
and define
$$
H_4=f_2^{-1} H_3f_2.
$$
It is easy to see that
$$
H_4(x)=
\begin{cases}
H_3(2x)/2 & \text{for }  x\in [0, 1/2),\\
x & \text{for }   x\in [1/2, 1).
\end{cases}
$$
Let
\begin{equation}\label{e2_22.04.21}
H_5=T_{1/2} H_4 T_{1/2}H_4.
\end{equation}
Observe that the graph of $H_5$ is built
from two scaled copies of $H_3$, i.e.
$$
H_5(x)=
\begin{cases}
H_3(2x)/2 & \text{for }  x\in [0, 1/2),\\
H_3(2x-1)/2+1/2 & \text{for }  x\in [1/2, 1).
\end{cases}
$$
Therefore, by (\ref{e1_22.04.21}) and (\ref{e2_22.04.21}), we finally obtain
\begin{equation}\label{e1_18.04.21}
T_{\beta+1/4}H_5 T_{-\beta-1/4}H_5=T_{2\beta}.
\end{equation}
Indeed, this is easy to see if we realize that (\ref{e1_18.04.21}) is simply equation (\ref{e1_22.04.21}) rewritten in the new coordinates $(x/2, y/2)$.
Plugging subsequently $H_5, H_4, H_3, H_2$ and $H_1$ into formula (\ref{e1_18.04.21}) we obtain the required equation with $g_1,\ldots, g_l\in\{f_1, f_2, f_1^{-1}, f_2^{-1}, T_{1/2}, T_{-1/2}, T_{1/4}, T_{-1/4}\}\subset G_{{\rm {triv}}}\cup {\rm T}$. Moreover, on the left hand side of this equation
there are eight more $T_{\beta}$'s than $T^{-1}_{\beta}$'s. On the other hand, on the right hand side we have  two $T_{\beta}$'s. Since $\beta\in (0, 10^{-3})$ was arbitrary and the functions $g_1,\ldots, g_l$ were independent of $\beta$, the proof of the lemma is complete. $\square$
\vskip2mm