**PREPRINT**

國立臺灣大學 數學系 預印本 Department of Mathematics, National Taiwan University

### www.math.ntu.edu.tw/ ~ mathlib/preprint/2012- 15.pdf

## Particle Swarm Optimization Techniques for Finding Optimal Mixture Designs

### Weichung Wang, Ray-Bing Chen, Chien-Chih Huang, and Weng Kee Wong

### July 5, 2012

### Particle Swarm Optimization Techniques for Finding Optimal Mixture Designs

### Weichung Wang

^{∗}

### , Ray-Bing Chen

^{†}

### , Chien-Chih Huang

^{‡}

### , and Weng Kee Wong

^{§}

### July 5, 2012

Abstract

Particle Swarm Optimization (PSO) is a meta-heuristic algorithm that has been shown to be successful in finding the optimum solution or close to the optimum for a wide variety of real and complicated op- timization problems in engineering and computer science. This paper adapts PSO methodology by first solving an optimization problem on the hypercube and then projecting the solution onto the q-simplex op- timization space to find different optimal designs for mixture models commonly used in agronomy, food science and pharmaceutical sci- ence. We show that PSO is also flexible in that it can be modified straightforwardly to find optimal designs for log contrast models and constrained mixture models. We conclude with a list of advantages of this simple and novel way for finding optimal mixture designs over current methods.

Keywords. A-, D- and L-optimal design; Efficiency; Equivalence theorem; Projection; Regression model.

### 1 Introduction

Mixture experiments are widely used in food processing, chemical, manu- facturing, agricultural, cosmetics and pharmaceutical industries. This list

∗Department of Mathematics, National Taiwan University, Taipei 106, Taiwan

†Department of Statistics, National Cheng-Kung University, Tainan 70101, Taiwan

‡Department of Mathematics, National Taiwan University, Taipei 106, Taiwan

§Department of Biostatistics, University of California, Los Angeles, CA 90095, USA

1

is growing as interest in mixture models expands and optimal designs are becoming more available thorough computer codes, software packages and interactive online websites. [10] provides an excellent introduction and broad coverage in mixture experiments. A comprehensive review of mixture models and their optimal designs can be found in [5] and a more recent review of research work in mixture experiments over the last 50 years is the work of Piepel in Chapter 12 of an edited volume by [23].

The aim of this paper is to introduce a modified version of a popular optimization technique already widely used in engineering and computer sci- ence research for finding optimal designs for mixture models. The technique has been around for more than ten years, but interestingly hasn’t had much impact in statistical methodology to date. Our experience reinforces the widely held findings that PSO is a very simple and powerful optimization tool. PSO requires no assumption of the objective function and has few and easy to work with parameters. PSO is intriguing because despite its lack of a firm theoretical basis, its repeated successes and increasing widespread use in various disciplines have even resulted in new journals that simply track its applications in various fields.

Section 2 provides statistical setup and briefly reviews general theory for optimal experimental designs before discussing mixture models and their op- timal designs. In Section 3, we propose a meta-heuristic method for finding a variety of optimal designs for different types of mixture experiments. The stochastic search is based on a modified version of the particle swarm opti- mization first proposed by [22]. In Section 4 we demonstrate this procedure is efficient for finding various optimal designs for different types of mixture models, including new models for which optimal designs are not known ana- lytically. We conclude in Section 5 with a summary of the advantages of the proposed PSO method over current methods for generating optimal designs.

### 2 Background

Our interest is in the general linear model given by

y = β^{0}f (x) + . (1)

Here y is the response variable, β is a t-dimensional column vector of unknown
coefficients, f (x) is a given vector of regression functions defined on a user-
defined compact multi-dimensional design space. The error is random noise
with zero mean and constant variance and we assume all errors are normally
and independently distributed. An approximate design ξ is defined by its
design points (x_{i}’s) and the proportions (p_{i}’s) of observations to be taken

at these points. Once the sample size n is fixed, either by cost or time
considerations, and an optimality criterion is given, the optimization problem
is to determine the number (k) of points needed and the values of x_{i}, p_{i}, i =
1, . . . , k subject to p_{1}+ . . . + p_{k} = 1. The implemented design takes roughly
np_{i} observations at x_{i}, i = 1 . . . , k subject to np_{1}+ . . . + np_{k} = n.

Following convention, the worth of a design is measured by its Fisher
information matrix which is obtained by taking the negative of the expec-
tation of the second derivative of the logarithm of the likelihood function
with respect to β. Given the above statistical model, this matrix is M(ξ) =
E_{ξ}(f (x)f (x)^{0}), which is inversely proportional to the variance-covariance ma-
trix of the estimated parameters β. Following convention, the design criterion
is formulated as a convex function of the information matrix. For example,
D-optimality for estimating model parameters seeks to minimize the gen-
eralized variance using the convex functional Φ(ξ) = −ln |M (ξ)|. Another
popular and useful criterion is L-optimality defined by Φ(ξ) = tr L M (ξ)^{−1}
and L is a user-selected matrix; if the goal is to minimize the average of the
variances of the estimated parameters, we set L = I, the identity matrix
whereupon L-optimality reduces to A-optimality. Alternatively, if we like
to estimate some average of the response over a user-selected region R, one
chooses L =R

Rf (x)f^{0}(x)µ(dx) and µ is a a selected measure over for R. This
corresponds to estimating the response surface over R with weights specified
by the measure µ with more important parts of R receiving a higher weight
assignment from µ. If there is equal interest over the region R, one chooses
µ to be the uniform measure on R. We obtain I-optimality when µ is the
uniform measure and R is the same as the design space.

When the design criterion is a convex function of the information matrix,
as it is for all the above criteria, [27], [25, 26] showed an effective way to check
whether a design is optimal for any model. In particular if Model (1) holds
and each of the following inequalities holds for all x in the design space, then
(a) a design ξ is D-optimal if and only if f (x)^{0}M^{−1}(ξ)f (x) ≤ t; (b) a design
ξ is L-optimal if and only if f (x)^{0}M^{−1}(ξ) L M^{−1}(ξ)f (x) ≤ trL M^{−1}(ξ).

The above results are frequently referred to as equivalence theorems or more informally as checking conditions and are derived from considerations of the Frechet derivatives of the convex functionals. When the regression model has one or two independent variables, the equivalence theorem can be easily applied to check the optimality of any design graphically. For instance, to check for D-optimality, we plot the function on the left hand side of the inequality (a) over the design space and determine whether the inequality is satisfied throughout the design space. If it is, the design ξ is D-optimal;

otherwise it is not. In the latter case, one can also obtain a lower bound on the efficiency of the design by examining the plot (without knowing the

3

optimum). Details are in standard design monograph such as [37].

### 2.1 Mixture Models

