SamplingProblem

SamplingProblem enables robust optimization over multiple system variants with shared controls. This is essential for designing pulses that perform well despite parameter uncertainty or variation.

When to Use

Use SamplingProblem when:

Your system parameters have uncertainty (e.g., qubit frequency drift)
You need pulses robust to fabrication variations
You want to optimize across an ensemble of similar systems
You're designing calibration-free gates

Key Design: Composition Pattern

SamplingProblemwraps an existing QuantumControlProblem and extends it with multiple system variants:

# Create base problem
qcp_base = SmoothPulseProblem(qtraj, N; Q=100.0)
solve!(qcp_base; max_iter=100)

# Create perturbed systems
sys_nominal = get_system(qcp_base)
sys_high = QuantumSystem(1.05 * H_drift, H_drives, drive_bounds)
sys_low = QuantumSystem(0.95 * H_drift, H_drives, drive_bounds)

# Robust optimization
qcp_robust = SamplingProblem(qcp_base, [sys_nominal, sys_high, sys_low])
solve!(qcp_robust; max_iter=100)

Constructor

SamplingProblem(
    qcp::QuantumControlProblem,
    systems::Vector{<:AbstractQuantumSystem};
    weights = fill(1.0, length(systems)),
    Q = 100.0,
    piccolo_options = PiccoloOptions()
)

Parameter Reference

Parameter	Type	Default	Description
`qcp`	`QuantumControlProblem`	required	Base problem providing trajectory structure
`systems`	`Vector{AbstractQuantumSystem}`	required	System variants to optimize over
`weights`	`Vector{Float64}`	`fill(1.0, length(systems))`	Relative importance of each system
`Q`	`Float64`	`100.0`	Infidelity weight (applied to all systems)
`piccolo_options`	`PiccoloOptions`	`PiccoloOptions()`	Solver options

Optimization Objective

SamplingProblem creates a weighted sum of objectives:

minimize: Σᵢ wᵢ × Qᵢ × (1 - Fᵢ) + regularization

Where:

wᵢ is the weight for system i
Qᵢ is the infidelity weight
Fᵢ is the fidelity for system i

All systems share the same control pulse, but each has its own state trajectory.

Examples

Robust to Frequency Drift

using Piccolo

# Nominal system
H_drift = 0.5 * PAULIS[:Z]
H_drives = [PAULIS[:X], PAULIS[:Y]]
sys_nominal = QuantumSystem(H_drift, H_drives, [1.0, 1.0])

# Create trajectory and base problem
T, N = 10.0, 100
times = collect(range(0, T, length = N))
pulse = ZeroOrderPulse(0.1 * randn(2, N), times)
qtraj = UnitaryTrajectory(sys_nominal, pulse, GATES[:X])

qcp_base = SmoothPulseProblem(qtraj, N; Q = 100.0)
cached_solve!(qcp_base, "sampling_base"; max_iter = 100)

# ±5% frequency variation
sys_high = QuantumSystem(1.05 * H_drift, H_drives, [1.0, 1.0])
sys_low = QuantumSystem(0.95 * H_drift, H_drives, [1.0, 1.0])

# Robust optimization
qcp_robust =
    SamplingProblem(qcp_base, [sys_nominal, sys_high, sys_low]; weights = [1.0, 1.0, 1.0])
cached_solve!(qcp_robust, "sampling_robust"; max_iter = 100)

    constructing SmoothPulseProblem for UnitaryTrajectory...
	applying timesteps_all_equal constraint: Δt
┌ Warning: Trajectory has timestep variable :Δt but no bounds on it.
│ Adding default lower bound of 0 to prevent negative timesteps.
│ 
│ Recommended: Add explicit bounds when creating the trajectory:
│   NamedTrajectory(...; Δt_bounds=(min, max))
│ Example:
│   NamedTrajectory(qtraj, N; Δt_bounds=(1e-3, 0.5))
│ 
│ Or use timesteps_all_equal=true in problem options to fix timesteps.
└ @ DirectTrajOpt.Problems ~/.julia/packages/DirectTrajOpt/4TeZc/src/problems.jl:65
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:    38354
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    38584

