Adaptive Synthesis Using Hybrid Genetic Algorithm and Particle Swarm Optimization for Reflectionless Filter With Lumped Elements

In this article, an adaptive synthesis based on the hybrid genetic algorithm and particle swarm optimization (HGAPSO) is proposed for reflectionless filter design with lumped capacitors, resistors, and inductors. The synthesis starts with a preset topology, where each branch of the topology represents a small passive network of lumped elements. The proposed HGAPSO is used to trim branches and obtain proper values of elements for a required filtering response. Focus on this synthesis model, the HGAPSO is embedded with local searching policies based on random coordinate and neighborhood search to improve its searching ability. Besides, a classifier-based strategy and a probabilistic method are introduced to accelerate convergence and boost iteration. Suitable topologies and component values are determined automatically by the HGAPSO to meet the specific filtering response. To predict the response accurately, the EM-simulated result of the corresponding layout and parasitic parameter models of lumped elements are considered during the fine-tuning. Based on the mechanisms mentioned above, four reflectionless bandpass filters (BPFs) are synthesized to validate the effectiveness of the proposed synthesis procedure. The fabricated filters exhibit good selectivity and low reflection coefficient in the measurement.


Adaptive Synthesis Using Hybrid Genetic Algorithm and Particle Swarm Optimization for Reflectionless
Filter With Lumped Elements Yifan Li , Graduate Student Member, IEEE, and Xun Luo , Senior Member, IEEE Abstract-In this article, an adaptive synthesis based on the hybrid genetic algorithm and particle swarm optimization (HGAPSO) is proposed for reflectionless filter design with lumped capacitors, resistors, and inductors.The synthesis starts with a preset topology, where each branch of the topology represents a small passive network of lumped elements.The proposed HGAPSO is used to trim branches and obtain proper values of elements for a required filtering response.Focus on this synthesis model, the HGAPSO is embedded with local searching policies based on random coordinate and neighborhood search to improve its searching ability.Besides, a classifier-based strategy and a probabilistic method are introduced to accelerate convergence and boost iteration.Suitable topologies and component values are determined automatically by the HGAPSO to meet the specific filtering response.To predict the response accurately, the EM-simulated result of the corresponding layout and parasitic parameter models of lumped elements are considered during the fine-tuning.Based on the mechanisms mentioned above, four reflectionless bandpass filters (BPFs) are synthesized to validate the effectiveness of the proposed synthesis procedure.The fabricated filters exhibit good selectivity and low reflection coefficient in the measurement.

I. INTRODUCTION
F ILTER plays a critical role in microwave systems.It is widely used to select the desired signals and suppress the interferences at unwanted frequencies.Conventional filters have relatively consummate design procedures based on lumped elements, transmission lines, or coupled resonators [1], [2], [3].However, such filters are usually reflective, where the unwanted out-of-band interferences are reflected back to the input terminals.Those reflected signals generate additional intermodulation products at the nonlinear front-end circuit (e.g., mixer).Meanwhile, the reactive input impedance of a reflective filter leads to system instability once the filter is connected with an amplifier [4], [5], [6], [7].Reflectionless filters using lumped components [4], [5] are developed to absorb the interferences at unwanted frequencies.Then, the performance of microwave systems can be improved when active modules are connected with reflectionless filters.Besides, the reflectionless filters replace the conventional ones with isolators or circulators in integrated systems, which improve system integration and reduce design cost [7], [8], [9].To improve the performance of reflectionless filters, two main methods are proposed.By using two networks with complementary frequency responses, the reflectionless operation can be obtained at one port [10], [11], [12], [13], [14].Such filters usually exhibit the merits of low reflection, high selectivity, and wide stopband.The other method is using symmetrical topologies to implement dual-port reflectionless filters [14], [15], [16], [17], [18], [19], [20], [21], [22].Such filters are featuring two-port reflectionless operation in a wide frequency range.Based on lumped elements, filters with ideal reflectionless response have been developed with specific topologies [23], [24].However, most of them require a number of inductors or transformers.The non-negligible parasitic parameters have a significant influence on the matching and degrade the reflectionless performance of manufactured filters.Meanwhile, transformers have larger sizes compared with other lumped components, which limits the placement of components and circuit connections on the PCB.To overcome the design challenges caused by parasitics and reduce the time consumption of finding suitable circuit topologies, an adaptive synthesis procedure is necessary.
Recently, adaptive optimization methods [25], [26], [27], [28] are introduced for filter design.In optimizing a filter with several design parameters and constraints, there are various solutions that can approach the desired requirements.Therefore, the target function has multiple local optimums, which makes it nonconvex.Up to now, global nonconvex optimization is still a great challenge and no algorithm can guarantee to find the best solution without traversing all possible solutions.Even so, heuristic methods, e.g., particle swarm optimization (PSO)-based algorithms [29], show remarkable performance in circuit synthesis and optimization [25], [28].In [25], PSO has a significant role to optimize the coupling matrix for designing coupled-resonator filters.Nevertheless, in high-dimensional solution space, particles are always sparse, which makes particles miss the global optimum and be congregated to a local optimum.To overcome this drawback, PSO and genetic algorithm (GA) are combined to form the hybrid genetic algorithm and particle optimization (HGAPSO) [30], [31], [32].The HGAPSO features the advantages of searching heuristically and reduces the probability of being trapped into local optimums.Nonetheless, once a time-consuming target function is used, PSO-based algorithms are slow in iteration.Introducing gradient-based local search is an intuitive method to accelerate PSO [33].However, in high-dimensional optimization, calculating gradient vectors still significantly decreases the speed of iteration.To avoid calculating gradients, algorithms using classifiers have been reported, which show high efficiency in high-dimensional optimization [34], [35].
The purpose of this work is to achieve the adaptive synthesis of proper topologies for reflectionless filters.In the implemented procedure, the inputs are constraints on the filtering response and range of component values, while the outputs are the synthesized topology and proper values of all components.The whole synthesis procedure is shown in Fig. 1.This procedure starts with a preset graph G 0 , and circuits are represented by trimming branches of G 0 .Constraints on G 0 can be applied to prevent complicated circuit structures and avoid unnecessary parasitics from the layout.After introducing the general branch model, the trimming procedure is realized by changing the variables' values.Such operations convert the synthesis into a high-dimensional nonconvex optimization.Focusing on this optimization, an enhanced HGAPSO is proposed, which embeds local search strategies and improves the heuristics ability of PSO-based algorithms.Meanwhile, a sorting-based classifier is introduced to accelerate convergence, and probabilistic methods are used for reducing unnecessary calculations.Once the preset graph G 0 is able to cover feasible topologies, the proposed HGAPSO can obtain the desired circuit.If the optimization could not get a satisfying performance after a long iteration, the optimization should be paused.Then, the graph is extended to improve its capability of expressing circuits and the optimization procedure is restarted.When the optimization using the proposed HGAPSO is finished, a suitable topology and corresponding values of all components are generated.After the layout is designed, a gradient-based algorithm is utilized for fine-tuning components' values.The influence of the parasitic parameters from the layout and components is considered to predict the response of practical filters accurately, which ensures the fabricated filters with expected characteristics.
The rest of this article is organized as follows.In Section II, the proposed general branch model and related network theory are introduced.Then, the procedure for constructing the target function which maps circuits into vector space is described.The details of the proposed HGAPSO algorithm and gradient-based fine-tuning method are introduced and compared with traditional algorithms in Section III.In Section IV, the schematic, layout, and response of synthesized reflectionless filters are presented and discussed.Those filters are fabricated and measured for verification.In Section V, the conclusions of this article are provided.

