Mosaic Based Optimization of NRD Guide Devices Using Binary Evolutionary Approaches and 2D-FVFEM

In this paper, efficient optimal design approaches based on mosaic optimization concept are developed for NRD guide devices to realize the high-performance compact millimeter-wave integrated circuit. Binary representation-based genetic algorithm, differential evolution algorithm, harmony search algorithm, firefly algorithm, and particle swarm optimization are developed to efficiently optimize the pixel pattern in the design region of NRD guide devices. To demonstrate the usefulness of these population-based optimizations, a comparative study based on the problem-solving success rate is conducted first. To carry out this study, four NRD circuit components are designed which include low crosstalk waveguide crossing, T-branch power splitter, bending waveguide, and frequency demultiplexer. The proposed optimal devices achieve high transmission efficiencies greater than 99.9%, 49.9%:49.9%, 99.9% at 60 GHz and 96.4%, 98.5% at 59 GHz and 61 GHz. In addition, the same NRD guide components except frequency demultiplexer are also designed at wideband operation and achieve broad bandwidth around 5 GHz, 4 GHz, and 3 GHz. In order to improve the computational efficiency, the originally developed two-dimensional full vectorial finite element method is employed for the numerical simulations. This paper demonstrates the detailed implementation procedure of developed evolutionary approaches for the material distribution in the design region of NRD guide devices, comparative study of developed optimization approaches, and proposed highly attributed NRD circuit components for the realization of NRD based high-performance compact millimeter-wave circuit.


