A Multi-Constrained Optimization Method for THz Backward Wave Oscillators

: The current design period for various backward wave oscillators (BWOs) is still at least several months. How to ﬁnd the best structure parameters with an efﬁcient and stable optimization method is a problem facing researchers in both scientiﬁc research and engineering work. In this paper, a non-randomized iterative optimization method is proposed. It applies orthogonal design methods to ﬁnd local solutions that can provide optimal ‘gradient direction’ for several successive next iteration steps. An evaluation function is designed to distinguish the better ones from the local solutions in the multi-constrained optimization of such BWOs. Optimizations from different starting points are performed separately for a global optimal solution. Two BWOs at different frequency ranges are optimized using the proposed method. The validity and stability of the method are veriﬁed. It is believed that the method can provide the global optimum and shorten the design period of THz BWOs.


Introduction
Terahertz radiation (0.3-3 THz) has great potential in a variety of applications; however, there is still a lack of compact and powerful THz radiation sources, which has become a bottleneck for the development of terahertz science and technologies [1][2][3]. Among different kinds of THz radiation source, backward wave oscillators (BWOs) are competitive in the 0.3-1 THz range and have attracted much attention in recent years [1,4].
To fully explore the potential of a THz BWO, an optimization of the slow-wave structure (SWS) is essential [5]. Typically, a large amount of calculation is required during the optimization of the geometrical parameters of the SWS. The optimization design often takes several months if 3-D full-electromagnetic particle-in-cell simulations are required for obtaining accurate output performances. Very large amounts of possible combinations of the geometrical parameters make the optimization design of such devices very timeconsuming. Sometimes, the design cycle can be reduced to several weeks when analytical models are applied for calculation of the output results [6,7]. However, it is still hard to find the global optimum under specific constraints on the working parameters such as beam current, beam voltage and so on. Thus, an efficient and stable global optimization method is necessary for the optimization design of THz BWOs and other Cerenkov devices working in linear regions, no matter that the working performances are obtained by theoretical calculation or numerical simulation.
As for the global optimization method, there is still no algorithm that can be proved mathematically to converge to the global optimum. To find the global optimum as far as possible, there are generally two types of method. The first type is the evolutionary algorithm [8], such as genetic algorithm (GA), particle swarm optimization (PSO) and differential evolution (DE). All evolutionary algorithms share the same procedures, i.e., 2 of 15 reproduction, mutation, recombination and selection, which are inspired by biological evolution. They also have different properties: EA converges more slowly than the others, but they have the best stability without error caused by noise; PSO has the lowest route length or number of generations; however, DE has the lowest calculation time due to the smaller requirement of the population size [9,10]. The second type is the local search algorithm, together with the multi-start strategy [11]. Local search is a heuristic method that moves from solution to solution in the space of candidate solutions (the search space) by applying local changes until a local optimal is found. Local search can handle most problems with continuity and weak nonlinearity [11]. Even in the optimization of some complex systems, there is a trend toward the second type of method for better performances [12]. Which one is better for the optimization of THz BWOs?
In previous studies on the optimization design of vacuum electronic devices, randomized optimization methods (the first type) such as genetic algorithms (GAs) were used in the optimization design of relativistic backward wave oscillators (RBWOs) [13][14][15]. Unlike non-randomized iteration optimization methods, randomized optimization methods choose a next move not only from the local solutions but also from randomly generated solutions. These 'strange candidates' are the key factor for successful optimization of RBWOs because, in these high-power Cerenkov devices, the energy conversion takes place in non-linear regions and small changes of geometrical and working parameters can result in significantly different output performances. Non-randomized optimization methods have low efficiency in the optimization of such devices operating in nonlinear regimes.
However, for the optimization of THz BWOs that operate in linear regions, nonrandomized iteration optimization becomes suitable for at least three reasons [7]: (i) small changes in geometrical parameters will not result in sudden and big changes in the output performance [13]; (ii) unlike the random solutions, local solutions will satisfy the constraints on the working parameters; (iii) it is time-consuming to calculate the working parameters for a candidate solution, so randomized optimization can be inefficient in finding candidate solutions. Therefore, the non-randomized optimization method (the second type) is more suitable and is applied in this paper for the optimization design of THz BWOs.
To improve the efficiency of local search optimization, the orthogonal design method is included in the proposed optimization method. Previous studies have shown that iteration optimization methods such as local search algorithms (LSAs) can be improved using orthogonal design [7,16]. Orthogonal design is applied for generating local solutions. It can decrease the amount of calculation by several orders of magnitude. However, it still requires lots of calculations to produce local candidates at each iteration step using the gradient descent method [11]. To deal with this problem, a new greedy algorithm is proposed based on the investigation of the optimization problem that provides rapid iteration and can save lots of calculations. As for the multiple constrains, a specifically designed evaluation function is proposed for choosing the best candidate solution effectively. It will be proved that the proposed new optimization method using local search algorithms can converge to the same global optimum from different initial solutions and is more efficient than the previous method presented by us in Ref. [7].
The following sections of this paper are arranged as follows: the optimization problem is briefly described in Section 2; the optimization method is presented in Section 3; the validations of the method and the case studies are shown in Section 4, followed by the conclusion in Section 5.

