Quantum Trajectory Optimization - Lecture 3

using Piccolo
using LinearAlgebra
using CairoMakie
using QuantumToolbox
using SparseArrays

const ⊗ = kron

kron (generic function with 73 methods)

I. Rotations

Goals

• Geodesics
• Controllability of LTI systems
• Dynamical Lie algebras and reachability

Comparing initialization and optimization when restricted to a feasible manifold, and how warm starts are enabled by direct control.

U_goal = PAULIS.X
T = 40
Δt = 0.1

# U_goal = exp(-im * H_eff * T * Δt)
H_eff = im * log(U_goal / T / Δt)

2×2 Matrix{ComplexF64}:
 -1.5708-1.38629im   1.5708+0.0im
  1.5708+0.0im      -1.5708-1.38629im

Exercise

• What happens if there is a drift operator? $H(t) = \Delta \omega Z + u(t) X$.
• Bonus: What about for an Embedded operator?

Controllability

• Quick check: What happens when we kick a system $x_{n+1} = A x_n + B u_n$?

\[\mathcal{C} = \begin{bmatrix} B & A B & A^2 B & \cdots & A^{n-1} B \end{bmatrix}\]

• Quick check: Why did we stop at $n-1$?

Example

• Test on a linear system in 2D. Recall our $F = ma$ system.

function continuous_to_discrete(A, B, h)
    # Construct augmented matrix for matrix exponential
    augmented_matrix = [
        A B;
        zeros(size(B, 2), size(A, 1)) zeros(size(B, 2), size(B, 2))
    ]

    # Compute matrix exponential
    exp_matrix = exp(augmented_matrix * h)

    # Extract discrete LTI system matrices
    A_h = exp_matrix[1:size(A, 1), 1:size(A, 2)]
    B_h = exp_matrix[1:size(A, 1), size(A, 2)+1:end]

    return A_h, B_h
end

# Extract discrete LTI system matrices
A_cts = [0.0 1.0; -1.0 -0.1]
B_cts = [0.0; 1.0]
h = 0.1  # Time step
A, B = continuous_to_discrete(A_cts, B_cts, h)

([0.9950207737420776 0.09933590957864684; -0.09933590957864684 0.9850871827842128], [0.004979226257922404; 0.09933590957864684;;])

• Let's create a 2D state.

\[z = \begin{bmatrix} x \\ \dot{x} \\ y \\ \dot{y} \end{bmatrix}\]

Axy = I(2) ⊗ A
Bxy = [B zeros(2); zeros(2) B]

C = hcat([Axy^n * Bxy for n in 0:size(Axy, 1)-1]...)

4×8 Matrix{Float64}:
 0.00497923  0.0         0.0148221  …  0.0        0.0336788  0.0
 0.0993359   0.0         0.0973599     0.0        0.0906016  0.0
 0.0         0.00497923  0.0           0.0244196  0.0        0.0336788
 0.0         0.0993359   0.0           0.0944356  0.0        0.0906016

rank(C)

4

Axy = I(2) ⊗ A

# Directly move the X particle only
Bxy = [[1; 0] zeros(2); [0; 1] zeros(2)]

C = hcat([Axy^n * Bxy for n in 0:size(Axy, 1)-1]...)

4×8 Matrix{Float64}:
 1.0  0.0   0.995021   0.0   0.980199  0.0   0.955779  0.0
 0.0  0.0  -0.0993359  0.0  -0.196696  0.0  -0.291131  0.0
 0.0  0.0   0.0993359  0.0   0.196696  0.0   0.291131  0.0
 1.0  0.0   0.985087   0.0   0.960529  0.0   0.926666  0.0

rank(C)

2

What about quantum systems? They are nonlinear.

• BCH

\[e^{X} e^{Y} = e^{X + Y + \tfrac{1}{2} [X, Y] + \tfrac{1}{12} ([X, [X, Y]] + [Y, [Y, X]])}\]

• Collect commutators, forming the Dynamical Lie algebra.

• Quick check: How can we test for linear dependence?

# Linearly dependent: QR decomposition has a zero on diagonal
H_drives = [PAULIS.X, PAULIS.Y]
# H_drives = [PAULIS.X, PAULIS.X]

M = stack(vec.(H_drives))
qr(M).R

2×2 Matrix{ComplexF64}:
 -1.41421+0.0im       0.0+0.0im
      0.0+0.0im  -1.41421+0.0im

• Quick check: What about systems with drift?

II. Demos

Bloch sphere