II. PROPOSED MODEL AND RELATED THEORY A. Expressing Circuits in Graph With Proposed General Branch Model
As mentioned in [3] and [36], a passive lumped circuit can be expressed using an undirected graph described by an incidence matrix.To avoid a variable number of components changing the dimensionality of the solution space, a predefined initial graph G 0 is used to fix the total number of components.G 0 can be a complete graph, e.g., the graph K 6 as shown in Fig. 2(a), to guarantee the ability to synthesize any graph with an equal or fewer number of nodes.In this article, for a graph of N nodes, port 1 is defined as node 1, port 2 is defined as node 2, and ground is defined as node N (i.e., the ground-node).For each pair of node i and node j, the corresponding branch b ij is composed of three parallel components, i.e., a resistor R ij , an inductor L ij , and a capacitor C ij , as shown in Fig. 2(b).Thus, the admittance of branch b ij is expressed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The value of each component is limited to ensure the synthesized result is manufacturable.R min , L min , and C min are the lower bound of possible prototype value for resistors, inductors, and capacitors, while R max , L max , and C max are the upper bound.The components' values smaller than the lower bound are set to positive infinitesimal.Meanwhile, the components' values are set to infinite once the upper bound is exceeded.Therefore, by changing the value of components, each branch can approach all situations summarized as follows.
1) Trimmed Branch: R ij → ∞ and L ij → ∞ and The number of required components grows with the number of branches.For a complete graph K n , the number of branches is calculated using the following: Thus, a massive number of branches are introduced once a large n is used.Meanwhile, variables of high dimensionality lead to a large amount of computation, which slows down the synthesis.Therefore, except for the complete graph K n , other graphs are used as G 0 .Besides, if the synthesized topology is planar (i.e., all components can be placed on the same side of the PCB, and no crossover is needed in the signal path), parasitics introduced by the interconnections are reduced.G 0 is initially set as a planar graph.However, for the PCB with metal layers at both sides, the whole bottom layer is used as the ground.To take advantage of the bottom layer, a type of quasi-planar graph is introduced.In the quasiplanar graph, the subgraphs without the ground node are planar graphs, and each node has an individual branch connecting the ground node.Such quasi-planar graphs are used as G 0 since it is acceptable to connect the ground plane using vias.Thus, parasitics from vias and long wires across multiple metal layers in the signal paths are avoided.An example of such a quasi-planar graph with 13 nodes (including the ports and ground node) is presented in Fig. 3.In this graph, all nodes are placed and connected in a grid.Then, diagonal branches are added to extend the flexibility in the synthesis.Compared to a complete graph of the same node number, such graphs have fewer branches, and the number of total variables is reduced.
Graphs with more nodes are needed to synthesize circuits with higher order and larger scale.Meanwhile, if the target filter is a bandpass filter (BPF) or a bandstop filter (BSF), and the desired response is symmetric about the center frequency, the response can be converted to a low-pass type.Thus, a low-pass filter (LPF) is firstly synthesized.Then, the desired BPF/BSF is converted from the synthesized LPF prototype [1].During the conversion, each inductor in the LPF is replaced with a series-connected LC pair, while each capacitor is replaced with a parallel-connected LC pair.Therefore, synthesizing LPF requires G 0 with a smaller scale compared to BPF/BSF, which reduces the total computation and improves the synthesis efficiency.

