Custom Optimization Algorithms

Our goal will be to optimize the coefficients of an aspheric singlet to minimize the RMS spot size. We start by first defining this singlet with all aspheric coefficients set to zero.

class AsphericSinglet(optic.Optic):
 """Aspheric singlet"""

 def __init__(self):
     super().__init__()

     # add surfaces
     self.add_surface(index=0, radius=np.inf, thickness=np.inf)
     self.add_surface(
         index=1,
         thickness=7,
         radius=20.0,
         is_stop=True,
         material="N-SF11",
         surface_type="even_asphere",
         conic=0.0,
         coefficients=[0, 0, 0],
     )
     self.add_surface(index=2, thickness=21.56201105)
     self.add_surface(index=3)

     # add aperture
     self.set_aperture(aperture_type="EPD", value=20.0)

     # add field
     self.set_field_type(field_type="angle")
     self.add_field(y=0)

     # add wavelength
     self.add_wavelength(value=0.55, is_primary=True)

Next, we need to define our optimization algorithm. This is a class that inherits from optimization.OptimizerGeneric. This class has two requirements:

• The constructor must accept the optimization problem as an argument. This is an instance of optimization.OptimizationProblem.
• The class must implement the optimize method.

The optimization algorithm is as follows:

• Get current variable values ("position") and the objective function value.

• Save the initial position in case we want to undo the optimization later (optional, but good practice).

• Generate a random step.

• Calculate the objective function at the new position.

• If the new position is better, accept the step.

• Continue for the maximum number of steps.

• Update the variables to the optimal position.

  pythonCopy

  class RandomWalkOptimizer(optimization.OptimizerGeneric):
   def __init__(self, problem):
       super().__init__(problem)
  
   def optimize(self, max_steps=100, delta=0.1, seed=42):
       # Set random seed
       np.random.seed(seed)
  
       # Get current position and objective function value
       current_position = [var.value for var in self.problem.variables]
       current_value = self._fun(current_position)
       num_variables = len(current_position)
  
       # save initial position to be able to revert
       self._x.append(current_position)
  
       # Save values of each iteration
       values = [current_value]
  
       for _ in range(max_steps):
           # Generate a random step
           random_step = np.random.randn(num_variables) * delta
           new_position = current_position + random_step
  
           # Calculate the objective function value at the new position
           new_value = self._fun(new_position)
           values.append(new_value)
  
           # If the new value is better, accept the step
           if new_value < current_value:
               current_position = new_position
               current_value = new_value
  
       # Update the variables with the best position
       for idvar, var in enumerate(self.problem.variables):
           var.update(current_position[idvar])
  
       return values

Let's create and draw the starting point lens.

lens = AsphericSinglet()
lens.draw(num_rays=5)

We now create the optimization problem, add our operand (RMS spot size target) and variables (3 aspheric coefficients).

problem = optimization.OptimizationProblem()

We print the optimization problem info to see the current status.

problem.info()

Let's now create an optimizer using our RandomWalkOptimizer class. We pass the optimization problem as an argument.

optimizer = RandomWalkOptimizer(problem)

Finally, we can run the optimization by calling the optimize. We choose to run the optimization for 1000 steps and use a delta value of 0.1. The delta value controls the size of the step in each iteration.

values = optimizer.optimize(max_steps=1000, delta=0.1)

On the author's machine, the optimization completes in 2 seconds. Let's print the problem information to see whether optimization was successful.

problem.info()

Indeed, the objective function value is minimized. We also draw the lens and generate the spot diagram to confirm.

lens.draw(num_rays=5)

Recall that we also saved the values of the objective function in each iteration. We can now use these values to plot the convergence of the objective function over the 1000 iterations.

plt.plot(values)
plt.xlabel("Iteration")
plt.ylabel("Objective Function Value")
plt.title("Convergence of Random Walk Optimizer")
plt.grid(True)
plt.show()

We can clearly see that the optimization converges within about 250 iterations. We can therefore conclude that the random walk optimizer is effective for this problem. For improved performance, we might consider adding a stopping condition in the future to avoid unnecessary iterations.