1. Introduction

The perfectly matched layer (PML) [1], a widely used method for truncating unbounded domains in numerical simulations, performs well in numerical simulations of wave propagation problems, which has been applied in the beam propagation method and the mode matching method in the optical propagation [2–4]. Meanwhile, the PML has been successfully applied to many other numerical methods including the the finite-difference-time-domain (FDTD) [5–8], the finite-element-method (FEM) [9,10], etc. In [5], Juntunen, Kantartzis and Tsiboukis presented a closed-form reflection coefficient for the PML, when realized using the FDTD algorithm, and they found that zero reflection can be obtained for isolated pairs of frequency and angle of incidence. In [6], Prokopidis, Kantartzis, and Tsiboukis provided a robust closed-form expression of the reflection coefficient concerning the discrete uniaxial FDTD/PML for the termination of lossy materials, and the proposed technique greatly diminishes artificial reflections and therefore achieves a significant accuracy. In [7], Richard and Nikolova used a problem independent approach to enhance the PML performance within the whole frequency band of excitation, in the presence of both evanescent and propagating fields. In [8], Laudani, Calvagna, Sorbello, Janner and Tramontana proposed a simple framework for a deeper understanding of approximation errors induced by discretization, leading to a closed form expression to compute numerical reflection of arbitrary PML profiles, which also features reliable prediction of the error. In [9], Pekel and Mittra modified and extended the PML to the frequency domain for FEM applications, and the governing equations require neither the splitting of the field components of interest, nor do they involve negative conductivity parameters.

It is well known that a PML can be regarded as a complex coordinate stretching [11] that turns the outgoing wave into an exponentially decaying solution. Theoretically, applying the PML can completely absorb outgoing waves of any angle and frequency, and also produce no reflection. However, due to the limited thickness of the PML, certain reflected waves will be generated.

Besides, the Helmholtz equation can be discretized by the finite difference method, and the corresponding average discrete reflectivity is computed after using a perfect matching layer to the edge of the solution region [12]. The optimal parameters of the PML can be chosen by minimizing an object function. The object function often adopts the discrete reflection coefficient [13], the sum of the weighted reflection coefficient [14,15], the reflection and computation time [16], the the average discrete reflectivity [17], or the sum of the reflection coefficient and the truncation error [18]. However, above methods limit the discrete degree of the equation system.

When the thickness of the PML is given, in order to reduce the reflection, the integral value of the absorption function should enlarge. In order to maintain the accuracy of numerical computation, the requirement to reduce the discrete step size should be satisfied, which easily contributes to large rounding error. The coefficient matrix of the discrete equation system will gradually become ill-conditioned, so that the rounding error will seriously impact the accuracy of numerical solution of the equation system.

In the present paper, we study the optimal parameters’ selection of the PML by minimizing both the average discrete reflectivity and the rounding error. The structure of the paper is as bellow. In section 2, we come up with two methods to find optimal PML parameters. One is used by minimizing the average discrete reflectance, based on the first, and the other adds rounding error to the model, which we are the first to propose. Section 3 demonstrates numerical simulations of light propagation using two optimal PML parameters. Finally, the conclusions are given in Section 4.

2. Mathematical models for optimal parameters of the perfectly matched layers

2.1 Reflectivity of a PML

Consider transverse electric waves in a two-dimensional (2-D) waveguiding structure shown in Fig. 1. The y-component of the electric field satisfies the Helmholtz equation. We have:

 

Fig. 1. An optical waveguide terminated by a perfectly matched layer (PML).

Download Full Size | PDF

For $D$ and $H$ satisfying $D<H<M$, we define the absorption function $\sigma (z)$ such that $\sigma (z)=0$ for $z \geq H$ and $\sigma (z)>0$ for $z<H$. Then, replace $z$ by $\hat {z}=\int _{0}^{z} s(\tau ) \,\mathbf {d}\tau$, where $s(z)=1+\mathrm {\mathbf i}\ \sigma (z)$ and i=$\sqrt {-1}$. Under this transformation, we have $\frac {\textstyle \partial }{\textstyle \partial z}\rightarrow \frac {\textstyle \partial }{\textstyle \partial \hat z}=\frac {\textstyle 1}{\textstyle 1+\mathrm {\mathbf i}\sigma (z)}\cdot \frac {\textstyle \partial }{\textstyle \partial z}$.

Then, the original Helmholtz equation is transformed to the following form:

