A Surrogate-Assisted Multiconcept Optimization Framework for Real-World Engineering Design Open Access

Design problems often have multiple conceptual solutions, referred to as concepts, typically represented using different variables. Ideally, designers would optimize across these concepts to identify the best-performing concept(s) and their corresponding design(s). However, existing optimization methods operate with a fixed set of variables, restricting their use in concept design. Multiconcept optimization (MCO) algorithms bridge this gap by allowing searches across multiple concept spaces. Yet, two key elements in MCO require further development for practical use: (1) efficient use of approximations to guide the search and (2) handling analysis failure during performance assessment. To address these challenges, we introduce an MCO framework that solves unconstrained and constrained single-objective optimization problems with a limited computing budget. The framework incorporates four surrogate-assisted optimization algorithms: predictor believer (PB), infeasibility preserved believer (IPB), and two novel approaches—enhanced constrained expected improvement (EcEI) and Bradley–Terry-based probabilistic sorting (BTPS). All these algorithms can handle analysis failure. For maximum flexibility in functional representation, the algorithms dynamically select surrogates from 23 classifiers and 41 regressors during the search. We demonstrate the framework on various analytical and practical examples, including single-concept constrained optimization problems (G-series), a modified G24 problem with analysis failure, a beam design problem involving six concepts, and a 3D shape-matching problem for coronary stent designs involving four concepts. Furthermore, we present the first parameterization scheme to represent single-helix stent designs. We believe that our contribution will enhance the adoption of optimization methods in concept design.

The engineering product design process comprises four phases—problem definition, concept design, preliminary design, and detailed design [1]. Concept design is a crucial phase where multiple concepts are assessed, and the most promising one is carried forward to the preliminary design phase for further development. It is estimated that the decisions in early design phases (concept and preliminary design) influence up to 70% of the product life-cycle cost [2]. While optimization methods are regularly used in preliminary design to improve performance, their use in the concept design remains limited. This limitation arises because existing optimization algorithms, in their native form, cannot simultaneously search across multiple concepts defined by different numbers and types of variables. In recent years, there has been significant interest in developing optimization methods capable of simultaneously searching across multiple concept spaces, and such methods are referred to as multiconcept optimization (MCO) algorithms.

MCO algorithms have been proposed to address both unconstrained and constrained single- [3,4] and multiobjective optimization problems [5–10]. For single-objective optimization, an unconstrained MCO algorithm was developed by modifying simulated annealing to probabilistically select a concept in each iteration based on the concepts’ past performance [4]. Multiobjective MCO algorithms aim to identify the overall Pareto front, referred to as s-Pareto (“set-based” Pareto) [5] or C-Pareto (“concept-based” Pareto) [6] solutions. Although the overall Pareto front/set can be generated by independently optimizing each concept and combining the results, this approach is computationally inefficient. Therefore, multiobjective MCO algorithms tend to concurrently evolve all participating concepts to yield the overall Pareto front [7,9,10]. In these approaches, traditional multiobjective evolutionary algorithms (MOEAs) such as NSGA-II [7], $ϵ$-MOEA [9], and differential evolution [10] have been adapted with modifications. The modifications include (1) concept-wise recombination, where parents from the same concept are only allowed to recombine to create offspring, and (2) modified environmental selection, where all solutions from all concepts are considered together and ranked. The environmental selection scheme was based on nondominated sorting of all solutions from all concepts, followed by concept-based crowding sort [7]. Similarly, Moshaiov and Snir [9] modified $ϵ$-MOEA by considering all solutions across concepts to decide if a candidate can be added to the population and/or the archive. To handle constraints, MCO algorithms have used the constraint domination principle [10]. Most studies used simple analytical benchmarks to illustrate the performance of MCO algorithms, with only a few reports utilizing practical examples such as truss design [5], lattice structure design [10], and robot path planning [8]. To facilitate further development of MCO algorithms, a test suite of analytical benchmark MCO problems with varying characteristics was recently introduced [11].

Despite significant advances, the application of current MCO algorithms on real-world optimization problems remains limited. In such settings, objectives and/or constraints are often evaluated via computationally expensive physics-based simulations or physical experiments, which severely restricts the number of allowable function evaluations. Most existing MCO techniques require a large number of evaluations, e.g., 22,964 and 73,000 solutions were evaluated for clustering a warehouse problem [4], 25,000 evaluations were required for a simple two-concept MCO [9], and a modest limit of 1548 evaluations was used for lattice structure optimization [10]. Surrogate-assisted optimization (SAO) is an effective approach for problems with limited computational budget. However, only a few SAO algorithms have been proposed for single- and multiobjective MCO problems. For single-objective MCO, three SAO strategies were recently introduced following a steady-state paradigm [3]: a believer, a Bayesian optimization (BO), and a tree-structure-based optimization. The believer strategy individually ran NSGA-II on surrogate models (Gaussian process or random forest). The best-predicted solution after SAO of all concepts was chosen for true evaluation and subsequently added to the surrogate training set. The BO strategy maximized an expected improvement (EI) acquisition function, thereby balancing the exploitation of promising regions and exploration of uncertain areas. The tree-structure-based optimization strategy used a multivariate tree-structured Parzons estimator (TPE) to predict a global rank of sampled solutions across concepts, selecting the highest ranked solution for evaluation. A comparison of these strategies on different MCO problems showed that the believer model often converged to local minima, BO performed well across various problems, and TPE showed the best exploration of the design space. For the multiobjective multiconcept SAO study [10], two surrogate-assisted MCO strategies were introduced alongside three direct evolutionary MCO methods. A lattice design problem was used as a case study with Kriging-based surrogates. The performance of the surrogate-assisted forms was suboptimal, largely due to surrogate inaccuracies at the steep sections of the Pareto front [10]. Overall, the current SAO algorithms for MCO have used only one or at most two predefined surrogate model types during optimization. Given the wide range of regressors and classifiers available in the modern machine learning domain, there is an opportunity to offer greater flexibility in function representation. A further complication in most real-life scenarios is that the simulations or experiments can fail to produce a numerical outcome for certain candidate designs. This issue of analysis failure in MCO or, more broadly, in SAO is a major challenge that requires further attention. Moreover, the effect of analysis failure on the performance of an optimization algorithm is amplified in settings with limited function evaluation budgets. Below, we provide a brief overview of standalone SAO algorithms (discussed in detail in Refs. [12–14]) and review in depth the existing approaches for dealing with analysis failure.

