The Optimal Antenna Layout for Maximum Ergodic Capacity of MISO Beamforming System

Abstract

The performance of the multiantenna systems can be deeply affected by the antenna layout, that is, the placement of the antenna elements in a given area. This paper focuses on the optimal antenna layout which can maximize the ergodic capacity of a multi-input single-output (MISO) beamforming system. Since the optimization problem cannot be solved analytically, we use simulated annealing (SA) algorithm to search for the optimal antenna layout. The results show that the optimal antenna layout depends highly on the power azimuth spectrum (PAS) distribution and power elevation spectrum (PES) distribution. In the three-dimensional (3D) space uniformly scattering environment, the optimal antenna layout will locate the antenna elements on the surface of a sphere. When the scattering is concentrated in a two-dimensional (2D) plane, the circular layout is the optimal and finally, when there is no scattering (single plane wave), the optimal layout is the nonuniform linear array.

1. Introduction

multiple-input multiple-output (MIMO) systems have received significant attention for their capability of increasing channel capacity and improving the transmission reliability [1, 2]. Beamforming (BF) is one particular way of utilizing multiple antenna elements [3, 4]. In a beamforming antenna array, the radiated electromagnetic field can be patterned to have high antenna gain in desired direction by controlling the phase and amplitude of the signal field at each antenna element. Antenna beamforming has found a diverse range of applications due to its simplicity and effectiveness [5–8].

For a practical system, it is typical that multiple antenna elements have to be installed in a limited space area. The constrained space area will inevitably introduce the correlation between antenna elements and such correlation will affect the performance of multiantenna system. There are already many published works discussing the relation between the antenna correlation and the performance of multiple antenna systems in terms of ergodic capacity, outage capacity, bit error rate, and so forth [9–14]. However, to our knowledge, there are relatively few reports [15–22] on the relation between the multiple antenna system performance and the antenna layout, that is, the placement of the antenna elements. It is well known that, given the power azimuth spectrum (PAS) distribution and power elevation spectrum (PES) distribution [23], the performance of the multiantenna system is completely determined by the antenna spatial correlation matrix with its entries determined by the relative distance of the antenna elements. Hence, it is clear that the performance of the multiple antenna system is a function of antenna layout. In case of 2-dimensional (2D) deployment of multiple antenna elements and uniform PAS, it is reported [19, 20] that the optimal antenna layout which can maximize MIMO capacity is the near circular array, but both of them considered a propagation model including PAS only. For short-range communications, [21] has investigated the optimal antenna placement of 2×2 MIMO taken into consideration the joint effects of path loss and phase difference. In [22], empirical results for the effects of antenna placement on vehicle-to-vehicle link performance in vehicular ad hoc networks have been reported.

This paper focuses on the 3-dimensional (3D) deployments of multiple transmit antenna elements. In particular, we are interested in the optimal antenna elements placement which can maximize the ergodic capacity of a multiple-input single-output (MISO) BF system. Since it is difficult (if not impossible) to find the optimal solution analytically, simulated annealing (SA) algorithm [24–26] is used to search for the optimal antenna layouts for different PAS and PES. The results show that the optimal antenna layout is highly dependent on the PAS and PES distribution. In the 3D space uniformly scattering environment, the layout with antenna elements located on a surface of sphere will obtain the best ergodic capacity. When the spherical power spectrum (SPS) [23] degrades to 2D distribution (with only PAS or PES), the circular layout is the optimal. While there is no angle spreading, the nonuniform linear array is the best.

The rest of the paper is organized as follows. The next section introduces the system model including the 3D spatial correlation model and the ergodic capacity expressed as the function of eigenvalues of antenna correlation matrix. In Section 3, we use the simulated annealing to search for the optimal antenna layouts for some representative PAS and PES. Section 4 concludes the paper.

Throughout the paper, the following notations are adopted. All the matrices and the vectors are denoted by boldfaced capital letters and boldfaced small letters, respectively. (⋅)^𝐻 and (⋅)^𝑇 denote the complex conjugate transpose and transpose operator, respectively. |⋅| is the absolute value and ‖⋅‖ is the two-norm operator. 𝐸{⋅} represents the expectation. CN(0,𝐀) denotes a complex circular-symmetric Gaussian random vector with zero mean and variance matrix 𝐀.

2. System Model

2.1. MISO Beamforming System

