Optimal Design of Aperture Illuminations for Microwave Power Transmission with Annular Collection Areas

Xun Li , Baoyan Duan , Yiqun Zhang , Yongxin Guo

Engineering ›› 2023, Vol. 30 ›› Issue (11) : 63 -74.

PDF (2836KB)
Engineering ›› 2023, Vol. 30 ›› Issue (11) :63 -74. DOI: 10.1016/j.eng.2023.07.016
Research
Article
Optimal Design of Aperture Illuminations for Microwave Power Transmission with Annular Collection Areas
Author information +
History +
PDF (2836KB)

Abstract

This work presents an optimal design method of antenna aperture illumination for microwave power transmission with an annular collection area. The objective is to maximize the ratio of the power radiated on the annular collection area to the total transmitted power. By formulating the aperture amplitude distribution through a summation of a special set of series, the optimal design problem can be reduced to finding the maximum ratio of two real quadratic forms. Based on the theory of matrices, the solution to the formulated optimization problem is to determine the largest characteristic value and its associated characteristic vector. To meet security requirements, the peak radiation levels outside the receiving area are considered to be extra constraints. A hybrid grey wolf optimizer and Nelder-Mead simplex method is developed to deal with this constrained optimization problem. In order to demonstrate the effectiveness of the proposed method, numerical experiments on continuous apertures are conducted; then, discrete arrays of isotropic elements are employed to validate the correctness of the optimized results. Finally, patch arrays are adopted to further verify the validity of the proposed method.

Graphical abstract

Keywords

Microwave power transmission / Beam collection efficiency / Ring-shaped beam / Annular collection area / Grey wolf optimizer / Nelder-Mead simplex method

Cite this article

Download citation ▾
Xun Li, Baoyan Duan, Yiqun Zhang, Yongxin Guo. Optimal Design of Aperture Illuminations for Microwave Power Transmission with Annular Collection Areas. Engineering, 2023, 30 (11) : 63-74 DOI:10.1016/j.eng.2023.07.016

登录浏览全文

4963

注册一个新账户 忘记密码

1. Introduction

Microwave power transmission (MPT) is a method to wirelessly deliver energy at microwave frequencies from a generation point to at least one receiver point. MPT has gained widespread attention since it was first proposed, because it eliminates the infrastructure needed to transmit power. Potential applications of MPT include powering unmanned aerial vehicles [1], charging electric vehicles [2], providing energy from one satellite to another [3], supplying energy to Internet of Things devices [4], and delivering power to forward operating bases [5]. In fact, most MPT developments have been driven by the progress of space solar power satellites (SSPSs). An SSPS [6] is a huge MPT system that collects solar power and converts it to direct current (DC) power in outer space, then transmits the DC power to earth via MPT technology. An MPT system mainly consists of a transmitting antenna and a rectenna. The transmitting antenna functions as a convex lens to focus the microwave beam on the rectenna, while the rectenna intercepts the incident microwave power and converts it back to DC power.

In an MPT system, the efficient transmission of microwave power to the target receiver is key. To evaluate this performance, beam collection efficiency (BCE), which is defined as the ratio of the power that impinges on the rectenna aperture to the total transmitted power, is usually adopted [7]. A great deal of work has been done to improve the BCE in both array antennas and continuous apertures [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19]. To be specific, the optimal array synthesis problem for maximizing the BCE of linear or planar arrays can be addressed by exploiting discrete prolate spheroidal sequences [8] or by solving a generalized eigenvalue problem [9], [10]. Numerical optimization methods such as linear programming [11], convex programming and compressive sensing [12], k-means clustering [13], contiguous partitioning [14], genetic algorithms [15], and particle swarm optimization (PSO) algorithms [16] have also been adopted to design arrays with high BCEs. For continuous apertures, stepped amplitude distribution [17] and isosceles trapezoidal distribution [18] have been proposed to reduce the transmitting antenna complexity while simultaneously ensuring a high BCE. Apart from BCE, other performance indexes such as the aperture power coefficient of the transmitting antenna [19] and the radiated power density variation on the rectenna [20] can be well addressed with multi-objective optimization techniques.

It is notable that the available antenna design approaches for MPT generally focus on the optimization design of spot beams with either circular or square shapes. However, special scenarios exist for which ring-shaped beams are required. For example, in some applications, it is necessary to direct the energy beam to the perimeter of an area rather than to its center. As a typical example, a radio telescope may be constructed in a lunar crater, with the power for the monitoring and control system around its rim being supplied by an MPT system [21]. As another example, it may be necessary to supply energy to a settlement surrounding a mountain. In addition, the use of an MPT system to charge wireless sensor networks around an active volcano holds great value for monitoring volcanic eruptions [22], [23], since a real-time monitoring and prediction system for detecting volcanoes is essential to save lives. Apart from MPT, ring-shaped beam antennas are considered to be good candidates for satellite communications and wireless local-area networks (WLANs) [24], [25]. Moreover, laser beams with a ring-shaped intensity distribution have many important applications in life science and technology [26], [27], [28], [29].

In this article, an optimal design method for antenna aperture illumination that generates a ring-shaped beam for MPT is proposed. The design goal is to achieve the maximum BCE. To address security concerns, constraints on peak radiation levels (PRLs) outside the annular collection area are considered. This paper is organized as follows. In Section 2, optimization designs of antenna illuminations for MPT with annular collection areas are formulated, first without (Section 2.1) and then with security concerns (Section 2.2). Next, a hybrid grey wolf optimizer (GWO) and Nelder-Mead (NM) simplex method is proposed in Section 2.3 to deal with the constrained optimization problem formulated in Section 2.2. To demonstrate the effectiveness of the proposed method, representative numerical experiments on continuous apertures are conducted in Section 3.1. Subsequently, to further confirm the validity of the optimized results, array antennas of isotropic elements and patch elements are analyzed in 3.2 Array antennas of isotropic elements for MPT with annular collection areas, 3.3 Array antennas of patch elements for MPT with annular collection regions. Finally, Section 4 concludes the paper.

2. Problem formulation and solution methods

2.1. Optimal aperture illumination design for MPT with an annular collection area

Fig. 1 shows an MPT system in which the transmitting antenna has a circular aperture with a radius of R t. The rectenna has an annular shape with inner and outer radii of R r 1 and R r 2, respectively. The transmitting antenna and rectenna are assumed to be aligned and separated by a distance of L in the Fresnel region. For simplicity, a circularly symmetric aperture distribution is considered for the transmitting antenna. Let E t ρ , ψ be the normalized aperture distribution, which can be written as follows:

