A Generalized SDP Multi-Objective Optimization Method for EM-Based Microwave Device Design

In this article, a generalized sequential domain patching (GSDP) method for efficient multi-objective optimization based on electromagnetics (EM) simulation is proposed. The GSDP method allowing fast searching for Pareto fronts for two and three objectives is elaborated in detail in this paper. The GSDP method is compared with the NSGA-II method using multi-objective problems in the DTLZ series, and the results show the GSDP method saved computational cost by more than 85% compared to NSGA-II method. A diversity comparison indicator (DCI) is used to evaluate approximate Pareto fronts. The comparison results show the diversity performance of GSDP is better than that of NSGA-II in most cases. We demonstrate the proposed GSDP method using a practical multi-objective design example of EM-based UWB antenna for IoT applications.


Introduction
Nowadays, full-wave electromagnetic (EM) analysis has become a necessary tool for designing microwave devices. Full-wave EM analysis methods include method of moment (MoM) [1] and the finite-difference time domain method (FDTD) [2], et al. At present, microwave device models can be analyzed by commercial full-wave EM-based simulation software packages such as CST and HFSS. The full-wave EM analysis process is computationally expensive, especially for complex models. When the model is designed using numerical optimization methods such as the quasi-Newton method [3], genetic algorithm [4] or other population-based metaheuristic methods, the cost may become prohibitive.
In the design process of microwave devices, such as multi-port component design, many design objectives need to be considered at the same time, such as return loss [5], insertion loss [6] and so on. In many cases, especially IoT component design, the miniaturization of devices is also an objective that a designer has to consider [7]. Therefore, the design process can be seen as a multi-objective problem that takes many design objectives into consideration. These objectives require simultaneous optimization and are usually in conflict with each other, which means that improving one objective may degrade other objectives [8].
The commercial EM-based simulation software packages (CST [9] and HFSS [10]) allow the inclusion of many optimization objectives. Traditionally, these objectives can be converted into a single objective by weighted sum and then optimized by popular methods such as the genetic algorithm and pattern search method [11]. However, this method can only find one of many possible solutions, and the solution found by the arbitrary setting of the weight of each objective may not be satisfactory. Currently, the most robust multi-objective optimization algorithms with good convergence are population-based metaheuristic multi-objective optimization methods (such as NSGA-II [12]) and multi-objective particle swarm optimization methods [13]. In the case of EM-based component design [14][15][16] (each function evaluation normally requires a full-wave EM analysis), these population-based metaheuristic methods are not applicable because of the large amount of objective function evaluation required and thus the very high computation cost. Multi-objective optimization can also be carried out by constructing surrogate models of EM-based simulation of microwave devices [17][18][19][20] and combining with optimization algorithms such as NSGA-II. The computational cost is saved, but the accuracy of the solutions may suffer.
The sequential domain patching (SDP) method [21] overcomes the problem of a large number of function evaluations (EM simulations) by means of searching small areas (patch) sequentially between two optimal solutions. Since the SDP method uses EM simulation directly for multi-objective optimization, the error is relatively small. However, the SDP method is only applicable to two-objective optimization problems.
In this work, a generalized SDP (GSDP) method is proposed to efficiently solve EM-based multi-objective optimization. The method is applicable to solve optimization problems of two or three objectives, and it has the potential to be expanded to more objectives. Compared with the existing algorithms such as NSGA-II [12], the GSDP greatly reduces the computational cost and ensures a good diversity of optimization results.

Problem Formulation
Let F k (x), k = 1, . . . , N obj , be the kth design objective (x stands for a vector of design parameters). The aim of the multi-objective optimization is to find an approximate Pareto set X so that for any x ∈ X, there is no y for which y ≺ x. Here, ≺ is the dominance relation defined as y ≺ x (y dominates over x) if F k (y) ≤ F k (x) for all k = 1, . . . , N obj , and F k (y) < F k (x) for at least one k [22]. A Pareto set in design space corresponds to a Pareto front [23] in the objective space. All the design solutions in X have the best possible trade-offs between the considered objectives.

Framework of GSDP Algorithm
The GSDP algorithm is designed to solve multi-objective problems of two or three objectives. The algorithm consists of four parts. The first step is to find all the extreme Pareto-optimal design solutions. Secondly, the Pareto set boundaries between every two extreme Pareto-optimal design solutions are found through an updated sequential domain patch (SDP [21]) method. If the problem is a two-objective optimization problem, the obtained boundary is defined as a "design profile". The profile is then refined to obtain the final approximate Pareto set.
If the problem is a three-objective optimization problem, the three Pareto set boundaries form a design profile. The profile is then filled with Pareto-optimal design solutions using the updated SDP method. These Pareto-optimal design solutions are refined to obtain the final approximate Pareto set. The flowchart of the GSDP is shown in Figure 1.

