MP-CE Method for Space-Filling Design in Constrained Space with Multiple Types of Factors

: Space-ﬁlling design selects points uniformly in the experimental space, bringing consider-able ﬂexibility to the complex-model-based and model-free data analysis. At present, space-ﬁlling designs mostly focus on regular spaces and continuous factors, with a lack of studies into the discrete factors and the constraints among factors. Most of the existing experimental design methods for qualitative factors are not applicable for discrete factors, since they ignore the potential order or spa-tial distance between discrete factors. This paper proposes a space-ﬁlling method, called maximum projection coordinate-exchange (MP-CE), taking into account both the diversity of factor types and the complexity of factor constraints. Speciﬁcally, the maximum projection criterion and distance criterion are introduced to capture the “bad” coordinates, and the coordinate-exchange and the optimization of experimental design are realized by solving one-dimensional constrained optimization problem. Meanwhile, by adding iterative perturbations to the traditional coordinate exchange process, the adjacent areas of the local optimal solution are explored and the optimum performances of the current optimal solution are retained, while the shortcomings of random restart are effectively avoided. Experiments in the regular space and constraint space, as well as experimental design for the terminal interception effectiveness of a missile defense system, show that the MP-CE method signiﬁcantly outperforms existing popular space-ﬁlling design methods in terms of space-projection properties, while yielding comparable or superior space-ﬁlling properties.


Introduction
The experimental design aims to strategically select a combination of experimental factors, so as to maximize benefits at the minimum computational cost. Space-filling designs are used to observe the response by placing a set of experimental points uniformly in the design space, which can provide insight into the performance of the experimental object in the whole experimental space and has high flexibility in exploring the relationship between the factors and the response. Typical space-filling design methods include Latin hypercube design (LHD) and uniform design. LHD is the standard method for uniformly placing points on each univariate dimension [1]. Johnson, Moore, and Ylvisaker proposed a method to optimize LHD using space-filling criterion, including the Maximin LHD (MmLHD) that maximizes the minimum distance between two points, and the Minimax LHD that minimizes the maximum distance between two points [2]. Kaitai Fang and Yuan Wang put forward a uniform design to find the test points uniformly distributed in the design space [3]. Joseph believes that the computer experiment output is deterministic, and any repeated design points will lead to a waste of computing resources. Therefore, a good design of computer experiment should be non-overlapping. LHD improves onedimensional projection properties but cannot guarantee good space-filling properties in higher-dimensional subspaces. Therefore, Joseph, Gul and Ba et al. proposed the maximum projection (MaxPro) design, which can simultaneously optimize the space-filling properties 2 of 13 of the design points with respect to all possible subsets of factors [4]. The test proves that it also has a good performance in the space-filling distance.
The current experimental design mainly focuses on the continuous numerical factors in the regular space, however, the actual experiments may involve discrete numerical factors and ordered qualitative factors, and there may be complex constraints between factors. For example, in the flight experiment of ship-to-air missile, the experimental factors include target speed, target height, target distance and so on. The target speed is generally a discrete numerical factor corresponding to the specific type of target, for example, the speed of a helicopter in a low-speed hovering state can be regarded as 0 m/s, while that of a subsonic anti-ship missile is around 300 m/s. The target height has three levels of high, medium, and low altitude, and the target distance has three levels of near, medium and far. They are both qualitative factors and their levels are in order. Furthermore, there are constraints between the target speed and the target height.
In a bid to solve the above problems of experimental design with discrete factors and qualitative factors, this paper proposes a coordinate exchange algorithm based on MaxPro criterion and distance criterion, namely the MP-CE method, which incorporates discrete factors for optimization, and calculates the probability of coordinate exchange based on distance properties and projection properties to avoid falling into a local optimal solution prematurely. Each time the local optimal solution is found, iterated perturbations are added to avoid the shortcomings of random restart of the coordinate exchange algorithm. As ordered qualitative factors can be transformed into discrete numerical factors, discrete numerical factors and ordered qualitative factors are collectively referred to as discrete factors in the following description.