E t ρ , ψ = g ρ e x p j ψ ρ

where g ρ and ψ ρ denote the aperture amplitude and phase distribution, respectively. The imaginary unit j is - 1 ρ = r / R t indicates the normalized radial distance, and r is the distance from the transmitting antenna center to another point on the transmitting aperture. To focus the transmitted beam at a distance of L in the Fresnel region (radiative near field), the transmitting antenna should be equipped with a phase distribution in the following form[30]:

ψ ρ = β ρ 2 R t 2 2 L

where β = 2 π / λ is the wavenumber and λ is the wavelength. This phase taper compensates for the phase difference due to the difference in the distances between each source point on the aperture and the focal point. Then, field contributions are added in phase at the focal point. In fact, it has been theoretically proven that, in the focal plane (i.e., the rectenna plane) near the axis, the electric field of the transmitting antenna has all the properties of the far field [30]. This conclusion is verified in Refs. [10], [31], [32]. Based on this fact, the radiation pattern of the circular transmitting antenna (E) is given by the following:

E ϑ = j β R t 2 e - j β L L F ϑ

where

F ϑ = 0 1 g ρ J 0 ϑ ρ ρ d ρ

and J 0 is the 0th-order Bessel function of the first kind, ϑ = β R t s i n θ , s i n θ = r ' / L , θ is the elevation angle, and r ' denotes the radial distance from the center of the beam at the rectenna site. Suppose the transmitting antenna radiates a peak microwave beam intensity of I t o ; then, the radiated beam intensity(I)on the rectenna plane can be calculated as follows:

I ϑ = I r ' = I  to  β 2 R t 4 L 2 F 2 ϑ

Referring to Fig. 1, the power confined in the annular collection area (Pr) is calculated by

P r = 0 2 π R r 1 R r 2 I r ' r ' d r ' d φ = 2 π I t o R t 2 ϑ 1 ϑ 2 F 2 ϑ ϑ d ϑ

where

ϑ i = β R t s i n θ i , i = 1 , 2

and θ 1 and θ 2 are the two angles of the annular area shown in Fig. 1. The total transmitted power P t can be calculated as follows:

P t = 2 π I t o R t 2 0 1 g 2 ρ ρ d ρ

Based on Eqs. (6), (8), the BCE of an MPT system with an annular collection area can be represented by

B C E = P r P t = ϑ 1 ϑ 2 F 2 ϑ ϑ d ϑ 0 1 g 2 ρ ρ d ρ

From Eqs. (4) and (9), it can be seen that, for a focused aperture, the antenna aperture amplitude g ρ plays a key role in determining the achievable BCE. To find the optimal g ρ that focuses the microwave beam on the annular collection area, it is helpful to express g ρ as a series [20]; that is,

g ρ = n = 1 N x n 1 - ρ 2 n - 1

where the n th basis function has the form 1 - ρ 2 n - 1, which is generally used to approximate aperture amplitudes that taper toward the edges of the apertures [33], and x n is the associated weight factor. It can be seen from Eq. (10) that various g ρ s can be formed by choosing a different truncation parameter N and weight factor x n. When N = 1, Eq. (10) can be reduced to a uniform aperture amplitude. If we let x = x 1 , , x N T and A = 1 , , 1 - ρ 2 N - 1 T, then Eq. (10) can be rewritten as follows:

g ρ = x T A

Accordingly, the denominator of Eq. (9) can be rewritten by

0 1 g 2 ρ ρ d ρ = x T B x

where

B = 0 1 A A T ρ d ρ

is an N × N matrix, and the(m, n)th element of B has a closed-form solution; that is,

B m n = 1 2 m + n - 1

A detailed derivation of Eq. (14) is given in Appendix A. In addition, by substituting Eq. (10) into Eq. (4), F ϑ can be rewritten as follows:

F ϑ = n = 1 N x n 0 1 1 - ρ 2 n - 1 J 0 ϑ ρ ρ d ρ

It is worth noting that the n th component of the integral in Eq. (15) also has a closed-form solution [33]; that is,

0 1 1 - ρ 2 n - 1 J 0 ϑ ρ ρ d ρ = 2 n - 1 n - 1 ! J n ϑ ϑ n (16) where J n denotes the n th-order Bessel function of the first kind, and the symbol "!" represents the factorial function. If we let C = J 1 ϑ ϑ , , 2 N - 1 N - 1 ! J N ϑ ϑ N T, then Eq. (15) can be rewritten as follows:

F ϑ = x T C

Accordingly, the numerator of Eq. (9) becomes

ϑ 1 ϑ 2 F 2 ϑ ϑ d ϑ = x T D x

where

D = ϑ 1 ϑ 2 C C T ϑ d ϑ

is an N × N matrix and the(m, n)th element of D is calculated by

D m n = ϑ 1 ϑ 2 2 m - 1 m - 1 ! J m ϑ ϑ m × 2 n - 1 n - 1 ! J n ϑ ϑ n ϑ d ϑ

Based on Eqs. (12), (18), the BCE formulated in Eq. (9) can be reduced to the ratio of two real quadratic forms:

B C E = x T D x x T B x

To maximize the achievable BCE, the optimal design variable vector x  opt , is given by

x  opt  = a r g m a x x x T D x x T B x

From Eqs. (14) and (20), it is clear that, for every m and n, B m n = B n m and D m n = D n m. Hence, B and D are two symmetric matrices. Also, for an arbitrary real vector x 0, we have x T D x > 0 and x T B x > 0, since they represent, respectively, the power radiated on the annular collection area and the total transmitted power. Thus, x T D x and x T B x are two positive definite quadratic forms. Based on the theory of matrices [34], the solution to Eq. (22) is to determine the largest characteristic value ω m a x and its associated characteristic vector; that is,

D x  opt  = ω m a x B x  opt 

The maximum BCE is equal to ω m a x, with the associated aperture amplitude distribution being g ρ = x  opt  T A. In this work, Eq. (23) is solved using the "eig" function in Matlab (MATLAB, USA), employing the Lapack package [35].

2.2. Optimal aperture illumination design for MPT with an annular collection area and with security constraints

