Direct constraint control for EM-based miniaturization of microwave passives

Handling constraints imposed on physical dimensions of microwave circuits has become an important design consideration over the recent years. It is primarily fostered by the needs of emerging application areas such as 5G mobile communications, internet of things, or wearable/implantable devices. The size of conventional passive components is determined by the guided wavelength, and its reduction requires topological modifications, e.g., transmission line folding, or utilization of compact cells capitalizing on the slow-wave phenomenon. The resulting miniaturized structures are geometrically complex and typically exhibit strong cross coupling effects, which cannot be adequately accounted for by analytical or equivalent network models. Consequently, electromagnetic (EM)-driven parameter tuning is necessary, which is computationally expensive. When the primary objective is size reduction, the optimization task becomes far more challenging due to the presence of constraints related to electrical performance figures (bandwidth, power split ratio, etc.), which are all costly to evaluate. A popular solution approach is to utilize penalty functions. Therein, possible violations of constraints degrade the primary objective, thereby enforcing their satisfaction. Yet, the appropriate setup of penalty coefficients is a non-trivial problem by itself, and is often associated to extra computational expenses. In this work, we propose an explicit approach to constraint handling, which is combined with the trust-region gradient-search procedure. In our technique, the decision about the adjustment of the search radius is determined based on the reliability of rendering the feasible region boundary by linear approximation models of the constraints. Comprehensive numerical experiments conducted using three miniaturized coupler structures demonstrate superiority of the presented method over the penalty function paradigm. Apart from the efficacy, its appealing features include algorithmic simplicity, and no need for tailoring the procedure for a particular circuit to be optimized.


