## Abstract

Five search algorithms from the literature of black-box optimization were implemented and applied to optical design problems. These algorithms are the Particle Swarm Optimization, Gravity Search Algorithm, Cuckoo Search, Covariance Matrix Adaptation Evolution Strategy and Nelder&Mead Simplex search. The performance of these search algorithms’ implementations was assessed using the BBOB2009 (Black Box Optimization Benchmark) benchmark suite. These algorithms were compared in the context of two optical case studies, one with conventional rotationally symmetric optics and one with freeform optics. A comparison was performed against a commercial optical design software. Finally we provided a simple restart scheme applicable in the workflow of an optical designer. To the best of our knowledge, this is the first in-depth quantitative comparison of optimization algorithms for optical design.

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

## Corrections

Thomas Houllier and Thierry Lépine, "Comparing optimization algorithms for conventional and freeform optical design: erratum," Opt. Express**27**, 28383-28383 (2019)

https://www.osapublishing.org/oe/abstract.cfm?uri=oe-27-20-28383

## 1. Introduction

Search algorithms play a crucial part in optical systems design and the literature shows longstanding efforts in the discipline [1–7]. In addition, the field of freeform optics, which commonly designates the study of optical systems with optical surfaces that do not have an axis of rotational symmetry, allows the design of very compact, very performant systems [8–12]. The use of freeform surfaces creates, as a general rule, optimization problems with more DOF (degrees of freedom) than in conventional optics. This makes the optimization process more difficult. The problem has prompted an interest in investigating various mathematical representations for freeform surfaces in their interplay with the search process [13–16]. We felt however that there was a need to compare quantitatively different search algorithms when applied to optical design examples, both conventional and freeform. We were encouraged by the good results of Menke’s research in applying Particle Swarm Optimization to optical design [17]. An investigation similar to our own, limited to conventional optics and with less detailed quantitative assessments can be found in [18].

Our aim is mainly to give optical designers an insight into what they could gain from using different search algorithms on their optimization problems. We would also like to show search algorithms specialists how their field is applied to optics, so that hopefully they will be inspired to make more recent, more performant advances available to the optical design community. In contrast with typical optical design methods where a starting point is required, we investigate algorithms which only need a defined search space. This brings the optical systems optimization closer to a fully automated process.

Other research, complementary to our own in the attempt to better the optimization of freeform optics, can be found in [19–22]. These provide much needed analytical and systematic methods for the design of freeform optics and focus on providing good starting points. In contrast, we take the more common raw raytracing approach of optical design software and focus on a more automatic optimization of optical systems by the use of various search algorithms.

In Section 2, we give an overview of where our work is situated in a typical optical design workflow. In Section 3, we describe the search algorithms we implemented using resources from the “Derivative-free optimization” and “Black-box optimization” literature. In Section 4, we then assessed their performance on a benchmark set from the literature. In Section 5, the performance of these algorithms was measured on two optical systems and compared with the performance of the search methods in the commercial optical design program OpticStudio. For this purpose, an independent raytracing tool was developed. Finally we give in Section 6 a practical template for a basic restart scheme and apply it to one of the optical systems.

## 2. Optical design workflow

A typical optical design project, like designing a photographic lens for example, can often be decomposed from the optical designer’s point of view in the way drawn in Fig. 1.

The MF (Merit Function) corresponds to a scalar performance score that the optical designer defines. The MF is a function of a vector of the optical system parameters. These parameters can be radii of curvature of lenses, lens thicknesses, geometrical positions of elements in the system, glass refraction indices, etc. The optical designer builds the MF to take into account basic specifications like the focal length of the optical system, the F-number, image quality etc. Advanced specifications like manufacturing tolerances or low levels of stray light are not commonly included in this “prototyping” MF, as they are too computationally costly. Each MF evaluation in the “prototyping” phase typically requires a raytracing simulation through the optical system of a few hundred rays.

Minimizing the MF means finding the set of parameters that gives the best performance according to the defined MF. Optical design software almost always includes an integrated optimization engine, often with different search methods. These search methods are not publicly disclosed in full details, at least not in commercial software.

We believe the prototyping phase to be crucial in the workflow of optical design work. In the present paper, we focus on the search algorithms used to minimize the MF.

## 3. Implemented search algorithms

We implemented 5 search algorithms from the derivative-free optimization literature. We also implemented a simple Monte Carlo search (abbrev. MC), which serves as a baseline for “worst” search algorithm. Our implementations may deviate from the canonical versions of each of the search algorithms, as numerous variations for each algorithm already exist in the literature.

Each algorithm has a set of tunable parameters, whose values can have a very significant impact on performance and must be adapted to the particular search problem the algorithm runs on. We describe each algorithm only sufficiently for the readers to understand general ideas and where the parameters come into play in our implementations. For the sake of conciseness, we refer the readers to the articles given as references for complete definitions of each algorithm.

#### 3.1. Particle Swarm Optimization

The Particle Swarm Optimization (PSO) involves a swarm of particles, each having a position in the search space and a corresponding function value [23–25]. Each particle remembers the best position it has personally visited (cognitive component). It also knows the best position currently occupied by an individual of the swarm (social component). We implemented the so-called *gbest* swarm topology [23], where the social component is constituted only by the current best particle in the swarm. We note that PSO has already been applied to optical design with success [17].

