A Hybrid Di ﬀ ractive Optical Element Design Algorithm Combining Particle Swarm Optimization and a Simulated Annealing Algorithm

: With the rapid development of computer hardware and the emergence of the parallel calculation of di ﬀ raction ﬁelds, a breakthrough has been made in terms of the limitation of the unacceptable amount of computational cost to design di ﬀ ractive optical elements (DOEs), and more accurate global search algorithms can be introduced to the design of complex DOEs and holographic projections instead of traditional iterative algorithms. In this paper, a hybrid algorithm which combines particle swarm optimization (PSO) with a simulated annealing (SA) algorithm is proposed for the designing of DOEs and projecting holographic images with less noise. PSO is used to reduce the invalid disturbance in SA, and SA can jump out from local extreme points to ﬁnd the global extreme points. Compared with the traditional Gerchberg–Saxton (GS) algorithm, the simulation and experimental results demonstrate that the proposed SA–PSO hybrid algorithm can improve uniformity by more than 10%.


Introduction
Diffractive optical elements (DOEs) have already played an important role in the fields of laser collimation [1,2], optical storage [3,4], projection systems [5][6][7], infrared focal plane arrays [8,9], bioimaging [10,11], terahertz [12], and so forth.However, the diffracted light pattern of the DOEs is not accurate due to the errors in design, fabrication and other factors.Among these factors, the impact of the design is particularly important.Basically, there are two methods for designing DOEs, one is a direct calculation based on single inverse diffraction, and the other is phase retrieval.The quality of the reconstruction is determined by the phase retrieval method being better than the single inverse diffraction method.The design accuracy of traditional phase retrieval algorithms, for example, the Gerchberg-Saxton (GS) algorithm [13,14] and its modified versions [15][16][17][18], and the Yang-Gu algorithm [19,20], is relatively low, because such algorithms tend to fall into local extreme points and cannot jump out, leading to a failure in terms of good design performance.For the DOEs designed by the above algorithm, the diffraction pattern is often poor in uniformity, with more noise.Furthermore, maskless lithography [21] and 3D printing based on photopolymer resins [22,23] need higher accurate and uniform exposure patterns as much as possible, which are formed by a liquid crystal on silicon chips (LCoS; work as dynamic DOEs).Hence, in the case that traditional design methods cannot meet the accuracy and uniformity requirement, it is necessary to explore new design methods.
The search algorithms, including local search algorithms and global search algorithms, theoretically can get a better solution, compared with the traditional iterative algorithms.However, the former, Appl.Sci.2020, 10, 5485 2 of 11 including the steepest descent [24], conjugate gradient [25], and direct binary search (DBS) [26] algorithms, tend to fall into local extreme points and cannot jump out.The latter mainly includes the simulated annealing (SA) algorithm [27], genetics algorithm (GA) [28][29][30], particle swarm optimization (PSO) [31][32][33], etc.Because the SA algorithm needs to search enough time at each synthetic temperature to reach the global optimal solution, the entire annealing process is extremely time consuming.The GA has fast and random search capabilities and potential parallel computing possibilities, however, the convergent effect is highly dependent on the initial parameters set by experience and the choice of initial values.Although the PSO has certain intelligence that can make it perform exceptionally well when solving continuous function optimization problems, when it is used to design DOEs, PSO is proven to be premature, i.e., it has a poor local search ability.Generally speaking, in order to improve the quality of holographic projection, a hybrid search algorithm is needed to devise the design of DOEs with better performances.
In this paper, we propose a hybrid algorithm combining PSO with SA, named the SA-PSO hybrid algorithm, which overcomes the time-consuming and low-efficiency random disturbances, and avoids falling into local optimal solutions.The simulation result shows that the performance has been improved compared with the traditional GS algorithm.
The rest of the paper is organized as follows.In Section 2, the SA-PSO hybrid algorithm for designing DOEs is proposed.We use a flow chart to illustrate how the algorithm works.In Section 3, we use the SA-PSO hybrid algorithm to shape a plane wave into a Gaussian beam, and the results show a better uniformity than the GS algorithm.Section 4 presents the experiments to verify the good performance of the SA-PSO algorithm in optimizing the uniformity of the diffraction pattern.Finally, the conclusions are presented in Section 5.

