## 3 Design Methodology
### Optimization Engine  
The optimization engine used in this work for solving Equations (4), (5), and (7) is a local gradient method [36]. The algorithm is embedded in a trust region (TR) framework [37]. The latter generates a series (\(i=1,2,3,\dots\) ) of local models around what is currently the best design \(\mathbf{x}^{(i)}\) in order to gradually approximate the optimal solution to be found. The TR-based optimization is defined as [37]  
\[\mathbf{x}^{(i+1)}=\arg\min_{\mathbf{x}\parallel\mathbf{x}-\mathbf{x}^{(i)}\parallel\leq r^{(i )}}\digamma\mathbf{F}\Big{(}\mathbf{G}^{(i)}(\mathbf{x})\Big{)} \tag{8}\]  
where \(\mathbf{G}^{(i)}\) is a first-order Taylor expansion model, given by [51]  
\[\mathbf{G}^{(i)}(\mathbf{x})=\mathbf{R}\Big{(}\mathbf{x}^{(i)}\Big{)}+\mathbf{J}\Big{(}\mathbf{x}^{(i) }\Big{)}\cdot\Big{(}\mathbf{x}-\mathbf{x}^{(i)}\Big{)} \tag{9}\]  
where \(\mathbf{J}(\mathbf{x}^{(i)})\) represents the Jacobian generated around \(\mathbf{x}^{(i)}\) using a large-step finite differentiation [51]. The TR radius \(r^{(i)}\) (i.e., the region around \(\mathbf{x}^{(i)}\) for which the model is considered valid) is adjusted based on the calculated gain ratio [51]:  
\[\rho=\frac{\mathbf{R}\Big{(}\mathbf{x}^{(i+1)}\Big{)}-\mathbf{R}\Big{(}\mathbf{x}^{(i)}\Big{)} }{\mathbf{G}^{(i)}\big{(}\mathbf{x}^{(i+1)}\big{)}-\mathbf{G}^{(i)}\big{(}\mathbf{x}^{(i+1)} \big{)}} \tag{10}\]  
The candidate design \(\mathbf{x}^{(i+1)}\) is accepted only if \(\rho>0\). Otherwise, the iteration is reset. The radius is updated according to the standard rules [37]. The initial value \(r^{(0)}\) is controlled by the nested optimization loop (i.e., for \(n=1\), the radius is set to \(r^{(0)}=1\), whereas for \(n>1\), \(r^{(0)}=1/d_{n}\)) (see Section 3.2). The algorithm is terminated either when the Euclidean distance between two consecutive designs or the radius are smaller than a specified threshold (i.e., \(\mid\mid\mathbf{x}^{(i+1)}-\mathbf{x}^{(i)}\mid\mid\leq\varepsilon\), or \(r^{(i+1)}\leq\varepsilon\) (here, \(\varepsilon=10^{-2}\))). The algorithm can be summarized as follows:  
1. Set \(i=0\) and \(\mathbf{x}^{(0)}\), then set \(r^{(0)}\) according to the iteration of the nested loop in Section 3.2;
2. Generate the Jacobian \(\mathbf{I}(\mathbf{x}^{(i)})\) and construct the local approximation model \(\mathbf{G}^{(i)}\);
3. Find \(\mathbf{x}^{(i+1)}\) by solving Equation (8);
4. Evaluate \(\mathbf{R}(\mathbf{x}^{(i+1)})\) and calculate \(\rho\);
5. If \(\rho>0.75\), set \(r^{(i+1)}=2r^{(i)}\). If \(\rho<0.25\), set \(r^{(i+1)}=r^{(i)}/3\);
6. If \(\rho<0\), go to 3; otherwise, go to 7;
7. If the termination condition is met, END; otherwise, set \(i=i+1\) and go to 2.  
Figure 5: Generation of the optimal designs between the existing solutions for \(K=2\), where: (**a**) \(k=1\) (before reindexing), (**b**) \(k=K\) (before reindexing), and (**c**) \(k=K\) (after reindexing). Note that the candidate designs (red squares) are obtained through approximation of the geometries optimized in previous iterations using \(d_{m}\) dimensions. The white and green triangles represent the optimized designs before and after reindexing, respectively. The final designs can be represented using a varying number of variables.  
Owing to utilization of the Taylor model, the computational cost of the algorithm is low, as it requires \(d_{n}+1\) simulations of the \(R\) model per successful iteration. Additional evaluations are required for the steps that do not improve the objective function value. The justification for decreasing the initial radius when \(n>1\) is that with each iteration of the nested loop, the goal of the algorithm is to exploit the region identified in the previous step. For more comprehensive discussion on TR optimization methods, see [37, 51].