Example - Optimization Using Reciprocal Radii

April 2025

This example demonstrates smooth surface transitions during optimization using the “reciprocal_radius” variable. This approach helps design optical systems without initially knowing whether each surface should be concave or convex.

The “reciprocal_radius” is defined as 1/radius (the inverse of the radius of curvature). You define initial conditions using the regular “radius” parameter, but the optimization operates on this proxy variable instead.

While the libary code attempts to handle edge cases, there are limitations. It’s best to avoid setting the initial radius to exactly zero (which is physically meaningless) or infinity. If you do create a system with a radius of ±infinity or if you manually set reciprocal_radius to zero, it is possible that the optimization stops prematurely. In this case, just re-run the optimization from where it stopped, using either the same algorithm or a different one.

import logging

import numpy as np

from optiland import analysis, optic, optimization

# Disable the warning "findfont: Font family 'cambria' not found." when running on Linux
logging.getLogger("matplotlib.font_manager").disabled = True

## 1. Definition of initial system

We’ll define a simple BK7 plano-concave singlet and optimize it to focus on the focal plane.

We’ll try two optimization approaches: - First using the standard “radius” variable - Then using the “reciprocal_radius” variable

class Singlet(optic.Optic):
    def __init__(self):
        super().__init__()

        # Define surfaces
        # object plane:
        self.surfaces.add(index=0, thickness=np.inf)
        # first surface, initially concave:
        self.surfaces.add(
            index=1, radius=-100, thickness=5, material="BK7", is_stop=True
        )
        # the second surface is flat:
        self.surfaces.add(index=2, radius=np.inf, thickness=45)
        # image plane:
        self.surfaces.add(index=3)

        # Define aperture
        self.set_aperture(aperture_type="EPD", value=10.0)

        # Define fields
        self.fields.set_type(field_type="angle")
        self.fields.add(y=0)
        self.fields.add(y=1)
        self.fields.add(y=2)

        # Define wavelengths
        # Uncomment the blue and red wavelengths to see chromatic aberration effects
        # self.wavelengths.add(value=0.4861)
        self.wavelengths.add(value=0.5876, is_primary=True)
        # self.wavelengths.add(value=0.6563)

Here’s our starting lens configuration:

lens = Singlet()
lens.draw()

(<Figure size 1000x400 with 1 Axes>, <Axes: xlabel='Z [mm]', ylabel='Y [mm]'>)

2. Optimization

Now we’ll define the optimization problem to achieve focus on a specific plane and minimize the RMS spot size.

problem = optimization.OptimizationProblem()

# Add requirement for spot size
for field in lens.fields.get_field_coords():
    input_data = {
        "optic": lens,
        "Hx": field[0],
        "Hy": field[1],
        "num_rays": 5,
        "wavelength": "all",
        "distribution": "hexapolar",
        "surface_number": 3,
    }
    problem.add_operand(
        operand_type="rms_spot_size", target=0, weight=1, input_data=input_data
    )

Approach 1: Optimizing with “radius”

First, let’s try optimizing using the standard “radius” variable. We don’t set limits here because we can’t easily specify that radius should be outside an interval (e.g., radius ≥ 10 or radius ≤ -10).

problem.add_variable(lens, "radius", surface_number=1)

Let’s check the initial problem configuration:

problem.info()

╒════╤════════════════════════╤═══════════════════╕
│    │   Merit Function Value │   Improvement (%) │
╞════╪════════════════════════╪═══════════════════╡
│  0 │                69.5755 │                 0 │
╘════╧════════════════════════╧═══════════════════╛
╒════╤════════════════╤══════════╤══════════════╤══════════════╤══════════╤═══════════════╤═════════╤═════════╤════════════════╕
│    │ Operand Type   │   Target │ Min. Bound   │ Max. Bound   │   Weight │   Eff. Weight │   Value │   Delta │   Contrib. [%] │
╞════╪════════════════╪══════════╪══════════════╪══════════════╪══════════╪═══════════════╪═════════╪═════════╪════════════════╡
│  0 │ rms spot size  │        0 │              │              │        1 │             1 │   4.815 │   4.815 │          33.32 │
│  1 │ rms spot size  │        0 │              │              │        1 │             1 │   4.816 │   4.816 │          33.33 │
│  2 │ rms spot size  │        0 │              │              │        1 │             1 │   4.817 │   4.817 │          33.35 │
╘════╧════════════════╧══════════╧══════════════╧══════════════╧══════════╧═══════════════╧═════════╧═════════╧════════════════╛
╒════╤═════════════════╤═══════════╤═════════╤══════════════╤══════════════╕
│    │ Variable Type   │   Surface │   Value │ Min. Bound   │ Max. Bound   │
╞════╪═════════════════╪═══════════╪═════════╪══════════════╪══════════════╡
│  0 │ radius          │         1 │    -100 │              │              │
╘════╧═════════════════╧═══════════╧═════════╧══════════════╧══════════════╛

Now we’ll run the optimization:

optimizer = optimization.OptimizerGeneric(problem)
res = optimizer.optimize(tol=1e-9)

Let’s examine the optimization results:

problem.info()
lens.draw()