The speed of each particle is updated at each iteration step in a loop by the formula [25]:

- • $\overrightarrow{{v}_{next}}$: The updated velocity vector of the particle.
- • $\overrightarrow{{v}_{last}}$: The velocity vector of the particle from the previous iteration. Also called inertia.
- • $ran{d}_{1},ran{d}_{2}$: Two different scalar random numbers from a uniform random number generator, from 0 to 1.
- • $\overrightarrow{{D}_{cog}}$: The distance vector between the particle and the best position it has personally visited.
- • $\overrightarrow{{D}_{soc}}$: The distance vector between the particle and the best particle in the swarm.
- • $\alpha ,\beta $: Two parameters for the algorithm. The recommended values when no information is available on the function to be minimized are $\alpha =0.7298$ and $\beta =1.4961$ [25].

Another parameter for the PSO is of course the number of particles we choose to use, which we will call ${n}_{parts}$.

#### 3.2. Gravity Search Algorithm

The Gravity Search Algorithm (GSA) [26,27], like the PSO, is a particle-based search method. Each particle is seen as a planet in an n-body gravitational system. A mass is assigned to each particle in relation to the function value at its position. To make the algorithm better suited to minimization, only the heaviest planets (with better function value) exert an attraction pull. The number of planets in this active set is${K}_{best}\left(t\right)$out of ${n}_{parts}$planets in total. Furthermore, the size of the “heaviest planets” set decreases linearly with time, with only one planet attracting all the others by the end of the optimization run. The gravitational constant $G\left(t\right)$of this fictitious universe is also dynamically decreased as a function of current iteration number.

We use the following update equations, as a function of $t$ the iteration number, with ${n}_{iter}$the total number of iterations in a given optimization run:

_{${G}_{0}$} and $\alpha $are parameters for our implementation of GSA.${K}_{best}\left(t\right)$decreases linearly with $t$ from ${n}_{parts}$ to $1$. The rest of the implementation follows closely the one outlined in [26].

#### 3.3. Cuckoo Search

The Cuckoo Search (CS) involves particles viewed as cuckoos [28–31], who jump from position to position in the search space using a random walk modeled by Levy flights and creating ‘nests’ in the process at the new positions. These nests have in turn a chance of being discovered by a virtual host bird with a fixed probability (thus mimicking the parasitic brooding behavior of cuckoos). This prompts another, different, random walk using a uniform probability density. We followed closely the MatLab implementation provided in [30].

For simplicity’s sake, let us consider each cuckoo as a particle that can probe around its position using two random walks: one being global, a Levy flight, applied to every particle, and the other being local and involving only a fraction of the particles. The position of a particle is updated only if it finds a better position while probing.

The position $\overrightarrow{{x}_{probe,L\xe9vy}}$ probed by a particle by (global) Levy flight from position $\overrightarrow{x}$ is:

And the position $\overrightarrow{{x}_{probe,walk}}$ probed by a particle via local random walk, with probability ${p}_{a}$ of happening is:

- • $\sigma $: A fixed parameter for Levy flight. $\sigma =0.6965745$.
- • $N\left(0,1\right)$: Random number drawn from a normal distribution with zero mean and unit standard deviation. Each occurrence in the equations is a different random number.
- • $rand$: Random number from a uniform distribution between 0 and 1.
- • $\overrightarrow{{x}_{best}}$: Position of the current best particle in the swarm.
- • $\overrightarrow{{x}_{shuff1}},\overrightarrow{{x}_{shuff2}}$: Positions of two particles in the swarm taken at random.

The parameters for this particular implementation of CS are then ${n}_{parts}$(the number of particles to use) and ${p}_{a}$.

#### 3.4. Covariance Matrix Adaptation Evolution Strategy

The Covariance Matrix Adaptation Evolution Strategy (CMAES) is a population-based search algorithm [32–35]. The general idea behind this algorithm is to draw generations of points at random positions using a probability density model. The probability density, which is a multivariate normal distribution, characterized by a mean and a covariance matrix, is recalibrated every few generations to reflect the results obtained by the previous generations of points and guide the population towards a minimum. The technical details of implementation are too involved to be explained here in a concise way, but the reader can refer to a very detailed guide published by Hansen on his algorithm [33], which we tried to follow as closely as possible. We use all the default parameters detailed in [33], leaving only $\sigma $, the initial step-size for mutation, as a control knob for the user.

#### 3.5. Nelder&Mead simplex search

The simplex search developed by Nelder and Mead (SPX) is one of the most common and better known search algorithm [36–38]. It is not stochastic, meaning that it has no random component and will always give the same result with a fixed set of starting conditions. The method involves a simplex, a geometrical shape which has $D+1$ vertices, $D$ being the dimensionality of the search space. The simplex search is based on a simple set of rules involving the function values at the positions of three points: the best, second worst and worst points of the simplex. A flowchart can be found in [36]. These rules give rise to four different possible geometric transformations of the simplex called reflection, expansion, contraction and multiple contraction (or general contraction or shrinking). The first three are homothetic transformations of the worst point in the simplex around the barycenter (non-weighted) of the rest of the points in the simplex. The new position $\overrightarrow{{x}_{new}}$ is obtained from the worst point position $\overrightarrow{{x}_{worst}}$ and barycenter position $\overrightarrow{{x}_{bary}}$ via:

The multiple contraction transformation shrinks all the points (here each noted $\overrightarrow{{x}_{last}}$) in the simplex towards the best point $\overrightarrow{{x}_{best}}$ by the ratio $\beta $ as such:

Our implementation follows closely the one from Numerical Recipes [38], and we kept the canonical transformation parameters for reflection, expansion and contraction of $\alpha =-1,2,0.5$respectively and$\beta =0.5$. Our implementation of simplex search without restart strategy has no user parameter.

#### 3.6. Other implementation concerns

Comments must be added on three issues, namely search starting point, search space bounds management and restart strategies.

Our search methods do not start at a given starting point defined by the user, as is usually the case in optical design software. The user only provides a minimum and maximum value for each variable, defining a search space boundary. We initialize every search method by drawing random points in the search space within the user-provided bounds, until we gather a sufficient number of feasible points to start the search. Feasible points are positions where the merit function can be evaluated. It is often the case in optical design that merit functions have unevaluable regions (take the case of rays missing elements entirely for example). For PSO, GSA and CS, we draw all the particles in the initial swarm at random, they are thus spread all over the search space. For CMAES and SPX, we draw a single initial random point and initialize the other particles at starting positions close to the initial point using a small fraction of the search space size in each dimension to spread apart the points while keeping them relatively close to the initial point.

The algorithms we present here generally do not have bound management methods in their canonical version. We implemented strict bound enforcing for each algorithm, we never evaluate the merit function outside of the user-provided bounds and set the eventual infringing points at the nearest bound instead. For CMAES, we used another method, we strongly penalized points falling outside of bounds by assigning them a huge constant value instead of evaluating the merit function. This is not ideal, but it is simple. A more thought out way of implementing bound enforcement for CMAES is discussed in [33], it was too involved for us to implement easily.

An “evaluation budget” is a maximum number of function evaluations a search method is allowed to perform before stopping. For simplicity’s sake and for ease of comparison, we did not provide any stop condition to our search methods other than the function evaluation budget. To highlight the effectiveness of restart strategies, we do however show throughout this paper the effect of adding a “stop and restart” condition to the simplex search. We designate by “restart strategies” the methods one can use to abort and restart an optimization run before having spent the whole evaluation budget, when an early convergence is detected for example. The remaining budget is used to find potentially better minima elsewhere in the search space. This stop condition, which is the one used in [38], is as follows:

- • ${y}_{hi}$: Value of the merit function at the highest (worst) point in the simplex.
- • ${y}_{lo}$: Value of the merit function at the lowest (best) point in the simplex.
- • $\epsilon $: Tiny numerical constant to avoid division by zero.
- • ${f}_{tol}$: Stop condition parameter. Note that it is independent of actual merit function value.

If the stop condition is reached, the search method aborts the current simplex run, saving the best result it had reached beforehand, and reinitializes a simplex randomly within the search space. This continues until the search method has spent its function evaluation budget.

## 4. Results on the BBOB2009 test suite

Given that we implemented our own versions of the search algorithms cited in the present paper, we needed to provide an objective point of reference for their performance. We chose to run the BBOB2009 (Black-Box Optimization Benchmarking) [39], which we re-implemented using the very useful online resources provided by the authors of this test suite [40].

The BBOB2009, consists of 24 different function types of $D$ dimensions [41], which have a search space of ${\left[-5;5\right]}^{D}$, a random minimum position $\overrightarrow{{x}_{opt}}$, and a random minimum value of ${f}_{opt}\in \left[-1000;1000\right]$ drawn from a Cauchy distribution. We draw 15 optimization problems for each function type and for each dimensionality $D$ considered, $D=\left[2,3,5,10,20,40\right]$, thus creating $24\times 15\times 6=2160$ independent optimization problems for our search algorithms to run on.

We run all the implemented search algorithms on each of the 2160 optimization problems and we save the best found minimum so far ${f}_{best}$ as a function of number of function evaluations for each problem. For the optimization problem $i$ having a global minimum ${f}_{opt,i}$, we compute $\Delta {f}_{i}={f}_{best,i}-{f}_{opt,i}$. The optimization problem is considered as solved when $\Delta {f}_{i}$ is below an arbitrary threshold$\Delta {f}_{thresh}$. A measure of performance for a search algorithm is then the proportion $R\in \left[0;1\right]$of the 2160 problems which have reached the threshold as a function of two parameters, first the number of function evaluations divided by $D$ the dimensionality of the problem ${n}_{eval,D}=\frac{{n}_{eval}}{D}$and secondly the threshold $\Delta {f}_{thresh}$. This performance metric can be thought of as a 2D surface, the altitude of the surface being the proportion $R$ of solved problems as a function of the two parameters ${n}_{eval,D}$ and $\Delta {f}_{thresh}$.