Total number of variables............................:     1587
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     1188
                     variables with only upper bounds:        0
Total number of equality constraints.................:     1386
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.0002783e+02 7.15e+00 2.32e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.6693747e+01 3.26e+00 3.94e+01  -1.4 5.19e+00   0.0 1.66e-01 5.45e-01f  1
   2  8.6059455e+01 2.24e+00 2.81e+01  -4.0 2.35e+00   0.4 9.23e-02 3.14e-01f  1
   3  4.9990637e+01 2.06e+00 2.63e+01  -4.0 4.92e+00  -0.1 2.62e-02 7.74e-02f  1
   4  5.9481147e-01 1.62e+00 1.98e+01  -1.2 5.47e+00    -  9.93e-02 2.16e-01f  1
⋮
  96  1.0053771e-03 1.57e-05 4.82e-03  -4.3 4.59e-02  -3.7 1.00e+00 1.00e+00h  1
  97  1.0046704e-03 1.03e-05 1.08e-03  -4.3 1.58e-02  -3.3 1.00e+00 1.00e+00h  1
  98  1.0087725e-03 3.04e-05 8.67e-03  -4.3 5.02e-02  -3.7 1.00e+00 1.00e+00h  1
  99  1.0076074e-03 2.32e-07 3.50e-04  -4.3 2.50e-03  -2.4 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  8.5765352e-04 1.82e-05 2.00e+02  -6.4 1.32e-02  -2.9 9.80e-01 9.97e-01h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   8.5765351876338023e-04    8.5765351876338023e-04
Dual infeasibility......:   1.9998315581268940e+02    1.9998315581268940e+02
Constraint violation....:   1.8227986546881514e-05    1.8227986546881514e-05
Variable bound violation:   9.9673094222652026e-09    9.9673094222652026e-09
Complementarity.........:   3.9665319148380616e-06    3.9665319148380616e-06
Overall NLP error.......:   1.9998315581268940e+02    1.9998315581268940e+02


Number of objective function evaluations             = 151
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 151
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 101
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 100
Total seconds in IPOPT                               = 19.072

EXIT: Maximum Number of Iterations Exceeded.
    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
