## Abstract

We introduce a new multi-objective genetic algorithm for wavefront shaping and realize controllable multi-point light focusing through scattering medium. Different from previous single-objective optimization genetic algorithms, our algorithm named Non-dominated Sorting Genetic Algorithm II based on hybrid optimization scheme (NSGA2-H) can make all focus points have uniform intensity while ensuring that their enhancement is as high as possible. We demonstrate the characteristics of NSGA2-H through simulations and experiments in amplitude optimization, analyze its optimization mechanisms and show its powerful optical control capability in uniform intensity focusing and even in customizable intensity focusing. This research will be expected to further promote future practical applications based on multi-point focusing of wavefront shaping, especially in optical trapping and optogenetics.

© 2019 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Wavefront shaping is a control technology of optical wave field. It can realize the refocusing of light through the scattering medium by modulating the amplitude or phase of the light wavefront. Since wavefront shaping was first demonstrated in experiment by Vellekoop and Mosk in 2007 [1], more and more people have participated in the research of this technology, which has made it develop rapidly and bring plenty of new applications, such as optical trapping [2], super-resolution imaging [3], endoscopy [4], cryptography [5], phototherapy [6], optogenetics [7] and optical communication [8].

Feedback-based optimization [9–21], transmission matrix measurement [22–27] and optical phase conjugation [28–33] are three classes of important methods in wavefront shaping. Over the years, many efforts have been made on single-point focusing with these methods. In contrast, there are fewer studies on multi-point focusing. In fact, multi-point focusing is more significant than single-point focusing in some specific research fields, such as optical trapping and optogenetics [2,7]. Unlike single-point focusing, not only the intensity (or enhancement) but also the uniformity (or relative standard deviation) of focus points needs to be taken into account on multi-point focusing. Therefore, multi-point focusing is essentially a more complex optimization process and it is susceptible to noise. Feedback-based optimization, especially genetic algorithm (GA), has stronger anti-noise ability than the other two methods, so it is more suitable for multi-point focusing.

In the past few years, people have tried to achieve multi-point focusing based on the GA [11] by changing discriminant. In 2012, Conkey *et al*. proposed a multi-point focusing discriminant as the fitness function of GA [34]. They linked the total intensity of focus points and the standard deviation of the intensities of focus spots through a fixed weight coefficient *M* (equal to the number of focus points) in the discriminant, and realized the reproduction of multi-point focusing patterns with uniform intensity. However, the discriminant with fixed weight coefficient *M* actually cannot ensure that the enhancement is as high as possible even if the intensity of focus points is uniform. Perhaps a suitable value can be used to replace the fixed weight coefficient to balance the intensity and uniformity of focus points for a better result, but it is not easy to find such a value in different situations by this combination method. After that, the Pearson Correlation Coefficient was proposed as the fitness function of GA and the optimized reproduction of the pre-designed patterns is also achieved [35–37]. However, this discriminant is easily influenced by the size of the background area. When all focus points are distributed in a large background area, the optimization effect will become very poor. Earlier this year, another discriminant was proposed to complete 3-point focusing [38]. But this discriminant will not make the intensity of focus points uniform if the intensity of focus points does not exceed the dynamic range of CCD.

In fact, multi-point focusing is a multi-objective optimization problem. The multi-objective optimization problem does not have only one optimal solution, but a global optimal solution set—Pareto-optimal solutions [39]. In the problem of multi-point focusing, we prefer to get the solutions corresponding to focus points with high intensity and good uniformity, i.e., the enhancement is as high as possible and the relative standard deviation is as low as possible. But these GAs mentioned above still belong to single-objective optimization algorithm, which cannot solve the multi-objective problem perfectly in a sense.

As one of the best multi-objective optimization algorithms, Non-dominated Sorting Genetic Algorithm II (NSGA2) has been applied in many fields [40–43]. In this paper, for the first time, we propose a new NSGA2 based on hybrid optimization scheme (NSGA2-H) for multi-point focusing of wavefront shaping. We demonstrate some characteristics of NSGA2-H through simulations and experiments in amplitude optimization, analyze its optimization mechanisms and show its outstanding performance in uniform intensity and customizable intensity focusing.

## 2. Principle and method

#### 2.1 The model of multi-point focusing of wavefront shaping

Principle of multi-point focusing of wavefront shaping is shown in Fig. 1. It is a closed-loop feedback system, which consists of four parts: input light, scattering medium, output light and iterative algorithm. Input light can be modulated by a digital micromirror device (DMD) (blue module in Fig. 1), and output light can be collected with a CCD camera (black module in Fig. 1). Iterative algorithm completed on PC is used to optimize input light, according to output light, for realizing multi-point focusing through scattering medium. Similar to the principle of single-point focusing based on the interference superposition of light, the whole optimization process of multi-point focusing is to use iterative algorithm to turn on the segments that are beneficial to the optimization objectives and turn off the detrimental segments for achieving the optimal multi-point focusing effect.

Our numerical simulations of wavefront shaping are based on the transmission matrix model, which is defined as follows [9]: ${E^{out}} = T \cdot {E^{in}}$. The input light field ${E^{in}}$ can be divided into *N* segments and it is a $N \times 1$ real matrix containing only 0 and 1 in amplitude optimization. The scattering medium is characterized by transmission matrix *T*, which is a $N \times N$ complex matrix and submits to the circular Gaussian random distribution. The output light field ${E^{out}}$ is a $N \times 1$ complex matrix, and the intensity of the m-th output channel can be expressed as ${I_m} = {|{E_m^{out}} |^2}$.