The results are synthesized in a way very similar to that adopted by the authors of the BBOB2009 [42], namely we draw two graphs, which are cuts of the 2D surface mentioned above:

- 1. Proportion of solved problems as a function of number of function evaluations per search space dimension with a fixed threshold distance to the global minimum value. This is in Fig. 2.
- 2. Proportion of solved problems as a function of threshold distance to the global minimum value at a fixed number of function evaluations. This is in Fig. 3.

A good algorithm on Fig. 2 will have a curve that is towards higher$R$values (more problems solved) as early as possible in terms of evaluation count${n}_{eval,D}$. On Fig. 3, a good algorithm will exhibit a curve closer to the top left of the diagram, meaning that it will solve many problems with a solution closer to the global minimum. The $\Delta {f}_{thresh}$value can be seen as a difficulty setting meaning “solve the problems to within $\Delta {f}_{thresh}$of the global minimum”.

The parameterization for all the algorithms we ran on the BBOB2009 can be found in Table 1. It is important for the readers to note that we do not pretend to compare in absolute terms the merits of the search algorithms we chose to implement, we merely state the performance of our own implementations, which could be faulty, with a particular set of search parameters, which can be ill-chosen for a given problem.

Our results show that our implementations of search algorithms do not perform very well against state of the art methods [42], which is not surprising given that we used very bare versions, with little attention to parameter tuning and no restart strategies. Despite the poor performance of our implementations against the state of the art, we will see that they nonetheless perform satisfactorily in optimizing optical systems when compared with a commercially available tool.

## 5. Application to two optical systems

We measure the performance of our implementations of search algorithms on a typical merit function for the design of optical systems. We use two different optical systems to perform the test. They are both eyepieces sharing a number of common specifications. We also compare the performance of the implemented algorithms against the search methods included in OpticStudio (Zemax). We will hereafter abbreviate ‘OpticStudio’ by ‘OS’. We used our own independent raytracing engine to compute the results of the implemented search algorithms.

#### 5.1. Simulation setup

We describe the merit function we used as a measure of performance to optimize in the two optical case studies.

### 5.1.1. Merit function definition

The MF we defined takes into account two typical optical quality criteria:

- • Spot size: Each bundle of rays coming from a given field at the entrance of the optical system must form the smallest possible spot on the image plane. This is an image quality criterion, if it is low then the image will be sharp, otherwise it will be blurry.
- • Image position: Each bundle of rays coming from a given field at the entrance of the optical system must focus at a given position on the image plane. This is a control over the focal length as well as the distortion in the optical system.
Our definition of the MF is then:

- • $Nfields$: The number of field points taken into account for the computation. We typically sample the object field space with up to a dozen points.
- • ${x}_{i},{y}_{i}$: The coordinates of impact of individual rays on the image plane.
- • $\overline{x},\overline{y}$: The coordinates of the centroid of the spot.
- • ${\overline{x}}_{trgt},{\overline{y}}_{trgt}$: Target coordinates on the image plane for the centroid of the ray bundle.

### 5.1.2. Raytracing implementation details

We have developed simple raytracing routines. They perform raytracing on optical systems with surfaces positioned arbitrarily in 3D space and support spherical as well as XY polynomials surface representations. They are limited to monochromatic raytracing, but support the sampling of multiple field points as well as arbitrary aperture stop sampling. We also limited ourselves to the case of aperture stops located in the object space of optical systems. For the two case studies in the present paper, we used an aperture stop sampling consisting of a uniform Cartesian grid cropped by the aperture shape, which is a disk.

#### 5.2. Characteristics of the case studies

We specify the optical characteristics of our two case studies as they would be in a typical optical design preliminary work. The two optical systems share the same set of basic target specifications but they differ greatly in shape and complexity.

### 5.2.1. Specifications for the two optical systems

We give the optical specifications for the two eyepieces we mean to optimize. The first eyepiece, drawn in Fig. 4, has 5 spherical lenses in an axio-symmetric configuration. This is a typical eyepiece for use with amateur telescopes for example. This system has 21 DOF and is comparatively easier to optimize: it is easy to reach realistic systems with good optical quality. The second eyepiece, drawn in Fig. 5, is a freeform prism with 2 optical surfaces described by fifth order XY polynomials (defined exactly as the surface type “Extended Polynomial” in OS). The whole prism is plano-symmetric. This system has 22 DOF. System 2 is inspired from the general shape of the system in [9]. The optical specifications are found in Table 2. The aperture stop diameter is the diameter of the eye pupil. The focal length is constrained by the field of view (FOV) in the eye space and size of the object. We give the design parameters specific to each system in Table 3. The raytracing for both system is done monochromatic and the refraction indices for the various glasses are indicated in Fig. 4 for system 1 and Fig. 5 for system 2. We note that it is preferable to use curvatures instead of radii of spheres for the specification of spherical surfaces. These quantities are inverse of each other. Indeed, radii give rise to unevaluable MF around their zero value, making a “hole” in the search space that hampers the performance of the search algorithms we implemented.

### 5.2.2. Obtaining results on a commercial software for comparison

We recreated the same minimization problems in OS. The surfaces of the optical system are modeled using the native functionalities of the software. The raytracing is also done using the native OS capabilities. The merit function is computed ray by ray using so-called “optimization operands” so as to mirror exactly Eqs. (9)–(11).