Δ = 0.2
qubit = QuantumSystem(Δ * PAULIS.Z, [PAULIS.X, PAULIS.Y])
ψ0 = ket_from_string("e+g", [2])
ψT = ket_from_bitstring("0")
T = 40
Δt = 0.1
prob = QuantumStateSmoothPulseProblem(qubit, ψ0, ψT, T, Δt)

DirectTrajOptProblem
   timesteps            = 40
   duration             = 3.900000000000002
   variable names       = (:ψ̃, :a, :da, :dda, :Δt)
   knot point dimension = 11

fig, = plot_bloch(prob.trajectory)
fig

plot(prob.trajectory, :a)

solve!(prob, max_iter=30, options=IpoptOptions(eval_hessian=true))

    initializing optimizer...
        applying constraint: timesteps all equal constraint
        applying constraint: initial value of ψ̃
        applying constraint: initial value of a
        applying constraint: final value of a
        applying constraint: bounds on a
        applying constraint: bounds on da
        applying constraint: bounds on dda
        applying constraint: bounds on Δt
This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:     1922
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     1685

Total number of variables............................:      432
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      196
                     variables with only upper bounds:        0
Total number of equality constraints.................:      351
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.3233215e-03 1.74e+00 2.10e+00   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  6.4022946e+01 1.19e-01 2.39e+02  -0.1 1.20e+00    -  6.65e-02 1.00e+00h  1
   2  1.0653966e+01 9.73e-02 1.83e+02   0.3 2.29e+00    -  3.58e-01 2.37e-01H  1
   3  5.8131184e+01 8.42e-03 1.23e+02  -0.1 7.05e-01    -  9.72e-01 1.00e+00f  1
   4  6.0778256e+01 1.43e-04 8.37e+00  -1.3 8.34e-02   2.0 1.00e+00 1.00e+00h  1
   5  5.2625001e+01 1.52e-04 1.74e+00  -2.7 5.21e-02   1.5 1.00e+00 1.00e+00f  1
   6  2.3049155e+01 1.59e-03 2.15e+00  -3.5 1.83e-01   1.0 1.00e+00 1.00e+00f  1
   7  2.0232837e+01 4.73e-03 1.60e+02  -4.0 3.15e-01   0.6 1.00e+00 1.00e+00f  1
   8  6.2438576e+00 2.16e-04 1.44e+02  -4.0 1.00e-01   0.1 1.00e+00 1.00e+00f  1
   9  4.3394200e+00 8.85e-05 4.90e+00  -4.0 4.67e-02   1.4 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  2.2189108e+00 1.04e-04 3.91e-01  -4.0 4.44e-02   0.9 1.00e+00 1.00e+00f  1
  11  2.3196951e-01 1.65e-04 1.81e-01  -4.0 6.19e-02   0.5 1.00e+00 1.00e+00f  1
  12  4.4477119e-01 1.35e-04 1.56e+02  -4.0 6.25e-02  -0.0 1.00e+00 1.00e+00h  1
  13  4.1626205e-01 1.15e-04 1.55e+02  -4.0 5.98e-02  -0.5 1.00e+00 1.00e+00f  1
  14  4.9305126e-02 1.25e-04 4.32e+01  -4.0 2.08e-01   0.8 1.00e+00 2.50e-01f  3
  15  7.5153165e-01 6.05e-06 1.59e+02  -4.0 1.63e-02   1.3 1.00e+00 1.00e+00h  1
  16  6.5791993e-01 1.39e-06 4.07e-01  -4.0 7.90e-03   1.7 1.00e+00 1.00e+00f  1
  17  4.5175637e-01 4.76e-06 2.42e-01  -4.0 1.47e-02   1.2 1.00e+00 1.00e+00f  1
  18  1.7435954e-01 1.31e-05 1.33e-01  -4.0 2.43e-02   0.7 1.00e+00 1.00e+00f  1
  19  3.1074795e-03 1.67e-05 1.61e+02  -4.0 2.28e-02   0.3 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  3.3996619e-02 1.38e-05 4.02e+01  -4.0 2.59e-02  -0.2 1.00e+00 2.50e-01h  3
  21  6.3215650e-02 1.62e-05 1.63e+02  -4.0 1.85e-02   0.2 1.00e+00 1.00e+00h  1
  22  1.7712228e-02 2.44e-06 1.63e+02  -4.0 9.75e-03  -0.3 1.00e+00 1.00e+00f  1
  23  1.4425159e-02 1.56e-07 6.41e-02  -4.0 1.97e-03   1.1 1.00e+00 1.00e+00h  1
  24  7.6569934e-03 3.09e-07 1.28e-02  -4.0 3.34e-03   0.6 1.00e+00 1.00e+00h  1
  25  1.6407745e-03 7.37e-07 1.64e+02  -4.0 5.28e-03   0.1 1.00e+00 1.00e+00h  1
  26  5.2036691e-03 6.91e-07 1.64e+02  -4.0 5.62e-03  -0.4 1.00e+00 1.00e+00h  1
  27  3.5876863e-03 1.02e-07 7.77e-02  -4.1 1.55e-03   1.0 1.00e+00 1.00e+00h  1
  28  1.9042651e-03 1.03e-07 6.13e-03  -4.0 2.02e-03   0.5 1.00e+00 1.00e+00h  1
  29  1.0888445e-03 1.85e-07 1.64e+02  -4.0 2.92e-03   0.0 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  1.3954065e-03 2.24e-03 1.64e+02  -4.0 5.09e-01    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   6.9770325428572251e-04    1.3954065085714450e-03
