Optimization of the Hardware Layer for IoT Systems Using a Trust-Region Method With Adaptive Forward Finite Differences

Trust-region (TR) algorithms represent a popular class of local optimization methods. Owing to straightforward setup and low computational cost, TR routines based on linear models determined using forward finite differences (FDs) are often utilized for performance tuning of microwave and antenna components incorporated within the Internet of Things (IoT) systems. Despite usefulness for design of complex structures, the performance of TR methods vastly depends on the quality of FD-based local models. The latter are normally identified from perturbations determined a priori using either rules-of-thumb, or as a result of manual tuning. In this work, a framework for the automatic determination of FD steps and their adjustment between the TR algorithm iterations is proposed. The method involves numerical optimization of perturbations so as to equalize the objective function changes w.r.t. the center design to the desirable precision. To maintain acceptable cost, the FD-tuning procedure is executed using the same approximation model as the one exploited in the course of the structure optimization. The proposed framework has been tested on a total of 12 design problems. Furthermore, the presented method has been thoroughly validated against TR algorithms with static, a priori selected perturbations. Numerical results indicate that the proposed framework provides up to 50% performance improvement (in terms of the optimized designs quality) compared to the state-of-the-art TR-based approaches. The usefulness of the proposed method for the real-world IoT systems has been implicitly demonstrated through the utilization of one of the optimized structures in a hardware layer of a real-time localization system.


I. INTRODUCTION
N UMERICAL optimization belongs to the most important tools in contemporary engineering. Modern problems, including the development of wireless hardware layer for Internet of Things (IoT) systems, are often characterized by many complex requirements that can be fulfilled only using intricate design solutions [1], [2], [3], [4]. Conventional approaches to circuits design-based on manual tuning of geometry intertwined with inspection of the performance characteristics-are impractical due to a combination of factors that include a large number of parameters, challenges related to maintaining reliable control over multiple objectives at a time, and high cost of electromagnetic (EM) simulations required for performance evaluation [3], [4]. From this perspective, neglecting the engineering insight in favor of unattended development of microwave components seems to be an interesting alternative [5], [6], [7], [8]. However, automatic design methods are largely based on population-based optimization algorithms coupled with EM solvers which make them infeasible for the development of complex structures [7], [8], [9], [10]. Different, and widely accepted, design approach involves synthesis and semi-manual adjustment of the structure geometry (including experience-driven topology changes) followed by its numerical optimization [3], [11], [12], [13]. From the design optimization perspective, synthesis narrows down the search space to a region of interest and determines a starting point for further tuning of the structure performance. Gradient-based methods are popular algorithms used for the development of microwave components in such a framework [9], [14], [15], [16]. Unfortunately, they often require dozens to hundreds of simulations to converge [9], [11], [17]. Having in mind that accurate performance evaluation of modern structures can be performed only using the expensive high-fidelity EM simulations, the application of conventional gradient-based algorithms to design optimization is numerically impractical.
The challenge associated with the high optimization cost of microwave/antenna structures EM models can be mitigated using surrogate-assisted methods [18], [19], [20], [21], [22], [23], [24], [25]. The latter ones shift the optimization burden from the high-fidelity EM model to a cheap surrogate constructed either from the corrected low-fidelity EM simulations [18], [19], [20], [21], or a response surface approximation model [22], [23], [24], [25]. The surrogate-based optimization (SBO) is performed in a loop where the desired solution is approximated using the auxiliary model and then verified and/or corrected based on the data obtained from scarce highfidelity EM simulations. Trust-region (TR) methods represent the most straightforward realizations of the SBO concept [11], [16], [17], [19]. Owing to low cost and simple implementation, TR routines find application in many research areas that include the development of microwave and antenna components, but also photonics and aerospace engineering [19]. A popular variant of the TR-based design exploits a firstorder Taylor-expansion model, constructed from the derivative data, for structure optimization [11], [16], [17], [26]. When the component is evaluated using adjoint-capable EM solvers, the derivative information required for the construction of the model can be obtained at a low computational overhead compared to the zero-order response [26]. Unfortunately, adjoint technology is supported only by a handful of commercial solvers [27], [28]. Alternatively, a finite-difference (FD) method can be used to obtain the derivatives. The latter are calculated based on the evaluation of the structure responses at a certain design of interest and small perturbations around it [16], [29], [30]. The practical challenge related to FD boils down to the selection of appropriate steps around the given design [30]. The perturbations determine the quality of the linear model and may substantially affect the optimization performance, both in terms of the obtained solution and the number of iterations required by the algorithm to converge. On the one hand, the steps should be large enough so that the numerical noise (i.e., slight performance inconsistency resulting from discretization imperfections of the EM model at closely located designs) does not affect the descent direction estimated using the derivative data. On the other hand, the steps have to be small so as to provide decent representation of the functional landscape in the vicinity of the given design. It is worth noting that the methods dedicated to mitigate the effects of numerical noise have been proposed in [31] and [32]. However, they are not widespread available in commercial EM solvers.
The problem related to selecting appropriate perturbations for FD-based models in TR optimization is often deemed as a "technical detail" that does not receive suitable attention in [17], [19], [33], and [34]. For instance, the discussion related to the determination of derivatives often boils down to mentioning that the first-order data have been determined using FD [33], or large-step FD (where large step refers to perturbations being large enough that they are not affected by the numerical noise) [34]. Similarly, in [17], the details concerning perturbations or their effects on algorithm performance are not provided. Although the problem of selecting suitable FD steps for TR-based methods seems not to receive the sufficient attention in the literature, the effects of FD on gradient-based methods (being the core of TR approaches discussed herein) performance are the subject of research [29], [30], [35], [36]. The generation of appropriate FD perturbations, often referred to as a step-size dilemma [30], involves balancing between errors caused by too small and too large steps [30], [35].
The FD perturbations are either setup as static, or adaptive [35]. The first category comprises the steps selected a priori (i.e., before the optimization process) as suitable fractions of the initial design, or the machine-precision of floatingpoint data [30]. Despite being useful for optimization [17], [34], [36], a priori defined FD might also produce inaccurate approximations of the structure response changes resulting in the increased number of iterations required by the algorithm to find a solution and/or cause its premature convergence [30], [35]. Adaptive FD methods incorporate a feedback mechanism that permits assessing the effects of perturbation size on change of the structure response [30]. Unfortunately, many of the adaptive methods require multiple model evaluations per parameter in order to determine appropriate step size [30], [37], [38]. Having in mind a large number of parameters used for the representation of modern wireless communication structures, as well as high-cost associated with their EM simulations, utilization of conventional adaptive methods for step-size estimation is impractical. Furthermore, although a number of advanced techniques oriented toward the identification of accurate derivatives has been proposed [30], they are either problem specific, expensive in terms of the number of model evaluations, or require access to the source code of the solver. Consequently, adaptive FD methods available in the literature are of little to no use for the design problems discussed here.
The above considerations indicate that, despite being the key tool in the design of modern wireless communication components, conventional numerical optimization incurs substantial (often prohibitive) computational cost. This challenge can be mitigated using SBO routines such as TR algorithms [16]. The core component of the TR framework is a local approximation model, which (due to low cost) is often based on FD [11], [33]. The problem of selecting appropriate FD perturbations does not receive sufficient attention in the context of engineering design using numerically efficient algorithms, whereas the gradient-based optimization either exploits the a priori (widely adopted to TR methods), or adaptively obtained steps. While the former cannot ensure that the model will represent the functional landscape with sufficient accuracy, the high cost of automated methods makes them unsuitable for the optimization of structures represented using high-fidelity EM models. In this context, the problem related to rapid estimation and adjustment of FD steps for TR-based optimization of modern wireless communication devices remains unsolved.
In this work, a TR framework that exploits linear approximation models with adaptively adjusted FD has been proposed. The method involves iterative tuning of the Taylor-expansion surrogate to yield a set of perturbed designs characterized by elevated changes of the objective function values w.r.t. the center design. The low cost of the procedure is ensued by the reuse of the Taylor surrogate for structure optimization and stepsize tuning. The quality of perturbations is validated through a comparison of the responses predicted by the model and obtained from EM simulations. The accepted designs are used to construct the surrogate for the consecutive TR iteration, whereas the rejected ones undergo reoptimization and further EM-based evaluations until the quality of FD steps is acceptable. The presented method has been validated using three test designs concerning a radiator for industrial scientific and medical (ISM) applications, an ultra-wideband (UWB) antenna, and a broadband rectifier for radio-frequency energy harvesting. The numerical results have been obtained for a total of 12 design problems. The proposed framework has been compared against the TR methods with a priori defined FD perturbations. The numerical results indicate that the presented approach is characterized by up to 50% higher performance compared to the state-of-the-art TR methods. Furthermore, the optimized UWB structure has been utilized in the hardware layer of a real-time localization system (RTLS) in order to implicitly demonstrate the usefulness of the proposed framework for real-world IoT applications.
The remainder of the work is organized as follows. Section II is devoted to the formulation of the design problem and discussion of the TR framework with emphasis on forward FDs. Section III focuses on formulation and explanation of the proposed mechanism for adaptive FD steps determination, as well as the discussion of the modified TR routine. Numerical results are presented in Section IV. A demonstration of the UWB radiator as a part of the UWB-RTLS is discussed in Section V, whereas Section VI summarizes this article and suggests possible directions for future work. Table I contains a glossary of the most important abbreviations and symbols used in this work.