B. Mapping Circuits Into a Vector Space: Element Denormalization, Circuit Solver, and Assessment
In this section, a vector ⃗ X is utilized to represent the components' values, where the dimensionality of the vector is D. And D is the same as the number of total components in G 0 .Note that the value of element x i in vector ⃗ X is limited from 0 to 1. Therefore, a circuit can be represented by G 0 and ⃗ X .Then, the target function includes the procedures of solving the circuit, assessing the response, and indicating the discrepancy between the response and constraints.For simplification, To assess the response of the circuits, the scattering matrix is obtained.Before calculating scattering parameters, ⃗ X is denormalized to get the actual values of components.In consideration of the finite practical component's value, a piecewise function is introduced.The actual value p(x i ) of component x i is expressed as where a and b represent the lower and upper bounds of x i , respectively.The exponential function is utilized due to the discrete and exponential distribution of available lumped components.Besides, to ensure p(x i ) is between the minimum and maximum component value, parameters in exponential function (i.e., k i and s i ) are determined as where V i,min is the available minimal component value and V i,max is the maximal one.When x i is less than a or more than b, the corresponding component is shorted or opened.Thus, m 1 should be small enough and m 2 should be large enough to ensure the calculated admittance is ideal.The reason for not setting m 1 and m 2 to zero and infinity is to avoid a singular matrix when calculating scattering parameters.
Meanwhile, when x i is in the range of [0, a) or (b, 1], the mapping function is still linear.Therefore, in the step of extending G 0 , if x i generated by the optimizer is close to 0 + or 1 − , then x i is used to determine whether the corresponding component should be removed.After the denormalization of ⃗ X , the scattering matrix of the corresponding circuit is solved using G 0 and denormalized component values.The first step of the solver is to calculate the primitive indefinite admittance matrix [Y p ] from the incidence matrix [A c ] of G 0 and the branch admittance matrix [Y b ].The matrix [A c ] is obtained from the topology of G 0 , which has the following form: where N and M are the number of nodes and the number of branches of G 0 , respectively.Each column of [A c ] represents a branch of G 0 .In the column representing branch b ij , the element of the ith row is 1 and the element of the jth row is −1 (e.g., the first column in ( 5) represents the branch b 13 with the direction from nodes 1 to 3).Note that in a passive reciprocal circuit, the predefined current direction in each column of The ith row and column element of [Y b ] is the admittance of the branch that the ith column of [A c ] indicates, and it is calculated using (1).Then, [Y p ] is calculated using following formula [3]: Since node N is the ground, the N th row and N th column of [Y p ] are directly removed to form the definite admittance matrix [Y d ].Then, the admittance matrix of this circuit is obtained by eliminating [Y d ] from node 3 to node (N − 1).This step is performed by applying the following [36]: Here, k denotes the node is eliminated.k begins from the last node (i.e., k = N − 1) to the third node (i.e., k = 3).In eliminating the kth node, ( 8) is applied to each element in [Y d ] other than the kth node by traversing i and j from 1 to (k − 1).Then, the kth row and column are removed, and k is decreased by 1.After k decreases to 3 (i.e., the third row and column are removed), the remaining [Y d ] with the first two rows and columns is the admittance matrix [Y ] of this circuit.The impedance matrix [Z ] is obtained by where [I ] is identity matrix, and Z 0 is characteristic impedance.
The scattering matrix [S] of each frequency ω is calculated from graph G 0 and vector ⃗ X , where the vector ⃗ X contains values of R, L, and C. The value of the target function f is calculated from [S] of each frequency.In constructing f , both the average and maximum distances of each frequency point are utilized.f is expressed as where f a,(i, j) denotes the average cost value calculated from scattering parameter S ij of all frequencies, and f m,(i, j) denotes the maximum cost value.f a,(i, j) is described as where r (x) is the ramp function defined as r (x) = max(0, x).Meanwhile, f m,(i, j) is described as In ( 11) and ( 12), the desired response is described by setting the upper bound S ij,UB (ω) and lower bound S ij,LB (ω) for each frequency point, as an example of constraints on reflectionless BPF shown in Fig. 4. When TZs are required, S 21,UB (ω) of the corresponding frequency points is set lower than nearby frequency points.The distance between each scattering parameter and the corresponding required value is passed into the ramp function.At a single frequency point ω 0 , the ramp function produces 0 only if the scattering parameter S ij (w 0 ) satisfies the requirements defined by S ij,LB (ω 0 ) and S ij,UB (ω 0 ).In (11), all ramp function outputs are summated, and the weighted averages are calculated.Quadratic operation is applied on the ramp function to make f tend to be convex, which can afford convergence in the optimization.Parameters of w a1 and w a2 should be configured to ensure f a,(i, j) contributes a larger proportion of f at the beginning of optimization.Thus, f a,(i, j) allows the target response to take shape, as illustrated in Fig. 5(a), where the optimizer needs to suppress all violations to decrease the value of f .However, (11) only calculates average distances.A small number of large violations would be ignored which allows unwanted ripples and peaks of the response to occur, as illustrated in Fig. 5(b).Equation ( 12) is introduced to suppress those unwanted ripples and peaks in the response.In ( 12), the highest ramp function output is selected which indicates the maximal violation.The weights of w m need to be configured to ensure that f m,(i, j) introduces less effect on f a,(i, j) when most of the requirements are unsatisfied.Meanwhile, f m,(i, j) can dominate f when only a few violations occur.Those parameters in (10) and (11) are set according to the specific filter performance.
Finally, circuits are represented by vectors defined in the solution space.The dimensionality of solution space is the same as the component number, and the definition domain is [0, 1] in each dimension.The target function f ( ⃗ X ) maps the vector ⃗ X to a single value.The goal is to search the solution space and find a suitable ⃗ X that makes f ( ⃗ X ) small enough or be zero.

