Globalized parametric optimization of microwave components by means of response features and inverse metamodels

Simulation-based optimization of geometry parameters is an inherent and important stage of microwave design process. To ensure reliability, the optimization process is normally carried out using full-wave electromagnetic (EM) simulation tools, which entails significant computational overhead. This becomes a serious bottleneck especially if global search is required (e.g., design of miniaturized structures, dimension scaling over broad ranges of operating frequencies, multi-modal problems, etc.). In pursuit of mitigating the high-cost issue, this paper proposes a novel algorithmic approach to rapid EM-driven global optimization of microwave components. Our methodology incorporates a response feature technology and inverse regression metamodels to enable fast identification of the promising parameter space regions, as well as to yield a good quality initial design, which only needs to be tuned using local routines. The presented technique is illustrated using three microstrip circuits optimized under challenging scenarios, and demonstrated to exhibit global search capability while maintaining low computational cost of the optimization process of only about one hundred of EM simulations of the structure at hand on the average. The performance is shown to be superior in terms of efficacy over both local algorithms and nature-inspired global methods.

Topological complexity of passive microwave components has been continuously increasing over the years 1,2 . This is a consequence of growing performance demands 3 , functionality requirements [4][5][6][7] , but also miniaturization trends 8 . In the latter case, techniques such as transmission line (TL) folding 9 or slow-wave phenomenon 10 , are often employed, leading to geometrically involved structures described by many parameters 11,12 . Circuittheory-based methods often turn out to be inadequate in describing the intricacies of such devices. Due to electromagnetic (EM) cross-coupling and similar effects, their reliable evaluation can only be realized through full-wave EM analysis.
As a result of involved interrelations between the circuit topology and its electrical characteristics, simultaneous optimization of geometry parameters by means of numerical algorithms becomes imperative to achieve the best possible performance of the structure. In fact, numerical optimization allows proper handling of several objectives and constraints over multi-dimensional parameter spaces. Yet, it is an expensive process as even local optimization involves a considerable number of system evaluations. In many cases, including multimodal tasks [13][14][15] , multi-criterial design 16 , or the lack of good starting points 17,18 , globalized search is necessary, which makes the optimization problem even more challenging.
Undoubtedly, the most popular global optimization methods today are population-based nature-inspired algorithms [19][20][21] . Their roots can be tracked back to late 1960s 22 , and eventually dominated global search practice since 2000s [23][24][25][26][27] . Over the last years, the number of nature-inspired algorithms has been growing tremendously (firefly algorithm 28 , harmony search 29 , and others [30][31][32][33][34][35][36] ). Population-based methods capitalize on exchanging information between the members of the candidate solution set [37][38][39] , but also generating new data using exploitative operators 40 . Avoiding local minima is facilitated by including randomness in various forms 41,42 . Typically, nature-inspired algorithms are straightforward to implement, yet, their computational efficiency is poor: a single optimization run may require from a few hundreds to many thousands of objective function evaluations, which becomes a serious problem when the system of interest is to be evaluated using full-wave EM simulation.