Problem Description
A THz BWO mainly consists of a cathode, an anode, an SWS, an output structure and a collector, as shown in Figure 1 [17]. The electromagnetic wave on the SWS has a phase velocity that is slowed down to be approximately equal to the drifting velocity of the electrons so that the beam-wave interaction can occur under the Cerenkov radiation scheme. A corrugated rectangular waveguide SWS (CRWSWS) is under observation and shown in Figure 2. phase velocity that is slowed down to be approximately equal to the drifting velocity of the electrons so that the beam-wave interaction can occur under the Cerenkov radiation scheme. A corrugated rectangular waveguide SWS (CRWSWS) is under observation and shown in Figure 2.   [7]. L is the grating period; s is the gap width between teeth; h is the tooth height; w is the waveguide width; b is the waveguide height; N is the SWS period number. The light gray part shows the position and size of the electron beam that has a beam thickness in y direction and a beam width in x direction. The height of the beam tunnel is b-h. Figure 3 shows the ω-kz diagram of the CRWSWS [7], which describes the relationship between the frequency ω and the wave-number kz. According to Fluquet's theorem, an electromagnetic mode propagating on a CRWSWS consists of infinite space harmonics so that its dispersion curves (as shown in Figure 3) repeatedly expand along the wave-number axis. The straight lines in Figure 3 correspond to the electron beams with different beam voltages. The slopes of the lines are proportional to the beam velocities, which can be derived from the beam voltages. When the beam voltage is given, the operating frequency can be determined by the intersection of the corresponding beam line and the dispersion curve of the lowest electromagnetic mode. As a result, the SWS determines the relationship between the beam voltage and the operating frequency. Before optimizing the SWS of the device, the operating frequency range and expected beam voltage are commonly given. phase velocity that is slowed down to be approximately equal to the drifting velocity of the electrons so that the beam-wave interaction can occur under the Cerenkov radiation scheme. A corrugated rectangular waveguide SWS (CRWSWS) is under observation and shown in Figure 2.   [7]. L is the grating period; s is the gap width between teeth; h is the tooth height; w is the waveguide width; b is the waveguide height; N is the SWS period number. The light gray part shows the position and size of the electron beam that has a beam thickness in y direction and a beam width in x direction. The height of the beam tunnel is b-h. Figure 3 shows the ω-kz diagram of the CRWSWS [7], which describes the relationship between the frequency ω and the wave-number kz. According to Fluquet's theorem, an electromagnetic mode propagating on a CRWSWS consists of infinite space harmonics so that its dispersion curves (as shown in Figure 3) repeatedly expand along the wave-number axis. The straight lines in Figure 3 correspond to the electron beams with different beam voltages. The slopes of the lines are proportional to the beam velocities, which can be derived from the beam voltages. When the beam voltage is given, the operating frequency can be determined by the intersection of the corresponding beam line and the dispersion curve of the lowest electromagnetic mode. As a result, the SWS determines the relationship between the beam voltage and the operating frequency. Before optimizing the SWS of the device, the operating frequency range and expected beam voltage are commonly given.  [7]. L is the grating period; s is the gap width between teeth; h is the tooth height; w is the waveguide width; b is the waveguide height; N is the SWS period number. The light gray part shows the position and size of the electron beam that has a beam thickness in y direction and a beam width in x direction. The height of the beam tunnel is b-h. Figure 3 shows the ω-k z diagram of the CRWSWS [7], which describes the relationship between the frequency ω and the wave-number k z . According to Fluquet's theorem, an electromagnetic mode propagating on a CRWSWS consists of infinite space harmonics so that its dispersion curves (as shown in Figure 3) repeatedly expand along the wave-number axis. The straight lines in Figure 3 correspond to the electron beams with different beam voltages. The slopes of the lines are proportional to the beam velocities, which can be derived from the beam voltages. When the beam voltage is given, the operating frequency can be determined by the intersection of the corresponding beam line and the dispersion curve of the lowest electromagnetic mode. As a result, the SWS determines the relationship between the beam voltage and the operating frequency. Before optimizing the SWS of the device, the operating frequency range and expected beam voltage are commonly given.  When the position of the electron beam and the geometrical parameters of the SWS are given, the working parameters such as the beam current and the beam voltage of the device can be determined with theoretical models [7]. As shown in a previous study, the beam voltage and beam current are constrained due to mode competition [7]. When the When the position of the electron beam and the geometrical parameters of the SWS are given, the working parameters such as the beam current and the beam voltage of the device can be determined with theoretical models [7]. As shown in a previous study, the beam voltage and beam current are constrained due to mode competition [7]. When the device operates at the maximum beam voltage, the corresponding beam line should intersect only with the first backward space harmonic (n = −1) of the TE x10 dispersion curve to avoid the mode competition between electromagnetic modes. This constraint on the beam voltage can be quantitatively described as [7] where U 1 is the beam voltage corresponding to the beam line (Line 2 in Figure 3) that intersects with the low-frequency cut-off point of the −1 harmonic of the TEx10 dispersion curve and U is the beam voltage. The beam current should be smaller than the minimum starting current of the second lowest oscillation mode in the operating frequency range to avoid mode competition between oscillation modes [7]. In summary, the specific constraints considered in the optimization problem and the optimization objective (the fourth item) are listed below:

1.
The oscillation frequency is centered at a given GHz, ranging from f min to f max ; 2.
The maximum beam voltage (denoted as V m ) at f max is~16 kV and ∆U < 0 (the beam line has no intersection with the first backward harmonic of higher-order electromagnetic modes); 3.
The current density (denoted as J) is~500 A/cm 2 and the beam current should not exceed the minimum starting current of the second-lowest oscillation mode in the operating frequency range; 4.
The minimum output power (denoted as P) in the operating frequency range is maximized and the acceptable relative error for each constraint is 10%.

Optimization Method
The proposed constrained optimization method is based on the local search algorithm, which optimizes a solution to a constrained optimization problem by moving from solution to solution in the search space until a local optimal solution is found. Composed of three processes, the method will be introduced step by step. The first process is to choose an initial solution as the start point; the second process is to make small moves successively toward a local optimum (LO) solution; the third process is to search for the global optimum.

Start Point
It is helpful to introduce a general mathematical form of constrained optimization problems for further discussion. A constrained optimization problem (COP) has the following form [18]: Find which optimizes the objective function. R is the feasible region defined by a set of m additional constraints (m > 0): where t j and T j are the reasonable value and maximum permissible relative error for the jth constraint, respectively; g j ( → s ) is the jth working parameter being constrained. These denotations will be used in the following description.
As for the first step, the initial solution in the optimization can be obtained by scaling up a presented structure [17]. This BWO has the oscillation frequencies from 900 GHz to 1100 GHz and the required frequency ranges from f min to f max GHz. Hence, the scaling factor is the ratio between the original and expected center frequencies and the initial solution can be obtained by where L 0 = 50 µm; s 0 = 25 µm; b 0 = 110 µm; h 0 = 60 µm; w 0 = 240 µm; N 0 = 170. The frequency of 1000 GHz is the center frequency of the operating frequency range from 900 GHz to 1100 GHz in Ref. [17].

Successive Iteration for Local Optimum
This subsection can be divided into three parts. The first is to generate local solutions; the second is to evaluate them; the third is to move toward the local optimum.