I. INTRODUCTION
T HE use of millimeter-and terahertz-wave bands is being actively explored to increase communication system capacity and meet modern communication requirements. Non-radiative dielectric (NRD) waveguide technology has received a lot of attention in recent years due to non-radiative and low loss nature. Several NRD guide components have been reported so far for microwave and millimeter-wave circuit system applications. At 50 GHz, the dispersion characteristics, transmission loss, coupling coefficient, and measurement setup with several NRD-based circuit components, including T-junction, right angle corner, 90°and 180°bends, and directional coupler, are reported [1]. An NRD-based Tjunction with a dielectric stub and thin metal patches with an output power of around -4 dB at 35 GHz has been proposed [2]. The coupling theory was used to analyze the losses in NRD bends (90°and 180°) at 50 GHz, with bending loss less than 0.3 dB has been presented [3]. Following a detailed examination of operational principles, propagating modes, losses, behavior of transmitting waves in bending structures, practical significance, and the confirmation of NRD as a lowloss guide in [1]- [3], the authors proposed several NRD circuit elements at an operating frequency of 35 GHz, including a matched terminator, directional coupler, circulator, beam lead diode, gun diode oscillator, transmitter, and receiver. At 50 GHz, an NRD guide filter, a four-way power divider, and an NRD leaky-wave antenna were fabricated [4].
In millimeter-wave and terahertz applications, a substrate integrated non-radiative dielectric waveguide (SINRD) for printed circuit boards has been developed [5]. The hybrid planar non-radiative waveguide is presented and confirmed as a building block for millimeter-wave circuits, with an experimental prototype based on hybrid technology that includes active and passive components. This novel technique will be greatly useful in the development of future millimeter-wave integrated circuit systems [6]. Using mode coupling theory, a technique for designing circular and racetrack-shaped NRD guiding ring resonators was successfully devised, and a low loss, small-sized ring resonator with radii less than 3.5 mm was built at 60 GHz with band rejection performance of more than 30 dB [7]. At a frequency of 77 GHz, a new sort of NRDbased directional coupler using two separate NRD guides interconnected with a bridge is presented and fabricated [8]. At 60 GHz, a first tunable liquid crystal filter based on NRD technology is presented, with fractional bandwidth and tunability of 1% and 2.5%, respectively [9]. Furthermore, over the tuning range, the filter insertion loss ranges between 4.9 and 6.2 dB.
Many additional research for the development of NRD waveguide theory and circuit components for millimeter application were presented in the literature. Some of them have already been addressed above, while the rest are listed in Table 1. [10]- [33]. In this comprehensive evaluation, all NRD circuit components are proposed without employing optimization approaches for NRD guiding structures. It is well understood that without optimization, a highperformance compact millimeter-wave integrated circuit is not viable. For this purpose, we developed several optimization approaches to optimize the design region of NRD guide devices. In this regard, our proposed research work is unique and distinct from the existing literature on NRD waveguide technology.
Several optimization approaches and variants based on size, shape, topology, and mosaic optimization concepts, such as genetic algorithms, particle swarm optimization, inverse design algorithms, and direct binary search algorithms, have been proposed for the efficient design of optical, dielectric, and photonic crystal guide devices [31]- [50]. However, a study about the development of optimal design approaches using the above-mentioned optimization concepts for NRD guides is insufficient and still has great potential in this domain. In the literature and to the authors best knowledge the development of binary evolutionary approaches based on mosaic optimization concepts, detailed implementation procedure, comparative study, and high performance NRD circuit components at single frequency operation and broadband operation, has not been proposed previously and it is also verified from Table 1.
In this paper, in order to add the valuable contributions in NRD platform, we develop several optimization approaches for the efficient design of NRD guide components with various functions to realize the high-performance compact millimeter-wave circuit as depicted in Fig. 1. For this purpose, digital material concept (mosaic-optimization) using binary evolutionary approaches is very useful. Using mosaicbased optimization a small number of variables are used to express the design region of NRD guide and those variables are optimized by five optimization approaches which include binary-representation-based genetic algorithm (GA), differential evolution algorithm (DEA), harmony search algorithm (HSA), firefly algorithm (FA), and particle swarm optimization (PSO). To find the most efficient ones, a comparative study based on problem-solving success rate (also known as convergence rate) is carried out first, using four numerical examples: low crosstalk waveguide crossing, T-branch power splitter, and bending waveguide at 60 GHz, as well as frequency demultiplexer at 59 GHz and 61 GHz. In addition, same NRD guide components except frequency demultiplexer are also designed by modifying the objective function for broadband operation. In modified objective function we used several frequency intervals within desired frequency band for each generation of optimization method. On other hand, to improve the design efficiency originally developed two-dimensional full-vectorial finite element method (2D-FVFEM) is used for the numerical simulation of NRD guide devices [51]. Our 2D-FVFEM is verified by 3D-FVFEM using several optimal and non-optimal NRD guide structures [18]- [20].
The paper is organized as follow: The brief discussion of NRD guide and structure representation using mosaic-based optimization are presented in section II. The evolutionary approaches for binary optimization and the brief discussion of BGA, BDEA, BHSA, BFA, and BPSO are described in section III. The simulation results of numerical design examples of NRD guides are described in section IV. Conclusion of this research work is in section V. This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.  Microstrip transition -Bandpass filter Planar NRD oscillator [7] 60 Ring resonator - [8] 77 Directional coupler - [9] 60 Liquid crystal filter - [10] 60 Mixer using bandpass filter - [11] 60 Flexible transmission line - [12] 60 Beam lead diode - [13] 77 Amplitude shift keying -Phase shift keying [14] 60 Transmitter and receiver - [15] 60 ASK Transceiver - [16] 60 High speed transceiver - [17] 60 Guide duplexer - [18] 60 90°bend DBSA

A. NRD GUIDE
NRD guide has been proposed as a millimeter waveguide by T. Yoneyama [1]. As illustrated in Fig. 2(a), its structure is made up of dielectric sandwiched between metal parallel plates with spacing less than λ/2 to realize the non-radiative nature. In NRD guide, propagation loss is relatively small because most of the electromagnetic energy flows through the dielectric and air and a very small amount of current flow in the metal plates. Usually, NRD guide support two orthogonal modes such as LSM 01 and LSE 01 mode. For LSM 01 , the electric field is parallel to the metal plate whereas in LSE 01 mode, the electric field is perpendicular to the plate as shown in Fig. 2(b). Both modes are non-radiative in nature and cannot propagate within air region sandwiched metal parallel plates and can be used as a non-radiative waveguide. LSM 01 mode whose polarization is parallel to the metal is usually used in NRD guide, in contrasted to microstrip line which supports polarization perpendicular to the metal. Furthermore, the dispersion relation of NRD guide is also calculated as shown in Fig. 2(c).