II. TRUST-REGION FRAMEWORK
The TR algorithm belongs to the class of SBO methods where the optimization burden is shifted from the high-fidelity EM model to an auxiliary surrogate. To make the paper selfcontained, a brief discussion of the TR framework is provided here along with an in-depth explanation of the FD role in the EM-driven design. The adaptive-step FD and its integration with TR are presented in Section III.

A. Design-Optimization Problem
Let R(x) denote the EM-based response of the structure under design obtained for the vector of input parameters x. The optimization problem can be formulated as the nonlinear task of the form where X represents a feasible region of the search space (determined by the lower and upper bounds l b and u b ), x * denotes the optimal design to be found, and U(x) = U(R(x)) is a scalar objective function. Due to a high cost of EM simulations and a large number of evaluations required to find the optimum solution, direct solving of (1) is impractical (especially when the design of complex, multiparameter structures is considered) [19], [39]. Instead, the surrogate-assisted design can be performed to approximate (1) at an acceptable computational cost.

B. Trust-Region Framework
Let y be the representation of x normalized to a unit-size hypercube Y = [0, 1] D (D is the search space dimensionality) and P(y) be the structure response [22]. Note that P(y) = R(x). Normalization provides an equal representation of the parameters regardless of their ranges in X, which yields a clearer description of the discussed methodology.
The normalized-input Taylor model G (i) is defined as where J denotes the Jacobian matrix of the first-order partial derivatives [30]. Preferably, the latter is obtained using EM simulations with adjoint sensitivity which ensure high accuracy and low numerical overhead as compared to the evaluation of the zero-order response [26]. In practice, the Jacobian is often approximated using the forward FD method [16], [30]. For more comprehensive discussion on the TR framework, see [16], [17], [26].