At $z=D$, we use a simple zero condition: $u=0$. The actual PML is the layer $D\le z\le H$ where Eq. (2) is different from the original Helmholtz equation.

For $H<z<M$, the solution to Eq. (2) is:

Denote $\theta$ as the angle between the $x$ axis and the wave vector $(-\alpha,\beta )$, where $\alpha =k_0n_0\sin \theta$ and $\beta =k_0n_0\cos \theta$, respectively. Obviously, $|R|$ depends on the angle of incidence $\theta$. In the special case, when the waves are almost parallel to the $x$ axis, $|R|$ is close to 1. For a fixed $\theta _0$, $|R|$ depends only on $\int _{D}^{H} \sigma (\tau )\,\mathbf {d}\tau$. Thus, we can make $|R|$ arbitrarily small by taking $\int _{D}^{H} \sigma (\tau )\,\mathbf {d}\tau$ large enough [17]. However, this is only a theoretical analysis. In fact, when the depth variable $z$ is discretized, the reflection coefficient depends on the actual distribution of the absorption function $\sigma (z)$ and it cannot be reduced indefinitely by increasing $\sigma (z)$[12].

For Eq. (2), consider a second-order finite difference approximation (central difference quotient) in the transverse direction $z$. The grid size $h= \varDelta z=z_j-z_{j-1}$, $z_0=D$, and $z_m=H=D+m\cdot h$. Then

When $H\leq z\leq M$, $\sigma (z)=0$, the coefficients are simplified to $a_j=c_j=1$ and $b_j=-2$. Therefore Eq. (7) becomes

For $H\leq z\leq M$, we have $\alpha ^2+\beta ^2=k_0^2n_0^2$, and

By using the previous Eqs. (9) and (10), combining with $\alpha ^2+\beta ^2=k_0^2n_0^2$, we can introduce

The above Eq. (11) can be used to solve the reflection coefficient. Consider a zero boundary condition: $u_0=0$, we write down Eq. (7) for $j=1,2,\ldots,m$ as follows.

From above, $\alpha$ clearly affects the value of the reflection coefficient $R$. In addition,

Then, $R$ can be considered as a function of $\theta$. When taking a fixed $\theta _0$, the corresponding reflection coefficient $R$ can be found. Notice that $R$ depends on $\alpha$, and the perfectly matched layer is the ideal material which we envision, so it is not proper to consider from a fixed $\theta _0$. Thus, we consider the average reflection coefficient module. In other words, the parameters in the absorption function should be selected appropriately so that the average reflection coefficient module is minimized for optimization purposes. The average reflection coefficient module is defined as

Since $|R(\theta )\vert$ is in the form of a column of points $\{ (\theta _i,|R(\theta _i)|),i=0,1,2,\ldots,n_{\theta }\}$, where $n_{\theta }$ is the number of dissection cells of $\theta$. Equation (19) is approximately computed by the Composite Simpson Integral Formula [19]:

2.2 Behavior of the difference equations for finding the discrete reflectivity

We speculate that $m$ will keep growing when $N=\lambda _0/(n_0h)$ obtains relatively large ($m=\frac {\delta }{h}=\frac {\delta }{\lambda _0/(n_0N)}=\frac { n_0N\delta }{\lambda _0}$, where $\delta$ is the thickness of the PML), where $N$ is represented as the number of grid points per wavelength. Then the order of the coefficient matrix (the coefficient matrix is denoted as $A$) in Eq. (16) becomes big enough, making it possible for the coefficient matrix $A$ to become ill-conditioned. For further verification, we study the condition number (the condition number is defined as $||A||_2\cdot ||A^{-1}||_2$) [20] of the coefficient matrix $A$ in Eq. (16) with respect to $N$.

For simplicity, we let $\lambda _0=n_0=1$, and

When $p=3$, $\theta = \frac {\textstyle \pi }{\textstyle 3}$, as shown in Fig. 2, the condition number of the coefficient matrix of Eq. (16) keeps increasing as $N$ goes up. For example, when $D=-3.25\mu \mathrm {m}$, $\delta = 0.25\mu \mathrm {m}$, $N=400$, and then $m=100$, $h=0.0025\mu \mathrm {m}$, the corresponding optimal parameters $S=S_1=396.5861$, under which the corresponding average discrete reflectivity $\overline {|R|}=2.34755\times 10^{-3}$, based on which, when we fix $\theta =\frac {\textstyle \pi }{\textstyle 3}$, the condition number of Eq. (16) is approximately $21545$. Obviously, the equation system in Eq. (16) is a sickness equation, so a small change in the rounding error will result in a deviation in the calculation of the average discrete reflectivity. Therefore, we can safely draw a conclusion that when $N$ is relatively large, the coefficient matrix is likely to be an ill-conditioned matrix, implying that despite a small rounding error, the calculated average discrete reflectivity will suffer greatly.

 