B. STRUCTURE REPRESENTATION OF NRD GUIDE USING MOSAIC-BASED OPTIMIZATION
Achieving the desired properties typically requires many design variables to optimize material distribution in the target area of the guide device. Several optimization approaches have been proposed to optimize these variables, using various optimization concepts such as shape, size, and topology optimization. In which, size optimization deals with dimensional parameters such as the length and width of the structure, whereas shape optimization optimizes the outer shape of the structure, and the topology optimization optimizes the material distribution within area of interest. On other hand, mosaic-based optimization is another way to efficiently optimize the guide devices. Mosaic structures show the presence of dielectric materials in the form of pixels at particular positions within the area of interest.
In this work, mosaic-based optimization using evolutionary approaches is considered to optimize the material distribution in the design region of NRD guide. Several design variables are employed to express material distribution in the design region. In order to create a mosaic-like structure, first discretize the design region into a grid pattern and then allocate a dielectric or air to each pixel (also known as a mosaic) in the form of 1 or 0, respectively as shown in Fig. 3. This is called the concept of digital materials. The degree of design freedom of our optimization approaches is determined by the number of pixels. A large number of pixels may make it difficult to fabricate a designed device. The pixel size is determined not by the convenience of the optimization method, but by the possibility of fabrication. A brief description of implementing the optimization approaches is given below in the next section.

III. EVOLUTIONARY APPROACHES FOR BINARY OPTIMIZATION
In evolutionary approaches, three processes are involved in the evolution of the population. The first step is initialization, which involves randomly generating an initial population of individuals based on a predetermined quantity. Determining an appropriate population size and individual representation technique that are more likely to offer feasible solutions are important. Another important parameter is the number of iterations. Both the population size and the number of iterations may have a significant impact on solution quality and time. After a few trials, an appropriate selection of these parameters can be determined based on the design problem resulting in a feasible solution. Individuals can be represented using a variety of ways, including binary representation, real value representation, and integer value representation.
The second step is to calculate the value of each individual using a problem-dependent objective function after the population has been randomly generated. The calculated value can be used to rank each individual for selection purposes. Whether the breaking criteria has been satisfied is judged after the evaluation of objective function. Different breaking criteria, such as static, dynamic, and hybrid, are used to stop the optimization process. When dynamic criteria is used, the optimization is executed until the desired characteristics is achieved, whereas the optimization is terminated at a fixed number of iterations in the case of static criteria. A combination of these two criteria is called a hybrid and is used in some cases. If no individual satisfies the desired property, the third stage of the evolutionary technique generates a new population by perturbation of solutions in the existing population. In evolutionary approaches, three processes are involved in the evolution of the population. The generic optimization procedure of our binary evolutionary approaches that we considered in this paper is shown in Fig. 4. In this paper, the size of population (Np = 64) is same for all evolutionary approaches. The starting search point is varies from one another due to the random generation of populations, resulting in differing optimal structures. The number of iterations required may be determined by the device function, whether it is designed for single-frequency operation or wideband operation. For single and broadband operation, we evaluated a maximum of 100 and 300 iterations for all devices, respectively. Binary representation is considered to represent an individuals structure because the digital material concept is employed to generate a mosaic-like structure. The optimization process is terminated by using static breaking criteria.