As mentioned above, there are two objectives for multi-point focusing optimization: the intensity and the uniformity of focus points. In this paper, we design two discriminants different from those in Ref. [34], such as Eq. (1) and Eq. (2):

*M*the number of focus points, ${I_{ref}}$ the average intensity of the output light field and ${\sigma _M}$ the standard deviation of the intensity of

*M*focus points. ${f_1}$ is equal to the total intensity of

*M*focus points divided by the average intensity of the output light field, and it is used to characterize the total enhancement of multi-point focusing optimization. Usually, the higher the enhancement, the higher the signal-to-noise ratio and the intensity of focus points. So we use the enhancement instead of the intensity to evaluate the first objective. The second objective ${f_2}$ is the relative standard deviation, also known as the coefficient of variation. It is obtained by dividing the standard deviation of the intensity of focus points by the average intensity of focus points, and it is used to characterize the discreteness of the intensity values of focus points. The smaller the coefficient of variation, the better the uniformity of focus points.

#### 2.2 NSGA2-H for multi-point focusing of wavefront shaping

GA is a heuristic global random search algorithm. It can search for the global optimal solution by eliminating worse solutions and retaining better solutions continuously based on an evolutionary strategy, which includes some key operations such as evaluation, selection, crossover, mutation and elitist preservation strategy. Among them, evaluation operation is to use the fitness function of optimization problem to obtain the fitness value of each solution. Then, an appropriate rank can be assigned to each solution according to their fitness values by sorting, and selection operation refers to the rank result to choose good solutions. Crossover and mutation operations use the result of selection to increase the diversity of solutions and avoid a local optimal solution. Finally, elitist preservation strategy can prevent the loss of good solutions and ensure the convergence of the whole optimization process. In wavefront shaping, solutions mentioned above can be expressed as input masks and obtaining the optimal mask (global optimal solution) is the ultimate goal of optimization.

There are two main differences between NSGA2-H and GA. One is the introduction of non-dominated sorting and crowding distance in the selection process. Unlike the single-objective sorting method in GA, non-dominated sorting in NSGA2-H is a special sorting method that can sort different solutions and assign each of them a reasonable rank according to two or more optimization objectives at the same time. For those solutions with the same rank, crowding distance is used as a new sorting indicator to ensure the uniform distribution of solutions. More detailed rules of these two operations can be found in Ref. [44]. The main idea of non-dominated sorting in multi-point focusing of wavefront shaping is that the larger the enhancement or the smaller the coefficient of variation, the lower the rank of corresponding mask. Then the mask with lower rank will be more easily selected for subsequent optimization, so that the whole algorithm will move towards the convergence direction of larger enhancement and smaller coefficient of variation.

The other one is the introduction of mixing factor *H* in the evaluation process for improving the enhancement of the algorithm as much as possible. In fact, mixing factor *H*, which divides the whole optimization into two stages, can be defined as the total algebra of the first stage optimization. In the first stage (${F_1} ={-} {f_1}$, ${F_2} = 100$, ${F_1}$ and ${F_2}$ are the fitness functions), which is actually a single-objective optimization process, we only consider enhancement and ignore the influence of coefficient of variation; in the second stage (${F_1} ={-} {f_1}$, ${F_2} = {f_2}$), as mentioned earlier, we consider the effects of enhancement and coefficient of variation comprehensively until the whole optimization is completed. Here the minus sign in front of ${f_1}$ is to turn it into a minimum problem for convenience. ${F_2}$ in the first stage is set to a fixed value of 100, which is large enough that the value of ${F_2}$ in the second stage is always smaller than it. This can ensure that the second stage optimization is effective. The concrete flowchart of NSGA2-H is depicted in Fig. 2. The detailed procedures of NSGA2-H in simulation are as follows:

**Step 1:**Initialize parameters and evaluate fitness values.- Step 1.1: Set $k = 0$. Assign population
*NP*, generation*NG*, segment*N*, crossover probability*P*, mutation probability_{c}*P*, mixing factor_{m}*H*, the number of focus points*M*and record the coordinate position of each focus point; generate a transmission matrix*T*and*NP*initial random masks. - Step 1.2: Calculate the corresponding output light field of each initial random mask by the transmission matrix model. Evaluate the fitness values of masks according to the fitness functions (${F_1} ={-} {f_1}$, ${F_2} = 100$).

**Step 2:**Identify the non-dominated rank and calculate crowding distance for each mask according to their fitness values.**Step 3:**Select*NP*parent masks by binary tournament selection.- Step 3.1: Randomly select two masks from
*NP*masks (which come from step 1.1 when $k = 0$ or step 7 when $k > 0$) each time. - Step 3.2: Compare their ranks, and the one with lower rank wins. If the ranks are the same, the one with larger crowding distance wins.
- Step 3.3: Repeat steps 3.1-3.2
*NP*times and get*NP*parent masks.