SAO algorithms are widely used for handling computationally expensive optimization problems with limited computing budgets. During SAO, approximators are typically used to predict the performance of the solutions as opposed to costly true evaluations. Thereafter, one or more promising solutions are selected for true evaluation, and the approximators are updated accordingly. In believer-based SAO approaches, only the predicted performance of a solution is used to identify the candidate solution for true evaluation. However, BO utilizes both the predicted performance and the uncertainty in predictions through an acquisition function to select candidate solutions for true evaluation. While believer-based methods have long been used due to their simplicity and ease of integration with a range of optimization algorithms, BO has seen increasing applications in recent years. BO Variants with different acquisition functions such as EI [15], knowledge gradient [16], entropy search [17], upper constraint bound [18], sampling-schemes like Thompson sampling [19], and hybrid acquisition function-sampling strategies such as predictive entropy search [20] have also been suggested in the literature. Originally developed for unconstrained problems, these schemes have been modified to tackle constrained optimization problems, resulting in approaches such as constrained EI (cEI) [21,22], constrained knowledge gradient [23], sampling-based methods like scalable constrained Bayesian optimization [24], and hybrid methods such as predictive entropy search with constraints [25]. In acquisition function-based methods, the objective function and the constraints are independently modeled, and the acquisition function balances exploitation and exploration. EI still remains the most popular form for BO [14], and its constrained counterpart, cEI, modifies the EI acquisition function by multiplying it by the probability of feasibility (PF). However, since the product of the probability of each constraint is used to compute the feasibility of a solution, the quantity quickly approaches zero, reducing its utility in ordering solutions, especially for problems with a large number of constraints. More elaborate lexicographic ordering schemes could potentially provide better guidance during evolution.

As highlighted above, optimization approaches need to deal with analysis failure, i.e., have mechanisms in place to deal with solutions that do not return a response/output upon evaluation. In literature, this phenomenon has been referred to using various terminologies such as infeasible solutions [26], simulator infeasible solutions [27], invalid vectors [28], hidden constraints [29–31], simulation failure [32–38], experimental failure [39], crash constraints [40,41], noncomputable design points [42], and unknown constraints [43]. In this work, we refer to all of these conditions as “analysis failure.” The specific cause depends on the domain or the analysis type. For instance, physics-based simulations often fail due to nonviable geometries/geometry interference/singularities in CAD constructions [26,35,44], poor quality of meshes, weak solver convergence [35,37,38,42,43], and hard simulation time limits [33]. In physical experiments, failure may arise due to the inability to perform desired objectives [41] or generate target material properties [39]. Effective management of such failures is crucial for maintaining an efficient optimization process. Two fundamental approaches have been reported to address analysis failure. The first relies on imputation, wherein failed solutions are assigned the worst performance value (also called floor padding) [39], or a penalty prescribed by a domain expert [32,39,42]. These penalties ensure that failed solutions are least preferred during the search. However, if surrogates are used to predict performance, artificially imputed values can distort the performance landscape, especially in the vicinity of the analysis failure region, thereby posing significant challenges to the underlying search. The second approach relies on using classifiers to predict the probability of analysis failure. Gaussian process (GP)-based classifiers have been used to differentiate between analysis success and failure solutions [33–35,37–40]). Support vector machine (SVM)-based classifiers have also been extensively used. These include SVM [42,43], canonical SVM [30], probabilistic SVM [45], and single class support vector domain descriptors [46]. Apart from GP- and SVM-based classifiers, there are also reports on the use of k-nearest neighbor [29,43], random forest [43], AdaBoost [43], decision trees [29], logistic regression [29], and ensemble of classifiers [27]. Typically, regressors and classifiers are retrained after evaluation of new solutions. The regressors are only trained using solutions that resulted in analysis success, whereas the classifiers are trained with solutions belonging to two classes (success and failure). In a recent work [33], the classifier training was only undertaken for a fixed number of iterations in the initial phase and then skipped later under the assumption that the classifier is successfully trained. Depending upon the application, significant variation can be observed in the proportion of solutions with analysis failure and thus the class imbalance needs to be managed carefully. Imputation-based schemes have been used in situations where the computational budget is extremely low [39,42], especially in the case of optimization with physical experiments in loop [39]. Apart from assessing such solutions via imputation or classifiers, the underlying SAO approach also plays an important role in the final outcome of the search. Since our focus is on settings with low computational budget, we limit our discussion to SAO algorithms where one or more solutions are likely to result in analysis failure. Both believer and BO approaches have been used in the SAO algorithms to manage analysis failure. In believer-based approaches, classifier prediction is often incorporated as a constraint [27,28]. In BO, various forms of acquisition functions such as EI [29,30,33,38–42], probability of improvement (PI) [30], lower confidence bound [30,45], and upper confidence bound [43] have been used. Such acquisition functions have been multiplied by the probability of analysis success derived from the classifier in Refs. [33,39,40] to deal with analysis failure. There are reports on subtle variations in the acquisition function, such as augmented EI [42] and variable fidelity EI [38]. Apart from directly altering the EI function, lexicographic sorting has also been altered where solutions are first sorted based on predicted analysis success followed by EI, lower constraint bound, or PI metrics within a differential evolution (DE) algorithm [30]. Overall, it is clear that the SAO algorithms for practical use need mechanisms to deal with analysis failure. Since believer and EI-based methods have their own set of advantages, a practically efficient SAO approach should ideally exploit the benefits offered by both. Apart from conventional GP-based regressors and classifiers used in believer or EI-based methods, there is a significant opportunity to exploit the whole suite of regressors/classifiers that have emerged from the machine learning domain.