A. BINARY GENETIC ALGORITHM (BGA)
John Holland proposed GA, a metaheuristic global searchoptimization method, in 1975 [52]. Conventional GA and its variants are used to solve a wide range of complex optimization problems. The basic goal of GA is to simulate natural selection and survival of the fittest individuals. In GA, the procedure starts with a randomly generated initial population. Each individual in population is a solution set or chromosome, which is encoded using digital material techniques. After the population is randomly generated, each individual is evaluated using a problem-dependent objective function and ordered from best to worst based on the calculated value. GA is a method for generating new populations that mimics natural selection of individuals, and it employs three biologically inspired operators: selection, crossover, and mutation. The developed optimization method considers two different selection techniques. To save the best individual from the entire population, first elite selection is used. On the other hand, rank selection, is used to select the best individuals for a new generation. Once the parent individuals x (n) p1 and x (n) p2 are selected, the crossover operator perturbs the selected individuals to produce new offspring as follows: where U (0, 1) is a random number between 0 and 1 and x i . After crossover, the newly generated offspring have a tendency to become highly similar to previously selected individuals, reducing population diversity and perhaps leading to population stagnation. The bit flip mutation operator is employed to inject diversity into the population to avoid stagnation. With this optimization approach, the mutation rate is 1%. The evolution of a randomly generated population is now complete. This search process is repeated until the number of iterations is reached. Finally, the best individual at the final iteration will be the optimal mosaic-like structure of the NRD guide.

B. BINARY DIFFERENTIAL EVOLUTION ALGORITHM (BDEA)
DEA is a metaheuristic optimization approach, and first proposed by Storn and price in 1995 for global optimization [53]. Usually, DEA performs better than GA because it explores the given search space more efficiently, when multi objectives are need to be optimized. The key reasons for differential evolution application in electromagnetic optimization are its reliability in dealing with multi-minima functional and its improved rate of convergence as compared to GA when applied to small scale real valued problem. It is critical to improve the candidate solution with each iteration, and it doesn't require explicit cost function gradients in the candidate solution search space. As a result, DEA can be applied to problems that are not even continuous or nondifferentiable. The most important feature of DEA is its unique mutation operation. The scaling parameter S and the population size N P are two critical control parameters in DEA. The key steps in the implementation of DEA are mutation, crossover, and selection. First, in binary strategy, a mutant individual is generated as follows: m,k is a mutant individual and may possible to take the value of other than binary number -1 or 2 if the values of randomly selected individuals are (x p2,k , and x (n) p3,k ) = (0, 0, 1) or (1, 1, 0) respectively. In order to binarize the mutant individual the value of x (n) m,k is set to be (0, 1) as follow.
The scaling parameter S is used to control the differential variations and it is set to be 0.75. The operation described above is called the mutation of DEA. The second step of DEA is crossover. In DEA, x t is the candidate of next generation by crossover between mutant individual x (n) is restored. The third step of DEA is to select the best individual between target individual and candidate individual for next round until the completion of iteration.

C. BINARY HARMONY SEARCH ALGORITHM (BHSA)
HSA is a population-based metaheuristic optimization approach proposed by Geem et al. in 2001 [54]. HSA is inspired by the improvisation process of experienced musicians. The key benefits of HSA are its clarity of execution, record of success, and ability to solve a variety of complex problems. Pitch adjustment rate, bandwidth, and harmony memory considering rate are the three main parameters that control exploitation and exploration in HSA. Through the selecting and tuning process, a new harmony is created in HSA. Each design variable is chosen from the harmony memory, which contains previous good harmonies. A random number, with a specified probability, tunes or replaces the selected harmony. The new harmony x new h is generated as follows: where C h = U (0, 1) and x m,k is k-th sound generated in considering harmony memory. R hmcr is the harmony memory considering rate (HMCR), and it is utilized to make the harmony memory more effective. If the HMCR value is very low, close to 0, only a few harmonies are used, and therefore the convergence rate is slow, whereas large rates exploit the harmonies, and the solution space is not fully explored, thus leading to an insufficient solution. The value of HMCR should be in the range of 0.7-0.95. VOLUME 4, 2016 In this optimization approach the R hmcr = 0.85 has been set. The pitch adjustment rate R par , is used to control the degree of adjustment. Because of the limitation in exploring only a small subspace of the entire subspace, a small R par with narrow bandwidth can slow down HSA convergence, whereas a large R par with a wide bandwidth can lead to optimization deviation in some optimal solutions. In this optimization approach, the pitch adjustment rate is set to be R par = 0.2.