A. Details of the Proposed HGAPSO
To obtain a suitable ⃗ X for the target function described in Section II-B, a specifically enhanced HGAPSO is proposed.This algorithm is based on the PSO algorithm.The strategy of PSO is that each particle moves toward its own best location and the global best location.This strategy gives the heuristic ability to PSO.To overcome the drawbacks of PSO in high-dimensional solution space, the mutation and selection parts of GA are introduced to construct the proposed HGAPSO.The mutation part is realized through random perturbation and classifier-assisted particle reinitialization.The selection part reduces unnecessary calculations by removing worse particles among congregated particles.The selection mechanism of the proposed HGAPSO is achieved through the lifespan and fitness of particles.Each particle has a parameter of lifespan with an initial value.In the iteration, the fitness of each particle is determined from the corresponding target function value.Then, lifespan is decreased according to fitness.Once the lifespan of a particle becomes negative, the particle is removed.Most of the removed particles are reinitialized in the optimal zone.The optimal zone is determined by the best particles using a classifier.Therefore, reinitialized particles can be regarded as variants of best particles, which form the mutation part and boosts the convergence of the proposed HGAPSO.Any type of classifier can be employed while a simple hyperrectangle classifier is used in this implementation.To avoid the optimal zone being very narrow and reduce the probability of converging to a local optimum, two policies are introduced.The first policy expands the upper and lower bounds of the hyperrectangle.The optimal zone is extended to prevent only the local optimum is covered.Another policy is allowing a portion of particles to be reinitialized outside the optimal zone.Thus, those particles explore solutions far from the particle swarm.In the experiment, such a strategy shows good performance in optimization while not taking up too much calculation.
The flowchart of the implemented procedure is shown in Fig. 6.Note that all variables and solution space are normalized and limited to [0, 1], while the target function is f ( ⃗ X ).The initialization part sets each element in vectors ⃗ C max to maximum and each element in ⃗ C min to minimum.Those two vectors describe the optimal zone of the classifier.Thus, the classifier initially covers the whole solution space of D dimensions.Here, D is the same as the component number.After that, the location of all particles, i.e., ⃗ X i , are initialized randomly in the solution space.Here, a function denoted as rand( ⃗ V ) is employed to generate a vector that has the same dimensionality as its input vector ⃗ V .The generated vector is located randomly inside the hyperrectangle defined by ⃗ V and coordinate origin.Then, the local best location ⃗ X LB,i of each particle is set to its initial location ⃗ X i .The lifespan L i of each particle is set to initial lifespan L init .The particle initialization process is executed in a loop which is carried for P times.Here, P is the number of particles.When all particles are initialized, the target function value of each particle is calculated.The location of the particle with minimum function value is selected as the global best point ⃗ X GB .Finally, the iteration counter is reset.This counter is employed to stop iteration once the number of total iterations reaches its limit of T .
After initialization, the iteration starts.In the normal iteration section, the subprocess shown in Fig. 7 controls the movement of each particle ⃗ X i .During the subprocess, the results using four movement strategies are calculated based on each particle's current location ⃗ X i .The first two movement policies, denoted as ⃗ X (1)   i and ⃗ X (2)  i , are based on the original PSO.Here, the particle moves toward the global best location ⃗ X GB with the weight of w g to location ⃗ X (1)  i .Then, it further moves toward local best location ⃗ X LB,i with the weight of w l from ⃗ X (1)   i to location ⃗ X (2)  i .The third strategy ⃗ X (3)   i introduces perturbation on current position ⃗ X i .In each dimension, the coordinate value of the particle is added with a random value.The scale of this random value is limited by the weight w r Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.multiplied by the length of the classifier's coverage of the corresponding dimension.Note that the gap in each dimension between ⃗ C max and ⃗ C min shrinks along with iteration, which makes the perturbation of ⃗ X (3)   i gradually become smaller.Thus, the resolution of random neighborhood search actualized by ⃗ X i increases with iteration.To calculate ⃗ X (4)  i , a function S( ⃗ X , D) is defined as follows: S( ⃗ X , D) randomly selects a dimension from all D dimensions.Then, S( ⃗ X , D) takes values in this dimension to generate a set of vectors, where values from other dimensions of ⃗ X are held.This step is noted as random coordinate search (RCS), which is a simplification of random direction search (RDS) [37], [38].The principle of RDS is that it randomly selects a direction vector, then take samples along the line defined by this vector.To balance the speed and search capability, in this embedded RCS, the direction vector is always parallel to a coordinate axis.In a single iteration, the best sample point from S( ⃗ X , D) is set as ⃗ X (4)  i .Note that the selected coordinate axis should be different for the next iteration.As shown in Fig. 8, it can help particles cross long-distance gaps between local optimums.After new positions of ⃗ X (1)  i -⃗ X (4)   i are calculated, the subprocess chooses the best one to replace the previous position of the particle ⃗ X i .If this new position is better than ⃗ X LB,i , ⃗ X LB,i is updated using ⃗ X i .When all particles are updated, the best ⃗ X LB,i is selected instead of the global best position ⃗ X GB .The subprocess is called independently for each particle which makes it suitable for parallel executing.Besides, to balance the speed of iteration and heuristic ability, a method inspired by simulated annealing (SA) is introduced in the subprocess of calculating ⃗ X (3)  i and ⃗ X (4)  i to accelerate iteration.In SA, the probability is utilized to determine whether to accept a new solution that is worse than the current solution.This probability decreases as the algorithm converges.In the proposed strategy, the probability depends on f ( ⃗ X GB ) and the maximum value of the target cost function f max .For the embedded RCS of calculating ⃗ X (4)  i , whether a single component is in the short zone or the open zone can be confirmed easily.Thus, it takes effect at the beginning when the topology is undetermined.In late iteration, better positions are hard to be found by RCS since the topology is stable.Meanwhile, since variables change significantly in the early iteration due to ⃗ X (4)  i , the random disturbance of ⃗ X i could not lead to contributions.If particles are trapped in a local optimum, random disturbance could afford particles to jump out to better locations nearby.Besides, because of the shrunk optimal zone, ⃗ X i searches the solution space finely when f ( ⃗ X GB ) is small.Therefore, the probability distribution for calculating ⃗ X (4)  i has a higher value when f ( ⃗ X GB )/ f max is near 1 other than 0. The probability for calculating ⃗ X (3)   i has the opposite distribution to ⃗ X (4)  i .To satisfy characteristics of such requirements, probability distribution functions using tanh are constructed as where parameters k 1 and k 2 are used to adjust the slope in the middle.The distributions of those functions are shown in Fig. 9.By introducing this probability-based strategy for deciding whether to calculate ⃗ X (3)   i and ⃗ X (4)  i , the amount of calculation for each iteration is nearly reduced by half.Meanwhile, the converging performance is almost unchanged.
After the normal iteration, the classifier needs to be updated, which covers optimal particles by defining the optimal zone.The optimal zone is used in subprocess 1 and the particlereinitialization, as an example shown in Fig. 10.Here, a simple classifier update strategy is introduced, which can be based on any sorting algorithm.The sorting function sorts target function values corresponding to each particle in order from lowest to highest.After all particles are sorted, the sorting function returns the function value list SV and index list SI of particles.Though sort is used in the classifier updating part, it can be carried out before updating ⃗ X GB in the normal iteration part.Then, ⃗ X GB is directly set as ⃗ X SI 1 if SV 1 is smaller than f ( ⃗ X GB ).After the sorting process, the best K particles are selected.From such K particles, the function max D returns the maximum value of each dimension to form the upper bound ⃗ C max , while the function min D returns the minimum value of each dimension to define the lower bound ⃗ C min .To avoid particles being aggregated to a local minimum, a vector ⃗ G is employed to extend the gap between ⃗ C max and ⃗ C min in each dimension.Here, the vector ⃗ G is generated by setting each element to (G/2), where the parameter G denotes the extending length in each dimension.After ⃗ C max and ⃗ C min are calculated, ⃗ C max is increased by ⃗ G, and ⃗ C min is decreased by ⃗ G.Note that each element of ⃗ C max and ⃗ C min is defined between 0 and 1.Thus, the value out of range is fixed to the maximum of 1 or the minimum of 0.
When ⃗ C max and ⃗ C min are updated, all particles are already sorted according to their value.The index list SI generated from the sorting function is utilized as fitness to calculate the lifespan of each particle.Here, a simple strategy is imported that each particle's lifespan L i is decreased by their sort index SI i (i.e., the lower the target function value of the particle, the less the lifespan is reduced.).Once L i becomes negative, the ith particle is reset.Such a procedure consists of two steps.The first step is resetting parameters, while the second step is executing the preoptimization.In the first step, the negative L i is reset to the initial value L init .Then, a random value between 0 and 1 is used to determine whether to reset ⃗ X i in the optimal zone.If the random value is larger than the probability p, and ⃗ X i is reset in the whole solution space.Otherwise, it is reset in the optimal zone defined by ⃗ C max and ⃗ C min .To help the reinitialized particles catch up with existing particles in iteration, preoptimization is applied to those particles in the second step.The preoptimization uses the same procedure by calling subprocess 1 with repeated W times, where W increases with I C .Here, I C is the abbreviation of IterationCount in Fig. 6.In this implementation, W is calculated as Since W is constructed using a logarithmic function, the preoptimization procedure does not slow down the iteration  ) is smaller than f ( ⃗ X GB ), ⃗ X GB should be updated using ⃗ X i .Finally, when the procedure described above is finished, an iteration is completed, and I C is increased by 1. Unless I C reaches its limitation T or f ( ⃗ X GB ) satisfies the requirement of f target , the iteration loop is continued.Once the iteration loop is stopped, the global best location ⃗ X GB is the final result.

