{'fraction_non_alphanumeric': 0.03544015487985423, 'fraction_numerical': 0.021295979956724748, 'mean_word_length': 4.323229873908827, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Design of Experiments for Nonlinear Response Surface Models

Many chemical and biological experiments involve multiple treatment factors and often it is convenient to fit a nonlinear model in these factors. This nonlinear model can be mechanistic, empirical or a hybrid of the two. Motivated by experiments in chemical engineering, we focus on D-optimal design for multifactor nonlinear response surfaces in general. In order to find and study optimal designs, we first implement conventional point and coordinate exchange algorithms. Next, we develop a novel multiphase optimisation method to construct D-optimal designs with improved properties. The benefits of this method are demonstrated by application to two experiments involving nonlinear regression models. The designs obtained are shown to be considerably more informative than designs obtained using traditional design optimality algorithms.


Nonlinear Multifactor Models
In chemical and biological studies, experimenters often wish to explore mechanisms relating controllable input variables (treatment factors) to observed outputs. After collecting data, they fit a statistical model describing the relationship between the levels of the treatment factors and the experimental responses. In this paper, we are interested in the estimation of the model parameters, in experiments in which a model is known or assumed at the design stage. A statistical model can be either mechanistic (theoretical) or empirical and illustrative examples of each are given in Section 2. Unlike common empirical models, mechanistic models tend to be nonlinear but more frugal in the use of parameters. For some discussion of mechanistic studies, see Box and Hunter (1965). The mechanistic model is often derived from scientific theories and formulas, though it is generally still an approximation to the true response function. It often allows more scientifically informative interpretations than some empirical models.
Although less often used, these benefits apply even when the unknown mechanism involves multiple treatment factors. In the event that there are multiple input factors, however, deriving a full mechanistic model can be very hard and therefore no such model might be available at the start of an experiment. Under these circumstances, it is advisable to resort to empirical modelling techniques, while still taking into account the known parts of the mechanism when building the model. These kinds of models can be referred to as hybrid models, incorporating both mechanistic and empirical features. Particularly for mechanistic models, it is recommended to use model-oriented optimal design of experiments, since no standard designs exist for such models. This approach requires advanced optimisation methods (Chaloner and Larntz, 1989;Gotwalt et al., 2009;Gilmour and Trinca, 2012;Goos and Mylona, 2018;Overstall and Woods, 2017). Little work has been done on designing multifactor experiments for hybrid models or multifactor nonlinear models in general. This is the focus of this paper. Optimal design for a few specific classes of nonlinear models have been studied in detail. In particular, some early work exists on two-factor inverse polynomial response surface models for agricultural experiments (Mead and Pike, 1975), as well as several studies with generalised linear models, e.g. Woods et al. (2006); Stufken and Yang (2012), including discrete choice models (Ruseckaite et al., 2017) and recent pharmacokinetic studies with nonlinear mixed models (Bogacka et al., 2017).
In Section 2, two experiments of the type which motivated the work presented are introduced. In Section 3, commonly used existing methods for constructing optimal designs are described. These methods are then built on in Section 4 to describe a novel continuous optimisation algorithm for finding optimal designs. The application of the new algorithm to the illustrative applications is described in Section 5. In Section 6, our continuous optimisation method is combined with traditional ones, to obtain a computationally efficient optimisation which produces the best known designs for these problems. This combined approach is applied to the two motivating applications in Section 7. Final comments and recommendations are made in Section 8.

Illustrative Applications
In order to see the need for improved design, we use two experiments which were actually run using central composite designs (CCDs) and are typical of the types of experiment performed in chemical engineering studies. In both cases, multifactor nonlinear models were useful for describing the data, though the experiments were not designed to be optimal for these models. The CCDs actually used give a benchmark against which improved designs can be compared.

Example 1: Multifactor Mechanistic Model
Our first motivating example is described in Box and Draper (1987, chapters 11-12), where both an empirical full second-order linear model and a mechanistic nonlinear model were fitted to the same data. A continuous stirred reactor was operated, for which the chemists wished to approximate the response surface and interpret the mechanism it may imply. The assumed mechanism, which was derived after the experiment was conducted, describes a two-step consecutive decomposition reaction of a chemical solution. The yield η (%) is the amount of the desirable product that has been formed, which is the primary response. The experiment involved three input factors: the rate R (L/min) at which the chemical flows into the reactor, the catalyst concentration C and the temperature ( • C). Yield data were obtained using a 24-run spherical CCD.
The full second-order linear model involving 10 parameters was fitted using the scaled variables The characteristics of this chemical reaction can be better understood if we fit the mechanistic model where X = 0.0028344 − 1/(T + 273) and ε is the error term. This model is derived from a number of differential equations and the Arrhenius equation for temperature effects. In this paper we will show how to plan a new experiment in order to precisely estimate the parameters θ 0 , θ ′ 0 , θ 1 , θ ′ 1 , θ 2 and θ ′ 2 , of the mechanistic model in Model (1).

Example 2: Hybrid Nonlinear Model
In the second motivating example, the experimenters' interest is in the enzymatic depolymerisation mechanism of a dextran substrate . In a stirred-cell membrane reactor, endodextranase is used as the enzyme activator, while the reactants are different kinds of oligodextrans. The treatment factors are the substrate concentration S (2.5-7.5% in weight/volume), the enzyme concentration E (0.625-62.5 Units ml −1 × S) and the transmembrane pressure P (200-400 kPa). The response we consider is the substrate conversion rate ξ (%). The data were obtained using an 18-run face-centred CCD.
No mechanistic model could be found to describe the response surface for ξ, so the alternative was to use an empirical approximation. In spite of the simplicity of deriving purely empirical polynomial models, however, it is worth considering a hybrid model which could have some mechanistic implications and so catalyse scientific reasoning. First, we can observe a negative effect from the experimental data: the higher the substrate concentration, the lower the conversion rate ξ. This suggests a smooth function of the form E(ξ/(100 − ξ)) = γ ′ 1 S/(γ ′ 2 + S), where E(ξ/(100 − ξ)) is the expectation of the transformed response and γ ′ 1 and γ ′ 2 are constants with γ ′ 1 ≥ 0 and γ ′ 2 ∈ (−2.5, 0). The parameter γ ′ 1 is then replaced with an empirical exponential function of the coded variables x E = log 10 (E/6.25) ∈ [−1, 1] and x P = (P − 300)/100 ∈ [−1, 1]. This results in the hybrid nonlinear model given by where ε is the error term and a 0 , a 1 , a 2 , a 3 , a 4 and a 5 are the six model parameters. For details of the model derivation from our conjectures, see APPENDIX A. In this paper, we show how to create optimal experimental designs for the model in Model (2).

Local D-Criterion for Nonlinear Models
To obtain an exact experimental design X , we need to choose a level for each of the v treatment factors for each of the n runs. For i = 1, 2, . . . , n, the ith response is where X i is the ith row of the n × v design matrix X and θ denotes the set of p parameters. Under ordinary least squares estimation, the errors ε i must be uncorrelated and come from identical distributions with zero mean and constant variance V(ε i ) = σ 2 .
If the function f (X i , θ) is linear in the parameters, the matrix form of the model is Y = F θ+ε. Here, Y is the column vector of the n responses, ε is the column vector of the n error terms, and F is the n × p model matrix based on the design X .
We assume the primary purpose of the experiment to be the precise estimation of the parameters of the model. This requires us to minimise the elements of the p × p variance-covariance matrix V(θ) = (F T F ) −1 σ 2 , where we assume, without loss of generality, that σ 2 = 1 when constructing designs and F T F is the Fisher information matrix. To ensure that X contains the maximal amount of information about θ, we can search for a design X that maximises the D-criterion function φ = log F T F ∈ (−∞, +∞). If φ A and φ B are the D-criterion function values of any two designs X A and X B respectively, the relative efficiency of X A with respect to To derive the D-criterion function of any nonlinear model (regardless of its functional form), the most widely used approach is to linearise f (X i , θ) in θ (see e.g. Atkinson et al. (2007, chapter 17)) by means of a first-order expansion about a selected centre θ 0 , giving where θ 1 , . . . , θ p are the elements of θ. The ith row of F is then and the corresponding information matrix is given by F T F . For a nonlinear model, φ is called the "local" D-criterion function, since its value is associated with the substitution θ = θ 0 . The choice of the prior θ 0 is generally based on similar experiments reported in the literature, the experimenters' expertise or even a "best guess" about the response surface. The design optimises φ when θ = θ 0 and usually is optimal or close to optimal when |θ − θ 0 | is small. As the difference in |θ − θ 0 | increases, the design might be suboptimal, though it is still valid in the sense that it gives asymptotically unbiased estimators of all the parameters. Because of the dependence of the Fisher information on the model, the parameter prior, and the number of experimental runs, the optimal designs for nonlinear models may be quite different from those for second-order or even higher-order linear models. Hence we need to find the optimal design for the proposed nonlinear model if feasible. To search for an exact design X comprising n×v coordinates or factor settings, we develop an algorithm that can derive and quickly integrate the local D-criterion function in an iterative search and optimisation. The algorithm should be able to work with any nonlinear parametric model, independent of model assumptions, and, up to computational limitations, any experimental size. The essential inputs to our algorithm are the model function f (X i , θ), the number of experimental runs n, the v-dimensional experimental region X and initial values of the parameters θ 0 . Fedorov (1972) built the theoretical framework for the earliest algorithm for locally D-optimal design (LDOD). The algorithm requires a discretisation of the continuous design region X , defined by the possible treatment factor levels, to obtain a set Ω of N candidate points. For instance, for a first-order linear model with three factors, two candidate levels should be defined for each factor (i.e. the maximum and minimum). In that case, there are 2 × 2 × 2 = 8 candidate points in total. A full second-order linear model demands at least one more candidate level for each factor, such that there would typically be 3 × 3 × 3 = 27 candidate points. Defining candidates is much harder for nonlinear models, since it is generally unknown what will be good candidate points. Having defined the candidate set Ω, we then apply discrete optimisation through numerous iterations, each of which can exchange at most one point of the design X with one candidate point. The more candidate points are included, the more reliable the optimisation is, but the higher the computational time.

Iterative Discrete Optimization
Despite the dependence between iterations, X can usually be globally optimised, though there is no guarantee of this. To speed up the computation, the modified Fedorov exchange algorithm (Cook and Nachtsheim, 1980) streamlines the iterative procedure of the Fedorov exchange algorithm, i.e. up to n points of X can be updated in each iteration and an exchange is executed whenever a clear improvement in the criterion function is achievable. A point exchange algorithm (PEA) we adapted for nonlinear models which, like all other algorithms we present in this paper, is implemented in the Matlab environment is as follows: 1 Generate an initial nonsingular n-run design X by means of n independent random draws (with replacement) from the candidate set Ω. After substituting θ = θ 0 , compute F and φ.
2 Set index i to 1 and Υ to 0.
3 Evaluate the relative improvement function d i for each possible point exchange between the ith row of X and one of the candidate points in Ω.
4 If the maximal d i exceeds the critical value, execute the corresponding substitution X i = X new,i in X . Update F and φ and set Υ to 1.
5 If i < n, set i to i + 1 and return to STEP 3.
6 If Υ = 1, return to STEP 2. Otherwise, save the current X and the maximal local D-criterion value φ.
8 Report the best of the designs found. That is the LDOD.
In the above algorithm, we use the indicator variable Υ, which takes the value 1 if at least one improvement has been made to X in the current iteration of the algorithm. In that case, a complete new iteration will be performed after finishing the current one. We also define the relative improvement function The iteration can take advantage of the updating function of |F T F | for the D-criterion (Fedorov, 1972). Sequentially maximising d i approximates maximising the D-criterion function φ defined in terms of X . A smaller critical value in STEP 4 would lead to more iterations and updates of X . As a result, the design X attained in STEP 6 may converge further towards (and get stuck in) a local optimum. We discuss in Section 4 how to determine the critical value.
In a PEA, each iteration is broken into n dependent steps, one for each row of X . Therefore, we should do a number of tries, i.e. run the algorithm for several random starting designs, to secure an efficient solution X d .
An alternative to a PEA is the coordinate exchange algorithm (CEA) of Meyer and Nachtsheim (1995). For discrete optimisation, this uses a unidimensional subset Ω k at each step of the iteration, which consists of the candidate coordinate levels of the kth factor. Compared with the outline above, in STEP 3, instead of updating the current row X i in one step of the iteration, at most one coordinate X ik shall be updated. Specifically, X ik is replaced by the best candidate coordinate from Ω k , for k = 1, 2, . . . , v. In the CEA, STEPS 3 and 4 of the PEA are replaced by 3a Set index k to 1.
3b The relative improvement function d i is evaluated for each possible coordinate exchange between the kth coordinate in the ith row of X and one of the candidate coordinates in Ω k .
4a If the maximal d i value exceeds the critical value, execute the corresponding substitution X ik = X new,ik in X , which maximises d i . Update F and φ and set Υ to 1.
4b If k < v, set k to k + 1 and return to STEP 3b.
Each iteration of the CEA takes vn steps to separately optimise all the coordinates of X . When the experiment involves many runs and factors, this approach reduces the computational time compared with the PEA. On the other hand, the local search for optimal coordinates is less extensive than that of the PEA, so that fewer designs are evaluated in each iteration.
For solving very complicated and high-dimensional LDOD problems, some advanced heuristic search methods, such as particle swarm optimisation (Wong et al., 2015;Phoa, 2017), have been used in the recent literature. These methods are able to learn from stochastic search and iteration, and tend to run faster than the PEA or the CEA, but often give suboptimal solutions (Blum and Roli, 2003;Sörensen, 2015;Ruseckaite et al., 2017) and can be very inefficient when users are unlucky with the choice of tuning constants. In low-dimensional cases, searching for LDODs is not too computationally expensive and we should use the PEA or the CEA for reliable solutions.

A Novel Continuous Optimisation Method
In this section, we present a new, improved, computing method which overcomes the main weaknesses of the PEA and CEA discussed in Section 3, namely that the candidate points are limited to a fixed set and must be predefined. These weaknesses can lead to lower efficiencies, especially for nonlinear multifactor experiments. It is often hard to decide on suitable discrete candidates for nonlinear models and it is generally useful to explore points outside any finite set Ω. One could use a denser Ω, with smaller meshes and narrower space between adjacent candidate points within X , but this would slow down the computation, and is not a fundamental solution to the problems of discrete optimisation.
For most experiments using linear models with sufficient numbers of runs, the traditional discrete optimisation methods work well because equally spaced levels are generally either optimal or very close to optimal. In contrast, when the model is nonlinear, the optimal factor levels are generally not equally spaced and it is difficult in advance to propose a good candidate set. Hence, a continuous optimisation of the relative improvement function d i over X in each iteration is more attractive. Below is a succinct outline of the continuous optimisation point exchange algorithm: STEPS 1-2 as in the PEA, but points of X are drawn from X instead of Ω.
3 Using the Nelder-Mead (or quasi-Newton) method, find the point in the design region X which numerically maximises d i .
STEPS 4-8 as in the PEA.
The main innovation in STEP 3 is to use a continuous optimisation search over the entire design region, rather than just a discrete set of candidate points, as is done in the traditional PEA. A second novel aspect is that algebraic computing is used to derive the updating formula for d i , based on the conventional updating formulae of Fedorov (1972). An alternative version of the new algorithm based on coordinate exchange can be constructed in a similar way, the difference being that the optimisation is then done one coordinate at a time in STEP 3 and X is updated more often in each iteration.
While the critical value in STEP 4 has little impact in the discrete optimisation over a small number of candidates, it can advantageously be set to be smaller in the continuous optimisation so as to allow for minor iterative improvements on X (i.e. a larger number of iterations on average). In both of our examples, 1.0001 is a reasonable critical value for the continuous optimisation, which allows for 4-6 iterations on average.
We choose the Nelder-Mead method (Nelder and Mead, 1965;Press, 2007) to maximise the function d i ∈ (0, +∞) and thus the D-criterion function φ at each exchange. This maximisation is subject to the boundary constraint X i ∈ X . The Nelder-Mead simplex method excels at optimising complex functions (Press, 2007). Because the continuous optimisation is still based on iterative search over the design region X , the solutions we find are guaranteed to be locally optimal only. In comparison, the optimal design algorithm in Gotwalt et al. (2009) is based on a traditional CEA combined with the less accurate nonlinear optimisation method of Brent (1973, chapter 5), which has been implemented in the commercial software package JMP. Brent's method is for one-dimensional optimisation of coordinates in X and thus infeasible for the PEA. In revisiting our first application in Section 7, we will compare our algorithm with the CEA in combination with Brent's method.
For problems in which the total number of coordinates in design X is small (e.g. vn < 10), Chaloner and Larntz (1989) introduced a method for direct optimisation of the entire design X . However, this approach is infeasible for large values of vn. Our continuous optimisation PEA and CEA can be considered to be a middle ground between this method and the traditional PEA/CEA.

Numerical Results and Comparisons
We now demonstrate how experiments can be optimally designed for nonlinear models by applying our improved PEA and CEA to the illustrative experiments in Section 2. To this end, we compare LDODs constructed using different approaches.

Example 1: Multifactor Mechanistic Model
For the experiment from Box and Draper (1987), we compare the results from our continuous PEA and CEA to those of the traditional algorithms that use discrete optimisation. The benchmark we use for the LDOD is a face-centred CCD with two replicates of the axial points and four centre points. This is a modification of the design actually used to fit a cubic region of experimentation. The design region is taken to be X = [1.5, 6] × [1, 4] × [70,90]. In APPENDIX B, additionally, we evaluate two alternative standard designs: 1) a spherical CCD with radius √ 2,   which is similar to the design used but restricted to the narrower design region; 2) a Box-Behnken design. These three types of standard designs are all commonly used and require fewer runs than a full 3 × 3 factorial. To evaluate the local D-criterion function for Model (1), the prior θ 0 is taken to be the nonlinear least squares estimateθ = {θ 0 = 5.90,θ ′ 0 = 1.15,θ 1 = 0.53,θ ′ 1 = −0.01,θ 2 = 15475,θ ′ 2 = 7489}, obtained from fitting (1) to the published experimental data. The criterion value of the reference design (i.e. the face-centred CCD) is φ = −52.7712.
Consider the full second-order polynomial model in terms of the scaled variables x 1 , x 2 and x 3 . This empirical model fits the data well, as shown in Box and Draper (1987, chapter 11). For this reason, we also computed a LDOD for that model, but this LDOD X emp , shown in Table 1, does not improve much on the reference CCD for parameter estimation in the mechanistic Model (1) (the criterion value of X emp is −51.0181).
Our interest is in the LDOD for Model (1), the mechanistic approximation to the response surface, so we would expect to find a better design than those intended for fitting the polynomial model. To ensure an efficient design, we perform τ = 100 independent tries with the critical value fixed at 1.0001. We use the factor levels in the reference CCD as the candidate coordinates, which are Ω 1 = {1.5, 3, 6}, Ω 2 = {1, 2, 4} and Ω 3 = {70, 80, 90} for our implementation of the traditional CEA. The full candidate set Ω for our implementation of the traditional PEA then consists of 3 3 = 27 points. To speed up the iterative continuous optimisation, in our modified PEA and CEA, we instructed Matlab to use parallel computation. Using both the PEA and the CEA, we find the LDOD X d in Table 2, the criterion value of which is −49.7321. With the PEA, the design in Table 2 was found in three of the 100 tries. With the CEA, the design in Table 2 was found only once. The best three designs are dissimilar and their criterion values are −49.7321, −49.7452 and −49.7452.
To compare this with the discrete optimisation over Ω, we also do 100 independent tries and use the same critical value 1.0001. In each step of the iteration, the current point (or coordinate) also acts as the initial vector (or value) in optimising the criterion function. When the continuous PEA is used, it produces the design in Table 3(a), which has a local D-criterion value of φ = −49.5528. The mean number of iterations is 4.4, so it does not take many steps to obtain a local maximum of φ. When we use the continuous CEA, the X d in Table 3(b) is found to be locally D-optimal. The maximal criterion value is −49.5573, which is worse than that of the design in Table 3(a). On average, the number of iterations is 5, so no clear difference between PEA and CEA is visible here.
When using continuous optimisation, exact replicate runs are rarer than near-replicates, as we   Table 3. This is because the design region X is continuous and implies an infinite number of feasible candidate runs. However, despite small differences in the later decimal places of the coordinates, it is clear what kind of coordinates we should use under the optimality criterion. Many of the runs are not far in distance from the candidates previously defined in the 3 × 3 × 3 set Ω. Hence the designs in Table 3 does not show much of an improvement over that in Table  2, which we found with the traditional discrete optimisation. With respect to Table 3(a) which consists of quasi-continuous coordinates (i.e. the real coordinate values are irrational numbers), the reference CCD is exp(−52.7712/6)/exp(−49.5528/6) ≈ 58.48% efficient and X emp is 78.33% efficient. However, the relative efficiency is as high as 97.06% for the design in Table 2, which is made up of candidate points from Ω. This shows that, even though the model is complex and nonlinear, the discrete optimisation can be effective provided a suitable candidate set is used.