Consider an MISO system shown in Figure 1(a) with 𝐿 transmit antennas and a single receive antenna, where 𝑑_max is the maximum allowable separation distance between any two antenna elements and 𝑑 is the distance between the transmitter and receiver. For most applications, 𝑑 is much larger than 𝑑_max. Suppose that all the antennas are mounted with the same orientation in 3D space, as shown in Figure 1(a); we can reasonably assume that the path loss, shadowing, and the antenna gain (due to the antenna pattern) are the same for all antenna elements.

(a) MISO system

(b) Electromagnetic waves (rays) departing from different directions. For each ray, 𝛼 denotes the elevation and 𝛽 denotes the azimuth

Figure 1 

System model.

Based on the principle of reciprocity, the characteristics of an antenna array will keep the same whatever it is used for transmitting or receiving. Consider the transmit antenna array shown in Figure 1(b); the electromagnetic waves (rays) will depart from the array in different directions, due to the scattering in the environment. For each ray, we use 𝛼 to denote the elevation and 𝛽 the azimuth. The total energy is distributed along these directions and the distribution is denoted, respectively, by 𝑃_𝑒(𝛼) and 𝑃_𝑎(𝛽), namely, the power azimuth spectrum (PAS) and the power elevation spectrum (PES) [14, 27].

The MISO channel is modeled as a 𝐿×1 vector 𝐡=(ℎ₁,ℎ₂,…,ℎ_𝐿)^𝑇, where ℎ_𝑙(𝑙=1,2,…,𝐿) is the channel coefficient from the 𝑙th transmit antenna to the receive antenna. Based on this assumption, we can model 𝐡 as a zero mean complex circular-symmetric Gaussian vector CN(0,𝐑), where 𝐑 is the transmit antenna correlation matrix defined as 𝐑=𝐸𝐡𝐡^𝐻. (1)

The received signal at the receiver is 𝑟=𝐡^𝑇𝐰𝑠+𝑛, (2) where 𝑛 is the zero mean complex additive Gaussian white noise (AWGN) with variance 𝑁₀ and 𝑠 is the scalar data symbol with 𝐸[|𝑠|²]=𝐸_𝑠. Before transmission on antenna 𝑙, the symbol is weighted by a complex number 𝑤_𝑙. The weights for all transmit antennas can be collected into a 𝐿×1 BF vector 𝐰=(𝑤₁,𝑤₂,…,𝑤_𝐿)^𝑇.

Assuming perfect channel state information (CSI) at the transmitter, then the optimal BF weighting vector is given by 𝐡𝐰=^𝐻‖𝐡‖, (3) and the instantaneous signal to noise ratio (SNR) at the output of receiver is given by 𝐸𝛾=^𝑠𝑁₀‖𝐡‖²=𝛾₀‖𝐡‖², (4) where 𝛾₀=𝐸_𝑠/𝑁₀.

2.2. Correlation between Antenna Elements

The elements of 𝐑 defined in (1) are the correlation coefficient of a pair of transmit antenna elements. With the classical multipath model [28–31], for any two transmit antennas located at 𝐩=(𝑥,𝑦,𝑧) and 𝐩^′=(𝑥^′,𝑦^′,𝑧^′), the 3D spatial correlation is given by 𝑅𝐩,𝐩^′=^𝜋/2_−𝜋/2^𝜋_−𝜋exp𝑗2𝜋𝜆𝑥−𝑥^′+cos𝛼cos𝛽𝑦−𝑦^′+cos𝛼sin𝛽𝑧−𝑧^′sin𝛼×𝑃_𝑒(𝛼)𝑃_𝑎𝛽d𝛼d𝛽, (5) where 𝜆 is the wavelength and 𝑃_𝑒(𝛼) and 𝑃_𝑎(𝛽) are the PAS and PES with 𝛼 and 𝛽 defined in Figure 1(b).

An antenna layout is a set of 𝐿 antenna positions: 𝐏={𝐩₁,𝐩₂,…,𝐩_𝐿} where 𝐩_𝑙=(𝑥_𝑙,𝑦_𝑙,𝑧_𝑙) is the 3D Cartesian coordinates of the 𝑙th antenna element. Given the PAS and PES, the antenna correlation matrix 𝐑 can be completely determined by 𝐏 through (5).

2.3. Ergodic Capacity

