During the past decades, order-of-addition (OofA) experiments are widely applied to biochemistry (Olsen et al. 1994), foot industries (Karim, Mccormick, and Kappagoda 2000), as well as chemical-related areas (Jourdain et al. 2009). In those mentioned experiments, several components are added into an apparatus sequentially rather than simultaneously, and different addition sequences of reactants yield different responses. However, as the number of components increases, the sample space suffers from a combinatorial explosion, and performing full design becomes unfeasible. Accordingly, our goal is to model for an OofA experiment and construct an informative and economical fractional design.

For the objective of optimizing and predicting the response, several statistical models are created for the OofA experiment: Van Nostrand (1995), Voelkel (2019), Peng, Mukerjee, and Lin (2019) proposed a PWO model with diminishing effects for pairwise factors; Mee (2020) added PWO factor interactions to account for sequencing effects not accounted for by pairwise main effects alone; Yang, Sun, and Xu (2020) proposed a component-position (CP) model for order-of-addition using categorical explanatory variables. In this article, we only focus on the first-order PWO model. However, the PWO model is also a regression model. Then, a family of criteria can be applied to find optimal designs under the PWO model, such as D-, A- and $M.S.$-optimal designs. The optimality proof indicates that a full PWO design with $m!$ distinct permutations of components is D-, A-, and $M.S.$-optimal, but the run size is extremely large. Taking $m=10$ as an example, there are $m!=3628800$ distinct permutations. Consequently, over three million runs of experiments should be implemented, which is impractical. Therefore, fractions of full PWO designs with a smaller number of runs are preferable.

Recently, four kinds of fractional PWO designs have been studied. Peng, Mukerjee, and Lin (2019) introduced a method for constructing optimal PWO designs. This method limits the run size to $m!/r!(2\le r\le m/2)$, which is also too many for experimenters to afford. For instance, if $m=10$, the method needs at least 30240 runs to be implemented. Yang, Sun, and Xu (2021) and Zhao, Lin, and Liu (2022) provided construction methods based on an orthogonal array, the resulting designs are component orthogonal arrays (COAs) and OofA orthogonal arrays (OofA-OAs), in which the run size is also inflexible. Zhao, Lin, and Liu (2021) provided a minimal-point design with $m(m-1)/2+1$ runs. The run size is small, but the efficiency is relatively low. However, theoretical constructions of these fractional PWO designs are highly dependent on run size. Winker, Chen, and Lin (2020) applied the threshold-accepting algorithm to construct the optimal designs (D-efficiency for application) based on the pairwise-order (PWO) model and the tapered PWO model, the designs obtained by threshold-accepting algorithm for $4\le m\le 30$ with $n=m(m-1)/2+1,m(m-1)+1,3m(m-1)/2+1,$ respectively, are provided for practical uses.

This article gives serious consideration to constructing D-, A-, and $M.S.$-optimal PWO designs. The main contributions of this article are as follows. First, we construct a new hybrid algorithm for generating the PWO design, which is flexible with no restriction on run size. Second, D-, A-, and $M.S.$-optimal PWO designs can be constructed using the proposed algorithm. Third, the constructed designs possess high efficiencies compared with the full PWO design, though with a small fraction of runs. Some optimal PWO designs and fully efficient PWO designs are listed in the tables.

This paper is organized as follows. We first introduce the PWO model in Section 2. Section 3 gives a review of Fedorov’s exchange algorithm for constructing the D-optimal designs. Then, this algorithm is modified and extended for constructing A- and $M.S.$-optimal designs. Some theoretical results on the information matrix and algorithm are also provided in Section 3. In Section 4, based on the exchange algorithm and particle swarm optimization (PSO) algorithm, a novel hybrid algorithm is proposed to achieve D-, A- and $M.S.$-optimal PWO designs. Some numerical results are given in Section 5. Finally, concluding remarks are provided in Section 6.