**Step 4:**Generate*NP*offspring masks by uniform crossover and single point mutation.- Step 4.1: Randomly divide
*NP*parent masks into*NP*/2 groups, and two masks in each group are ${P_1}$ and ${P_2}$. - Step 4.2: Generate a random binary mask
*S*and calculate the new masks by uniform crossover: ${C_1} = S \cdot {P_1} + (1 - S) \cdot {P_2}$, ${C_2} = S \cdot {P_2} + (1 - S) \cdot {P_1}$. - Step 4.3: Obtain two offspring masks by single point mutation of ${C_1}$ and ${C_2}$.
- Step 4.4: Repeat steps 4.2-4.3
*NP*/2 times and get*NP*offspring masks.

**Step 5:**Evaluate the fitness values of masks for*NP*offspring masks (if $k < H$, ${F_1} ={-} {f_1}$, ${F_2} = 100$; else, ${F_1} ={-} {f_1}$, ${F_2} = {f_2}$).**Step 6:**Mix*NP*masks with*NP*offspring masks, and execute step 2 for 2·*NP*mixed masks.**Step 7:**Execute elitist preservation strategy: Retain*NP*new offspring masks with lower rank and larger crowding distance (rank has a high priority) from the 2·*NP*mixed masks. Set $k = k + 1$.**Step 8:**If $k ={=} NG$, go to step 9; otherwise, return to step 3.**Step 9:**Plot the Pareto-optimal front figure according to the fitness values of*NP*new offspring masks (X-axis ${f_1} ={-} {F_1}$, Y-axis ${f_2} = {F_2}$). Output*NP*new offspring masks, which are Pareto-optimal solutions.

## 3. Experimental setup

Our experimental setup for multi-point focusing is shown in Fig. 3. A laser beam (532 nm, 30 mw) is expanded by an objective lens OL1 (25×, NA = 0.4), and the middle part of the input light passes through an aperture A (20 mm) and a linear polarizer P for covering the entire mask area on the DMD (TI, DLP6500FYE). The reflected light from DMD is narrowed by lenses L1 (f = 175 mm) and L2 (f = 50.8 mm), and a spatial filter SF is used to eliminate the redundant diffraction orders. After passing through OL2 (Plan 10×, NA = 0.25), a scattering sample S (Edmund ground glass diffuser, 120 grit, 1.6 mm thickness) and OL3 (Plan 20×, NA = 0.4), the output light is finally collected by a 12-bit CCD camera (AVT, Manta G-031B) with 656×492 pixels. The focal plane of the OL3 is located on the rear surface of the sample, which is the actual position of multi-point focusing. In particular, all components in the blue dotted line need to be covered with a black box to avoid interference from ambient light. The DAQ (NI USB-6009) controlled by PC is used to ensure the synchronization of timing sequence between DMD and CCD. In the software aspects, except that the CCD is controlled by LabVIEW, the rest are run on MATLAB, including the control of trigger signal of DAQ, the display of mask on DMD and the implementation of NSGA2-H. The actual optimization speed of our system can reach 11 frames per second, i.e. 11 measurements per second. The main limitations are the data transmission rate of mask and the acquisition frame rate of CCD.

## 4. Results and analysis

#### 4.1 Effects of mixing factor H on multi-point focusing

Mixing factor *H* is an important parameter for NSGA2-H. Here we will illustrate the effects of *H* on multi-point focusing through simulations and experiments of five-point focusing in amplitude optimization. The other parameters of NSGA2-H: population *NP* = 20, generation *NG* = 400, segment *N* = 3600, crossover probability *P _{c}* = 0.9, mutation probability

*P*= 1/

_{m}*N*.

In Figs. 4(a)–4(c), we compare three different values of *H* (0.3NG, 0.4NG, 0.5NG) for revealing the role of *H* in NSGA2-H. As shown in Fig. 4(a), the X-axis is the total enhancement of five focus points, and the Y-axis is their coefficient of variation. We can find that the Pareto-optimal front of *H* = 160 is obviously better than that of the other two values both in simulation and experiment. Figures 4(b) and 4(c) clearly show the evolution of enhancement and coefficient of variation. For *H* = 200, the minimum coefficient of variation exceeds 20%, which leads to a poor uniformity of focus points in Figs. 4(f) and 4(i). Although the minimum coefficient of variation of *H* = 120 is less than 1%, its maximum enhancement is the smallest of the three values. So that the results of simulation and experiment show more serious background noise in Figs. 4(d) and 4(g). In contrast, the focusing effects of *H* = 160 in Figs. 4(e) and 4(h) are more satisfactory. In other words, an appropriate value of *H*, which makes an effective trade-off between enhancement and coefficient of variation, can ensure a higher enhancement with acceptable uniformity. Because of the disturbance of laser and the noise of CCD, the enhancement of the experiment is lower than that of the simulation. Fortunately, the effects of *H* show a consistent trend in experiment and simulation, and we can easily choose the appropriate *H* value through several attempts. In Appendix A, we make a more detailed discussion on the mixed factor *H* based on the full text, including the limit values and the appropriate value of *H*.

#### 4.2 Comparison of arithmetic mean (AM) and geometric mean (GM) discriminants

As mentioned above, an appropriate *H* can improve the enhancement of the algorithm to a certain extent. In this section, we will replace the original discriminant (Eq. (1)) with a new discriminant (Eq. (3)) to further improve the enhancement of NSGA2-H. This modification scheme is inspired by arithmetic-geometric mean inequality. In fact, the original discriminant can be called as arithmetic mean (AM) discriminant, and the new one is geometric mean (GM) discriminant.

*NG*= 600. In Figs. 5(a)–5(c),