The MCO may appear similar to optimization using categorical variables, such as the mixed-variable BO. However, there is a key distinction between MCO and optimization with categorical variables. In MCO, the design variables (their numbers, types, and bounds) for each concept typically have their own distinct variable spaces; while for mixed-variable optimization problems involving categorical variables, the search space for the remaining variables (i.e., real and discrete variables) typically remains unchanged for a given categorical variable. For mixed-variable optimization problems with categorical variables, BO modified for categorical spaces is a suitable approach that considers uncertainty in prediction along with mean prediction. For example, Zhang et al. [47] recently introduced an approach to perform BO for problems involving a combination of categorical and continuous variables via the use of GP model with Con Cat kernel. However, there are still a number of improvements that can still be incorporated to better deal with real-world problems. These include (1) the ability to deal with analysis failure, (2) the incorporation of multiple types of regressors/classifiers with dynamic selection among them to offer greater flexibility in functional representation, (3) the use of multiple surrogate-assisted search strategies by capitalizing on the strengths of different formulations as opposed to BO driven by a single acquisition function, and finally, (4) offering a provision of independent definition of continuous variables for each concept with specific variable bounds as opposed to modeling it as a mixed-variable optimization problem. Such an approach would be suitable in settings such as MCO where each concept is independently defined using its own set of variables with different meanings and bounds and demands a computationally efficient approach to operate within the realm of limited function evaluation budgets.

In this article, we present a framework for efficient optimization across multiple concepts with a limited computational budget. The key elements of the framework are outlined below:

1. The approach relies on a steady-state model, i.e., a single solution is evaluated in each iteration. The framework is embedded with four different SAO strategies relying on a predictor believer (PB), an invisibility preserved believer (IPB), a novel enhanced constrained expected improvement (EcEI), and a novel Bradley–Terry-based probabilistic sorting (BTPS) scheme. Analysis failure is managed using a classifier, while the objective and constraints are modeled using independent regressors. To offer maximum flexibility in representation, $23$ different classifiers and $41$ regressors have been used in the framework.

2. To deal with multiple constraints and analysis failure, a novel EcEI optimization strategy has been introduced. EcEI adopts a lexicographic ordering based on the traditional cEI, probability of analysis success, PF, constraint violation, and objective function value to order solutions which in turn guides the environment selection step of the evolutionary search process.

3. A novel BTPS-SAO strategy relying on probabilistic sorting based on Bradley–Terry model is included in the framework. The same BTPS scheme is also used to select one solution for true evaluation from multiple candidate solutions identified via PB, IPB, EcEI, and BTPS strategies.

In addition to the optimization framework, we present a shape parameterization scheme for single-helix coronary stent designs, which has never been attempted before [48]. The geometry is defined using five design variables via a combination of PARSEC [49], ellipse, and straight line segments. We demonstrate the principle of MCO using a 3D shape-matching problem where a search is conducted across four realistic coronary stent concepts.

The article is organized in five sections. Section 1 provided the introduction and background of MCO along with relevant literature pertaining to SAO and methods to handle analysis failure. Section 2 describes the proposed MCO framework, highlighting the novelties in the current work, including the consideration of analysis failure and selection of best surrogate models for SAO, EcEI, and BTPS sorting algorithms. The performance of the approach is first assessed using 13 single-concept, single-objective constrained optimization problems from the G-series benchmark [50] and compared with state-of-the-art GPSAF algorithm [51]. We then present a modified G24 problem involving analysis failure and illustrate the working principle of the approach. Thereafter, a multiconcept beam design optimization problem is presented in Sec. 3. Finally, in Sec. 4, we present a shape parameterization scheme for single-helix coronary stent design and conduct multiconcept optimization using a 3D shape-matching problem where the search is conducted across four realistic coronary stent design concepts. We conclude the findings and discuss the limitations of the current work in Sec. 5.

In this section, we present a framework for multiconcept design optimization using a steady-state paradigm. The overall architecture of the framework is presented in Fig. 1, consisting of four major steps: initialization and evaluation of solutions, construction of surrogates, generation of candidates via optimization, and candidate selection for evaluation. The first step is performed only once at the start of the optimization process, whereas subsequent steps are repeated iteratively until the total evaluation budget is exhausted.

The first step involves initialization of solutions for each concept. The framework offers sampling based on Latin hypercube sampling (LHS) and Random Sampling. All generated solutions are rounded to five decimal places, and the integer solutions are rounded to the nearest integer prior to evaluation. In this study, we focus on single-objective, unconstrained and constrained optimization problems.

The initialized solutions for each concept are evaluated, i.e., the objective function and the constraint functions (if they exist in the case of constrained optimization) are computed. We assume that numerical techniques such as finite element method (FEM) or computational fluid dynamics (CFD) are used to evaluate the objectives and constraints of the problem, where evaluations may fail to produce a result (i.e., an analysis failure). The framework is designed to handle such situations and a not a number (NaN) value is recorded as an output against the performance of the solution for the objective/constraint functions. For a concept defined by $D$ variables, the regressors downstream in the framework require at least $⌈D/0.8⌉$ solutions with valid objective/constraint values (assuming an 80–20 train–test split). If the number of valid initial solutions is less than $⌈D/0.8⌉$⁠, the initialized solutions for the concept are discarded, and the initialization step is retried with a new random seed. This process is repeated up to 10 times, and if no case produces the required number of valid solutions, the optimization process is terminated. After the solutions of each concept are successfully initialized and evaluated, they are stored in the archive. Thereafter, whenever a candidate solution undergoes true evaluation, the information is added to the archive.

The iterative part of the optimization process starts from the surrogate generation phase. Independent surrogates are created for every objective/constraint for each concept. Two types of surrogates are constructed for each concept, i.e., a classifier that partitions the solution space into regions of analysis success/failure and a regressor for each objective/constraint function. The surrogate generation process is presented in Fig. 2 and the details are outlined below.