Determining Extreme Pareto-Optimal Designs
The GSDP algorithm begins by obtaining the extreme Pareto-optimal designs [24]. Here, the extreme pareto-optimal designs are a set of optimal solutions for all single objectives. The optimal solution x * k for kth objective, k = 1, . . . , N obj , is defined as follows:

Determining Extreme Pareto-Optimal Designs
The GSDP algorithm begins by obtaining the extreme Pareto-optimal designs [24]. Here, the extreme pareto-optimal designs are a set of optimal solutions for all single objectives. The optimal solution x * k for kth objective, k = 1, …, Nobj, is defined as follows:

Determining Design Profile
After all the extreme Pareto-optimal designs are obtained, we use an updated SDP method to obtain the Pareto set boundaries between every two extreme Pareto-optimal designs. There is only one Pareto set boundary which can be obtained when it is a two-objective problem, and this boundary forms a so-called design profile. When it is a three-objective problem, all three Pareto set boundaries are linked together to form a design profile. This process of finding a Pareto set boundary is equivalent to solving a two-objective optimization problem.
Suppose a two-objective optimization problem is defined as follows: where x is the n-dimensional design parameter vector. * 1 x and * 2 x are the single objective optimal solutions of F1 and F2, respectively,

Determining Design Profile
After all the extreme Pareto-optimal designs are obtained, we use an updated SDP method to obtain the Pareto set boundaries between every two extreme Pareto-optimal designs. There is only one Pareto set boundary which can be obtained when it is a two-objective problem, and this boundary forms a so-called design profile. When it is a three-objective problem, all three Pareto set boundaries are linked together to form a design profile. This process of finding a Pareto set boundary is equivalent to solving a two-objective optimization problem.
Suppose a two-objective optimization problem is defined as follows: where x is the n-dimensional design parameter vector. x * 1 and x * 2 are the single objective optimal solutions of F 1 and F 2 , respectively, The method to find the Pareto set boundary of the two objectives is described in detail as follows: 1. Determine the design space; the lower bound vector lb and upper bound vector ub can be defined as follows: where lb j and ub j are the jth variables of vectors lb and ub respectively; x * 1 j and x * 2 j are the jth variables of vectors x * 1 and x * 2 respectively. only) and select the one that leads to the largest improvement with regard to the first objective F 1 as the current point x 2 c . 6. If x 1 c or x 2 c exceeds the design space, the process is stopped. If both x 1 c and x 2 c are in the design space, return to step 4.
In step 2, the perturbation size vector d can be obtained by the following formula: where d j ( j = 1, . . . , n) means the size of the perturbations in the jth direction. The value of K is a user-defined constant (typically between 5 and 20).
As an illustration, the above method is applied to solve a two-objective optimization problem: where the design parameter dimension n = 2 and constant K = 15; the process to find the Pareto set boundary is shown in Figure 2. Corresponding to the Pareto set boundary for the two-objective problem in Figure 2, the Pareto front boundary in the objective space is shown in Figure 3.    Since (8) is a two-objective optimization problem, the design profile is formed by the Pareto set boundary (shown as red dots in Figure 2b). Its corresponding Pareto front boundary (shown in Figure  3) is called the "objective profile". According the GSDP algorithm flowchart (shown in Figure 1), all the solutions in the design profile need to be refined to obtain the final approximate Pareto set. Its corresponding approximate Pareto front is shown in Figure 4. The refinement process is described later in Section 2.6. Since (8) is a two-objective optimization problem, the design profile is formed by the Pareto set boundary (shown as red dots in Figure 2b). Its corresponding Pareto front boundary (shown in Figure 3) is called the "objective profile". According the GSDP algorithm flowchart (shown in Figure 1), all the solutions in the design profile need to be refined to obtain the final approximate Pareto set. Its corresponding approximate Pareto front is shown in Figure 4. The refinement process is described later in Section 2.6.  If the problem is a three-objective optimization problem, the extreme Pareto-optimal designs X12, X13, and X23 are obtained first, where Xij is the Pareto set boundary of the ith and jth objectives found by updated SDP, and Fij is the Pareto front boundary corresponding to Xij. F12, F13 and F23 form the objective profile for three objectives. The conceptual results for the obtained design profile and objective profile for three objectives are illustrated in Figure 5. If the problem is a three-objective optimization problem, the extreme Pareto-optimal designs X 12 , X 13 , and X 23 are obtained first, where X ij is the Pareto set boundary of the ith and jth objectives found by updated SDP, and F ij is the Pareto front boundary corresponding to X ij . F 12 , F 13 and F 23 form the objective profile for three objectives. The conceptual results for the obtained design profile and objective profile for three objectives are illustrated in Figure 5.

