Optimization

We are currently in the early stages of implementing optimization for lens surfaces.

We have had success using Optim.jl and ForwardDiff.jl among other packages, though ForwardDiff.jl gradient and hessian pre-compilation can be very slow for complex systems.

Other optimization packages should integrate easily with the system too: JuMP.jl, Ipopt.jl and NLOpt.jl are some options.

Merit functions must currently be implemented by the user, this is quite straight-forward as trace returns the vast majority of information that could be needed in the form of a LensTrace which hits the detector (or nothing if the ray doesn't reach the detector).

There are some helper functions implemented for AxisymmetricOpticalSystems which can make optimization of basic systems much easier:

OpticSim.Optimization — Module

Optimization interface consists of two functions optimizationvariables and updateoptimizationvariables. optimizationvariables packs variables to be optimized into a vector. updateoptimizationvariables receives a vector of variables and creates a new optical system with the variable values.

source

OpticSim.Optimization.optimizationvariables — Function

optimizationvariables(a::AxisymmetricOpticalSystem{T}) -> Vector{T}

Pack variables that have been marked to be optimized into a vector in a form suitable for the optimizer. Variables are marked for optimization by having a true value in the :OptimizeName column, where Name can be Radius, Thickness or Conic.

source

OpticSim.Optimization.updateoptimizationvariables — Function

updateoptimizationvariables(a::AxisymmetricOpticalSystem{T}, optimizationvariables::Vector{S}) -> AxisymmetricOpticalSystem{S}

Creates a new optical system with updated variables corresponding to the optimization variables.

source

It is of course possible to write your own optimization loop for more complex (i.e. non-AxisymmetricOpticalSystem) systems and this should work without issue.

Example

using OpticSim
using ForwardDiff
using Optim

function objective(a::AbstractVector{T}, b::AxisymmetricOpticalSystem{T}, samples::Int = 3) where {T}
    # RMSE spot size
    system = Optimization.updateoptimizationvariables(b, a)
    # distribute rays evenly across entrance pupil using HexapolarField. The latest version of OpticSim no longer supports HexapolarField. 
    field = HexapolarField(system, collimated = true, samples = samples)
    error = zero(T)
    hits = 0
    for r in field
        traceres = trace(system, r, test = true)
        if traceres !== nothing # ignore rays which miss
            hitpoint = point(traceres)
            if abs(hitpoint[1]) > eps(T) && abs(hitpoint[2]) > eps(T)
                dist_to_axis = hitpoint[1]^2 + hitpoint[2]^2
                error += dist_to_axis
            end
            hits += 1
        end
    end
    if hits > 0
        error = sqrt(error / hits)
    end
    # if hits == 0 returns 0 - not ideal!
    return error
end

start, lower, upper = Optimization.optimizationvariables(system)
optimobjective = arg -> objective(arg, system)
gcfg = ForwardDiff.GradientConfig(optimobjective, start, ForwardDiff.Chunk{1}()) # speed up ForwardDiff significantly
hcfg = ForwardDiff.HessianConfig(optimobjective, start, ForwardDiff.Chunk{1}())
g! = (G, x) -> ForwardDiff.gradient!(G, optimobjective, x, gcfg)
h! = (H, x) -> ForwardDiff.hessian!(H, optimobjective, x, hcfg)
df = TwiceDifferentiable(optimobjective, g!, h!, start)
dfc = TwiceDifferentiableConstraints(lower, upper) # constrain the optimization to avoid e.g. thickness < 0
res = optimize(df, dfc, start, algo, Optim.Options(show_trace = true, iterations = 100, allow_f_increases = true))
final = Optim.minimizer(res)
new_system = Optimization.updateoptimizationvariables(system, final)