Quickstart Guide

This guide shows you how to set up and solve a quantum optimal control problem in Piccolo.jl. We'll synthesize a single-qubit X gate.

The Problem

We want to find control pulses that implement an X gate on a single qubit with system Hamiltonian:

\[H(t) = \frac{\omega}{2} \sigma_z + u_1(t) \sigma_x + u_2(t) \sigma_y\]

using Piccolo

Step 1: Define the Quantum System

First, we define our quantum system by specifying the drift Hamiltonian (always-on), the drive Hamiltonians (controllable), and the bounds on control amplitudes.

# Drift Hamiltonian: qubit frequency term
H_drift = 0.5 * PAULIS[:Z]

# Drive Hamiltonians: X and Y controls
H_drives = [PAULIS[:X], PAULIS[:Y]]

# Maximum control amplitudes
drive_bounds = [1.0, 1.0]

# Create the quantum system
sys = QuantumSystem(H_drift, H_drives, drive_bounds)

QuantumSystem: levels = 2, n_drives = 2

Step 2: Create an Initial Pulse

We need an initial guess for the control pulse. ZeroOrderPulse represents piecewise constant controls.

# Time parameters
T = 10.0   # Total gate duration
N = 100    # Number of timesteps

# Create time vector
times = collect(range(0, T, length = N))

# Random initial controls (scaled by drive bounds)
initial_controls = 0.1 * randn(2, N)

# Create the pulse
pulse = ZeroOrderPulse(initial_controls, times)

ZeroOrderPulse
  drives: 2
  duration: 10.0

Step 3: Define the Goal via a Trajectory

A UnitaryTrajectory combines the system, pulse, and target gate.

# Target: X gate
U_goal = GATES[:X]

# Create the trajectory
qtraj = UnitaryTrajectory(sys, pulse, U_goal)

UnitaryTrajectory{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.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}(QuantumSystem: levels = 2, n_drives = 2, ZeroOrderPulse(Number of drives = 2, T = 10.0), ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im], ComplexF64[0.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im], ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im;;; 0.9985507115505804 - 0.049975842857885096im 0.002900417523501793 + 0.01976054580883007im; -0.002900417523501794 + 0.019760545808830077im 0.9985507115505803 + 0.04997584285788509im;;; 0.9947973767169338 - 0.09971897507199891im 0.008812316409221893 + 0.018885136144634548im; -0.008812316409221897 + 0.018885136144634548im 0.9947973767169335 + 0.0997189750719989im;;; … ;;; 0.16196027596370674 + 0.967208298726967im 0.11033988836472247 - 0.1615613967478452im; -0.11033988836472253 - 0.16156139674784506im 0.1619602759637095 - 0.9672082987269693im;;; 0.21006722280917137 + 0.9602163727160283im 0.09378462754264098 - 0.15830578986888094im; -0.09378462754264116 - 0.15830578986888086im 0.21006722280917436 - 0.9602163727160311im;;; 0.25851090814435507 + 0.9487763430584681im 0.08758421906895261 - 0.15913694037862183im; -0.08758421906895278 - 0.15913694037862164im 0.2585109081443581 - 0.9487763430584705im])

Step 4: Set Up the Optimization Problem

SmoothPulseProblem creates the optimization problem with:

Fidelity objective (weight Q)
Regularization for smooth controls (weight R)
Derivative bounds for control smoothness

qcp = SmoothPulseProblem(
    qtraj,
    N;
    Q = 100.0,       # Fidelity weight
    R = 1e-2,        # Regularization weight
    ddu_bound = 1.0,  # Control acceleration bound
)

QuantumControlProblem{UnitaryTrajectory{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.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}}
  System: QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}
  Goal: Matrix{ComplexF64}
  Trajectory: 100 knots
  State: Ũ⃗
  Controls: u

Step 5: Solve!

solve!(qcp; max_iter = 20, verbose = false, print_level = 1)

Step 6: Analyze Results

After solving, we can check the fidelity and examine the optimized controls.

# Check final fidelity
fidelity(qcp)

0.9999995867273045

Access the trajectory and check the final unitary:

traj = get_trajectory(qcp)
U_final = iso_vec_to_operator(traj[:Ũ⃗][:, end])
round.(U_final, digits = 3)

2×2 Matrix{ComplexF64}:
    0.0+0.0im  0.001-1.0im
 -0.001-1.0im    0.0-0.0im

Visualization

Piccolo provides specialized plotting functions for quantum trajectories:

using CairoMakie

# Plot the unitary evolution (state populations over time)
fig = plot_unitary_populations(traj)

Saving Your Results