B. Comparison of Heuristic Algorithms
To demonstrate the advantages of the proposed HGAPSO, the function of (10) with the same parameters for synthesizing reflectionless filters is applied to different algorithms for comparison.The results of convergence are shown in Fig. 11.Each algorithm is executed 50 times while the convergence curve of the best result is selected.The proposed HGAPSO is denoted as "HGAPSOCRS" in the legend of Fig. 11.The "C" suffix stands for "classification."The "R" suffix stands for "random neighborhood search" which is the step of calculating ⃗ X (3) .The "S" suffix stands for RCS which is the step of calculating ⃗ X (4) .In these tests, the same hyperparameters are used for all algorithms, except for the gradient descent (GD) and RDS, which have no concept of particles.The hyperparameters used are listed in Table I.The details of selecting hyper-parameters are presented in the Appendix.Besides, since GD tends to fall into local optimums, a preliminary procedure of randomly generating initial positions is carried out.Then, the best position of 1000 initial positions is selected to start GD and RDS.As illustrated in Fig. 11, particles of PSO can be trapped into local optimums in high-dimensional optimization.Compared to PSO, HGAPSO has the ability to prevent particle aggregation, but particles still lack searching ability in such a high-dimensional solution space.When classification is introduced, convergence is boosted.However, classification cannot improve searching ability.Meanwhile, the embedded optimizations of random neighborhood search and RCS ameliorate the heuristic ability significantly.In Fig. 11(b), the target function of synthesizing wideband reflectionless BPF is used, which has more sample points and tighter passband constraints.In both comparisons, the proposed HGAPSO shows the ability of fast convergence and better final results.

C. Design Layout and Fine-Tune
The topology of the synthesized filter is determined by referring to ⃗ X GB generated from the proposed HGAPSO and its corresponding graph G 0 .Then, the layout is designed along with substrate characteristics, components' package, and circuit terminations.In this step, all normalized variables are denormalized to actual component values.When the circuit is synthesized directly (i.e., the cost function describes a bandpass response instead of a low-pass response), the LPF prototype converting formulas from [1] is utilized which are listed below where g Li , g Ci , and g Ri are prototype values of the inductor, capacitor, and resistor, respectively.Z 0 is the characteristic impedance of the transmission line connected to the terminals, while ω c is the desired center frequency.
If the target is a BPF but the response is described in a low-pass response, the LPF is synthesized.The low-passto-bandpass transformation is performed to convert the synthesized LPF into BPF.Here, each inductor is connected in series with a new capacitor, while each capacitor is paralleled with a new inductor.The prototype value of the added component in each pair is the reciprocal value of the original component value [1].After all components in the prototype are converted, ( 15) is applied to calculate the actual values.Besides, if the target is a BSF or high-pass filter, similar procedures described in [1] can be applied for converting.
During the circuit implementation, the influence of the layout is considered.The scattering matrix of the layout without lumped components is obtained through EM simulation.Then, the indefinite impedance matrix is calculated by using ( 9) in reverse as follows: The admittance matrix of the layout (i.e., [Y l ]) is the inverse matrix of the impedance matrix [Z ].Then, the admittance matrix of each component is added to the corresponding rows and columns of [Y l ] to form the indefinite admittance matrix [Y p ].After [Y p ] is formed, the method described in (8) and ( 9) is utilized to calculate scattering matrix [S] of the circuit.The assessment procedure in (10) is applied to finetune components' values considering the effect of the layout.Besides, the frequency in the target function is changed to the denormalized frequency.
The parasitic parameters of the packaged lumped components are considered once self-resonator frequencies are near the passband frequency.The models of parasitic parameters employed in this work are shown in Fig. 12.For inductors, C p is the parasitic capacitance of the package, which is fixed for a specific inductor series.R s is the metal resistor that is approximately linear to inductance L. It is expressed as a function R s (L).Thus, the admittance of an inductor with parasitic parameters is calculated as For capacitors, L x and R x are the parasitic inductance and resistance of the metal, respectively, which are approximately linear to 1/C.Therefore, functions L x (C) and R x (C) are introduced to model parasitics of The admittance of a capacitor with parasitic parameters is calculated as Parasitic parameters are obtained by referring to the self-resonator frequency and Q from datasheets.To assist the optimization, linear spline interpolation is used, which is expressed as follows: where x i and x i+1 are two adjacent component values, s i and s i+1 are the corresponding parasitic parameters, and Ŝi (x) is the predicted parasitic value, respectively.The domain of interpolated function Ŝi (x) is [x i , x i+1 ].During the finetuning, discrete parasitic parameters are modeled as piecewise continuous functions using (19).After considering the parasitic parameters of both layout and lumped components, reflectionless characteristics deteriorate.Since the components' values generated from the HGAPSO are near the optimal solution, it is effective to use a gradient-based method for further optimization.The flowchart is shown in Fig. 13.Parameter D 0 limits the minimum length of the gradient vector.Once the length is small enough, it means that the shape of the function f is flat enough near ⃗ X .Thus, ⃗ X is returned as the final value.Parameters step 0 and step min control the sweep step in the direction of the gradient.During the sweeping, if f ( ⃗ X * ) is smaller than f ( ⃗ X cur ), ⃗ X cur is updated using ⃗ X * , then ⃗ X will be replaced by ⃗ X cur after the sweeping loop.When step is too large, i.e., the first attempt strides across all regions where the value of f is lower than f ( ⃗ X cur ), the branch of t = t + step is not triggered.Therefore, t is equal to step which triggers the branch of reducing parameter step by half.Once step reaches the minimum limitation of step min (i.e., current location ⃗ X is close enough to the optimal solution), the iteration is stopped.Besides, the variable count that functions as an iteration counter are employed.If count reaches its limit count max , the iteration is terminated, and ⃗ X is the final result.