Clearly optimal designs are dependent on model assumptions and for non- linear models they depend on the nominal values of the model parameters as well, see for example, [8], where they considered using a D-optimal design to estimate parameters in an inverse quadratic with 4 parameters in agronomy.

In mixture problems, it appears much of the recent design work for mixture models focuses on constructing designs robust to model mis-specification.

For example, [18] is one of the recent papers to study robustness properties of optimal designs for mixture models when there is uncertainty in the model assumptions. Specifically, they investigated how to find A-optimal designs for mixture experiments that are robust to the linear and quadratic models proposed by [38]. In addition, [16] and, [17] found D- and A-optimal designs for linear log contrast and quadratic log contrast models for experiments with mixtures, respectively. Some new criteria were proposed in [29] and [36]. The former advocated a trace criterion to estimate the best propor- tions for the ingredients or components and the latter explored a minimax criterion to estimate the response surface in a mixture experiment, including using a deficiency criterion to measure the goodness of a mixture experiment.

In both papers, the model was a quadratic polynomial in several factors over the simplex region.

We assume our mixture experiments have q factors x1, x2, . . . , xq defined
on the regular q simplex S^{q−1} = {x^{0} = (x_{1}, x_{2}, . . . , x_{q}) ∈ [0, 1]^{q} : Pq

i=1x_{i} =
1}. Some common mixture models used in practice are Scheff´e’s polynomials
of degree n. If denotes random error, the simplest is an additive polynomial
mixture model when n = 1 and f (x)^{0} = (x_{1}, x_{2}, . . . , x_{q}) given by

y = β^{0}x =

q

X

i=1

β_{i}x_{i}+ , (2)

When f (x)^{0} = (x_{1}, x_{2}, . . . , x_{q}, x_{1}x_{2}, x_{1}x_{3}, . . . , x_{q−1}x_{q}) and n = 2 the second
degree Scheff´e’s polynomial mixture model is

y =

q

X

i=1

βixi+ X

1≤i<j≤q

βijxixj+ . (3)

This is an example of a Scheff´e quadratic canonical polynomial models widely used in blending experiments in engineering, agriculture, biological and the medical sciences. In the notation of model (1), t = q for model (2) and

t = q + q(q + 1)/2 for model (3). More generally, the Scheff´e polynomial of order n for a q-component mixture model is

y =

q

X

i=1

θ_{i}x_{i}+ X

1≤i<j≤q

φ_{ij}x_{i}x_{j} + · · · + X

1≤i1<···<in≤q

φ_{i}_{1}···inx_{i}_{1}· · · x_{i}_{n}+ . (4)

[2, 3] proposed a flexible modeling techniques for studying mixture experi- ments with additive effects using Becker polynomials. They are also useful for modeling the mean response when it also depends linearly on the to- tal amount A used in the experiment and all components in the regression function are homogenous of degree 1. Becker’s models include

y =

q

X

i=1

θ_{i}x_{i}+X

i<j

φ_{ij}min(x_{i}, x_{j}) + · · · + φ_{1,2,···,q}min(x_{1}, · · · , x_{q}) + , (5)

y =

q

X

i=1

θ_{i}x_{i}+ X

1≤i<j≤q

φ_{ij}x_{i}x_{j}

x_{i}+ x_{j} + · · · + φ_{1,2,···,q}x_{1}x_{2}· · · , x_{q}

(x_{1}+ · · · + x_{q})^{q−1} + , (6)

y =

q

X

i=1

θ_{i}x_{i}+ X

1≤i<j≤q

φ_{ij}(x_{i}x_{j})^{1/2}+ · · · + φ_{1,2,···,q}(x_{1}x_{2}· · · x_{q})^{1/q} + . (7)
In metallurgy, when there are q = 2 ingredients in a mixture experiment,
[20] found a certain type of polynomials is useful for modeling the response
and he called them Kasatkin’s polynomials. Such a polynomial of n^{th} order
has the form:

y = θ_{1}x_{1}+ θ_{2}x_{2}+

n−2

X

i=0

φ_{i}x_{1}x_{2}(x_{1}− x_{2})^{i}+ . (8)

Among these models, Scheff´e’s models are the most widely used. Fur- ther details on rationale and specific applications of Becker’s and Kasatkin’s models can be found in [10], [20], [39], etc. Interestingly, these papers allude to D-optimal designs for Kasatkin’s polynomial models but we were unable to find them in English publications. We used our proposed method and generate D-optimal designs for some models in Table 4.

### 2.2 Optimal Designs

We first briefly review known optimal designs for mixture models from the literature. Analytical descriptions of D-optimal designs for different orders of Scheff´e polynomial models were reported in [38], [24], [40], [14], [15], [42]

5

and, [28]. Formulae for A- and integrated or I-optimal designs are available for a much smaller class of models.

[38] found the theoretical A- and D-optimal designs for the first order
linear models with q factors over S^{q−1}. Both A- and D-optimal designs
coincide and are equally supported on the q vertices of the simplex given
by (1, 0, . . . , 0), . . . , (0, . . . , 0, 1). The A-optimal design for the quadratic
mixture model with q ≥ 4 was found by [42] where they showed that the
A-optimal design is the weighted {q, 2} simplex-centroid design. It has a
combined weight of (4q − 3)^{1/2}/(q(4q − 3)^{1/2}+ 2q(q − 1)) equally distributed
among support points of the form (1, 0, . . . , 0), . . . , (0, . . . , 0, 1), and a com-
bined weight of 4r_{1}/(4q − 3)^{1/2} equally distributed among points of the form
(1/2, 1/2, 0, . . . , 0), . . . , (0, . . . , 0, 1/2, 1/2). When q = 3, they numerically
identified the A-optimal as the weighted {q, 3} simplex-centroid design with
(r_{1}, r_{2}, r_{3}) = (0.1418, 0.1873, 0.0128), where r_{1}, r_{2} are as before and r_{3} is now
the weight at each of the support point of the form (1/3, 1/3, 1/3, 0 . . . , 0), . . . ,
(0, . . . , 0, 1/3, 1/3, 1/3). Table 1 shows these optimal designs and the {q, 2}

simplex-centroid designs which are A- or D-optimal for the quadratic models.

We next consider two third-degree polynomial models for mixture studies.

The first one is the cubic model without 3-way effect polynomial mixture model, i.e.

E(y) =

q

X

i=1

β_{i}x_{i}+ X

1≤i<j≤q

β_{ij}x_{i}x_{j} + X

1≤i<j≤q+1

γ_{ij}x_{i}x_{j}(x_{i}− x_{j}). (9)

The D-optimal design was found by [30] to be equally supported at the
the following design points: C_{1}^{q+1} points given by x_{i} = 1, x_{j} = 0, i 6= j, i =
1, . . . , q, 2C_{2}^{q+1} points given by x_{i} = 1 − x_{j} = ^{1}_{2}(1 −^{√}^{1}

5), i 6= j, i, j = 1, . . . , q,
and x_{k} = 0, k 6= i, j. The second third-degree polynomial model of interest is
the full cubic model with all 2-factor and 3-factor interactions:

E(y) =

q

X

i=1

β_{i}x_{i}+ X

1≤i<j≤q

β_{ij}x_{i}x_{j}+ X

1≤i<j≤q

γ_{ij}x_{i}x_{j}(x_{i}−x_{j})+ X

1≤i<j<k≤q+1

β_{ijk}x_{i}x_{j}x_{k}.
(10)
The D-optimal design was also found by [31] and is equally supported at
the following design points: C_{1}^{q+1} points given by x_{i} = 1, x_{j} = 0, i 6= j, i =
1, . . . , q, 2C_{2}^{q+1} points given by x_{i} = 1 − x_{j} = ^{1}_{2}(1 −^{√}^{1}

5), i 6= j, i, j = 1, . . . , q,
x_{k} = 0, k 6= i, j, and C_{3}^{q+1} points given by x_{i} = x_{j} = x_{k} = 1/3, x_{l} = 0,
l 6= i, j, k; i, j, k = 1, 2, . . . , q + 1.

An analytical description of the optimal design is clearly desirable but there are also obvious limitations of such a theoretical approach. This is because (i) formulae exist only for relatively simple problems and even when

they exist, they are complicated and are obtained only after long and tedious mathematical calculation that at times require exploiting the mathematical properties of the specific regression functions, (ii) frequently artificial as- sumptions are imposed to arrive at a closed form solution and (iii) the math- ematical derivation typically works only for one specific model and quickly breaks down for a sub-model or a slight generalization of the model. An example of the latter situation is finding D-optimal designs for the above two cubic models with one or two missing interaction terms.

For these reasons, search algorithms are more practical approaches to find optimal designs than one that relies on mathematics alone. Statistical software packages like SAS, JMP and Design-Expert have some routines for finding selected types of optimal designs for some mixture models, however, the scope still seems limiting. Our experience is that most statistical software packages do not have the flexibility of generating an optimal design for a user- specified model directly, numerically or otherwise.

The next section describes a PSO-based algorithm with great potential for finding many types of optimal designs quickly for a wide variety of mixture models, including optimal designs for sub-models in (10) for which analytical results remain elusive.

### 3 Particle Swarm Optimization with Projec- tion Capabilities

Given the mixture model and the optimality criterion, PSO begins its search for the optimal mixture design with a population of randomly generated candidate designs or particles that cover the design space. Each design se- quentially adapts its movement toward where it believes is the optimum and does so with an velocity that depends on its current location and locations that other particles believe is the optimum. Usually these movements are governed by two key equations with a few parameters that includes random vectors that are responsible for the stochastic movement of the designs.

For our problem, we first use a projection function to optimize over the hypercube instead of the simplex. This simplifies the optimization problem and speeds up the computation process. To fix ideas, suppose the given mixture model has q factors and we wish to find a k-point optimal design.

Let m = k × (q + 1) and let Ξ = [0, 1]^{m} denote the m-dimensional hypercube.

Define the m × 1 vector, ˜ξ = (x^{0}_{1}, . . . , x^{0}_{k}, p^{0})^{0} ∈ Ξ, where xi is a q × 1 vector
in [0, 1]^{q}, i = 1, . . . , k, and p ∈ [0, 1]^{k}. Let Ξ^{∗} = Ξ \ { ˜ξ = (x^{0}_{1}, . . . , x^{0}_{k}, p^{0})^{0} ∈
Ξ | 1^{0}_{k}· p = 0 or 1^{0}_{q}· x_{i} = 0 for some i}. To transform ˜ξ into a proper design

7

ξ, we define the projection function P : Ξ^{∗} −→ (S^{q−1})^{k}× S^{k−1} by

P ( ˜ξ) = ( x^{0}_{1}

(1^{0}_{q} x1), . . . , x^{0}_{k}

(1^{0}_{q} xk), p^{0}

(1^{0}_{k} p))^{0}. (11)
The projection function P is invariant in the sense that P ◦ P ( ˜ξ) = P ( ˜ξ) and
the design ξ has support on x^{0}_{i}/(1^{0}_{q} xi), i = 1, . . . , k and the components in
p^{0}/(1^{0}_{k} p) are the corresponding weights. The notation ξ = P (˜ξ) signifies
that the design ξ is transformed from ˜ξ via the projection P .

Our particle swarm optimization is based on the projection function P
in Eq. (11) as follows. We first initialize a random population of n candi-
dates of k-point designs from Ξ^{∗}. We define two notion at each stage of the
iteration: let (i) ˜ξ_{i}^{pbest} denote the personal best position for the i^{th} particle,
i.e. the design ˜ξ_{i}^{pbest} provides the optimal value for the criterion among all
the positions that the i^{th} particle has ever visited, and (ii) let ˜ξ^{gbest} denote
the global best position, i.e. ξ^{gbest} provides the optimal value for the crite-
rion among all the positions that all of the particles have ever visited. The
strategy for the i^{th} particle, ˜ξ_{i} at the t^{th} iteration is as follows:

• Generate a new velocity v_{i}^{t} to reach to next position given by v_{i}^{t} =
wtv^{t−1}_{i} +c11( ˜ξi

gbest

− ˜ξ_{i}^{t−1})+c22( ˜ξ_{i}^{pbest}− ˜ξ_{i}^{t−1}), where v^{t−1}_{i} is the velocity
generated at the t − 1 iteration, w_{t} is the inertia weight, c_{1}, c_{2} are two
pre-specified positive constants, and _{1}, _{2} are m × 1 uniform random
vectors.

• The next location for the i^{th} particle is

ξ˜^{t}_{i} = ˜ξ_{i}^{t−1}+ χv^{t}_{i}, (12)
where χ is a pre-specified positive constant. If ˜ξ_{i}^{t}is not in Ξ^{∗}, we project
ξ˜_{i}^{t} to a location closest to the boundary of Ξ^{∗}.

• Obtain the current design ξ^{t} by projecting ˜ξ_{i}^{t} to the simplex using i.e.

ξ_{i}^{t}= P ( ˜ξ_{i}^{t}).

• Evaluate the objective function at the new design ξ_{i}^{t}.

• Update the current best for each particle ξ_{i}^{pbest} and ˜ξ_{i}^{pbest}.

• Update the inertia weight wt+1 = g(wt), where g is a user-selected monotonic decreasing function.

After updating all particles, ˜ξ_{i}^{t}, we identify ˜ξ^{gbest} and ξ^{gbest}, and repeat the
procedure. The procedure terminates and reports ξ^{gbest} as our “best” design
after a pre-specified maximal number of iterations is reached or when the cri-
terion value does not change much according to some user-specified tolerance
level.