[ Info: Loaded cached trajectory from sampling_base_573ffb2.jld2
    constructing SamplingProblem...
	using 3 systems
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:   132752
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:   117116

Total number of variables............................:     2775
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     2576
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2475
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  5.9703455e+00 9.21e-03 1.00e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.2458275e+01 4.34e-02 1.25e+01  -1.4 8.61e-01   0.0 2.93e-01 3.37e-01h  1
   2  2.5326163e+01 1.44e-02 2.15e+01  -0.8 2.52e-01   1.3 9.71e-01 1.00e+00h  1
   3  2.1415266e+01 2.69e-04 2.61e+00  -1.8 1.79e-02   1.8 9.99e-01 1.00e+00f  1
   4  2.1164636e+01 4.53e-05 7.67e-01  -3.1 1.72e-02   1.3 1.00e+00 1.00e+00f  1
⋮
  96  7.6061830e-03 3.37e-05 2.39e-02  -4.0 1.46e-02  -2.0 1.00e+00 1.00e+00h  1
  97  7.5622363e-03 8.17e-06 3.55e-03  -4.0 5.48e-03  -1.6 1.00e+00 1.00e+00h  1
  98  7.4818316e-03 3.69e-05 1.94e-02  -4.0 1.54e-02  -2.0 1.00e+00 1.00e+00h  1
  99  7.4146603e-03 7.08e-06 6.38e-03  -4.0 5.66e-03  -1.6 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  7.3035617e-03 3.98e-05 3.41e-02  -4.0 1.63e-02  -2.1 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   7.3035616565949816e-03    7.3035616565949816e-03
Dual infeasibility......:   3.4090622132435744e-02    3.4090622132435744e-02
Constraint violation....:   3.9771404570052883e-05    3.9771404570052883e-05
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.0000000000000002e-04    1.0000000000000002e-04
Overall NLP error.......:   3.4090622132435744e-02    3.4090622132435744e-02


Number of objective function evaluations             = 105
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 105
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 101
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 100
Total seconds in IPOPT                               = 142.548

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
        applying constraint: initial value of Ũ⃗1
        applying constraint: initial value of Ũ⃗2
        applying constraint: initial value of Ũ⃗3
        applying constraint: bounds on Ũ⃗1
        applying constraint: bounds on Ũ⃗2
        applying constraint: bounds on Ũ⃗3
        applying constraint: bounds on u
        applying constraint: bounds on Δt
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
[ Info: Loaded cached trajectory from sampling_robust_573ffb2.jld2

Weighted Sampling

Prioritize certain parameter values:

qcp_weighted = SamplingProblem(
    qcp_base,
    [sys_nominal, sys_high, sys_low];
    weights = [2.0, 1.0, 1.0],  ## Nominal weighted 2x
)
cached_solve!(qcp_weighted, "sampling_weighted"; max_iter = 100)

    constructing SamplingProblem...
	using 3 systems
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:   132752
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:   117116

Total number of variables............................:     2775
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     2576
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2475
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  7.9634503e+00 4.54e-02 2.02e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.0192814e+00 2.59e-02 6.19e+00  -1.1 8.27e-02   2.0 9.91e-01 4.38e-01h  1
   2  4.3404151e+00 1.11e-03 1.75e+01  -1.3 5.49e-02   1.5 9.95e-01 1.00e+00f  1
   3  9.9741186e-01 1.06e-04 1.91e+00  -2.6 1.05e-02   1.9 1.00e+00 1.00e+00f  1
   4  2.2472025e-01 3.24e-04 6.75e-01  -3.9 2.65e-02   1.5 1.00e+00 1.00e+00f  1
⋮
  96  6.6432556e-03 1.15e-05 1.00e-02  -4.1 8.54e-03  -1.8 1.00e+00 1.00e+00h  1
  97  6.6883938e-03 2.59e-03 4.13e-01  -4.1 7.59e-02  -2.3 1.00e+00 1.00e+00h  1
  98  6.6433482e-03 3.21e-06 5.40e-03  -4.1 2.97e-03  -0.0 1.00e+00 1.00e+00h  1
  99  6.6442048e-03 3.57e-08 1.34e-04  -4.1 4.31e-04  -0.5 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  6.6586489e-03 5.43e-07 3.88e-02  -4.0 1.61e-03  -1.0 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   3.3330055735692117e-03    6.6586488960568404e-03
Dual infeasibility......:   3.8777317747391543e-02    7.7468980567653487e-02
Constraint violation....:   5.4290098439047085e-07    5.4290098439047085e-07
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.0073676194822144e-04    2.0125100721645598e-04
Overall NLP error.......:   3.8777317747391543e-02    7.7468980567653487e-02


Number of objective function evaluations             = 102
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 102
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 101
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 100
Total seconds in IPOPT                               = 145.884

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
        applying constraint: initial value of Ũ⃗1
        applying constraint: initial value of Ũ⃗2
        applying constraint: initial value of Ũ⃗3
        applying constraint: bounds on Ũ⃗1
        applying constraint: bounds on Ũ⃗2
        applying constraint: bounds on Ũ⃗3
        applying constraint: bounds on u
        applying constraint: bounds on Δt
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
[ Info: Loaded cached trajectory from sampling_weighted_573ffb2.jld2

Dense Parameter Sampling

For smooth performance across a parameter range:

scales = range(0.9, 1.1, length = 3)
systems = [QuantumSystem(s * H_drift, H_drives, [1.0, 1.0]) for s in scales]

qcp_dense = SamplingProblem(qcp_base, systems)
cached_solve!(qcp_dense, "sampling_dense"; max_iter = 100)

    constructing SamplingProblem...
	using 3 systems
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:   132752
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:   117116

Total number of variables............................:     2775
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     2576
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2475
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  5.9738213e+00 5.37e-02 3.17e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  5.4149480e+00 3.60e-02 1.38e+01  -0.8 1.58e-01   2.0 9.82e-01 3.32e-01h  1
   2  1.4886952e+01 1.17e-02 9.34e+01  -0.4 1.42e-01   1.5 8.30e-01 1.00e+00f  1
   3  1.0529015e+00 2.08e-04 3.92e+00  -1.4 3.20e-02   1.9 1.00e+00 9.94e-01f  1
   4  6.0950319e-01 2.52e-04 1.45e+01  -2.4 2.02e-02   1.5 9.99e-01 9.12e-01f  1
⋮
  96  5.3256805e-03 2.92e-06 8.04e-03  -4.1 1.69e-03  -0.9 1.00e+00 1.00e+00h  1
  97  5.4351440e-03 8.14e-07 1.08e-01  -4.1 3.63e-03  -1.4 1.00e+00 1.00e+00H  1
  98  5.4280269e-03 6.85e-07 9.43e-03  -4.1 1.36e-03  -0.9 1.00e+00 1.00e+00h  1
  99  5.4244892e-03 3.54e-06 4.70e-03  -4.1 4.04e-03  -1.4 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  5.4143237e-03 2.54e-05 3.27e-02  -4.1 1.18e-02  -1.9 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   5.4143236729862479e-03    5.4143236729862479e-03
Dual infeasibility......:   3.2727470432375567e-02    3.2727470432375567e-02
Constraint violation....:   2.5360231117622315e-05    2.5360231117622315e-05
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   8.0000531851293414e-05    8.0000531851293414e-05
Overall NLP error.......:   3.2727470432375567e-02    3.2727470432375567e-02


Number of objective function evaluations             = 106
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 106
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 101
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 100
Total seconds in IPOPT                               = 148.478

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
        applying constraint: initial value of Ũ⃗1
        applying constraint: initial value of Ũ⃗2
        applying constraint: initial value of Ũ⃗3
        applying constraint: bounds on Ũ⃗1
        applying constraint: bounds on Ũ⃗2
        applying constraint: bounds on Ũ⃗3
        applying constraint: bounds on u
        applying constraint: bounds on Δt
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
[ Info: Loaded cached trajectory from sampling_dense_573ffb2.jld2

With Time Optimization

Chain with MinimumTimeProblem:

# Base problem with free time
qcp_base = SmoothPulseProblem(qtraj, N; Q=100.0, Δt_bounds=(0.05, 0.5))
solve!(qcp_base; max_iter=100)

# Add robustness
qcp_robust = SamplingProblem(qcp_base, [sys_nominal, sys_high, sys_low])
solve!(qcp_robust; max_iter=100)

# Minimize time
qcp_mintime = MinimumTimeProblem(qcp_robust; final_fidelity=0.95)
solve!(qcp_mintime; max_iter=100)

Trajectory Structure

SamplingProblem creates a SamplingTrajectory internally with:

Variable	Description
`:u`	Shared control values
`:Ũ⃗1`, `:Ũ⃗2`, ...	State for each system (unitary case)
or `:ψ̃1`, `:ψ̃2`, ...	State for each system (ket case)

Note: Derivative variables (:du, :ddu) from the base problem are not carried over. The robustness is achieved through multiple dynamics integrators, one per system.

Difference from MultiKetTrajectory

These serve different purposes:

Feature	`SamplingProblem`	`MultiKetTrajectory`
Systems	Multiple different systems	Single system
States	Same goal across systems	Different initial/goal pairs
Use case	Parameter uncertainty	Multi-state gates

SamplingProblem: "Same gate, different systems" MultiKetTrajectory: "Same system, different state transfers"

Evaluating Robustness

After solving, check fidelity across the parameter range:

# Sample more densely for evaluation
eval_scales = range(0.8, 1.2, length=21)
eval_systems = [QuantumSystem(s * H_drift, H_drives, drive_bounds) for s in eval_scales]

# Get optimized pulse
optimized_pulse = get_pulse(qcp_robust.qtraj)

# Evaluate
for (s, sys) in zip(eval_scales, eval_systems)
    qtraj_eval = UnitaryTrajectory(sys, optimized_pulse, GATES[:X])
    fid = fidelity(qtraj_eval)
    println("Scale $s: Fidelity = $fid")
end