To create surrogates, data is retrieved from the archive for each of the concepts. To build the classifier, the data is labeled into two classes: analysis success (label 0) and analysis failure (label 1). The design variables for the concept are used as input features of the classifier. The features are scaled using z-score normalization, and an 80–20 train–test split (with stratification) from scikit-learn library [52] is used to partition the scaled dataset into training and test datasets. All $40$ classifiers within the scikit-learn library were considered as options. However, seven ensemble classifiers (ClassifierChain, MultiOutputClassifier, OneVsRestClassifier, OneVsOneClassifier, OutputCodeClassifier, VotingClassifier, and StackingClassifier) require a base classifier to be specified and, therefore, do not successfully fit the train data with default hyperparameters. Additionally, seven classifiers (LinearSVC, NearestCentroid, Perceptron, RidgeClassifier, RidgeCV, SGDClassifier, and PassiveAggressiveClassifier), which do not provide probability estimation of the classes have been excluded. Finally, three naive Bayes classifiers which require either categorical or nonnegative data (CategoricalNB, ComplementNB, and MultinomialNB) have also been excluded. Thus, the remaining 23 classifiers have been used in the framework. Each of the classifiers is first trained with default hyperparameters and followed up with hyperparameter tuning using 10 RandomizedSearchCV attempts. The hyperparameter range for each tuned classifier is provided in Supplemental Material A available in the Supplemental Materials on the ASME Digital Collection. All classifiers with default and tuned hyperparameters are assessed using the accuracy score on the test data, and the classifier with the highest accuracy score is selected. In the case of a tie between multiple classifiers having the same highest accuracy score, the choice is based on their names following an alphabetical order. Classifiers are only built for a concept if there are at least two solutions with analysis failure in the archive. If there are less than two solutions with analysis failures, a custom classifier is used that suggests 100% probability of analysis success for any given input. In the case where the total number of solutions in the archive is less than 10, a 50–50 train–test split is used so that both classes are represented in the train and test datasets.

While classifiers are used to predict the analysis success/failures of potential solutions, regressors are used to predict the objective and constraint functions. For each concept, independent regressors are trained for the objective and each of the constraint functions. Regressor construction for a function (objective or any constraint) begins by retrieving all data of analysis success after z-score normalization. The scaled data are split into 80–20 for generating training and test datasets. All 55 scikit-learn regressors have been considered for inclusion. However, 14 of them cannot be used with their default hyperparameters and thus have been excluded from the framework. These include four ensemble regressors (VotingRegressor, StackingRegressor, MultiOutputRegressor, and RegressorChain) that require a base classifier to be specified, two regressors (GammaRegressor and PoissonRegressor) that require only positive features (our z-score normalized inputs contain both positive and negative values), four others (MultiTaskElasticNet, MultiTaskElasticNetCV, MultiTaskLasso, and MultiTaskLassoCV) that are limited to the cases with multiple targets, and three more (CCA, PLSCanonical, and QuantileRegressor) that had conflicts with scikit-learn default hyperparameters and finally the IsotonicRegression model that expects only one input. We have also excluded RadiusNeighborsRegressor from the list due to the excessive time required for predictions. The remaining 40 scikit-learn regression models are first trained using default hyperparameters which are followed up with hyperparameter tuning using 10 RandomizedSearchCV attempts. The hyperparameter range for each tuned regressor is provided in Supplemental Material A.

In addition to the above, a python clone of design and analysis of computer experiments (DACE) regressor [53,54] with Gaussian correlation function, linear regression model, and optimized theta parameter was also added to the list of regressors. The performance of the regressors (scikit-learn default hyperparameter regressors, tuned regressors, and optimized DACE model) is based on the mean squared error on the test data. The model with the least mean squared error was selected as the regressor for the objective and each constraint function. During the prediction phase, the algorithm uses bootstrapping to estimate standard deviation if the selected regressor does not inherently provide the error estimates. For such regressors, bootstrapping was used with 100 models as described in our earlier work [26].

In the framework, candidate solutions are generated using different SAO strategies for each concept. The overall process is outlined in Fig. 3. Four SAO strategies are deployed to generate a total of six candidate solutions for each concept. The four SAO strategies include PB, IPB, EcEI, and BTPS. The details of each of these strategies are outlined in the following subsections.

PB refers to a scheme that assumes the predictors are accurate and the optimization algorithm believes their predictions. The optimization process starts with 10D randomly initialized solutions. The trained regressors for each concept are used to assess these initialized solutions, and the corresponding classifier is used to predict their analysis success/failure. The initial parent population is sorted using a modified feasibility first ranking method, wherein the solutions with predicted analysis success are placed at the top, followed by the ones that are predicted to yield analysis failure. The solutions in each of the above blocks are ordered using feasibility first sorting. Differential Evolution with best/1/bin technique [55] is used to generate 10D offspring solutions from the sorted parent population. The offspring population is assessed using surrogates and combined with the parent population. The combined population is sorted using the same feasibility first sorting mechanism described earlier and the process of evolution is carried over for 10D generations. Finally, the top solution is selected as the first candidate solution for each concept. The process is repeated for all the concepts under consideration.

Population-based optimization algorithms are known to benefit in terms of convergence if marginally infeasible solutions are maintained in the population during the course of the search. IPB strategy is specifically aimed at exploiting this aspect wherein the top 2D solutions in a population are marginally infeasible solutions. The optimization process starts with 10D randomly initialized solutions, assessment of their analysis success/failure using the classifier, and evaluation of their performance through the predicted objective and constraints using the regressor. Solutions that are predicted to be analysis successful and infeasible are ranked using nondominated sort based on their objective function value and the sum of constraint violations (CV). The same process of nondominated sort is repeated for the block of solutions that are predicted to be failed in analysis and infeasible. This block of solutions is placed below the first block, and the top 2D solutions from this collection are carried forward as marginally infeasible solutions. To ensure the active participation of these solutions in recombination, they are placed at the top of the parent population. All the remaining solutions are ordered based on the feasibility first sorting method described earlier for PB. The population is allowed to evolve over 10D generations. Thereafter, three candidate solutions are selected from the final population. The first two are the top two solutions from the final population, i.e., one which has the least nonzero-sum of CV and the least objective function value. The third is the best feasible solution from the final population based on the feasibility first sorting described earlier.