Experimental Design with Discrete Factors
The space-filling designs with qualitative factors include the slice Latin hypercube design [5][6][7], the marginally coupled design [8,9], the clustering filling design and so on [10,11]. The optimal design strategy for experiments with discrete factors differs from that for experiments with continuous factors, as the level of discrete factors may be repeated, and it also differs from that for experiments with qualitative factors, because the difference between qualitative factors cannot be measured by distance. However, the current main strategy for experimental designs with discrete factors is to ignore the order between factors and treat them as disordered qualitative factors. Another strategy is to treat the discrete numerical factors as continuous factors and then round each of them to the nearest discrete value. The problem with this method is that the final design may not be optimal, or even worse than the local optimal solution. Joseph first extends its MaxPro design for experiments with continuous, discrete, and qualitative factors, which can achieve good space-filling properties in the whole design space and all possible sub dimensional projections [12]. The design matrix is randomly initialized by generating a random LHD and folding each value of the factors to the closest given discrete level. and the simulated annealing algorithm is used to iteratively search the optimal solution by randomly exchanging two coordinates in a column of the design matrix. In fact, this method does not change the distribution of discrete levels in the initial design, thus it is not suitable for experiments with constraints between factors.

Experimental Design in Constrained Space
The typical methods of constructing space-filling designs in hypercubes are usually not suitable for constrained experimental spaces. Although it is possible to create design in unconstrained space and remove unacceptable points, this may result in unexplored areas of the design space and require additional work to ensure that the number of design points matches the requirements of the experiment. Draguic, Santner, and Dean introduced the CoNcaD algorithm [13], which is used for a space-filling design in bounded nonrectangular regions, focusing on maximum and low-dimensional properties. Beal et al. (2014) improved the WSP algorithm to construct a space-filling design, which can be applied to experimental designs with specific constraints and can increase the observation density of specific regions of interest [14]. Lulu Kang (2019) provided a random coordinate exchange algorithm for constructing experimental designs with space-filling criterion in arbitrary regular or irregular constrained spaces, which improves the optimization efficiency of coordinate exchange [15]. Lekivetz and Jones (2015) proposed a fast flexible space-filling (FFF) design [16], which clusters random points in the design space to generate design. Subsequently, the FFF algorithm was extended to allow qualitative factors to be added to the design, while maintaining the flexibility of quickly creating space-filling designs in rectangular and non-rectangular spaces [10]. However, none of these methods propose an appropriate measure of space-filling properties for designs with discrete factors.

The MP-CE Method
Assuming there are p 1 continuous factors, namely X 1 , . . . , X p1 , and p 2 discrete factors, namely X p1+1 , . . . , X p , and p 1 + p 2 = p. The number of the discrete levels of the l-th discrete factor X l (l = p 1 + 1, . . . , p), is denoted as m l . There are s constraints between factors, denoted by the vector cons = (constraint 1 , . . . , constraint s ), then the value space of factors can be denoted as R p cons . Our task is to construct a design matrix D n×p of the input factors, where the i-th row x i = (x i1 , x i2 , . . . , x ip ) represents the i-th design point, so that these n design points representing n experiments are uniformly filled in the design space.
In the literature [17][18][19][20] it was found that the ideal space-filling design should possess high separation distance, low fill distance and excellent projection distance. The separation distance refers to the minimum distance between any two experimental points. Designs with high separation distance aim to maximize the following criterion: The fill distance refers to the maximum distance between a given position in the design space and the nearest test point. Designs with low filling distance aim to minimize the following criterion: max x∈ R cons p min There is almost no contradiction between the above two distance criteria. For projection distance, the MaxPro design with multiple factors proposed by Joseph can achieve good projection properties in the whole design space and all possible sub dimensional spaces. Referring to the literature [12], this method minimizes the following criterion: Expanding the above ideas, this paper proposes a new space-filling design method, namely MP-CE, which applies the MaxPro criterion and distance criterion to guide coordinate exchange, making it suitable for experimental designs with discrete factors and constraints between factors.
MP-CE selects a "bed" coordinate of the current design by applying the coordinate selection criterion, exchanges it into a better coordinate to improve the current design criterion to generate a new design, which reduces the high-dimensional optimization problem to one-dimensional optimization, and it overcomes the problem that the algorithm is easy to fall into local optimization by applying iterated perturbation, as shown in Figure 1. Next, we will describe coordinate selection, coordinate exchange and iterated perturbation.
MP-CE selects a "bed" coordinate of the current design by applying the coordinate selection criterion, exchanges it into a better coordinate to improve the current design criterion to generate a new design, which reduces the high-dimensional optimization problem to one-dimensional optimization, and it overcomes the problem that the algorithm is easy to fall into local optimization by applying iterated perturbation, as shown in Figure  1. Next, we will describe coordinate selection, coordinate exchange and iterated perturbation.