Dual infeasibility......:   1.6411747563608097e+02    3.2823495127216194e+02
Constraint violation....:   2.2379176476689661e-03    2.2379176476689661e-03
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   9.9999990995489434e-05    1.9999998199097887e-04
Overall NLP error.......:   1.6411747563608097e+02    3.2823495127216194e+02


Number of objective function evaluations             = 38
Number of objective gradient evaluations             = 31
Number of equality constraint evaluations            = 38
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 31
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 30
Total seconds in IPOPT                               = 14.817

EXIT: Maximum Number of Iterations Exceeded.

fig, = plot_bloch(prob.trajectory)
fig

plot(prob.trajectory, :a)

CZ gate

\[H(t) = \sum_j \tfrac{\eta}{2} a_j^\dagger a_j^\dagger a_j a_j + \Delta(t) (a_1^\dagger a_1) + g(t) (a_1^\dagger a_2 + a_1 a_2^\dagger)\]

n_levels = 2
a = lift_operator(annihilate(n_levels), 1, 2)
b = lift_operator(annihilate(n_levels), 2, 2)
η = -0.3

H_drift = η/2 * (a'a'*a*a + b'b'*b*b)
H_drives = [a'a', (a'b + a*b')]
transmons = QuantumSystem(H_drift, H_drives)

U_goal = EmbeddedOperator(
    GATES.CZ,
    get_subspace_indices([1:2, 1:2], [n_levels, n_levels]),
    [n_levels, n_levels]
)

PiccoloQuantumObjects.EmbeddedOperators.EmbeddedOperator{ComplexF64}(ComplexF64[1.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 1.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im 1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im -1.0 + 0.0im], [1, 2, 3, 4], [2, 2])

commutator(A, B) = A * B - B * A

commutator (generic function with 1 method)

A = H_drives[1]
B = H_drives[2]
commutator(A, B) |> sparse

4×4 SparseArrays.SparseMatrixCSC{ComplexF64, Int64} with 0 stored entries:
     ⋅          ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅          ⋅    
     ⋅          ⋅          ⋅          ⋅    

try
    is_reachable(U_goal, transmons)
catch e
    println(e)
end

operator algebra depth = [1]
ArgumentError("`stack` on an empty collection is not allowed")

Continuous-variable quantum computing

\[H(t) = \Omega(t) a^\dagger + \Omega(t) a + \kappa(t) a^2 + \kappa^*(t) (a^\dagger)^2\]

n_levels = 5