The instantaneous capacity, in nats per symbol, for MISO BF system is given by 𝐶=ln1+𝛾, (6) where 𝛾 is defined in (4). The ergodic capacity [31] is defined as the expectation of (6), which is 𝐶=𝐸ln1+𝛾=^∞₀[ln(1+𝑥)]𝑓_𝛾(𝑥)d𝑥, (7) where 𝑓_𝛾(𝑥) is the probability density function (p.d.f.) of 𝛾.

If the correlation matrix 𝐑 is known, then 𝑓_𝛾(𝑥) is given by [32, 33] as 𝑓_𝛾(𝑥)=^𝐿_𝑙=1𝜎_𝑙^𝐿−2∏_𝑘≠𝑙𝜎_𝑘−𝜎_𝑙−𝑥exp𝛾₀𝜎_𝑙, (8) where 0<𝜎₁<⋯<𝜎_𝑙<⋯<𝜎_𝐿 are the distinct eigenvalues of 𝐑. Substituting (8) into (7), the ergodic capacity is given as 𝐶=^𝐿_𝑙=1𝜎_𝑙^𝐿−2∏_𝑘≠𝑙𝜎_𝑘−𝜎_𝑙^∞₀()−𝑥ln1+𝑥exp𝛾₀𝜎_𝑙=d𝑥_𝐿𝑙=1𝜎^𝐿_𝑙∏_𝑘≠𝑙𝜎_𝑘−𝜎_𝑙1exp𝛾₀𝜎_𝑙𝐸₁1𝛾₀𝜎_𝑙, (9) where 𝐸₁(𝑥) is the exponential integral function, defined as 𝐸₁(𝑥)=^∞_𝑥𝑒^−𝑡𝑡d𝑡. (10) Note that (9) is valid only when the eigenvalues 𝜎₁,𝜎₂,…,𝜎_𝐿 are distinct. However, it can be easily verified that the discontinuous point of 𝐶 at 𝜎_𝑘≠𝜎_𝑚,𝑘≠𝑚 is the removable discontinuity; that is, lim_{𝜎𝑘→𝜎+𝑚}𝐶=lim_{𝜎𝑘→𝜎−𝑚}𝐶. (11) Hence, in case that 𝐑 has identical eigenvalues, the right-hand side of (9) can be evaluated with the left or right limits around the discontinuity.

From (9), we know that 𝐶 is a function of 𝐑. On the other hand, 𝐑 is given by (5). Hence, for fixed PAS and PES, 𝐶 is completely determined by the antenna layout 𝐏={𝐩₁,𝐩₂,…,𝐩_𝐿}. Our objective is to find the optimal 𝐏 which can maximize 𝐶 under the diameter constraint 𝑑_max. The optimization problem can be formally stated as 𝐏^opt=argmax_{𝐏∈𝐑𝐿×3}𝐶(𝐏)s.t.𝐩_𝑙−𝐩_𝑙′≤𝑑_max,∀1≤𝑙,𝑙^′≤𝐿,𝑙≠𝑙^′. (12) The constraint condition |𝐩_𝑙−𝐩_𝑙′|≤𝑑_max implies that the antenna elements should be placed within a sphere of diameter 𝑑_max.

3. The Optimal Antenna Layout

The optimization problem shown in (12) is too difficult to be solved directly, so in this section, we use simulated annealing (SA) [24–26] to search for the optimal transmit antenna layout. A brief description on SA and the related SA parameters used in this paper are presented in the Appendix.

We consider the 8×1 MISO BF system. With 𝐿=8 antennas, some of the regular antenna placement patterns are illustrated in Figure 2 for reference.

Figure 2 

Reference regular antenna layouts.

The optimal antenna layout will certainly depend on the angle spreading of the radio wave. For simplicity, we consider the uniform distribution PAS and PES, which are defined as 𝑃_𝑒(𝛼)=12Δ_𝛼𝛼₀−Δ_𝛼≤𝛼≤𝛼₀+Δ_𝛼𝑃_𝑎𝛽=12Δ_𝛽𝛽₀−Δ_𝛽≤𝛽≤𝛽₀+Δ_𝛽, (13) where Δ_𝛼 and Δ_𝛽 are the range of 𝛼 and 𝛽⋅𝛼₀ and 𝛽₀ denote the center of angle of departure (AoD). Without loss of generality, we set 𝛼₀=𝛽₀=0. Substituting (13) into (5), the spatial correlation of transmit antenna is given as 𝑅𝐩,𝐩^′=14Δ_𝛼Δ_𝛽×^Δ𝛼_Δ𝛼^Δ𝛽_Δ𝛽exp𝑗2𝜋𝜆𝑥−𝑥^′+cos𝛼cos𝛽𝑦−𝑦^′+cos𝛼sin𝛽𝑧−𝑧^′sin𝛼d𝛼d𝛽. (14)