C. Jacobian Approximation Using Finite Differences
The Jacobian, which is a crucial component of the model (4), has the form where P d (y (i) ) represents the first-order derivative of the design y (i) about dth (d = 1, 2, . . . , D) parameter. As already mentioned, numerical derivatives are not widespread available in commercial EM solvers. Instead, first-order responses can be approximated using forward FD as [30] P d (1) The parameter h d represents perturbation size w.r.t. the dth variable of y. The vector of perturbations and the standard basis in the dth dimension are h = [h 1 · · · h d · · · h D ] T and e d = [0 · · · 1 · · · 0] T , whereas h d = h • e d (here, " • " denotes component-wise multiplication). In practical FD realizations, the residual O(h d ) is often neglected [30]. Assuming the existence of analytical derivatives P d (y (i) ), then O(h d ) = P d (y (i) ) -P (1) d (y (i) ) which indicates that-with O(h d ) set to zero-FD only approximates the first-order data, i.e., The approximation accuracy-and hence the quality of the model (4)-depends on the composition of the perturbations vector h. In other words, the set of parameters h used in (6) implicitly affects the TR framework performance.
The Jacobian can also be estimated using the backward, central, or higher order FDs [38]. The residual of backward and forward methods is proportional to h d resulting in their similar performance [40]. The central FD is [30] P d (2) In (8), the residual is proportional to h 2 d which provides a better estimation of the analytical derivative for small perturbation steps [40]. On the other hand, central FD it is characterized by twofold higher cost compared to the forward one [30], [40]. Higher order methods are associated with even greater computational burden and thus are not considered here [38].
The problem concerning the determination of the appropriate perturbations-often referred to as a step-size-dilemma-is conceptually illustrated in Fig. 1. It boils down to balancing the FD size so as to minimize both the truncation (too large step) and round-off (too small perturbation) errors [30]. Although the evaluation of the structure for a range of steps with varying sizes can provide a general understanding of their close-to-optimal values, the computational cost of the process is often too high for the EM-driven design [30], [37], [38]. Another important challenge pertinent to the optimization of contemporary structures involves the numerical noise associated with the EM simulations [31], [32]. Conceptual illustration of its effects on the FD-based model quality-shown in Fig. 2-indicates that the step size in each dimension has to be large enough so as to accurately track the objective function changes regardless of its local "oscillations." At the same time, the steps have to be small so as to ensure that they capture local variations of the performance features.
Despite the significance of the mentioned problems for the TR-optimization performance, FD perturbations are the  most often determined based on the rule-of-thumb approaches [35], [39]. The main problems here include the lack of knowledge on appropriate step sizes, as well as changes of the optimum perturbations in the course of the optimization process (i.e., from one design to another). Furthermore, a priori selection of perturbations does not provide feedback on their effects on the response which might hinder the determination of a good design solution. As illustrated in Fig. 3, FD step sizes may considerably affect the performance of the TR-based design, both in terms of the solutions quality and numerical cost (expressed in the number of EM simulations for successful TR steps). Clearly, the FD set-up through execution of a few optimizations with different sets of a priori selected perturbations might produce a satisfactory solution, yet is numerically impractical.
A more-rigorous, feedback-enhanced estimation of perturbations may yield improved quality of the Jacobian, yet at the expense of increased cost compared to the rule-of-thumbbased FD steps [30], [35]. A rudimentary feedback method involves manual adjustment of perturbations intertwined with a visual inspection of their effects on changes of the structure response [21], [23], [33]. On the other hand, manual FD tuning is laborious, limited in accuracy, and cannot guarantee improved optimization performance compared to methods that exploit models based on conventional FD steps. More advanced, adaptive techniques can accurately estimate FD perturbations. However, they require a large number of EM model evaluations per design parameter to determine the optimal FD steps [30], [38]. Consequently, they are numerically impractical when the evaluation of the structure response involves EM simulations.
Finally, for the sake of low cost, the a priori FD perturbations are not updated in between the TR iterations. Instead, the presumed validity of the model (4) is controlled using the TR radius. Notwithstanding, modern structures are characterized by complex, nonlinear, and multimodal functional landscapes [11]. Therefore, the step sizes adjusted for the given design may incorrectly represent the response changes at another one. In extreme-yet not uncommon-situations, lack of control over step size leads to inconsistency between signs of the structure response changes obtained by the linear and EM models. As a consequence of violating this (mildest possible) restriction, the descent direction predicted by (4) counters the one obtained from the EM simulations (in the affected dimensions) [35]. The conceptual illustration of the problem is shown in Fig. 2, where the linear model identified based on the parameters distorted by the numerical noise features incorrect slope w.r.t. the objective function changes.

III. FINITE DIFFERENCES WITH ADAPTIVE STEP SIZE
In this section, a method for automated determination and in-optimization tuning of FD steps for EM-driven optimization using the TR algorithm is presented. Numerical validation and discussion of the approach are provided in Section IV.