In the above PSO steps, the constant c_{1} is the cognitive learning factor
and the constant c_{2} is the social learning factor. Following [21], the conven-
tion is to set c_{1} = c_{2} = 2. The constant χ is the constriction factor and
usually is set to be 1. The function g is a pre-selected non-increasing func-
tion to control the time-varying inertia weight w_{t}. In our work, we set the
function g to be the linear decreasing function such that w_{t} varies from 0.9
to 0.4.

The key advantage of this modified PSO is that it operates on the simple
hypercube first before it projects any non-feasible point into (S^{q−1})^{k}× S^{k−1}.
This simplifies and makes the computation more efficient. If we had directly
implemented PSO to search for the optimal mixture design and worked with
the constraints in the simplex, our experience is that some of the sequentially
generated particles “flew” outside the simplex and the subsequent work re-
quired to ignore them or bring them back to the simplex can complicate
and prolong the search for the best design considerably. This proposed PSO
based on the projection is called ProjPSO.

### 4 Numerical Results

We now demonstrate our ProjPSO algorithm is an effective way to generate different types of optimal designs for various mean functions in a mixture experiment. For comparison sake, the examples we work with include results already reported in the literature and also new results where analytical for- mulae for the optimal designs are not available. In the latter case, we used an equivalence theorem to confirm its optimality and if the PSO-generated de- sign is not optimal, we assessed its proximity to the optimal (without knowing the optimum) via an efficiency lower bound obtained from the equivalence theorem. Details in [37].

### 4.1 ProjPSO codes

We implemented the ProjPSO algorithm written in MATLAB codes and they are available from the authors upon request. Our PC has a 2.67GHz Intel(R) Core(TM) i7 CPU. We always start with a modest size of the particles and a modest number of the iterations and increase them when the dimensionality

9

of the experimental region gets larger or the model has more parameters.

We follow convention and use default values for the other parameters in
PSO algorithm, for example, setting c_{1} and c_{2} both equal to 2 in ProjPSO.

In almost all our examples, PSO was able to generate designs which were optimal or very close the theoretical optimal designs after a few minutes of CPU time. For example, to find the D-optimal design for the full cubic model with 3 factors, we implemented the codes using 1024 particles and 200 iterations and it took ProjPSO around 2.5 minutes to produce the analytical D-optimal design.

For each mixture problem, the main generic part of the code remains intact and only the information matrix and the criterion need to be changed in the code. For all our examples, we first searched among all designs with design points equal to the number of parameters in the model using the ProjPSO algorithm. Our guiding principle was larger number of particles or larger number of iterations for more complex models. This is because the time required to generate the optimal design is usually fast and the differences in additional computational time required when the number of particles or iterations are increased is usually not long. For instance, in the examples below, the number of particles we chose to generate the optimal designs for the linear Scheff´e polynomial models were 64, 128 and 256 for q = 3, 4 and 5 and the corresponding number of iterations used were 200, 400 and 800.

We applied ProjPSO to find the optimal designs when there are q factors
and q = 3, 4 and 5. The generated A- and D-optimal designs for Scheff´e’s
linear mixture models are all numerically the same as the theoretical A-
and D-optimal designs reported in Table 1. Figure 1 is an illustrative PSO-
generated plot that shows the particles representing the support points of
different designs at various stages in the search when q = 3. The sub-figures
display the support points of the 64 particles at the 1^{st}, 5^{th}, 10^{th}, 20^{th}, 30^{th}
and the 40^{th} iterations and show the very rapid convergence of the ProjPSO
procedure to the optimum. In this case, ProjPSO only takes less than 3
seconds for 40 iterations. When we applied the ProjPSO algorithm to find
D-optimal designs for the Scheff´e’s quadratic mixture model (3), we also
obtained the {q, 2} simplex-centroid design in Table 1 and shown by Kiefer
to be D-optimal [24].

As specific examples, we applied the ProjPSO algorithm to search for these A-optimal designs using 1024 particles 600 and 400 numbers of iter- ations for q = 3 and 4 respectively. The PSO-generated designs for these two cases were found to be identical to the theoretical optimal designs for all practical purposes.

We also applied ProjPSO to models with more complex structures. Con- sider the full cubic model with q = 3 and the cubic model does not contain

Figure 1: The movement of particles in the PSO search for the D-optimal
design for the linear mixture model with q = 3 factors. Each subfigure
displays the PSO-generated mixture design at a particular iteration. At the
40^{th} iteration, ProjPSO seems to have converged in 4 seconds of CPU time
to the D-optimal design equally supported at the vertices.

11

Figure 2: The plot of the directional derivative of the PSO-generated design confirms that ξP SO is I-optimal design for cubic Scheff´e model.

3-way effects. The corresponding best designs generated by ProjPSO are all the same these reported in [30, 31]. Here the number of iterations in ProjPSO is set to be 400, but the number of particles are 128 and 512 for the cubic model without 3-way effect and full cubic model respectively. In addition to A- and D-optimal criteria, we also study I-optimal design for the Scheff´e’s quadratic and cubic mixture models. To generate the numerical I-optimal designs, we set 1024 particles and iterate ProjPSO 400 times. The results are shown in Table 2. The equivalence plot of the I-optimal design for the cubic Scheff´e model is shown in Figure 2.

We also used ProjPSO to determine optimal designs for several submodels (or incomplete models) obtained by deleting a few interaction terms from the full cubic model. As far as we know, theoretical optimal designs for these submodels or incomplete (IC) models remain unknown. We used 1024

Figure 3: The plot of the directional derivative of the PSO-generated design
confirms that ξ_{P SO} is D-optimal design for Katkasin polynomial model with
order 5.

13

Figure 4: The plot of the directional derivative of the PSO-generated design
confirms that ξ_{P SO} is D-optimal design for IC Model 1.

particles and iterated 400 times. To ensure the generated designs are D- optimal for the submodels, we used an equivalence theorem for each submodel and checked for its optimality. Figure 3 is an example of graphical version of equivalence theorem for the submodel IC Model 1 in Table 3 and since the 3-dimensional plot is bounded above by 0, optimality of the reported corresponding design in Table 3 is confirmed.

In addition to the Scheff´e type models, we also apply ProjPSO to generate D-optimal designs for the other mixture models to demonstrate the flexibility of our ProjPSO. Here three Becker’s models, (5) to (7), with q = 3 and Kasatkin’s polynomials, (8), with 3, 4, and 5 orders are considered. In our ProjPSO code, we only need to change the regressor set-up according to the target model. The numerical best designs for all 6 models are satisfied the corresponding equivalence theorem. The equivalence plot for the D-optimal design of Kasatkin’s polynomial with order 5 is shown in Figure 3.