IV. DESIGN EXAMPLE
In this section, four synthesized reflectionless filters are presented to demonstrate the effectiveness of the proposed procedure.The grid-type quasi-planar graphs, which have an example of 13 nodes illustrated in Fig. 3, are used as the final G 0 for all results.The layouts are designed according to the topologies synthesized by the proposed HGAPSO.Then, the component values of each filter are denormalized and fine-tuned with the corresponding layout using gradientbased optimization.Note that parasitics of both layouts and components are considered in the fine-tuning.The closest practical value of each component is selected to fabricate the filter, where the inductors are from the muRata LQW15AN series and the capacitors are from the muRata GJM1555 series.

A. Reflectionless Filter With Deep Stopband Attenuation
The first synthesized result is a reflectionless BPF with deep stopband attenuation.This design example is also used to validate that topologies of fewer components can be found by the proposed procedure, while passband characteristics are maintained.By inputting a synthesized topology as G 0 to the synthesis procedure, a simpler topology is obtained if G 0 can be simplified.This procedure is repeated multiple times until G 0 can no longer be further simplified.To demonstrate this functionality, the three-stage reflectionless filter from [4] is used as the reference of response.The passband and TZs are constrained as the reference, while the stopband rejection is set to 55 dB, which is 10 dB deeper than the reference.Besides, the topology is also constrained symmetrically due to the reference also having a symmetrical topology.The reflection coefficient is limited to −30 dB.All hyper-parameters used in synthesizing this filter are listed in Table I.The schematic of the final synthesized filter is shown in Fig. 14.This filter only needs 15 inductors, 15 capacitors, and 4 resistors.Compared to the referenced filter, three inductors, three capacitors, and two resistors are saved.The calculated responses from schematics are presented in Fig. 15.The synthesized filter has a passband that is nearly identical to the reference.The TZs are added at the same frequency, while the insertion peak in the stopband is 58.3 dB.
In fabricating this filter, the center frequency is set to 433 MHz.The designed layout is shown in Fig. 16.The PCB is implemented using the dielectric substrate of RO4350B (i.e., typical ϵ r of 3.66 and a thickness of 0.762 mm).Values of practical components are presented in Fig. 14.Photographs of the fabricated filter are shown in Fig. 17.The measured, simulated, and calculated responses are shown in Fig. 18.In Fig. 18(a), the normalized frequencies are predicted using the following: where is the normalized frequency, FBW is the fractional bandwidth (FBW) and is set to 1, f is the original frequency, and f 0 is the center frequency for normalization, respectively [2].Note that in calculating the result with parasitics, both layout and components are considered.The post-layout simulation is performed using the component models provided by muRata.As shown in Fig. 18, the measured insertion loss at the center frequency is 1.70 dB.The 0.5-dB passband is 383-517 MHz (FBW of 29.8%).For the upper stopband starting from 674 MHz, the rejection is better than 52 dB.The measured peak reflection coefficient is −17.6 dB near the passband, and this performance is sustained up to 1.65 GHz.The measured return loss is larger than 10 dB in the range from dc to 3.04 GHz, which is roughly seven times the center frequency.The comparison between the measured responses of this filter and the referenced filter is presented in Table II.

B. Wideband Reflectionless BPF
The second synthesized result is a reflectionless BPF with a flat and wider passband.Meanwhile, multiple TZs are also required in stopbands.The constraints similar to Fig. 4 are used, while 0.5-dB FBW is set to 30%, the rejection in stopbands is set to 30 dB, and the minimum return loss of both ports is assigned to 20 dB.Besides, values of inductors and capacitors are limited to avoid the use of tiny capacitance and huge inductance in synthesis.In this case, the upper bound of inductors is set to 47 nH, and the lower bound of capacitors is set to 2 pF.After 3000 iterations, the topology is obtained.The synthesized schematic and corresponding layout are shown in Figs.19 and 20, respectively.The PCB is implemented using Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.This filter is fabricated and its photographs are presented in Fig. 21.The calculated, simulated, and measured results are shown in Fig. 22. Fig. 22 exhibits that the measured insertion loss at center frequency is 1.24 dB.The 0.5-dB passband is 372-499 MHz (FBW of 29.2%), while the 3-dB passband is 337-540 MHz (FBW of 46.3%).The stopband is from 590 MHz to 5.65 GHz, while the rejection is better

C. Reflectionless Wideband Filter With Equiripple Passband
To further verify that filters with wider passband can be synthesized, another result is a reflectionless BPF with a wider and equiripple passband.The same hyper-parameters and initial graphs used in synthesizing the previous filter are utilized.Similar boundary constraints in Fig. 4 are applied with an additional constraint of the 0.5-dB ripple in the passband.The 3-dB FBW is set to 100%.The upper bound of inductors is set to 91 nH and the lower bound of capacitors is set to 1 pF.After 3000 iterations, a suitable topology is obtained.The schematic and corresponding layout are shown in Figs.23 and 24, respectively.The PCB is implemented using the dielectric substrate of RT5880 with a thickness of 0.508 mm.The component values of the prototype and practical filter are presented in Fig. 23.
This filter is fabricated and its photographs are shown in Fig. 25.As shown in Fig. 26, the measured insertion loss at center frequency is 0.89 dB.The 3-dB passband is 247-709 MHz (FBW of 96.7%).For the upper stopband starting from 759 MHz, the rejection is better than   3.11 GHz, the measured return loss is larger than 15 dB.The peak of return loss is controlled to 10 dB by gradient-based fine-tuning, which extends the frequency range of 10-dB return loss up to 5.59 GHz.A good agreement among the calculated, simulated, and measured results is achieved.