Fig. 2. Relation of the condition number of the coefficient matrix of Eq. (16) to $N$.

Download Full Size | PDF

2.3 Round-off error of the difference equation

In the previous section, we approximate the second-order derivative using the second-order difference formula, which leads to errors containing both truncation errors and rounding errors. Since the truncation error (about $O(h^2)$) is reduced with the discrete step size $h$, when the parameter $S$ is large enough, the optimal step size $h$ is so small at this time that the truncation error would be negligible. In this paper, we mainly consider the effect of rounding errors.

Denote $\mathrm {fl}(*)$ as the floating-point input of $*$ in the computer, and $t$ as the floating-point number of the commonly used binary computer, generally $t=32$. For any number $a$, we have $\mathrm {fl}(a)=a(1+\iota )$, where $\iota$ is real number and $|\iota \vert \leq \mu$($\mu =2^{-t}$) [19].

For the $j$-th difference equation of the linear system, namely Eq. (16), we obtain its round-off error, denoted as $Error_j$, as follows:

Similarly, we can obtain the upper bounds of round-off errors of other difference equations as follows.

For the $j-$th equation of the linear system (16), we obtain its upper bound of the round-off error, denoted as $Error(\theta )_{j}$, $j\in \{1,m-1,m\}$ as follows:

To sum up, we can get the upper bound vector of rounding error: $[ Error(\theta )_1,\ldots,Error(\theta )_m]^T$ of Eq. (16). We let $Error(\theta )=\displaystyle \sum _{j = 1}^{m} Error(\theta )_j$. In the process of minimizing the average discrete reflectivity, we discretize $\theta$, where then use the compound Simpson integral formula to solve $\overline {|R|}$, and the discrete fraction is $n_{\theta }$. Then, let the total rounding error be:

Obviously, with the continuous increasing of the discrete lattice number $m$, the total rounding error of Eq. (16) will become larger and larger, and thus affect the accuracy of the average discrete reflectivity $\overline {|R|}$. From Eq. (21) and Eq. (22), we can see $|R|$ depends only on $S|H-D\vert$. Furthermore, we give two models to find the optimal parameters of PMLs in the bellow.

Model (I): PML parameter selection based on minimizing average discrete reflectivity [17], that is,

3. Numerical simulations

3.1 Choice of optimal perfectly matched layer by Model I

In this section, we use the Model I to find the optimal parameters. In this discrete case, the choice of $\sigma _0$ is important. We study the convergence of the variation of $\overline {|R|}$ with $n_{\theta }$. Take a PML with thickness of five cells $(m=5)$ as an example, where the grid size $h$ is $\frac {\textstyle 1}{\textstyle 20}$ of the wavelength in the medium $(h=0.05\mu \mathrm {m})$. Precisely, we have $\lambda _0 =1\mu \mathrm {m}$, $k_0=\frac {\textstyle 2\pi }{\textstyle \lambda _0}$, and $n_0=1$. To further determine the value of $n_{\theta }$ in the actual calculation of (20), we should fix $p=3$, and $S=25$. The variation of $\overline {|R|}$ with $n_{\theta }$ is shown in Fig. 3. From Fig. 3, we notice that, when $n_{\theta } \geq 70$, $\overline {|R|}$ is gradually stable ($|\ \overline {|R|}^{( n_{\theta }+1)}-\overline {|R|}^{(n_{\theta })}\ |<10^{-5}$, $n_{\theta }>70$). Therefore, we can fix $n_{\theta }=70$. By using Eq. (16), we can further solve $\overline {|R|}$ corresponding to different powers $p$ and parameters $S$. The optimal parameters $S$ corresponding to different powers $p$ are solved below by Model I ($p$ is given).

 

Fig. 3. The variation of $\overline { |R|}$ with $n_{\theta }$ $(m=5)$.

Download Full Size | PDF