Algorithm
The cost function of the SA-PSO hybrid algorithm is the distance between the target light intensity and the intensity diffracted by the calculated DOE, as defined by Equation (1).
where I m,n is the intensity diffracted by the calculated DOE, I m,n,ex is the expected target light intensity, and M and N represent the number of rows and columns of the matrix, respectively.During the annealing process, the initial solution and the new solution after random perturbation are compared by using an acceptance probability term.According to the Metropolis principle [34], the acceptance probability P t is related to the current synthetic temperature t.
where f i is the cost function value of the solution before the disturbance, and f j is the cost function value of the new solution.It should be noted that, the "solution" represents the phase distribution of the DOE.The light intensity distribution on the target plane is calculated from a forward diffraction from the phase distribution of the DOE, and then f can be obtained by using Equation (1).The diffraction model can be the angular spectrum, the Fresnel or the Fraunhofer, depending on the distance from the DOE to the target plane [35].By using the Metropolis principle, it is more likely to accept worse solutions at higher synthetic temperatures, but almost only better solutions can be accepted at lower synthetic temperatures.We set a disturbance coefficient la to control the disturbance's severity.
where α is the linear attenuation coefficient and β represents the smallest disturbance value.Through adjusting the values of α and β, intense disturbances at higher synthetic temperatures and slow disturbances at lower synthetic temperatures can be obtained.In this way, the initial particle swarm at synthetic temperature t is obtained, and the size of particle swarm is 30.The choice of α and β depends on the initial synthetic temperature and the ending synthetic temperature.If the linear change strategy in Equation ( 3) is used, the values of α and β can be calculated by the method of undetermined coefficients.
In the PSO process, the cost function in Equation ( 1) is still used to evaluate the initialized population, and then the population is updated according to the following rules: where ω is the variable speed update rate, which denotes that the new generation of particle swarms is affected by the overall speed of the previous generation.The value of ω is a constant between 0.8 and 1.2 based on experience.P id is the optimal solution of single particles and P gd is the optimal solution of particle swarms.C 1 and C 2 are constants between 1 and 2, which respectively represent the influence weight of P id and P gd on V id (n + 1), the larger value, with the greater influence on V id (n + 1).φ 1 and φ 2 are random numbers between 0 and 1, randomly changed in every updating in PSO.X id (n) represents the solution, i.e., the phase distribution of the DOE obtained after the n-th update.V id (n + 1) represents the amount of change when X id (n) is updated.
The disturbance coefficient la is also introduced to ensure that the first several evolutions can shuttle quickly in the solution space, while later evolutions can search finely in the solution space at a slower speed.It should be pointed out that the random perturbation strategy is where the number of random disturbance points decreases linearly with a decreasing synthetic temperature.This strategy is similar to the annealing process, which is fast at the beginning of the search, but becomes stable at the end of the search.
In theory, C1 and C2 will affect the local search ability and global search ability of the particle swarm algorithm; the higher the number of particles, the larger the search range, and the more likely it is to get the global optimal solution.However, the results of searching algorithms are uncertain, and it is impossible to obtain the best parameter value strictly by a mathematical formula.Therefore, we have explored the law through many calculations, and given the parameter values that meet our requirements based on experience.
We use a flow chart to illustrate how the algorithm works, as shown in Figure 1.The algorithm stops running when one of the following two conditions happens.The first is that the calculation value reaches the threshold close to the target value, and the second is that the number of loops reaches the maximum.

Simulation
The performance of the SA-PSO hybrid algorithm is compared with the traditional GS algorithm.For the diffraction process from DOEs to the target plane, the Fresnel diffraction equation or the Fraunhofer diffraction equation [35] can be used to calculate the intensity distribution on the target plane.The goal of the simulation is to shape a plane wave into a Gaussian beam.The diffraction model used in the process is Fresnel diffraction and the numerical algorithm used in the simulation is a single fast Fourier transform (SFFT).The intensity distribution on the target plane is shown in Figure 2a.The specific parameters used in the simulation are shown in Table 1.

Simulation
The performance of the SA-PSO hybrid algorithm is compared with the traditional GS algorithm.For the diffraction process from DOEs to the target plane, the Fresnel diffraction equation or the Fraunhofer diffraction equation [35] can be used to calculate the intensity distribution on the target plane.The goal of the simulation is to shape a plane wave into a Gaussian beam.The diffraction model used in the process is Fresnel diffraction and the numerical algorithm used in the simulation is a single fast Fourier transform (SFFT).The intensity distribution on the target plane is shown in Figure 2a.The specific parameters used in the simulation are shown in Table 1.Two DOEs are calculated by using the GS algorithm and the SA-PSO hybrid algorithm, respectively.The diffracted intensity distributions on the target plane are calculated, as shown in Figure 2b,c.The shaping effects are compared in Figure 2d.The deviation to the target by the GS algorithm and the SA-PSO hybrid algorithm are shown in Figure 2e,f.It can be seen that the SA-PSO algorithm has a better shaping effect than the GS algorithm, especially in the central high-light intensity area.The cost function values of the light intensity by the GS algorithm and the SA-PSO hybrid algorithm are 9.30 and 8.09, respectively.So, the relative reduction by the SA-PSO hybrid algorithm is 13.00%.

