# 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.OptimizationModule

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.optimizationvariablesFunction
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

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)