D. BINARY PARTICLE SWARM OPTIMIZATION (BPSO)
PSO is an iterative optimization approach that takes its inspiration from the animal kingdom, specifically a group of animals moving in search space for common objective. Dr. Eberhart and Dr. Kennedy were the first to suggest PSO [55]. In PSO, a swarm of particles is the first to be randomly initialized. Each particle has its own unique position and velocity. Each particle changes its position continuously inside a specified search space with velocity to find the optimal location. The particles movement is determined by its memory path and iteration with other particles. The best position of each particle and the best position of the swarm are updated once all particles have reached the new place. PSO finds the best solution by collaborating and sharing information among swarm particles. The following are the steps in the PSO procedure: 1. Generate an initial swarm with initial position and velocity. 2. Evaluate objective function of each particle. 3. Update both the individual and global best positions. 4. Update each particle velocity and position.
The optimization process is repeated until the termination criterion is met. The following equations are used in PSO to update the position and velocity of the particles.
where y are the position and velocity of i-th particle respectively, while y (n) p and y (n) g are the individual and global best positions of n-th iteration respectively. W is the inertial weight parameter with a value 0.8 that adjusts the ability to search locally and globally. The initial cognitive and social cognitive coefficients, C 1 and C 2 can adjust the balance between individual and global optimum positions. r 1 and r 2 are two random values between 0 and 1 which are used to simulate the random components of swarm behavior.
To obtain x (n+1) i , the following binarization scheme is employed: where ∆ is a coefficient introduced to avoid being trapped in a local solution and is set to be ∆ = 0.125.

E. BINARY FIREFLY ALGORITHM (BFA)
FA is a swarm-based metaheuristic optimization algorithm first proposed by Xin She Yang in 2008 [56]. FA is inspired by the movement of fireflies as they interact based on their flashing light. Regardless of their gender, particles are attracted to others because of their attractiveness. Attractiveness is determined by brightness and distance between them. The optimization process of FA consists of three steps. First, a particles called fireflies are randomly generated. Each firefly has a brightness level. The brightness is used to calculate the value of each firefly by using objective function. For minimization problem, the fireflies with small value are more brighter than higher ones. Once brightness of fireflies is evaluated, then fireflies follow those ones who are rich in brightness. The brightest firefly conducts a local search by moving randomly in its surroundings. The position of firefly y (n) i is updated by equation as follow: The attraction to the brighter neighbor is calculated using the second term of the above equation where, β (0,j) is attractive force of the j-th particle and r ij is a distance between i-th and j-th particles. u(ξ) is the unit step function. y The light absorption coefficient of medium is γ = 1/ √ L and L is used to adjust the search range of particles. The third term in the above equation is for randomization, where α is the scaling factor used to adjust the random step length, δ is the damping factor and ε is a random vector. The values of δ and α are 0.99 and 0.25, respectively. Using the above equation, the firefly's position is updated before the end of iterations.