Globalized microwave optimization using feature-based inverse metamodels
The purpose of this section is to introduce the optimization technique discussed in the work. Our approach employs inverse regression surrogates established using pre-selected random parameter vectors, and the characteristic points of EM simulated responses of the microwave component under design. Weakly-nonlinear relationship between these feature points and geometry parameters enables global search capability at low computational cost. The inverse model is used to render a good starting point, which is subsequently tuned by means of a local (here, gradient-based) procedure.
Formulation of EM-driven design task. The simulation-driven design problem is formulated here as a nonlinear minimization task of the form where U is a scalar objective function, and F t = [F t.1 … F t.K ] T is a target vector of operating parameters. The objective function quantifies the quality of the design based on EM-simulated responses of the microwave component at hand, which are most often scattering parameters S kl (x,f), where k and l denote the corresponding ports of the circuit, x is a vector of designable parameters, and f is the frequency.
For the sake of example, let us consider a microwave coupler, which is to operate at the frequency f 0 so that its matching and isolation characteristics, |S 11 | and |S 41 |, are minimized at f 0 , and the power split ratio d S (x,f 0 ) =|S 21 (x,f 0 )| -|S 31 (x,f 0 )| reaches a target value K P (e.g., 0 dB for equal power split). In this case, the operating parameter vector would be F t = [f 0 K P ] T , whereas the objective function may be defined as (1) x * = arg min x U(x, F t ) www.nature.com/scientificreports/ where the second term is a penalty function enforcing the required power split ratio, with β being the penalty factor controlling the contribution of the penalty term to the overall objective function. Another example is to design a dual-band coupler for a substrate characterized by a specific relative permittivity ε r , so that the circuit minimizes both |S 11 | and |S 41 | responses at the two operating frequencies f 0.1 and f 0.2 , while providing equal power split at these frequencies. Here, the operating parameter vector is F t = [f 0.1 f 0.2 ε r ] T , and the objective function may be defined as Other design scenarios can be treated in a similar manner. This includes cases where the operating bandwidth of the device is handled explicitly (e.g., bandwidth enhancement, etc.).
Design quality evaluation using response features. Adjustment of geometry parameters is a necessary step of microwave design process. It aims at improving the performance parameters, and, as explained in "Formulation of EM-driven design task", it can be formulated as an optimization task, normally solved at the level of EM analysis to ensure reliability. In this work, we address globalized optimization. It is often required, either due to the lack of good starting point, or the presence of multiple local optima, some of which may fail to satisfy the prescribed performance requirements. A representative situation is a design of miniaturized structures where conventional transmission lines are replaced by CMRC or similar unit cells 10 . Therein, the relationships between geometry parameters and electrical characteristics of the cell are generally complex 12 , which renders an identification of a reasonable initial design a difficult problem. Similar issues may arise when re-designing a given circuit for operating frequencies or substrate that are away from those at the current design.
Global exploration of the parameter space is a daunting task due to nonlinearity of system characteristics (both as a function of geometry parameters and frequency), but also dimensionality issues. While direct EMdriven global search using, e.g., nature-inspired algorithms, is almost always computationally prohibitive (unless the computational model is relatively cheap to evaluate), the aforementioned reasons also hinder utilization of surrogate-assisted procedures, as rendering reliable metamodels is rarely feasible beyond a few parameters and within narrow ranges thereof. Figure 1 illustrates several situations where a local (e.g., gradient-based) search may fail due to the lack of a good initial design or the necessity of re-designing the structure for operating conditions that are distant from those at the current design. The coupler of Fig. 1a is a compact microstrip rat-race coupler, described by six independent geometry parameters x = [l 1 l 2 l 3 d w w 1 ] T , further details pertaining the circuit can be found in Table 2. The example is based on one of the miniaturized coupler structures considered as verification cases in "Numerical verification". Figure 1. Miniaturized microstrip coupler and its scattering parameters versus frequency: (a) coupler geometry, (b) S-parameters at selected random designs within the assumed parameter space. The vertical lines mark the target operating frequency (here, 1.6 GHz). Local search carried out using the objective function such as (2) would fail when starting from most of the shown designs, due to severe misalignment between the target and the actual operating conditions. www.nature.com/scientificreports/ The issues discussed above may be mitigated by processing information extracted from the EM-simulated circuit responses in the form of the characteristic points (or response features), which is a foundation of featurebased optimization (FBO) technology 61 . FBO explores the fact that despite intrinsic nonlinearity of the system responses (cf. Fig. 1b), the relationships between the characteristic point coordinates and geometry/material parameters is much less nonlinear, which allows for obtaining a considerable amount of information about the system using a limited amount of EM simulation data, as shown in Fig. 2. Figure 2a shows the selection of the response features corresponding to minimum of matching and isolation characteristics as well as the power split ratio of miniaturized microwave coupler of Fig. 1a.
The power split K P =|S 21 (f 0 )| -|S 31 (f 0 )| is evaluated at the approximate operating frequency f 0 of the circuit. The frequency f 0 is assessed as the average of the |S 11 | and |S 41 | minima. Note that some of the feature points may not exist depending on the particular parameter vector (e.g., the operating frequency located outside the simulation frequency range). The relationship between the operating frequency and power split ratio and the a b Figure 2. Miniaturized coupler of Fig. 1a: (a) response features: minima of |S 11 | and |S 41 | and power split ratio K P (o); K P is evaluated at the frequency (approximate operating frequency f 0 of the circuit) being the average of the said minima (thick vertical line); (b) relationship between f 0 and K P and the three selected geometry parameters; the circles mark coupler designs and the gray points denote the regression model. www.nature.com/scientificreports/ three selected geometry parameters is presented in Fig. 2b, where the specific coupler designs are shown along with the regression model of the form a 0 + a 1 exp(a 2 f 0 + a 3 K P ), which represents the trends between the operating parameters and the circuit dimensions. The definition of the feature points depends on the particular shape of the circuit responses and on the formulation of the design task. These could be simply the frequency and level locations of the resonances 61 , local minima/maxima of the pass-band part of the return loss 62 , or points defining a circuit bandwidth, power split, etc. 64 . In this paper, utilization of characteristic points is one of the foundations of the presented optimization technique, specifically, at its first stage ("Globalized optimization with inverse regression models"). However, the feature points are used here primarily to estimate the actual operating conditions of the circuit, rather than directly (as in FBO 61 ), which was graphically illustrated in Fig. 2a.
Globalized optimization with inverse regression models. As announced in "Design quality evaluation using response features", this paper capitalizes on a weakly-nonlinear dependence of the geometry variables of the circuit on its operating parameters (frequency, bandwidth, power split) in order to explore the parameter space in a computationally-efficient manner. Examples of such relationships can be found in Fig. 2b for a representative miniaturized microstrip coupler. By weakly nonlinear we mean the type of relation that is usually monotonic (e.g., of the exponential type, close to proportional or inversely proportional).
This sort of relationship holds for many practical circuits, especially when the figures of interest are operating frequencies/bandwidth, or material parameters (substrate permittivity/height), but also other parameters (e.g., power split ratios for couplers). In particular, operating parameters are in typically monotonic relations with certain major parameters, despite the fact that the entire frequency characteristics may be strongly nonlinear function of frequency.
The information necessary to estimate the mentioned dependencies is acquired using randomly generated parameter vectors (observables). Some of these may be of good quality from the perspective of the assumed performance requirements, whereas others may be poor and need to be rejected. The observable quality is evaluated using the feature points extracted from the simulated scattering parameters, and comparing the utility metrics obtained this way with the design targets. The subset of the best observables is then employed to identify an inverse regression model. The latter serves two purposes: (1) to find the promising parameter space region, and (2) to generate infill points for refining the inverse surrogate. The globalized search process is executed iteratively, with a single infill point rendered per iteration, and the worsts observables replaced by those being closer to the target. The details of the procedure are explained in the remaining part of this section.
First, we introduce the notation used throughout: ] T -a vector of the operating parameters at the design x (e.g., centre frequency, bandwidth, power split ratio), extracted from the EM simulated circuit responses. As mentioned before, the particular operating parameters are estimated based on the feature points (cf. Fig. 2a); e.g., the operating frequency of the coupler can be estimated as the average of the frequencies corresponding to the minimum of the matching and isolation characteristics. If some of the parameters cannot be extracted (e.g., some of the relevant feature points cannot be distinguished or are allocated outside the frequency range of simulation), we assign ] T -a vector of auxiliary coefficients reflecting the design quality and corresponding to the entries of the vector F(x). For example, if the objective is to reduce the level of |S 11 | and |S 41 | at the operating frequency, the corresponding l k can be the average of |S 11 | and |S 41 | at their respective minima: the lower value indicates that the design x is of higher quality. Similarly, if the operating parameter is a power split ratio, the corresponding l k might be a deviation from the estimated power split and its target value. If some of the entries of L(x) cannot be extracted, we assign L(x) = [0 … 0] T ; • D(F,F t )-a function quantifying the misalignment between the target vector F t (cf. "Formulation of EMdriven design task") and the operating parameter vector F; in this work, we use • D accept -user-defined control parameter utilized to terminate the global search stage of the optimization process, i.e., we assume that the current design is sufficiently close to the target if D(F,F t ) ≤ D accept .
Before providing a rigorous formulation of the global search process, the following outline is discussed to clarify the operation and the meaning of the specific steps: 1. Observable generation obtain a set of parameter vectors x (j) , j = 1, …, N, generated randomly over the assumed space X (most often, an interval defined by the lower and upper parameter bounds), typically, using a uniform probability distribution. The vectors are generated as long as necessary to obtain N designs for which ||F(x (j) )||> 0, j = 1, …, N. 2. 2. Inverse model construction Using the set of triples {F(x (j) ), L(x (j) ), x (j) } j = 1,…,N , identify an inverse regression model r I (F) with the values in X; the model quantifies the dependence between the operating and geometry parameters of the circuit. The analytical formulation of r I will be discussed later in the section; 3. Design prediction Employ the inverse model r I to identify a candidate parameter vector x tmp = r I (F t ), where F t is the vector of target operating parameters (cf. "Formulation of EM-driven design task"). If ||F(x tmp )||> 0 and D(F(x tmp ),F t ) < max{j = 1, …, N : D(F(x (j) ),F t )}, replace the vector realizing the above maximum by x tmp , and reset r I . www.nature.com/scientificreports/ The second and the third step are iterated in the attempt to find a design that is sufficiently close to the target. More specifically, the procedure is terminated if D(F(x tmp ),F t ) < D max (a user-defined acceptance threshold). The next stage is local optimization as described in "Local optimization procedure". It can be noted that the procedure generates random parameter vectors until a sufficient number of designs are found for which clearly defined feature points can be extracted. This is followed by constructing the inverse model (see Fig. 2b for a graphical illustration), which is then used as a predictor to yield a location of the design for which the operating parameters are possibly close to the target F t . If the candidate design is of sufficient quality (according to function D(F,F t )), it replaces the worst of the existing base vectors.
Because the design replacement in the dataset {x (j) } is governed by the proximity function D(F,F t ), over time, the inverse model will be focused on the region containing designs that exhibit low values thereof. For the same reason, the local accuracy of the surrogate will gradually improve. A graphical illustration of the procedure can be found in Fig. 3. The acquisition of the observables is carried out as follows: only the designs with their corresponding operating parameters (centre frequency, power split ratio) within the region of interest and the simulation frequency range will contribute to a construction of the inverse model r I (see Fig. 3a).
The inverse regression model r I (F) is the fundamental component of the proposed optimization procedure. As outlined above, it is constructed using the triples {F(x (j) ), L(x (j) ), x (j) } j = 1,…,N . The analytical form of r I can be simple because the dependence between the geometry parameters and the operating conditions of the circuit is typically only weakly nonlinear. Nevertheless, the model has to have a sufficient flexibility to account for the fact that the aforementioned dependence may be close to inverse proportionality for certain parameters. Having this in mind, the following form has been assumed: The surrogate is identified by solving "Good" designs (operating parameters within the frequency range of interest) "Poor" designs (operating parameters outside the frequency range of interest) X Random observables in parameter space Here, it is assumed that the components l k are normalized, and can only assume the values from the interval [0, 1], with zero corresponding to the highest-quality design (with respect to the kth operating parameter), and one corresponding to the lowest-quality design. In the example considered before, with the objective being to reduce the level of |S 11 | and |S 41 | at the operating frequency, the corresponding l k could be selected as the average of |S 11 | and |S 41 | at their respective minima. Clearly, good design corresponds to l k close to zero (low reflection and high isolation), whereas poor design would be associated with l k closer to one. Overall, the inverse regression model r I is essentially a trend function that approximates the observable set {x (j) } in the weighted L-square sense (cf. (2)). The reason for introducing the weights w k is to discriminate between the low-and high-quality observables so that the latter affect r I to a greater extent.
It should also be noted that, in general, construction of the inverse models may be hindered by non-uniqueness issues (e.g., Refs. 71,72 ). Notwithstanding, for typical microwave passive components, the operating conditions (e.g., operating frequency, bandwidth, etc.) are mainly dependent on specific geometry parameters controlling electrical lengths of its parts, therefore, the relationship between designable parameters and operating conditions is usually monotonic. Although it might not be so for compact structures (especially those utilizing slow-wave phenomenon), the mentioned major parameters enforce monotonicity. Furthermore, the inverse model is established as a regression surrogate, so it does not follow exactly its training data but only the trend. Consequently, it accommodates possible non-uniqueness due to the parameters of minor importance (from the point of view of the trends). Finally, it should be noted that identification of the inverse model is a weighted regression task (cf. (5)), with less weight put on "poor" observables, thereby, extracting the trends from the best points only. Figure 3b provides a graphical illustration of the inverse model for the microstrip coupler of Fig. 1a. Whereas the conceptual illustration of the global search process is presented in Fig. 3c,d. In the first iteration of this search (cf. Fig. 3c), the infill point predicted by r I replaces the worst observable and the model r I is updated. Figure 3d presents the last iteration, in which the observables are concentrated near the target operating parameters, hence the updated inverse model is capable of yielding the design which is sufficiently close to the target. Hence, the procedure may be terminated and followed by a local optimization ("Local optimization procedure"). Observe that Fig. 3b-d refer to a single geometry parameter x; the same scheme is applied to all parameters simultaneously.
The operating flow of the global search process as proposed in this work has been summarized in Fig. 4 in the form of a pseudocode. Therein, Steps 1 through 4 correspond to identification of N observables with their operating parameters being within the prescribed ranges, in particular, the frequency-related parameters being within the range of circuit simulation. These parameter sets are employed in Step 5 to construct the inverse model r I . The remaining steps describe utilization of r I for generating a candidate design x tmp , its evaluation, and insertion into the observable pool (provided it is of sufficient quality). These are followed by rebuilding the inverse model.
The termination condition is D(F(x (0) ),F t ) ≤ D accept , i.e., identification of a design, which is sufficiently close to the target F t , which then becomes a starting point x (0) for local optimization (cf. "Local optimization procedure"). If such a design cannot be found, the procedure is terminated upon exceeding its computational budget, which results in returning the best design found so far.
Local optimization procedure. Global search procedure, as described in "Globalized optimization with inverse regression models", yields a design x (0) satisfying the condition D(F(x (0) ),F t ) ≤ D accept , with the threshold D accept assigned to make sure that the operating parameters at x (0) are sufficiently close to F t to make the target attainable by means of local optimization. In this work, it is realized using the trust-region (TR) gradient-based algorithm with numerical derivatives 65 . The TR algorithm produces a series of approximations to the optimum design x * , denoted as x (i) , i = 0, 1, … The subsequent iteration points are obtained as where the objective function U L takes the same form as the function U (cf. (1)); however, it is computed using of the first-order Taylor expansion model G (i) (x,f) of the system responses established at the current point x (i) . The linear model is defined, for the S-parameter S kl , as The gradients in (8) are estimated using finite differentiation. The trust region in (7) is an interval [x (i) -d (i) , x (i) + d (i) ]. The size vector d (i) is adjusted using the standard TR rules 65 . The candidate vector x (i+1) is accepted if it reduces the objective function value at the EM simulation level, i.e., if U(x (i+1) ,F t ) < U(x (i) ,F t ). Otherwise, it is rejected and the iteration is repeated with reduced d (i) .
The algorithm termination is determined by the convergence in argument ||x (i+1) -x (i) ||< ε, or diminishing the TR size, i.e., ||d (i) ||< ε (whichever occurs first). Here, we use ε = 10 -3 . In order to reduce the computational cost of the optimization process, finite differentiation (normally entailing n additional EM analysis of the circuit (5) p j.0 p j.1 ... p j.K+1 = arg min Global optimization framework. The operating flow of the globalized optimization framework discussed in this work is summarized below (see also Fig. 5). Its two main components are the global and local optimization procedures formulated in "Globalized optimization with inverse regression models" and "Local optimization procedure", respectively. The framework uses the control parameters gathered in Table 1. Note that only the first two parameters (primary parameters of Table 1), i.e., N and D accept , are specific to the proposed technique, whereas the remaining ones are conventional (in terms of numerical optimization routines). Parameter N (i.e., the number of observables utilized for inverse model setup) should be of the same order as design space dimensionality but also take into account the number of design objectives in order to properly account for the geometry of the set comprising high-quality designs. An appropriate value of D accept is problem-specific and should be set to make the target operating parameters attainable by means of a local algorithm (at the local optimization stage). In practice, operating parameters most often relate to the operating frequency (or frequencies) of a component under design. In such a case, it suffices to set D accept equal to approximately half of the intended operating bandwidth, which, in turn, typically ranges from a dozen to few dozen percent of the operating frequency. Whereas the parameters N max.k , k = 1, 2, 3, should be set up with some margin in order to ensure the algorithm termination due to convergence rather than exceeding the budget. Thus, N max.1 and N max.2 should be of one order of magnitude larger than the design space dimensionality, and N max.3 a few times higher. In short, the operating flow of the entire optimization process can be described using the following three stages: 1. Input argument setup: • Target operating frequencies f t , • Objective function U, • Parameter space X; 2. Global search: Obtain initial design x (0) by performing the algorithm of "Globalized optimization with inverse regression models"; 3. Local optimization: Find the final design x * using the TR algorithm of "Local optimization procedure". www.nature.com/scientificreports/ The flow diagram of the process can be found in Fig. 5, where the first stage is broken down into several components, whereas the local refinement is represented by a single block. It should be emphasized that in this work, the global search procedure does not involve any direct global optimization algorithm. Instead, identification of the approximate allocation of the globally optimum solution is obtained from the inverse model. Its evaluation at the target objectives directly brings us to the appropriate portions of the parameter space. The key factor behind the efficacy of this approach is that the objective space is normally of considerably lower dimension than  www.nature.com/scientificreports/ the parameter space, which allows us to construct a relatively accurate inverse model using a limited number of observables. These operating principles are the major differences between the methodology proposed in this work and the majority of approaches to global optimization of expensive simulation models reported in the literature.

Numerical verification.
Here, the globalized optimization framework introduced in "Globalized microwave optimization using feature-based inverse metamodels" is validated and its performance is demonstrated using several examples of microstrip components. These include two miniaturized couplers and a dual-band power divider. The verification results are supplemented by comparisons with multiple-start local optimization (to validate the need for global search) and a state-of-the-art population-based metaheuristic algorithm (to corroborate the efficacy of the presented approach).
The geometries of the test structures are introduced in "Case studies" along with the formulations of the respective design problems. "Setup and results" showcases the numerical results obtained using the proposed and the benchmark methods. "Discussion" provides a summary and discussions thereof.
Case studies. Numerical verification of the presented optimization procedure has been carried out using three microstrip circuits shown in Fig. 6. We have: • Circuit I: a compact microstrip rat-race coupler (RRC1) 67 , implemented on RO4003 substrate (ε r = 3.38, h = 0.762 mm). The adjustable geometry parameters are x = [l 1 l 2 l 3 d w w 1 ] T ; the remaining parameters are d 1 = d +|w -w 1 |, d = 1.0, w 0 = 1.7, and l 0 = 15 fixed (dimensions in mm). For this circuit, the goal is to minimize the matching and isolation characteristics, |S 11 | and |S 41 |, at the intended operating frequency f 0 , while ensuring a required power split ratio K P . The objective function is defined as in (2) ("Formulation of EM-driven design task"). • Circuit II: a compact rat-race coupler (RRC2) using a defected microstrip structure (meander spurline) within a folded transmission line, implemented on 0.15-mm-thick substrate 68 . The adjustable parameters are x = [L 1 b r g h fr s l fr ] T . The dimensions are in mm except for the relative quantities denoted using the subscript r; these parameters are unitless. The following relationships hold: , and h f = s + (w 0 -s)h fr ; dW = dL = 10 mm. The input line width w 0 is computed for a given substrate permittivity ε r so as to ensure 50 Ω input impedance. The design goal is-for a given the substrate permittivity ε r -to minimize the matching and isolation characteristics, |S 11 | and |S 41 |, at the intended operating frequency f 0 , while ensuring equal power split. The objective function is defined as in (2) but with K P = 0 dB. • Circuit III: a dual-band equal-split power divider (PD) 69 , implemented on AD250 substrate (ε r = 2.5, h = 0.81 mm). The adjustable geometry parameters are x = [l 1 l 2 l 3 l 4 l 5 s w 2 ] T (dimensions in mm); w 1 = 2.2 mm and g = 1 mm are fixed. The design goal is to simultaneously minimize the input matching |S 11 |, output matching |S 22 |, |S 33 |, as well as isolation |S 23 | simultaneously at two operating frequencies f 1 and f 2 . The objective function is similar to (3) but the equal power split condition is not directly handled in the optimization process as it is implied by the structure symmetry.
The computational models for all three circuits are implemented in CST Microwave Studio, and simulated using the time-domain solver. Table 2 provides information about the specific design tasks considered (i.e., the target values of the operating parameters), as well as the lower and upper bounds for parameters. It should be noted that parameter ranges are very wide with the average ratio of the upper-to-lower bound being 11.7, 4.6, and 10.3 for Circuit I, II, and III, respectively. The reason for selecting such broad ranges of geometry parameters was to ensure that the considered design tasks are sufficiently challenging, as well as to emulate the scenario under which the designer does not have a clear indication of what a good initial design should be, thereby to  Table 3 shows the experimental setup for the optimization framework proposed in this work and the benchmark methods. These include particle swarm optimization, employed as a representative population-based metaheuristic, as well as the trust-region gradient-based algorithm (as described in "Local optimization procedure") with random initial designs. The TR procedure is considered to demonstrate that local search is insufficient for the considered design task and often fails when the initial design is away from the target. The numerical results are provided in Tables 4, 5 and 6 for Circuit I, II, and III respectively. The scattering parameter responses at the designs obtained using the proposed framework for selected algorithm runs can be found in Figs. 7, 8 and 9.

Setup and results.
The breakdown of the computational cost of the proposed algorithm is the following (averaging over ten algorithm runs):  Table 2. Target operating parameters and parameter spaces for circuits I through III. a The circuit is to be optimized for a specific substrate of given relative permittivity ε r .  Table 3. Experimental setup: Proposed optimization framework and the benchmark.

Control parameters Termination condition Comments
Inverse

Discussion
The results gathered in "Setup and results" allow us to draw several conclusions concerning the presented optimization strategy, not only in terms of its efficacy and the computational complexity but also how its performance compares to the benchmark methods. These are the main points: • The presented approach does exhibit a global search capability, which is corroborated by the fact that satisfactory designs have been found in all runs of the algorithm (ten per case). Also, as illustrated in Figs. 7, 8 and 9, the parameter vectors x(0) found by the global search stage are of good quality (in particular, with the operating frequency being close to the target), which makes local optimization sufficient. On the contrary, the benchmark local search from random initial designs often fails (at about fifty percent of the cases), and its performance is highly dependent on the initial design quality. Consequently, the average value of the optimized objective function is considerably worse than for the proposed method. Population-based metaheuristics (here, PSO) performs much better, but its computational complexity is high. It can also be observed that there is a noticeable difference in the design quality produced after 50 and 100 PSO iterations. This indicates that computational budget set at 500 EM simulations is insufficient for this method. • The quality of the designs obtained using the presented technique surpasses that obtained using both local search and PSO. Although local optimization may produce a design that is of similar quality, it only happens when the initial points was sufficiently close to the target. Clearly, for the algorithm presented in this work, the problem of initial design has been eliminated, which makes it considerably more robust. • In terms of computational efficiency, the presented approach compares favorably with the population-based algorithm as the average running cost is only about 128, 88, and 99 EM analyses of the system for Circuit I, II, and III, respectively. At the same time, our methodology is only slightly more expensive than local optimization (by 23 percent on the average across all three circuits considered). • This is because the cost of the first (global) search stage is low, only 40, 22, and 31 EM simulations on the average for Circuit I, II, and III, respectively. This level of efficiency is a result of combining the response feature technology with inverse modeling, in particular, establishing the model over low-dimensional oper-   www.nature.com/scientificreports/ ating condition space. The latter requires a limited number of samples to identify relationships between the operating frequency, power split, etc., in a reliable manner.
As demonstrated, the presented technique offers both global search capability and computational efficiency. Both make it a low-cost alternative to mainstream global optimization methods, especially nature-inspired algorithms. This is particularly the case if the parameter space for the problem is set up in a reasonable way (e.g., not excessively large), and the likelihood that electrical characteristics of a randomly-generated design are not overly distorted is not excessively low.
A practical limitation (or, inconvenience) is that the feature-based approximations of the operating conditions need to be extracted from the EM simulated circuit responses, which is normally realized on case-to-case basis, although the respective implementations are transferrable within the same type of circuit responses (coupler, power divider, etc.). Automation of this process will be addressed elsewhere.

Conclusion
In this paper, a simple and reliable procedure for computationally-efficient globalized design optimization of passive microwave circuit has been presented. The fundamental component of the algorithm is an inverse regression model constructed using information extracted from a pre-selected subset of randomly generated parameter vectors and the corresponding EM-simulated circuit characteristics. Inverse model rendition involves the response feature technology, which exploits weakly nonlinear relationships between the geometry and operating parameters of the system at hand. Numerical verification of the proposed procedure has been carried out using three microstrip components, including two miniaturized rat-race couplers, and a dual-band power divider. In each case, ten independent algorithm runs were executed to validate the efficacy of the method as well as repeatability of solutions. Satisfactory designs have been found in all executions, which was not the case for the benchmark algorithms, especially multiple-start local search, where inferior designs were produced for a considerable number of optimization runs. At the same time, the computational cost of our framework is significantly lower than that of state-of-the-art global optimizers (here, PSO). As a matter of fact, it is comparable to the cost of local gradient-based optimization. The methodology discussed in this work might be an attractive alternative to conventional global search methods, particularly nature-inspired algorithms, but also hybrid methods incorporating forward surrogate modelling methods. The major advantages include low computational complexity, global search capability, as well as the improved immunity to dimensionality and parameter range issues. Although demonstrated for microwave passives, the proposed algorithm might be generalized to other type of circuits, including amplifiers or mixers; however, this requires additional investigation. While extending the applicability of the presented approach, one should keep in mind that in the case of other microwave components, the major limitation factor might be that the assumption concerning weakly-nonlinear relation between the figures of interest and designable parameters may no longer hold. In the latter case, utilization of the inverse surrogates (at the global search stage) would not be as efficient as demonstrated for the passive circuits. This will be considered in the future work. On the other hand, the non-uniqueness issues may be detrimental to the algorithm performance as well.