*H*= 180 and 500 are the appropriate values for NSGA2-H with AM and GM discriminants, respectively. We can find that NSGA2-H with GM discriminant owns higher enhancement than with AM discriminant in 600 generations from Fig. 5(b). In Fig. 5(c), the final coefficients of variation of both methods are reduced to 1% ∼ 2%, but NSGA2-H with GM discriminant has a much lower coefficient of variation (about 10%) than with AM discriminant (more than 60%) at the end of the first stage of optimization. This means GM discriminant can not only improve the enhancement but also reduce the coefficient of variation as much as possible. Thus the value of

*H*for NSGA2-H with GM discriminant can be much larger than with AM discriminant so that it can obtain higher enhancement under similar coefficient of variation. Figure 5(a) shows that NSGA2-H with GM discriminant has better Pareto-optimal fronts. The focusing effects of simulation (Figs. 5(d) and 5(e)) and experiment (Figs. 5(f) and 5(g)) also prove this point. In Appendix B, we analyze the optimization mechanisms of AM and GM discriminants in NSGA2-H, and further explain the advantage of GM discriminant in multi-point focusing.

#### 4.3 Uniform intensity focusing for more points and customizable intensity focusing

In this section, we will demonstrate the capability of NSGA2-H (with GM discriminant) in uniform intensity focusing for more points and customizable intensity focusing through experiments. Firstly, we tried the uniform intensity focusing of 20-point. Figure 6(a) shows the image captured by the CCD before optimization. Here our optimization target is a Golden Spiral with 20 points. After 32000 measurements (*NP* = 80, *NG* = 400, *H* = 350), it took 48.5 minutes and we got a good focusing effect in Fig. 6(b). Its total enhancement is 509, and the coefficient of variation is about 5%. Secondly, we tested customizable intensity focusing of 9-point. In the process of optimization, the pre-designed intensity ratio of focus points, which is 1:2:4, should be added to the calculation of fitness values. After 20000 measurements (*NP* = 20, *NG* = 1000, *H* = 400), it took 30.3 minutes and we also got the expected result in Fig. 6(c). The total enhancement is 395, and the intensity of these focus points is shown in ascending order in Fig. 6(d).

In practice, focusing of more points and more complex customizable intensity optimization can also be achieved if we continue to increase the number of measurements in NSGA2-H. However, the overall focusing effect will be weakened because the average enhancement (${{{f_{AVG}} = {f_1}} \mathord{\left/ {\vphantom {{{f_{AVG}} = {f_1}} M}} \right.} M}$) decreases with the increase of the number of focus points *M*. When the average enhancement is too small, the intensity of focus point is easily submerged in noise, which may lead to the failure of optimization. An effective method is to increase the average enhancement ${f_{AVG}}$ by increasing the total enhancement ${f_1}$. The theoretical formula of total enhancement ${f_1}$ is defined as follows [45]:

*N*is the segment, $\alpha \in [0,1]$ depends on the type of light modulation and the value of $\alpha$ in phase optimization is 2-5 times that in amplitude optimization. Therefore, we can increase segment

*N*or use phase optimization instead of amplitude optimization to achieve a better optimization effect for more focus points in future work.

## 5. Conclusion and outlook

In conclusion, we have proposed a new multi-objective genetic algorithm NSGA2-H for multi-point light focusing in wavefront shaping. Through the analysis of simulation and experimental results in amplitude optimization, we have clarified that NSGA2-H with the appropriate value of mixing factor *H* can ensure a higher enhancement with acceptable uniformity. Additionally, the mechanism of *H* in NSGA2-H was illustrated and the empirical range of *H* was given in Appendix A, which would help us select the appropriate value of *H* quickly to make better use of the algorithm for multi-point focusing. Then, we used GM discriminant instead of AM discriminant, and the comparison results confirmed that GM discriminant could further improve the enhancement of NSGA2-H. The corresponding optimization mechanisms of these two discriminants also were explained in Appendix B. The final experiments on 20-point uniform intensity focusing and 9-point customizable intensity focusing have fully demonstrated the powerful multi-point focusing capability of NSGA2-H. Additionally, NSGA2-H also has great potential for focusing of more points and more complex customizable intensity optimization.

For more complex conditions, multi-point focusing optimization of wavefront shaping based on NSGA2-H is attractive and promising. For instance, it may be used to achieve multi-point focusing with specific intensity distribution in different positions in three-dimensional space for optical trapping and optogenetics, or realize color focusing with multi-wavelength light for fluorescence imaging. These are potential applications of multi-point focusing, and we may only need to modify the objective functions of NSGA2-H slightly to meet different requirements. But there are still many limitations in these applications at present, especially how to further improve the optimization speed of the whole system and how to effectively collect the feedback signals of focus points in different spatial locations or different wavelengths at the same time. We believe that these limitations will be overcome with the deepening of research, and the improvement of optical control capability in multi-point focusing will certainly bring more valuable applications for wavefront shaping.

## Appendix A: Discussion on the limit values and the appropriate value of mixing factor *H*

In Section 4.1 of the text, we can find that mixing factor *H* plays an important role in NSGA2-H. In fact, someone may wonder why the mixing factor *H* has such a great impact on the algorithm, and what the actual mechanism of *H* is in the algorithm. In order to better explain this point, we will discuss the two limit values of *H* and compare them with the appropriate value of *H* though the simulations of five-point focusing in amplitude optimization.