Five propagation scenarios are considered in this paper, as shown in Table 1. Scenario I is the 3D whole space uniformly scattering environment. In Scenario II, the elevation of the departing rays is compressed to a very narrow domain, the scattering range converges like a 2D thin plane. Scenario III represents the narrow spreading situation where the angle spreading of departing rays for both azimuth and elevation has been compressed in the ±10° range. In this case, the departing rays converge like a thin beam. Scenario IV and Scenario V are the degraded cases of Scenario II and Scenario III.

Table 1 

PAS and PES for different propagation scenarios.

3.1. The Optimal Layouts for Scenario I

Figure 3 shows the optimal antenna layout for Scenario I with diameter constraint as 𝑑_max=0.2𝜆,  𝑑_max=0.6𝜆, and 𝑑_max=𝜆, respectively. For each 𝑑_max value, two subfigures are presented which are the two different side views of the same antenna placement layout. An interesting observation is that the optimal antenna layout will always place the antenna elements on the surface of the sphere.

(a) 𝑑max=0.2𝜆

(b) 𝑑max=0.6𝜆

(c) 𝑑max=𝜆

Figure 3 

The optimal layout for Scenario I.

In Figure 4, we compare the ergodic capacity of the optimal layout with that of regular layouts listed in Figure 2. We can see that the cube and the cone layout are near optimal. Specifically, when 𝑑_max≤0.4𝜆, the ergodic capacity of cube layout is almost the same as that of optimal layout and the cone layout is very close to optimal. When 𝑑_max>0.4𝜆, the cone layout is more close to the optimal layout than the cube layout.

Figure 4 

Capacity difference between the optimal layout and the reference regular layouts for Scenario I with 𝛾0=20dB.

3.2. The Optimal Layouts in Scenario II

Figure 5 shows the optimal antenna layout for Scenario II where the angle spread is uniform in the azimuth plane and is narrow in elevation (constrained in ±10° in elevation plane). Comparing Figures 5 and 3, we can see that the angle spread of the departing rays has serious impact on the optimal antenna layout. Since the rays mainly depart from the horizontal plane, the optimal layout is essentially a ring with its normal line pointing in the 𝑧-direction. Figure 5(b) is an exception where 7 antennas have formed a ring in horizontal plane while the 8th antenna is placed at the norm line of the ring.

(a) 𝑑max=0.2𝜆

(b) 𝑑max=0.6𝜆

(c) 𝑑max=𝜆

Figure 5 

The optimal layout for Scenario II.

If we force the 8th antenna in Figure 5(b) at the norm line of the ring to the center of the horizontal plane ring and keep the positions of the rest 7 antennas, we will arrive at a layout similar to the center-circle layout shown in Figure 2. Figure 6 shows the capacity difference between the layout in Figure 5(b) and the center-circle layout. The results show that the capacity difference is very slight, since the angle spread in elevation is narrow (constrained in ±10° in elevation plane). In other words, although the optimal layout suggested by SA is Figure 5(b) for 𝑑_max=0.6𝜆, which is quite different with Figures 5(a) and 5(c) where the optimal solution is a plane, we know from Figure 6 that forcing the highest point to the 𝑥-𝑦 plane will lead to almost no loss. What is implicit here is that, with a very narrow angel spreading at elevation, the optimization can focus only on the 𝑥-𝑦 plane. In fact, we will see later that if the angle spread in elevation is even smaller as to zero, the 𝑧-axial coordinate will no more affect the performance.

Figure 6 

Capacity difference between the layout in Figure 5(b) and the corresponding center-circle layout. 𝛾0=20dB.