A. Effect of Perturbations on Objective Function Change
As already explained, conventional FD mechanisms (cf. Section II-C) provide little-to-no means for analyzing the effects of perturbations on the quality of the Jacobian. The validation of the selected FD steps (if any) is typically performed by visualizing variations of the structure response compared to the given reference design. On the other hand, from the objective function standpoint, the relevance of such analysis might be limited. As conceptually illustrated in Fig. 4, seemingly noticeable change of the structure response might not be important in terms of the defined performance requirements. Conversely, visually small variation of the response may contribute to the evident change of the objective function value. Another important remark is that optimization algorithms are normally set-up to analyze and adjust the objective function (i.e., a representation of the selected response features) rather than the entire characteristic of the structure at hand. From this perspective, analysis of the FD step-size effects on the defined design criteria seems natural and straightforward approach for selecting appropriate perturbations of individual parameters. In other words, examination of the objective function changes allows maintaining systematic control over perturbation sizes in each dimension d while ensuring that all relevant components of the structure response (objective-function-wise) are accounted for.

B. Precision of Input Parameters and Structure Response
The precision of simulation results obtained from complex computational models is often significantly lower compared to the machine precision of contemporary computers. The reasons include the accumulation of errors in the course of calculating the structure response based on its discretized model [35], [40], but also the numerical noise [32]. The undesirable consequence for the TR optimization is that one needs to obtain a possibly accurate approximation of derivatives while having data responses characterized by a limited precision. Similarly, in the case of input parameters, the number of significant digits (i.e., the machine precision used to represent the floating-point data) processed by commercial solvers is much lower than the machine precision. In practice, many EM solvers support single-precision floatingpoint data [27], [28]. As already indicated in Section II-C, the data accuracy is normally lower than expected from the obtained number of significant digits [35]. Here, the problem Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.
is accounted for by assuming that the minimum acceptable perturbation size does not exceed where ε A = 10 −7 corresponds to the single-precision floatingpoint data. In fact, due to numerical noise, the accuracy of the responses obtained for h min might be insufficient for reliable FD. Nonetheless, the value represents the lower bound on the (nonnormalized-cf. Section II-B) step size permitted by the discussed algorithm.

C. Optimization of FD Perturbations
and u h = 0.1·1 T be the lower and upper bounds on the FD step size that defines the search space H. Here, s(·) is the design normalization function (defined w.r.t. l b and u b -see Section II-A), 1 is the D-dimensional vector of ones, and h min = h min ·1 T . Then, let y (i) be the center design for FD obtained from (2). Given the availability of the model (4), the vector of perturbations can be approximated as The minimization in (10) is performed using a MATLABbased fmincon function with the standard set of termination conditions [43]. The objective function is of the form where w d is the coefficient that activates the optimization of the dth perturbation (i.e., w d = 1 when h d is to be adjusted and w d = 0 for the optimized and EM-validated perturbations; cf. Section III-D) and k denotes the desired number of significant digits for the function [35] Here, U obj represents the design objective for the structure under design, (4). The minimization process (10) yields the set of FD steps h t that equalize the variations of the objective function responses obtained for the perturbed designs around y (i) to the selected number of significant digits k. Note that the use of log 10 (·) in (12) permits adjusting (11) so as to equalize the difference between y (i) and the perturbed designs (in terms of the objective function U obj ) according to the selected number of digits.
The values from the vector h t yield by solving (10) might reach the lower and/or upper bounds on the perturbations, indicating that the predicted objective function U obj changes are either too high or low. Here, the problem is addressed through embedding the step-size optimization problem into a heuristic routine that iteratively adjusts k to mitigate the effects of hitting the bounds on the optimization. The algorithm is below. l , set k = min{k + k stp ; k max } and go to step 7. 7) If j < j max , set j = j + 1 and go to step 2; otherwise, END. The presented routine counts the number of h t variables that attain the bounds and increase/decrease k accordingly to minimize the U obj for the perturbed designs. Here, n l and n u represent the number of parameters from the h t vector that reach the lower and upper bounds, respectively. The coefficient k stp is user-defined (here, k stp = 0.3); k max = 4 which roughly reflects half of the digits that represent the machine precision of the objective function responses and j max = 3. Fig. 5 shows a comparison of the objective function responses and perturbation sizes before and after the optimization. One should emphasize that the computational cost of the algorithm is negligible as the vector h t is adjusted using the Taylor-expansion model obtained in iteration i -1 of the TR algorithm.