D. Dual-Band Reflectionless BPF
To synthesize a reflectionless dual-band BPF, an efficient method is using the low-pass-to-bandpass transformation twice.Compared with synthesizing the filter directly, which needs a topology of a considerably large scale, only a small-scale prototype is required when transformations are applied.This section provides an example of synthesizing dual-band BPF using such operations.Firstly, an LPF prototype is synthesized, and the schematic is presented in Fig. 27(a).The passband insertion loss is limited to 1 dB, the stopband rejection is set to 30 dB, and a TZ is added to the stopband.The max reflection coefficient is constrained to   each capacitor is converted into the network of Fig. 27(c).After that, the schematic presented in Fig. 27(d) is obtained.Note that in Fig. 27(b) and (c), g is the prototype value of the corresponding component, ω c is the center frequency of the transformation and 1 and 2 are the FBW in the first and second transformations, respectively.1 only affects the bandwidths of the passbands, while 2 affects both the bandwidths and the center frequencies of the two passbands.If the center frequencies of the two passbands are ω 1 and ω 2 , and ω 1 < ω 2 , the value of ω c is calculated as At frequencies of ω 1 and ω 2 , the original frequency before the second low-pass-to-bandpass transformation should satisfy the equation | | = 1.Thus, 2 is derived as In this case, the center frequencies are set to 433 and 915 MHz.By using (21), center frequency f c is 629.4MHz, and ω c is calculated by 2π f c .Meanwhile, 2 is calculated from (22), and the value is 0.766.The calculated values of all components are labeled in blue in Fig. 27(d).Note that C 1 and C 5 in Fig. 27(a) can be simplified into a single capacitor.However, these components are kept in the transformation to provide more flexibility in fine-tuning after introducing parasitics.
The design layout of this reflectionless dual-band BPF is shown in Fig. 29.All values of components labeled in green in Fig. 27(d) are used in fabricating this filter, and the photographs are presented in Fig. 30.This PCB is implemented using the dielectric substrate of RO4350B with a thickness of 0.762 mm.The calculated, measured, and simulated results are shown in Fig. 31.In the measured response, the insertion loss at the center frequency of the lower passband is 1.70 dB, while it is 1.41 dB in the upper passband.The 3-dB passbands are 356-503 MHz (FBW of 34.2%) and 805-1119 MHz (FBW of 32.6%).The measured maximum peaks of insertion loss of the lower, middle, and higher stopbands near passbands are 25.2, 25.0, and 25.8 dB, respectively.The maximum reflection coefficient near passbands is −11.6 dB.

V. CONCLUSION
In this article, an adaptive synthesis procedure using lumped components is proposed and applied for synthesizing reflectionless filters.The algorithm can automatically select the suitable topology of given constraints through the specifically enhanced HGAPSO.The proposed HGAPSO is embedded with random neighborhoods and coordinates the search for improving the heuristics.Besides, a sorting-based classification is introduced for accelerating convergence, and a probabilistic method is used for boosting iteration.Components' values are fine-tuned using gradient-based optimization with parasitics of lumped components and layout being considered.To verify the effectiveness of the proposed algorithm, four reflectionless filters are synthesized and implemented.The results calculated with parasitics, the simulated results, and the measured results are matched accurately.The synthesized filters fulfill the design goals, which proves the effectiveness of the proposed methodology.

APPENDIX A FACTORS AFFECTING THE CONVERGENCE
OF THE PROPOSED HGAPSO Under a properly-configured hyperparameter, the proposed procedure is adaptive to select suitable topologies for synthesizing circuits with different response constraints.In this section, several factors affecting the convergence of the proposed HGAPSO are investigated by comparison.Each comparison is based on the setup in synthesizing the reflectionless BPF shown in Fig. 19.Note that the reflection coefficient is constrained to −15 dB for faster convergence during the benchmark.Every convergence curve is the median result of all tests.
All hyper-parameters configured in synthesizing the flat-passband BPF are used as a reference in testing the effect of different hyper-parameters.The PSO parameters w g and w l are fixed to 0.024 and 0.016, respectively.Fig. 32 shows the effect of hyper-parameters introduced by the proposed HGAPSO.The convergence curves of different total number of particles (P) are shown in Fig. 32(a).It can be seen that more particles provide greater search capability, but the time consumption of each iteration is increased.The lifespan of each particle (L init ) is set proportionally related to parameter P. Since the lifespan of a particle is reduced by its rank in each iteration, the factor of L init to P is the minimum searching time of each particle.As shown in Fig. 32(b), L init = 120• P is used to achieve the high particle efficiency during the convergence.In Fig. 32(c), the optimal weight factor w r locates in the range of 0.6-0.7,which provides the proper searching range for random neighborhood search.The convergence curves of different values for the parameter G are exhibited in Fig. 32(d).This parameter defines the minimum gap between the upper bound and the lower bound in each dimension of the classifier.Setting G between 0.3 and 0.6 would now significantly affect the convergence.The population of optimal particles K is set proportionally related to parameter P. Using particles of the top 2% to update the classifier and define the optimal zone can obtain the best convergence performance, as shown in Fig. 32(e).According to the results in Fig. 32(f), the best probability of regenerating new particles in the optimal zone is around 0.7.
Meanwhile, constraints are factors that affect the convergence of the synthesis.Fig. 33(a) shows the effect of limiting the value boundaries (i.e., V min and V max ).Limiting the range of components can avoid extreme values and reduce the performance degradation caused by the parasitics.However, once prototype values are limited to a very narrow range, the convergence is slow.Another factor is the constraint on the reflection coefficient.As shown in Fig. 33(b), the convergence is slower under tighter constraints.Fig. 33(c) shows the effect of the graph scale on the convergence.In Fig. 33(c), E(G 0 ) and C(G 0 ) denote the number of edges (i.e., the number of branches) and the number of components of the graph G 0 , respectively.If the scale of the graph used in the synthesis is larger than the graph that it actually needs, the synthesis requires more time to find a suitable topology.But a graph of a larger scale has more subgraphs, which ensures more circuits be obtained in the synthesis.
In short, setting hyper-parameters to unsuitable values will result in slower convergence.However, as shown in Fig. 32, slightly deviated parameters do not deteriorate the convergence significantly.If the synthesis constraints are changed while hyper-parameters are not configured to optimal values, suitable results can still be obtained after more iterations.