To enforce strict bounds on the search space, we resorted to the same type of method as we applied on the CMAES, namely we penalized heavily the MF variables that wander outside the user-defined bounds. We have$M{F}_{with\_bounds}=MF+bound\_penalty$, the search algorithms in OS ran on$M{F}_{with\_bounds}$, but we always report $MF$ since it is the true measure of the optical systems’ fitness. This may very well have impacted the performance exhibited by the OS search routines. However, we consider the comparison is still fair if the design problem the optical designer sets out to solve is “Find the best optical system according to the MF with system variables strictly within the bounds”.

### 5.2.3. Validation of the raytracing merit function

To assess the exactness of the comparison between our raytracing routines and OS, we compared the values of the MF for 1000 random points within the search space bounds for both of our case studies. The MF is evaluable at these 1000 positions. This test ensures the exactness of our raytracing routines and the concordance between the MF definitions in our routines and in OS.

For each of these 1000 points, we computed the difference $\left|M{F}_{author}-M{F}_{zmx}\right|$ between the two results. We report the maximum and mean values of this quantity across our test set in Table 4. The errors are negligible when compared to the lowest MF values found, the lowest MFs for System 1 being in the order of ${10}^{-5}$ and for System 2 in the order of ${10}^{-4}$ (see Section 5.3).

#### 5.3. Results of running the implemented algorithms on the two case studies

We launch, for each case study, 100 independent optimization runs of *at least* 5000 function evaluations for each search algorithm, including those in OS. The OS searches were started with a fixed list of 100 random, evaluable points within the user-provided bounds. We ran two types of search routines provided by OS, which are called: Local optimization Orthogonal Descent (ZLOD) and Hammer OD (ZHOD). The Hammer is a global optimizer. A DLS (Damped Least Squares) mode also exists in OS. However it could not be brought to perform at least 5000 MF evaluations for each optimization run, hence the comparison with other search methods would not have been fair to OS. We report that the performance of the DLS mode appears to be worse than that of the OD mode for both our case studies in test runs that are not reported here. OS provides another search routine called “Global Optimization”. We did not test it due to practicality issues. This routine carries out the optimization over many different files, making the results difficult to track.

The results can be found in Fig. 6 and Fig. 7 for the 5 lenses eyepiece (System 1) and in Fig. 8 and Fig. 9 for the freeform prism eyepiece (System 2). The results are formatted in the same manner as those in the BBOB2009 benchmark, except for the y axis of each graph. Rather than expressing the ratio of successful number of runs (reaching below$\Delta {f}_{thresh}$), we give the raw number of successful runs. The total number of launched runs is 100 for each search algorithm. The SPX with restart can have a higher number of successful runs, since it has the ability to restart as many times as it wants within its 5000 evaluations budget. We count each convergence of the SPX method with restart as a complete run since it provides each time a viable optical system. The parameterization of the search algorithms for the two case studies is given in Table 5.

Since the global minima of the MF for the two case studies are unknown, we use the raw MF value as criterion for$\Delta {f}_{thresh}$. Namely, we compute $\Delta f={f}_{best}-{f}_{opt}$with${f}_{opt}=0$, ${f}_{best}$being the minimum MF value found so far in a given optimization run.

Note that whenever an algorithm does not appear in a graph, it means that its performance was too low. Let us analyze the results from the perspective of an optical designer. On Fig. 6, a designer running 100 optimization runs on System 1 would find 50 optical systems with a MF below ${10}^{-2}$ in about 50 evaluations per problem dimension (~1050 MF evaluations in total) for the SPX with restart versus about 90 evaluations per dimension (~1890 in total) for the ZLOD. This means that, for the particular problem of obtaining 50 optical systems of a given quality out of a 100 tries, the SPX with restart can be seen as a little less than twice as fast as the ZLOD. Likewise on Fig. 8, using the PSO over the ZHOD to obtain 50 optical systems with a MF below ${10}^{-2}$ would be about 2 times faster (56 evaluations per dimension for PSO versus 106.5 evaluations per dimension for ZHOD).

Note that, for every example presented here, the computational cost of evaluating the MF was always hugely predominant over the computational cost of managing the searches. This means that the time it takes a designer to run the optimizations is directly proportional to the number of MF evaluations needed.

Now let us look at the ability of the algorithms to reach low global minima, *ie* find better optical systems for the optical designer. If we were to perform a hundred runs of 5000 MF evaluations for each algorithm and look at the quality of the 5 best systems we obtained out of 100, for System 1 on Fig. 7 we would get 5 systems with a MF better than ${4.10}^{-5}$using SPX and 5 systems better than ${2.10}^{-4}$with ZLOD, the 5 best systems of the SPX runs are 5 times better in terms of MF than the 5 best systems of the ZLOD runs. On Fig. 9, for System 2, we can compare the performance of the 5 best systems for PSO and ZHOD. The 5 best systems for PSO are below $8\times {10}^{-4}$ in MF value, while they are below $1.6\times {10}^{-4}$ for ZHOD, the ZHOD shows a MF value 5 times better than the PSO.