Filling Design Profile
In the case of three design objectives ( 1 f , 2 f and 3 f ), the corresponding extreme optimal designs are * 1 x , * 2 x and * 3 x , respectively. The design profile found following Section 2.4 is filled by the process described below.
The process finds the appropriate design set

Filling Design Profile
In the case of three design objectives ( f 1 , f 2 and f 3 ), the corresponding extreme optimal designs are x * 1 , x * 2 and x * 3 , respectively. The design profile found following Section 2.4 is filled by the process described below.
The process finds the appropriate design set X i between x * i and its opposite Pareto set boundary X jk (j, k = 1, 2, 3 and j, k i) for all i = 1, 2, 3. The design sets X i (i = 1, 2, 3) and the design profile are combined to form the initial approximate Pareto set S. In other words, the design profile is filled by the three similar steps shown in the flowchart in Figure 6.
Here, we only describe the process of obtaining the solution set X 1 between x * 1 and X 23 for ease of reading. The processes to find X 2 and X 3 are similar.

1.
Let X i 23 be the ith solution in X 23 , i = 1.

2.
If i > N, the process stops, where N is the number of solutions in X 23 .

3.
Determine the design subspace, the lower bound vector lb and upper bound vector ub can be defined as follows: where lb j and ub j are the jth variables of vectors lb and ub, respectively, and x * 1 j and X i The weighting factors a and b in step 6 are defined as follows: where dl is the distance between X i 23 and x * 2 in the design space, and dr is the distance between X i 23 and x * 3 . When x * 1 and X i 23 are taken as starting points to locate set X 1 by the process mentioned above, the relative location of X i 23 in X 23 affects the selection of the perturbed design. As shown in Figure 7a, if dl < dr, the selection is biased towards x * 3 and if dl > dr, the selection is biased towards x * 2 . Therefore, the weighting factors a and b are necessary to balance the bias. The process to connect x * 1 and X i 23 is illustrated as Figure 7.
When the filling design profile is complete, all the obtained paths and design profiles obtained in Section 2.4 together constitute the initial approximate Pareto set. Here, we only describe the process of obtaining the solution set 1 X between * 1 x and 23 X for ease of reading. The processes to find 2 X and 3 X are similar. 3. Determine the design subspace, the lower bound vector lb and upper bound vector ub can be defined as follows: Obtain the solution set between and Obtain the solution set between and Obtain the solution set between and Obtain the initial approximate Pareto set

Pareto Set Refinement
In the above method, the initial approximate Pareto set and the corresponding Pareto front are obtained, so there are still many dominated solutions in the result. These dominated solutions are deleted to ensure the results are all non-dominated solutions (the final approximation Pareto set). 23 3 1 23 1 X mentioned above, the relative location of 23 i X in 23 X affects the selection of the perturbed design.
As shown in Figure 7a, if dl < dr, the selection is biased towards * 3 x and if dl > dr, the selection is biased towards * 2 x . Therefore, the weighting factors a and b are necessary to balance the bias. The process to connect * 1 x and 23 i X is illustrated as Figure 7. When the filling design profile is complete, all the obtained paths and design profiles obtained in Section 2.4 together constitute the initial approximate Pareto set.

Evaluation
The Diversity Comparison Indicator (DCI) [25] is used to assess the diversity of the obtained Pareto front with respect to the other Pareto fronts. In this method, all the concerned Pareto fronts are compared together and the dominated solutions are removed. The remaining solutions are put into a grid environment. DCI considers the total contribution of each Pareto front to all the nonempty grid cells. The DCI method can only be applied to the assessment of relative diversity of a Pareto front against others rather than providing an absolute measure of distribution for a single Pareto front. The DCI value is a number in the range [0, 1], and the higher the DCI value, the better the distribution of this Pareto front.

Illustration and Verification Examples
In this section, we use 4 DTLZ three-objective functions to verify the GSDP algorithm.