Coordinate Selection
The coordinates to be exchanged, also known as 'bad' coordinates, are selected according to the projection criterion and distance criterion. Define the projection coefficient of the l-th dimension between two points xi and xj as: The steps of selecting coordinates are as follows: (1) From the row vector of the design matrix, [x1, x2,…, xn] T , a pair of "bad" design points (xa, xb) are randomly selected with a probability proportional to the following formula:

Coordinate Selection
The coordinates to be exchanged, also known as 'bad' coordinates, are selected according to the projection criterion and distance criterion. Define the projection coefficient of the l-th dimension between two points x i and x j as: The steps of selecting coordinates are as follows: (1) From the row vector of the design matrix, [x 1 , x 2 , . . . , x n ] T , a pair of "bad" design points (x a , x b ) are randomly selected with a probability proportional to the following formula: From the perspective of numerical stability, when p is large, the following formula can be used instead: (2) Calculate the sum of the projection coefficients of x a , x b and other design points respectively, and determine the row coordinate index i (i = a or i = b) with a probability proportional to the following formula: Mathematics 2021, 9, 3314 5 of 13 (3) Finally, in the pair of "bad" design points (x a , x b ), the column coordinate index j (j = 1, . . . , p) is determined with a probability proportional to the following formula:

Coordinate Exchange
According to the above steps, the coordinate to be improved is selected, which is x ij , and it will be optimized next. First of all, keep the remaining coordinates x il (l = j) unchanged, which belongs to point x i , and substitute the values of the remaining coordinates into the set of constraints cons = (constraint 1 , . . . , constraint s ), then solve the constraint problem for x ij , and get the updated range of x ij , which is denoted as low x ij , up x ij . Solve the one-dimensional optimization problem in the new range and use the obtained optimal solution x ij * to generate new design D* iteratively.
For each iteration to solve the one-dimensional optimization problem, the sum of the projection coefficients between other points remains unchanged except for the projection coefficients involved in point x i . The sum of the projection coefficients between the other points is expressed by the following formula.
The projection coefficients between point x i and other points remain unchanged except the j-th dimension, and the projection coefficients matrix is expressed as follows: where the column vector pro when the coordinate x ij is updated, the projection coefficients between point x i and other points in the j-th dimension will all be updated accordingly, which can be represented by the following column vector: After Formulas (10) and (11) are combined, the sum of the projection coefficients between point x i and other points is calculated, and it is updated as follows: According to Formula (3), combining Formulas (9) and (12), our objective function is to minimize Formula (13): To simplify the calculation, the one-dimensional optimization problem of coordinate exchange can be simplified as follows: For the implementation of our method, formula (14) is calculated using the R package 'nloptr' (Ypma, 2020). The entire coordinate exchange algorithm repeats the coordinateexchange optimization and compares the projection properties according to formula (13). If the new design matrix D* obtained after the coordinate exchange is better than the current optimal design matrix D, the coordinate exchange is implemented, and the current optimal design is updated. Algorithm 1 summarizes the detailed steps of coordinate exchange.  (3) i.pair ← sampl a pair of "bed" points with Equation (6) as the probability i ← sampl the "bed" point with Equation (7) as the probability j ← sampl the "bed" coordinate with Equation (8)  x ij * ← use nloptr with Equation (14) and [low xij ,up xij ] as input else x ij * ← use Equation (14) and [low xij ,up xij ] pro(x il ,x jl )* ← use Equation (12) opt(D*) ← use Equation (13) if opt(D*) < opt opt ← opt(D*) x ij ← x ij * D ← D* Output: D:design matrix opt:objective function value.

Iterated Perturbation
The local search optimization algorithm may be affected by the initial design to a certain extent. Generally, the problem of easily falling into local optimal solution is overcome by repeating the search (i.e., restarting the search with a new set of initial experimental points). Studies have found that the restart strategy merely expands the scope of search and increases the chance of avoiding local optimal solution. In addition, the strategy lacks stability and is not applicable to large-scale experimental designs [21,22]. Therefore, we add iterated perturbation to MP-CE, and set the perturbation operator as the coordinate. The iterated perturbation strategy is shown in Figure 1. Given the current optimal solution D, it is perturbed to jump out of the local area, and the perturbed solution is denoted as D'. Then the MP-CE algorithm is applied to D', and a new current optimal solution is obtained. The local search and perturbation are carried out repeatedly to return to the optimal design after the number of iterations is reached.
Through experimental research, Cuervo found that the size of the perturbation operator that leads to the best algorithm performance is PERT_SIZE = 10% [23]. Drawing reference from its experimental conclusion, this paper sets the scale of coordinate perturbation to 10% × (n × p). Algorithm 2 shows the detailed steps of iterated perturbation.