Local Solutions
Orthogonal array or orthogonal design is applied. An index-one, strength-two orthogonal array which has v 2 rows and k columns (denoted as A v,k ) is introduced. Here, v is the number of levels; k is the number of parameters. Each row in an orthogonal array is a combination of levels alike the combination of the geometrical parameters of the SWS. The following is an example of A 3,4 : Define the orthogonal array as A v,k = [a i,j ], where the jth parameter in the ith row has a level of a i,j and a i,j ∈ {1,2, . . . ,v}, the sample taking from the local solutions can be described as where → s 0 = (p 1 , . . . , p k ) is a given solution; → s i is the ith local solution obtained by changing the given solution with the ith row of an orthogonal array; δ is a factor used for changing the parameters, for example δ = 1%, and v is an odd positive integer. Equation (5) can provide a sample of v 2 solutions in the local region near a given solution. However, the total number of possible combinations is k v . It can be observed from Figure 4 that the optional local solutions in this sample are uniformly distributed. This uniformity makes sure that better local solutions can always be found by the orthogonal design method.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 7 Equation (7) is given conservatively to guarantee that more than one local solu can pass the evaluation. It is reasonable to choose the local solution that has the la non-zero evaluated value without worrying about the satisfaction of the constraints.

Iterations toward Local Optimum by Greedy Algorithm
Being different from typical local search algorithms, the proposed optimiz method does not solve the best local solution for the next move at each iteration ste greedy algorithm is applied after finding the best local solution at the start. The be cal solution is not only considered as the next move, but also considered as a 'gra direction' along which successive moves are made without evaluating local solutio each step until the last move itself does not satisfy the constraints any more.
The successive moves, which reduce a lot of solving and evaluating calculat can be defined using a geometric series such as is the ith solution in the serie has been defined in Equation (4), α1, ..., αk are the levels in the orthogonal array c sponding to the best local solution, imax is the maximum number of iterations. The geometrical parameters of the successively variant solutions increase o crease exponentially with a slow rate. Figure 5 shows the corresponding working rameters as a function of the serial number that are calculated by Pierce's theory [7] monotonicity of each working parameter remains within several moves, but it cou broken as the ratio for the geometrical series may not influence the working param monotonically. Thus, a greedy strategy is let α1, ..., αk be constant during the succe iteration following Equation (7), which means that the local solutions will not be pled at each step for reducing the total calculation cost.
imax controls the number of iterations for each 'direction', i.e., combination of α αk. It is expected that the output performances and the working parameters can be proved monotonically during the imax times of iterations. Thus, this successive iter process should stop once the reached solution provides the output performances an working parameters that broke their previous monotonic variation. The last solu must be a bad choice for further iteration, and the second to last solution in the s

Evaluation of Local Solutions
As for the evaluation function, unit step functions are used to describe the constraints. The evaluation function is defined as and H is the unit step function. It can be observed from Equation (6) that a solution will be eliminated if the evaluation function value is equal to zero, which means that at least one of the H functions is equal to zero. Half of the optional solutions, whose working parameters are farther to the expectation values, are eliminated by each constraint, considering that the local samples are obtained by orthogonal design and distribute uniformly. The remaining solutions tend to satisfy the constraints since their working parameters are closer to the reasonable values. The reasonable value of v can be given as Equation (7) is given conservatively to guarantee that more than one local solution can pass the evaluation. It is reasonable to choose the local solution that has the largest non-zero evaluated value without worrying about the satisfaction of the constraints.