Now, we introduce the Van Nostrand PWO model. Suppose there are m components denoted as $1,\dots ,m$. Any treatment in the OofA experiment corresponds to a permutation of $1,\dots ,m$, denoted as $\alpha$, and the first-order PWO model can be expressed as follow where each $z_{jk}(\alpha )$ is a PWO indicator between j and k, For an n-point PWO design, let Y be the n-dimensional response vector, Z be the design matrix with $\left(\genfrac{}{}{0pt}{}{m}{2}\right)$ columns corresponding to PWO indicators $z_{12},z_{13},\dots ,z_{(m-1)m}$, and $\beta =({\beta }_{12},{\beta }_{13},\dots ,$ ${\beta }_{(m-1)m}{)}^{\prime }$, where ${}^{\prime }$ denotes the transpose. Then, the first-order PWO model can be written as or where $X={(\mathbf{1}Z)}_{n\times p}$ with $p=\left(\genfrac{}{}{0pt}{}{m}{2}\right)+1$ is the model matrix, $\tilde{\beta }={({\beta }_0,{\beta }^{\prime })}^{\prime }$ represents the parameter of interest, and the random error $ϵ\sim N(\mathbf{0},{\sigma }^{2}I)$. Mee (2020) extended the PWO model to the high-order case. Here, we only consider a first-order PWO model. The proposed algorithms also apply to a higher-order PWO model.

Furthermore, we refer to $\tilde{M}=X^{\prime }X/n$ as the information matrix of an n-point PWO design. Under the PWO model (3), the variance-covariance matrix of the least squares estimator of $\beta$ is proportional to $\tilde{M}$. Hence, it is desirable to maximize the matrix $\tilde{M}$ under some criteria. The popular criteria include the D-criterion $det{(\tilde{M})}^{1/p}$, the A-criterion $tr({\tilde{M}}^{-1})$, the $M.S.$-criterion $tr({\tilde{M}}^2)$ (see the reference Atwood 1969). Note that $tr({\tilde{M}}^{-1})$ is interpreted as $+\infty$ for singular $X^{\prime }X$. Let $X_f$ be the full PWO design and the corresponding information matrix be ${\tilde{M}}_f=X_f^{\prime }{X}_f/n$. For clarity, we take $m=3$ as an example to illustrate the characteristics of the full PWO design under D-, A- and $M.S.$-criteria. The levels of PWO factors in the full PWO design with 3 components are as follows.

From this, we obtain and $det({\tilde{M}}_f)=16/27,$ $tr({\tilde{M}}_f^{-1})=11/2$ and $tr({\tilde{M}}_f^2)=14/3$.

Theoretical constructions on optimal designs are always complicated; hence, computer algorithms are applied for constructing approximate and exact optimal designs in the literature. The exchange algorithm is one of the popular computer algorithms for constructing optimal designs for the cases with the design points being selected from a finite design space. Fedorov (1972) first proposed an exchange algorithm for generating D-optimal designs. This algorithm chooses n points to include in the design from a finite set of possible points called candidate points. It starts with nonsingular n-point designs and then adds and deletes one observation to achieve increases in the determinant. After that, some improved implementations are proposed based upon Fedorov’s exchange algorithm, such as the Kiefer round-off algorithm, the Mitchell algorithm, the Wan Schalkwyk algorithm, the combined Fedorov, the Wynn–Mitchell algorithm, and so on; see the references Mitchell (1992), Nguyen and Miller (1992).

Before introducing the new algorithm, we add some details of the PSO algorithm. PSO is a population-based stochastic algorithm for optimization. Each population member is described as a particle that moves around a search space testing new criterion values. All particles survive from the beginning of a trial until the end, and their interactions result in iterative improvement of the quality of the problem solutions over time. The most common type of implementation defines the particles’ behavior as adjusting toward each of its personal best position (local-best) and global-best position so that its trajectory shifts to new regions of the search space and the particles gradually cluster around the optima. To generate optimal experimental designs, a particularly challenging task is to redefine the particle designs’ movement toward its personal local-best design and global-best design. A review of some recent applications of PSO and its variants to tackle various types of efficient experimental design is Chen, Chen, and Wang (2022). Since finding optimal PWO designs for the OofA experiment is to solve a discrete optimization problem, we utilize an update procedure for the particle designs that is similar to the modified PSO algorithms in Chen et al. (2014) and Phoa et al. (2016). Each particle design relates to its personal local-best design which is derived by exchange procedures starting from itself. During each iteration, the current particle design is adjusted toward its personal local-best design as well as the global-best design by exchanging points with each other.

