SmoothPulseProblem

SmoothPulseProblem is the most commonly used problem template in Piccolo.jl. It sets up trajectory optimization with piecewise constant controls (ZeroOrderPulse) where control smoothness is enforced through discrete derivative variables.

When to Use

Use SmoothPulseProblem when:

You want piecewise constant control pulses
You need smooth transitions between control values
You're starting a new optimization (not warm-starting from a previous solution)
You want fine control over derivative bounds

Pulse Requirement

SmoothPulseProblem requires a trajectory with a ZeroOrderPulse:

pulse = ZeroOrderPulse(controls, times)
qtraj = UnitaryTrajectory(sys, pulse, U_goal)
qcp = SmoothPulseProblem(qtraj, N)  # Works

Using a spline pulse will result in an error directing you to SplinePulseProblem.

Constructor

SmoothPulseProblem(
    qtraj::AbstractQuantumTrajectory{<:ZeroOrderPulse},
    N::Int;
    kwargs...
)

Parameter Reference

Required Parameters

Parameter	Type	Description
`qtraj`	`AbstractQuantumTrajectory{ZeroOrderPulse}`	Quantum trajectory containing system, pulse, and goal
`N`	`Int`	Number of discretization timesteps

Objective Weights

Parameter	Type	Default	Description
`Q`	`Float64`	`100.0`	Weight on infidelity objective. Higher values prioritize achieving target fidelity.
`R`	`Float64`	`1e-2`	Base regularization weight applied to all derivative terms.
`R_u`	`Float64` or `Vector{Float64}`	`R`	Regularization on control values. Can be per-drive.
`R_du`	`Float64` or `Vector{Float64}`	`R`	Regularization on first derivative (control jumps).
`R_ddu`	`Float64` or `Vector{Float64}`	`R`	Regularization on second derivative (control acceleration).

Bounds

Parameter	Type	Default	Description
`du_bound`	`Float64`	`Inf`	Maximum allowed control jump between timesteps (all drives).
`ddu_bound`	`Float64`	`1.0`	Maximum allowed control acceleration.
`Δt_bounds`	`Tuple{Float64, Float64}`	`nothing`	Time-step bounds `(min, max)` for free-time optimization. Required for `MinimumTimeProblem`.

Advanced Options

Parameter	Type	Default	Description
`integrator`	`AbstractIntegrator`	`nothing`	Custom integrator. If `nothing`, uses `BilinearIntegrator`.
`global_names`	`Vector{Symbol}`	`nothing`	Names of global parameters to optimize (requires custom integrator).
`global_bounds`	`Dict{Symbol, ...}`	`nothing`	Bounds on global variables. Values can be `Float64` (symmetric ±) or `Tuple{Float64, Float64}`.
`constraints`	`Vector{AbstractConstraint}`	`[]`	Additional constraints to add to the problem.
`piccolo_options`	`PiccoloOptions`	`PiccoloOptions()`	Solver and leakage options.
`free_phase`	`Bool`	`false`	Optimize a per-subsystem frame phase alongside the pulse. Only available for `MultiKetTrajectory`; requires `subsystem_levels`.
`initial_phases`	`Vector{Float64}`	`nothing`	Initial values for the per-subsystem phase variables when `free_phase=true` (`MultiKetTrajectory` only).

Examples

Basic Gate Synthesis

using Piccolo

# Define system
H_drift = PAULIS[:Z]
H_drives = [PAULIS[:X], PAULIS[:Y]]
sys = QuantumSystem(H_drift, H_drives, [1.0, 1.0])

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

# Solve
qcp = SmoothPulseProblem(qtraj, N; Q = 100.0, R = 1e-2)
solve!(qcp; max_iter = 100)

fidelity(qcp)

0.999999229699887

With Derivative Bounds

Constrain how fast controls can change:

qcp = SmoothPulseProblem(
    qtraj,
    N;
    Q = 100.0,
    du_bound = 0.5,    ## Limit control jumps
    ddu_bound = 0.1,    ## Limit control acceleration
)

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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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}, Matrix{ComplexF64}}}
  System: QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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

Enabling Free Time for MinimumTimeProblem

To later use MinimumTimeProblem, give bounds on variable timesteps:

qcp = SmoothPulseProblem(
    qtraj,
    N;
    Q = 100.0,
    Δt_bounds = (0.01, 0.5),  ## Timesteps can vary from 0.01 to 0.5
)
solve!(qcp; max_iter = 100)

# Now can minimize time
qcp_mintime = MinimumTimeProblem(qcp; final_fidelity = 0.99)
solve!(qcp_mintime; max_iter = 100)

    constructing SmoothPulseProblem for UnitaryTrajectory...
	applying timesteps_all_equal constraint: Δt
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:        0
                variables with lower and upper bounds:     1288
                     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.9935660e+00 1.37e-01 3.91e-01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  3.0686251e+00 4.12e-04 7.01e-01  -1.2 1.33e-01    -  9.16e-01 1.00e+00h  1
   2  3.4425026e-02 2.02e-04 2.06e-01  -2.9 5.06e-02    -  9.71e-01 9.73e-01f  1
   3  3.8221340e-03 5.25e-06 3.11e-02  -3.8 7.84e-02  -2.0 1.00e+00 9.75e-01h  1
   4  3.4109603e-03 1.79e-05 9.88e-01  -4.0 2.12e-01  -2.5 1.00e+00 1.00e+00h  1
⋮
  96  3.1344359e-03 2.46e-05 1.55e+02  -3.7 4.91e-03   0.5 1.00e+00 3.34e-01h  1
  97  3.6207417e-03 2.82e-08 1.83e+02  -4.0 8.27e-04   0.0 1.43e-01 1.00e+00f  1
  98  3.5014467e-03 2.56e-05 6.59e+01  -3.8 7.66e-02    -  1.00e+00 4.13e-01h  2
  99  3.6043272e-03 2.71e-09 1.94e-03  -3.8 2.62e-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  3.6008210e-03 1.48e-08 6.48e-04  -3.8 5.33e-04  -0.9 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   3.6008209952797005e-03    3.6008209952797005e-03
Dual infeasibility......:   6.4842964282854371e-04    6.4842964282854371e-04
Constraint violation....:   1.4831264688597301e-08    1.4831264688597301e-08
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.6350900156167264e-04    1.6350900156167264e-04
Overall NLP error.......:   6.4842964282854371e-04    6.4842964282854371e-04


Number of objective function evaluations             = 225
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 225
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                               = 20.816

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 Δt
        applying constraint: bounds on du
        applying constraint: bounds on ddu
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
[ Info: Loaded cached trajectory from smooth_pulse_free_time_573ffb2.jld2
    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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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}, Matrix{ComplexF64}}}...
	final fidelity constraint: 0.99
	minimum-time weight D: 100.0
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.:        8
Number of nonzeros in Lagrangian Hessian.............:    38584

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

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  1.8784727e+03 1.00e-02 1.87e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.8746751e+03 9.50e-05 1.85e+00  -2.1 1.85e-02   2.0 9.87e-01 9.90e-01f  1
   2  1.8691722e+03 6.02e-05 1.84e+00  -3.0 5.52e-02   1.5 1.00e+00 1.00e+00f  1
   3  1.8530806e+03 6.64e-04 1.82e+00  -1.9 1.64e-01   1.0 1.00e+00 1.00e+00f  1
   4  1.8472034e+03 6.29e-05 1.85e+00  -2.1 6.25e-02   1.5 1.00e+00 1.00e+00f  1
⋮
  96  1.0925341e+03 1.55e-05 5.72e+00  -3.2 1.57e+00    -  1.77e-01 2.45e-01f  1
  97  1.0925026e+03 1.70e-05 2.06e+01  -3.2 2.11e+00    -  4.33e-01 2.24e-01f  1
  98  1.0924628e+03 1.32e-05 1.92e+01  -4.0 6.95e-01  -1.7 4.75e-01 2.69e-01f  1
  99  1.0924595e+03 1.26e-05 2.29e+01  -3.6 2.93e-01    -  1.00e+00 4.53e-02h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  1.0924243e+03 7.84e-06 4.18e+00  -4.0 2.56e-02    -  1.00e+00 3.78e-01f  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   1.0924128553821258e+03    1.0924243302993173e+03