APPENDIX B TIME CONSUMING AND THE ACCELERATION OF PARALLELIZATION
The proposed HGAPSO is implemented in C programming language, with a simple built-in thread pool for multithreading.The platform is built by an Intel-D1581 with Ubuntu 20.04, and the program is compiled by GCC 10.3.0.The synthesis of the flat-passband BPF (Fig. 19) is used as a benchmark.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.The test result is shown in Fig. 34.The speedup ratio S p and efficiency E p are calculated by following: where T p is the average time consumption of each iteration using p threads.Because the cost value of each particle can be calculated individually, the speedup ratio of processing particles should be ideally linear.However, due to the additional process of submitting thread work, thread-switching in the system, and nonparallelized parts (i.e., sorting the cost value, calculating life span, and reinitializing particles), the speedup ratio is below the ideal value.Therefore, the efficiency of acceleration gradually decreases along with the number of threads.Because the platform only has 16 CPU cores, the average time consumption is slightly longer when using 16 threads compared with 15 threads.Even so, the algorithm is benefited from parallelization.Besides, when 16 threads are used, the program dominates 15 MB of RAM in total (from the "RES" column in the "htop" processes viewer).In this test, the speedup ratio can reach up to 13 when 15 cores are used.

Fig. 2 .
Fig. 2. (a) Complete graph K 6 that is used as an example of G 0 .(b) Circuits model for each branch in a graph.

Fig. 3 .
Fig. 3. Example of a quasi-planar graph that has 12 nodes aligned in a grid and an extra node utilized as the ground node.(The subgraphs connected by branches drawn in solid lines are planar graphs.Branches drawn in dashed lines are ground-connecting branches.) [A c ] has no effect on the result.In matrix [Y b ], all elements are zero except diagonal ones.[Y b ] has the following form:

Fig. 4 .
Fig. 4. Example of bandpass response and its constraints described by the target function.

Fig. 5 .
Fig. 5. Examples of two different sampling methods in different situations.(a) Average error of all sampling points dominates the early stage of optimization.(b) Maximum error of all sampling points dominates the late stage of optimization.

Fig. 8 .
Fig.8.Example for illustrating the effectiveness of the embedded RCS.(The range of X i → 0 and X j → 0 denotes a shorted resistor, a shorted inductor, or an opened capacitor, while the range of X i → ∞ and X j → ∞ denotes an opened resistor, an opened inductor, or a shorted capacitor.)

Fig. 9 .
Fig. 9. Distributions of the constructed probability functions (k 1 and k 2 are set to 5).

Fig. 10 .
Fig. 10.Example showing one step iteration of five particles in Ackley function of two dimensions.In this example, the classifier is used to cover the best two particles.

Fig. 11 .
Fig. 11.Convergence comparison of different algorithms using target functions in synthesizing reflectionless filter.(a) Target of a small-scale reflectionless BPF.(b) Target of a reflectionless wideband BPF.

Fig.
Fig. Schematic of the synthesized reflectionless BPF with symmetrical topology.

Fig. 15 .
Fig. 15.Calculated response of the synthesized filter and the three-stage BPF from [4].

Fig. 18 .
Fig. 18.Calculated, simulated, and measured results of the synthesized reflectionless BPF with symmetrical topology.(a) Response with frequency normalization using the center frequency of 433 MHz.(b) Response without frequency normalization.

Fig. 19 .
Fig. 19.Schematic of the synthesized reflectionless wideband BPF with flat-passband.(The middle line of each label, i.e., the blue line, is the corresponding prototype value.The bottom line of each label, i.e., the green line, is the value of the component used in fabricating the filter.)TABLEII COMPARISON BETWEEN MEASURED RESPONSES OF LUMPED BPFS

Fig. 22 .
Fig. 22. Calculated, simulated, and measured results of the synthesized reflectionless wideband BPF.(a) Response with frequency normalization using the center frequency of 433 MHz.(b) Response without frequency normalization.than 30 dB.The measured reflection coefficient is lower than −15 dB from dc to 4.55 GHz, which is more than ten times the center frequency.The calculated, simulated, and measured results show good agreement.

Fig. 26 .
Fig. 26.Calculated, simulated, and measured results of the synthesized reflectionless wideband BPF with equiripple passband.(a) Response with frequency normalization using the center frequency of 433 MHz.(b) Response without frequency normalization.

Fig. 27 .
Fig. 27.Schematic of the synthesized LPF prototype, the low-pass-to-dual-band transformation, and the obtained dual-band BPF.(a) Schematic of the reflectionless LPF prototype.(b) Network converted from an inductor in the transformation.(c) Network converted from a capacitor in the transformation.(d) Schematic of the reflectionless dual-band BPF transformed from the prototype.

Fig. 31 .
Fig. 31.Calculated, simulated, and measured results of the reflectionless dual-band BPF.(a) Response with frequency normalization using the center frequency of 629.4 MHz.(b) Response without frequency normalization.

Fig. 32 .
Fig. 32.Comparisons of the effect of different hyperparameters on the convergence.(a) Number of total particles P. (b) Factor of the lifespan of each particle L init to P. (c) Factor w r of the max step of random neighborhood search to the range of the classifier in each dimension.(d) Value G of the minimum gap between the upper and the lower bound of the classifier in each dimension.(e) Factor of K -P, where K is the number of optimal particles in generating the optimal zone.(f) Probability p of generating reinitialized particles in the optimal zone.

Fig. 33 .
Fig. 33.Convergence curves of different constraints in synthesis.(a) Different limits on the range of prototype value.(b) Different constraints on reflection coefficient.(c) Initial graphs of different scales.

Fig. 34 .
Fig. 34.Acceleration obtained by multithreading.(a) Average time consumption of each iteration versus the number of threads.(b) Speedup ratio and efficiency versus the number of threads.

TABLE I LIST
OF HYPER-PARAMETERS USED IN THE TEST significantly.Besides, during the preoptimization, if f ( ⃗ X i