Iterations toward Local Optimum by Greedy Algorithm
Being different from typical local search algorithms, the proposed optimization method does not solve the best local solution for the next move at each iteration step. A greedy algorithm is applied after finding the best local solution at the start. The best local solution is not only considered as the next move, but also considered as a 'gradient direction' along which successive moves are made without evaluating local solutions at each step until the last move itself does not satisfy the constraints any more.
The successive moves, which reduce a lot of solving and evaluating calculations, can be defined using a geometric series such as where S ex is an example of series of solutions, S ex,i is the ith solution in the series, → s 1 has been defined in Equation (4), α 1 , . . . , α k are the levels in the orthogonal array corresponding to the best local solution, i max is the maximum number of iterations.
The geometrical parameters of the successively variant solutions increase or decrease exponentially with a slow rate. Figure 5 shows the corresponding working parameters as a function of the serial number that are calculated by Pierce's theory [7]. The monotonicity of each working parameter remains within several moves, but it could be broken as the ratio for the geometrical series may not influence the working parameters monotonically. Thus, a greedy strategy is let α 1 , . . . , α k be constant during the successive iteration following Equation (7), which means that the local solutions will not be sampled at each step for reducing the total calculation cost.
Appl. Sci. 2022, 12, x FOR PEER REVIEW will be suitable as a new initial solution for next successive moves. The search fo optimum will end if no better solution can pass the evaluation. Figure 6 shows the flow chart of the optimization for local optimum. A local will terminate if all the local solutions, having not been arrived at, have the same ated value of zero. The result should satisfy the constraints mentioned in Section 2  i max controls the number of iterations for each 'direction', i.e., combination of α 1 , . . . , α k . It is expected that the output performances and the working parameters can be improved monotonically during the i max times of iterations. Thus, this successive iteration process should stop once the reached solution provides the output performances and the working parameters that broke their previous monotonic variation. The last solution must be a bad choice for further iteration, and the second to last solution in the series will be suitable as a new initial solution for next successive moves. The search for local optimum will end if no better solution can pass the evaluation. Figure 6 shows the flow chart of the optimization for local optimum. A local search will terminate if all the local solutions, having not been arrived at, have the same evaluated value of zero. The result should satisfy the constraints mentioned in Section 2.
ated value of zero. The result should satisfy the constraints mentioned in Section

Optimization for Global Optimum
As for the global search, more than one local optimum could be found. The zation method for searching local optimum from given starting point has bee lished. If enough starting points can be provided, local optimums will be fou these starting points, respectively. Then, the global optimum is obtained from th optimums.
Apparently, to find few but enough starting points is important for searc global optimum effectively. It is always difficult to try all possible starting poi

Optimization for Global Optimum
As for the global search, more than one local optimum could be found. The optimization method for searching local optimum from given starting point has been established. If enough starting points can be provided, local optimums will be found from these starting points, respectively. Then, the global optimum is obtained from these local optimums.
Apparently, to find few but enough starting points is important for searching the global optimum effectively. It is always difficult to try all possible starting points. One reasonable method is to find new starting points based on the first-obtained local optimum (FOLO).
Let the FOLO be denoted as → s op = (q 1 , q 2 . . . , q k ). The local solutions near the FOLO can be given using Equation (5) as Successive moves from the FOLO can be performed for new starting points. From Equation (9), the successive moves, i.e., series of solutions, can be written as where → S i,j is the ith series; j max,i is defined as the last solution in the ith series, which breaks the specific constraints and can be solved during the iteration processes; → χ i is the 'moving direction' for the ith series; other variables have been defined; (q 1 , . . . , q k ) is → s op .
It can be known from Equation (10) that the successive moves have v 2 different 'directions' corresponding to the local solutions near the FOLO. The uniformity of these directions is guaranteed by the uniformity of these local solutions that are obtained by orthogonal design as described in Section 3.2.1. Therefore, at the ends of the successive moves from the FOLO along these 'directions', v 2 new starting points can be found which are fully-spreading and located near the boundary of the parameter space of feasible solutions. It should be mentioned that all these 'directions' are not the same as the 'gradient direction'. Therefore, the greedy procedure to find new starting points is different from the iteration procedure in a gradient descent algorithm. Figure 7 shows the flow chart of the global search algorithm. To avoid an endless loop, the maximum number of LOs is limited by k max = 1000. However, the search for the global optimum will commonly end due to the failure to find new LOs.
It should be mentioned that it is very time-consuming to search the global optimum from all v 2 new starting points. The search, in practice, typically needs only several new starting points because the iterations from different beginnings always end at an approximately same solution. An example will be described in detail in Section 4.2. Appl. Sci. 2022, 12, x FOR PEER REVIEW 10 of 16 It should be mentioned that it is very time-consuming to search the global optimum from all v 2 new starting points. The search, in practice, typically needs only several new starting points because the iterations from different beginnings always end at an approximately same solution. An example will be described in detail in Section 4.2.