We note that the SPX that includes a restart strategy is mostly preferable to the SPX without restart strategy. This is not surprising given that a simplex cannot get out of a local minimum. Once the simplex is trapped in a local minimum, it is useless to let it run further. In consideration of this result, the implementation of restart strategies for the other algorithms, detecting similar “stalled” system states, seems very desirable.

To synthesize the present results, for System 1, it is preferable to use the SPX with restart, the SPX or the PSO over the rest of the algorithms. For System 2, the ZHOD and ZLOD are almost always preferable, and the PSO could be used as a complementary method to obtain a greater result variety.

The optical performance of the best found systems in each case study is satisfactory. The best results of System 1 have RMS spot radii of about 5-10 µm and maximum distortion of about 3.5%. For System 2, these values become 10-15 µm and 2.5%. System 2 could for example be used in virtual or augmented reality applications, as it is very small and has a shape compatible with injection molding. The comparison of both case studies highlights the advantages of freeform optics for miniaturization.

The performance of search algorithms depends heavily on the type of problem being optimized. We see this when comparing the relative performance of the search algorithms on the two present case studies. For example, SPX performs very well comparatively to the other methods on System 1, but its performance is underwhelming on System 2. This phenomenon is known as the No Free Lunch theorem of optimization [43]. It states, in layman’s terms, that no search method can outperform another search method on all MFs. In fact, if we compare search method A and search method B, there are always as many MFs for which A is better than B as there are for which B is better than A. It would be helpful to the optical design community to study the performance of search methods depending on the types of optical systems and/or on MF definitions.

## 6. A basic restart scheme for practical use

Using pure search methods facilitates comparisons for the purpose of publishing the present results. However, for practical purposes, there is a great performance gain in using even simple restart schemes, as suggested in [37]. It allows a more user-friendly interface with the search algorithms, as the optical designer can specify parameters like the desired local algorithm (out of the ones in Section 3), the number of best systems to save and display at the end of the optimization process, the evaluation budget of the runs (which can be thought of as the time spent in each run), the number of runs etc. An advanced prototyping scheme, closer to the state of the art in general engineering, can be found in [44].

The flowchart for our simple strategy is found in Fig. 10. The general idea is that we first run a “local optimizer”, taken from the list of pure search methods, until a stop criterion is met or an evaluation budget is spent. Then we take the result outputted by the local optimizer and run it through a simplex optimizer to reach the bottom of the particular valley we have landed in. We compensate in this way the tendency of some optimizers to hover about the final result. Then we compare the final result with a buffer of past results and save it only if the MF value is good enough and if the system is sufficiently different from the results already stored in the buffer (according to a distance metric between optical systems). The so-called “Global method inputs” are:

- •
**Merit function**: The function for which we want to have a list of minima. - •
**MF variables bounds**: The limits of the search space. The optimizers strictly enforce the limits of the search space. - •
**Global budget**: Global merit function evaluation budget for the whole search. We stop the global search once this budget is exhausted. - •
**Local budget**: Local merit function evaluation budget for one local optimization round. The local optimizer must stop when its budget is exhausted or when another optimizer-specific stop criterion is reached. - •
**Size of results buffer**: This is how many solutions we wish to keep track of and present to the user at the end of the global search. Setting this to 1 makes the global search output only the best global minimum it found. - •
**Distance metric threshold**: Threshold below which two points in the solution space are considered the same according to some distance metric and thus not stored twice in the solution buffer. The simplest such metric is the Euclidian norm. But it is not suitable for ill-conditioned problems, and is unsuitable for large dimensionality. A more advanced distance metric between optical systems could be devised by the optical designer on a case by case basis. - •
**Other local optimizer specific inputs:**Stop criteria, population size, algorithm parameters etc.

The output of the method is a list of the best optical systems found that are at the same time sufficiently different from one another according to a provided distance metric. An example of such a list of results running on the MF of System 1 is provided in Fig. 11. We think such a list of results is very useful to the optical designer, as it allows exploration of the search space in preliminary design phases. This is a feature that also exists in commercial software. It is called “Global Optimization” in OS.

## 7. Conclusion

We demonstrated useful performance improvements for the search algorithms we implemented against commercial software on optimizing two optical case studies. These performance gains exist despite the poor performance of these search algorithms on a standard benchmark suite when compared to the state of the art. This seems to indicate that further improvements could be made in search algorithms for optical design by applying state of the art search methods. In the future, we think massive improvements could be made by applying state of the art machine learning methods and artificial intelligence to the problem of designing optical systems, particularly in the area of structural changes to optical systems, *eg* adding or removing elements from a system.

We have applied a rigid comparison method between search algorithms (the methodology of the BBOB2009) for given optical systems with given MF definitions. This framework can be applied in future works in order to perform quantitative assessments of the performance of other search algorithms in comparison with the present work.

We have also shown on these two case studies that one should not expect a given search algorithm to outperform all the others for every problem and that it would be beneficial to select the correct algorithm for each problem. We gave an example of a simple restart scheme that can be used in practice by optical designers for global optimization purposes. In combination with internal raytracing software tools that we have developed, we have effectively created a simple yet functional optical design software adapted to higher dimensionality and difficult optimization problems, which include freeform optical design.