╒════╤════════════════════════╤═══════════════════╕
│    │   Merit Function Value │   Improvement (%) │
╞════╪════════════════════════╪═══════════════════╡
│  0 │                44.5212 │           36.0102 │
╘════╧════════════════════════╧═══════════════════╛
╒════╤════════════════╤══════════╤══════════════╤══════════════╤══════════╤═══════════════╤═════════╤═════════╤════════════════╕
│    │ Operand Type   │   Target │ Min. Bound   │ Max. Bound   │   Weight │   Eff. Weight │   Value │   Delta │   Contrib. [%] │
╞════╪════════════════╪══════════╪══════════════╪══════════════╪══════════╪═══════════════╪═════════╪═════════╪════════════════╡
│  0 │ rms spot size  │        0 │              │              │        1 │             1 │   3.852 │   3.852 │          33.33 │
│  1 │ rms spot size  │        0 │              │              │        1 │             1 │   3.852 │   3.852 │          33.33 │
│  2 │ rms spot size  │        0 │              │              │        1 │             1 │   3.852 │   3.852 │          33.33 │
╘════╧════════════════╧══════════╧══════════════╧══════════════╧══════════╧═══════════════╧═════════╧═════════╧════════════════╛
╒════╤═════════════════╤═══════════╤═════════╤══════════════╤══════════════╕
│    │ Variable Type   │   Surface │   Value │ Min. Bound   │ Max. Bound   │
╞════╪═════════════════╪═══════════╪═════════╪══════════════╪══════════════╡
│  0 │ radius          │         1 │ -141387 │              │              │
╘════╧═════════════════╧═══════════╧═════════╧══════════════╧══════════════╛

(<Figure size 1000x400 with 1 Axes>, <Axes: xlabel='Z [mm]', ylabel='Y [mm]'>)

The optimization struggles to converge to a good solution. This happens because there’s no continuous path from a concave to a convex surface when using the “radius” directly.

Approach 2: Optimizing with “reciprocal_radius”

lens = Singlet()

problem = optimization.OptimizationProblem()

# Add requirement for spot size
for field in lens.fields.get_field_coords():
    input_data = {
        "optic": lens,
        "Hx": field[0],
        "Hy": field[1],
        "num_rays": 5,
        "wavelength": "all",
        "distribution": "hexapolar",
        "surface_number": 3,
    }
    problem.add_operand(
        operand_type="rms_spot_size", target=0, weight=1, input_data=input_data
    )

# Use reciprocal_radius instead of radius for optimization
problem.add_variable(lens, "reciprocal_radius", surface_number=1)

Let’s run the optimization:

optimizer = optimization.OptimizerGeneric(problem)
res = optimizer.optimize(tol=1e-9)

And check the results:

problem.info()
lens.draw()

╒════╤════════════════════════╤═══════════════════╕
│    │   Merit Function Value │   Improvement (%) │
╞════╪════════════════════════╪═══════════════════╡
│  0 │            0.000455027 │           99.9993 │
╘════╧════════════════════════╧═══════════════════╛
╒════╤════════════════╤══════════╤══════════════╤══════════════╤══════════╤═══════════════╤═════════╤═════════╤════════════════╕
│    │ Operand Type   │   Target │ Min. Bound   │ Max. Bound   │   Weight │   Eff. Weight │   Value │   Delta │   Contrib. [%] │
╞════╪════════════════╪══════════╪══════════════╪══════════════╪══════════╪═══════════════╪═════════╪═════════╪════════════════╡
│  0 │ rms spot size  │        0 │              │              │        1 │             1 │   0.012 │   0.012 │          32.21 │
│  1 │ rms spot size  │        0 │              │              │        1 │             1 │   0.012 │   0.012 │          31.48 │
│  2 │ rms spot size  │        0 │              │              │        1 │             1 │   0.013 │   0.013 │          36.31 │
╘════╧════════════════╧══════════╧══════════════╧══════════════╧══════════╧═══════════════╧═════════╧═════════╧════════════════╛
╒════╤═══════════════════╤═══════════╤═══════════╤══════════════╤══════════════╕
│    │ Variable Type     │   Surface │     Value │ Min. Bound   │ Max. Bound   │
╞════╪═══════════════════╪═══════════╪═══════════╪══════════════╪══════════════╡
│  0 │ reciprocal_radius │         1 │ 0.0396763 │              │              │
╘════╧═══════════════════╧═══════════╧═══════════╧══════════════╧══════════════╛

(<Figure size 1000x400 with 1 Axes>, <Axes: xlabel='Z [mm]', ylabel='Y [mm]'>)

Conclusions

• When using the “radius” variable directly, the optimizer often struggles find solutions that require changing from concave to convex or vice versa.

• The “reciprocal_radius” variable allows for smooth transitions between concave and convex surfaces, producing better optimization results in this specific example.

• Initial system conditions are still specified with “radius” in the lens definition, but the optimization works on “reciprocal_radius” internally.

• For fine tuning, after surfaces have reached the correct concavity/convexity, there is no problem in switching back to “radius” for the final optimization steps.

• If the optimization stops prematurely due to edge cases (e.g., initial radius equal to ±infinity), simply re-run the optimization from where it stopped, using either the same algorithm or a different one.