Comparison Method
We chose several popular space-filling experimental design methods for comparison. The following is a brief introduction to these comparison methods.
(1) uniform design Kaitai Fang and Yuan Wang jointly put forward the uniform design in 1978, using the uniform distribution theory in number theory to select n experimental points, and applying number theory to ensure the points were uniformly distributed within the integration range, and ensure the distribution points were sufficiently close to various values of the integrand function. The uniform design is obtained through the "space-filling: uniform design" in the DOE of JMP software.
(2) MmLHD In the LHD, each factor has as many levels as designed, and these levels are uniformly spaced from the lower limit to the upper limit of the factor. The MmLHD optimizes LHD by maximizing the minimum distance between design points. The MmLHD is obtained through the maximinLHS function in R packet "lhs" (Carnell, 2020).

(3) Fast flexible space-filling (FFF)
The FFF design (Lekivetz and Jones, 2015) [16], which is a space-filling design based on the Fast Ward clustering algorithm. First obtain a large sample U(χ), then use the Ward's minimum variance criterion (Ward Jr, 1963) for systematic clustering to form n clusters, and finally obtain the final design points by using the MaxPro optimality criterion. The FFF design is obtained through "space-filling: FFF Design: MaxPro criterion "in the DOE of JMP software.

(4) Maximum projection
The maximum projection designs with quantitative and qualitative factors (Max-ProQQ), which is a maximum projection design with quantitative and qualitative factors. First, the design matrix is randomly initialized by generating a random LHD and folding n levels to the nearest given discrete values of the factor. Use the simulated annealing algorithm (Morris and Mitchell 1995) and MaxPro criterion to optimize the initial design, and iteratively optimize in the design space by randomly selecting two coordinates in a column in the design matrix to exchange. The MaxProQQ design is obtained through the MaxProQQ function in R packet "MaxPro" (Ba, 2018).

Implementation Details
To implement the MP-CE method, we set the number of internal loop searches to 100 and the number of external perturbations to 20. We used the continuous optimization algorithm nloptr to optimize the design matrix to find the local optimal solution of coordinate-exchange. This function can be obtained in the R package "nloptr". The scale of the coordinate perturbation is set to 10% × (n × p).
For the initial design matrix, the continuous factors are uniformly distributed using the LHD, and the discrete factor levels are allocated to each design point by random sampling, and the feasibility of each design point is tested according to the constraints. The comparison experiments, MmLHD and MaxProQQ, adopt the same initial design as MP-CE. The number of iterations of MmLHD and MaxProQQ is set to 1000. The number of restarts for uniform design and FFF design is set to 20.
In Section 3, we have pointed out the properties that the ideal space-filling design should have, so in the experimental analysis of this paper, we use the MaxPro index to evaluate the space projection of the design, see formula (3), and use the index φ p based on distance to evaluate the space-filling. The ∅ p is defined as follows: where p is a positive integer and d ij represents the distance between any two design points in the design matrix D n×m , that is, the row vectors x i and x j .

Space-Filling Design in [0,1] p
In this section, we compare the uniform design, MmLHD, FFF design, MaxProQQ design, and our MP-CE design. We conducted experiments on several different input configurations, including p = 2, p = 6, p = 10, and n = 25, n = 50, n = 100. The experiments on these input configurations led to basically similar conclusions. Table 1 takes n = 50 as an example to list the performance of different experimental design methods in 2, 6 and 10-dimensional space. The results produced by the best performer and the best baseline in each column are boldfaced and underlined, respectively. First of all, it can be seen that MP-CE outperforms other space-filling design methods, i.e., uniform design, MmLHD, FFF and MaxProQQ in terms of φ p and MaxPro. When the dimension of the experimental space changes from 2 to 10, the improvements of MP-CE over the best baseline are increased. In detail, the improvements are −0.21%, 5.09%, and 3.05% in terms of φ p , 1.47%, 4.47% and 4.26% in terms of MaxPro when p = 2, p = 6, p = 10, respectively. It shows the superiority of this method in the space-filling of high-dimensional. Since MP-CE is designed based on the MaxPro criterion, the improvements of MP-CE in terms of MaxPro are larger than that of φ p . In addition, both MaxProQQ and MP-CE start with the same initial design, and MP-CE outperforms MaxProQQ in terms of φ p and MaxPro, indicating that the perturbed coordinate-exchange algorithm also has certain advantages compared with the simulated annealing algorithm. The initial design also has a certain impact on the final effect of MP-CE, especially when the number of design points is small or the dimension is low, it can be overcome by increasing the number of iterations.