## 8. Technical software details

Software implementations were done using Steel Bank Common Lisp (SBCL) [45] version 1.4.4, which is a Common Lisp compiler. We used OpticStudio Professional 18.7 for System 1 and Premium 18.9 for System 2, belonging to Sophia Engineering. We ask our readers to avoid treating the present results as a rigorous benchmark of the commercial software.

## Funding

Sophia Engineering; Thales Alenia Space; Association Nationale de la Recherche et de la Technologie (ANRT).

## Acknowledgements

We wish to thank H. Benard and F. Keller at Thales Alenia Space (Cannes, France) for their part in the PhD work. We also thank G. Cassar and M. Guilhem at Sophia Engineering for proofreading this article. We thank especially M. Guilhem for the many hours of discussion about raytracing: We also thank the reviewers of this paper for their insightful remarks.

## References

**1. **E. Glatzel and R. Wilson, “Adaptive automatic correction in optical design,” Appl. Opt. **7**(2), 265–276 (1968). [CrossRef] [PubMed]

**2. **M. J. Kidger, “Use of the Levenberg-Marquardt (damped least-squares) optimization method in lens design,” Opt. Eng. **32**, 1731–1740 (1993). [CrossRef]

**3. **D. Shafer, “Global optimization in optical design,” Comput. Phys. **8**, 188–195 (1994). [CrossRef]

**4. **K. Höschel and L. Vasudevan, “Genetic algorithms for lens design: a review,” J. Opt. **48**(1), 134–144 (2019). [CrossRef]

**5. **S. Banerjee and L. Hazra, “Experiments with a genetic algorithm for structural design of cemented doublets with prespecified aberration targets,” Appl. Opt. **40**(34), 6265–6273 (2001). [CrossRef] [PubMed]

**6. **L. Hazra and B. Saswatee, “Genetic algorithm in the structural design of Cooke triplet lenses,” Proc. SPIE **3737**, 172–180 (1999). [CrossRef]

**7. **X. Cheng, W. Yongtian, H. Qun, and I. Masaki, “Global and local optimization for optical systems,” Optik (Stuttg.) **117**, 111–117 (2006). [CrossRef]

**8. **J. P. Rolland, K. Fuerschbach, G. E. Davis, and K. P. Thompson, “Pamplemousse: The optical design, fabrication, and assembly of a three-mirror freeform imaging telescope,” Proc. SPIE **9293**, 92930L (2014).

**9. **Z. Shen, J. Yu, Z. Song, L. Chen, Q. Yuan, Z. Gao, S. Pei, B. Liu, and J. Ye, “Customized design and efficient fabrication of two freeform aluminum mirrors by single point diamond turning technique,” Appl. Opt. **58**(9), 2269–2276 (2019). [CrossRef] [PubMed]

**10. **Q. Meng, H. Wang, W. Liang, Z. Yan, and B. Wang, “Design of off-axis three-mirror systems with ultrawide field of view based on an expansion process of surface freeform and field of view,” Appl. Opt. **58**(3), 609–615 (2019). [CrossRef] [PubMed]

**11. **Z. Li, X. Liu, F. Fang, X. Zhang, Z. Zeng, L. Zhu, and N. Yan, “Integrated manufacture of a freeform off-axis multi-reflective imaging system without optical alignment,” Opt. Express **26**(6), 7625–7637 (2018). [CrossRef] [PubMed]

**12. **J. Zhu, H. Wei, Z. Xiaodong, and J. Guofan, “Design of a low F-number freeform off-axis three-mirror system with rectangular field-of-view,” J. Opt. **17**(1), 015605 (2015). [CrossRef]

**13. **A. Yabe, “Representation of freeform surfaces suitable for optimization,” Appl. Opt. **51**(15), 3054–3058 (2012). [CrossRef] [PubMed]

**14. **A. Broemel, H. Gross, D. Ochse, U. Lippmann, C. Ma, Y. Zhong, and M. Oleszko, “Performance comparison of polynomial representations for optimizing optical freeform systems,” Proc. SPIE **9626**, 96260W (2015).

**15. **A. Broemel, C. Liu, Y. Zhong, Y. Zhang, and H. Gross, “Freeform surface descriptions. Part II: application benchmark,” Adv. Opt. Technol. **6**(15), 337–347 (2017).

**16. **M. I. Nikolic, B. Pablo, N. Bharathwaj, A. G. Dejan, L. Jayao, and M. Juan Carlos, “Optical design through optimization for rectangular apertures using freeform orthogonal polynomials: a case study,” Opt. Eng. **55**(7), 071204 (2016).

**17. **C. Menke, “Application of particle swarm optimization to the automatic design of optical systems,” Proc. SPIE **10690**, 106901A (2018). [CrossRef]

**18. **D. Vasiljevic, *Classical and Evolutionary Algorithms in the Optimization of Optical Systems* (Springer Science & Business Media, 2012).

**19. **A. Bauer, E. M. Schiesser, and J. P. Rolland, “Starting geometry creation and design method for freeform optics,” Nat. Commun. **9**(1), 1756 (2018). [CrossRef] [PubMed]