Unlike the previous SAO strategies that rely on “believer” principles, the EcEI method uses predictor uncertainty along with the mean prediction. The maximization of EI has long been suggested and used in practice, often coupled with PF to deal with constrained optimization problems. The most widely adopted form relies on the maximization of $EI×PF$⁠, where EI is computed based on the predicted objective function value and its uncertainty along with the PF (for constrained problems). The probability of a solution being feasible is computed as a product of the probability of satisfying each constraint. In the event there are a large number of constraints, the PF term is often close to 0, complicating the ordering of solutions. The proposed EcEI method is a step towards resolving the “ties” that may occur in the process of ordering the solutions. The solutions are sorted based on the decreasing value of the product of $EI×PS×PF$⁠, where EI is the expected improvement of the objective function value, PS is the probability of analysis success, and PF denotes the probability of feasibility of the solution. In the case of ties, the following lexicographic sorting order is enforced where preference order is based on higher PS, higher PF, lower sum of CV, and finally, lower predicted objective function value (in a minimization sense). The best objective function value used in computing EI is the least objective function value of a feasible solution across all concepts. In the event, there are no feasible solutions evaluated so far, the expected improvement of the objective function for all solutions is set to 1. The process of ordering is used during the course of evolution over 10D generations and the top solution is selected as a candidate.

The fourth SAO strategy relies on the proposed BTPS. The BTPS strategy uses uncertainty in the objective function prediction, constraint prediction, and the probability of analysis success to rank solutions during optimization. The algorithm initially computes pairwise comparisons between all solutions to be ranked and then uses the Bradley–Terry model [56] to rank them.

We use Eqs. (2)–(8) to evaluate $PA≻B′$⁠, which is finally substituted into Eq. (1) to evaluate the probability (⁠$PA≻B$⁠), i.e., the measure of solution $A$ being better than solution $B$⁠. Equation (1) is used to perform a pairwise comparison of all possible solutions in the population and generate a score matrix $M$⁠, where $Mij$ is the probability $Pi≻j$ of a solution $i$ better than solution $j$⁠. We use the Bradley–Terry model [56] to estimate the relative strength of each solution and the decreasing order of these relative strengths provides the rank of the solutions.

Having described the principle of ranking solutions, we move on with its use in the course of SAO. The five candidate solutions obtained from the PB, IPB, and EcEI (described in the earlier sections) are used to seed the initial population of 10D solutions. Solutions are ranked based on the BTPS scheme discussed above and allowed to evolve over 50 generations with the DE/best/1/bin scheme. The top solution in the final generation is selected as the sixth candidate solution. The sixth candidate solution can be one of the five seeded solutions or a novel solution that emerged through the BTPS optimization process.

The SAO approaches for candidate generation outlined in the previous section deliver six candidate solutions per concept, and we need to identify one of those solutions for true evaluation. Assuming there are $N$ concepts, $6N$ candidate solutions are ranked using the BTPS sorting. Within the ranked candidates, the solutions that have already been evaluated are removed from the candidate list. A solution is designated as already evaluated if its Euclidean distance from any design point in the archive (within the same concept) is less than $1×10−4$ in the normalized variable space. The top solution from the resultant list is selected for true evaluation.

If all $6N$ candidate solutions exist in the archive (i.e., they have been evaluated in the past), a roulette wheel selection is undertaken to choose the concept to be evaluated. The concepts are ranked on the basis of their best-evaluated solution, which decides the size of the pie. Once the concept is selected, a single solution is generated using the DE/best/1/bin scheme from the archive of solutions sorted using the feasibility first principle.

The selected solution undergoes true evaluation, and the information is added to the archive. The above cycle is repeated until only $N$ evaluations are left from the total evaluation budget, where $N$ is the total number of concepts. For the final $N$ evaluations, we perform surrogate-assisted direct COBYLA optimization [57] on each concept. The optimization aims to minimize the objective function with constraints retrieved from the trained regressors for the concept. Additionally, the classification output defining the analysis success/failure is considered an additional constraint during the COBYLA optimization process, with analysis failure being considered an infeasible outcome. Each concept is optimized once using the surrogate-assisted COBYLA and the resulting solution is truly evaluated and added to the archive.

Although our approach is designed to deal with MCO, we first observe its performance on single-concept problems using 13 single-objective, inequality-constrained test problems from G-series [50]. Out of a total of 24 problems in the G-series, 13 problems contain only inequality constraints, which have been used in this study. The detailed description of these G-series problems, including the problem definition, constraints, variable bounds, and the best-known optimized solutions, is provided in Supplemental Material B. The problems contain a mix of linear and nonlinear constraints, a wide range of feasible regions, varying numbers of overall constraints (2–38), and active constraints (0–6). We have only used the DACE regressor to reduce overall computation time and benchmarked the performance against the GPSAF-ISRES algorithm using the same evaluation budget [51]. The overall statistics are presented based on 11 independent runs with different random seeds.

The statistical comparison between MCO and GPSAF using the Wilcoxon test is presented in Supplementary Table 1. The individual performance of both approaches is also presented in Supplementary Tables 2 and 3, respectively. The results indicate that the proposed approach offers statistically better results than GPSAF in two problems, equivalent results in two, and losses against GPSAF in nine problems. However, its important to acknowledge that GPSAF cannot deal with solutions having analysis failure or MCO problems. The surrogate of GPSAF operates as a predictor believer, and such approaches have been known to perform well across a range of functions. Additionally, we have incorporated the absolute coefficient of variability (⁠$|standarddeviation/mean|$⁠) metric to depict the overall variability in the MCO framework performance due to initialization and the algorithm’s own stochastic nature. The mean absolute coefficient of variability for G-series problems was 0.094, indicating minor variation in the performance on G-series problems. We observe from Supplementary Table S2 that most G-series problems showed an absolute coefficient of variability of less than 0.1. However, four problems, G2, G7, G9, and G19 showed a higher absolute coefficient of variability with the values of $0.15$⁠, $0.6$⁠, $0.14$⁠, and $0.27$⁠, indicating a larger impact of initialization and the stochastic nature of the algorithm on the performance. This indicated that the algorithm has not converged to its global minimum within the allocated evaluation budget for these problems.