Now, we introduce a new hybrid algorithm called Ex-PSO algorithm, which combines the single-point exchange algorithm and PSO algorithm for generating D-, A-, and $M.S.$-optimal designs. The single-point exchange algorithm is used for generating the local-best design concerning each particle design. The PSO algorithm ensures the particle designs gradually cluster around the optimal PWO design. To avoid singularity, the technique proposed in Remark 2 is used; hence, a parameter $\theta$ with a small value is involved in this algorithm.

Since the Ex-PSO algorithm involves a set of parameters denoted as $s,t_{ex},t_{pso},\theta ,c_1,c_2$, we also refer to it as Ex-PSO $(m,n;s,t_{ex},t_{pso},\theta ,c_1,c_2)$ for generating optimal PWO design with m components and n runs. For clarity, we create a programming chart to illustrate the steps of Ex-PSO $(m,n;s,t_{ex},t_{pso},\theta ,c_1,c_2)$. Further, we explain the optimization process and the uses of these parameters as follows. Denote ${\mathcal{L}}_k$s and $\mathcal{G}$ as the local-best designs and the global-best design, respectively. These designs are updated during each iteration of the Ex-PSO algorithm. Each local-best design ${\mathcal{L}}_k$ is derived from the current particle design $D_k$ via a fixed number of iterations of the single-point exchange procedure, denoted as $t_{ex}$. In addition, the global-best design $\mathcal{G}$ is the optimal local-best design that maximizes $\varphi ({\mathcal{L}}_k)$. And the number of iterations of the PSO algorithm is denoted as $t_{pso}$. Meanwhile, two parameters are used to control the PSO behavior of the Ex-PSO algorithm: $c_1$ and $c_2$, which account for the velocities at which each current design drifts toward the corresponding local-best and global-best design. More specifically, during each iteration of the PSO algorithm, we randomly exchange $c_1$ points from the difference set $D_k\setminus {\mathcal{L}}_k$ with $c_1$ points from ${\mathcal{L}}_k\setminus D_k$ and then randomly exchange $c_2$ points from the difference set $D_k\setminus \mathcal{G}$ with $c_2$ points from $\mathcal{G}\setminus D_k$. This procedure corresponds to the “Update $D_k$ by PSO” box in the programming chart.

Finally, we note that the time complexity for each iteration is $O(n^3)$. However, the search procedure of PSO also relates to the full design which is a combinatorial large set with $m!$ points, so Ex-PSO algorithm may be time-consuming for large m. The Ex-PSO algorithm implemented in MATLAB running on Intel(R) Core(TM) i7-8550U GHz with 8 GB Memory. Take the case of $m=6,n=16$ for example, it takes 30.42 seconds for running Ex-PSO $(6,16;10,20,100,0.005,1,1)$.

In this section, we illustrate the performances of the obtained designs constructed by the Ex-PSO algorithm. For brevity, the generated designs are denoted as Ex-PSO-D, Ex-PSO-A and Ex-PSO-M.S. designs that respectively correspond to the objective functions (4)–(6) which are considered in the exchange algorithm. Numerical simulations show that these designs are powerful for fitting PWO models in terms of the D-, A- and $M.S.$-efficiencies. The efficiencies are derived from comparison with the full PWO design since the information matrix of the full PWO design has been proven to be universally optimal. Therefore, we have the D-, A- and $M.S.$-efficiencies that calculate $\frac{det{(\tilde{M})}^{1/p}}{det{({\tilde{M}}_f)}^{1/p}}$, $\frac{tr({\tilde{M}}_f^{-1})}{tr({\tilde{M}}^{-1})}$ and $\frac{tr({\tilde{M}}_f^2)}{tr({\tilde{M}}^2)}$, respectively, where $\tilde{M}$ and ${\tilde{M}}_f$ are the information matrices of the obtained design and full PWO design respectively, and p is the number of the columns of the model matrix X.