Table1:DifferenttypesofPSO-generateddesignsmatchedtheoreticalresultsfordifferentScheff´epolynomialmodels. modeltypecriteriontypeofsupportpointsweightscomment 1Linear(q=3)D(1,0,0)1/3×3Scheff´e(1958) 2Linear(q=4)D(1,0,0,0)1/4×4Scheff´e(1958) 3Linear(q=5)D(1,0,0,0,0)1/5×5Scheff´e(1958) 4Linear(q=3)A(1,0,0)1/3×3Scheff´e(1958) 5Linear(q=4)A(1,0,0,0)1/4×4Scheff´e(1958) 6Linear(q=5)A(1,0,0,0,0)1/5×5Scheff´e(1958) 7Quadratic(q=3)D(1,0,0)1/6×3Kiefer(1961) (1/2,1/2,0)1/6×3 8Quadratic(q=4)D(1,0,0,0)1/10×4Kiefer(1961) (1/2,1/2,0,0)1/10×6 9Quadratic(q=3)A(1,0,0)0.142×3YuandGuan(1993) (1/2,1/2,0)0.187×3 (1/3,1/3,1/3)0.0128×1 10Quadratic(q=4)A(1,0,0,0)0.094×4YuandGuan(1993) (1/2,1/2,0)0.104×6 15

Table2:PSO-generatedDandI-optimaldesignsmatchedtheoreticaldesignsforvariousScheff´epolynomialmodels (cont’d). modeltypecriteriontypeofsupportpointsweightscomment 11CubicwithoutD(1,0,0)1/9×3Mikaeili(1989) 3-wayeffects(0.2764,0.7236,0)1/9×3 (q=3)(0.7236,0.2764,0)1/9×3 12FullCubic(q=3)D(1,0,0)1/10×3Mikaeili(1993) (0.2764,0.7236,0)1/10×3 (0.7236,0.2764,0)1/10×3 (1/3,1/3,1/3)1/10×1 13Quadratic(q=3)I(1,0,0)0.1002×3LiuandNeudecker(1995) (1/2,1/2,0)0.2016×3 (1/3,1/3,1/3)0.0949×1 14Cubic(q=3)I(1,0,0)0.0925×3LiuandNeudecker(1995) (1/2,1/2,0)0.1483×3 (1/3,1/3,1/3)0.2776×1

Table3:PSO-generatedD-optimaldesignsfornew3-factorsubmodelsfromScheff´epolynomialmodels. modeltypesupportpointsweightscomment 15

P iβixi+β13x1x3+β23x2x3(1,0,0),(0,1,0)0.0833×2new +P i<jγijxixj(xi−xj)+β123x1x2x3(0,0,1)0.1111×1 (0.2764,0.7236,0),(0.2764,0,0.7236)0.1111×4 (0.2113,0.7887,0),(0.7887,0.2113,0)0.0833×2 (0.3333,0.3333,0.3333)0.1111×1 16

P iβixi+β23x2x3(1,0,0)0.0938×3new +P i<jγijxixj(xi−xj)+β123x1x2x3(0,0.2764,0.7236),(0,0.7236,0.2764)0.1250×2 (0.7887,0.2113,0),(0.7887,0,0.2113)0.0937×2 (0.2113,0.7887,0),(0.2113,0,0.7887)0.0937×2 (0.3333,0.3333,0.3333)0.1250×1 17 P iβixi+

P i<jβijxixj(1,0,0)0.1250×3new +γ13x1x3(x1−x3)+γ23x2x3(x2−x3)(0,0.2764,0.7236),(0.2764,0,0.7236)0.1250×2 (0.7236,0,0.2764),(0,0.7236,0.2764)0.1250×2 (0.5,0.5,0))0.1250×1 18

P iβixi+

P i<jβijxixj(1,0,0)0.1111×3new +γ13x1x3(x1−x3)+γ23x2x3(x2−x3)+β123x1x2x3(0,0.2764,0.7236),(0.2764,0,0.7236)0.1111×2 (0.7236,0,0.2764),(0,0.7236,0.2764)0.1111×2 (0.5,0.5,0))0.1111×1 (0.3333,0.3333,0.3333)0.1111×1

17

Table4:PSO-generatedD-optimaldesignsforBeckerandKasatkin’smixturemodelswith3factors. modeltypesupportpointsweightscomment 19BeckerModel1(1,0,0)0.1429×3Becker(1968,1978) P iβixi+P i<jβij(xixj)(1/2) (0.5,0.5,0)0.1429×3 +β123(x1x2x3)(1/3) (0.3333,0.3333,0.3333)0.1429×1 20BeckerModel2(1,0,0)0.1429×3Becker(1968,1978) P iβixi+β12x1x2/(x1+x2)+β13x1x3/(x1+x3)(0.5,0.5,0)0.1429×3 +β23x2x3/(x2+x3)+β123x1x2x3(0.3333,0.3333,0.3333)0.1429×1 21BeckerModel3(1,0,0)0.1429×3Becker(1968,1978) P iβixi+β12min{x1,x2}+β13min{x1,x3}(0.5,0.5,0)0.1429×3 +β23min{x2,x3}+β123min{x1,x2,x3}(0.3333,0.3333,0.3333)0.1429×1 22Kasatkin3rdOrderModel(1,0,0),(0,1,0)1/4×2Kasatkin(1974) P2 i=1θixi+P1 i=0φix1x2(x1−x2)i (0.2764,0.7236,0),(0.7236,0.2764,0)1/4×2 23Kasatkin4thOrderModel(1,0,0),(0,1,0)1/5×2Kasatkin(1974) P2 i=1θixi+P2 i=0φix1x2(x1−x2)i (0.1727,0.8273,0),(0.8273,0.1727,0)1/5×2 (0.5,0.5,0)1/5 24Kasatkin5thOrderModel(1,0,0),(0,1,0)1/6×2Kasatkin(1974) P2 i=1θixi+P3 i=0φix1x2(x1−x2)i (0.1175,0.8825,0),(0.8825,0.1175,0)1/6×2 (0.3574,0.6426,0),(0.6426,0.3574,0)1/6×2

### 4.2 Constrained Mixture Experiments

4.2.1 The linear log contrast models

This subsection shows ProjPSO is flexible and can be directly modified to find optimal designs for related mixture models. Consider the linear log contrast model proposed by [1] for studying a certain type of mixture experiments.

Such models continue to be of interest and were recently studied by [16] and, [17] where they found exact D- and A-optimal designs for linear log contrast and quadratic log contrast models. [4] found the continuous D-optimal design for the log contrast model,

E(y) = β_{0}+

q−1

X

i=1

β_{i}log(x_{i}/x_{q}).

To ensure a D-optimal design exists, additional constraints on all the factors
are required. One common way to do this is to select a constant δ ∈ (0, 1)
with the conditions δ ≤ x_{i}/x_{j} ≤ 1/δ, for all 1 ≤ i, j ≤ q as added constraints
on the design region S^{q−1}.