To demonstrate the full functionality of the proposed approach in a single-concept setting, we modify the G24 two-variable, two-constraint problem. We use constraint $g1$ to represent analysis success/failure, wherein any true evaluation that violates the $g1$ constraint would result in the NaN outcome for the objective function value and the remaining constraint $g2$⁠. Additionally, we allow all possible regressors of the framework to be selected as opposed to only DACE used in the earlier example. We limit the total function evaluation budget to $300$ and use 2D evaluations for the initial LHS. We report the overall statistics of $11$ runs with different random seeds, and the graphs presented correspond to the median result.

The statistical results are presented in Table 1. We observe that the worst objective function value obtained from 11 runs is 0.005% higher than the ground truth. Additionally, the results show a very low absolute coefficient of variability (⁠$4×10−05$⁠), indicating an overall satisfactory outcome of the optimization exercise. To better visualize the search behavior in the variable space, we present regions that correspond to feasible, infeasible, and analysis failure in Fig. 4(a), along with the solutions sampled during the course of optimization. The analysis failure, infeasible, and feasible region constitute 20.48%, 35.39%, and 44.13% of the total design space, indicating a well-balanced test problem with a significant proportion of all three constituent regions. The samples are also labeled, i.e., whether it was suggested by PB, IPB, EcEI, or BTPS. The optimum lies on the boundary of all three regions, necessitating the need for an accurate classifier and accurate regressors for constraint and the objective. We observe that a large number of solutions were evaluated near the optimal location. Additionally, the solutions from all four constituent algorithms, i.e., PB, IPB, EcEI, and BTPS, were selected during the SAO process for true evaluation. Since BTPS starts with an initial population seeded by the best solutions of PB, IPB, and EcEI, a solution is labeled as BTPS only if it is different from the seeded candidate solutions. One can observe from Fig. 4(b) that the three constituent surrogates (classifier, constraint regressor, and objective regressor) show high accuracy in the region of interest.

With the basic working principles illustrated through single-concept constrained optimization benchmarks and the problem with analysis failure, we advance to MCO problems. We use a modified version of the multiconcept beam design optimization problem presented in a recent work [3]. The objective is to minimize the overall weight of a cantilever beam with a length 400 mm, Young’s modulus $200×103N/mm2$⁠, density $7820×10−6kg/mm3$⁠, subject to the end load of $600N$⁠, with a constraint that the maximum deflection of the beam is limited to 5 mm. The beam can have six different cross sections, each represented using a different number of design variables and variable bounds, as shown in Fig. 5. We solve the MCO problem with a limited computation budget of 480 function evaluations. We use 20D evaluations for the initial LHS and utilize the remaining 240 evaluations for optimization. We report the overall statistics of 11 runs with different random seeds, and the presented graphs correspond to the median run.

The statistical results for the MCO-beam problem are presented in Table 2. In this problem, the I-beam concept corresponds to the optimal solution. We observe that the proposed approach obtained a median objective function value of 264.72, which is 0.1% higher than the ground truth that has the objective function value of 264.43 (at location x51: 28.85 mm, x52: 1 mm, x53: 1 mm). This indicates that the algorithm was able to identify the correct concept along with its variables in the majority of runs. We also see a high absolute coefficient of variability of 0.38. This high variability is attributed to the framework’s inability to identify the I-beam as the best design for a few random seeds, leading to higher objective function values. For the median case, we observe from the convergence plot in Fig. 6(a) that the approach identified I beam as the best concept even when it had the worst initial performance among all concepts during the LHS phase. The optimum solution (objective function value: 264.72) is achieved within the first 20 optimization iterations. It should be noted that the achieved solution is already close to the global optimum solution (objective function value: 264.43). As the global optimum is identified in the initial phase of the optimization process, further reduction during the later phase is not expected and, therefore, not observed from the convergence plot. Additionally, the approach optimized other concepts once the minimum for the I section was identified. We also observe from the cumulative evaluation plot in Fig. 6(b) that although the algorithm performed the highest sampling of the I section, it also heavily sampled other promising concepts (H and square hollow), indicating the dynamic concept selection ability of the approach. The above results substantiate the ability of the approach to deal with MCO problems.

In the previous section, we demonstrated the performance of our approach using a simple multiconcept problem involving beam design. Here, we present a practical optimization problem focused on 3D shape matching across four different stent designs: two independent ring (IR1, IR2) stents and two helical stents (HS1, HS2). Before introducing the 3D shape-matching problem, we present the methodology for the geometric representation of a realistic single-helix HS2 stent which has never been attempted before. The IR1 and HS1 designs were parameterized using 17 and 18 independent variables, respectively, with a detailed description available in our earlier work [26]. The IR2 stent resembles an old Cypher (Johnson & Johnson Inc.) design which has been used in past literature [58]. We also parameterized the IR2 stent using 9 independent variables and the detailed parameterization scheme of the IR2 stent is provided in Supplemental Material C.

The geometry of the single-helix HS2 design is parameterized using five independent variables, as presented in Fig. 7. The stent structure consists of four primary components - the main helix, transition helix-1, transition helix-2, and the end helix. Three independent variables, i.e., crimped stent radius (⁠$CR$⁠), stent length (⁠$SL$⁠), and stent strut radius (⁠$SR$⁠)) define the overall characteristics of the entire stent. The fourth variable, $Nn$ (number of unit struts) defines the total number of unit struts in the main helix, and the fifth variable $Ry$ (axial semiaxis) defines the axial semiaxis length of the crown. The entire stent structure consists of a single strand of circular wire extending from one end to the other. The baseline values of these variables are listed in Table 3.