D. Surrogate Model Reset Using Adaptive Perturbations
As already explained, the vector h t is obtained through optimization of the linear model G (i−1) , i.e., the one used for generation of the y (i) design. Consequently, h t is just an approximation of the optimal FD perturbations. For the sake of notational clarity, let G (i.l) denote the Jacobian used to generate the h (i.l) = h t vector in the ith TR algorithm iteration with l being the step count of the presented surrogate model reset procedure. Note that for l = 0 the model G (i.l) = G (i−1) . As discussed in Section II-C, the linear model is obtained using the forward FD method, where the set . . D of designs is evaluated by the EM solver and utilized to construct the Taylor-expansion candidate G (i.l+1) . Next, the responses of G (i.l+1) are validated against the ones obtained from the G (i.l) in terms of the signs. It should be noted that, due to the interpolative nature of the surrogate, the responses G (i.l+1) at the designs from H (i.l) set are the same as the ones determined from the EM simulations [i.e., G (i.l+1) (H (i.l) ) = P(H (i.l) )]. Validation of the predicted and obtained objective function responses in terms of the sign represents the mild acceptance threshold for the predicted perturbations and allows for the identification of parameters that might require further refinement. Let w = [w 1 · · · w d · · · w D ] T represent the vector of coefficients used to activate the optimization of the specific parameters from h (i.l) ; the dth element of w (d = 1, . . . , D) enables/disables individual components of the objective function (11) when set to 1/0 and is calculated as follows: The values of w obtained from (13) measure the quality of the objective function changes predicted by G (i.l) compared to G (i.l+1) in terms of the individual perturbations (cf. Section III-C). In other words, nonzero elements of w indicate that the predicted and the obtained objective function changes w.r.t. design parameters with their corresponding indices are opposite (e.g., one increases whereas the other decreases) and hence violate the acceptance threshold. In such a case, the vector h (i.l) is used as a starting point for reoptimization of the finite-difference perturbations using the G (i.l+1) model. Note that, in this step, only the parameters activated by the nonzero components of w are refined (as explained in Section III-C). Next, the obtained h (i.l+1) = h t vector is used to update the linear model. The procedure continues until ||w|| = 0, i.e., the predicted and obtained changes of the objective functions are the same otherwise, go to step 6. 6) Set l = l + 1 and go to step 2. As indicated above, the proposed model reset procedure embeds the routine of Section III-C. It should be reiterated that the first approximation of the perturbations (l = 0; w = 1 T ) is performed using the Jacobian calculated around the y (i−1) design. From this perspective, the model predictions might not provide an accurate representation (even sign-wise) of the actual objective function changes around y (i) . On the other hand, reuse of G (i−1) allows maintaining a relatively low cost of the presented FD steps adjustment method. Normally, the number of algorithm iterations l required to reset the surrogate is no greater than 3. Fig. 6 illustrates the predicted and obtained effects of the perturbations on the objective function changes before and after the refinement, as well as the obtained finite-difference step sizes.

E. TR Optimization With Adaptive Step Size Adjustment
The proposed FD step-size adjustment algorithm is embedded within the TR-optimization routine of Section II. Given the availability of the starting point y (0) and the vector of (userdefined) initial perturbations h (0) , the entire procedure can be summarized as follows (see Fig. 7 for a block diagram).
2) Construct G (i) model based on EM responses obtained at the H (i) test designs and set y (i+1) = y (i) , h (i+1) = h (i) , and i = i + 1; set δ (i) = 1. 3) Generate G (i) model using the method of Section III-D and its nested subroutines of Section III-C. 4) Find y (i+1) by solving (2) with U = U obj , and obtain ρ (i+1) . 5) If ||y (i+1) -y (i) || < ε, set y * = y (i) and END; otherwise, go to step 6. 6) Adjust δ (i+1) according to ρ (i+1) as in Section II-B. 7) If ρ (i+1) < 0, set δ (i) = δ (i+1) and go to step 4; otherwise, go to step 8. 8) If termination condition is met, set y * = y (i+1) and END; otherwise set i = i + 1 and go to step 3. Upon optimization, the design y * is rescaled back to the search space ranges determined by l b and u b vectors (cf. Section II) to obtain the final solution x * to the problem. The computational cost of the procedure corresponds to about D(χ + 1) + 1 EM simulations for the first TR algorithm iteration and χ D + 1 EM evaluations for the remaining (successful) steps. The iterations that do not improve the objective function require additional EM simulations. The parameter χ = 1.2 represents the average computational overhead necessary to ensure that the signs of the predicted and obtained objective function changes-determined through the adaptive perturbations adjustment procedure-are the same.

IV. NUMERICAL RESULTS
In this section, the performance of the proposed TR-based optimization algorithm with adaptive step-size adjustment is demonstrated based on three case studies that include the design of: 1) an ISM antenna; 2) a UWB radiator; and 3) a broadband rectifier for RF energy harvesting. The structures have been used to generate a total of 12 designs (including a set of ISM antennas obtained through optimization from various starting points). The proposed design procedure has been thoroughly benchmarked in terms of the cost and performance against the start-of-the-art TR algorithms that do not involve in-between-iterations tuning of FD steps. The initial number of significant digits and perturbations vector for optimization of all considered designs are set to k init = 1 and h (0) = 0.005·1 T , respectively. The rules concerning the update of the radius and TR algorithm termination are specified in Section II-B. All the computations have been performed on a 16-core Intel Xeon machine with 32-GB RAM. A case study related to the application of the UWB antenna in a hardware layer of a real-world IoT system is presented in Section V.