Case Studies of THz BWOs
The proposed optimization method is implemented via MATLAB and applied for the optimization of an SWS shown in Figure 2 with constraints given in Section 2. Two cases are considered with different operation frequency ranges. However, other theories or mathematical tools applied are the same.
The analytical models presented by McVey [19] and Johnson [6] are used to calculate the working parameters of the device. A seven-level and six-factor orthogonal array is used to sample local solutions. The number of the constraints m is three and the number of levels of the orthogonal array is set to be seven, which is determined by Equation (7). The position of the electron beam is obtained using the empirical formulas suggested by P.C. Panda [20] and K.T. Nguyen [21]. The cathode is simplified by an emitting surface and the beam width and beam thickness will be one third of the waveguide width and two fifths of the height of the beam tunnel, respectively, considering the distribution of the electric field [17].

Case Studies of THz BWOs
The proposed optimization method is implemented via MATLAB and applied for the optimization of an SWS shown in Figure 2 with constraints given in Section 2. Two cases are considered with different operation frequency ranges. However, other theories or mathematical tools applied are the same.
The analytical models presented by McVey [19] and Johnson [6] are used to calculate the working parameters of the device. A seven-level and six-factor orthogonal array is used to sample local solutions. The number of the constraints m is three and the number of levels of the orthogonal array is set to be seven, which is determined by Equation (7). The position of the electron beam is obtained using the empirical formulas suggested by P.C. Panda [20] and K.T. Nguyen [21]. The cathode is simplified by an emitting surface and the beam width and beam thickness will be one third of the waveguide width and two fifths of the height of the beam tunnel, respectively, considering the distribution of the electric field [17].

First Case
For the first case, the operation frequency range is set to be 900-1100 GHz for comparison with the results in Ref. [22]. The conductivity of the SWS is 5.96 × 10 7 S/m, which is the same as that used in Ref. [22]. The optimized geometrical parameters are L = 52 µm, s = 22 µm, h = 61 µm, w = 215 µm, b = 124 µm, and N = 217. The current density is 520 A/cm 2 , the beam voltage ranges from 8 to 16 kV, and the guiding magnetic field is 1 T. The dispersion curves of the optimized SWS are shown in Figure 8. The 16-kV beam voltage line has no intersection point with the dispersion curves of higher-order modes as shown in Figure 8. The optimized geometrical parameters are the same as the optimized parameters in a previous paper [7] because (i) two methods converge to the same solution and (ii) the parameters are rounded to be integer numbers in the units of µm.
For the first case, the operation frequency range is set to be 900-1100 GHz for com parison with the results in Ref. [22]. The conductivity of the SWS is 5.96 × 10 7 S/m, whic is the same as that used in Ref. [22]. The optimized geometrical parameters are L = 5 μm, s = 22 μm, h = 61 μm, w = 215 μm, b = 124 μm, and N = 217. The current density i 520 A/cm 2 , the beam voltage ranges from 8 to 16 kV, and the guiding magnetic field is T. The dispersion curves of the optimized SWS are shown in Figure 8. The 16-kV beam voltage line has no intersection point with the dispersion curves of higher-order mode as shown in Figure 8. The optimized geometrical parameters are the same as the opti mized parameters in a previous paper [7] because (i) two methods converge to the sam solution and (ii) the parameters are rounded to be integer numbers in the units of μm. The starting currents of the lowest and second lowest oscillation modes are show in Figure 9a. The beam current can be calculated as 0.4(b-h)×0.33w×520 A/cm 2 = 9 mA. I can be observed that the starting currents of the second lowest mode are always large than 9 mA in the operation frequency range. The output powers of the optimized BWO are compared with the results in Ref. [22] in Figure 9b. The output powers are simulate using CST Particle Studio, and the time signal and its spectrum are shown in Figure 10 The optimized BWO has acceptable output powers that are close to those published i Ref. [22]. Under the same beam voltage range and beam current, the optimized BWO satisfies all constraints, delivers a stable output signal at correct frequencies and pro vides comparable output powers with the BWO in [22]. Thus, the case under discussio demonstrates the validity of the proposed method.
It should be mentioned that the optimization results in Figures 8-10 are close bu different from those in Ref. [22]. Note that the second constraint requires the beam lin not to intersect with the dispersion curve of TEx20 mode in the −1-harmonic region however, in Ref. [22], the beam line is permitted to do this. The starting currents of the lowest and second lowest oscillation modes are shown in Figure 9a. The beam current can be calculated as 0.4(b − h) × 0.33w × 520 A/cm 2 = 9 mA. It can be observed that the starting currents of the second lowest mode are always larger than 9 mA in the operation frequency range. The output powers of the optimized BWO are compared with the results in Ref. [22] in Figure 9b. The output powers are simulated using CST Particle Studio, and the time signal and its spectrum are shown in Figure 10. The optimized BWO has acceptable output powers that are close to those published in Ref. [22]. Under the same beam voltage range and beam current, the optimized BWO satisfies all constraints, delivers a stable output signal at correct frequencies and provides comparable output powers with the BWO in [22]. Thus, the case under discussion demonstrates the validity of the proposed method.