We model the main helix geometry by repeating the unit strut, with each new unit shifted by $AAX$ and $AAY$ distance in the circumferential and axial directions, respectively. The unit strut comprises five unit strut components, with each being a sequence of a bottom crown, a short strut, a top crown, and a tall strut. The bottom and top crowns are half-ellipses with circumferential and axial semiaxis lengths $Rx$ and $Ry$⁠, respectively. The short strut is a straight vertical line segment of length $SSy$⁠, and the tall strut is an inclined line segment with length $TSy$ and width $TSx$⁠. We use PARSEC connections [49] between crowns and tall strut to maintain $c0$ and $c1$ continuity. The difference in the strut length results in the helical nature of the stent and generates the helix angle (⁠$α$⁠). This stent does not have any separate connectors and the axial connection is provided by one top and one bottom connector crown within the unit strut. A connector crown has the same circumferential semiaxis length (⁠$Rx$⁠) but a larger axial semiaxis length (⁠$Ry′$⁠) than the other crowns. Each bottom connector crown is joined with the top connector crown of the next unit strut. We model the end helix by maintaining a 0-deg helix angle at the end side and a transition helix angle (⁠$β$⁠) for the first six strut components at the helix side. The next three strut components link the end helix to the transition helix-1 with the end helix angle (⁠$θ$⁠) at the end side and the helix angle (⁠$α$⁠) at the helix side. The length of the first short strut (⁠$Sy0$⁠) is the final dependent variable required for the overall construction of the helix. The other lengths of the end helix struts, transition helix-1 struts, and transition helix-2 struts are derived from the $Sy0$⁠, $SSy$⁠, $α$⁠, $β$⁠, and $θ$ angles. All standard and connector crowns across the end helix, transition helices, and main helix have the same dimensions.

The dependent variables $TSy$⁠, $SSy$⁠, $TSx$⁠, $Rx$⁠, $Sy0$⁠, $θ$⁠, and $Ry′$ are analytically obtained by solving Eqs. (12)–(18). In addition to these seven, we still need the values for three other dependent variables, $α$⁠, $SLT$⁠, and $β$⁠. We perform the repair procedure as defined in our earlier parameterization of the double helix HS1 [26]. The three dependent variables are obtained by solving a single-objective constrained optimization problem. We set the initial values based on the contemporary stent designs to be [23.5 deg, 6.7 mm, and 11.75 deg] and minimize the absolute difference between $Sy0$ and $SSy$⁠, subject to constraints listed in Eqs. (19)–(21). The variables are allowed to vary between $±50%$ of the initial values.

Similar to IR1 and HS1 parameterizations, we created a 3D stent geometry builder by coupling matlab 2023b (Mathworks Inc.) programming interface with Solidworks 2024 API (Dassault Systemes SE) through python 3.9 and Pywin32. The 3D stent geometry builder would be fed with independent design variables (rounded to five decimal places) as inputs, and it would, in turn, generate a neutral parasolid CAD file and a stereolithography file of the resultant stent geometry. Building a stent geometry takes approximately 434 s on a workstation with Intel(R) Xeon(R) Gold 6226R CPU @ 2.90 GHz, 2.89 GHz processors, and 128 GB RAM.

Shape matching is often performed as the first step before undertaking a comprehensive design optimization exercise involving FEM analysis or CFD tools. A shape-matching optimization exercise aims to identify a geometry that resembles a given target geometry. We solve a multiconcept shape-matching optimization problem where the target geometry is the baseline design of IR1 presented in Ref. [26]. The stent geometry builder created the target stent mesh with $28,697$ points. The objective function for shape matching was set as the norm of the distances between each point in the target design geometry and its nearest point in the candidate design. The candidate design can be from any of the four different stent designs. The variables were allowed to vary $±10%$ of their baseline design values. For the IR1 stent, the number of connectors was allowed a variation of 33% to accommodate discrete values. Similarly, for the IR2 stent, a higher 50% variation was allowed for the connector $ρ1$ and $ρ2$ variables to allow significant changes in the oscillating behavior of the stent’s connector design. The optimization process was started with 20D LHS evaluations per concept, which meant 340 (IR1), 180 (IR2), 360 (HS1), and 100 (HS2) designs were evaluated during the initialization phase. The total evaluation budget was set to 1180 evaluations which meant the algorithm could evaluate 200 additional designs beyond the initialization phase. It is noteworthy that the problem is challenging since analysis failure is fairly common. Within the total LHS evaluations, 21%, 7%, 43%, and 35% evaluations for IR1, IR2, HS1, and HS2 stents resulted in geometry construction failure.