a = annihilate(n_levels)
X = a + a'
Y = -im * (a - a')
X2 = a^2 + (a')^2
Y2 = -im * (a^2 - (a')^2)

Ω = 1.0
κ = 0.1

sys = QuantumSystem([X, Y, X2, Y2])

# Displacement and squeezing operators
function displacement(α)
    return exp(α * a' - conj(α) * a)
end

function squeezing(r)
    return exp((r / 2) * (a^2 - (a')^2))
end

# Initial states
ψ0 = I(n_levels)[:, 1] .+ 0.0im
ψα = displacement(im * 1.5) * ψ0;
ψs = squeezing(0.5) * displacement(0.5 + 0.5im) * ψ0;

5-element Vector{ComplexF64}:
   0.7274368279485289 + 0.08379539462600397im
   0.2732317136965428 + 0.36403586011179717im
 -0.26906111979151676 + 0.1634807537631869im
  -0.2884595663616072 - 0.15919291236653724im
  0.08267502401506288 - 0.20525595962886536im

Coherent state

T = 40
Δt = 0.2
prob = QuantumStateSmoothPulseProblem(sys, ψ0, ψα, T, Δt, dda_bound=0.1)
solve!(prob, max_iter=30, options=IpoptOptions(eval_hessian=true))

    constructing QuantumStateSmoothPulseProblem...
	using integrator: typeof(QuantumCollocation.QuantumIntegrators.KetIntegrator)
	using 1 initial state(s)
	applying timesteps_all_equal constraint: Δt
    initializing optimizer...
        applying constraint: timesteps all equal constraint
        applying constraint: initial value of ψ̃
        applying constraint: initial value of a
        applying constraint: final value of a
        applying constraint: bounds on a
        applying constraint: bounds on da
        applying constraint: bounds on dda
        applying constraint: bounds on Δt
This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:     7418
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:     4419

Total number of variables............................:      902
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      352
                     variables with only upper bounds:        0
Total number of equality constraints.................:      741
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.0828501e-02 1.94e+00 1.07e+01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  8.9882803e+00 1.60e+00 1.77e+02  -0.0 1.09e+00   2.0 9.82e-02 1.72e-01f  1
   2  1.1620073e+01 9.49e-01 3.88e+02   0.0 8.81e-01   2.4 1.00e+00 4.08e-01f  1
   3  8.8190617e+00 4.53e-01 4.51e+02   0.6 5.25e-01   2.9 1.00e+00 5.22e-01h  1
   4  1.2837611e+01 1.71e-01 7.24e+02   0.8 2.74e-01   3.3 1.00e+00 6.20e-01h  1
   5  1.4902916e+01 9.17e-02 1.04e+03  -4.0 6.92e-01   2.8 1.94e-01 4.64e-01h  1
   6  2.1102296e+01 1.06e-02 4.48e+02   0.3 1.24e-01   3.2 1.00e+00 1.00e+00f  1
   7  2.0395859e+01 4.33e-04 4.35e+01  -0.4 3.91e-02   2.7 1.00e+00 1.00e+00f  1
   8  1.6807694e+01 4.59e-04 7.15e+00  -1.4 3.68e-02   2.3 1.00e+00 1.00e+00f  1
   9  1.0411555e+01 1.43e-03 7.49e+00  -2.4 6.69e-02   1.8 9.98e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  9.8346830e+00 1.14e-04 4.75e+00  -3.6 2.45e-02   2.2 1.00e+00 1.00e+00f  1
  11  6.4962736e+00 6.36e-04 5.02e+00  -4.0 5.12e-02   1.7 1.00e+00 1.00e+00f  1
  12  5.8589910e+00 9.67e-05 3.47e+00  -4.0 1.95e-02   2.2 1.00e+00 1.00e+00f  1
  13  3.6994579e+00 4.64e-04 3.52e+00  -4.0 3.98e-02   1.7 1.00e+00 1.00e+00f  1
  14  3.2997432e+00 6.58e-05 2.25e+00  -4.0 1.40e-02   2.1 1.00e+00 1.00e+00f  1
  15  2.6449569e+00 1.01e-04 2.28e+00  -4.0 2.77e-02   1.6 1.00e+00 4.94e-01f  1
  16  2.2379577e+00 4.71e-05 1.68e+00  -4.0 1.14e-02   2.1 1.00e+00 1.00e+00f  1
  17  2.1055139e+00 6.65e-06 1.38e+00  -4.0 4.05e-03   2.5 1.00e+00 1.00e+00f  1
  18  1.7781658e+00 2.98e-05 1.38e+00  -4.0 1.01e-02   2.0 1.00e+00 8.37e-01f  1
  19  1.6602575e+00 5.87e-06 1.12e+00  -4.0 3.65e-03   2.4 1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  1.3303777e+00 3.54e-05 1.12e+00  -4.0 9.00e-03   2.0 1.00e+00 1.00e+00f  1
  21  1.2432955e+00 4.82e-06 8.65e-01  -4.0 3.16e-03   2.4 1.00e+00 1.00e+00f  1
  22  9.8426438e-01 2.67e-05 8.57e-01  -4.0 7.64e-03   1.9 1.00e+00 1.00e+00f  1
  23  9.1833490e-01 3.51e-06 6.53e-01  -4.0 2.65e-03   2.3 1.00e+00 1.00e+00f  1
  24  7.2040600e-01 1.81e-05 6.45e-01  -4.0 6.25e-03   1.9 1.00e+00 1.00e+00f  1
  25  6.7178547e-01 2.29e-06 4.83e-01  -4.0 2.28e-03   2.3 1.00e+00 1.00e+00f  1
  26  5.2399137e-01 1.09e-05 4.81e-01  -4.0 5.52e-03   1.8 1.00e+00 1.00e+00f  1
  27  4.8881958e-01 1.51e-06 3.63e-01  -4.0 1.98e-03   2.2 1.00e+00 1.00e+00h  1
  28  3.8004068e-01 1.06e-05 4.50e-01  -4.0 4.70e-03   1.8 1.00e+00 1.00e+00f  1
  29  3.5486985e-01 1.45e-06 2.88e-01  -4.0 1.66e-03   2.2 1.00e+00 1.00e+00h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  2.7563720e-01 9.95e-06 4.09e-01  -4.0 4.02e-03   1.7 1.00e+00 1.00e+00f  1

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   2.5438859752659493e-01    2.7563720249999463e-01
Dual infeasibility......:   4.0923429159578362e-01    4.4341686852035972e-01
Constraint violation....:   9.9538102384723077e-06    9.9538102384723077e-06
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   9.9999999999767593e-05    1.0835281344343183e-04
Overall NLP error.......:   4.0923429159578362e-01    4.4341686852035972e-01


Number of objective function evaluations             = 31
Number of objective gradient evaluations             = 31
Number of equality constraint evaluations            = 31
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 31
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 30
Total seconds in IPOPT                               = 18.706

EXIT: Maximum Number of Iterations Exceeded.

rollout_fidelity(prob.trajectory, sys)

0.9971892582932071

fig, = plot_wigner(prob.trajectory, 1)
fig

fig, = plot_wigner(prob.trajectory, prob.trajectory.T)
fig

fig, = plot_wigner(QuantumObject(ψα))
fig

plot(
    prob.trajectory, [:a],
    transformations = [:ψ̃ => ψ̃ -> abs2.(iso_to_ket(ψ̃)),],
    use_autolimits=true,
    transformation_titles = ["Population"],
)

Minimum Time

min_prob = QuantumStateMinimumTimeProblem(prob, ψα)

DirectTrajOptProblem
   timesteps            = 40
   duration             = 15.023833275928611
   variable names       = (:ψ̃, :a, :da, :dda, :Δt)
   knot point dimension = 23

solve!(min_prob, max_iter=30)

    initializing optimizer...
        applying constraint: timesteps all equal constraint
        applying constraint: initial value of ψ̃
        applying constraint: initial value of a
        applying constraint: final value of a
        applying constraint: bounds on a
        applying constraint: bounds on da
        applying constraint: bounds on dda
        applying constraint: bounds on Δt
This is Ipopt version 3.14.17, running with linear solver MUMPS 5.8.0.

Number of nonzeros in equality constraint Jacobian...:     7408
Number of nonzeros in inequality constraint Jacobian.:       10
Number of nonzeros in Lagrangian Hessian.............:     4474

Total number of variables............................:      902
                     variables with only lower bounds:        0
                variables with lower and upper bounds:      352
                     variables with only upper bounds:        0
Total number of equality constraints.................:      741
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.5026590e+03 2.75e-03 6.19e+01   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  8.4192840e+02 5.34e-02 2.03e+02  -0.1 1.02e+00   2.0 1.00e+00 2.36e-01f  1
   2  8.0608345e+02 4.67e-02 8.02e+02   0.3 1.79e-01   2.4 1.00e+00 1.25e-01f  1
   3  4.8321797e+02 1.96e-02 5.12e+02   0.5 3.22e-01   1.9 1.00e+00 6.07e-01f  1
   4  6.3434653e+02 6.16e-03 1.17e+02   0.4 1.96e-01    -  9.96e-01 9.35e-01h  1
   5  6.0445164e+02 1.93e-03 3.39e+02  -0.0 7.65e-02    -  9.96e-01 7.33e-01f  1
   6  5.7033787e+02 3.81e-04 7.20e+01  -0.7 4.45e-02    -  1.00e+00 1.00e+00f  1
   7  5.6167328e+02 1.48e-04 7.99e+02  -1.2 5.33e-02    -  1.00e+00 7.04e-01f  1
   8  5.5559917e+02 2.98e-04 4.97e+02  -1.7 7.30e-02    -  1.00e+00 8.93e-01f  1
   9  5.5697132e+02 2.17e-04 7.38e+03  -1.1 2.34e-01    -  1.00e+00 2.42e-01f  2
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  5.5699160e+02 1.87e-05 1.42e+03  -1.6 2.98e-02    -  1.00e+00 9.16e-01h  1
  11  5.5701072e+02 1.81e-05 6.13e+04  -4.0 1.56e-01    -  3.16e-01 3.11e-02h  3
  12  5.5750652e+02 2.40e-05 5.67e+04  -1.4 1.88e+00    -  1.60e-01 8.45e-02h  1
  13  5.5773174e+02 2.23e-05 5.09e+04  -1.8 1.12e+00    -  1.33e-01 6.15e-02h  2
  14  5.5795339e+02 2.05e-05 3.78e+04  -1.8 9.28e-01    -  4.07e-01 7.43e-02h  2
  15  5.5778635e+02 1.53e-04 2.85e+04  -2.7 3.00e-01    -  2.87e-01 2.67e-01f  1
  16  5.5782765e+02 1.52e-04 1.58e+04  -2.0 1.34e+00    -  3.29e-02 8.08e-03h  4
  17  5.5817749e+02 5.73e-05 2.79e+03  -2.0 2.40e-01   1.5 6.46e-01 8.19e-01h  1
  18  5.5787366e+02 7.14e-05 3.66e+03  -2.0 5.01e+00   1.0 1.76e-03 5.60e-03f  1
  19  5.5811494e+02 5.84e-05 2.31e+03  -2.0 4.86e-01   1.4 2.46e-01 2.09e-01h  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20  5.5824377e+02 5.73e-05 1.34e+03  -2.0 7.28e-02   0.9 4.83e-01 4.05e-01h  1
  21  5.5839143e+02 3.32e-05 3.13e+04  -2.0 4.38e-01    -  5.50e-01 3.01e-01H  1
  22  5.5875637e+02 2.43e-05 1.68e+04  -2.0 2.30e-01    -  3.56e-01 2.48e-01h  2
  23  5.5911825e+02 2.30e-05 1.31e+04  -2.0 2.08e-01    -  4.20e-01 3.12e-01h  1
  24  5.5938348e+02 1.45e-05 1.33e+04  -2.0 9.67e-02    -  6.28e-01 5.00e-01h  2
  25  5.5923090e+02 8.91e-06 3.44e+04  -2.5 6.22e-02    -  8.46e-01 4.56e-01f  2
  26  5.5895532e+02 1.78e-05 1.29e+04  -2.9 5.00e-02   0.5 1.00e+00 8.54e-01f  1
  27  5.5885723e+02 4.12e-06 3.49e+03  -3.7 2.06e-02   0.9 1.00e+00 9.49e-01f  1
  28  5.5884963e+02 2.60e-06 1.26e+05  -4.0 4.67e-02   0.4 1.00e+00 3.76e-01h  1
  29  5.5882103e+02 3.16e-06 1.80e+05  -3.9 1.75e-01  -0.1 1.00e+00 1.45e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30  5.5880918e+02 3.72e-06 1.65e+05  -2.8 2.50e-01  -0.5 1.00e+00 8.48e-02f  1

Number of Iterations....: 30

                                   (scaled)                 (unscaled)
Objective...............:   5.1797514152728161e+02    5.5880918342810162e+02
Dual infeasibility......:   1.6484330731854338e+05    1.7783856129592101e+05
Constraint violation....:   3.7211417886099741e-06    3.7211417886099741e-06
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   2.1853615288838571e-03    2.3576422757470412e-03
Overall NLP error.......:   3.6446764822853936e+02    1.7783856129592101e+05


Number of objective function evaluations             = 56
Number of objective gradient evaluations             = 31
Number of equality constraint evaluations            = 56
Number of inequality constraint evaluations          = 56
Number of equality constraint Jacobian evaluations   = 31
Number of inequality constraint Jacobian evaluations = 31
Number of Lagrangian Hessian evaluations             = 30
Total seconds in IPOPT                               = 10.301

EXIT: Maximum Number of Iterations Exceeded.

fig, = plot_wigner(min_prob.trajectory, min_prob.trajectory.T)
fig

plot(
    min_prob.trajectory, [:a],
    transformations = [:ψ̃ => ψ̃ -> abs2.(iso_to_ket(ψ̃)),],
    use_autolimits=true,
    transformation_titles = ["Population"],
)

Squeezed states

Further exercises and explorations with squeezed states can be implemented here.