IV. DESIGN EXAMPLES OF NRD GUIDE
In this paper, we proposed the concept of realizing a highperformance millimeter-wave integrated circuit using basic NRD guide components as shown in Fig. 1. Evolutionary approaches have been developed to obtain the optimal design of the NRD guides. To demonstrate the usefulness of these optimization approaches, we considered four numerical examples of NRD guide devices. Except symmetrical conditions in design region, size of design region, position and number of output ports, the other basic geometrical parameters are the same in all devices. In our design examples, we assume the spacing between metal plates is a = 2.2 mm, length and width of the dielectric strip are l = 10 mm and b = 2 mm, respectively. The relative permittivity of dielectric and air is ε r = 2.2 and ε air = 1.0 respectively. The size of design region in waveguide crossing, T-branch power splitter and 6 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3180732 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ bending waveguide is W x = W y = 6 mm and it is discretized into 144 pixels. In case of frequency demultiplexer W x = 6 mm and W y = 16 mm with total 384 pixels in the design region. The size of each pixel is 0.5 mm × 0.5 mm that is quarter of the dielectric strip. In order to make the fabrication feasible, pixel size is determined by considering the operating wavelength, guide width and fabrication sensitivity. The computational domain is surrounded by a perfectly matched layer (PML) with a thickness of d PML = 5 mm. LSM 01 mode at 60 GHz is incident at input port 1 in first three examples. In case of frequency demultiplexer wave at two frequencies 59 GHz and 61 GHz are incident, and it is selectively output to port 2 and 3 respectively. In order to achieve the desired properties, the design region of size W x × W y is optimized.

A. LOW CROSSTALK WAVEGUIDE CROSSING
First, we consider low crosstalk NRD waveguide crossing as shown in Fig. 5. The geometrical parameters and incident conditions are explained above. In order to make the optimization process easy and faster, a unique one-two structural symmetric conditions are applied in the design region indicated by red line, to optimize the half design region except whole as shown in Fig. 5. The following objective function is used at 60 GHz to achieve the maximum transmission efficiency at output port 3 when LSM 01 mode is incident at input port 1.
The pixel pattern is optimized by using five different optimization that are described in section III. The convergence behavior of these optimization approaches at 60 GHz is shown in Fig. 6(a). We can see that all optimization approaches show good convergence in the range of 60 to 80 number of iterations. The performance of BHSA and BDEA is slightly better, but other optimizations also achieved sufficient device performances. Figure 7 shows the optimized structure, propagation field, and frequency characteristics of BHSA which is the best one in case of waveguide crossing at 60 GHz operation. The frequency characteristics of waveguide crossing using other optimizations are shown in Fig. 8. Due to random generation of initial population in each optimization the optimal structures are different to each other. Although, optimized structures are different but overall device performance is almost same for each design as shown in Fig. 8. All optimized waveguide crossings achieve high transmission power more than 0.95. In Fig. 7, frequency analysis by BHSA relatively show broad bandwidth characteristics. By considering this clue, we modified our objective function to realize the broadband operation possible. In order to achieve the broader bandwidth we used the following modified objective function. The updated objective function analyzes the waveguide structure using three different frequencies at a time for each generation of optimization approaches. Usually, a large number of iterations are required to optimize the structure for broadband operation because multiple frequencies involve in the objective function. In Fig. 6(b) we can see that how efficient our developed optimization approaches that converge the solution at broadband operation with almost same number of iterations used for single frequency operation. The optimized structure at broadband operation, propagation fields at considered frequencies in the objective function and frequency characteristic by BDEA are shown in Fig. 9. Furthermore, the performance detail of proposed NRD waveguide crossing is summarized in Table 2. By using modified objective function, all optimization approaches achieve high broadband property around 5 GHz in the frequency range of 58 GHz-63 GHz.