Since the locally optimal solutions to $\underset {S\ge 0}{\min }\overline {|R|}$ are hard to obtain, over the entire domain of positive real numbers, the genetic algorithm is applied, which is more convenient than the traditional numerical solution method with randomly taken points. As given in Fig. 4, when $p=3$, and $h=0.05\mu \mathrm {m}$, the best parameters $S$ corresponding to $\underset {S\ge 0}{\min }\overline {|R|}$ are concentrated in the interval $[20,30]$, which indicates that the solution to this problem optimized with the genetic algorithm is more stable, and the accuracy of the solution is improved.

 

Fig. 4. Distribution of population individual sizes (parameters $S$) after the last iteration $(p=3)$.

Download Full Size | PDF

The average discrete reflectivity $\overline {|R|}$ and the optimal values of $S$ are listed in Table 1 for $p=0,1,2,3,4,5$. Table 2 is the results of [17] on the selection of optimal PML parameters. Comparing Tables 1 and 2, we see that our results calculated by the genetic algorithm, and the other results listed in the Ref. [17], are basically consistent. The results shows that when $p=2,3,4,5$, the average reflection coefficient $\overline {|R|}$ is quite similar. In order to further analyze the details of the PML under different $p$, we plot $|R|$ against $\theta$ in Fig. 5. From Table 1 and Fig. 5, it appears that when $p=3$, the overall effect of the absorption function performs the best. Also, when $\theta$ takes different values, the respective results for $p=2,3,4,5$ vary widely. The quadratic profile gives the lowest reflectivity for large angle ($\theta >0.1\frac {\textstyle \pi }{\textstyle 2}$). When $\theta < 0.05\frac {\textstyle \pi }{\textstyle 2}$, the cubic profile enjoys the lowest reflectivity and the quintic profile is the worst. Fig. 6 demonstrates the discrete reflectivities for $p=2,3,4$ and $h=0.05\mu \mathrm {m}$ in [17]. Comparing Figs. 5 and 6, the effect of PML absorption under the PML parameters selected by the genetic algorithm is basically consistent with the results in the reference [17].

Table 1. Average reflectivity $\overline {|R|}$ and optical scaling parameter $S$ for PML

View Table | View all tables in this article

Table 2. Results of the selection of optimal PML parameters in [17]

View Table | View all tables in this article

 

Fig. 5. Discrete reflectivities for $p=2,3,4,5$, using the optimal values of $S$ given in Table 1.

Download Full Size | PDF

 

Fig. 6. Discrete reflectivities for $p=2,3,4$ in [17].

Download Full Size | PDF

To illustrate the dependence on the scaling parameter $S$, we fix $p=2$, and plot the discrete reflectivity $|R(\theta )|$ for $S=16.0177$, $S=18.0177$, $S=20.0177$ and $S=22.0177$ in Fig. 7. We notice that the four curves show the almost same trend. When $\theta <0.14\frac {\textstyle \pi }{\textstyle 2}$, the value of $|R(\theta )|$ decreases as the value of $S$ increases. However, when $\theta >0.15\frac {\textstyle \pi }{\textstyle 2}$, $|R(\theta )|$ decreases with smaller value $S$. The optimal value of $S$ (i.e., $S=18.0177$) gives a reasonable compromise, as it minimizes the average reflectivity. Figure 8 shows the figure of the average discrete reflectance of the quadratic profile for $S=16$, $S=18$ and $S=20$ in [17].The comparison shows that the experimental results are basically consistent.

 

Fig. 7. Discrete reflectivity for $p=2$ and $S=16.01778$, $S=18.0177$, $S=20.0177$ and $S=22.0177$.

Download Full Size | PDF

 

Fig. 8. Discrete reflectivity for $p = 2$ and $S = 16$, $S = 18$, and $S = 20$ in [17].

Download Full Size | PDF

It is clear from the above proof that the choice of the optimal parameter $S$ is only related to the degree of dispersion $m$ and $N=\frac {\textstyle \lambda _0}{\textstyle n_0h}$. We use the genetic algorithm to calculate optimal parameters $S$ when $5<N<25$, $5<m<15$, $p=3$. As shown in Fig. 9, the optimal parameter $S$ does not always increase respect to $m$. The curves for $m=6$ and $m=9$ share some unexpected behavior. However, the optimal parameter $S$ always increases respect to $N$.

 

Fig. 9. Dependence of the optimal scaling parameter $S$ on $N$ for the cubic profile and for different $m$.

Download Full Size | PDF