The optimal design of aperture distribution for MPT in terms of maximizing the BCE is given in Section 2.1. It should be noted, however, that this design fails to deal with the PRLs outside the receiving area. In practice, in many MPT applications, a large amount of power is transferred. Thus, microwave radiation safety must be taken into account. To address this issue, extra constraints on the PRLs outside the receiving area must be considered. In this case, the optimization design of the aperture distribution can be transformed into a constrained optimization problem:

m i n x f x = - B C E x
s.t.   P R L 1 = 20 l o g 10 m a x ϑ ϑ 1 F ϑ m a x ϑ F ϑ C 1
P R L 2 = 20 l o g 10 m a x ϑ ϑ 2 + Δ ϑ F ϑ m a x ϑ F ϑ C 2
- 1 x n 1 n = 1 , , N

where f is the objective function, x = x 1 , , x N T is the design variable vector used to define the shape of g ρ. The objective is to maximize the BCE. To change the maximization problem to a standard minimization problem, the objective function is multiplied by -1, as shown in Eq. (24). Eq. (25) is a constraint used to ensure that the PRL at the edge of or inside region 1 (shown in Fig. 1), denoted by P R L 1, is below C 1 dB. Similarly, Eq. (26) is a constraint used to ensure that the PRL at the edge of or outside the exclusion zone [36], denoted by P R L 2, is below C 2 d B.

2.3. Solution strategy to the optimization problem formulated in Section 2.2

It is clear that Eqs. (24), (25), (26) represent a constrained nonlinear optimization problem. To handle this problem, we first convert the constrained optimization problem to an unconstrained one using a penalty method:

f x = - B C E x + K × m a x 0 , P R L 1 + C 1 , m a x 0 , P R L 2 + C 2

where K is a penalty factor and is set to 10 6. It can be seen that, when P R L 1 or P R L 2 is larger than the desired value C 1 or C 2, Eq. (28) yields a very large value. When the constraints in Eqs. (25) and (26) are met, the optimization will go on to search for the maximum BCE. In order to solve Eq. (28), a hybrid GWO and NM optimization method (GWO-NM) is proposed that combines the advantages of GWO [37] and NM [38]. Details of the GWO-NM algorithm are provided below.

2.3.1. Grey wolf optimizer

GWO [37] is a swarm intelligence optimization technique that imitates the social hierarchy and group hunting behavior of grey wolves in nature. In GWO, all "wolves" (i.e., solutions) are divided into four kinds based on their fitness, simulating the social structure of wild wolves. The best solution is defined as the alpha; the second and third best solutions are named the beta and delta, respectively; and the remaining solutions are called the omega. In a hunt, three group hunting strategies are employed-namely, searching for prey, encircling prey, and attacking prey. To simulate the encircling behavior of grey wolves, the following two equations are used:

X t + 1 = X P t - E G
G = F X P t - X t

where X t + 1 and X t are vectors of dimension N 1 that are used to indicate the positions of a wolf (candidate solutions) at the t + 1 th and t th iterations, respectively, and N 1 is the number of design variables. X P t is an N 1 -dimensional vector denoting the position of the prey (potential optimal solution). E and F are two N 1 -dimensional vectors, defined as follows:

E = 2 r 1 a - a
F = 2 r 2

where r 1 and r 2 are random vectors of dimension N 1 whose components are within the interval 0 , 1 . a is an N 1 -dimensional vector whose elements decrease linearly from 2 to 0 throughout the iteration:

a = 2 1 - t T

where t and T are the current iteration and the maximum number of iterations, respectively. By adjusting E and F, a solution, X t, can adjust its position with respect to X P t in an N 1 -dimensional search space to mimic the encircling behavior of grey wolves.

In a wolf pack, the hunt is led by the pack leaders, who are assumed to have better knowledge of the position of the prey. Other wolves then follow the pack leaders to approach the prey. This group hunting mechanism can be represented by the following:

G α = F 1 X α - X , G β = F 2 X β - X , G δ = F 3 X δ - X  
X 1 = X α - E 1 G α , X 2 = X β - E 2 G β , X 3 = X δ - E 3 G δ
X t + 1 = X 1 + X 2 + X 3 / 3

where X α , X β , X δ, and X are the positions of alpha, beta, delta, and omega in the t th iteration, respectively. X t + 1 is the updated position of a wolf in the t + 1 th iteration. For a better understanding of GWO, interested readers are directed to Ref. [37].

GWO has been shown to have good global search ability [37] but poor local search capability. In addition, similar to other population-based optimization techniques, GWO has a low convergence speed compared with gradient-based optimization methods. To address this issue, the NM simplex algorithm is incorporated into GWO.

2.3.2. NM simplex algorithm

The NM simplex algorithm [38] is a derivative-free optimization method with a powerful local search ability. It has a very fast convergence speed and has been widely used in unconstrained optimization problems. For an optimization problem with N 1 design variables, the NM simplex algorithm starts by forming a simplex Δ with N 1 + 1 initial vertices, x l l = 1 , , N 1 + 1, each of which represents a candidate solution. A simplex is a geometrical objective generated by N 1 + 1 points in an N 1 -dimensional space. For example, a triangle is a simplex in two-dimensional (2D) space, and a tetrahedron is a simplex in three-dimensional (3D) space.

In the NM algorithm, all the N 1 + 1 vertices x l l = 1 , , N 1 + 1 are sorted in ascending order based on their objective function values; that is,

f x 1 f x 2 f x N 1 + 1

The vertex yielding the minimum objective function value, x 1, is named the best vertex. Similarly, x N 1 + 1 is referred to as the worst vertex. At each iteration, a new simplex is formed by replacing the worst vertex with a newly generated vertex or by shrinking the old simplex while keeping its best vertex unchanged. Four operations (reflection, expansion, contraction, and shrinking) are employed to shape the simplex, each of which is associated with a parameter: γ (reflection), δ (expansion), ε (contraction), and σ (shrinking). These parameters are selected to satisfy γ > 0 , δ > 1, δ > γ , 0 < ε < 1, and 0 < σ < 1 [38]. The k th iteration of the NM algorithm is given below [38].

(1) Evaluate the objective function. Evaluate the objective function f at the N 1 + 1 vertices x l k l = 1 , , N 1 + 1 of the simplex Δ k, and sort them so that Eq. (37) holds.

(2) Perform the reflection operation. Obtain the reflection point x r k :

x r k = x k ¯ + γ x k ¯ - x N + 1 k

where x k ¯ is the centroid of the first N 1 best vertices

x k ¯ = l = 1 N 1 x l k / N 1