Design Case 1: DTLZ1
The DTLZ1 [26] function is as follows: where x i means the ith variable of n-dimensional vector parameter x.
Here, parameter x has a dimension of n = 6. The built-in function in MATLAB fmincon is used with the initial optimization points [0, 0, 0, 1, 1, 1] to obtain the three single-objective optimal solutions, x * 1 = [0.0008, 0.0008, 0.5000, 0.5000, 0.5000, 0 The distribution performances of the obtained two Pareto fronts are evaluated by the Diversity Comparison Indicator (DCI) [25]. In this case, the DCI value of P1 is 0.7667, while the DCI value of P2 is 0.2333. The distribution performances of the obtained two Pareto fronts are evaluated by the Diversity Comparison Indicator (DCI) [25]. In this case, the DCI value of P1 is 0.7667, while the DCI value of P2 is 0.2333.
As demonstrated in this example, GSDP is much more efficient than NSGA-II. Moreover, the distribution of the approximate Pareto front of GSDP is relatively better than that of NSGA-II.

Design Case 2: DTLZ2
The DTLZ2 [26] function is as follows: where x i means the ith variable of the n-dimensional vector parameter x. At n = 9, the built-in function in MATLAB fmincon is used with the initial optimization points [0, 0, 0, 1, 1, 0, 0, 1, 0] to obtain the three single-objective optimal solutions x * As a comparison, NSGA-II is used to obtain 833 solutions of the approximate Pareto front. Each objective function is evaluated 268,940 times. GSDP saves about 97% of computational cost compared to NSGA-II. The approximate Pareto fronts P1 and P2 obtained by the GSDP and NSGA-II, respectively, are shown in Figure 9.
After combining P1 and P2 and deleting the dominated solutions, all 833 solutions in the P1 are retained, while only 454 solutions in the P2 are retained. The DCI values of P1 and P2 are 0.4840 and 0.5426, respectively. Although the DCI of P1 is slightly less than that of P2, the GSDP has absolute advantages in terms of computational cost.

Design Case Validation Summary
In addition to the above test results, the GSDP method is tested using more three-objective functions with variables of various dimensions. We summarize the results in Table 1. P1 and P2 represent the approximate Pareto fronts obtained by the GSDP and NSGA-II, respectively. P1 and P2 have the same number of solutions. It can be seen that for a variety of problems and variable dimensions, the GSDP method performs better than or similar to NSGA-II in terms of distribution, but GSDP saves more than 85% in computational cost for all the testing cases. After combining P1 and P2 and deleting the dominated solutions, all 833 solutions in the P1 are retained, while only 454 solutions in the P2 are retained. The DCI values of P1 and P2 are 0.4840 and 0.5426, respectively. Although the DCI of P1 is slightly less than that of P2, the GSDP has absolute advantages in terms of computational cost.

UWB Antenna Multi-Objective Design Example
The proposed GSDP method is applied to the multi-objective optimization of a microwave device: a UWB antenna [27]. The antenna is implemented on an FR4 substrate (ε = 4.3, h = 1.55 mm, tanδ = 0.02). Here, f 1 (the first objective) is used to minimize |S 11 |(dB) between 4 GHz and 10 GHz, while f 2 (the second objective) is used to minimize the difference of the realized gain between 4 GHz and 10 GHz and f 3 (the third objective) is used to minimize the UWB antenna size. The antenna structure and design parameter vector x = [a, b, d, kL, ds, dW, dWs, L 0 , Lg, Ls, Ws] are shown in Figure 10.

UWB Antenna Multi-Objective Design Example
The proposed GSDP method is applied to the multi-objective optimization of a microwave device: a UWB antenna [27]. The antenna is implemented on an FR4 substrate (ε = 4.3, h = 1.55 mm, tanδ = 0.02). Here, 1 f (the first objective) is used to minimize 11  Here, 0 0.731 W = is fixed. The first step of the GSDP method is to find the three extreme optimal design solutions using single-objective optimization (here, a genetic algorithm is used): * 1  Figure 11. Here, W 0 = 0.731 is fixed. The first step of the GSDP method is to find the three extreme optimal design solutions using single-objective optimization (here, a genetic algorithm is used):  Figure 11. Some of the solutions in the Pareto front shown in Figure 11 may not applicable to the design of UWB antennae. For example, the 1 f ( 11 S ) value is greater than -7 dB. These solutions are discarded.
Then, the final UWB antenna multi-objective design results are shown in Figure 12. Eight solutions ( (1) x , (2) x ,…, (8) x ) are selected according to customized specifications (e.g., all are below −9 dB between 4 GHz and 10 GHz, all the realized gain (dBi) variation are below 2.6 between 4 GHz and 10 GHz, and sizes are below 400 mm 2 ). Their objective values and parameter values are shown in Table  2. The 11 (dB) S responses (reflection coefficient) and the realized gain (dBi) are shown in Figure 13. A total of 3466 CST simulations are used to optimize the UWB antenna, which includes 1600 simulations to obtain the extreme optimal solutions (single-objective optimization process). This takes about 48.14 h on an Intel(R) Core i5-6500 processor with 8 GB RAM. If the NSGA-II method were used for the three-objective optimization of this UWB antenna (11 design variables), the computational cost is estimated at more than 100,000 CST EM simulations (see DTLZ2 (11) in Table  1), which is prohibitive. Some of the solutions in the Pareto front shown in Figure 11 may not applicable to the design of UWB antennae. For example, the f 1 (|S 11 |) value is greater than −7 dB. These solutions are discarded. Then, the final UWB antenna multi-objective design results are shown in Figure 12. Eight solutions (x (1) , x (2) , . . . , x (8) ) are selected according to customized specifications (e.g., all are below −9 dB between 4 GHz and 10 GHz, all the realized gain (dBi) variation are below 2.6 between 4 GHz and 10 GHz, and sizes are below 400 mm 2 ). Their objective values and parameter values are shown in Table 2. The S 11 (dB) responses (reflection coefficient) and the realized gain (dBi) are shown in Figure 13.
A total of 3466 CST simulations are used to optimize the UWB antenna, which includes 1600 simulations to obtain the extreme optimal solutions (single-objective optimization process). This takes about 48.14 h on an Intel(R) Core i5-6500 processor with 8 GB RAM. If the NSGA-II method were used for the three-objective optimization of this UWB antenna (11 design variables), the computational cost is estimated at more than 100,000 CST EM simulations (see DTLZ2 (11) in Table  1), which is prohibitive.    (1) x (2) x (3) x (4) x (5) x (6) x (7) x (8) Figure 13. The selected design results ( (1) x , ( 2) x ,…, (8) x ref. Table 2) on the approximate Pareto set obtained by GSDP: (a) All 11 (dB) S responses are below -9 dB between 4 GHz and 10 GHz; (b) All the realized gain (dBi) variationd are below 2.6 between 4 GHz and 10 GHz. Figure 13. The selected design results (x (1) , x (2) , . . . , x (8) ref. Table 2) on the approximate Pareto set obtained by GSDP: (a) All S 11 (dB) responses are below −9 dB between 4 GHz and 10 GHz; (b) All the realized gain (dBi) variationd are below 2.6 between 4 GHz and 10 GHz. Table 2. Objective values and parameter values of selected designs found by GSDP.
x (1) x (2) x (3) x (4) x (5) x (6) x (7) x (8)  A total of 3466 CST simulations are used to optimize the UWB antenna, which includes 1600 simulations to obtain the extreme optimal solutions (single-objective optimization process). This takes about 48.14 h on an Intel(R) Core i5-6500 processor with 8 GB RAM. If the NSGA-II method were used for the three-objective optimization of this UWB antenna (11 design variables), the computational cost is estimated at more than 100,000 CST EM simulations (see DTLZ2 (11) in Table 1), which is prohibitive.

Conclusions
An efficient multi-objective optimization method (GSDP) for EM-based microwave device design is presented. The GSDP method generalizes the SDP method to solve multi-objective problems of both two and three objectives. The GSDP method includes four parts: (i) determining the extreme Pareto-optimal designs; (ii) determining the design profile; (iii) filling the design profile; and (iv) refining the Pareto set. Each part is described in detail. The GSDP method is compared to NSGA-II in terms of computational efficiency and the performances of Pareto fronts. The GSDP method saves more than 85% in computational cost compared to NSGA-II with similar or better performance in terms of Pareto front distribution (DCI). A UWB antenna multi-objective EM-based design example demonstrates the efficiency of GSDP in finding the approximate Pareto set (front). Selected designs from the approximate Pareto set satisfy all three customized specifications.
Author Contributions: Y.L., Q.S.C. and S.K. devised the research and wrote the paper; Y.L. and Q.S.C. designed algorithms and the experiments; Y.L. and Q.S.C. improved the English. All authors have read and approved the final manuscript.