As an illustration, consider the log contrast model with q = 3. [4] showed that for a given δ, the D-optimal design has 3 points and is supported equally at (1/(1 + 2δ), δ/(1 + 2δ), δ/(1 + 2δ)), (δ/(1 + 2δ), 1/(1 + 2δ), δ/(1 + 2δ)), (δ/(1 + 2δ), δ/(1 + 2δ), 1/(1 + 2δ)), or (1/(2 + δ), 1/(2 + δ), δ/(2 + δ)), (δ/(2 + δ), 1/(2 + δ), 1/(2 + δ)), (1/(2 + δ), δ/(2 + δ), 1/(2 + δ)).

To find the optimal design using ProjPSO, we redefined the regressors as
log(x_{i}/x_{q}) and also amended the projection operator in ProjPSO so that
it projects into the right space that includes the additional constraints,
δ ≤ x_{i}/x_{j} ≤ 1/δ for all i. We used ProjPSO to find the D-optimal de-
signs when δ = 0.145 and 0.2. Using a flock size of 1024 and 100 number of
iterations, ProjPSO took approximately 11 seconds of CPU time to generate
the D-optimal designs below, which also agree with the result in [4]. For each
of these two δ’s, there are two optimal designs equally supported at 3 points.

For δ = 0.145, one set of support points is x_{1}^{0} = (0.1124, 0.1124, 0.7752),
x_{2}^{0} = (0.7752, 0.1124, 0.1124) and x_{3}^{0} = (0.1124, 0.1124, 0.7752) and the
other set is x_{1}^{0} = (0.4662, 0.4662, 0.0676), x_{2}^{0} = (0.0676, 0.4662, 0.4662)
and x_{3}^{0} = (0.4662, 0.0676, 0.4662). For δ = 0.2, one set of support points
is x_{1}^{0} = (0.7143, 0.1429, 0.1429), x_{2}^{0} = (0.1429, 0.7143, 0.1429) and x_{3}^{0} =
(0.1429, 0.1429, 0.7143) and the other set is x_{1}^{0} = (0.0909, 0.4545, 0.4545),
x_{2}^{0} = (0.4545, 0.0909, 0.4545) and x_{3}^{0} = (0.4545, 0.4545, 0.0909).

It is interesting to note here that this is one situation where we did not
obtain good results using the default values c_{1} = c_{2} = 2 in the PSO al-
gorithm. This may be due to the smaller design space resulting from the

19

several constraints. Our experience suggests that setting c_{1} = c_{2} = 0.5 seems
to work well for log contrast models.

4.2.2 Another constrained mixture problem

Our last example concerns mixture experiments with constraints on the com-
ponents imposed physically by the problem. Because of practical or cost
considerations, upper or lower bound constraints are imposed on some of
the x_{i}’s with user-specified constants L_{i} and U_{i}, such that L_{i} ≤ x_{i} ≤ U_{i},
i = 1, 2, ..., q. A commonly cited example is in the making of a 3-component
fruit juice comprising watermelon, pineapple and orange. The limits on the
first component are .4 and .8, the limits on the second component are .1 and
.5 and the limits on the third component are .1 and .3. The design question
is how to blend these components so that the drink tastes best subject to
some cost constraints. In this case, the mixture has a larger requirement
on the watermelon component because it is the cheapest among the three.

Examples where mixture experiments have constraints on the components abound in pharmaceutical problems as well. For instance in tablet formu- lations, typically a D-optimal design is sought in the constrained mixture design with limits imposed on the various ingredients, see for example [13], [19], [34] and [35].

ProjPSO algorithm can also be modified directly to generate optimal
designs for mixture experiments with constraints on its components. We
discuss one such application due to space consideration. [7] employed a
mixture design to study amphiphilic cyclodextrin nano particles. There were
three variables z_{1}, z_{2} and z_{3} and their constraints were 0.4 ≤ z_{1} ≤ 0.7;

0 ≤ z_{2} ≤ 0.6 and 0 ≤ z_{3} ≤ 0.6. Working with pseudo-variables, [11] searched
for an optimal design inside the pseudo-variable simplex region defined by x_{1},
x_{2} and x_{3} with the following constraints: 0 ≤ x_{1} ≤ 0.5; 0 ≤ x_{2}, x_{3} ≤ 1 and
P3

i=1x_{i} = 1. Because this region for the pseudo-variables is not the regular
simplex, we modified our ProjPSO algorithm to ensure that all design points
are projected correctly into the proper region.

Choisnard et al. (2005) used the full quadratic model in their work with- out explanation but because this model contains the intercept, along with the constraint that the components sum to unity, the information matrix of any design for this model is singular. Accordingly we worked with an illustrative cubic mixture model without the 3-way effect, i.e.

E(y) =

q

X

i=1

β_{i}x_{i}+ X

1≤i<j≤q

β_{ij}x_{i}x_{j}+ X

1≤i<j≤q+1

γ_{ij}x_{i}x_{j}(x_{i}− x_{j}). (13)

Figure 5: The plot of the directional derivative of the PSO-generated de-
sign confirms that ξ_{P SO} is D-optimal design on the restricted design space
P3

i=1x_{i} = 1, 0 ≤ x_{1} ≤ 0.5 and 0 ≤ x_{2}, x_{3} ≤ 1.

There are two constraints for the variables, x_{i}’s. One is P3

i=1x_{i} = 1 and
another is x1 ≤ 0.5. To find the optimal design, we modified ProjPSO by
including the additional constraint and found a 9-point optimal design using
1024 particles and 1000 iterations. This design ξ_{pso,D}^{3} is equally supported at
the 9 support points xi0 = (x1, x2, x3) shown as columns in this matrix:

0.5000 0.3645 0.0000 0.2135 0.5000 0.0000 0.0000 0.2135 0.0000 0.5000 0.3178 0.0000 0.7865 0.0000 1.0000 0.2764 0.0000 0.7236 0.0000 0.3178 1.0000 0.0000 0.5000 0.0000 0.7236 0.7865 0.2764

Figure 5 is a multi-dimensional plot of the directional derivative of this
generated design ξ_{pso,D}^{3} , which is f (x)^{0}M^{−1}(ξ_{pso,D}^{3} )f (x) − 9. The plot shows
the derivative is always bounded above by 0 with equality at the support
points and so confirms the D-optimality of this design.

21

### 5 Conclusions

There are several common methods for finding optimal designs for mixture experiments. In this section, we mention some of them and discuss compar- ative advantages of PSO techniques. Early search algorithms that have been proposed for finding optimal designs are quite comprehensively reviewed in [9]. Roughly they include: Dykstra’s method [12], Wynn-Mitchell’s method [33, 41], DETMAX [32] and modified Fedorov’s method. Monographs on op- timal design usually devote a chapter on algorithmic construction of optimal designs.