Miniaturization of microwave passives with direct constraint control
Here, we introduce the procedure for miniaturization of microwave components proposed in this work. "Simulation-based size reduction: problem formulation" provides the formulation of the miniaturization task. In "Explicit constraint handling: the concept", we discuss the concept of direct control of CPU-heavy constraints within trust-region gradient-based algorithm. The technical details of controlling the tolerance levels for constraint violation as well as decision-making process that adjusts search radius are elaborated on in "Explicit constraint handling: constraint-related gain ratios". Finally, "Explicit constraint handling: optimization algorithm" summarizes the entire procedure.
Simulation-based size reduction: problem formulation. We denote by x = [x 1 ⋯ x n ] T a vector of adjustable parameters of the circuit under design. For passive components, these are normally geometry parameters (circuit dimensions). Let A(x) be the circuit size, e.g., its footprint area. The objective is to reduce the size as much as possible, i.e., to solve where x * is the optimum parameter vector to be identified, whereas U(x) is the objective function. In the case of miniaturization, we have U(x) = A(x). The problem (Eq. 1) is subject to constraint, which can be of or equality type In this work, we will only consider inequality constraints, which are the most common. Also, an equality constraint η k (x) = 0 can be represented in an inequality form as |η k (x)|≤ 0.
Let us consider an example of a microwave coupler, which is to be miniaturized while satisfying the following conditions.
Thus, for the exemplary coupler, we may formulate the miniaturization task as subject to constraints For the sake of illustration, let us consider another example of an optimization task, which is but oriented towards improving selected electrical performance figures rather than size reduction. Assume that the goal is to minimize the maximum reflection within the frequency range of interest of an impedance matching transformer. In this case, the merit function is defined as In Eq. (7), |S 11 (x,f)| stands for the circuit reflection, whereas F refers to the frequency range of interest.
Explicit constraint handling: the concept. Evaluation of constraints imposed on electrical characteristics of microwave components is computationally expensive: their values are obtained by post-processing EM simulation data. This is troublesome from the point of view of numerical optimization procedures, as local methods typically require constraint gradients. Unless adjoint sensitivities are available 49 , estimation of these requires finite differentiation, and the constraints may not be differentiable due to their very formulation as minimax functions (cf. "Simulation-based size reduction: problem formulation"). Also, EM simulation results may contain a certain level of numerical noise, being a result of adaptive meshing techniques, or specific termination criteria used by the EM solvers. As mentioned before, a common mitigation method is a penalty function approach 47 , where the cost function is defined through aggregation of the main objective (size reduction) and contributions from constraint violations, appropriately scaled using weighting factors (penalty coefficients). Although conceptually attractive, optimum setup of the coefficient values is generally an intricate task, often associated with preparatory optimization runs.
This paper offers an alternative approach to constraint handling, which is an explicit method. It employs linear approximation models of the system response, therefore, a natural choice for the underlying optimization algorithm is the trust-region (TR) framework 50 . The standard TR procedure yields a series of approximate solutions x (i) , i = 0, 1, …, that converge to x * . The new vector x (i+1) is obtained by solving where U L (i) is a first-order Taylor model of the scalar merit function U. The solution to Eq. (8) is restricted to the vicinity of x (i) determined by the size parameter δ (i) . Additionally, we have inequality constraints of the form of Eq. (2); equality constraints are not considered for the sake of simplicity, cf. "Simulation-based size reduction: problem formulation". The TR radius δ (i) is adaptively adjusted using the standard TR rules (e.g. 50 ).
The problem (Eq. 8) is solved using Matlab's fmincon algorithm, which implements the Sequential Quadratic Programming (SQP) procedure, one of the state-of-the-art procedures for constrained continuous optimization. The SQP procedure directly handles geometry constraints (here, the TR condition ||x − x (i) || ≤ δ (i) ), whereas the constraints related to electrical performance figures are controlled using the explicit approach being the subject of this paper. The details are explained in the remaining part of this section.
If size reduction is of interest, the evaluation of the merit function incurs negligible costs: the structure size can be obtained directly from the system geometry parameters, i.e., the vector x. On the other hand, maintaining solution feasibility becomes problematic due to expensive constraints. In this work, constraint control is achieved by the incorporation of linearized models γ L.k of the constraints γ k , k = 1, …, n γ , and adaptive adjustment of the trust-region size parameter δ (i) . The decision-making process governing the latter involves quantification of the reliability of γ L.k in predicting the feasibility status in the course of the optimization process.
In the following, we will denote as r(x) the vector of EM-simulated outputs (e.g., S-parameters) of the circuit of interest. A first-order Taylor model r L (i) (x) of the response r(x), established at the design x (i) , is defined where the Jacobian matrix of the component response at the design x (i) is denoted as J(x (i) ). In most cases, it is estimated by means of finite differentiation. For the purpose of subsequent considerations, we will explicitly indicate that the constraints are functions of the circuit characteristics, i.e., we have γ k = γ k (r(x)). Then, the linear model γ L.k is defined as www.nature.com/scientificreports/ When solving the trust-region sub-problem (Eq. 2), the exact constraints γ k (r(x)) will be replaced by their linearized versions (Eq. 10). The accuracy of representing γ k by γ L.k depends on a particular location in the parameter space and on the trust-region size parameter δ (i) , which is because ||γ k (r(x)) − γ L.k (x)|| is proportional to ||x − x (i) || 2 (for sufficiently small design relocations). Consequently, a proper updating procedure for δ (i) is essential. In particular, maintaining small values of the TR radius improves the alignment between γ L.k (x) and γ k (r(x)), whereas increasing it allows for increased-size steps in the design space while solving (Eq. 8). The adjustment of δ (i) should take into account the solution feasibility predictions according to γ L.k (x), but also the actual feasibility status (as verified by EM simulation). At the generic level, the adaptation scheme is arranged the same way as for conventional TR algorithms 50 , i.e., In Eq. (11), the coefficients m inc and m dec are used for incrementing and decrementing the TR region size, respectively, whereas θ inc and θ dec represent appropriate threshold. In our approach, we employ their typical values, i.e., we have m inc = 2 and m inc = 3, as well as θ inc = 0.75 and θ dec = 0.25 50 . As mentioned earlier, these are the typical values used in the TR algorithms. According to classical theory (e.g. 50 ), the specific values of coefficients are not critical for the algorithm operation.
However, the decision-making factor θ in Eq. (11) is the gain ratio pertinent to the constraints. It compares the actual alteration of constraint violations to those predicted by the linear model for subsequent iterations (specifically, the (i + 1)th versus the ith one). The modification coefficients as well as the thresholds in Eq. (11) mimic the conventional rules of the TR frameworks (cf. 51 ).
In the case of multiple constraints, the coefficient ρ is generalized to account for the worst-case situation over the entire set g k , k = 1, …, n g . We have Definition of the factors θ k for each γ k (x), k = 1, …, n γ , is pivotal to the successful operation of the proposed optimization procedure. It will be discussed at length in the next section.
Explicit constraint handling: constraint-related gain ratios. This section elaborates on the definition and evaluation of the constraint-based gain ratios θ k , utilized to control the trust region size as discussed in "Explicit constraint handling: the concept". In the following, we will denote by Γ k (i) the threshold for accepting the violation of γ k (x) at the iteration i. The threshold is iteration dependent for the reasons explained at the end of the section. At this point, we will outline the rules for computing the ratios θ k utilized in decision-making process that governs the search radius adjustments: • Rule 1: If γ k (r(x (i) )) > Γ k (i) , i.e., the constraint violation before executing the (i + 1)th iteration exceeds the acceptance threshold, then i.e., the constraint violation is at the acceptable level, then • Rule 3: If θ k < 0 (as computed using Eqs. (13) or (14)) but γ k (r(x (i+1) )) ≤ Γ k (i) , i.e., EM-evaluated constraint violation is acceptable, then the value of θ k is overwritten to θ k = 0.5.
The above rules serve for two purposes. On the one hand, one needs to impose penalty on insufficient prediction accuracy of γ L.k if the constraint violation is large (Rule 1), or the feasibility condition has not been improved in the case of minor constraint infringement (Rule 2). On the other hand, the conditions (Eqs. 13 and 14) are employed to promote sufficient prediction of γ k by γ L.k (cf. Eq. (13)), or relocation of the design towards the feasible region (cf. Eq. (14)). The role of Rule 3 is to overwrite the previous ones if the EM-evaluated constraint violation γ k (r(x (i+1) )) at the candidate design x (i+1) is within the acceptance threshold. Rule 3 has been introduced to prevent erratic operation for the designs residing in the vicinity of the feasible region boundary, in particular, near-zero constraint infringements (either positive or negative) in any given iteration. The graphical illustration of acceptable and insufficient evaluation of the constraint γ k by the linearized model γ L.k has been provided in Fig. 1.
In the remaining part of this section we discuss the acceptance thresholds Γ k (i) . At the early stages of the optimization process (far from convergence), it is advantageous to relax the acceptance thresholds for constraint violation in order to facilitate identification of small-size designs. However, when close to convergence, the thresholds should be tightened to ensure more precise control over constraints. In practice, this is realized by adjusting the threshold values based on the convergence status of the optimization process. Let Γ k.max be a user-defined maximum violation level. Further, let ε be a small positive number determining the algorithm www.nature.com/scientificreports/ termination. In this work, the termination condition is formulated as ||x (i+1) − x (i) ||< ε or δ (i) < ε. For a given iteration i, let us define the convergence indicator Note that ψ (i) is large when the optimization process is launched, and it is reduced to unity upon convergence. We also have the threshold , it is initially equal to Γ k.max , and gradually diminished to αΓ k.max upon convergence. Here, α is assumes small values greater than zero, e.g., 0.1 or 0.01. This value is not of key importance for the operation of the procedure.
Explicit constraint handling: optimization algorithm. Figure 2 shows the operating flowchart of the proposed size reduction procedure with explicit handling of design constraints. Apart from the termination threshold, the algorithm only contains the following control parameters: the threshold Γ k.max (maximum tolerance for constraint violation), and the scaling coefficient α. These parameters are not critical for the algorithm performance. As a matter of fact, we will keep these values fixed for all verification case studies considered in "Demonstration examples". The acceptance of the candidate design x (i+1) produced by solving (Eq. 8) is based on the standard TR rules, i.e., it is accepted if the decision-making factor θ is positive, and rejected otherwise. In the latter case, the iteration is repeated with a reduced search region. Note that θ > 0 if either the violation of the constraint has been reduced to a sufficient extent, or the design was relocated to the feasible region.