Evaluate f x r k. If f x r k < f x 1 k, go to step 3; if f x 1 k f x r k < f x N 1 k, replace x N 1 + 1 k with x r k and go to step 7 ; if f x N 1 k f x r k < f x N 1 + 1 k, go to step 4 ; otherwise, go to step 5.

(3) Perform an extension operation. Obtain the extension point x e k :

x e k = x k ¯ + δ x r k - x k ¯ (40) Evaluate f x e k. If f x e k < f x r k, replace x N 1 + 1 k with x e k and go to step 7; otherwise, replace x N 1 + 1 k with x r k and go to step 7.

(4) Perform an outside contraction. Obtain the outside contraction point x o c k :

x o c k = x k ¯ + ε x r k - x k ¯

Evaluate f x o c k. If f x o c k < f x r k, replace x N 1 + 1 k with x o c k and go to step 7; otherwise, go to step 6.

(5) Perform an inside contraction. Obtain the inside contraction point x i c k :

x i c k = x k ¯ - ε x k ¯ - x N 1 + 1 k

Evaluate f x i c k. If f x i c k < f x N 1 + 1 k, replace x N 1 + 1 k with x i c k and go to step 7; otherwise, go to step 6.

(6) Perform a shrinking operation. All the vertices of the simplex except the best one, that is, x l k l = 2 , , N 1 + 1, are replaced by new vertices:

x l k = x 1 k + σ x l k - x 1 k

(7) Determine whether the stopping condition has been met. Sort the vertices of the new simplex, x l k l = 1 , , N 1 + 1, such that Eq. (37) holds. If the maximum number of function evaluations (MNFE) is reached, or if f x 1 k - f x N 1 + 1 k ϵ, where ϵ is a user-defined small predetermined tolerance, then stop the algorithm; otherwise, k = k + 1, in which case, go to step 2.

For clarity, Fig. 2 shows the effects of reflection, expansion, contraction, and shrinking for a simplex in 2D space, with the coefficients γ = 1.0 , δ = 2.0 , ε = 0.5, and σ = 0.5. With successive iterations, the simplex gradually converges to the optimal point. As suggested in Ref. [38], the coefficients of reflection, expansion, contraction, and shrinking are set to the following:

γ = 1 , δ = 1 + 2 N 1 , ε = 0.75 - 2 2 N 1 , σ = 1 - 1 N 1

NM is a very efficient local search method. However, the obtained result is extremely sensitive to the initial points. Thus, the initial points should be carefully selected.

2.3.3. The GWO-NM algorithm

As described above, GWO has good exploration ability but possesses the drawbacks of poor local search ability and slow convergence speed. While the NM simplex method has a good exploitation capability and a fast convergence speed, the result obtained is highly dependent on the initial solutions. To make full use of the advantages of these two algorithms, a hybrid GWO-NM algorithm is proposed. There are two stages in our GWO-NM: a coarse global search stage and an intensive local search stage. In the first stage (the coarse global search stage), GWO is adopted to explore the search space globally and quickly find a promising search space. In the second stage (the intensive local search stage), the NM simplex method is utilized to find a high-quality solution by performing an intensive local search based on the promising solutions found by GWO.

In the GWO-NM, the initial vertices x l l = 1 , , N 1 + 1 of the simplex are generated in two ways. First, we execute GWO N 1 times, and the best solutions are stored as the N 1 initial vertices x l l = 1 , , N 1. In order to speed up the convergence speed of the GWO-NM, the N 1 + 1 th initial vertex of the simplex is obtained by solving Eq. (23)-that is, x N 1 + 1 = x  opt . With these N 1 + 1 initial solutions, an N 1 -dimensional simplex is formed, and the process of NM is triggered. A diagram of the GWO-NM algorithm is provided in Fig. 3.

3. Numerical analysis and discussion

In this section, the optimal design of antenna aperture distributions for MPT with annular collection areas with or without security constraints in the Fresnel region will be conducted and discussed. First, optimization designs of continuous apertures with different receiving areas will be presented. With a quadratic phase taper, the radiation pattern in the Fresnel region can be approximated to the far-field pattern [10], [30]. Based on this fact, for simplicity, array antennas of isotropic elements working in the far-field region are employed to verify the validity of the optimized results. The excitation coefficients are obtained by sampling the optimized continuous distributions. Finally, the proposed method is further validated with patch arrays working in the far-field region, while considering the mutual coupling effect. For all the numerical experiments, the working frequency is set to 5.8 GHz, which lies within the atmospheric window.

3.1. Continuous aperture illumination designs for MPT with annular collection areas

Recall the antenna amplitude distribution $g(ρ)$ in Eq. (10), which is formed by a summation of a series of the form 1 - ρ 2 n - 1 and order of N. To find the optimal g ρ that maximizes the achievable BCE, the truncation parameter N should be determined. For this purpose, a design of g ρ for maximizing the BCE without constraint is considered. Here, ϑ 1 and ϑ 2 in Eq. (9) are respectively set to 3 and 9, and N is first set to 4. By solving Eq. (23), the maximum BCE is found to be 96.047%. The optimal design variable vector, x  opt , is given in Table 1. Then, we increase N from 4 to { 5 , 6 , 7 , 8 , 9 , 10 }, respectively. For each N, Eq. (23) is solved; the optimized BCEs and the associated optimal design variable vectors are summarized in Table 1. It is found that, when N is larger than 7, the optimized BCE tends to be stable. Therefore, in this article, N is chosen as 8.

Fig.4(a) plots the optimized $g(ρ)$ (solid blue curve), which has apeak value in the antenna center, then drops to below zero, andfinally increases to above zero. It should be noted that there is anarea in which $g(ρ) < 0 $, which means that the antenna should befed 180° out of phase. With this unusual amplitude distribution,the power radiated by the transmitting antenna can be mainlyfocused on the annular collection area,since a very high BCE ofabout 97.59% is obtained. Fig. 4(b) shows the associated normal-ized radiation pattern (solid blue curve). Although a very highBCE is achieved,the PRL in region 1 $\(\vartheta \le \vartheta_9 \)$ is rather high(PRL1 = — 6.44 dB), which appears at the edge of region 1. Such ahigh PRL may cause interference to the electronic equipment orpose a hazard to people nearby.