A. Planar Antenna for ISM Applications
The first considered structure is a planar quasi-patch antenna for ISM-band applications shown in Fig. 8 [39]. The structure is implemented on a Taconic RF-35 dielectric substrate (ε r = 3.5, h = 0.762 mm, and tanδ = 0.0018). It consists of a driven element in the form of a quasi-patch radiator with an inset feed, loaded with the monopole component. The ground plane length below the patch is adjustable. The antenna The objective function for the structure design is U obj (x) = max{R(x)} fL≤f ≤fH , where R(x) = P(y) (see Section II-B) represents the antenna reflection (in dB) within the f L = 5 GHz to f H = 6 GHz frequency range of interest, and f is the frequency sweep. The structure has been optimized using the proposed TR algorithm with adaptively adjusted FD steps. A total of ten test cases comprising the designs randomly generated within the box defined by the l b and u b bounds have been considered. The initial and optimized parameter vectors are gathered in Tables II and III, respectively, whereas responses of the selected final designs are shown in Fig. 9. It should be emphasized that the optimized characteristics are virtually the same. The average cost of the antenna optimization-for the considered test cases-corresponds to 86.2 EM simulations (2.8 h of CPU-time) per design.  The routine of Section III-E has been benchmarked in terms of performance and cost against the TR-based optimization approaches with a priori determined FD steps. The perturbations have been selected as (i) a fraction of the initial design (normalized to the unit box), and (ii) a square root of the machine precision as in (9), and (iii) a result of the manual FD adjustments followed by visual inspection of the structure responses at the initial design. The numerical tests have been performed using the set of starting points in Table II. The optimization results, gathered in Table IV, indicate that the proposed method generates the designs that exhibit (on average) from 2.8 to 3.36 dB (median varies from 2.8 to 4.69 dB) lower reflection levels compared to the benchmark algorithms with a priori FD steps. At the same time, the average CPUtime resulting from the application of the methodology is from 11.5% to 24.5% higher (median: from 10.9% to 15.1%). The increased cost stems from the computational overhead related to the optimization of perturbations. It should be noted that for all but two test cases the proposed method obtained either the best solution or its cost was lower (from 9.6% to 41.9%) compared to the benchmark techniques (i.e., for x 4

and x 7 ).
The cost corresponding to optimization of x 1 and x 10 using methods (ii) and (iii) was lower compared to the proposed approach. Furthermore, the resulting designs are characterized by 0.55 and 1.7 dB lower in-band reflection levels. Slightly deteriorated numerical performance of the proposed algorithm for the x 1 and x 10 designs results from locally reduced sensitivity of the objective function to step-size adjustment which results in reaching the bounds for some perturbations. Low correlation between the variable and function changes induces adjustment of k for the affected dimensions which increases the design cost (cf. Section III-C).
The efficiency of the presented framework has also been compared-for the x 1 , x 2 , x 6 , and x 9 designs-against the TR method with central FD where perturbations are defined a priori as the fraction of the starting point (iv). The results gathered in Table IV indicate that the solutions based on central FD are characterized by improved performance compared to the ones obtained using the forward method, yet at the expense over twofold higher cost (on average). The latter stems from doubling the number of EM simulations required to estimate the derivative data, but also the improved quality of the linear models that may increase the number of successful (hence more expensive) iterations in the course of the optimization process. Regardless of the performance gains, the results obtained using (iv) are still inferior compared to the ones generated by the proposed algorithm. Furthermore, the average cost of the presented method is 28.9% lower than for the central FD-based optimization.
For the sake of comparison, the adaptive adjustment of FD perturbations has been performed-for the design x (0) 1using a gradient tuning algorithm of [38]. The computational cost of estimating the derivative w.r.t. the first and second antenna parameters correspond to 62 (2.11 h of CPU-time) and 53 (1.81 h) EM simulations, respectively. At the same time, the average cost of determining the acceptable perturbation size using the algorithm of Section III-D is around 2.4 (295 s of CPU-time) and 1.4 (172 s) EM evaluations per parameter for the first and the remaining (successful) iterations of the TR algorithm. To put that into perspective, the cost of the presented method is two orders of magnitude lower (22-to 44-fold) compared to the benchmark algorithm. Prohibitively high numerical overhead renders the gradient tuning method useless for the class of design problems considered in this work.
It should be emphasized that the presented algorithm does not involve any manual tuning or decision making in terms of FD steps other than the determination of h (0) and k init parameters. The adjustment of perturbations-in the course of the optimization process-is performed without any insight of the designer. Consequently, the setup of the method is substantially simpler compared to the manual tuning of FD. At the same time, it offers notable improvement of the optimized designs performance (on average) compared to the benchmark techniques.