In Figure 7, we compare the ergodic capacity of the optimal layout under Scenario II with that of regular layouts listed in Figure 2. Unsurprisingly, the uniform circle layout and the center circle layout are near optimal. Specifically, when 𝑑_max≤0.5𝜆, the ergodic capacity of the uniform circle layout is almost the same as that of optimal layout. When 𝑑_max>0.5𝜆, the center circle layout is more close to the optimal layout than the uniform circle layout.

Figure 7 

Capacity difference between the optimal layout and the reference regular layouts for Scenario II. 𝛾0=20dB.

3.3. The Optimal Layouts in Scenario III

Figure 8 shows the optimal antenna layout for Scenario III where the angle spread is narrow in both azimuth and elevation (constrained in ±10° in azimuth and elevation plane). Comparing Figures 8, 5, and 3, we can see that the optimal antenna layout dimension is highly dependent on the power azimuth spectrum (PAS) distribution and power elevation spectrum (PES) distribution. In this case, the departing rays converge to a thin beam which is approximately a single plane wave. The optimal layout in this case is different with all the results in the two scenarios mentioned previously, the antennas no longer distribute on the sphere surface. The eight antennas are some like eight nodes on a knotted rope, but each of them has its own vibration amplitude and these amplitudes are irregular.

(a) 𝑑max=0.2𝜆

(b) 𝑑max=0.6𝜆

(c) 𝑑max=𝜆

Figure 8 

The optimal layout for Scenario III.

In Figure 9, we compare the ergodic capacity of the optimal layout under Scenario III with that of regular layouts listed in Figure 2. We can see that when the angle spreading is very narrow both in elevation and azimuth, the uniform linear array is the one closet to the optimal, among the layouts listed in Figure 2.

Figure 9 

Capacity difference between the optimal layout and the reference regular layouts for Scenario III. 𝛾0=20dB.

3.4. The Optimal Layouts in Scenario IV

Scenario IV is the degraded case of Scenario II. The range of PES has been further compressed by setting Δ_𝛼 in (13) to zero, resulting in that the 𝑃_𝑒(𝛼) becomes the Dirac impulse function and the double integral in (14) simplifies to a single integral as shown in (15). The 𝑧-axial coordinate disappears from (15) implying that the 𝑧-axial coordinate will no more affect the antenna spatial correlation. The 3D antenna layout optimization problem is then reduced to a 2D optimization problem. With a constraint that the antenna elements must be placed within a circle of diameter 𝑑_max: 𝑅𝐩,𝐩^′=12Δ_𝛽^Δ𝛽_Δ𝛽exp𝑗2𝜋𝜆𝑥−𝑥^′+cos𝛽𝑦−𝑦^′sin𝛽d𝛽. (15)

Figure 10 shows the optimal antenna coordinates for Scenario IV which is the output of SA. The points in red, blue or black denote, respectively, the optimal antenna locations for different 𝑑_max constraints. In case that 𝑑_max=0.2𝜆 and 𝑑_max=𝜆, the antennas are almost uniformly distributed on a circle while in the case of 𝑑_max=0.6𝜆, one antenna is placed at the center of the circle and the rests are nonuniformly distributed on the circle.

Figure 10 

The optimal layout for Scenario IV with various 𝑑max.

3.5. The Optimal Layouts in Scenario V

Scenario V is the degraded case with both Δ_𝛼 and Δ_𝛽 set to zero. Now the double integral in (14) disappears and it reduces to (16). This indicates that the optimal layout does not depend on 𝑧-axial and 𝑦-axial coordinates and the original 3D antenna layout optimization problem is reduced to a 1D optimization problem with a constraint that all antennas should be placed in a interval of length 𝑑_max: 𝑅𝐩,𝐩^′=exp𝑗2𝜋𝜆𝑥−𝑥^′. (16)

Figure 11 is the optimal layout obtained through SA. The points in red, blue, or black correspond to three difference 𝑑_max constraints. We can see that the optimal layout in this case is the nonuniform linear array.

Figure 11 

Optimal layout for Scenario V with various 𝑑max.

In summary, in the three-dimensional (3D) space uniformly scattering environment, the optimal antenna layout will locate the antenna elements on the surface of a sphere. For practical consideration, the cube layout is a good choice for 𝑑_max≤0.4𝜆 and the cone layout is better for 𝑑_max>0.4𝜆. When the scattering is concentrated in a two-dimensional (2D) plane, the circular layout is the optimal. For practical consideration, we may consider the uniform circular layout for 𝑑_max≤0.5𝜆 and the center circle layout for 𝑑_max>0.5𝜆. While the scattering is very narrow in both elevation and azimuth, the optimal layout is nearly a nonuniform linear array. The uniform linear layout is near optimal particularly when 𝑑_max>0.4𝜆.