For security concerns, constraints on the PRL outside the annular receiving area must be considered. In the first set of design cases, suppose that the PRL outside the exclusion region is required to be below -20 dB; that is, C 2 in Eq. (26) is set to - 20 d B, and the exclusion zone Δ ϑ is set to 1. Different PRLs in region 1 will be considered, with the goal of maximizing the BCE. To deal with this constrained optimization problem, the proposed GWO-NM is adopted. For the GWO algorithm in the first stage of the GWO-NM, the population size(nPop)is set to n P o p = 20 and the maximum number of iterations is set to T = 200, which are both very small for a population-based optimization algorithm. For the NM simplex algorithm in the second stage of the GWO-NM, MNFE is set to MNFE = n P o p × T, and ϵ is set to 1 × 10 - 6.

As the first design case, the PRL in region 1 is required to be below - 18 d B C 1 = - 18 d B. With the GWO-NM, the optimized BCE is 93.09 %, and the optimal design variable vector, x  opt , is given in Table 2. We then gradually decrease C 1 and solve the constrained optimization problem in sequence. This process is stopped when no feasible solution can be found. The optimized BCEs and the corresponding design variable vectors are presented in Table 2. For clarity, the relationship between C 1 and the associated optimized BCE is plotted in Fig. 5. Clearly, there is a linear relationship between the achievable BCE and C 1. The suppression limit of C 1 is about - 29 d B, and the associated BCE is reduced to 89.25 %. For the purposes of comparison and clarity, only the optimized g ρ s obtained for C 1 = - 20 , - 25, and -29dB are plotted in Fig. 4(a). The associated radiation patterns are depicted in Fig. 4(b), where the left vertical black dotted line denotes the inner edge of the receiving area ϑ = ϑ 1, while the right vertical black dotted line indicates the outer edge of the exclusion area ϑ = ϑ 2 + Δ ϑ. From Fig. 4(b), it can be seen that the PRL in region 1 is suppressed at the expense of raising other sidelobes outside the exclusion zone. It should be noted that the area outside the exclusion region is much larger than that of region 1. Thus, as C 1 is suppressed, much more power will be distributed outside the exclusion region. Hence, BCE decreases as C 1 is reduced. In addition, it is found that the directions of the main beams are shifted. However, the radiated power can be well confined in the annular collection area, which is confirmed by the high BCEs obtained, as shown in Table 2.

In the second set of design cases, ϑ 1 and ϑ 2 in Eq. (9) are set to 4 and 10, respectively. Solving Eq. (23) shows that the optimized BCE is as high as 97.27 %, and the optimal design variable vector is x  opt  = [ 0.0051 , 0.0135 , - 0.0363 , 0.0057 , - 0.3694 , 0.6381 , - 0.5634 , 0.3707 ] T. Fig. 6(a) plots the optimized g ρ. Similar to the first design case, the optimized g ρ has the maximum value in the antenna center, then decreases to below zero, and finally increases to above zero. Fig. 6(b) plots the associated radiation pattern, and a relatively high PRL P R L 1 = - 10.67 d B is observed at the edge of region 1.

To reduce the PRL while simultaneously achieving a high BCE, the GWO-NM algorithm is employed to handle the constrained optimization problem formulated in Section 2.2. Here, C 2 in Eq. (26) is fixed at - 20 d B, and Δ ϑ is set to 1 ; in Eq. (25), different C 1 values are considered. The parameters of the GWO-NM algorithm are set to be the same as in the first set of design cases. Table 3 summarizes the optimized BCEs and the associated optimal design variable vectors concerning different C 1 values. To be specific, when C 1 is set to - 18 d B, the BCE drop can be neglected 0.42 % compared with the unconstrained optimal one (96.85% vs 97.27%). The suppression limit of C 1 is about - 22 d B. In this case, the optimized BCE is still high enough (BCE = 95.28 % ). For comparison, three constrained optimized g ρ s with C 1 = - 18 , - 20, and - 22 d B are plotted in Fig. 6(a). It can be seen that all the optimized g ρ s have similar shapes. Fig. 6(b) displays the associated radiation patterns, clearly showing that the radiation patterns strictly meet the design constraints.

To investigate the performance of the GWO-NM, a comparison study was carried out on the GWO-NM versus GWO and PSO in dealing with the constrained optimization problem. Due to space limitations, only one design case is considered. For this problem, ϑ 1 and ϑ 2 in Eq. (9) are respectively set to 3 and 9, and C 1 and C 2 are respectively set to -18 and -20. In the first stage of the GWO-NM algorithm, the population size of GWO is set to n P o p = 20, and the maximum number of iterations is set to T = 200. For the NM simplex algorithm in the second stage of the GWO-NM, the MNFE is set to MNFE = n P o p × T = 4000. For GWO and PSO, the population size is set to n P o p = 100, and the maximum number of iterations is set to T = 1000. Other parameters for the two algorithms are set the same as those shown in Table 1 in Ref. [39]. To obtain statistical results, all three algorithms are independently run five times. The optimized results and associated computation time are summarized in Table 4. The central processing unit (CPU) adopted for the numerical simulations was an Intel ® Xeon ® E-2224G at 3.5 G H z with 32 G B random access memory (RAM). The numerical analysis software was Matlab R2018a [35].

It is clear that both the GWO and PSO consume much more time than the GWO-NM. Moreover, they have lower success rates, at 40% and 20%, respectively. Here, the success rate refers to the ratio of the number of an algorithm to successfully find a feasible solution (PRL1 < −18 dB and PRL2 < −20 dB) to the total number of trials. The proposed GWO-NM has the capacity to find a stable optimal solution in each independent run. From this analysis, it can be concluded that the GWO-NM has both a good exploitation capability and a fast convergence speed compared with GWO and PSO.

3.2. Array antennas of isotropic elements for MPT with annular collection areas

Section 3.1 presented the optimization design of continuous aperture distributions for MPT with annular collection areas with and without constraints. From a practical standpoint, it is difficult (or sometimes impossible) to design a continuous aperture antenna with optimized aperture distribution. However, optimized continuous distributions can serve as references when designing array antennas of arbitrary sizes, since the array excitation coefficients can be easily determined by sampling the continuous distributions. In order to illustrate this point and to demonstrate the validity of the optimized results, numerical experiments on planar arrays with different aperture sizes and receiving regions in the far field were conducted and are described in this section. The antenna elements are considered to be isotropic sources.

In the first set of numerical experiments, the optimized g ρ s obtained for ϑ 1 = 3 and ϑ 2 = 9 without constraint and with the constraints C 1 = - 20 d B and C 2 = - 20 d B are used. The diameter of the circular transmitting array, D t, is assumed to vary from 5 λ to { 10 λ , 15 λ , 20 λ , 25 λ , 30 λ }. The receiving angles θ i i = 1 , 2 shown in Fig. 1 can be calculated as follows:

θ i = s i n - 1 2 ϑ i β D t , i = 1 , 2

The circular transmitting array is positioned to be positioned along a rectangular grid in the xoy plane, and the inter-element spacing is d x = d y = 0.5 λ in both directions. The method of forming the circular array is as follows. First, construct a square array with side length D t. For this square array, there are P P = D t / d x rows and Q Q = D t / d y columns of radiating elements. Then, calculate the distance r p q from the array center to the element located at the p th p 1 , P row and q th q 1 , Q column; that is,

r p q = p - P + 1 2 d x 2 + q - Q + 1 2 d y 2

Finally, remove the elements whose r p q > D t / 2 from the square array, forming a circular array with a diameter of D t. For clarity, the upper right quadrant of a circular array of D t = 10 λ is plotted in Fig. 7.

Once the circular array is available, the excitation coefficient of the(p, q)th element can be obtained using I p q = g ρ p q, where ρ p q = 2 r p q / D t denotes the normalized radial distance. After obtaining the excitation coefficient of each radiating element, the associated BCE can be calculated as follows [9]:

B C E = P r P t = 0 2 π θ 1 θ 2 A F u , v 2 s i n θ d θ d φ 0 2 π 0 π A F u , v 2 s i n θ d θ d φ

where A F u , v is the array factor; u = s i n θ c o s φ , v = s i n θ s i n φ are direction cosines; and φ is the azimuth angle.

Similar to the continuous apertures, for array antennas, we are concerned with the achievable BCE and the associated PRL in region 1 P R L 1 and outside the exclusion region P R L 2. The detailed performance indexes for a transmitting array with different diameters are calculated and presented in Table 5. Here, The subscripts "U" and "C" indicate the results associated with the optimized g ρ without constraint and with the constraints C 1 = - 20 d B and C 2 = - 20 d B. As expected, the achievable BCE increases as the array diameter increases for both the unconstrained and constrained cases. The optimized BCEs for the continuous apertures are the theoretical limits of discrete arrays with the same ϑ 1 and ϑ 2. In addition, it is observed that, for the unconstrained case, when the diameter of the array is larger than 10 λ, the resulting BCE becomes very close to that of the continuous aperture. Although a very high BCE can be obtained, the resulting P R L 1 is relatively high (about - 7 d B ), which is not desirable. For the constrained cases, the associated P R L 1 values are all reduced to below - 20 d B. These low P R L 1 values are achieved at the expense of small BCE drops, which is consistent with the results obtained for the continuous case. For the P R L 2 value, as the array size increases, it becomes close to - 20 d B. For clarity and due to space limitations, only 3 D radiation patterns for the array with D t = 10 λ are shown in Fig. 8. The two circles marked by dash-dotted lines refer to the inner and outer boundaries of the annular collection area. It can be clearly seen that, in both cases, the radiated power can be mostly focused on the receiving region, which is numerically confirmed by the high BCEs obtained B C E U = 97.574 % and B C E C = 90.206 % ). In addition, comparing Fig. 8(b) with Fig. 8(a) shows that the P R L 1 is greatly reduced for the array with excitation coefficients obtained by sampling the optimized g ρ with the constraints C 1 = - 20 d B and C 2 = - 20 d B.

In the second set of numerical experiments, the optimized g ρ s obtained for ϑ 1 = 4 and ϑ 2 = 10 without constraint and with the constraints C 1 = - 20 d B and C 2 = - 20 d B are used to obtain the array excitation coefficients. The diameter of the transmitting array is also varied from 5 λ to 30 λ. The performance indexes for the B C E , P R L 1, and P R L 2 are calculated and presented in Table 6. It is notable that, in this set of test cases, as the array size increases, the obtained BCE first increases and then fluctuates to approach that of the continuous case. This phenomenon may be caused by the discretization error of g ρ. Furthermore, it is found that there is little difference in the BCE for the same array with and without the PRL constraints, which is consistent with the continuous case. In addition, for those arrays D t 10 λ whose excitation coefficients are obtained by sampling the optimized g ρ with the constraints C 1 = - 20 d B and C 2 = - 20 d B, the associated P R L 1 s and P R L 2 s are very close to - 20 d B. Due to space limitations, only 3D power patterns for the array with D t = 10 λ are plotted and shown in Fig. 9. It can be seen that most of the radiated power is confined in the annular collection area for both cases. With the optimized constrained g ρ, the P R L 1 is reduced from -11.35 to - 20.03 d B, which is accompanied by a small reduced B C E B C E C = 96.117 % vs B C E U = 96.889 % ).

3.3. Array antennas of patch elements for MPT with annular collection regions

In this subsection, the validity of the proposed method is verified by using real antenna elements. For simplicity, a patch antenna (shown in Fig. 10(a)) is used. The ground plane size of the patch element is 25.86mm×25.86mm, and a 1.00 mm thick Rogers 5880 substrate is used. The square patch has a side length of 16.25 mm and is centered at the middle of the ground plane. The feed position is shown in Fig. 10(a). This patch element is designed to resonate at the center frequency of 5.8 GHz. It should be noted that the radiation pattern of an isolated element is different from the one embedded in an array, due to the mutual coupling effect. Indeed, the element patterns are all different in the array; thus, it is often difficult to model the exact mutual coupling effect among elements.

One possible way to model the mutual coupling effect is to use an embedded element pattern (EEP) [40]. The EEP refers to the pattern of a single antenna element embedded in a finite array. To obtain the EEP of the patch element, a 5×5 patch array (shown in Fig. 10(b)) is adopted; only the center element is excited, with all other elements being terminated in 50 Ω. A commercial full-wave simulator, Ansys HFSS (ANSYS, USA), is adopted to extract the radiation pattern of the embedded element. Fig. 10(c) plots the 3D gain pattern of the embedded element, which incorporates the effect of coupling with the neighboring elements. Under the hypothesis that most element patterns are the same as the EEP (ignoring the edge effects), the radiation pattern of the patch array can be approximated as follows:

F u , v = A F u , v × E E P u , v