Infeasible region
Feasible region The first constraint enforces equal power split ratio within 0.1 dB tolerance, whereas the second ensures that the circuit impedance matching and port isolation are better or equal − 20 dB within the operating band. Table 1 provides essential data about all three structures.
Experimental setup and results. The proposed optimization procedure has been applied to the circuits of Fig. 3. In each case, the initial design (the last row of Table 1) was obtained by optimizing the respective circuits to improve the matching and isolation within the operating bandwidth F, subject to equal power split constraint. This means, in particular, that the starting points are feasible from the point of view of both constraints γ 1 and γ 2 (cf. "Verification circuits"). The termination threshold is set to ε = 10 −3 , the acceptance threshold are chosen to be Γ 1.max = 1 dB, Γ 2.max = 0.3 dB, and α = 0.1.

Initial design x (0)
Evaluate response r(x (i) ) and gradients J(x (i) ) Set iteration index i = 0 Prepare linear models γ L.k , k = 1,...,n γ Constraints γ k , k = 1,...,n γ Compute convergence factor ψ (i) and thresholds Γ k The benchmark algorithm is run for all combinations of the penalty coefficients β 1 ∈ {10, 100, 1000, 10,000}, β 2 ∈ {10, 100, 1000, 10,000}. This is to illustrate the fact that optimum performance of the algorithm requires identification of the appropriate setup of the penalty terms, and sub-optimal setup leads to inferior constraint control or miniaturization rates.
The numerical results have been gathered in Tables 2, 3 and 4. These include the achieved footprint area of the respective circuits, as well as constraint violations at the final design. Figures 4, 5 and 6 show the circuit characteristics at the initial and the optimized designs for all circuits, along with the history of the circuit footprint area and violation of constraints during the optimization run.
Discussion. The results presented in "Experimental setup and results" allow us to formulate a number of remarks concerning performance of the proposed optimization procedure with explicit handling of design constraints. Furthermore, our methodology can be conclusively compared with the benchmark methods employing the penalty function approach. The observations are as follows.  Design parameters Other parameters L = 2dL + L s , L s = 4w 1 + 4 g + s + l a + l b , W = 2dL + W s , W s = 4w 1 + 4 g + s + 2w a , l 1 = l b l 1r , w 2 = w a w 2r , w 3 = w 3r w a , and w 4 = w 4r w a , w c = 1.9 mm www.nature.com/scientificreports/ • The proposed algorithm performs consistently for all considered verification circuits. In particular, it enables a satisfactory control of both design constraints. There is no violation for power split ratio constraint observed, whereas the maximum violation for the matching/isolation constraint is only 0.5 dB for Circuit II (design scenario II); however, it should be noted that the acceptance threshold is − 20 dB. • The performance of benchmark methods is highly dependent on the penalty coefficient setup. Among sixteen combination of parameters, only a few lead to satisfactory results in terms of both ensuring good miniaturization ratio and sufficiently precise constraint handling. For Circuit I we have about three of such 'good' setups, for Circuit II there is only one (per design scenario), whereas for Circuit III about three.   Table 3. Optimization results for Circuit II. www.nature.com/scientificreports/ • For the particular setups ensuring good performance of the benchmark procedure, the obtained circuit sizes are comparable to those obtained by the proposed algorithm (which, on the other hand, does not require any setup or tailoring to the task at hand). • In the case of Circuit II, for most combinations of penalty coefficients featuring β 1 ≥ 10 2 , the optimization process becomes stuck at the early stages of the optimization process, leaving large feasibility margin for the second constraint. This is due to the fact that a large value of the first penalty coefficient along with a small margin for the power split ratio constraint (only 0.1 dB), makes the problem numerically challenging. More specifically, the objective function becomes highly nonlinear near the feasible region boundary, which hinders exploration of that region and leads to a premature convergence of the process. A similar effect can be observed for Circuit III, although it is less pronounced. • The proposed algorithm turns out to be less prone to the aforementioned issues due to the adaptive adjustment of the acceptance thresholds governed by the convergence status of the algorithm (cf. Eqs. (15,16)). • The average costs of rendering the optimal designs by the proposed approach equal: 110, 80 and 55 full-wave simulations, for Circuit I, II and III, respectively. Whereas in the case of the benchmark procedure the corresponding costs are: 115, 57 and 45 EM simulations. Therefore, the proposed procedure is around 20% more expensive in terms of the number of EM analyses necessary for the algorithm to converge. Yet, given that in our approach there is virtually no need for tailoring the algorithm to render a satisfactory design meeting the design specifications, this additional cost seems to be justifiable. This is because tuning the optimization procedure to ensure its satisfactory operation normally entails additional computational expenses (e.g., for adjusting penalty coefficients in the case of implicit methods). Furthermore, the primary purpose of the presented technique was to improve the precision of controlling design constraints, and miniaturization rate, both of which have been conclusively demonstrated.
Given a large combined number of circuits, design scenarios, as well as penalty coefficient setups involved in this verification study, the observations summarized above should be categorized as conclusive. Overall, the performance of the presented procedure can be considered competitive over the benchmark (implicit) methods, both with respect to the accuracy of constraint handling and achievable miniaturization rates. The important advantages of the proposed algorithm include easy implementation and no need for adjusting any control parameters. The latter normally incurs extra computational expenses and may require a certain level of experience pertinent to optimization methods.