The number of PSO particles (s), the maximum iteration counts of the single-point exchange algorithm and PSO algorithm ($t_{ex},t_{pso}$), and the numbers of pairs of exchanging points with which each particle design drifts toward the local-best and global-best design $(c_1,c_2)$, control the optimization process of Ex-PSO algorithm. It seems reasonable that these parameters should be increased for large n and m. In our test of searching for optimal PWO designs with $m=4,5,6,7$ components, we recommend that these parameters set at $s=10$, $t_{ex}=20$, $t_{pso}=100$, $c_1=c_2=1$. Furthermore, we recommend setting the maximum iteration counts of the exchange algorithm and PSO at $t_{ex}=20$ and $t_{pso}=100$, respectively, which achieves high computational efficiency. Further, to demonstrate the performance of such a set of parameters, we randomly run the algorithm Ex-PSO $(m,n;10,20,100,0.005,1,1)$ one hundred times for generating the Ex-PSO-A designs, because the exchange algorithm seems inefficient under A-optimal criterion, as shown in Table 2. Therefore, one hundred Ex-PSO-A designs with $m=4,5,6,7$ components and $n=m(m-1)$ runs are generated, and Figure 3 highlights that all Ex-PSO-A designs reach at least $93\%$ of the efficiency of the full PWO design. For the cases with large m, the settings on maximum iterations, $t_{ex}$ and $t_{pso}$ may not be enough, but the Ex-PSO algorithm still returns approximate optimal PWO designs; see Table 4, Table 5 and Table 6 in the following part.

To illustrate the advantages of the obtained designs for fitting PWO model, we compare the Ex-PSO-D designs for 4 components and 12 runs with the optimal PWO design in Peng, Mukerjee, and Lin (2019).

Furthermore, for the OofA experiment with 4 components, we generated optimal PWO designs with 7 to 23 runs using the Ex-PSO algorithm. Figure 4 shows the efficiencies of these designs. All the obtained designs with $n\ge 12$ reach at least $95\%$ efficiency of the full PWO design, though with less than one-fifth of the runs. Especially for the cases with $n=12$, the design attains the same efficiency as the full PWO design.

Furthermore, we compare four types of fractional PWO designs, which are COA, and the corresponding designs obtained by the threshold-accepting algorithm (Winker, Chen, and Lin 2020), the Federov’s exchange algorithm (which iteratively optimizes a delta function of the $x_i$ and x where $x_i$ is in the design and x is not, see reference to Section 3.3 in Fedorov 1972) and the Ex-PSO algorithm, denoted as $D_{coa}$, $D_{ta}$, $D_{ex}$ and $D_{ex-pso}$, respectively. $D_{ta}$ is the best result obtained over repeated runs of the threshold-accepting algorithm with up to $10000000$ iterations, $D_{ex}$ is generated by the optFederov function (implemented in the R library AlgDesign) with $nRepeats=5$, and $D_{ex-pso}$ is the best result obtained over five repeated runs of the Ex-PSO algorithm with $t_{ex}=20$ and $t_{pso}=100$. The optimal PWO design constructed by Peng, Mukerjee, and Lin (2019) which serves as a benchmark for evaluating fractional PWO designs is also listed here and denoted as $D_{peng}$. In addition, the new hybrid algorithm needs an exhaustive search over the design space during the single-point exchange procedure, and it can be computationally expensive if m is large. Hence, we only report designs with $n=m(m-1)/2+1,m(m-1),m!/r!(r=\lfloor m/2\rfloor )$ where $4\le m\le 7$. Nevertheless, given the tremendous growth in computational resources available, it is feasible to conduct the Ex-PSO algorithm for constructing designs with $m>7$.

Table 4, Table 5 and Table 6 exhibit the values of $det{(\tilde{M})}^{1/p}$, $tr({\tilde{M}}^{-1})$, $tr({\tilde{M}}^2)$, and D-, A- and $M.S.$-efficiency (in parentheses) for the corresponding designs. Note that the larger the value of $det{(\tilde{M})}^{1/p}$ is, the better, while smaller values of $tr({\tilde{M}}^{-1})$ and $tr({\tilde{M}}^2)$ are better. For any number of components, $D_{mp}$ is not unique and the corresponding $tr({\tilde{M}}^{-1})$ or $tr({\tilde{M}}^2)$ is by no means a fixed value. Hence, $D_{mp}$ is not listed in Table 5 and Table 6. From the tables, we can find that $D_{ex-pso}$ reaches a higher efficiency than the other types of designs under the PWO model in most cases. Further, we report the best PWO designs with $n=m(m-1)/2+1,m(m-1)$ and $4\le m\le 7$ under the D-, A- and $M.S.$-optimal criteria in the supplementary material.