Accordingly, the BCE for the patch arrays can be calculated by replacing A F u , v in Eq. (47) by F u , v in Eq. (48). In this set of test cases, the diameter of the circular patch array is fixed at 10 λ, and the elements are half-wavelength spaced in both the x and y directions. The array consists of 316 patch elements. Two sets of excitation coefficients are used, which are obtained by sampling the optimized g ρ s ϑ 1 = 4 and ϑ 2 = 10 without constraint and with the constraints C 1 = - 20 d B and C 2 = - 20 d B. The associated performance parameters for the B C E , P R L 1, and P R L 2 are calculated and listed in Table 7. For the unconstrained optimized solution, the PRL outside the exclusion region is about - 22 d B ; however, the PRL in region 1 is relatively high(-12.61dB). For the constrained optimized solution, the PRLs outside the exclusion region and in region 1 are both below - 20 d B. The corresponding normalized power patterns are shown in Fig. 11.

In order to further study the validity of the proposed method, the patch array was simulated by means of Ansys HFSS. Fig. 12 plots the normalized radiation patterns, with the associated B C E , P R L 1, and P R L 2 summarized in Table 7. A comparison of Fig. 12(a) with Fig. 11(a) and Fig. 9(a), and of Fig. 12(b) with Fig. 11(b) and Fig. 9(b), reveals that only small deviations between them are noticeable. This is confirmed by the BCE and PRL values obtained. In fact, the BCE and PRL values of the patch arrays are slightly improved compared with the ideal arrays of isotropic elements. Through this study, the validity of the proposed method is further confirmed.

Although the effectiveness of the proposed method is verified by arrays working in the far field, the obtained results can also be applied to the radiative near field when a quadratic phase taper is employed. Finally, it must be stressed that the proposed method also applies to MPT applications in which circular collection areas are of interest. This is because a circular collection area is only a special case of an annular collection area where the inner radius of the annular collection area is zero. Therefore, there is no doubt that the proposed method will find various applications in practical engineering.

4. Conclusions

In this paper, an optimal design method of antenna aperture illumination used for MPT with an annular collection area with and without security concerns is proposed. The design goal is to achieve the maximum BCE. After formulating the aperture illumination by means of a summation of a special set of series, the unconstrained optimal design problem is revealed to be finding the maximum ratio of two real quadratic forms. The problem can then be solved mathematically. To meet security requirements, constraints on the PRLs outside the annular collection area are considered. A hybrid GWO-NM method is proposed to deal with the constrained optimization problem and is demonstrated to quickly find the optimal solutions. With the proposed method, continuous aperture distributions yielding the maximum BCE with/without extra constraints can be achieved. Then, array antennas with arbitrary sizes can be easily designed. Notably, the proposed method is also applicable to MPT applications in which circular collection areas are of interest, since a circular collection area is only a special case of an annular collection area.

Acknowledgments

This work was supported in part by the National Key Research and Development Program of China (2021YFB3900300), in part by the National Natural Science Foundation of China (62201416), in part by the Fundamental Research Funds for the Central Universities (QTZX23070), in part by the Qin Chuang Yuan High-Level Innovative and Entrepreneurial Talents Project (QCYRCXM-2022-314), and in part by Singapore Ministry of Education Academic Research Fund Tier 1.

References

[1]

Satoru S, Nguyen DH, Nishioka Y, Shimamura K, Mori K, Yokota S.The logistics system by rotary wing unmanned aerial vehicle with 28 GHz microwave power transmission. In:Proceedings of IEEE Wireless Power Transfer Conference (WPTC); 2019 Jun 18-21; London, UK; 2019.

[2]

Shinohara N. Wireless power transmission progress for electric vehicle in Japan. In: Proceedings of 2013 IEEE Radio and Wireless Symposium; 2013 Jan 20-23; Austin, TX, USA; 2013.

[3]

C. Bergsrud, J. Straub. A space-to-space microwave wireless power transmission experiential mission using small satellites. Acta Astronaut, 103 ( 2013), pp. 193-203

[4]

L. Sun, L. Wan, K. Liu, X. Wang. Cooperative-evolution-based WPT resource allocation for large-scale cognitive industrial IoT. IEEE Trans Industr Inform, 16 (8) ( 2020), pp. 5401-5411 DOI: 10.1109/tii.2019.2961659

[5]

C.T. Rodenbeck, P.I. Jaffe, B.H. Strassner II, P.E. Hausgen, J.O. McSpadden, H. Kazemi, et al.. Microwave and millimeter wave power beaming. IEEE J Microw, 1 (1) ( 2021), pp. 229-259 DOI: 10.1109/jmw.2020.3033992

[6]

X. Li, B. Duan, L. Song, Y. Yang, Y. Zhang, D. Wang. A new concept of space solar power satellite. Acta Astronaut, 136 ( 2017), pp. 182-189

[7]

X. Li, K.M. Luk, B. Duan. Aperture illumination designs for microwave wireless power transmission with constraints on edge tapers using bezier curves. IEEE Trans Antennas Propag, 67 (2) ( 2019), pp. 1380-1385 DOI: 10.1109/tap.2018.2884850

[8]

S. Prasad. On an index for array optimization and the discrete prolate spheroidal functions. IEEE Trans Antennas Propag, AP-30 (5) ( 1982), pp. 1021-1023

[9]

G. Oliveri, L. Poli, A. Massa. Maximum efficiency beam synthesis of radiating planar arrays for wireless power transmission. IEEE Trans Antennas Propag, 61 (5) ( 2013), pp. 2490-2499

[10]

S. Kojima, T. Mitani, N. Shinohara. Array optimization for maximum beam collection efficiency to an arbitrary receiving plane in the near field. IEEE Open J Antennas Propag, 2 ( 2021), pp. 95-103 DOI: 10.1109/ojap.2020.3044443

[11]

A.F. Morabito, A.R. Laganà, T. Isernia. Optimizing power transmission in given target areas in the presence of protection requirements. IEEE Antennas Wirel Propag Lett, 14 ( 2015), pp. 44-47

[12]

A.F. Morabito. Synthesis of maximum-efficiency beam arrays via convex programming and compressive sensing. IEEE Antennas Wirel Propag Lett, 16 ( 2017), pp. 2404-2407

[13]

X. Li, B. Duan, L. Song. Design of clustered planar arrays for microwave wireless power transmission. IEEE Trans Antennas Propag, 67 (1) ( 2019), pp. 606-611

[14]

Rocca P, Oliveri G, Massa A. Innovative array designs for wireless power transmission. In: Proceedings of IEEE MTT-S International Microwave Workshop Series on Innovative Wireless Power Transmission: Technologies, Systems, and Applications; 2011 May 12-13; Kyoto, Japan; 2011.