Experiment
In the experiment, we load the calculated 256-step phase distribution on an LCoS (liquid crystal on silicon) as a DOE.The parameters in the calculation and experiment are shown in Table 2.The target intensity pattern is a rectangle with a size of 282 × 150, as shown in Figure 3.In order to reduce the influence of the zero-order light on the experimental results, we use a lens to converge the zero-order light in the center, while placing the rectangle in the upper left corner to be far away from the zero-order light; this requires the use of the Fraunhofer diffraction equation to calculate the target plane light intensity distribution.Two 256-step phase distributions of the DOE are calculated by the GS algorithm and the SA-PSO hybrid algorithms, and then the simulated Fraunhofer diffraction results on the target plane are shown in Figure 4a,b, respectively.The calculation by the GS algorithm takes 20 s, and the one by the SA-PSO takes 1 h.The cost function values of the diffracted pattern in Figure 4a,b are 100.33 and 89.71, respectively, so the relative reduction by the SA-PSO algorithm is 10.59%.

Experiment
In the experiment, we load the calculated 256-step phase distribution on an LCoS (liquid crystal on silicon) as a DOE.The parameters in the calculation and experiment are shown in Table 2.The target intensity pattern is a rectangle with a size of 282 × 150, as shown in Figure 3.In order to reduce the influence of the zero-order light on the experimental results, we use a lens to converge the zeroorder light in the center, while placing the rectangle in the upper left corner to be far away from the zero-order light; this requires the use of the Fraunhofer diffraction equation to calculate the target plane light intensity distribution.Two 256-step phase distributions of the DOE are calculated by the GS algorithm and the SA-PSO hybrid algorithms, and then the simulated Fraunhofer diffraction results on the target plane are shown in Figure 4a,b, respectively.The calculation by the GS algorithm takes 20 s, and the one by the SA-PSO takes 1 h.The cost function values of the diffracted pattern in Figure 4a,b are 100.33 and 89.71, respectively, so the relative reduction by the SA-PSO algorithm is 10.59%.We also compared the convergence speed of the SA-PSO algorithm and SA algorithm.It can be seen from Figure 1 that, at each synthetic temperature, the particle swarm search is performed first, and then the optimal solution is used as the initial value to perform 100 simulated annealing searches, and finally the search at this synthetic temperature ends.Therefore, after removing the particle swarm search process, the simulated annealing process is obtained.Both the two algorithms start the search from the same initial value (obtained by the GS algorithm).In both We also compared the convergence speed of the SA-PSO algorithm and SA algorithm.It can be seen from Figure 1 that, at each synthetic temperature, the particle swarm search is performed first, and then the optimal solution is used as the initial value to perform 100 simulated annealing searches, and finally the search at this synthetic temperature ends.Therefore, after removing the particle swarm search process, the simulated annealing process is obtained.Both the two algorithms start the search from the same initial value (obtained by the GS algorithm).In both the SA and SA-PSO algorithms, at a synthetic temperature, we record the cost function value when the simulated annealing process begins.The entire simulated annealing process experiences 31 synthetic temperatures.Each of the two algorithms is performed five times.The changes of the cost function values during the searching processes in the SA-PSO and the SA algorithm are shown in Figure 5a,b, respectively.However, because the SA-PSO algorithm has an extra particle swarm search process than the SA algorithm, its program running time is longer than the SA algorithm.The program running times of the SA-PSO and the SA algorithms are 1 h and 15 min, respectively.
The schematic of the experimental setup is shown in Figure 6.A Holoeye Pluto 2 NIR-11 LCoS is employed to load the DOEs designed by the GS algorithm and the SA-PSO algorithm.The charge coupled device (CCD) is placed on the focal plane of the focusing lens, i.e., the Fraunhofer plane of the DOEs.The light intensity is adjusted by a continuously neutral density filter, to ensure that none of the pixels on the CCD are overexposed.The diffracted patterns are captured by the CCD and are normalized to 0~1.The pseudo color images of the normalized patterns are shown in Figure 7.Because we use the calculation result of the GS algorithm as the initial value of the SA-PSO algorithm, The initial cost function value is 99.93, and this value becomes 94.93 and 89.94 by a reduction of 5% and 10%, respectively.In order to compare the convergence speed of the two algorithms, the y = 94.93 and y = 89.94lines are used as the baseline.In each figure, we pick out the curve with the fastest convergence speed and marked the data points that just passed the baseline.It can be calculated that the calculation numbers, to achieve a reduction of 5%, reduce by 29.7% with the SA-PSO algorithm, and that value for a 10% reduction is even larger.Therefore, the convergence speed of the SA-PSO is faster, and the final cost function value of the SA-PSO is also smaller than SA.