Conclusion
The purpose of this work was to propose a novel procedure for simulation-based miniaturization of microwave passives. Our approach involves direct control of constraints imposed on electrical performance figures of the circuit under design. Linear approximation models of the constraint functions are employed to make predictions concerning solution feasibility. Appropriately quantified quality of these predictions is utilized in the decisionmaking process that controls the search region size within the trust region framework. Furthermore, the constraint violation tolerance thresholds are governed by the convergence indicators of the optimization process in order to foster more aggressive size reduction at the early stages of the optimization process. Comprehensive numerical verification involving three microstrip couplers and six design scenarios demonstrate superior performance of the proposed technique as compared to the benchmark methods employing a penalty function approach for implicit constraint handling. Its major advantages include competitive size reduction ratios, accuracy in controlling constraint violation levels, consistency of results obtained for a variety of problems, straightforward implementation, as well as no need for tailoring the procedure to handle a particular microwave structure. The last feature is particularly important in practical applications: tuning the optimization procedure to ensure satisfactory operation (e.g., setting up penalty coefficients for implicit methods) normally entails additional computational expenses and may require optimization-related knowhow lacking by many microwave engineers. www.nature.com/scientificreports/ www.nature.com/scientificreports/

Data availability
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request. Contact person: anna.dabrowska@pg.edu.pl.