The results of 11 independent runs with different random seeds are presented in Table 4. We compare the results obtained from our proposed approach with single concept (SC)-PB and EI shape-matching methods reported in our earlier work [26]. For the same additional function evaluation budget of 200, the proposed MCO approach resulted in the objective function value of 328.9, an improvement of 57.36% over SC-PB and 21.05% over SC-EI. This indicates that an optimal balance between the exploration and exploitation components is achieved by the proposed MCO framework. Furthermore, the proposed framework also obtained the best overall mean, median, and minimum objective function value over the SC-PB and SC-EI approaches. This demonstrates the successful combination of constituent PB and EI algorithms in the MCO framework with improved minimum objective function value (PB algorithm feature) and mean/median objective function values (EI algorithm feature). We observe a consistent reduction in the objective function value in the convergence plot in Fig. 8(a), even in the later stages of the optimization process, with the minimum objective function value of 328.9 achieved after approximately 170 optimization iterations. This consistent decline is observed as the global optimum objective function value (0) was not achieved within the provided evaluation budget. Additionally, the limited reduction in the objective function value of other stents indicates the framework’s ability to identify and optimize the promising concept. This is also evident from the pie chart in Fig. 8(b), which shows that IR1 stent received the maximum sampling of 45.5% during the complete optimization process (including the concept-wise initialization phase). The computational time comparison of the MCO framework with SC-PB and SC-EI approaches for a typical run is presented in Supplementary Table 5. The MCO approach exhibits higher wall clock time, surrogate creation time, and optimization time than the SC-PB and SC-EI approaches. The higher surrogate creation time is due to fitting 23 classifiers and 41 regressors with hyperparameter tuning in MCO for each of the 4 concepts, as opposed to 6 classifiers and 8 regressors for a single concept in SC-PB and SC-EI. Similarly, the MCO framework shows higher optimization time because four SAOs are performed across four concepts compared to one SAO for a single concept in SC-PB and SC-EI approaches. Among the four SAO techniques, believer approaches (PB and IPB) have the least computational cost (1859 s and 1872 s, respectively) followed by EcEI (12,291 s) due to bootstrapping, and BTPS (52,826 s) due to expensive generation of score matrix based on pairwise comparisons. The comparison between the target geometry and the best designs obtained from SC-PB, SC-EI, and the proposed MCO framework is shown in Fig. 9. The SC-PB design was trapped in local minima as shown by the four-connector design as compared to the three-connector design of the target. It also showed a longer length than the target. The results from SC-EI and the proposed MCO approach closely resemble the target design, with the proposed approach demonstrating a closer match with the target in strut dimensions and the overall length of the stent. Although in the G-series constrained test problems presented earlier, PB performed significantly better as a part of GPSAF [51], they are known to be susceptible to being trapped in local optimum as evidenced in this example. Finally, we performed independent baseline optimization for each concept with an equal function evaluation budget (50 each) and combined them to compare the results with the proposed MCO approach with its adaptive budget allocation mechanism. The detailed results are presented in Supplementary Table 6. The MCO framework achieved improved median (45.8% lower), minimum (⁠$45.7%$ lower), and mean (34.7% lower) objective function values than the combined optimization. This highlights the superior performance of the proposed MCO framework as compared to the simpler strategies, such as combining independent single-concept optimization with equal function evaluation budgets.

In this article, we introduced a framework to deal with unconstrained/constrained, single-objective optimization problems involving computationally expensive analysis, wherein a search is conducted across multiple concepts defined using its own set of variables. Since most optimization methods operate with a fixed set of variables, they cannot be directly used to solve such problems. Four SAO algorithms have been used in the framework, PB and IPB, along with two novel algorithms, i.e., EcEI and BTPS. Each of these algorithms independently evolve a population of solutions and are equipped to handle analysis failure. Their search strategies exploit various principles, e.g., (1) relying on surrogates as accurate approximators (PB), (2) preserving marginally infeasible solutions (IPB) during the course of the search to achieve convergence to the constrained optima, (3) utilizing uncertainty information of the surrogates via expected improvement and lexicographic sorting (EcEI), and (4) employing probabilistic model to rank and evolve solutions (BTPS). To exploit the benefits of various functional representations, the algorithms dynamically choose surrogates from a suite of 23 classifiers and 41 regressors. In every step, a single solution is evaluated which conforms to a steady-state paradigm.

Several analytical and practical examples have been used for illustration and benchmarking. First, we observed the performance of the approach on 13 single-concept G-series benchmark problems and compared it with the results obtained from the state-of-the-art GPSAF algorithm. Second, the ability of the approach to solve optimization problems involving analysis failure was demonstrated using the modified G24 problem. The approach successfully handled analysis failure, generated high-accuracy surrogate models near the region of interest, and successfully identified the optimal solution. Having established the performance on single-concept optimization problems using the above two examples, we focused on a cantilever beam design problem involving six concepts, wherein different concepts correspond to different cross-sectional shapes of the beam. Our proposed approach dynamically allocated search effort across the concept spaces and quickly identified the best concept along with the corresponding design. To demonstrate the performance in a more complex MCO problem with analysis failure, we presented a 3D shape-matching problem involving four coronary stent concepts. For this task, we used existing geometry representation schemes for three of the stents, i.e, IR1, IR2, and HS1. Notably, we introduced, for the first time, a geometry parameterization scheme for single-helix stent designs using a combination of PARSEC, ellipse, and straight line segments. The results from the 3D shape matching were objectively compared with single-concept [believer (PB), and BO (EI)]-based methods to substantiate the benefits of the proposed approach. It is important to take note that while we have used “shape matching” as the objective in this study, other domain-specific objectives derived from CFD or FEM simulations could be easily used.

There are a few areas for further development. The first is to extend the current capability of the approach to deal with multiple objectives. This would certainly increase the uptake of the approach for practical applications. The second is in the context of shape and topology design where a more radical advancement would be to step out of the family of concepts to uncover novel shape/topologies using deep neural network-based representations coupled with innovative optimization strategies. These representations can generate completely free-form geometries, unrestricted by any parameterization scheme. The optimization of these representations can theoretically allow the creation of any possible geometry. However, there are multiple challenges that require further research for the concept free design optimization. These include using domain-specific objectives instead of a shape-matching metric. While shape matching can be easily swapped with the domain objective in the presented MCO framework, this swap is nontrivial in design optimization using neural networks as the underlying Adam or stochastic gradient descent optimizers require the objectives to be auto-differentiable for managing a large number of neural network weights. Since the domain-specific objectives are obtained using computationally expensive physics-based simulations, generating them in an auto-differentiable setting is not straightforward. Similarly, these underlying optimizers are designed to deal with unconstrained single-objective optimization and significant modifications will be required to deal with the constraints and multiple objectives. Third, there is an inherent dependency of the optimization results on the initialization and stochastic nature of the population-based optimization algorithms. The development of novel algorithms to minimize this dependence is an ongoing topic of research in the field of evolutionary computation. The authors are currently pursuing some of these directions.

A.K. gratefully acknowledges the support from the Commonwealth Government through the Australian Government Research Training Program Scholarship. S.B. would like to acknowledge the support from the National Heart Foundation Vanguard grant.

NJ declares grant and honoraria support from Abbott Vascular. The other authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The stent 3D geometry builder codes for IR2 and HS2 stent are publicly available (https://doi.org/10.5281/zenodo.15117426, GitHub repository: https://github.com/ankushkapoor2003/stent_geometry_builder.git).