Second Case
For the second case, the same structure of SWS and same constraints listed in Sec. II have been applied again. The operation frequency range changes to 800-900 GHz, corresponding to a THz window in air. Following Equation (4) The conductivity of the SWS is also 5.96 × 10 7 S/m. The optimized geometrical parameters are L = 61 μm, s = 33 μm, h = 79 μm, w = 250 μm, b = 130 μm and N = 147. The current density is 543 A/cm 2 , the beam voltage ranges from 8 to 15.8 kV, and the guiding magnetic field is 1 T. The dispersion curves of the optimized SWS are shown in Figure  11a. The 16 kV beam voltage line has no intersection point with the dispersion curves of higher-order modes as shown in Figure 11a. The output powers of the optimized BWO are compared with the theoretical results as shown in Figure 11b. The theoretical output powers can estimate the simulated output powers with acceptable accuracy. The simulated time signal at the beam voltage of 12 kV and its spectrum are shown in Figure 12. The optimized BWO delivers stable output signal and the maximum efficiency is 0.52 mW/9.2 mA/15.8 kV = 0.35%, which is close to the maximum efficient of 0.34% in Ref. [22]. Considering that the beam size in cross-plane, and the current density and the beam voltage are the same, it is reasonable to obtain a similar total efficiency for the optimized BWO in 800-900 GHz. It should be mentioned that the optimization results in Figures 8-10 are close but different from those in Ref. [22]. Note that the second constraint requires the beam line not to intersect with the dispersion curve of TE x20 mode in the −1-harmonic region; however, in Ref. [22], the beam line is permitted to do this.