Save the optimized pulse so you can reload it later without re-solving. save writes the pulse under the JLD2 key "pulse"; load_pulse reads it back as the original pulse type:

optimized_pulse = get_pulse(qcp.qtraj)
save("quickstart_pulse.jld2", optimized_pulse)

Reload in any script with:

saved_pulse = load_pulse("quickstart_pulse.jld2")
qtraj = UnitaryTrajectory(sys, saved_pulse, GATES[:X])

To bundle a pulse with metadata (fidelity, gate name, system config), use jldsave with keyword arguments:

using JLD2
jldsave("quickstart_pulse.jld2"; pulse=optimized_pulse, fidelity=fidelity(qcp))

See the Saving and Loading Pulses guide for warm-starting and other patterns.

Minimum Time Optimization

Now let's find the shortest gate duration that achieves 99% fidelity.

First, we need to create a new problem with variable timesteps enabled:

# Create problem with free-time optimization
qcp_free = SmoothPulseProblem(
    qtraj,
    N;
    Q = 100.0,
    R = 1e-2,
    ddu_bound = 1.0,
    Δt_bounds = (0.01, 0.5),  # Enable variable timesteps
)
solve!(
    qcp_free;
    max_iter = 20,
    verbose = false,
    print_level = 1,
)

# Convert to minimum time problem
qcp_mintime = MinimumTimeProblem(qcp_free; final_fidelity = 0.99)
solve!(
    qcp_mintime;
    max_iter = 20,
    verbose = false,
    print_level = 1,
)

    constructing SmoothPulseProblem for UnitaryTrajectory...
    constructing MinimumTimeProblem from QuantumControlProblem{UnitaryTrajectory{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.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Base.Pairs{Symbol, Vector{Float64}, Nothing, @NamedTuple{tstops::Vector{Float64}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqLinear.MagnusAdapt4, OrdinaryDiffEqCore.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.AutoSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#53#54"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#55#56"{Vector{SparseArrays.SparseMatrixCSC{Float64, Int64}}, Int64, SparseArrays.SparseMatrixCSC{Float64, Int64}}, @NamedTuple{}, Vector{DriftTerm{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Returns{Float64}, Returns{Float64}}}, SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}}, Val{()}}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Piccolo.Quantum.Rollouts.PiccoloRolloutSystem{Union{Int64, AbstractVector{Int64}, CartesianIndex, CartesianIndices}}, Nothing, Nothing}, Vector{Matrix{ComplexF64}}, Vector{Float64}, Vector{Vector{Matrix{ComplexF64}}}, Nothing, OrdinaryDiffEqLinear.MagnusAdapt4Cache{Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Matrix{ComplexF64}, Nothing}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}}}...
	final fidelity constraint: 0.99
	minimum-time weight D: 100.0

Compare durations:

initial_duration = sum(get_timesteps(get_trajectory(qcp_free)))
minimum_duration = sum(get_timesteps(get_trajectory(qcp_mintime)))

initial_duration

12.429085117242698

minimum_duration

5.605073009405118

fidelity(qcp_mintime)

0.9908208074426424

Plot the time-optimal solution:

fig_mintime = plot_unitary_populations(get_trajectory(qcp_mintime))

State Preparation

Instead of synthesizing a full unitary gate, you can prepare a specific quantum state using KetTrajectory:

ψ_init = ComplexF64[1.0, 0.0]  # |0⟩
ψ_goal = ComplexF64[0.0, 1.0]  # |1⟩
qcp_state = SmoothPulseProblem(KetTrajectory(sys, pulse, ψ_init, ψ_goal), N)
solve!(
    qcp_state;
    max_iter = 20,
    verbose = false,
    print_level = 1,
)
fidelity(qcp_state)

0.994225519370145

Robust Control

To optimize a pulse that works across parameter variations (e.g., uncertain qubit frequency), use SamplingProblem:

# Perturbed systems: ±10% drift Hamiltonian
sys_low = QuantumSystem(0.9 * H_drift, H_drives, drive_bounds)
sys_high = QuantumSystem(1.1 * H_drift, H_drives, drive_bounds)

# Start from a nominal solution, then add robustness
qcp_robust = SamplingProblem(qcp, [sys_low, sys, sys_high])
solve!(
    qcp_robust;
    max_iter = 20,
    verbose = false,
    print_level = 1,
)
fidelity(qcp_robust)

3-element Vector{Float64}:
 0.9970676593183584
 0.9999785687241922
 0.9971921002966301

Next Steps

See Concepts for the mathematical formulation
Learn about different Problem Templates
Read Saving and Loading Pulses to organize and reuse your results
Explore Tutorials for more complex examples

This page was generated using Literate.jl.