Constructing two-level $Q_B$-optimal screening designs using mixed-integer programming and heuristic algorithms

Abstract

Two-level screening designs are widely applied in manufacturing industry to identify influential factors of a system. These designs have each factor at two levels and are traditionally constructed using standard algorithms, which rely on a pre-specified linear model. Since the assumed model may depart from the truth, two-level $Q_B$-optimal designs have been developed to provide efficient parameter estimates for several potential models. These designs also have an overarching goal that models that are more likely to be the best for explaining the data are estimated more efficiently than the rest. However, there is no effective algorithm for constructing them. This article proposes two methods: a mixed-integer programming algorithm that guarantees convergence to the two-level $Q_B$-optimal designs; and, a heuristic algorithm that employs a novel formula to find good designs in short computing times. Using numerical experiments, we show that our mixed-integer programming algorithm is attractive to find small optimal designs, and our heuristic algorithm is the most computationally-effective approach to construct both small and large designs, when compared to benchmark heuristic algorithms.

Appendices

Appendix A: Proofs

Proof of Lemma 1

Let $\textbf{D} = (D_{ij})$ where $D_{ij} = \pm 1$. We define

$$\begin{aligned} s_{i_1 i_2 \ldots i_k} = \sum _{l = 1}^{n} D_{l i_1} D_{l i_2} \cdots D_{l i_k}, \end{aligned}$$

where $i_r \in \{1, \ldots , m\}$. Note that when all elements in $(i_1, i_2, \ldots , i_k)$ are different, $s_{i_1 i_2 \ldots i_k}$ is a $J_k$-characteristic of $\textbf{D}$ (Tang 2001). The generalized word count of length k of $\textbf{D}$ can be calculated as follows:

$$\begin{aligned} B_k = \frac{1}{n^2 } \sum _{i_1< \cdots < i_k} s^{2}_{i_1 \ldots i_k}. \end{aligned}$$

From Lemma 1 of Butler (2003b), we have that

$$\begin{aligned} E_k = \frac{1}{n^2 } \sum ^{m}_{i_1 = 1} \sum ^{m}_{i_2 = 1} \cdots \sum ^{m}_{i_k = 1} s^{2}_{i_1 i_2 \ldots i_k}. \end{aligned}$$

(8)

We will show that the power moments $E_1$, $E_2$, $E_3$ and $E_4$, can be expressed as linear combinations of the generalized word counts $B_1$, $B_2$, $B_3$ and $B_4$. The first power moment is obvious since $E_1 = B_1$. For the second power moment, we have that

$$\begin{aligned} E_2 = \frac{1}{n^2} \sum _{i_1 = 1}^{m} s^2_{i_1 i_1} + \frac{2}{n^2} \sum _{i_1 = 1}^{m-1} \sum _{i_2 = i_1 + 1}^{m} s^{2}_{i_1 i_2} = m + 2 B_2. \end{aligned}$$

The third power moment $E_3$ in Equation (8) can expressed as a sum of $n^{-2} s^2_{i_1, i_2, i_3}$ over three types of triplets $(i_1, i_2, i_3)$. The first type involves triplets in which all elements are equal. The sum of all $n^{-2} s^2_{i_1, i_2, i_3}$ over this type equals $B_1$. The second type involves triplets in which all elements are distinct. The sum of all $n^{-2} s^2_{i_1, i_2, i_3}$ over the second type equals $6 B_3$, where the factor 6 is due to the fact that there are 6 permutations of the elements in $(i_1, i_2, i_3)$. The third type of triplets is such that exactly two elements are equal. The sum of all $n^{-2} s^2_{i_1, i_2, i_3}$ over this type is $3(m-1)B_1$. Therefore, we have that $E_3 = (3\,m - 2)B_1 + 6B_3$.

The fourth power moment $E_4$ in Equation (8) can expressed as a sum of $n^{-2} s^2_{i_1, i_2, i_3, i_4}$ over five types of quadruplets $(i_1, i_2, i_3, i_4)$. These types are summarized as follows: (1) All elements are equal; (2) all elements are distinct; (3) exactly two elements are equal; (4) exactly three elements are equal; and, (5) each element is equal to exactly one other element. By a counting argument, we can show that the sum of the $n^{-2} s^2_{i_1, i_2, i_3, i_4}$ over the first, second, third, fourth and fifth type of quadruplets equals m, $24 B_4$, $12(m-2)B_2$, $8 B_2$ and $3\,m(m-1)$, respectively. The sum of these totals shows that $E_4 = 24B_4 + 4(3m - 4) B_2 + m(3m-2)$. $\square$