Second Case
For the second case, the same structure of SWS and same constraints listed in Sec. II have been applied again. The operation frequency range changes to 800-900 GHz, corresponding to a THz window in air. Following Equation (4), the beginning can be obtained → s case2,1 = (L 0 , s 0 , b 0 , h 0 , w 0 , N 0 ) 1000 850 (11) The conductivity of the SWS is also 5.96 × 10 7 S/m. The optimized geometrical parameters are L = 61 µm, s = 33 µm, h = 79 µm, w = 250 µm, b = 130 µm and N = 147. The current density is 543 A/cm 2 , the beam voltage ranges from 8 to 15.8 kV, and the guiding magnetic field is 1 T. The dispersion curves of the optimized SWS are shown in Figure 11a. The 16 kV beam voltage line has no intersection point with the dispersion curves of higher-order modes as shown in Figure 11a. The output powers of the optimized BWO are compared with the theoretical results as shown in Figure 11b. The theoretical output powers can estimate the simulated output powers with acceptable accuracy. The simulated time signal at the beam voltage of 12 kV and its spectrum are shown in Figure 12. The optimized BWO delivers stable output signal and the maximum efficiency is 0.52 mW/9.2 mA/15.8 kV = 0.35%, which is close to the maximum efficient of 0.34% in Ref. [22]. Considering that the beam size in cross-plane, and the current density and the beam voltage are the same, it is reasonable to obtain a similar total efficiency for the optimized BWO in 800-900 GHz.
To check whether the optimized structure is the global optimum, two more optimizations are performed. The beginning geometrical parameters are The optimization results from three different beginnings are listed in Table 1. It can be observed that the three results are close to each other although they correspond to quite different beginnings. This demonstrates the validity of the proposed global optimization method. mization method.
It takes only 16 steps for the searching algorithm to reach a LO as shown in Figure  13. Figure 13 shows that the optimization process has three durations (≤11, [12][13][14][15][16] corresponding to three different gradient directions. The initial solution has out-ranged working parameters and the optimization pulls the working parameters back within the constraints successfully. As relative errors are allowed, the optimized working parameters may exceed the expectation values.   quite different beginnings. This demonstrates the validity of the proposed global optimization method. It takes only 16 steps for the searching algorithm to reach a LO as shown in Figure  13. Figure 13 shows that the optimization process has three durations (≤11, [12][13][14][15][16] corresponding to three different gradient directions. The initial solution has out-ranged working parameters and the optimization pulls the working parameters back within the constraints successfully. As relative errors are allowed, the optimized working parameters may exceed the expectation values.    It takes only 16 steps for the searching algorithm to reach a LO as shown in Figure 13. Figure 13 shows that the optimization process has three durations (≤11, [12][13][14][15][16] corresponding to three different gradient directions. The initial solution has out-ranged working parameters and the optimization pulls the working parameters back within the constraints successfully. As relative errors are allowed, the optimized working parameters may exceed the expectation values.  A comparison is made between the previous optimization method in Ref. [7] and the proposed method in this paper. In Ref. [7], the evaluation function takes a common form where Ul is the upper limit of the beam voltage and Jl is the upper limit of the beam current. The optimization in Ref. [7] takes 160 iteration steps to stop. However, it takes only 16 steps for the proposed optimization method to stop. From Equation (13), it can be found that the evaluation function in Ref. [7] is too strict so that the working parameters tend to converge exactly at the upper limits. However, the best candidate solution may not be chosen due to the deviations of the working parameters, which results in the low efficiency of the method in Ref. [7]. The minimum output power curve in Figure 13 is only auxiliary, which shows the lowest output power from all local solutions for each iteration step. Note that these calculations for local solutions are not necessary for the optimization. The curve is drawn here to show the tendency of the variation of the output power. The curve goes down because the beam current and the beam voltage have to be smaller to satisfy the constraints.

Conclusions
Mode competition restricts the beam voltages and beam currents of THz BWOs. This makes the optimization of this kind of device a multi-constrained problem. At the same time, it requires no random samples during the optimization process.
A new local search algorithm combined with greedy algorithm is proposed for the optimization of THz BWOs. Iterations are made in a successive way, holding a fixed 'direction' to save calculation consumption. Local solutions are generated uniformly by orthogonal design and evaluated by a newly designed formula that is suitable for the multi-constrained optimization problem. A comparison is made between the previous optimization method in Ref. [7] and the proposed method in this paper. In Ref. [7], the evaluation function takes a common form F(U, J, P) = Pe −|U−U l |/U l −|J−J l |/J l (13) where U l is the upper limit of the beam voltage and J l is the upper limit of the beam current. The optimization in Ref. [7] takes 160 iteration steps to stop. However, it takes only 16 steps for the proposed optimization method to stop. From Equation (13), it can be found that the evaluation function in Ref. [7] is too strict so that the working parameters tend to converge exactly at the upper limits. However, the best candidate solution may not be chosen due to the deviations of the working parameters, which results in the low efficiency of the method in Ref. [7]. The minimum output power curve in Figure 13 is only auxiliary, which shows the lowest output power from all local solutions for each iteration step. Note that these calculations for local solutions are not necessary for the optimization. The curve is drawn here to show the tendency of the variation of the output power. The curve goes down because the beam current and the beam voltage have to be smaller to satisfy the constraints.

Conclusions
Mode competition restricts the beam voltages and beam currents of THz BWOs. This makes the optimization of this kind of device a multi-constrained problem. At the same time, it requires no random samples during the optimization process.
A new local search algorithm combined with greedy algorithm is proposed for the optimization of THz BWOs. Iterations are made in a successive way, holding a fixed 'direction' to save calculation consumption. Local solutions are generated uniformly by orthogonal design and evaluated by a newly designed formula that is suitable for the multi-constrained optimization problem.
Two cases are investigated to verify the proposed method. The optimized BWOs have reasonable output performances and satisfy the constraint well. In the second case, it takes only 16 moves to obtain a local optimum. Three searches from different beginnings coincidentally find three local optimal solutions that are close to each other. The results demonstrate that the proposed optimization method can find the global optimum effectively, efficiently and stably.