Cost function value
However, because the SA-PSO algorithm has an extra particle swarm search process than the SA algorithm, its program running time is longer than the SA algorithm.The program running times of the SA-PSO and the SA algorithms are 1 h and 15 min, respectively.
The schematic of the experimental setup is shown in Figure 6.A Holoeye Pluto 2 NIR-11 LCoS is employed to load the DOEs designed by the GS algorithm and the SA-PSO algorithm.The charge coupled device (CCD) is placed on the focal plane of the focusing lens, i.e., the Fraunhofer plane of the DOEs.The light intensity is adjusted by a continuously neutral density filter, to ensure that none of the pixels on the CCD are overexposed.The diffracted patterns are captured by the CCD and are normalized to 0~1.The pseudo color images of the normalized patterns are shown in Figure 7.Because we use the calculation result of the GS algorithm as the initial value of the SA-PSO algorithm, the speckle distributions of the GS algorithm and the SA-PSO algorithm are quite similar.Each intensity distribution in Figure 7 is normalized.Because there is no truth value in the experiment, the standard deviation of each result is calculated to evaluate the quality.The standard deviation is expressed as: where I m,n is the intensity diffracted by the calculated DOE, I is the mean value of the target light intensity, and M and N represent the number of rows and columns of the matrix, respectively.The standard deviation of Figure 7a,b is 0.0440 and 0.0383, respectively.The relative reduction by the SA-PSO hybrid algorithm is 13.0%, which is consistent with the simulation result.
where , is the intensity diffracted by the calculated DOE, ̅ is the mean value of the target light intensity, and M and N represent the number of rows and columns of the matrix, respectively.
The standard deviation of Figure 7a,b is 0.0440 and 0.0383, respectively.The relative reduction by the SA-PSO hybrid algorithm is 13.0%, which is consistent with the simulation result.To show the feasibility of the SA-PSO algorithm on display, we use this algorithm to design a DOE for a 2D Rubik's cube, and compare the experimental reconstruction result with those by the GS algorithm and the single inverse diffraction method, as shown in Figure 8.To show the feasibility of the SA-PSO algorithm on display, we use this algorithm to design a DOE for a 2D Rubik's cube, and compare the experimental reconstruction result with those by the GS algorithm and the single inverse diffraction method, as shown in Figure 8.To show the feasibility of the SA-PSO algorithm on display, we use this algorithm to design a DOE for a 2D Rubik's cube, and compare the experimental reconstruction result with those by the GS algorithm and the single inverse diffraction method, as shown in Figure 8.Because the GS algorithm is an iteration method with many forward and backward diffractions, the reconstruction quality is better than that of the single inverse diffraction method, as shown in Figure 8c,d.The SA-PSO algorithm uses the result of the GS algorithm as the initial search value, so the imaging quality is slightly better than that of the GS algorithm, as shown in Figure 8b,c.

Conclusions
In this paper, we propose a hybrid algorithm for the design of diffractive optical elements, named the SA-PSO hybrid algorithm.PSO is used to reduce the invalid disturbance in SA, and SA can jump out from local extreme points to find the global extreme points.The simulation result of a 256 × 256 shaping demonstrates that the SA-PSO hybrid algorithm improves the accuracy of diffraction patterns.Because of the normalization operation to the captured intensity distributions, the standard deviation of each result is calculated to evaluate the quality in the experiment.The experiment results also show that the SA-PSO hybrid algorithm can improve the uniformity by more than 10%.
The performance of the SA-PSO algorithm depends on the number of searches and the size of the particle swarm.Therefore, increasing these two parameters can effectively improve the performance of the algorithm, but this will lead to an increase in computing time.The proposed algorithm with a better performance will be practical in the rapid designing of the computer hardware.Because the GS algorithm is an iteration method with many forward and backward diffractions, the reconstruction quality is better than that of the single inverse diffraction method, as shown in Figure 8c,d.The SA-PSO algorithm uses the result of the GS algorithm as the initial search value, so the imaging quality is slightly better than that of the GS algorithm, as shown in Figure 8b,c.