Space-Filling Design in Constrained Spaces
This section will discuss the performance of the MP-CE design in constrained spaces. The experimental setting is the same as above, considering a constrained space of p = 2, Figure 2 plots the design results of n = 10 and n = 100. FFF, which performs well in the regular space and is used for comparison with MP-CE. Figure 3 compares these two methods in terms of the space projection and space-filling properties with different run sizes. As the measurements are different with different run sizes, Figure 3b,d show the performance improvement of MP-CE over FFF in terms of ratio, in a bid to make the results more intuitive. The following two points can be observed from Figures 2 and 3. On the one hand, in the constrained space of p = 2, FFF and MP-CE have very close performance in terms of space-filling and space-projection, which confirms the applicability of our method in constrained spaces. On the other hand, compared with FFF, MP-CE provides optimal maximin and maximum projection design for each run size. In particular, MP-CE performs better as the size of the design increases. certain advantages compared with the simulated annealing algorithm. The initial design also has a certain impact on the final effect of MP-CE, especially when the number of design points is small or the dimension is low, it can be overcome by increasing the number of iterations.

Space-Filling Design in Constrained Spaces
This section will discuss the performance of the MP-CE design in constrained spaces. The experimental setting is the same as above, considering a constrained space of p = 2, Figure 2 plots the design results of n = 10 and n = 100. FFF, which performs well in the regular space and is used for comparison with MP-CE. Figure 3 compares these two methods in terms of the space projection and space-filling properties with different run sizes. As the measurements are different with different run sizes, Figure 3b,d show the performance improvement of MP-CE over FFF in terms of ratio, in a bid to make the results more intuitive. The following two points can be observed from Figures 2 and 3. On the one hand, in the constrained space of p = 2, FFF and MP-CE have very close performance in terms of space-filling and space-projection, which confirms the applicability of our method in constrained spaces. On the other hand, compared with FFF, MP-CE provides optimal maximin and maximum projection design for each run size. In particular, MP-CE performs better as the size of the design increases.  For the above experimental design in a 2-dimensional constrained space, a discrete factor was added, with the value of 0, 0.5, 1, and the constraints between the factors remain unchanged.  For the above experimental design in a 2-dimensional constrained space, a discrete factor was added, with the value of 0, 0.5, 1, and the constraints between the factors remain unchanged. For the above experimental design in a 2-dimensional constrained space, a discrete factor was added, with the value of 0, 0.5, 1, and the constraints between the factors remain unchanged.
For the purpose of comparing with MP-CE method, we adopted the "custom designs" platform of JMP, which can define factor types and factor constraints. Enter two continuous factors, specify constraints on the design space, and impose constraints on the third factor to make it infinitely approximate to the level of discrete value. Then select the FFF design and set the optimality criterion to MaxPro. Set n = 27. In Figure 4, we made the 2-dimensional mapping discrete diagram to show the results from the two experimental design methods. The circle, triangle, and plus sign in the Figure represent x 3 = 0, 0.5, 1, respectively. It can be seen from the diagram that the space-filling performance of the FFF is obviously inferior to that of MP-CE.
For the purpose of comparing with MP-CE method, we adopted the "custom designs" platform of JMP, which can define factor types and factor constraints. Enter two continuous factors, specify constraints on the design space, and impose constraints on the third factor to make it infinitely approximate to the level of discrete value. Then select the FFF design and set the optimality criterion to MaxPro. Set n = 27. In Figure 4, we made the 2-dimensional mapping discrete diagram to show the results from the two experimental design methods. The circle, triangle, and plus sign in the Figure represent x3 = 0, 0.5, 1, respectively. It can be seen from the diagram that the space-filling performance of the FFF is obviously inferior to that of MP-CE.

