Constant amplitude optimization in GRAPE

[mabuni1998]

Is it possible to have a GRAPE optimization where one of the controls being optimized is just a flat constant while having other controls that are time-dependent?

For example, I could define the following constant function structure:

using ComponentArrays
using QuantumPropagators.Controls: ParameterizedFunction, get_parameters
#Constant function definition ----------------------------------------------------
struct ConstantFunction <: ParameterizedFunction
    parameters::ComponentVector{Float64,Vector{Float64},Tuple{Axis{(A=1,)}}}
end
function ConstantFunction(; A=1.0)
    return ConstantFunction(ComponentVector(; A))
end
function (control::ConstantFunction)(t)
    @unpack A = control.parameters
    return A
end

Then if I do:

using QuantumOptics
using QuantumControl.Amplitudes: ShapedAmplitude
using QuantumControl.Functionals: J_T_sm
using QuantumPropagators: ExpProp

b = FockBasis(1)
sm = destroy(b)
x = (sm + dagger(sm))
sigma_z = dagger(sm) * sm - sm * dagger(sm)

tlist = collect(range(0, 5, step=0.05))
control_field(t) = QuantumControl.Shapes.flattop(t, T=5, t_rise=0.3, func=:blackman);
constant_ampl = ConstantFunction(A=5)

H = hamiltonian((sigma_z.data, constant_ampl), (x.data, control_field))

psi_init = (fockstate(b, 0).data)
psi_target = fockstate(b, 1).data
trajectories = [Trajectory(initial_state=psi_init, generator=H, target_state=psi_target)];

problem = ControlProblem(
    trajectories=trajectories,
    tlist=tlist,
    iter_stop=500,
    prop_method=ExpProp,
    pulse_options=Dict(),
    J_T=J_T_sm,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end,
);

propagator = QuantumControl.init_prop_trajectory(trajectories[1], tlist; psi_init,method=ExpProp)

using GRAPE
datadir(names...) = joinpath(@__DIR__, names...)
opt_result_LBFGSB = @optimize_or_load(
    datadir("TLS", "GRAPE_opt_result.jld2"),
    problem;
    method=GRAPE,
    prop_method=ExpProp,
);

And plot it...

using PyPlot
pygui(true)
fig,ax = subplots(1,1)
ax.plot(opt_result_LBFGSB.tlist, opt_result_LBFGSB.optimized_controls[1], label="Optimized control 1")
ax.plot(opt_result_LBFGSB.tlist, opt_result_LBFGSB.optimized_controls[2], label="Optimized control 2")
ax.set_xlabel("Time")
ax.set_ylabel("Control amplitude")
ax.legend()
savefig("optimized_controls.png")

Here, we see that the amplitude that should have been a constant has a time dependence. I think it is also explicitly mentioned that parameterized functions should not be used together with GRAPE, so this is perhaps expected. Is it simply not possible to combine a constant and a time-varying field in GRAPE, or should it be defined differently?

[goerz]

This does not seem like it would be easily possible. GRAPE by definition optimizes piecewise-constant pulses on a fixed time grid. The implementation in GRAPE.jl calculates the gradient vector (gradient of J w.r.t. the field value in each interval of the time grid) and then feeds that to a general solver, LBFGSB by default.

You might be able to figure out something with constraints. LBFGSB allows for box constraints, for example. If you found a solver that allows you to express the constraint that all values should be the same, that would hypothetically solve the problem. I'm not sure that there is such a solver.

You could also try something with running costs: basically, add a smoothness term with a very large weight (penalizing ϵₙ - ϵₙ₋₁ ≠ 0). That would only be approximate, of course, so I'm not sure how far that will get you.

In general, I would approach this in one of two ways:

1. Start out with treating both fields as parameterized pulses. For one field, it would just be a constant (so 1 parameter), the other one could use something more flexible (something like a CRAB parameterization, for example). Once you get a result from that, do a second optimization where you don't touch the constant field at all and only optimize the second field with GRAPE.

   The big caveat of course is that ParameterizedQuantumControl is very unfinished. It has little if any documentation, and I just haven't implemented any of the machinery for calculating gradients. So you'd be limited to gradient-free optimization (which might be okay: after all, the original CRAB has always been gradient-free, simply using Nelder-Mead)

2. Alternate between optimizing the constant field and the non-constant field: Minimize the fidelity with a univariate optimization of the constant amplitude, and then do a few iterations of GRAPE. Repeat until you converge.

Of course, approach (2) also combines with approach (1). You can just iterate over the two optimizations.

It might also be possible that one could derive GRAPE from scratch for a constant field… but if that yields a numerical scheme, you would have to implement that scheme completely by hand (GRAPE.jl wouldn't help you at all). But it might be worth thinking about.