Statistical software packages like SAS and JMP typically have a few menus for finding optimal designs and some have specific codes for gen- erating optimal designs for several mixture models. The types of optimal designs available in conventional packages are usually for Scheff´e polynomial models and frequently limited to D- or A-optimal designs and sometimes also for I-optimal designs. Different packages employ different methods for finding optimal design; for instance, SAS uses the exchange coordinate type algorithms which requires a candidate point set be pre-specified, and JMP uses the candidate-free exchange algorithm for finding the optimal design.

When a model is not available, it is not always clear how to generate the desired optimal design or if such an option is even possible. For instance, we were unable to find a statistical package capable of generating the D-optimal designs in Table 4 directly.

There is a R package called AlgDesign that generates optimal designs for several mixture models. This package uses the Federov exchange algo- rithm under the name optFederov, and claims to calculate the approximate designs for D-, A- and I-criteria. Optimal designs for mixture experiments are obtained using the function gen.mixture and details are freely available at http://cran.r-project.org/web/packages/AlgDesign/AlgDesign.pdf. The function “optFederov” in AlgDesign quits when no more profitable exchanges are possible.

We now compare PSO performance with AlgDesign which seems like the
most appropriate program to use since it finds optimal continuous designs
for different models and criteria in a mixture experimental problem. Results
found from AlgDeisgn and our ProjPSO algorithm were basically the same
but we observed optimal designs found from the latter are sometimes slightly
better in terms of the criterion value. For example, for the full cubic model
with three factors, the optimal design ξAD−D found by AlgDesign has 33 de-
sign points whereas the one found by ProjPSO ξ_{P SO−D}has 10 points. The rel-
ative D-efficiency of the two designs is {det(M (ξ_{AD−D}))/det(M (ξ_{P SO−D}))}^{1/10} =
0.9985. As another example, for the quadratic model with four factors, Al-

gDesign produced ξ_{AD−A}, a 25-point A-optimal design and ProjPSO pro-
duced a design ξ_{P SO−A} with only 10 points. The relative A-efficiency of the
two designs is trace M (ξ_{P SO−A})^{−1}/trace M (ξ_{AD−A})^{−1} = 0.9668. In either
case, the PSO-generated design wins. A possible explanation for the discrep-
ancy in the results is that AlgDesign does not use an equivalence theorem to
verify optimality or incorporate it as part of the stopping criterion.

To implement optFederov, AlgDesign requires a candidate set of points to be pre-specified first. The grid set we chose was 100 levels uniformly spread out for each factor. This common requirement in AlgDesign and several other packages means that the optimal design found depends on this initial grid set of points. PSO works on the a continuous domain and differentiates itself from this and other algorithms by not requiring the user to specify a candidate set of points at the onset. We view this feature of PSO a distinct advantage over its competitors.

Other advantages of our ProjPSO algorithm over current methods are (1) our experience is that the time required to generate the optimal design is gen- erally a lot faster than many of the current methods; a detailed comparison will be reported in another paper; (2) we have applied ProjPSO successfully to find optimal designs for models not covered in the standard software pack- ages, suggesting that it is versatile and can generate optimal designs for a broader range of models; for instance, ProjPSO finds the I-optimal design for the cubic Scheff´e model with three factors quickly and Figure 2 confirms its I-optimality. We were not able to find current packages that will produce such an optimal design; (3) in the few cases we had examined, ProjPSO always produced an optimal design with fewer points and higher efficiencies compared with other methods; this can be an advantage if taking observa- tions at a new “site” is expensive or laborious; (4) the ProjPSO algorithm is easy to build after downloading the basic PSO codes from various websites and then making changes to them; it is also freely available to interested reader by writing to Dr. Chen, and (5) many current methods proceed by adding a point to the current design sequentially, necessitating the user to collapse and redistribute the accumulated points and the cumulative weights to selected sites as the iteration progresses. PSO does not require such an ad-hoc procedure, which can be cumbersome.

PSO also compares favorably with other heuristic algorithms in one major aspect. For example, in genetic algorithms (GA), all tuning parameters have to be carefully selected before the algorithm works well. Our experience to date is that tuning parameters in PSO seems easy to use and are not sensitive for solving optimal design problems. Frequently, the default values for the tuning parameters work well for a broad range of design problems.

For instance, we always set c_{1} = c_{2} = 2 in our work here and when we applied
23