We conclude this section with some numerical results on constructions of fractional PWO designs which have the same correlation structure as the full PWO design. Since these designs are 100% efficient under diverse design criteria including the D-, A-, $M.S.$-optimal criteria, we call them fully efficient PWO designs. Using the Ex-PSO algorithm, we find the following results.

For saving space, some selected fully efficient PWO designs with minimized runs for $m=4,5,6,7$ are exhibited in Table 7, other fully efficient PWO designs and the MATLAB codes for the Ex-PSO algorithm are available upon request.

In this article, we have presented a new hybrid algorithm for searching optimal PWO designs for OofA experiments. This algorithm starts with a set of randomly selected particle designs. It combines the single-point exchange algorithm for generating the local-best design related to each particle and the PSO algorithm for clustering the particle designs gradually around the global optimal design. We have addressed two points for this stochastic optimization technique: (i) the selection of initial particle designs; and (ii) the definition of the particle designs’ movement toward its personal local-best or global-best design. As the sample space is a large combinatorial set of $m!$ points, it’s satisfying to randomly select points for initial designs. And the movement is defined as randomly exchanging points to adjust the particle design toward the target design. Our generated design has an appealing efficiency in comparison with the full PWO design, though with only a small fraction of runs. And it achieves better efficiency compared with the existing designs. Also, we have reported optimal PWO designs and some selected 100% efficient PWO designs for application.

To illustrate the effectiveness of our model and designs, we provide a real-life OofA experiment. The data is from the four-drug data in Table 3 of Yang, Sun, and Xu (2021). The first-order PWO model fits using the data from the fully efficient PWO design with 12 runs (runs {1,2,6,7,10,11,16,18,20,21,22,24} of four-drug data) as well as the full design with 24 runs (i) with an adjusted R square of 0.93 (vs. 0.77), (ii) an RMSE of 2.19 (vs. 3.53) based on 5 (vs. 17) degrees of freedom. For further comparison, we consider variable selection to simplify the first-order PWO model. We start with a model with intercept only and perform stepwise regression with the AIC criterion. By using the data of fully efficient design, the stepwise regression leads to the following PWO model: $y=44.22-3.82z_{23}+6.02z_{34}+2.79z_{12}-2.96z_{14}+1.37z_{13}$, which has an adjusted R square of 0.95, an RMSE of 2.00 with 6 degrees of freedom. When using the data of full design, the stepwise regression leads to the following PWO model: $y=45.22-3.91z_{12}+3.86z_{23}-1.59z_{13}$, which has an adjusted R square of 0.79, an RMSE of 3.43 with 20 degrees of freedom. Thus, the first-order PWO model is sufficient for the four-drug OofA experiment. Meanwhile, the fully efficient design has an appealing efficiency, though with only a half run of the full design. However, in some applications, we are persuaded that the first-order PWO model cannot account for all systematic effects caused by component ordering, thus we shall consider the second-order PWO model to obtain a lower error. This approach can be referred to Mee (2020).

In addition, from the viewpoint of model robustness, it is interest to examine how our designs behave under high-order PWO models. From the perspective of model parsimony, we only consider the second-order PWO model with only a few interactions of two PWO factors sharing a common component being entertained, and recommend using the fully efficient PWO designs for analysis. We omit the details to save space but note that, even under such augmented models the fully efficient PWO designs tend to perform well, particularly under the D-, $M.S.$-criteria and often under the A-criterion. For example, the fully efficient design for $m=6,n=24$, has D-, A-, $M.S.$-efficiencies of at least 90%, under such augmented model. This allows further flexibility in model selection with a provision for penalty for the additional parameters that correspond to two-factor interactions. As before, one remains assured of relatively high design efficiency under the model that one arrives at.

As pointed out by reviewers, our work can be extended using alternation optimization techniques, such as the memetic optimization algorithm based on local searches (Abbasi-khazaei and Rezvani 2022). We conclude with the hope that the present endeavor will generate further interest in optimal designs for OofA experiments and related topics.