Case Study
Sections 5.1 and 5.2 discussed the performance of MP-CE in rule space and constraint space. This section will discuss the applicability and performance of MP-CE in the constrained space with multiple types of factors, in combination with an experimental design for the terminal interception effectiveness of the missile defence system.
At the terminal of the missile-target intersection, when the missile attacks the target, the kill probability is associated with many factors, including the accuracy of the missile guidance system, the miss distance, the positional relationship when the missile and target meet, the initiation trigger performance of the warhead, the explosion power, and the anti-damage ability of the target, etc. In one word, it is the result of randomly combined multiple factors. Table 2 lists the factors that may be related to missile interception effectiveness. The assessment of missile intercepting low-speed hovering helicopter targets, subsonic aircraft targets, subsonic missile targets, supersonic maneuverable missile targets, and supersonic anti-ship missile targets correspond to the following target speeds in Table 2. Regarding the constraint between the factors, when the target speed is 0, the target height can only correspond to ultra-low altitude and hollow altitude.

Case Study
Sections 5.1 and 5.2 discussed the performance of MP-CE in rule space and constraint space. This section will discuss the applicability and performance of MP-CE in the constrained space with multiple types of factors, in combination with an experimental design for the terminal interception effectiveness of the missile defence system.
At the terminal of the missile-target intersection, when the missile attacks the target, the kill probability is associated with many factors, including the accuracy of the missile guidance system, the miss distance, the positional relationship when the missile and target meet, the initiation trigger performance of the warhead, the explosion power, and the antidamage ability of the target, etc. In one word, it is the result of randomly combined multiple factors. Table 2 lists the factors that may be related to missile interception effectiveness. The assessment of missile intercepting low-speed hovering helicopter targets, subsonic aircraft targets, subsonic missile targets, supersonic maneuverable missile targets, and supersonic anti-ship missile targets correspond to the following target speeds in Table 2. Regarding the constraint between the factors, when the target speed is 0, the target height can only correspond to ultra-low altitude and hollow altitude. First of all, the experimental factors are normalized. The encounter distance is mapped to X 1 = [0,1], the target speed is converted to discrete values X 2 = 0, 0.25, 0.3, 0.8, 1, the target height is converted to discrete values X 3 = 0, 0.7, 1, and the target maneuver overload is converted to discrete values X 4 = 0, 0.2, 0.4, 0.6, 0.8, 1. Constraint between factors is when X 2 = 0, X 3 = 0 or 0.7.
Additionally, we used the FFF design of JMP as the baseline, set the optimality criterion to MaxPro, entered 4 continuous factors, and imposed constraints on the second, third, and fourth factors to make it infinitely approximate the value of each discrete level. Taking n = 10 as an example, the final experimental design results returned by FFF and MP-CE are shown in Table 3. Experimental designs with different number of points were performed, as shown in Figure 5, leading to similar conclusions as described above. For the experimental design with multiple types of factors and complex constraints, MP-CE performs similarly to FFF in terms of space-filling distance but performs much better in terms of projection properties.

Conclusions
This paper proposes the MP-CE method to solve the problem of experimental design with multiple types of factors and complex constraints between factors. A coordinate selection strategy based on the projection property and distance property of design points is designed. By assigning probability of coordinate selection to overcome the shortcomings of randomness and easy to fall into local optimum and by adding perturbation to each local search of coordinate exchange, a new solution is explored in the areas where the exchange of coordinates may have an improvement, so as to retain the good features and attributes of the existing solution and avoid the shortcomings of random restart. The experimental results and analysis show that the MP-CE method is effective for experimental designs with multiple types of factors in complex constraints spaces, coupled with good performance in terms of space-filling properties.

Conclusions
This paper proposes the MP-CE method to solve the problem of experimental design with multiple types of factors and complex constraints between factors. A coordinate selection strategy based on the projection property and distance property of design points is designed. By assigning probability of coordinate selection to overcome the shortcomings of randomness and easy to fall into local optimum and by adding perturbation to each local search of coordinate exchange, a new solution is explored in the areas where the exchange of coordinates may have an improvement, so as to retain the good features and attributes of the existing solution and avoid the shortcomings of random restart. The experimental results and analysis show that the MP-CE method is effective for experimental designs with multiple types of factors in complex constraints spaces, coupled with good performance in terms of space-filling properties.