PSO to find minimax optimal designs for nonlinear models in [6]. The only
exception we encountered so far is that for log contrast model where setting
c_{1} = c_{2} = 0.5 seems to work better than setting them equal to 2. Our
experience this far suggests that for finding an optimal continuous design,
only two parameters in the PSO algorithm may require changes; the flock
size and the number of iterations. A larger size of randomly generated flock
of birds covers a broader range of the search space and so is suitable for
more complex and high dimensional problems. A larger number of iterations
minimizes the chance of early termination and allows PSO additional time
to find the optimum, which it usually does not need for solving our design
problems. Our typical value for a flock size is 256 and 512 if the model is
more complex. A typical iteration number that we used is 300. In addition, a
further advantage of PSO over GA is that there are fewer parameters in PSO
to adjust than in GA (http://www.swarmintelligence.org/tutorials.php).

The numerous optimal designs found here on a regular or irregular sim- plex for mixture experiments and our earlier success using PSO to find min- imax optimal designs for several nonlinear models further reinforce the great potential of PSO as a novel, easy and powerful way to generate optimal ex- perimental designs. An immediate application is to follow up work in [43] to further modify the ProjPSO to search for multiple-objective optimal designs for mixture models. On a longer term goal, our aim is to create an interac- tive web site for practitioners to generate tailor-made optimal designs where users can select the optimality criterion and input a fully specify the mixture model for their problem and let ProjPSO does its job.

### Acknowledgements

This manuscript originated at The Sir Isaac Newton Institute at Cambridge, England when Wong was a visiting fellow and a scientific adviser for a six month workshop in the design and analysis of experiment at the Institute.

He would like to thank the Institute for the support during his repeated visits there in the second half of 2011. The research of Chen and Wang were partially supported by the National Science Council of Taiwan. Wang was partially supported by the Taida Institute of Mathematical Sciences. Chen was also partially supported by the Mathematics Division of the National Center for Theoretical Sciences (South) in Taiwan.

### References

[1] J. Aitchison and J. Bacon-Shone. Log contrast models for experiments with mixtures. Biometrika, 71:323–330, 1984.

[2] N. G. Becker. Models for response of a mixture. J. Roy. Statist. Soc.

Ser. B, 31(2):107–112, 1968.

[3] N. G. Becker. Models and designs for experiments with mixtures. Aus- tral. J. Statist., 20(3):195–208, 1978.

[4] L. Y. Chan. Optimal design for a linear log contrast model for experi- ments with mixtures. J. Statist. Plann. Inf., 20:105–113, 1988.

[5] L. Y. Chan. Optimal design for experiments with mixtures: a survey.

Commun. Statist.-Theor. Meth., 29:342–373, 2000.

[6] R. B. Chen, S. P. Chang, W. Wang, and W. K. Wong. Optimal ex- perimental designs via particle swarm optimization methods. Preprint, 2012.

[7] L. Choisnard, G´eze A., Bigan M., J. L. Putaux, and Wouessidjewe D.

Efficient size control of amphiphilic cyclodextrin nanoparticles through a statistical mixture design methodology. J. Pharm Parmaceut Sci, 8:593–600, 2005.

[8] J. M. Cobby, P. F. Chapman, and D. J. Pike. Design of experiments for estmating inverse quadratic polynomial responses. Biometricsi, 42:659–

664, 1986.

[9] R. D. Cook and C. J. Nachtsheim. A comparison of algorithms for constructing exact D-optimal designs. Technometrics, 22:315–323, 1980.

[10] J. A. Cornell. Experiments with Mixtures: Designs, Models, and the Analysis of Mixture Data (2nd ed.). Wiley, New York, 1990.

[11] R. B. Crosier. The geometry of constrained mixture experiments. Tech- nometrics, 28(2):95–102, 1986.

[12] O. Dykstra. The augmentation of experimental data to maximize

|X^{0}X|^{−1}. Technometrics, 13:682–688, 1971.

[13] Y. El-Malah, S. Nazzal, and N. M. Khanfar. D-optimal mixture design:

optimization of ternary matrix blends for controlled zero-order drug re- lease from oral dosage forms. Drug Dev Ind Pharm., 32:1027–1218, 2006.

25

[14] Z. Galil and J. Kiefer. Comparison of simplex designs for quadratic mixture models. Technometrics, 19(4):445–453, 1977.

[15] Q. He and Y. Guan. Note on simplex-centroid design of degree 3: D- optimality. J. Northeast Uni. Technology, 15(5):504–507, 1990.

[16] M. K. Huang and M. N. L. Huang. D_{s}-optimal designs for quadratic log
contrast model for experiments with mixtures. Commun. Statist.-Theor.

Meth., 38:1607–1621, 2009.

[17] M. K. Huang and M. N. L. Huang. φ_{p}-optimal designs for a linear log
contrast model for experiments with mixtures. Metrika, 70:239–256,
2009.

[18] M. N. L. Huang, H. L. Hsu, C. J. Chou, and T. Klein. Model-robust D- and A-optimal designs for mixture experiments. Statistica Sinica, 19:1055–1075, 2009.

[19] X. Jin, Y. Zhang, L. Xiao, and Z. Zhao. Optimization of extended zero- order release gliclazide tablets using D-optimal mixture design. Yaku- gaku Zasshi, 182:1475–1483, 2008.

[20] O. G. Kasatkin. On the construction of D-optimal design on a sim- plex (in russian). In Application of Mathematical Methods for Muti- component Systems Investigation, pages 43–51. Metallurgia, Moscow, 1974.

[21] J. Kennedy. The Particle Swarm: Social Adaptation of Knowledge. In Evolutionary Computation, 1997., IEEE International Conference on, pages 303–308. IEEE, 1997.

[22] J. Kennedy and R. Eberhart. Particle Swarm Optimization. In Neu- ral Networks, 1995. Proceedings., IEEE International Conference on, volume 4, pages 1942–1948. IEEE, 1995.

[23] A. I. Khuri. Linear model methodology. Chapman and Hall, London, 2009.

[24] J. Kiefer. Optimal designs in regression problems II. Annals of Mathe- matical Statistics, 32:298–325, 1961.

[25] J. Kiefer. General equivalence theory for optimum designs (approximate theory). Annals of Statistics, 2(5):849–879, 1974.

[26] J. Kiefer. Optimal design: Variation in structure and performance under change of criterion. Biometrika, 62(2):277–288, 1975.

[27] J. Kiefer and J. Wolfowitz. The equivalence of two extremum problems.

Canad. J. Math., 12:363–366, 1960.

[28] S. Liu and H. Neudecker. A V -optimal design for scheff´e polynomial model. Stat. &. Prob. Letters, 23:253–258, 1995.

[29] N. K. Mandal and M. Pal. Optimum mixture design using deficiency criterion. Communications in Statistics-Theory and Methods, 37:1565–

1575, 2008.

[30] F. Mikaeili. D-optimum design for cubic without 3-way effect on the simplex. Journal of Statistical Planning and Inference, 21:107–115, 1989.

[31] F. Mikaeili. D-optimum design for full cubic on q-simplex. Journal of Statistical Planning and Inference, 35:121–130, 1993.

[32] T. J. Mitchell. An algorithm for the construction of D-optimal experi- mental designs. Technometrics, 16:203–211, 1974.

[33] T. J. Mitchell and JR. Miller, F. L. Use of design repair to construct de- signs for special linear models. In Math. Div. Ann. Progr. Rept. (ORNL- 4661), pages 130–131. Oak Ridge National Laboratory, Oak Ridge, TN., 1970.

[34] T. Nahata and T. R. Saini. D-optimal designing and optimization of long acting microsphere-based injectable formulation of aripiprazole. Drug Dev Ind Pharm., 34:668–675, 2008.

[35] T. Nahata and T. R. Saini. Formulation optimization of long-acting depot injection of aripiprazole by using d-optimal mixture design. PDA J Pharm Sci Technol., 63:113–122, 2009.

[36] M. Pal and N. K. Mandal. Minimax designs for optimum mixtures.

Statistics & Probability Letters, 78:608–615, 2008.

[37] A. Pazman. Foundations of Optimum Experimental Design. D. Reidel Publishing Company, 1986.

[38] H. Scheff´e. Experiments with mixtures. J. Royal Statist. Soc., Ser. B, 20:344–360, 1958.

27

[39] N. N. Sobolev and T. A. Chemleva. Construction of simplicial-linear mdoels of composition-property diagrams. Industrial Laboratory USSR, 42:103–108, 1976.

[40] H. Uranisi. Optimal design for the special cubic regression on the q- simplex. In Mathematical Reports, volume 1, pages 7–12. General Edu- cation Department, Kyushu University, 1964.

[41] H. P. Wynn. Results in the theory and construction of D-optimum experimental designs. J. Roy. Statis. Soc. Ser. B., 34:133–147, 1972.

[42] G. Yu and Y. Guan. A-optimal design of parameter estimation for mixture models of two degree. Journal of Northeast University of Tech- nology, 14(3):307–311, 1993.

[43] C. Zhang, W. K. Wong, and H. Peng. Dual-objective optimal designs for mixture experiments. Australian and New Zealand Journal of Statistics, 2012. In press.