GitHub

Working with Solutions

This guide covers how to solve problems, extract data from solutions, evaluate performance, and save/load results.

using QuantumCollocation
using PiccoloQuantumObjects
using NamedTrajectories

Solving Problems

Once you've created a problem template, solving is straightforward with the solve! function:

system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y], [1.0, 1.0])
U_goal = EmbeddedOperator(GATES.H, system)
T = 10.0 # time duration
qtraj = UnitaryTrajectory(system, U_goal, T)

N = 51 # number of timesteps
prob = SmoothPulseProblem(qtraj, N)
QuantumControlProblem{PiccoloQuantumObjects.QuantumTrajectories.UnitaryTrajectory{PiccoloQuantumObjects.Pulses.ZeroOrderPulse{DataInterpolations.ConstantInterpolation{Matrix{Float64}, Vector{Float64}, Vector{Union{}}, Float64}}, SciMLBase.ODESolution{ComplexF64, 3, Vector{Matrix{ComplexF64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, SciMLBase.ODEProblem{Matrix{ComplexF64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{PiccoloQuantumObjects.Rollouts.var"#update!#_construct_operator##2"{PiccoloQuantumObjects.QuantumSystems.QuantumSystem{PiccoloQuantumObjects.QuantumSystems.var"#26#27"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, PiccoloQuantumObjects.QuantumSystems.var"#28#29"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, PiccoloQuantumObjects.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}, saveat::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusGL4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{PiccoloQuantumObjects.Rollouts.var"#update!#_construct_operator##2"{PiccoloQuantumObjects.QuantumSystems.QuantumSystem{PiccoloQuantumObjects.QuantumSystems.var"#26#27"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, PiccoloQuantumObjects.QuantumSystems.var"#28#29"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, PiccoloQuantumObjects.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusGL4Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, PiccoloQuantumObjects.EmbeddedOperators.EmbeddedOperator{ComplexF64}}}
  System: PiccoloQuantumObjects.QuantumSystems.QuantumSystem{PiccoloQuantumObjects.QuantumSystems.var"#26#27"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, PiccoloQuantumObjects.QuantumSystems.var"#28#29"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}}
  Goal: PiccoloQuantumObjects.EmbeddedOperators.EmbeddedOperator{ComplexF64}
  Trajectory: 51 knots
  State: Ũ⃗
  Controls: u

The solve! function accepts several key options:

solve!(prob;
    max_iter = 100,        # Maximum optimization iterations
    verbose = true,         # Print convergence information
    print_level = 1         # Ipopt output level (0-5)
)
    initializing optimizer...
        applying constraint: timesteps all equal constraint
        applying constraint: initial value of Ũ⃗
        applying constraint: initial value of u
        applying constraint: final value of u
        applying constraint: bounds on Ũ⃗
        applying constraint: bounds on u
        applying constraint: bounds on du
        applying constraint: bounds on ddu
        applying constraint: bounds on Δt
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

Understanding Convergence

Ipopt reports several key metrics:

  • Objective: Current cost function value
  • inf_pr: Constraint violation (primal infidelity) - should go to ~0
  • inf_du: Dual infidelity - measure of optimality
  • lg(mu): Log of barrier parameter
  • alphadu/alphapr: Step sizes

Successful convergence typically shows inf_pr < 1e-6 and status Optimal Solution Found.

Extracting Data from Solutions

Accessing Controls

Control pulses are stored in the trajectory with automatic naming:

u = prob.trajectory.u       # Control values [n_drives × N]
du = prob.trajectory.du     # First derivatives
ddu = prob.trajectory.ddu   # Second derivatives

println("Control shape: ", size(u))
println("Number of drives: ", size(u, 1))
println("Number of timesteps: ", size(u, 2))
Control shape: (2, 51)
Number of drives: 2
Number of timesteps: 51

Access individual drive controls:

u_drive_1 = u[1, :]   # First drive over time
51-element Vector{Float64}:
 0.0
 0.09102972396310519
 0.1683407407610345
 0.22596245578236104
 0.26220656442967155
 0.2778324893716056
 0.27519111197603585
 0.25758291699716895
 0.22860014092405484
 0.19166799775224544
 ⋮
 0.20730834935883788
 0.20598862060908002
 0.19709453853529052
 0.18042351020245914
 0.1560649265242895
 0.12450022846037355
 0.0867252119063542
 0.04439604362165721
 0.0

Accessing States

For unitary problems:

Ũ⃗ = prob.trajectory.Ũ⃗  # Vectorized unitary [N² × T]
8×51 view(reshape(view(::Vector{Float64}, :), 16, 51), 1:8, :) with eltype Float64:
 1.0   0.999779      0.998895    0.996209   …   0.0214745    1.65558e-5
 0.0  -4.23077e-32  -0.0092729  -0.0270897     -0.00825884   4.5671e-6
 0.0  -0.0210169    -0.0420183  -0.0629547     -0.709828    -0.707091
 0.0  -3.48142e-33  -0.0189377  -0.0537388     -0.703971    -0.707095
 0.0   9.17594e-33   0.0092729   0.0270897      0.00825884  -4.5671e-6
 1.0   0.999779      0.998895    0.996209   …   0.0214745    1.65558e-5
 0.0   5.63116e-34  -0.0189377  -0.0537388     -0.703971    -0.707095
 0.0   0.0210169     0.0420183   0.0629547      0.709828     0.707091

The unitary is stored in "isovec" format (vectorized). To get the actual unitary matrix at timestep k:

using LinearAlgebra
k = N  # Final timestep
U_k = iso_vec_to_operator(Ũ⃗[:, k])
println("Final unitary dimensions: ", size(U_k))
Final unitary dimensions: (2, 2)

For ket (state) problems, use: ψ̃ = prob.trajectory.ψ̃ # Vectorized state [2N × T]

Time Grid

Access timestep information:

Δt_vec = prob.trajectory.Δt  # Timestep durations
1×51 view(reshape(view(::Vector{Float64}, :), 16, 51), 9:9, :) with eltype Float64:
 0.210184  0.210184  0.210184  0.210184  …  0.210184  0.210184  0.210184

Calculate total duration:

duration = get_duration(prob.trajectory)
println("Total gate time: ", duration, " (arbitrary units)")
Total gate time: 10.509209626401601 (arbitrary units)

For minimum time problems, timesteps vary: minprob = UnitaryMinimumTimeProblem(prob, Ugoal) solve!(minprob, maxiter=100) Δtoptimized = minprob.trajectory.Δt # Variable timesteps

Evaluating Solutions

Fidelity Calculations

Direct fidelity - Compare final state to goal:

U_final = iso_vec_to_operator(prob.trajectory.Ũ⃗[:, end])
fid_direct = unitary_fidelity(U_final, U_goal.operator)
println("Direct fidelity: ", fid_direct)
Direct fidelity: 0.9999607437739154

Rollout fidelity - Simulate dynamics forward:

fid_rollout = unitary_rollout_fidelity(prob.trajectory, system)
println("Rollout fidelity: ", fid_rollout)
Rollout fidelity: 0.9997937552347185

The rollout fidelity is more accurate as it accounts for actual dynamics, while direct fidelity only checks the final point.

For Embedded Operators (Multilevel Systems)

When working with subspaces (e.g., qubit in transmon): op = EmbeddedOperator(:X, system) probembedded = UnitarySmoothPulseProblem(system, op, N) solve!(probembedded, max_iter=100)

Evaluate fidelity only in computational subspace

fidsubspace = unitaryrolloutfidelity( probembedded.trajectory, system; subspace = op.subspace )

Leakage Evaluation

For multilevel systems, check population in leakage levels: This requires analyzing state populations during evolution See the Multilevel Transmon example for details

Constraint Violations

Check if solution satisfies all constraints:

  • Dynamics constraints: Compare rollout vs direct fidelity
  • Bound constraints: Verify controls within system.drive_bounds
  • Derivative constraints: Check |du|, |ddu| within bounds
println("Max control value: ", maximum(abs.(u)))
println("Max control derivative: ", maximum(abs.(du)))
Max control value: 0.2778324893716056
Max control derivative: 0.43312082255200207

Saving and Loading

Save a Trajectory

using JLD2 saveobject("mysolution.jld2", prob.trajectory)

Load and Reuse

Load trajectory for warm-starting: trajloaded = loadobject("my_solution.jld2")

Use as initial guess for new problem

probrefined = UnitarySmoothPulseProblem( system, Ugoal, T, Δt; uguess = trajloaded.u, piccolo_options = PiccoloOptions(verbose=false) )

Save Entire Problem

To save all problem information including constraints and objectives: saveobject("myproblem.jld2", prob)

Post-Processing

Resampling Trajectories

Change time discretization while preserving control shape: Tnew = 101 # More timesteps trajresampled = resample(prob.trajectory, T_new)

Extracting for Experiments

Prepare control pulses for hardware:

using PiccoloPlots  # For visualization
using CairoMakie

Plot controls

fig = plot(prob.trajectory)
Example block output

save("controls.png", fig)

Extract control data for export

control_data = Dict(
    "time" => cumsum([0; prob.trajectory.Δt[:]]),
    "drive_1_real" => u[1, :],
    "drive_2_real" => u[2, :],
    "duration" => duration
)
Dict{String, Any} with 4 entries:
  "duration"     => 10.5092
  "time"         => [0.0, 0.210184, 0.420368, 0.630553, 0.840737, 1.05092, 1.26…
  "drive_1_real" => [0.0, 0.0910297, 0.168341, 0.225962, 0.262207, 0.277832, 0.…
  "drive_2_real" => [0.0, -0.0422436, -0.0797889, -0.110932, -0.135205, -0.1528…

Save to CSV or other format for AWG

using CSV, DataFrames df = DataFrame(controldata) CSV.write("pulsesequence.csv", df)

Pulse Filtering

Apply smoothing to reduce high-frequency noise: using DSP for i in 1:size(u, 1) u_filtered[i, :] = filtfilt(responsetype, u[i, :]) end

Best Practices

Starting a new optimization:

  1. Begin with coarse discretization (small T)
  2. Use relaxed bounds and convergence criteria
  3. Refine solution incrementally
  4. Use previous solution as warm start

Debugging poor convergence:

  1. Check inf_pr - high values indicate constraint violations
  2. Verify system Hamiltonian is correct
  3. Try looser derivative bounds (dubound, ddubound)
  4. Increase regularization weights (R_u, R_du, R_ddu)
  5. Use piccolo_options.bound_state=true for better numerics

Improving solutions:

  1. Increase T (more timesteps = finer control)
  2. Add derivative constraints for smoother pulses
  3. Use minimum time optimization for fastest gates
  4. Apply leakage constraints for multilevel systems
  5. Use sampling problems for robust control

This page was generated using Literate.jl.