Table 5 Average computing times in seconds and their standard deviations for 10 optimizations performed by the heuristic algorithms

Proof of Corollary 1

Let $\textbf{d}_j$ denote the $m \times 1$ vector involving the j-th row of $\textbf{D}$. Without loss of generality, $\textbf{D} = [\textbf{d}_j, \textbf{D}^{T}_{-j}]^{T}$ where $\textbf{D}_{-j}$ is the $(n-1) \times m$ design matrix excluding the j-th row. We then have that

$$\begin{aligned} \textbf{T} = \left( \begin{array}{cc} m &{} \textbf{d}_{j}^{T}\textbf{D}^{T}_{-j} \\ \textbf{D}_{-j} \textbf{d}_j &{} \textbf{D}_{-j} \textbf{D}_{-j}^{T} \end{array}\right) . \end{aligned}$$

Let $\textbf{r} = \textbf{D}_{-j} \textbf{d}_j$ and $\textbf{S} = \textbf{D}_{-j} \textbf{D}_{-j}^{T}$. The k-th power moment of $\textbf{T}$ is

$$\begin{aligned} n^2 E_k = m^k + 2 \sum _{u=1}^{n-1} r^{k}_{u} + \sum _{u=1}^{n-1} \sum _{v=1}^{n-1} S^{k}_{uv}, \end{aligned}$$

(9)

where $r_u$ is the u-th element of $\textbf{r}$ and $S_{uv}$ is the element in the u-th row and the v-th column of $\textbf{E}$. For Case II, substituting the first four power moments given by Equation (9) in the $Q_B$ criterion in Theorem 1 gives

$$\begin{aligned} \frac{1}{n}\left\{ w_0 + \sum _{k=1}^{4} \frac{w_k}{n^2} \left( m^k + 2 \sum _{u=1}^{n-1} r^{k}_{u} \right) + \sum _{k=1}^{4} \frac{w_k}{n^2} \left( \sum _{u=1}^{n-1} \sum _{v=1}^{n-1} S^{k}_{uv}\right) \right\} . \end{aligned}$$

The contribution of the j-th row of $\textbf{D}$ to its $Q_B$ criterion value is the second term in the expression above, where we have that $\sum _{u=1}^{n-1} r^{k}_{u} = \sum _{i\ne j} T^{k}_{ij}$. The expression for Case I is found similarly. $\square$

Appendix B: Computing times

Our goal here is to compare the computing times required for a completed optimization by the heuristic algorithms. An optimization of the PBCE algorithm involves its standard settings, that is, a perturbation size ($\alpha$) of 0.1, a maximum number of perturbation without improvement (M) equal to 100, and 5 restarts of the whole procedure. An optimization of the CE, RCP and PE algorithms involves 1000 iterations. In this way, we mimic executions of the algorithms conducted by a standard user.

We consider three sets of design problems. The first set has all combinations of numbers of runs and numbers of factors in Table 1 with $\pi _1$ = 0.625, while the second set has the 7-factor designs in Table 4. The third set of problems involves the 11-factor designs with 20 and 24 runs in Table 4.

Table 5 shows the computing times required for a completed optimization by the heuristic algorithms. For each design problem, the table gives the averages and standard deviations of 10 optimizations. Preliminary experiments revealed that a single optimization of the PE algorithm with 1000 iterations is computationally more demanding than the other algorithms, especially for designs with more than nine factors. For this reason, Table 5 gives the computing times required for optimizations of the PE algorithm with 10 iterations only. Using these computing times, we estimate the average time required by an optimization of the PE algorithm with 1000 iterations.

Clearly, the PBCE algorithm outperforms the CE and RCP algorithms in terms of the computing time. For the 9-, 11-, 13- and 17-factor designs in Table 5, the PBCE algorithm even takes less computing time than the PE algorithm with 10 iterations. For the 17-factor 18-run design, the PE algorithm was computationally infeasible because its 10 iterations took longer than an hour. In contrast, the PBCE algorithm required 15 s to solve the optimization problem in this case.

For the 6- and 7-factor designs, the average computing times of the PE algorithm in Table 5 should be multiplied by 100, so as to obtain the times required by optimizations with 1000 iterations. Therefore, the computing times for the PE algorithm are around 9.3 s for six factors, and between 227.9 and 532.9 s for seven factors. Clearly, these computing times are much larger than those of the PBCE algorithm in the table.

In conclusion, a completed optimization of the PBCE algorithm is computationally less expensive than one of the benchmark heuristic algorithms.