Since *H* represents the total algebra of the first stage optimization, its value ranges from 0 to *NG*. Here, *NG* is the total algebra of the whole algorithm. It is obvious that 0 and *NG* are two limit values of *H*. When *H = *0, NSGA2-H which contains only the second stage optimization (${F_1} ={-} {f_1}$, ${F_2} = {f_2}$) can be transformed into NSGA2. When *H* = *NG*, NSGA2-H only carries out the first stage optimization (${F_1} ={-} {f_1}$, ${F_2} = 100$) and it can be regarded as GA, that is, the traditional single-objective optimization genetic algorithm. Figures 7(a) and 7(b) show evolution curves of enhancement and coefficient of variation of three algorithms based on AM discriminant. An appropriate value of *H* is 400 for NSGA2-H. The other relevant parameters: population *NP* = 20, generation *NG* = 1000, segment *N* = 3600, crossover probability *P _{c}* = 0.9, mutation probability

*P*= 1/

_{m}*N*.

In Fig. 7(a), we can find that the evolution speed of GA is obviously faster than that of NSGA2. For example, the enhancement of GA in the 200th generation is the same as that of NSGA2 in the 1000th generation. The huge difference in optimization speed can be attributed to different optimization mechanisms of the two algorithms, especially the unique sorting mechanism of NSGA2. Under the non-dominated sorting mechanism of NSGA2, the masks with high enhancement or low coefficient of variation have high probability to be selected. However, in the initial stage of optimization, a few of masks with low coefficient of variation and even low enhancement may also be selected to participate in the subsequent optimization. This weakens the enhancement of offspring masks to a certain extent, so the evolution speed of NSGA2 is much slower than that of GA and the enhancement of NSGA2 is less than half of that of GA at the end of optimization. Nevertheless, the coefficient of variation of NSGA2 can reach about 0.1% when the optimization is completed in Fig. 7(b). On the contrary, the actual coefficient of variation of GA in the whole optimization process is not less than 40% (the coefficient of variation is set to a fixed value of 100%, so the actual value is not shown in Fig. 7(b)). To make full use of the advantages of GA and NSGA2, we introduce the mixed factor *H* and mix these two algorithms into NSGA2-H. *H* represents the total algebra of GA (the first stage optimization). The first stage optimization takes advantage of the speed advantage of GA to quickly increase the enhancement to a certain value. Under the condition of high enhancement, NSGA2 can be used to make the intensity of focus points more uniform in the second stage. In Figs. 7(a) and 7(b), we can see that this successful hybrid optimization scheme make NSGA2-H perform well. When the optimization is completed, its enhancement is more than 500 and its coefficient of variation also converges to about 0.1%.

In addition, similar numerical simulations also have been done for the GM discriminant proposed in Section 4.2 of the text. The appropriate value of *H* for NSGA2-H is 900. As shown in Figs. 8(a) and 8(b), the comprehensive performance of NSGA2-H is also better than that of NSGA2 and GA. But it is worth noting that even GA, its final coefficient of variation can also reach a small value about 2% (the actual value is not shown in Fig. 8(b)) due to the characteristics of GM discriminant. Therefore, it takes only a few generations for NSGA2-H to make the coefficient of variation converge from 2% to 0.1% in the second stage optimization. So the value of *H* and enhancement in NSGA2-H with GM discriminant can be larger than that with AM discriminant.

As mentioned in the text, the appropriate value of *H* is an effective trade-off between enhancement and coefficient of variation. The final enhancement can increase with the increase of the value of *H*, but an excessive *H* may lead to inadequate generation in the second stage optimization, resulting in uneven intensity of focus points. Thus, how to determine the appropriate value of *H* is a key issue. Although it is difficult to find the appropriate value of *H* from a theoretical perspective because of the complex optimization mechanism of NSGA2-H, we can analyze the law of appropriate value of *H* for different numbers of focus points in simulation and obtain the empirical range of appropriate value of *H* as a reference for practical application. For this purpose, we tested six groups of simulations based on AM and GM discriminants with different numbers of focus points. Each group was repeated five times and its appropriate value of *H* was determined in advance by several attempts. The final results are shown in Table 1. To ensure the fairness of comparison, the other relevant parameters of these simulations are the same as those above.

From Table 1, we can find that whether AM discriminant or GM discriminant, appropriate values of *H* can ensure that their coefficients of variation are close and acceptable at the same number of focus points. As the number of focus points increases from 5 to 20, the appropriate value of *H* based on AM discriminant decreases from 400 to 300, while the appropriate value of *H* based on GM discriminant decreases from 900 to 800. This is because when the number of focus points increases, the degree of interaction between different focus points will become more intense, and the complexity of optimization will increase dramatically. At this time, more generations need to be provided for the second stage optimization to ensure the uniformity of focus points. Thus, the appropriate value of *H* will decrease when *NG* does not change. Here, the empirical range of *H* is 0.3 ∼ 0.4*NG* (AM) or 0.8 ∼ 0.9*NG* (GM) for 5 ∼ 20 focus points. As a reference, it can help us lock in the appropriate value of *H* more quickly in experiment. For more focus points and even complex customizable intensity optimization mentioned in Section 4.3 of the text, *H* can be a lower value than the empirical range to ensure the optimization effect.

## Appendix B: Optimization mechanisms of arithmetic mean (AM) and geometric mean (GM) discriminants in NSGA2-H