[15]

N. Anselmi, A. Polo, M.A. Hannan, M. Salucci, P. Rocca. Maximum BCE synthesis of domino-tiled planar arrays for far-field wireless power transmission. J Electromagn Wave, 34 (17) ( 2020), pp. 2349-2370

[16]

X. Li, B. Duan, J. Zhou, L. Song, Y. Zhang. Planar array synthesis for optimal microwave power transmission with multiple constraints. IEEE Antennas Wirel Propag Lett, 16 ( 2017), pp. 70-73

[17]

X. Li, B. Duan, L. Song, Y. Zhang, W. Xu. Study of stepped amplitude distribution taper for microwave power transmission for SSPS. IEEE Trans Antennas Propag, 65 (10) ( 2017), pp. 5396-5405

[18]

A.K.M. Baki, N. Shinohara, H. Matsumoto, K. Hashimoto, T. Mitani. Study of isosceles trapezoidal edge tapered phased array antenna for solar power station/satellite. IEICE Trans Commun, E90-B (4) ( 2007), pp. 968-977 DOI: 10.1093/ietcom/e90-b.4.968

[19]

X. Li, Y. Guo. Multiobjective optimization design of aperture illuminations for microwave power transmission via multiobjective grey wolf optimizer. IEEE Trans Antennas Propag, 68 (8) ( 2020), pp. 6265-6276 DOI: 10.1109/tap.2020.2981736

[20]

X. Li, K. Luk, B. Duan. Multiobjective optimal antenna synthesis for microwave wireless power transmission. IEEE Trans Antennas Propag, 67 (4) ( 2019), pp. 2739-2744 DOI: 10.1109/tap.2019.2893312

[21]

Potter SD. Specialized phased-array antenna patterns for wireless power and information transmission. In: Proceedings of Space Manufacturing 10 Pathways to the High Frontier; 1995 May 4-7; Princeton, NJ, USA; 1995.

[22]

N. Takabayashi, N. Shinohara, T. Mitani, M. Furukawa, T. Fujiwara. Rectification improvement with flat-topped beams on 2.45-GHz rectenna arrays. IEEE Trans Microw Theory Tech, 68 (3) ( 2020), pp. 1151-1163 DOI: 10.1109/tmtt.2019.2951098

[23]

Prasad D, Hassan A, Verma DK, Sarangi P, Singh S. Disaster management system using wireless sensor network: a review. In:Proceedings of 2021 International Conference on Computational Intelligence and Computing Applications (ICCICA); 2021 Nov 26-27; Nagpur, India; 2021.

[24]

S. Son, S. Jeon, C. Kim, W. Hwang. GA-based design of multi-ring arrays with omnidirectional conical beam pattern. IEEE Trans Antennas Propag, 58 (5) ( 2010), pp. 1527-1535

[25]

D. Hua, S. Qi, W. Wu, D. Fang. Synthesis of conical beam array antenna with concentric loop configuration using element-level pattern diversity technique. IEEE Trans Antennas Propag, 66 (11) ( 2018), pp. 6397-6402 DOI: 10.1109/tap.2018.2867041

[26]

I. Manek, Y.B. Ovchinnicov, R. Grimm. Generation of a hollow laser beam for atom trapping using an axicon. Opt Commun, 147 (1) ( 1998), pp. 67-70

[27]

G. Roosen, C. Imbert. The TEM*01 mode laser beam—a powerful tool for optical levitation of various types of spheres. Opt Commun, 26 (3) ( 1978), pp. 432-436

[28]

B. Shao, S.C. Esener, J.M. Nascimento, E.L. Botvinick, M.W. Berns. Dynamically adjustable annular laser trapping based on axicons. Appl Opt, 45 (25) ( 2006), pp. 6421-6428

[29]

J.F. Guan, Z. Shen, X. Ni, J. Lu, J. Wang, B. Xu. Numerical simulation of the ultrasonic waves generated by ring-shaped laser illumination patterns. Opt Laser Technol, 39 (6) ( 2007), pp. 1281-1287

[30]

J.W. Sherman. Properties of focused apertures in the Fresnel region. IEEE Trans Antennas Propag, 10 (4) ( 1962), pp. 399-408

[31]

S. Karimkashi, A.A. Kishk. Focused microstrip array antenna using a Dolph-Chebyshev near-field design. IEEE Trans Antennas Propag, 57 (12) ( 2009), pp. 3813-3820

[32]

A. Buffi, A.A. Serra, P. Nepa, H.T. Chou, G. Manara. A focused planar microstrip array for 2.4 GHz RFID readers. IEEE Trans Antennas Propag, 58 (5) ( 2010), pp. 1536-1544

[33]

C.A. Balanis. Antenna theory: analysis and design. ( 3rd ed.), Wiley, New York City ( 2005)

[34]

F.R. Gantamacher. The theory of matrices. Chelsea, New York City ( 1959)

[35]

Version 9.4 ( R2018a), MathWorks. Natick: MATLAB. 2018.

[36]

S.D. Potter. Optimization of microwave power transmission from solar power satellites [dissertation]. New York University, New York City ( 1993)

[37]

S. Mirjalili, S.M. Mirjalili, A. Lewis. Grey wolf optimizer. Adv Eng Softw, 69 ( 2014), pp. 46-61

[38]

F. Gao, L. Han. Implementing the Nelder-Mead simplex algorithm with adaptive parameters. Comput Optim Appl, 51 (1) ( 2012), pp. 259-277 DOI: 10.1007/s10589-010-9329-3

[39]

X. Li, Y.X. Guo. Grey wolf optimizer for antenna optimization designs: continuous, binary, single-objective, and multiobjective implementations. IEEE Antennas Propag Mag, 64 (6) ( 2022), pp. 29-40

[40]

R.J. Mailloux.Phased array antenna handbook. (2nd ed.), Artech House, Norwood ( 2005)

Funding

the National Key Research and Development Program of China(2021YFB3900300)

the National Natural Science Foundation of China(62201416)

the Fundamental Research Funds for the Central Universities(QTZX23070)

the Qin Chuang Yuan High-Level Innovative and Entrepreneurial Talents Project(QCYRCXM-2022-314)

Singapore Ministry of Education Academic Research Fund Tier 1.

PDF (2836KB)

Supplementary files

Supplementary Material

4624

Accesses

0

Citation

Detail

Sections
Recommended

/