B. UWB Monopole for IoT-Based Precise Localization
The second structure is a UWB monopole shown in Fig. 10 [41]. The antenna is designed on a Rogers RO4003C substrate (ε r = 3.55, h = 0.813 mm, and tanδ = 0.0027). It comprises a radiator represented using a spline curve which is excited through a 50-microstrip line with the edge-mounted The design process involves optimization of the structure for minimization of: 1) reflection, i.e., F 1 (x) = max{R(x)} fL≤f ≤fH (in dB) within the f L = 3.1 GHz to f H = 10.6-GHz band and 2) footprint, i.e., F 2 (x) = A·B. The objective function is given as where The function (14) enables the minimization of the structure footprint along with the in-band reflection while ensuring that the latter is maintained around the target level F t.1 = -10 dB. The coefficients β = 1000 and F t.2 = 180 control the shape of the objective function landscape so as to ensure its smoothness for close-to-optimal designs. The structure has been optimized using the algorithm summarized in Section III-E. The initial design x (0) = [10 6 16 0.8 The antenna size and in-band reflection are 169.9 mm 2 and −9.44 dB, respectively, which represent almost 6-dB improvement compared to the initial design with no increase of the occupied area. Antenna reflection characteristics at the initial and final designs are shown in Fig. 11.
The presented algorithm has been compared in terms of performance and cost against the TR routines with a priori selected perturbations (cf. Section IV-A). The results gathered in Table V indicate that the proposed method generates the smallest antenna with the best in-band performance (by up to 1.5 dB). It is worth emphasizing that the design obtained using the presented approach offers a slightly reduced size of the structure compared to the starting point (170 mm 2 ), while the solutions generated using the benchmark methods are characterized by 5%-8% larger dimensions. Note that size reduction stays in conflict with maintaining high antenna performance [19], [39]. In the course of the optimization, the minimization is performed close to the edge of the feasible search space region which challenges the performance of local-search routines [11], [19]. Hence, the improvement of the antenna performance by up to 1.5 dB, as compared to the designs generated using methods with a priori selected FD steps is not negligible. However, the improved performance has been achieved at the expense of increased computational cost.

C. Rectifier for Wireless Energy Harvesting
The last considered structure is a broadband rectifier with the dual converter dedicated to harvest the RF energy shown in Fig. 12 [42]. The circuit is implemented on the Rogers RO4003C substrate (ε r = 3.55, h = 0.813 mm, and tanδ = 0.0027) and consists of a diplexer with two Dicksonpump converters with a common output. The RF-to-DC part of the rectifier comprises SMS7630 Schottky diodes, as well as C 1 = 100 pF and C 2 = 10 nF capacitors. It is connected to the high-and low-frequency transformers (arranged in a diplexer configuration) through the L 1 = 15 nH and L 2 = 3. The design goals are: 1) minimization of the rectifier reflection R 1 (x) expressed in dB and 2) maximization of its RF-to-DC conversion efficiency R 2 (x), both within the f L = 0.5 GHz to f H = 3.5 GHz band of interest. The objective function is given as where F 1 (x) = max{R 1 (x)} fL≤f ≤fH and the objective function coefficients are defined as The efficiency-related performance figure is expressed as F 2 (x) = -min{|R 2 (x)|} fL≤f ≤fH , whereas the threshold values  Fig. 13. The in-band reflection and efficiency of the optimized design are -9.9 dB and 30.1% (41.4% on average), respectively, which represents 6.5 dB and 13% improvement, compared to the initial design.
Numerical comparisons against the benchmark methods gathered in Table VI indicate that the proposed algorithm generated the best design (in terms of the objective function values). As it comes to the considered performance figures, the solution obtained using the presented method offers the lowest in-band reflection, while its worst case in-band efficiency is slightly (up to ∼0.15%) than lower compared to the designs obtained using algorithms (i) and (iii). It should be noted, however, that the computational cost of the algorithm is competitive compared to the benchmark methods-with the exception of (ii) which produced inferior design. Also, as for the example of Section IV-B, the approach that involves manual adjustment of the performance offers second to best objective function response for the optimized design.  numerical cost. The results indicate that for all of the considered test problems the proposed method generates the best solutions (the ISM antenna data points are based on the average of all considered test cases; cf. Table IV). The methods with a priori selected PE perturbations offer solutions that are characterized by 6.2%-48.7% worse performance. At the same time, the computational cost of the presented algorithm is from 6% upto 75.4% higher compared to the cheapest benchmark techniques. It is worth noting that, for the rectifier circuit, the cost of the presented approach is around 28% lower compared to the method (iii). Furthermore, the results show that the presented framework exhibits especially high cost compared to the solutions obtained using algorithm (ii). On the other hand, the designs generated by the latter are inferior by a large margin in terms of performance (from 28.1% to 48.7%, respectively). The design cost using algorithms (i) and (iii) is up to 47.2% lower compared to the presented method. Consequently, the numerical experiments indicate that the proposed framework is capable of obtaining competitive gains in terms of the optimized designs performance compared to the benchmark TR algorithm with a priori FD, yet at the expense of increased cost and complexity (implementation-wise) of the method.

D. Discussion of the Results
It is worth noting that the presented framework provides the means for dynamic estimation of the FD step sizes at a substantially lower cost compared to conventional approaches that exploit multiple simulations of the structure at hand per parameter to estimate acceptable perturbations [30], [37], [38]. To put that into perspective, for the Antenna of Section IV-A, the combined cost of adjusting FD steps (for just two out of six design parameters) using the method [38] amounts to 116 EM simulations (3.96 h of CPU-time). Consequently, it is from 4% to 61% higher than the optimization cost of all but two (x 1 and x 3 ) designs using the presented TR routine (cf . Table IV).
Finally, due to the local nature, the performance of the proposed TR algorithm (in terms of the optimized designs quality) depends on the selected starting points [11], [16], [17], [21]. On the other hand, the method proved to be useful for handling the problems with up to 38 dimensions. Such design tasks are considered prohibitively expensive for handling by conventional global optimization algorithms such as population-based metaheuristics [19], [44].

V. APPLICATION CASE STUDY
In this section, the applicability of the UWB monopole optimized using the presented TR algorithm as a part of the IoT, RTLS is demonstrated. The experimental validation of the structure far-field and electrical properties is followed by the concise discussion of the IoT-based indoor localization, and comparative studies concerning the positioning performance of the considered UWB-RTLS when equipped with the prototype and stock (i.e., delivered with the system) radiators.

A. Experimental Validation of the Spline UWB Antenna
The optimized radiator of Section IV-B has been fabricated and measured. Photographs of the manufactured prototype and comparison of the reflection characteristics obtained from the EM simulations and measurements are shown in Fig. 15. The in-band performance of both responses is similar. The discrepancy between the in-band maxima is only 0.65 dB. It is worth noting that slightly different shapes of the characteristics result from the current flow on the outer shell of the measurement cables, which alter the electrical size of the structure (as seen from the ports of the vector network analyzer), [45], [46]. The effect is calibration-independent.
A comparison of the radiation patterns obtained for the antenna in the x-y plane is shown in Fig. 16. The agreement between the responses is acceptable. Note that, although the radiation properties have not been considered as the design objective, small dimensions of the radiator promote its omnidirectional behavior [19]. Slight deterioration of the radiation performance around -90 • (more pronounced for the 3.5-GHz frequency) is due to the L-shaped stub which acts as a reflector for the far-field radiation.
It should be noted that, within the frequency range from 5 to 6.5 GHz, the antenna features close-to-omnidirectional radiation characteristics, low reflection, and high resemblance between the simulated and experimentally obtained responses. Mentioned performance figures make it useful for application in the hardware layer of the UWB localization system.