B. T-BRANCH POWER SPLITTER
Next, we considered second numerical example NRD Tbranch power splitter as shown in Fig. 10. The geometric parameters are the same as in the previous example, except for the number of output ports and the condition for structural symmetry. In order to split the transmission power equally at output ports, one-two symmetric condition in rectangle form is applied indicated by red line as shown in Fig. 10. Due to considered symmetry, pixels in the upper half of the design region are optimized. The following objective function is used to achieve the maximum transmission power at 60 GHz.
This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication.       where |S 21 | 2 and |S 41 | 2 are transmission powers from port 1 to port 2 and port 1 to port 4 respectively. Figure 11  to 60 number of iterations. However, BGA requires more number of iterations for better convergence in comparison to BDEA, BHSA and BFA but it also achieves a desired device performance. The optimal structure, propagating field at 60 GHz and frequency characteristics of BHSA obtained by 2D-FVFEM are shown in Fig. 12, and it depicts that optimal structure can efficiently split the power equally into two output ports. The frequency characteristics of T-branch guide using other optimizations are shown in Fig. 13 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3180732 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ all optimization approaches for broadband operation by using objective function as follow: (f 1,2,3,4,5 = 58, 59, 60, 61, 62 GHz) In order to ensure the broad bandwidth, more close intervals within interested frequency band are considered in the objective function. Figure 11(b) show the convergence behavior at broadband operation and it shows almost same behavior and convergence order as in 60 GHz operation. In this design, BDEA, BHSA, and BGA achieved desired property around 25, 60 and 250 number of iterations. Due to structural limitations of T-branch guide, it is quite challenging for each optimization approach to achieve a desired property. Therefore, the performance of BPSO and BFA is not satisfactory. Successful optimizations achieved same device performance with broad bandwidth around 4 GHz in the frequency range of 58 GHz-62 GHz, but BHSA and BDEA are more efficient than BGA. The best optimal results obtained by BDEA are shown in Fig. 14. Furthermore, the performance detail of proposed NRD T-branch power splitter is summarized in Table 3.

C. BENDING WAVEGUIDE
Next, we considered the third numerical example of an NRD bending waveguide as shown in Fig. 15. Except number of output ports, geometrical parameters and structural symmetrical conditions are the same as in the example of waveguide crossing. For maximum transmission, the following objective function is used at a single operating frequency.
where |S 21 | 2 is the transmission power at port 2. Figure 16 (a) shows the convergence analysis of optimization approaches at 60 GHz in which BDEA, BHSA, BGA, and BFA achieving the best convergence results not more than 60 iterations but VOLUME 4, 2016 as usual GA require more iterations confirmed in the previous example in both operation. Because of the sharp bending structure, it is challenging for each optimization approach to satisfy the desired properties. Figure 17 shows the optimal results of a bending waveguide at 60 GHz using BHSA. The frequency characteristics of bending waveguide using other optimizations are also shown in Fig. 18. As shown in Fig. 17   For broadband operation, we used the objective function as follow. The above objective function analyzes the guide structure using three frequencies 59, 60, 61 GHz. Figure 16(b) depicts the convergence behavior of a bending waveguide at broadband operation, with almost same convergence order to that of the preceding example of a broadband T-branch, where BDEA and BHSA are more efficient than BGA. BDEA, BHSA, and BGA achieved broad bandwidth about 3 GHz in this design example, while BPSO and BFA achieved 2.3 GHz and 2 GHz respectively. The best optimal results at broadband operation using BDEA are shown in Fig. 19. Furthermore, the performance detail of proposed NRD bending waveguide is summarized in Table 4.