4. Conclusions

This paper investigates the problem of arranging the transmit antennas to attain the maximum ergodic capacity in an MISO beamforming system. Since the optimization problem cannot be solved analytically, we use simulated annealing (SA) algorithm to search for the optimal antenna layout. The results show that the optimal antenna layout dimension is highly dependent on the power azimuth spectrum (PAS) distribution and power elevation spectrum (PES) distribution.

Appendix

Simulated annealing [24–26] is a technique for combinatorial optimization problems, such as minimizing functions of many variables. For our purposes, a combinatorial optimization problem is one in which we seek to find some configurations of antenna layout 𝐏={𝐩_𝟏,𝐩_𝟐,…,𝐩_𝐿} in (12) that minimizes the 𝑓(𝐏)=−𝐶 in (9), which is defined as the cost or objective function of SA in this paper. SA starts with a random or feasible initial solution and moves step by step towards a solution 𝐒^opt giving hopefully the minimum (or close to the minimum) of the cost function.

When the configuration of antenna layout 𝐏={𝐩_𝟏,𝐩_𝟐,…,𝐩_𝐿} is at a feasible solution, we attempt some small random perturbation to the configuration that yields a nearby solution. This process can continue starting from the new configuration until no further improvements are obtained, at which point the process terminates. This strategy seems reasonable, but it has a serious problem: it is easily trapped in local minimum, solutions that look good in some small neighborhood of the cost surface but are not necessarily the global optimum. To avoid this problem, the acceptance function is used to decide the acceptance probability.

The probability of accepting a move that causes an increase Δ𝑓 in 𝑓(𝐏)=−𝐶 is determined by the acceptance function. Like most researchers, we use the Metropolis acceptance function [24] exp(−Δ𝑓/𝑇), where 𝑇 is a control parameter corresponding to the temperature in analogy with physical annealing. Because Δ𝑓>0 and 𝑇>0, exp(−Δ𝑓/𝑇) is always smaller than 1. We generate a random number 𝑐 uniformly distributed on the interval [0,1] and compare exp(−Δ𝑓/𝑇) with 𝑐. If exp(−Δ𝑓/𝑇)>𝑎, then the move is accepted; otherwise, the move is rejected.

The initial temperature 𝑇₀ is set to be 100, and the ending temperature 𝑇_end equals 0.0001. The cooling schedule will be the simplest possible: 𝑇_new=𝑏𝑇_old,0<𝑏<1, where the initial temperature and cooling rate 𝑏 are determined empirically to give good results and we let 𝑏 be 0.9 in this paper. At each temperature, we perform 100 moves. The stopping criterion is to terminate annealing when the cost improvement seen across five successive temperatures or the temperature is cooler than the ending temperature. The key SA parameters in this paper are listed in Table 2 and the detailed process of SA is shown in Algorithm 1.

Table 2 

Parameters for SA algorithms.

(1) Generate a random initial configuration of antennas 𝐏𝟎, and 𝐏opt=𝐏𝟎, compute the cost function value 𝑓(𝐏0); (2) Set the initial temperature 𝑇(0)=𝑇0 and the iterations 𝑖=1; (3) Do while 𝑇(𝑖)>𝑇end  (1) for inner loop 𝑚=1∼𝐾  (2) Do some small random perturbation to the 𝐏opt that yields a nearby solution 𝐏new. Compute 𝑓(𝐏new)  and Δ𝑓=𝑓(𝐏new)−𝑓(𝐏opt)  (3) if Δ𝑓≤0, 𝐏opt=𝐏new;  (4) if Δ𝑓>0, 𝑢=exp(−Δ𝑓/𝑇(𝑖)); generate a random variable 𝑐=random[0,1], if 𝑐<𝑢, 𝐏opt=𝐏new.  (5) end for (4) 𝑖=𝑖+1; (5) 𝑇(𝑖+1)=𝑇(𝑖); (6) End Do while (7) Output 𝐏opt

Algorithm 1 

Process of SA algorithms.

Acknowledgment

This paper was supported by National Natural Science Foundation of China (no. 61072059).