Example 2: Hybrid Nonlinear Model
For the experiment in Mountzouris et al. (1999), we assume Model (2). To determine θ 0 , we fit that model to the dataset obtained from the 18-run face-centred CCD (  (2). As an empirical approximation to the response surface, again, one can fit a full second-order linear model in terms of the scaled factors x S , x E and x P , where It is straightforward to find the LDOD X emp for the empirical model over discrete candidates. If we evaluate X emp under Model (2), the local D-criterion function value is 34.4783, which is clearly better than that from the CCD. It can be shown that Model (2) gives a better approximation to the data than the polynomial model. The critical value 1.0001 is used in our search for an 18-run LDOD X to maximise the local D-criterion function for Model (2). For discrete optimisation with the PEA, the candidate set Ω of 3 3 points is based on the factor levels in the CCD in Table 4. The best three solutions obtained from 100 tries of that algorithm are similar. For each of them, the criterion function value is φ = 38.8433. Under the CEA, the largest three values obtained are 38.7313, 38.6514 and 38.6514, so the PEA works better than the CEA in this case.
For the continuous optimisation, each exchange in the iterative procedure starts at multiple initial points, in addition to the current design point X i . We use a coarse set of 2 3 initial points {3.75, 6.25} × {1.9764, 19.764} × {250, 350} forming a cube within the design region X halfway between the centre and the edges of the region. To maximise the criterion function in v = 3 factors, at each iteration, there are in total (1 + 2 3 ) = 9 initial starts under the PEA and three for unidimensional maximisation under the CEA. The designs found using the continuous PEA and CEA are shown in Table 5(a) and (b) respectively. We can see little difference between these two solutions, which both have a D-criterion value of 41.2246 to four decimal places.
Here the continuous optimisation results in a bigger improvement than in the previous example. Compared with the designs produced by continuous optimisation, shown in Table 5, the design obtained by discrete optimisation is exp(38.8433/6)/exp(41.2246/6) ≈ 67.24% efficient. The design for the empirical model X emp is 32.49% efficient and the reference CCD is only 20.63% efficient. This demonstrates that, if a nonlinear multifactor model is to be used, a design for this model must be sought, rather than relying on standard designs. It is also clear that the greatest benefit can only be realised by using a continuous optimisation algorithm.

A Novel Multiphase Optimisation Algorithm
The continuous optimisation algorithm demands many tries to find a highly efficient design. This is due to the fact that random sampling to construct a starting design lacks effectiveness, because many initial designs cannot lead to D-optimal designs. In order to improve the effectiveness of the search, we introduce a refinement to the new continuous optimisation method. In most situations, there is insufficient information about the appropriate candidate points that would eventually compose the LDOD. However, it can be expected that the more efficient are the initial designs, the more efficient the final design will be under the D-criterion.
A new idea is to perform a swift discrete optimisation, over a small, coarse, candidate set, to convert each random initial design into an intermediate solution X d1 , to obtain a reasonably efficient starting design for the continuous optmisation. To ensure a fast generation of the interim solutions in Phase 2 and to ensure that we get several different intermediate designs, we increase the critical value for this optimisation, e.g. to 1.1. This is Phase 1 of the multiphase algorithm. Next, the distinct X d1 shall act as new starting designs for the continuous optimisation in Phase 2, the computation of which is more intensive.
Practical limitations will influence how we can set the levels of the factors. To take account of them, in Phase 3, the closest distance of a factor indicates the minimum space between two feasible coordinates next to each other, the levels of which must be distinguishable to experimenters. As we have set the v closest distances, the continuous variable space X is redefined to be a discrete set of candidate levels. We then adjust the continuous design to use only levels from this discrete set. Our multiphase optimisation PEA is as follows: 1.1 Phase 1: Compute the LDOD over discrete set Ω, using STEPS 1-7 in the traditional PEA.

Refer to the intermediate solutions obtained in
Phase 1 as X d1 .
2.1 Phase 2: Compute the LDOD over X , using either the continuous PEA or the contiuous CEA, for each X d1 acting as the starting design X .
2.2 Refer to the revised solutions as X d2 , the most efficient of which is selected to be X d .
3.1 Phase 3: Find the LDOD for refinement. For k = 1, 2, . . . , v, define the closest distance for the kth factor. Points with all factors within the closest distance will be considered quasi-replicates. Group the quasi-replicates in X d , so that the n runs of X = X d are divided into n * homogeneous clusters. Sort the rows of X with respect to the clusters and factor levels. Reset k to 1.
3.2 Set both i and i * to 1 as the starting values.
3.3 Pick out the i * th cluster of X , corresponding to n i * quasi-replicate points, which will be replaced by true replicates. Let the maximum value in the cluster of the kth factor be X max,k and the minimum be X min,k .
3.4 Create a provisional subset of candidate coordinates Ω k , on the basis of the closest distance, the interval [X min,k , X max,k ], and X . Search over Ω k and let X new,k be the candidate that maximises φ. Use this substitute for each point of the i * th cluster. Set i to i + n i * . After substitution, there are n i * identical values of the new coordinate.
3.5 Unless i = n + 1, set i * to i * + 1 and return to STEP 3.3.
3.6 Unless k = v, set k to k + 1 and return to STEP 3.2.
3.7 Create a new candidate set Ω, consisting of the n * distinct design points of the current X . Perform a discrete optimisation over Ω, the final solution X * d of which is the LDOD.
In this improved algorithm, STEP 3.3 creates the variable space for factor k in the i * th cluster, and then STEP 3.4 uses a variant of coordinate exchange to optimise the common factor level. The aim is to find the best level for the kth factor. Suppose the closest distance is 0.01 unit for the flow rate in the Box and Draper example. In Table 3b, the 13th and 14th runs are quasi-replicates. The corresponding two flow rates are within the interval [3.4081, 3.4102] ⊂ [3.40, 3.42]. To select their factor levels as in STEP 3.4, the candidate subset should be Ω 1 = {3.40, 3.41, 3.42}. In the final solution X * d , n * is equal to the number of distinct design points and n i * indicates the number of exact replicates at each point, for i * = 1, 2, . . . , n * .
The CEA can improve the continuous optimisation if the number of factors v is large. Otherwise, the PEA is found to be more reliable in our examples.
In the multiphase algorithm, the new continuous optimisation method is used to complement the basic discrete optimisation over an imperfect candidate set. Phase 3 aims to convert the quasi-replicate runs of X d into exact replicates. Donev and Atkinson (1988) proposed another algorithm for adjusting an already efficient design, though this is not specifically to create replicates.

Example 1: Multifactor Mechanistic Model
For the Box and Draper example, the new multiphase optimisation algorithm can start with the same 3 × 3 × 3 candidate set Ω, based on the factor levels from the reference CCD. This candidate set is suitable for an initial discrete optimisation and also facilitates Phases 2 and 3. With 30 random tries and a critical value of 1.1 in Phase 1 of the multiphase algorithm, 30 distinct intermediate solutions X d1 are obtained using the traditional PEA. Next, in Phase 2, the continuous optimisation starts from each of the designs in that set using a critical value 1.0001 in Phase 2. Typical interim solutions after Phase 2 are shown in Table 6 With the CEA, we also tested Brent's unidimensional method instead of the Nelder-Mead method. The CEA takes the maximum of 30 iterations in all but one case. The mean criterion value is −50.0487 and the largest three values are −49.8285, −49.8633 and −49.8766. Brent's method is faster for a single iteration, but it fails to deliver a stable and efficient solution after Phase 2. In our demonstration on a Windows desktop computer with an Intel Core i7 Processor, the total required computing time for 30 tries is 8 seconds with the PEA, 10 with the CEA using the Nelder-Mead algorithm, and 15 with the CEA using Brent's method. Thus, our improved algorithm takes just over half of the computing time of the best existing method, though in this case the difference is practically unimportant.
We can choose either interim solution X d in Table 6 for Phase 3. The closest distance is set to be 0.1 unit for R, 0.1 unit for C and 1 unit for T . After the reallocation of experimental runs described in STEP 3.7 (i.e. recalculating the optimal numbers of replicates for each distinct design point), the final solution X * d in Table 7 is obtained. This design has only eight  Table 7: 24-Run LDOD for Model (1) produced by the multiphase algorithm distinct design points, to estimate the eight parameters. The distinct design points are unequally replicated. We also note that only two levels of the catalyst concentration are used. The new maximal criterion value is −49.5116 which is better than all values obtained previously. More importantly, now we have fully achieved the objective to accurately identify the factor levels of every distinct design point (at desirable precision) as well as the corresponding number of replicates. A standard experimental design (e.g. CCD or Box-Behnken design) requires 13-15 distinct design points here, as illustrated in APPENDIX B, while a full factorial design requires at least 3 3 = 27 distinct points. By using the multiphase optimisation algorithm, in addition to achieving higher efficiency and more reliable parameter estimation, we are able to simplify the initially constructed design, making it more practical and more efficient. In this example, there are only n * = 8 support points in X * d , in contrast to the 11 distinct design points for the initial design in Table 2 which was found using the discrete optimisation.
As a side note, even if we use Brent's method for CEA (i.e. unidimensional continuous optimisation), it is feasible to obtain Table 7. We can, for instance, increase the maximum number of iterations and reduce the critical value (currently at 1.0001) to allow for more optimisations and achieve a stable D-optimality. More effectively, we can carry out STEP 3.7 of Phase 3 for every Phase 2 design X d2 , as suggested by Gotwalt et al. (2009) and implemented in JMP. Using JMP 13, with 100 tries using Brent's method, we found a design which is only slightly worse than that in Table 7 and the same up to rounding of factor levels, though it takes longer to find.

Example 2: Hybrid Nonlinear Model
For the Mountzouris et al. example, after the initial discrete optimisation (Phase 1), there are just four distinct designs X d1 that can be used as inputs for Phase 2, each of which is made up of 18 runs selected from the 24-run candidate set Ω. At the end of Phase 2, with the PEA,  Table 8, the criterion value of which is also 41.2246. This is similar to the D-criterion value of the design in Table 5, so the design does not improve the D-efficiency, but it simplifies the search without sacrificing D-efficiency. The new design has nine distinct design points, seven of which are replicated twice and one is replicated three times.

Discussion and Recommendations
Empirical polynomial response surfaces are widely used and good designs for estimating them are readily available. However, as noted in Box and Draper (1987), experimenters might benefit from deriving a mechanistic model or an empirical nonlinear model, which can better approximate the response surface and facilitate the scientific interpretation of results. The optimal design of multifactor experiments for estimating nonlinear models has been rather neglected and the traditional design method via discrete optimisation is not always successful. To solve this problem, we have developed a multiphase optimisation algorithm. Searching for such an optimal design is no longer challenging with the multiphase optimisation method. When this algorithm is applied to both illustrative examples, we easily obtain a locally D-optimal design. These locally D-optimal designs have a simple structure that is different from those of standard experimental designs. It would therefore be interesting to apply the methodology in this paper to study different classes of nonlinear experiments.
To promote such designs, there are two challenges in practice: (i) to find a suitable form of the nonlinear model; and (ii) to choose a set of realistic values as the parameter prior. On the first challenge, if we cannot assume a mechanistic model, often it is still possible to use a hybrid nonlinear model reflecting some characteristics of the underlying mechanism. This is demonstrated in Example 2. As to the second challenge, we should base the locally D-optimal design on the most realistic parameter values we are able to assume. Some reference data, either from the literature or earlier experiments, are helpful. Otherwise, values obtained from the literature on related systems can be useful, or in the absence of anything else experimenters have to use their knowledge, experience and intuition to suggest prior values. It is important to note that, even if these values are inaccurate, the design chosen will still be valid and usually quite efficient, even though it will not in general be optimal for the unknown true parameter values.
Our continuous optimisation method is robust in updating the factor levels in X , as we demonstrate in both examples. It takes the whole design space X into account, rather than a finite candidate set Ω. Hence, we avoid the risk of using an improperly specified candidate set. Further, the multiphase optimisation algorithm circumvents ineffective tries and iterations, in order to find a locally D-optimal design with less computational effort.
The most important message from our work is that, if a multifactor nonlinear model is to be used to gain scientific insight from experimental data analysis, it is important to tailor the experimental design to that model. We have provided the tools needed to do this effectively and we recommend their use to experimenters in chemical engineering, biochemistry, and related areas. The enormous potential to harvest data which are informative about plausible mechanisms, rather than just giving rough predictions more than compensates for the extra efforts taken to acquire designs.
Under an extreme scenario when even the maximal substrate concentration S = 7.5 (% weight/volume) is far too low (relative to the fixed concentration of endodextranase), the reaction rate will be high and stable over time. The conversion is then close to 100% at the measurement time (after the reaction almost ceases), when the instantaneous substrate concentration is near zero and treated as independent of the initial concentration. The expected conversion rate is E(ξ) → 100(S − a)/S, where a is a nonnegative constant near zero. In contrast, if even the minimal substrate concentration level 2.5 is too high for the concentration of endodextranase, the reaction will be slow and its rate can decrease quickly over time. The concentration of the converted substrate is almost independent of the initial substrate concentration, so that E(ξ) → 100b/S, where b is a nonnegative constant.
A tradeoff must be found between these two impractical scenarios, where we assume an ideal solution (or ideal mixture). When the substrate concentration increases, ξ shall decrease. It is also more realistic to consider the function of ξ to be concave rather than convex. As an approximation to the observed response surface in Mountzouris et al. (1999), we consider a nonlinear function in parameters γ 1 ≥ 0 and γ 2 ∈ (−2.5, 0): E(ξ i ) = γ 1 S i γ 2 + S i , for i = 1, 2, . . . , n.
To use a transformed dependent variable, we derive the statistical model which incorporates an error ε ′ i and two constants γ ′ 1 ≥ 0 and γ ′ 2 ∈ (−2.5, 0). When the substrate concentration is fixed, we use a second-order linear function to explain the transformed conversion rate in the scaled variables x E and x P . With interaction x E x P excluded because of insignificance in this case, the final model is where ε ′′ i is the error and a ′ 0 , a ′ 1 , a ′ 2 , a ′ 3 and a ′ 4 are unknown parameters. We combine (4) and (5) to write the nonlinear multifactor model as ξ i 100 − ξ i = exp(a 0 + a 1 x E,i + a 2 x P,i + a 3 x 2 E,i + a 4 x 2 P,i )S i a 5 + S i + ε i , for i = 1, 2, . . . , n, where ε i is the error and a 0 , a 1 , a 2 , a 3 , a 4 and a 5 are the parameters. This model fits the reference data well, so we can use it to approximate the unknown response surface.

B Standard Designs for Example 1
Tables B1-B3 below show three standard experimental designs we may construct for Example 1, each of which consists of 24 independent runs. For each design, we define four exact replicate runs at the centre (3, 2, 80) of the variable space X . These designs are used in this paper as benchmarks in evaluating the relative efficiency of the LDODs found using different algorithms in Sections 5 and 7. Particularly, the face-centred CCD is the reference design we choose and study in Section 5, the local D-criterion value of which is −52.7712. Note that the spherical CCD (Table B2) is even less efficient in this case, with a criterion value of −54.6880. This is because of the suboptimal factor levels it uses. The modified Box-Behnken design (BBD) ( Table  Table B1: 24-Run Reference Face-Centred CCD R C T R C T R C T R C T 1.5 1 70 3 1 80 3 2 80 6 1 70