In Section 4.2 of the text, we can find that GM discriminant is better than AM discriminant in multi-point focusing. A key issue is why the intensity of focus points based on GM discriminant is much more uniform than that based on AM at the end of the first stage of optimization. In this appendix, we will explain it by clarifying the relationship between AM and GM discriminants and analyzing their optimization mechanisms in NSGA2-H.

AM and GM discriminants can be expressed as ${f_{AM}} = {{\sum\limits_{m = 1}^M {{I_m}} } \mathord{\left/ {\vphantom {{\sum\limits_{m = 1}^M {{I_m}} } {{I_{ref}}}}} \right. } {{I_{ref}}}}$ and ${f_{GM}} = {{\sqrt[M]{{\prod\limits_{m = 1}^M {{I_m}} }} \bullet M} \mathord{\left/ {\vphantom {{\sqrt[M]{{\prod\limits_{m = 1}^M {{I_m}} }} \bullet M} {{I_{ref}}}}} \right. } {{I_{ref}}}}$, respectively. They are subject to the following relationship:

Equation (B3) is actually an arithmetic-geometric mean inequality. If and only if ${I_1} = {I_2} = \cdots = {I_m}$, the equal sign holds.

In the process of multi-point focusing optimization of NSGA2-H, the main difference between ${f_{AM}}$ and ${f_{GM}}$ is reflected in the first stage of optimization, in which only enhancement is taken as the optimization objective. According to the selection mechanism of genetic algorithm, a better solution with higher enhancement will have a higher probability to be selected, so that the searched solutions will approach the optimal solution after a certain time of iteration optimization. Assuming that all feasible solutions have the chance to be retrieved randomly in the first stage of optimization. Set $\sum\limits_{m = 1}^M {{I_m}} = A$ (A is an arbitrary constant), then ${f_{AM}} = {A \mathord{\left/ {\vphantom {A {{I_{ref}}}}} \right.} {{I_{ref}}}}$. According to the properties of the arithmetic-geometric mean inequality in Eq. (B1), the maximum value of ${f_{GM}}$, which is ${A \mathord{\left/ {\vphantom {A {{I_{ref}}}}} \right.} {{I_{ref}}}}$, can be obtained if and only if ${I_1} = {I_2} = \cdots = {I_m}$. In other words, the smaller the difference of ${I_1} \sim {I_m}$, the larger the corresponding enhancement. Therefore, for GM discriminant, these solutions with small differences between ${I_1} \sim {I_m}$ will have a high probability to be selected. However, all different feasible solutions with random differences of ${I_1} \sim {I_m}$ for AM discriminant have equal chance to be selected because of the same enhancement ${A \mathord{\left/ {\vphantom {A {{I_{ref}}}}} \right.} {{I_{ref}}}}$. Next, we will see the difference between AM and GM discriminants more intuitively through a simple example. Set $M = 2$, $\sum\limits_{m = 1}^M {{I_m}} = 10$ (${I_1},{I_2} \in [{0,10} ]$), ${I_{ref}} = 1$, then ${f_{AM}} = {I_1} + {I_2}$, ${f_{GM}} = \sqrt {{I_1}\cdot {I_2}} \times 2$. Figure 9 shows the relationship between ${I_1},{I_2}$ and ${f_{}}$ with these two discriminants.

From Fig. 9, we can clearly see that when the sum of the intensity of two focus points is a fixed value, the curves of ${f_{AM}}$ and ${f_{GM}}$ are a straight line and a parabola, respectively. To make the enhancement *f* as large as possible, there are countless combinations of feasible solutions for ${f_{AM}}$ (e.g., ${I_1} = 0.5,{I_2} = 9.5$ or ${I_1} = 3,{I_2} = 7$), but only one perfect feasible solution for ${f_{GM}}$ (${I_1} = {I_2} = 5$). For AM discriminant, the value of *f* is fixed so that all feasible solutions will be selected randomly with same probability, which makes the uniformity of focus points poor. However, for GM discriminant, those feasible solutions, which have higher *f* with smaller difference between ${I_1}$ and ${I_2}$, will be selected with high probability, so the uniformity of focus points will be much better. Of course, in the actual optimization process, it is generally difficult for ${f_{GM}}$ to converge to a perfect feasible solution. Not all feasible solutions have the chance to be retrieved under the limited number of measurements on the one hand, on the other hand, the solutions distributed near the perfect feasible solution also have a high probability to be selected, which results in that the optimal solutions obtained in more cases are only distributed near the perfect feasible solution. Nevertheless, the intensity of focus points based on GM discriminant is much more uniform than that based on AM discriminant.

For more focus points, we cannot explain it graphically (4-D or more), but the optimization mechanisms of AM and GM discriminants in NSGA2-H are the same as described above.

## Funding

Key Technologies Research and Development Program (2017YFB1104500); National Natural Science Foundation of China (21627813); Natural Science Foundation of Beijing Municipality (7182091); China-Japan Friendship Hospital (PYBZ1801).

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **I. M. Vellekoop and A. P. Mosk, “Focusing coherent light through opaque strongly scattering media,” Opt. Lett. **32**(16), 2309 (2007). [CrossRef]

**2. **T. Čižmár, M. Mazilu, and K. Dholakia, “In situ wavefront correction and its application to micromanipulation,” Nat. Photonics **4**(6), 388–394 (2010). [CrossRef]

