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 PiccoloStep 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 = 2Step 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.0Step 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.0029004175235017933 + 0.01976054580883007im; -0.0029004175235017937 + 0.019760545808830077im 0.9985507115505802 + 0.04997584285788508im;;; 0.9947973767169337 - 0.09971897507199891im 0.008812316409221891 + 0.018885136144634544im; -0.008812316409221893 + 0.018885136144634548im 0.9947973767169334 + 0.09971897507199888im;;; … ;;; 0.16196027596370707 + 0.9672082987269651im 0.11033988836472255 - 0.16156139674784528im; -0.11033988836472326 - 0.16156139674784525im 0.16196027596370863 - 0.9672082987269711im;;; 0.2100672228091717 + 0.9602163727160264im 0.09378462754264102 - 0.15830578986888102im; -0.0937846275426419 - 0.1583057898688811im 0.21006722280917356 - 0.960216372716033im;;; 0.25851090814435534 + 0.948776343058466im 0.08758421906895265 - 0.15913694037862186im; -0.08758421906895344 - 0.15913694037862183im 0.2585109081443573 - 0.9487763430584721im])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 · BilinearIntegrator, DerivativeIntegrator, DerivativeIntegrator
│
├─ System
│ dim=2 drives=2
│
├─ Trajectory
│ N=100 T=10.000 Δt∈[0, Inf]
│ Ũ⃗ (8) ±[1.0, 1.0, 1.0, … (8 total)] ✓ state
│ Δt (1) [0.0, Inf] ✓ timestep
│ t (1) · state
│ u (2) ±[1.0, 1.0] ✓ control
│ du (2) · control
│ ddu (2) ±[1.0, 1.0] ✓ control
│
├─ Goal
│ U_goal (2×2)
│
├─ Objective total = 104.8 @ current x
│ KnotPointObjective w=1 97.47
│ QuadraticRegularizer(:u) w=1 1.147e-04
│ QuadraticRegularizer(:du) w=1 0.02373
│ QuadraticRegularizer(:ddu) w=1 7.331
│ NullObjective w=1 0
│
├─ Constraints 1/13 violated at x₀
│ [dyn] BilinearIntegrator ✗ (‖c‖∞ = 0.01424)
│ [dyn] DerivativeIntegrator ✓ (‖c‖∞ = 2.776e-17)
│ [dyn] DerivativeIntegrator ✓ (‖c‖∞ = 4.441e-16)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [eq] TimeConsistencyConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│
└─ Status
variables: 1600 (1200 bounded)
equality: 117617
inequality: 0
F (raw) = 0.025325
Hint: show_problem(qcp; detail=:full) for pulse plot + sparsity
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.9999995867273076Access 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.0imVisualization
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 [UnitaryTrajectory]
QuantumControlProblem
├─ UnitaryTrajectory · ZeroOrderPulse · BilinearIntegrator, DerivativeIntegrator, DerivativeIntegrator
│
├─ System
│ dim=2 drives=2
│
├─ Trajectory
│ N=100 T=11.737 Δt∈[0.01, 0.5]
│ Ũ⃗ (8) ±[1.0, 1.0, 1.0, … (8 total)] ✓ state
│ Δt (1) [0.01, 0.5] ✓ timestep
│ t (1) · state
│ u (2) ±[1.0, 1.0] ✓ control
│ du (2) · control
│ ddu (2) ±[1.0, 1.0] ✓ control
│
├─ Goal
│ U_goal (2×2)
│
├─ Objective total = 1.282e-03 @ current x
│ KnotPointObjective w=1 4.133e-05
│ QuadraticRegularizer(:u) w=1 3.177e-04
│ QuadraticRegularizer(:du) w=1 3.065e-04
│ QuadraticRegularizer(:ddu) w=1 6.160e-04
│ NullObjective w=1 0
│
├─ Constraints 1/13 violated at x₀
│ [dyn] BilinearIntegrator ✗ (‖c‖∞ = 8.163e-05)
│ [dyn] DerivativeIntegrator ✓ (‖c‖∞ = 3.469e-18)
│ [dyn] DerivativeIntegrator ✓ (‖c‖∞ = 6.939e-18)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [eq] TimeConsistencyConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│
└─ Status
variables: 1600 (1300 bounded)
equality: 117617
inequality: 0
F (raw) = 1.000000
Hint: show_problem(qcp; detail=:full) for pulse plot + sparsity
constructing MinimumTimeProblem [from UnitaryTrajectory]
final fidelity ≥ 0.99
min-time weight D = 100.0
QuantumControlProblem
├─ UnitaryTrajectory · ZeroOrderPulse · BilinearIntegrator, DerivativeIntegrator, DerivativeIntegrator
│
├─ System
│ dim=2 drives=2
│
├─ Trajectory
│ N=100 T=12.305 Δt∈[0.01, 0.5]
│ Ũ⃗ (8) ±[1.0, 1.0, 1.0, … (8 total)] ✓ state
│ Δt (1) [0.01, 0.5] ✓ timestep
│ t (1) · state
│ u (2) ±[1.0, 1.0] ✓ control
│ du (2) · control
│ ddu (2) ±[1.0, 1.0] ✓ control
│
├─ Goal
│ U_goal (2×2)
│
├─ Objective total = 1230 @ current x
│ KnotPointObjective w=1 2.088e-04
│ QuadraticRegularizer(:u) w=1 3.206e-04
│ QuadraticRegularizer(:du) w=1 2.944e-04
│ QuadraticRegularizer(:ddu) w=1 5.478e-04
│ NullObjective w=1 0
│ MinimumTimeObjective w=1 1230
│
├─ Constraints 4/14 violated at x₀
│ [dyn] BilinearIntegrator ✗ (‖c‖∞ = 2.011e-06)
│ [dyn] DerivativeIntegrator ✗ (‖c‖∞ = 3.089e-07)
│ [dyn] DerivativeIntegrator ✗ (‖c‖∞ = 1.901e-06)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [bnd] BoundsConstraint ✓
│ [eq] TimeConsistencyConstraint ✓ (no eval)
│ [eq] EqualityConstraint ✓ (no eval)
│ [ineq] Terminal fidelity bound on :Ũ⃗ ✗ (viol = 9.998e-03)
│
└─ Status
variables: 1600 (1300 bounded)
equality: 117617
inequality: 1
F (raw) = 0.999998
Hint: show_problem(qcp; detail=:full) for pulse plot + sparsityCompare durations:
initial_duration = sum(get_timesteps(get_trajectory(qcp_free)))
minimum_duration = sum(get_timesteps(get_trajectory(qcp_mintime)))
initial_duration12.429085117242698minimum_duration5.605073009405118fidelity(qcp_mintime)0.9908208074426408Plot 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.9942255193701437Robust 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.9970676593183564
0.9999785687241878
0.9971921002966296Next 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.