Conclusions
In this paper, we propose a hybrid algorithm for the design of diffractive optical elements, named the SA-PSO hybrid algorithm.PSO is used to reduce the invalid disturbance in SA, and SA can jump out from local extreme points to find the global extreme points.The simulation result of a 256 × 256 shaping demonstrates that the SA-PSO hybrid algorithm improves the accuracy of diffraction patterns.Because of the normalization operation to the captured intensity distributions, the standard deviation of each result is calculated to evaluate the quality in the experiment.The experiment results also show that the SA-PSO hybrid algorithm can improve the uniformity by more than 10%.
The performance of the SA-PSO algorithm depends on the number of searches and the size of the particle swarm.Therefore, increasing these two parameters can effectively improve the performance of the algorithm, but this will lead to an increase in computing time.The proposed algorithm with a better performance will be practical in the rapid designing of the computer hardware.

Figure 1 .
Figure 1.The flow chart of the simulated annealing-particle swarm optimization (SA-PSO) hybrid algorithm.

Figure 1 .
Figure 1.The flow chart of the simulated annealing-particle swarm optimization (SA-PSO) hybrid algorithm.

FigureFigure 2 .Figure 2 .
Figure2b,c.The shaping effects are compared in Figure2d.The deviation to the target by the GS algorithm and the SA-PSO hybrid algorithm are shown in Figure2e,f.It can be seen that the SA-PSO algorithm has a better shaping effect than the GS algorithm, especially in the central high-light intensity area.The cost function values of the light intensity by the GS algorithm and the SA-PSO hybrid algorithm are 9.30 and 8.09, respectively.So, the relative reduction by the SA-PSO hybrid algorithm is 13.00%.

Figure 4 .
Figure 4.The simulated diffraction results of the diffractive optical elements (DOEs) designed by (a) the GS algorithm and (b) the SA-PSO algorithm.

Figure 5 .
Figure 5.The convergence process of the SA-PSO hybrid algorithm (a) and the simulated annealing (SA) algorithm (b).

Figure 5 .
Figure 5.The convergence process of the SA-PSO hybrid algorithm (a) and the simulated annealing (SA) algorithm (b).

Figure 6 .
Figure 6.Schematic of the experimental setup.

Figure 7 .
Figure 7.The diffracted intensity patterns obtained by the CCD.(a) The GS algorithm and (b) the SA-PSO algorithm.

Figure 6 .
Figure 6.Schematic of the experimental setup.
Appl.Sci.2020, 7, x FOR PEER REVIEW 9 of 12where , is the intensity diffracted by the calculated DOE, ̅ is the mean value of the target light intensity, and M and N represent the number of rows and columns of the matrix, respectively.The standard deviation of Figure7a,b is 0.0440 and 0.0383, respectively.The relative reduction by the SA-PSO hybrid algorithm is 13.0%, which is consistent with the simulation result.

Figure 6 .
Figure 6.Schematic of the experimental setup.

Figure 7 .
Figure 7.The diffracted intensity patterns obtained by the CCD.(a) The GS algorithm and (b) the SA-PSO algorithm.

Figure 7 .
Figure 7.The diffracted intensity patterns obtained by the CCD.(a) The GS algorithm and (b) the SA-PSO algorithm.

Figure 8 .
Figure 8. Two-dimensional Rubik's Cube patterns obtained by the CCD, using different algorithms.(a) Target pattern, (b) the SA-PSO algorithm, (c) the GS algorithm and (d) the single inverse diffraction method.

Figure 8 .
Figure 8. Two-dimensional Rubik's Cube patterns obtained by the CCD, using different algorithms.(a) Target pattern, (b) the SA-PSO algorithm, (c) the GS algorithm and (d) the single inverse diffraction method.

Table 1 .
Parameters in the simulation EndSearch for optimal solution jk with PSO, then replace i with jk Disturb i randomly, generate particle swarm jn ?Evaluate i by using the cost function f(i) Replace i with jm as a new solution Output the result Reduce the temperature tThe initial solution i the temperature Ts

Table 1 .
Parameters in the simulation

Table 2 .
Parameters in simulation.

Table 2 .
Parameters in simulation