B. UWB-RTLS
Consider the UWB-RTLS dedicated to operate in a monitorbased localization (MBL) architecture [47], [48]. The system consists of-precisely synchronized in time [49]-reference nodes (so-called anchors) that register the arrival time of the pulses received from the transmitters (so-called tags) [47]. The latter ones are to be localized in the area covered by the system. The positions of tags are computed by the localization server using the time difference of arrival method [47]. The technique exploits the time delay between the reception of the signal by individual anchors (the location of which is known) to determine the position of the transmitter through solving the system of hyperbola equations. The discussed MBL architecture can be briefly summarized as follows: the tag sends a pulse signal, which is received by the anchors at the time instances that depend on the transmitter-receivers distances. The registered time signatures are then transmitted to the localization engine, which calculates the position of the tag [50]. A conceptual illustration of the system components is shown in Fig. 17.
It is worth noting that, due to challenging propagation conditions, UWB radio technology is particularly attractive for indoor positioning [47], [50]. In contrast to narrowband systems [47], [51], the transmission of short pulses in UWB increases its immunity to interferences and aids in distinguishing the line-of-sight signals from the reflected ones. Furthermore, UWB localization is supported by standardization, availability of the specialized integrated circuits, as well as dedicated system architectures [47], [48], [49], [50], [51]. The mentioned aspects are important for the development of reliable in-door positioning systems capable of obtaining both high precision and accuracy [47], [48], [49], [50], [51].

C. UWB-RTLS With Spline-Based Monopole Radiators
The UWB-RTLS system shown in Fig. 17 consists of five anchors installed in a rectangular test site (840 cm × 450 cm × 305 cm) at the height of 289 cm and two tags. The nodes are connected to the localization engine using an Ethernet-based local area network [48]. The experiment considered here involves the analysis of the positioning accuracy η and the (worst-case) precision σ within the monitored area. Note that the localization of the tag is represented in two dimensions [47], [48]. The performance metrics are expressed as For the considered experiments, the anchors are equipped with: (i) commercially available antennas and (ii) the manufactured prototypes of the UWB structure. The tests are performed at ten reference points uniformly distributed at the edges of the rectangle with known dimensions (cf. Fig. 17). For the anchors in configuration (i), the measurement campaign involved sequential determination of the tag location at the reference points. The position data have been obtained as a result of two sets of measurements for each tag which amount to a total of 40 experiments (N = 10 and M = 4). The reason for using two transmitters is to mitigate the potential bias resulting from slight physical differences between specific realizations of their transceiver sections. The same sequence of 40 measurements has also been performed for the anchors in (ii) configuration.
The accuracy and precision data obtained in setups (i) and (ii) are shown in  With the averaged discrepancies in the range of only 1 cm, the localization performance of the system (cf. Fig. 17) equipped with stock and prototype antennas is virtually the same. Although certain differences between precision and accuracy-attributed to distinct farfield performance of the considered radiators-can be noticed at specific points (cf. Table VII), they are within the acceptable ranges (especially, given the capabilities of UWB-RTLS) [48], [50]. Furthermore, the optimized prototype antenna is smaller compared to the stock radiator, which could be leveraged in the future for miniaturization of the anchors. The obtained results implicitly demonstrate the usefulness of the proposed TR-based optimization algorithm for real-world IoT applications.

VI. CONCLUSION
In this work, a TR framework that enables adaptive adjustment of FD perturbations has been proposed. The method involves (in-between TR iterations) optimization of the FD steps in order to elevate their effects on changes of the objective function values w.r.t. the center design. The refined perturbations are then used for construction of the first-order Taylor-expansion surrogates exploited in the course of the TR-based circuits optimization. Acceptable computational cost of the approach is maintained by reuse of the existing linear models for adjustment of the perturbations. The presented approach has been thoroughly validated based on three different structures for IoT applications, i.e., quasi-patch antenna, UWB monopole, and rectifier for RF energy harvesting. A total of 12 test cases have been considered. The method has been benchmarked against the state-of-the-art TR algorithms with the static, a priori selected perturbations. The numerical results indicate that, for the considered test cases, the presented algorithm offers up to almost 50% improvement of the optimized designs performance, yet at a higher cost compared to the reference methods. Notwithstanding, the CPU overhead of the approach is low compared to conventional algorithms that require dozens of EM simulations per parameter to adaptively determine the optimal FD step.
Usefulness of the proposed algorithm for real-world IoT systems has been implicitly demonstrated through the utilization of the optimized spline-based UWB monopole antenna prototypes in the hardware layer of the real-time system for indoor positioning. The results of comparative studies proved that, for the considered test setup, the substitution of commercially available antennas with prototype structures has negligible effects on the localization performance. At the same time, small dimensions of the spline-based radiators could be leveraged in order to develop miniaturized nodes of the RTLS.
Future work will focus on the reduction of the method numerical cost through automatic preselection of the variables characterized by the largest sensitivity of the objective function to the FD step changes. Application of the approach to other research areas that heavily rely on the use of computationally expensive simulation models, as well as embedding the algorithm into a surrogate-assisted design framework, with FD steps estimated based on the low-fidelity EM simulations, will also be considered.