It can be expected that the optimal parameter $S$ goes up as $N$ increasing, and thus the computational effort of the computer increases. Furthermore, we investigate the optimal parameter $S$ and the average discrete reflectivity $|R(\theta )|$ versus $N$ when fixing the thickness of the PML. In Figs. 10 and 11, the variation of the optimal parameters $S$ are consistent with our expectation, although $|R|$ decreases with $N$ increasing. $|R|$ gradually tends to $0.00234$, and the decreasing rate becomes small.

 

Fig. 10. Dependence of the optimal scaling parameter $S$ on $N$ for the cubic profile $(m=5)$.

Download Full Size | PDF

 

Fig. 11. Dependence of $|R|$ on $N$ for the cubic profile $(m=5)$.

Download Full Size | PDF

3.2 Choice of optimal perfectly matched layer by Model II

In this section, we use the Model II to find the optimal parameters. The general approach to solve the multi-objective problem in Eq. (25) is to transform it into a single-objective problem. Using a simple weighting method here, the objective function becomes

Table 3. Comparison of the results of two different PML parameter selection methods

View Table | View all tables in this article

From Table 3, when $m$ is small ($m\leq 40$), the solution of Eq. (16) is more stable due to the lower discretization, so that $S_1$ and $S_2$ are almost the same even after considering the rounding error. When $m$ is larger ($m>40$), the condition number of the coefficient matrix of Eq. (16) becomes progressively larger due to the progressively larger discretization, resulting in rounding errors which begin to have an impedance on the numerical results of Eq. (16), and thus the results of $S_1$ and $S_2$ are different.

In the next section we compare the absorption effects of two different PML parameters in the actual propagation situation.

3.3 PML actual effect simulation

(A) Wave propagation with PML Let $M=0\mu \mathrm {m}$, $H=-3\mu \mathrm {m}$, $D=-3.25\mu \mathrm {m}$,and $L=5\mu \mathrm {m}$. More precisely, we have $\lambda _0=1\mu \mathrm {m}$, $k_0=\frac {2\pi }{\lambda _0}$, and $n_0=1$. Thus, we have:

In the following, we study the propagation of waves from $x=0$ to $x=L$ (TE case) with the following initial conditions, boundary conditions and the PML interface conditions:

For weakly range-dependent waveguides, the Helmholtz equation is often approximated by the following one-way Helmholtz equation [21]:

In particular, when $k_0$ and $n_0$ are $x$-independent, the one-way equation $u_x= \mathrm {\mathbf i}\sqrt {\partial _z^2+k_0^2n_0^2}u$ is equivalent to the Helmholtz equation. For $x\in [0,L]$, we discretize the $x$-axis: $x_1=0$, $x_p=p\Delta x$, where $p=1,2,\ldots, \tilde {N}$ and $\Delta x=\frac {\textstyle L}{\textstyle \tilde {N}}$. Consider a step from $x_p$ to $x_{p+1}$, the solution to the one-way Helmholtz equation can be approximated by

The above stepping algorithm involves the problem of second-order partial derivative $\sqrt {\partial _z^2+k_0^2n_0^2}$. Since our problem involves the PML with two different discrete steps, we discretize the $z$-axis in segments.

When $z\in [D,H]$, the discrete step is $h_1=\frac {\textstyle H-D}{\textstyle m}$, and step number is $m$, so $z_q=D+qH_1\in [D,H]$, $q=1,2,\ldots,m$.

When $z\in [H,M]$, the discrete step is $h_2=\frac {\textstyle M-H}{\textstyle n}$, and step number is $n$, so $z_q=H+(q-m)h_2\in [H,M]$, $q=m+1,m+2,\ldots,m+n$.

From Eq. (28), we have

Except for the special difference formula of the second derivative at the point $u(x,z_m)$ on the PML boundary, the second derivative of other points can be directly solved by the second-order difference formula:

For $u(x,z_m)$, we have

Plug Eq. (32) into Eq. (37), we obtain

In summary, on the weird $[x_p,x_{p+1}]$, operator $\partial_{\hat z}^{2}$ can be approximated as discrete matrix $D_2$:

Therefore, after plugging Eq. (39) into Eq. (29), the propagation result of the wave from $x=0$ to $x=L$ in the waveguide can be obtained.

(B) Exact solution of wave propagation

In order to compare the propagation effect of the PML, we use the separation parameter method to find the exact solution of the wave propagation in vacuum by expanding the depth $|D|$ to $|\tilde {D}|$. Making use of the following conditions:

More specifically,

Next, we take $L=5\mu \mathrm {m}$ and plot the propagation results for $\tilde {D}=-25.25\mu \mathrm {m}$ and $\tilde {D}=-30.25\mu \mathrm {m}$, respectively, as shown in Fig. 12. It can be seen that the actual solutions almost coincide in interval $[-3,0]$ as $\tilde {D}=-25.25\mu \mathrm {m}$ decreases to $\tilde {D}=-30.25\mu \mathrm {m}$ and the average relative error between the two curves is $5.69\times 10^{-4}$. Therefore, the approximate solution Eq. (40) with $\tilde {D}=-25.25\mu \mathrm {m}$ is used as the “exact solution" (called as reference solution).

(C) Comparison between the approximate and the actual solutions

 

Fig. 12. “Exact solutions” corresponding to $\tilde {D}=-25.25\mu \mathrm {m}$ and $\tilde {D}=-30.25\mu \mathrm {m}$.

Download Full Size | PDF

Consider the case with $\Delta x=0.01\mu \mathrm {m}$ and $n=600$, the incident wave is $u(0,z)=u_0(z)=\sin (\pi z)\sin (\frac {\textstyle 4}{\textstyle 13}\pi z)$. When $m=100$, the corresponding PML parameters $S=S_1=396.5861$ and $S=S_2=1203.64$ for Eqs. (24) and (26), respectively. We use Eq. (40) as the “exact solution” (called as the reference solution). The concrete comparison of the module of optical propagation solutions are shown in Fig. 13, from which the module of the corresponding propagation solutions for two different PML parameters are pretty close to the “exact solution" in the non-PML region $[-3,0]$. To further compare the two different PML parameters, we average the relative error of the value of the wave function at each discrete point between the approximate solution and the exact solution. Denote the average relative error of the module at $S=S_1=396.5861$ be $E_1$ and the average relative error of the module at $S=S_2=1203.64$ be $E_2$. $E_1=0.00469$ and $E_2=0.00338$ can be obtained by calculation. It can be seen from the comparison that the PML absorption effect of the PML parameter selection considering the rounding error will be better. In the meanwhile, the relative errors of the propagated and exact solutions for two different PML parameters when the discretization degree $m\geq 50$ are shown in Table 4. From Table 4, when the degree of discretization ($m$) is large enough, the rounding error does affect the selection of the best PML parameters. Therefore, PML absorption effect will be better when choosing PML parameters which take into account rounding errors.

Table 4. Relative errors of the modules of the approximate solutions at $x=L$ with two different PML’s parameters

View Table | View all tables in this article

 

Fig. 13. Comparison of the modules between two approximate solutions and “exact solution” for different PML parameters at $x=L=5\mu \mathrm {m}$($m=100$).

Download Full Size | PDF

4. Conclusion

This paper mainly focuses on the optimal design of perfect matching layer based on two-dimensional Helmholtz equation and on the parameter optimization of absorption function $\sigma (z)$. The equation is discretized by finite difference method, and then the corresponding average discrete reflectivity and the rounding error are calculated after adding a perfect matching layer on the edge of the region. The optimal parameters $S$ are obtained by minimizing the average discrete reflectivity $\overline {|R|}$ and the rounding error. Then, the propagation effect of PML is simulated by the solution of one-way equation, and then compared with the exact solution of Helmholtz equation. The research shows that when the degree of dispersion $m$ is relatively large, the condition number of equation coefficient matrix used to solve the average discrete reflectivity becomes very large, which shows that the increase of the degree of dispersion has affected the calculation accuracy of the average discrete reflectance. Given the grid size, the free-space wavelength, the reflective index, and the thickness of the PML, we can solve the optimal parameter $S$ by minimizing the objective function $0.7*f_1(S)+0.3*f_2(S)$. The results show that when the degree of dispersion $m$ is small enough, the results of $\displaystyle \min _{S\ge 0}\{0.7*f_1(S)+0.3*f_2(S)\}$ and $\displaystyle \min _{S\ge 0}f_1(S)$ are almost the same. When the degree of dispersion $m$ is large, the result of $\displaystyle \min _{S\ge 0}\{0.7*f_1(S)+0.3*f_2(S)\}$ is obviously better than that of $\displaystyle \min _{S\ge 0}f_1(S)$. The further research on the choice of weights is in progress.