**3. **E. G. van Putten, D. Akbulut, J. Bertolotti, W. L. Vos, A. Lagendijk, and A. P. Mosk, “Scattering Lens Resolves Sub-100 nm Structures with Visible Light,” Phys. Rev. Lett. **106**(19), 193905 (2011). [CrossRef]

**4. **Y. Choi, C. Yoon, M. Kim, T. D. Yang, C. Fang-Yen, R. R. Dasari, K. J. Lee, and W. Choi, “Scanner-Free and Wide-Field Endoscopic Imaging by Using a Single Multimode Optical Fiber,” Phys. Rev. Lett. **109**(20), 203901 (2012). [CrossRef]

**5. **S. A. Goorden, M. Horstmann, A. P. Mosk, B. Škorić, and P. W. H. Pinkse, “Quantum-secure authentication of a physical unclonable key,” Optica **1**(6), 421 (2014). [CrossRef]

**6. **F. Wang, H. He, H. Zhuang, X. Xie, Z. Yang, Z. Cai, H. Gu, and J. Zhou, “Controlled light field concentration through turbid biological membrane for phototherapy,” Biomed. Opt. Express **6**(6), 2237 (2015). [CrossRef]

**7. **H. Ruan, J. Brake, J. E. Robinson, Y. Liu, M. Jang, C. Xiao, C. Zhou, V. Gradinaru, and C. Yang, “Deep tissue optical focusing and optogenetic modulation with time-reversed ultrasonically encoded light,” Sci. Adv. **3**(12), eaao5520 (2017). [CrossRef]

**8. **Z. Cao, X. Zhang, G. Osnabrugge, J. Li, I. M. Vellekoop, and A. M. J. Koonen, “Reconfigurable beam system for non-line-of-sight free-space optical communication,” Light: Sci. Appl. **8**(1), 69 (2019). [CrossRef]

**9. **I. M. Vellekoop and A. P. Mosk, “Phase control algorithms for focusing light through turbid media,” Opt. Commun. **281**(11), 3071–3080 (2008). [CrossRef]

**10. **D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, “Focusing light through random photonic media by binary amplitude modulation,” Opt. Express **19**(5), 4017 (2011). [CrossRef]

**11. **D. B. Conkey, A. N. Brown, A. M. Caravaca-Aguirre, and R. Piestun, “Genetic algorithm optimization for focusing through turbid media in noisy environments,” Opt. Express **20**(5), 4840 (2012). [CrossRef]

**12. **B. R. Anderson, P. Price, R. Gunawidjaja, and H. Eilers, “Microgenetic optimization algorithm for optimal wavefront shaping,” Appl. Opt. **54**(6), 1485 (2015). [CrossRef]

**13. **R. Li, T. Peng, Y. Liang, Y. Yang, B. Yao, X. Yu, J. Min, M. Lei, S. Yan, C. Zhang, and T. Ye, “Interleaved segment correction achieves higher improvement factors in using genetic algorithm to optimize light focusing through scattering media,” J. Opt. **19**(10), 105602 (2017). [CrossRef]

**14. **Y. Wu, X. Zhang, and H. Yan, “Focusing light through scattering media using the harmony search algorithm for phase optimization of wavefront shaping,” Optik **158**, 558–564 (2018). [CrossRef]

**15. **Z. Fayyaz, N. Mohammadian, F. Salimi, A. Fatima, M. R. R. Tabar, and M. R. N. Avanaki, “Simulated annealing optimization in wavefront shaping controlled transmission,” Appl. Opt. **57**(21), 6233 (2018). [CrossRef]

**16. **L. Fang, X. Zhang, H. Zuo, and L. Pang, “Focusing light through random scattering media by four-element division algorithm,” Opt. Commun. **407**, 301–310 (2018). [CrossRef]

**17. **L. Fang, H. Zuo, Z. Yang, X. Zhang, L. Pang, W. Li, Y. He, X. Yang, and Y. Wang, “Particle swarm optimization to focus coherent light through disordered media,” Appl. Phys. B: Lasers Opt. **124**(8), 155 (2018). [CrossRef]

**18. **Z. Wu, J. Luo, Y. Feng, X. Guo, Y. Shen, and Z. Li, “Controlling 1550-nm light through a multimode fiber using a Hadamard encoding algorithm,” Opt. Express **27**(4), 5570 (2019). [CrossRef]

**19. **G. Osnabrugge, L. V. Amitonova, and I. M. Vellekoop, “Blind focusing through strongly scattering media using wavefront shaping with nonlinear feedback,” Opt. Express **27**(8), 11673 (2019). [CrossRef]

**20. **Z. Fayyaz, N. Mohammadian, M. Reza Rahimi Tabar, R. Manwar, and K. Avanaki, “A comparative study of optimization algorithms for wavefront shaping,” J. Innovative Opt. Health Sci. **12**(04), 1942002 (2019). [CrossRef]

**21. **D. Wu, J. Luo, Z. Li, and Y. Shen, “A thorough study on genetic algorithms in feedback-based wavefront shaping,” J. Innovative Opt. Health Sci. **12**(04), 1942004 (2019). [CrossRef]

**22. **S. M. Popoff, G. Lerosey, R. Carminati, M. Fink, A. C. Boccara, and S. Gigan, “Measuring the Transmission Matrix in Optics: An Approach to the Study and Control of Light Propagation in Disordered Media,” Phys. Rev. Lett. **104**(10), 100601 (2010). [CrossRef]