Dual infeasibility......:   4.1842001571365017e+00    4.1842441087891018e+00
Constraint violation....:   7.8426719529733901e-06    7.8426719529733901e-06
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.9949974748657987e-03    1.9950184307074557e-03
Overall NLP error.......:   4.1842001571365017e+00    4.1842441087891018e+00


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

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 Δt
        applying constraint: bounds on du
        applying constraint: bounds on ddu
        applying constraint: time consistency constraint (t_{k+1} = t_k + Δt_k)
        applying constraint: initial time t₁ = 0
[ Info: Loaded cached trajectory from smooth_pulse_mintime_573ffb2.jld2

Per-Drive Regularization

Apply different regularization to different control channels:

qcp = SmoothPulseProblem(
    qtraj,
    N;
    Q = 100.0,
    R_u = [1e-3, 1e-2],     ## Less regularization on drive 1
    R_du = [1e-2, 1e-1],    ## Different smoothness weights
    R_ddu = [1e-2, 1e-1],
)

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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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.FullSpecialize, SciMLOperators.MatrixOperator{ComplexF64, Matrix{ComplexF64}, SciMLOperators.FilterKwargs{Nothing, Val{()}}, SciMLOperators.FilterKwargs{Piccolo.Quantum.Rollouts.var"#update!#_construct_operator##2"{QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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}, Matrix{ComplexF64}}}
  System: QuantumSystem{Piccolo.Quantum.QuantumSystems.var"#37#38"{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}, Vector{SparseArrays.SparseMatrixCSC{ComplexF64, Int64}}, Int64}, Piccolo.Quantum.QuantumSystems.var"#39#40"{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

With Leakage Suppression

For multilevel systems, use PiccoloOptions to enable leakage handling, in this example a 3-level transmon system:

sys = TransmonSystem(levels = 3, δ = 0.2, drive_bounds = [0.2, 0.2])

# Embedded X gate (only affects computational subspace)
U_goal = EmbeddedOperator(:X, sys)

pulse = ZeroOrderPulse(0.1 * randn(2, N), times)
qtraj = UnitaryTrajectory(sys, pulse, U_goal)

# Enable leakage suppression
opts = PiccoloOptions(leakage_constraint = true, leakage_constraint_value = 1e-3)

qcp = SmoothPulseProblem(qtraj, N; Q = 100.0, piccolo_options = opts)
solve!(qcp; max_iter = 100)

    constructing SmoothPulseProblem for UnitaryTrajectory...
	applying leakage suppression: Ũ⃗ < 0.001
	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/gniE1/src/problems.jl:66
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:   113244
Number of nonzeros in inequality constraint Jacobian.:     9598
Number of nonzeros in Lagrangian Hessian.............:   101039

Total number of variables............................:     2577
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     2178
                     variables with only upper bounds:        0
Total number of equality constraints.................:     2376
Total number of inequality constraints...............:      800
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:      800

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  4.6506424e+01 5.78e+00 2.35e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  4.0342378e+01 5.20e+00 1.38e+01  -0.8 3.26e+00    -  4.55e-02 1.01e-01f  1
   2  3.9661117e+01 5.11e+00 1.64e+01  -1.2 2.92e+00    -  1.00e-01 1.60e-02h  1
   3  3.8454794e+01 4.97e+00 2.00e+02  -0.9 2.88e+00    -  9.04e-02 2.72e-02f  1
   4  3.8179168e+01 4.91e+00 5.39e+02  -1.3 2.83e+00    -  4.31e-02 1.23e-02h  1
⋮
  96  6.3507075e+01 1.22e+00 8.32e+04   0.4 3.04e+01    -  4.82e-03 1.40e-02f  1
  97  6.3466042e+01 1.22e+00 8.15e+04   0.5 1.20e+02   0.9 1.51e-03 7.86e-04f  1
  98  6.3433892e+01 1.21e+00 7.51e+04  -1.5 5.46e+00   2.3 1.49e-02 6.44e-03h  1
  99  6.3247786e+01 1.20e+00 7.54e+04   0.9 8.01e+00   2.7 1.21e-02 1.31e-02f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  6.3240072e+01 1.18e+00 6.12e+04  -0.4 1.17e+00   3.1 5.30e-02 1.45e-02h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   6.3240071511482711e+01    6.3240071511482711e+01
Dual infeasibility......:   6.1193617937769312e+04    6.1193617937769312e+04
Constraint violation....:   1.1808700636995419e+00    1.1808700636995419e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   6.0281175158191864e+00    6.0281175158191864e+00
Overall NLP error.......:   1.6110732330078197e+03    6.1193617937769312e+04


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

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 smooth_pulse_leakage_573ffb2.jld2

State Transfer

Works with KetTrajectory for state preparation. Here we go back to the 2-level system:

sys = QuantumSystem(H_drift, H_drives, [1.0, 1.0])

ψ_init = ComplexF64[1, 0]  ## |0⟩
ψ_goal = ComplexF64[0, 1]  ## |1⟩

pulse = ZeroOrderPulse(0.1 * randn(2, N), times)
qtraj = KetTrajectory(sys, pulse, ψ_init, ψ_goal)

qcp = SmoothPulseProblem(qtraj, N; Q = 100.0)
solve!(qcp; max_iter = 100)

    constructing SmoothPulseProblem for KetTrajectory...
	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/gniE1/src/problems.jl:66
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.

Number of nonzeros in equality constraint Jacobian...:    19430
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:    21862

Total number of variables............................:     1191
                     variables with only lower bounds:      100
                variables with lower and upper bounds:      792
                     variables with only upper bounds:        0
Total number of equality constraints.................:      990
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  9.9799936e+01 6.27e+00 2.02e-01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.9636265e+01 2.76e+00 2.43e+01  -0.3 3.76e+00    -  2.65e-01 5.59e-01f  1
   2  9.8041976e+01 7.04e-04 2.31e+02  -0.7 1.68e+00   2.0 9.61e-01 1.00e+00f  1
   3  9.8216263e+01 1.42e-06 4.45e+02  -4.0 5.84e-03   2.4 6.84e-01 1.00e+00h  1
   4  9.7428169e+01 7.71e-04 4.47e+02   0.6 4.06e+01    -  4.98e-02 1.97e-02f  1
⋮
  96  3.4068479e-03 4.90e-09 3.47e+02  -4.1 9.09e-04  -0.2 1.00e+00 1.00e+00h  1
  97  3.3625349e-03 4.43e-09 5.37e-02  -4.1 3.23e-04   1.1 1.00e+00 1.00e+00h  1
  98  3.3510194e-03 2.06e-09 9.35e-04  -4.1 2.08e-04   0.7 1.00e+00 1.00e+00h  1
  99  3.3525189e-03 1.85e-09 3.47e+02  -4.1 2.30e-04   0.2 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  3.3487504e-03 8.66e-10 3.47e+02  -4.1 9.47e-04  -0.3 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   3.3487504378023137e-03    3.3487504378023137e-03
Dual infeasibility......:   3.4703911823096945e+02    3.4703911823096945e+02
Constraint violation....:   8.6591400716429234e-10    8.6591400716429234e-10
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   7.9999989709547719e-05    7.9999989709547719e-05
Overall NLP error.......:   3.4703911823096945e+02    3.4703911823096945e+02


Number of objective function evaluations             = 175
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 175
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                               = 7.260

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 smooth_pulse_state_transfer_573ffb2.jld2

Multiple State Transfers

Use MultiKetTrajectory for gates defined by state mappings:

ψ0, ψ1 = ComplexF64[1, 0], ComplexF64[0, 1]

# X gate: |0⟩ → |1⟩ and |1⟩ → |0⟩
qtraj = MultiKetTrajectory(sys, pulse, [ψ0, ψ1], [ψ1, ψ0])

qcp = SmoothPulseProblem(qtraj, N; Q = 100.0)
solve!(qcp; max_iter = 100)

    constructing SmoothPulseProblem for MultiKetTrajectory (2 states)...
	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/gniE1/src/problems.jl:66
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  9.9984809e+01 6.27e+00 6.28e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.9642622e+01 2.62e+00 3.50e+01  -0.3 3.75e+00    -  2.90e-01 5.82e-01f  1
   2  9.8853658e+01 6.52e-04 2.21e+02  -0.7 1.59e+00   2.0 9.74e-01 1.00e+00f  1
   3  9.9019463e+01 4.23e-06 4.23e+02  -4.0 8.44e-03   2.4 7.26e-01 1.00e+00h  1
   4  9.7820646e+01 4.98e-04 4.00e+02  -0.1 1.00e+01    -  1.91e-01 7.37e-02f  1
⋮
  96  2.7046004e-03 5.54e-04 2.00e+02  -6.6 1.11e-01    -  8.09e-01 1.00e+00h  1
  97  5.0601758e-03 2.05e-03 9.59e+01  -3.1 2.32e-01  -2.4 4.79e-01 1.00e+00f  1
  98  3.4377805e-03 2.02e-03 2.01e+02  -3.5 5.19e-02  -1.9 1.00e+00 1.68e-02h  1
  99  4.2549974e-03 4.35e-04 1.14e+03  -3.5 6.64e-02  -2.4 6.89e-02 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  4.3639865e-03 1.13e-04 4.94e+00  -3.5 4.36e-02  -2.9 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   4.3639865058668834e-03    4.3639865058668834e-03
Dual infeasibility......:   4.9360555880880685e+00    4.9360555880880685e+00
Constraint violation....:   1.1285948842998383e-04    1.1285948842998383e-04
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   3.1087072410498060e-04    3.1087072410498060e-04
Overall NLP error.......:   4.9360555880880685e+00    4.9360555880880685e+00


Number of objective function evaluations             = 221
Number of objective gradient evaluations             = 101
Number of equality constraint evaluations            = 221
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                               = 22.134

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
        applying constraint: timesteps all equal constraint
        applying constraint: initial value of ψ̃1
        applying constraint: initial value of ψ̃2
        applying constraint: initial value of u
        applying constraint: final value of u
        applying constraint: bounds on ψ̃1
        applying constraint: bounds on ψ̃2
        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 smooth_pulse_multi_ket_573ffb2.jld2

Free-Phase Optimization (MultiKetTrajectory)

When the goal is defined up to a per-subsystem frame rotation (common for multi-level transmons), use free_phase=true. This is supported for the MultiKetTrajectory dispatch. Piccolo adds one phase variable per subsystem and applies a number-operator rotation $e^{i\theta_j \hat{n}_j}$ to the goal states — level $s$ of subsystem $j$ acquires phase $s \theta_j$. For qubits this is a Z-gate rotation; for multi-level systems it generalizes to diag(1, e^{iθ}, e^{2iθ}, ...).

Use initial_phases to warm-start the phases (e.g. from a previous solve):

sys_fp = QuantumSystem(PAULIS[:Z], [PAULIS[:X], PAULIS[:Y]], [1.0, 1.0])
ψ0_fp, ψ1_fp = ComplexF64[1, 0], ComplexF64[0, 1]
pulse_fp = ZeroOrderPulse(0.1 * randn(2, N), times)
qtraj_fp = MultiKetTrajectory(sys_fp, pulse_fp, [ψ0_fp, ψ1_fp], [ψ1_fp, ψ0_fp])

qcp_fp = SmoothPulseProblem(
    qtraj_fp,
    N;
    Q = 100.0,
    free_phase = true,
    subsystem_levels = [2],      ## one 2-level subsystem (qubit)
    initial_phases = [0.1],      ## warm-start the phase variable
)
solve!(qcp_fp; max_iter = 100)

    constructing SmoothPulseProblem for MultiKetTrajectory (2 states)...
    free_phase=true: added phase variables [:φ_1]
	applying timesteps_all_equal constraint: Δt
    added GlobalBoundsConstraint for :φ_1 with bounds ([-6.283185307179586], [6.283185307179586])
┌ 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/gniE1/src/problems.jl:66

******************************************************************************
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
******************************************************************************

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.............:    38593

Total number of variables............................:     1588
                     variables with only lower bounds:      100
                variables with lower and upper bounds:     1189
                     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.0002528e+02 5.71e+00 4.84e-02   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  9.8546388e+01 1.00e+00 7.28e+01  -0.4 3.43e+00   0.0 2.00e-01 8.25e-01f  1
   2  9.8157346e+01 1.64e-04 8.99e+01  -1.5 5.93e-01   1.3 1.40e-01 1.00e+00h  1
   3  9.5082350e+01 1.17e-04 3.97e+00  -1.0 3.39e-01   0.9 9.64e-01 1.00e+00f  1
   4  9.2504894e+01 5.45e-05 4.32e+00  -1.1 1.29e-01   1.3 9.50e-01 1.00e+00f  1
⋮
  96  4.5271296e-03 5.30e-13 5.68e-04  -4.0 1.32e-05   1.6 1.00e+00 1.00e+00h  1
  97  4.5269521e-03 4.75e-12 5.71e-04  -4.0 3.98e-05   1.2 1.00e+00 1.00e+00h  1
  98  4.5264512e-03 4.26e-11 5.73e-04  -4.0 1.20e-04   0.7 1.00e+00 1.00e+00h  1
  99  4.5250359e-03 3.80e-10 5.75e-04  -4.0 3.60e-04   0.2 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
 100  4.5208115e-03 3.38e-09 5.70e-04  -4.0 1.07e-03  -0.3 1.00e+00 1.00e+00h  1

Number of Iterations....: 100

                                   (scaled)                 (unscaled)
Objective...............:   4.5208114578137422e-03    4.5208114578137422e-03
Dual infeasibility......:   5.7045166390267072e-04    5.7045166390267072e-04
Constraint violation....:   3.3835301227469472e-09    3.3835301227469472e-09
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   1.0000000000000000e-04    1.0000000000000000e-04
Overall NLP error.......:   5.7045166390267072e-04    5.7045166390267072e-04


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                               = 29.461

EXIT: Maximum Number of Iterations Exceeded.
    initializing optimizer...
      building evaluator: 4 integrators, 0 nonlinear constraints
      dynamics constraints: 1188, nonlinear constraints: 0
        integrator 1 jacobian structure: 0.138s
        integrator 2 jacobian structure: 0.003s
        integrator 3 jacobian structure: 0.051s
        integrator 4 jacobian structure: 0.001s
      jacobian structure: 38016 nonzeros (0.21s)
        integrator 1 hessian structure: 0.038s
        integrator 2 hessian structure: 0.005s
        integrator 3 hessian structure: 0.026s
        integrator 4 hessian structure: 0.016s
        computing objective hessian structure (CompositeObjective)...
          sub-objective 1 (GlobalKnotPointObjective): 0.294s
          sub-objective 2 (QuadraticRegularizer): 0.186s
          sub-objective 3 (QuadraticRegularizer): 0.001s
          sub-objective 4 (QuadraticRegularizer): 0.001s
          sub-objective 5 (NullObjective): 0.004s
        objective hessian structure: 0.52s
      hessian structure: 38953 nonzeros (0.604s)
      linear index maps built (0.105s)
      evaluator ready (total: 1.011s)
    evaluator created (1.271s)
    NL constraint bounds extracted (0.017s)
    NLP block data built (0.0s)
    Ipopt optimizer configured (0.006s)
    variables set (0.076s)
        applying constraint: timesteps all equal constraint
        applying constraint: bounds constraint on global variable φ_1
        applying constraint: initial value of ψ̃1
        applying constraint: initial value of ψ̃2
        applying constraint: initial value of u
        applying constraint: final value of u
        applying constraint: bounds on ψ̃1
        applying constraint: bounds on ψ̃2
        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
    linear constraints added: 14 (0.447s)
    optimizer initialization complete (total: 1.817s)
[ Info: Loaded cached trajectory from smooth_pulse_free_phase_83da686.jld2

How It Works

SmoothPulseProblem internally:

Adds derivative variables: Creates :du (first derivative) and :ddu (second derivative) variables alongside controls :u
Sets up integrators: Uses DerivativeIntegrator to enforce consistency:
- du[k] = (u[k+1] - u[k]) / Δt
- ddu[k] = (du[k+1] - du[k]) / Δt
Configures objectives:
- Infidelity objective with weight Q
- Quadratic regularization on u, du, ddu with weights R_u, R_du, R_ddu
Applies bounds: Enforces du_bound and ddu_bound as hard constraints

The derivative variables act as auxiliary optimization variables that encourage smooth control transitions without requiring explicit smoothness constraints on the original controls.