**20. **J. C. Papa, J. M. Howard, and J. P. Rolland, “Starting point designs for freeform four-mirror systems,” Opt. Eng. **57**(10), 101705 (2018). [CrossRef]

**21. **T. Yang, G. F. Jin, and J. Zhu, “Automated design of freeform imaging systems,” Light Sci. Appl. **6**(10), e17081 (2017). [CrossRef] [PubMed]

**22. **L. I. Jun, W. Huang, and F. Hongjie, “A novel method for finding the initial structure parameters of optical systems via a genetic algorithm,” Opt. Commun. **361**, 28–35 (2016). [CrossRef]

**23. **R. Eberhart and J. Kennedy, “A new optimizer using particle swarm theory,” Proceedings of the Sixth International Symposium on Micro Machine and Human Science (1995), pp. 39–43. [CrossRef]

**24. **C. Leboucher, S. Hyo-Sang, C. Rachid, L. M. Stéphane, S. Patrick, F. Mathias, T. Antonios, and K. Alexandre, “An Enhanced Particle Swarm Optimization Method Integrated With Evolutionary Game Theory,” IEEE Trans. Games **10**, 221–230 (2018). [CrossRef]

**25. **T. Zeugmann, P., Poupart and J. Kennedy, “Particle swarm optimization,” in *Encyclopedia of Machine Learning*, C. Sammut and G. I. Webb, ed. (Springer Science & Business Media, 2011).

**26. **E. Rashedi, N.-P. Hossein, and S. Saeid, “GSA: a gravitational search algorithm,” Inf. Sci. **179**(13), 2232–2248 (2009). [CrossRef]

**27. **S. Özyön, Y. Celal, and T. Hasan, “Incremental gravitational search algorithm for high-dimensional benchmark functions,” Neural Comput. Appl. **29**, 1–25 (2018).

**28. **X.-S. Yang and D. Suash, “Engineering optimisation by cuckoo search,” Math. Model. Num. Opt. **1**(4), 330–343 (2010). [CrossRef]

**29. **A. H. Gandomi, X.-S. Yan, and A. H. Alavi, “Cuckoo search algorithm: a metaheuristic approach to solve structural optimization problems,” Eng. Comput. **29**, 17–35 (2013). [CrossRef]

**30. **X.-S. Yang and D. Suash, “Cuckoo search: recent advances and applications,” Neural Comput. Appl. **24**, 169–174 (2014). [CrossRef]

**31. **G.-G. Wang, D. Suash, A. H. Gandomi, Z. Zhaojun, and A. H. Alavi, “Chaotic cuckoo search,” Soft Comput. **20**, 3349–3362 (2016). [CrossRef] [PubMed]

**32. **N. Hansen and A. Ostermeier, “Completely derandomized self-adaptation in evolution strategies,” Evol. Comput. **9**(2), 159–195 (2001). [CrossRef] [PubMed]

**33. **“The CMA evolution strategy: A tutorial,” 2016, https://arxiv.org/pdf/1604.00772.pdf.

**34. **N. Hansen, “Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed,” Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference: Late Breaking Papers (2009), pp. 2389–2396. [CrossRef]

**35. **A. Auger and H. Nikolaus, “A restart CMA evolution strategy with increasing population size,” 2005 IEEE Congress on Evolutionary Computation (IEEE, 2005), pp. 1769–1776. [CrossRef]

**36. **J. A. Nelder and R. Mead, “A simplex method for function minimization,” Comput. J. **7**, 308–313 (1965). [CrossRef]

**37. **A. Luersen, L. R. Marco, Rodolphe, and G. Frédéric, “A constrained, globalized, and bounded Nelder–Mead method for engineering optimization,” Struct. Multidis. Optim. **27**, 43–54 (2004).

**38. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes 3rd edition: The Art of Scientific Computing* (Cambridge University, 2007).

**39. **N. Hansen, F. Steffen, R. Raymond, and A. Auger, “Real-parameter black-box optimization benchmarking 2010: Experimental setup,” Research report RR-7215, INRIA, 2010.

**40. **COCO, (COmparing Continuous Optimisers) homepage, (INRIA, 2019), http://coco.gforge.inria.fr/.

**41. **N. Hansen, F. Steffen, R. Raymond, and A. Auger, “Real-parameter black-box optimization benchmarking 2009: Noiseless functions definitions,” Research report RR-6829, INRIA, 2009.

**42. **N. Hansen, A. Auger, R. Ros, S. Finck, and P. Pošík, “Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009,” Proceedings of the 12th Annual Conference Companion on Genetic and Evolutionary Computation (2010). pp. 1689–1696. [CrossRef]

**43. **D. H. Wolpert and W. G. Macready, “No free lunch theorems for optimization,” IEEE Trans. Evol. Comput. **1**, 67–82 (1997). [CrossRef]

**44. **A. Hagg, A. Asteroth, and T. Bäck, “Prototype discovery using quality-diversity,” International Conference on Parallel Problem Solving from Nature (2018), pp. 500–511. [CrossRef]

**45. **Steel Bank Common Lisp (SBCL) homepage, (2019), http://www.sbcl.org/.