**23. **Y. Choi, T. D. Yang, C. Fang-Yen, P. Kang, K. J. Lee, R. R. Dasari, M. S. Feld, and W. Choi, “Overcoming the Diffraction Limit Using Multiple Light Scattering in a Highly Disordered Medium,” Phys. Rev. Lett. **107**(2), 023902 (2011). [CrossRef]

**24. **D. B. Conkey, A. M. Caravaca-Aguirre, and R. Piestun, “High-speed scattering medium characterization with application to focusing light through turbid media,” Opt. Express **20**(2), 1733 (2012). [CrossRef]

**25. **T. Chaigne, O. Katz, A. C. Boccara, M. Fink, E. Bossy, and S. Gigan, “Controlling light in scattering media non-invasively using the photoacoustic transmission matrix,” Nat. Photonics **8**(1), 58–64 (2014). [CrossRef]

**26. **X. Tao, D. Bodington, M. Reinig, and J. Kubby, “High-speed scanning interferometric focusing by fast measurement of binary transmission matrix for channel demixing,” Opt. Express **23**(11), 14168 (2015). [CrossRef]

**27. **H. Yu, K. Lee, and Y. Park, “Ultrahigh enhancement of light focusing through disordered media controlled by mega-pixel modes,” Opt. Express **25**(7), 8036 (2017). [CrossRef]

**28. **M. Cui and C. Yang, “Implementation of a digital optical phase conjugation system and its application to study the robustness of turbidity suppression by phase conjugation,” Opt. Express **18**(4), 3444 (2010). [CrossRef]

**29. **I. M. Vellekoop, M. Cui, and C. Yang, “Digital optical phase conjugation of fluorescence in turbid tissue,” Appl. Phys. Lett. **101**(8), 081108 (2012). [CrossRef]

**30. **D. Wang, E. H. Zhou, J. Brake, H. Ruan, M. Jang, and C. Yang, “Focusing through dynamic tissue with millisecond digital optical phase conjugation,” Optica **2**(8), 728 (2015). [CrossRef]

**31. **Y. Shen, Y. Liu, C. Ma, and L. V. Wang, “Focusing light through scattering media by full-polarization digital optical phase conjugation,” Opt. Lett. **41**(6), 1130 (2016). [CrossRef]

**32. **Y. Shen, Y. Liu, C. Ma, and L. V. Wang, “Focusing light through biological tissue and tissue-mimicking phantoms up to 9.6 cm in thickness with digital optical phase conjugation,” J. Biomed. Opt. **21**(8), 085001 (2016). [CrossRef]

**33. **Y. Liu, C. Ma, Y. Shen, J. Shi, and L. V. Wang, “Focusing light inside dynamic scattering media with millisecond digital optical phase conjugation,” Optica **4**(2), 280 (2017). [CrossRef]

**34. **D. B. Conkey and R. Piestun, “Color image projection through a strongly scattering wall,” Opt. Express **20**(25), 27312 (2012). [CrossRef]

**35. **H. He, Y. Guan, and J. Zhou, “Image restoration through thin turbid layers by correlation with a known object,” Opt. Express **21**(10), 12539 (2013). [CrossRef]

**36. **L. Wan, Z. Chen, H. Huang, and J. Pu, “Focusing light into desired patterns through turbid media by feedback-based wavefront shaping,” Appl. Phys. B: Lasers Opt. **122**(7), 204 (2016). [CrossRef]

**37. **B. R. Anderson, R. Gunawidjaja, and H. Eilers, “Initial tamper tests of novel tamper-indicating optical physical unclonable functions,” Appl. Opt. **56**(10), 2863 (2017). [CrossRef]

**38. **T. Peng, R. Li, S. An, X. Yu, M. Zhou, C. Bai, Y. Liang, M. Lei, C. Zhang, B. Yao, and P. Zhang, “Real-time optical manipulation of particles through turbid media,” Opt. Express **27**(4), 4858 (2019). [CrossRef]

**39. **A. Konak, D. W. Coit, and A. E. Smith, “Multi-objective optimization using genetic algorithms: A tutorial,” Reliab. Eng. Syst. Saf. **91**(9), 992–1007 (2006). [CrossRef]

**40. **E. G. Bekele and J. W. Nicklow, “Multi-objective automatic calibration of SWAT using NSGA-II,” J. Hydrol. **341**(3-4), 165–176 (2007). [CrossRef]

**41. **G. Vilcot and J.-C. Billaut, “A tabu search and a genetic algorithm for solving a bicriteria general job shop scheduling problem,” Eur. J. Oper. Res. **190**(2), 398–411 (2008). [CrossRef]

**42. **S. Kannan, S. Baskar, J. D. McCalley, and P. Murugan, “Application of NSGA-II Algorithm to Generation Expansion Planning,” IEEE Trans. Power Syst. **24**(1), 454–461 (2009). [CrossRef]

**43. **S. Honda, T. Igarashi, and Y. Narita, “Multi-objective optimization of curvilinear fiber shapes for laminated composite plates by using NSGA-II,” Composites, Part B **45**(1), 1071–1078 (2013). [CrossRef]

**44. **K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput. **6**(2), 182–197 (2002). [CrossRef]

**45. **I. M. Vellekoop, “Feedback-based wavefront shaping,” Opt. Express **23**(9), 12189 (2015). [CrossRef]