D. FREQUENCY DEMULTIPLEXER
Finally, we consider another interesting example of NRD frequency demultiplexer whose initial structure, placement of output ports, and incident conditions are different from those in the prior examples as shown in Fig. 20. It has one input and two output ports, as well as a design region of size 6 mm × 16 mm between them. The design region is divided 10 VOLUME 4, 2016 This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and  into 12 × 32 pixels, with each of size is 0.5 mm × 0.5 mm.
The gap between two output waveguides is set to be d = 8 mm to avoid the coupling effect. To get the maximum power of a particular frequency at a desired output port, we ignore structural symmetry conditions in this case. Two frequencies f 1 = 59 GHz and f 2 = 61 GHz with LSM 01 mode are incident at input port 1 and are separated into two different output ports. The following objective function is used to separate the input frequencies.
Minimize  2 and port 3 for 59 GHz and 61 GHz respectively. After going through a design region optimized by using binary evolutionary approaches, two input frequencies are divided into two different output ports. The convergence trend of optimizations for designed device at 59 GHz and 61 GHz is shown in Fig. 21. In comparison to the preceding devices, more iterations are required for optimum convergence to obtain the desired properties due to complicated device functionality, multiple frequencies, and large size of design region. As demonstrated in Fig. 21, all optimization algorithms have a good enough convergence rate of less than 300 iterations. Figure 22 shows the best optimal results of the NRD frequency demultiplexer device using BFA. The frequency characteristics of the designed device is shown in Fig. 22 (a), with maximum transmission power at 59 GHz and 61 GHz, as indicated by the straight magenta line. At 59 GHz and 61 GHz, the transmission power is 0.964 and 0.985 with minimal crosstalk level 0.03 and 0.006 respectively. Figure  22(b) and (c) indicate that the propagation fields of 59 GHz and 61 GHz are fully guided into output waveguides without any coupling effect. The main purpose of the demultiplexer is to achieve multi-path interference in the design region. For 59 GHz operation, constructive interference occurs at output port 2 and destructive interference occurs at output port 3, and vice versa for 61 GHz operation. Figure 23 shows the frequency characteristics analysis of designed device utilizing various optimization approaches, which similarly obtained the desired transmission property of more than 0.95. Furthermore, the performance detail of proposed NRD frequency demultiplexer is summarized in Table 5.
The developed evolutionary approaches have a large degree of design freedom and require a small number of variables to efficiently optimize the pixel pattern in the design region. They result in high computation efficiency and converge the solution very quickly within a few iterations. This has been observed in proposed NRD guide devices. In addition, no matter how complex the device functionality and structure is or how large the design region is, developed evolutionary techniques are totally trustworthy for any type of NRD guide at single frequency operation. In spite of several practical advantages that are mentioned above, developed evolutionary approaches have some limitations too, the most important one being the careful selection of hyper parameters. Inappropriate parameter selection can result in population stagnation and premature convergence, rendering a solution unfeasible. Two of our five developed optimization approaches, BPSO and BFA, are not completely reliable in broadband operation, which is a only limitation. On the other hand, BDEA, BHSA, and BGA, are extremely efficient in both operations.  In order to investigate the stability performance of our developed evolutionary approaches, we considered the NRD T-branch power splitter which is the most challenging device for developed optimizations has been proven above. The device configuration, initialization and other hyper parameters of each optimization algorithm are same. The stability analysis is conducted by optimizing the considered device through ten attempts as shown in Fig. 24. We can see that BHSA and BDEA proves its high stability characteristics in each attempt but BHSA is more efficient than BDEA. The value of objective function in BGA is less than 0.2 with 100 iterations and may possible to achieve almost 0 with more iteration as obtained in BDEA and BHSA. The stability performance of BPSO and BFA is not impressive as achieved in rest of the optimizations. However, we can consider both of them for other simple NRD guide devices such as bending and crossing waveguide. Overall this analysis state that the BHSA is best one in term of stability, computational efficiency and reliability. In this paper, we developed a set of evolutionary approaches using digital material concept for the efficient optimization of NRD guide devices. To reduce the computational efforts, originally developed 2D-FVFEM is used for the numerical simulations. Binary representation based GA, DEA, HSA, FA, and PSO are developed to efficiently optimize the material distribution in the design region of NRD guide devices. Four NRD circuit components including low crosstalk waveguide crossing, T-branch power splitter, bending waveguide and frequency demultiplexer are designed, and achieved high transmission throughput greater than 99.9%, 49.9%:49.9%, 99.9% at 60 GHz and 96.4%, 98.5% at 59 GHz and 61 GHz respectively. Furthermore, same devices except frequency demultiplexer are designed at broadband operation using modified objective function. Through comparative study we found that developed evolutionary approaches are highly efficient for single frequency operation. However, it is quite challenging for broadband operation because of several frequencies in the objective function and structural limitation. We have achieved quite impressive results from proposed devices with minimal reflection and crosstalk level due to highly efficient optimization approaches. The proposed NRD guide devices may have a potential to realize a high-performance NRD based millimeter wave circuit. However, some improvements may possible to make optimization approaches more sophisticated for complex device functionalities. We are considering to design more interesting and efficient devices such as NRD-based isolators and circulators with attractive structural symmetrical conditions using developed evolutionary approaches.