diff --git a/Project.toml b/Project.toml index 939f1d8f..8aff5bdd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "QuantumCollocation" uuid = "0dc23a59-5ffb-49af-b6bd-932a8ae77adf" authors = ["Aaron Trowbridge and contributors"] -version = "0.7.0" +version = "0.7.2" [deps] DirectTrajOpt = "c823fa1f-8872-4af5-b810-2b9b72bbbf56" @@ -21,14 +21,14 @@ TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe" TrajectoryIndexingUtils = "6dad8b7f-dd9a-4c28-9b70-85b9a079bfc8" [compat] -DirectTrajOpt = "0.1.0" +DirectTrajOpt = "0.2" Distributions = "0.25" ExponentialAction = "0.2" Interpolations = "0.15" JLD2 = "0.5" LinearAlgebra = "1.10, 1.11" NamedTrajectories = "0.3" -PiccoloQuantumObjects = "0.3" +PiccoloQuantumObjects = "0.4" Random = "1.10, 1.11" Reexport = "1.2" SparseArrays = "1.10, 1.11" diff --git a/README.md b/README.md index 56e918f8..154bc062 100644 --- a/README.md +++ b/README.md @@ -67,30 +67,7 @@ ``` where $\mathbf{Z}$ is a trajectory containing states and controls, from [NamedTrajectories.jl](https://github.com/harmoniqs/NamedTrajectories.jl). -### Problem Templates - -*Problem Templates* are reusable design patterns for setting up and solving common quantum control problems. - -For example, a *UnitarySmoothPulseProblem* is tasked with generating a *pulse* sequence $a_{1:T-1}$ in orderd to minimize infidelity, subject to constraints from the Schroedinger equation, -```math - \begin{aligned} - \arg \min_{\mathbf{Z}}\quad & |1 - \mathcal{F}(U_T, U_\text{goal})| \\ - \nonumber \text{s.t.} - \qquad & U_{t+1} = \exp\{- i H(a_t) \Delta t_t \} U_t, \quad \forall\, t \\ - \end{aligned} -``` -while a *UnitaryMinimumTimeProblem* minimizes time and constrains fidelity, -```math - \begin{aligned} - \arg \min_{\mathbf{Z}}\quad & \sum_{t=1}^T \Delta t_t \\ - \qquad & U_{t+1} = \exp\{- i H(a_t) \Delta t_t \} U_t, \quad \forall\, t \\ - \nonumber & \mathcal{F}(U_T, U_\text{goal}) \ge 0.9999 - \end{aligned} -``` - -In each case, the dynamics between *knot points* $(U_t, a_t)$ and $(U_{t+1}, a_{t+1})$ are enforced as constraints on the states, which are free variables in the solver; this optimization framework is called *direct collocation*. For details of our implementation please see our award-winning IEEE QCE 2023 paper, [Direct Collocation for Quantum Optimal Control](https://arxiv.org/abs/2305.03261). If you use QuantumCollocation.jl in your work, please cite :raised_hands:! - -Problem templates give the user the ability to add other constraints and objective functions to this problem and solve it efficiently using [Ipopt.jl](https://github.com/jump-dev/Ipopt.jl) and [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) under the hood. +For details of our implementation please see our IEEE QCE 2023 paper, [Direct Collocation for Quantum Optimal Control](https://arxiv.org/abs/2305.03261). If you use QuantumCollocation.jl in your work, please cite :raised_hands:! ## Installation diff --git a/docs/literate/man/ket_problem_templates.jl b/docs/literate/man/ket_problem_templates.jl new file mode 100644 index 00000000..9d2283ba --- /dev/null +++ b/docs/literate/man/ket_problem_templates.jl @@ -0,0 +1,97 @@ +# ```@meta +# CollapsedDocStrings = true +# ``` +using NamedTrajectories +using PiccoloQuantumObjects +using QuantumCollocation + +# ----- + +#= +## Quantum State Smooth Pulse Problem + +```@docs; canonical = false +QuantumStateSmoothPulseProblem +``` + +Each problem starts with a `QuantumSystem` object, which is used to define the system's +Hamiltonian and control operators. The goal is to find a control pulse that drives the +intial state, `ψ_init`, to a target state, `ψ_goal`. +=# + +# _define the quantum system_ +system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) +ψ_init = Vector{ComplexF64}([1.0, 0.0]) +ψ_goal = Vector{ComplexF64}([0.0, 1.0]) +T = 51 +Δt = 0.2 + +# _create the smooth pulse problem_ +state_prob = QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, T, Δt); + +# _check the fidelity before solving_ +println("Before: ", rollout_fidelity(state_prob.trajectory, system)) + +# _solve the problem_ +solve!(state_prob, max_iter=100, verbose=true, print_level=1); + +# _check the fidelity after solving_ +println("After: ", rollout_fidelity(state_prob.trajectory, system)) + +# _extract the control pulses_ +state_prob.trajectory.a |> size + +# ----- + +#= +## Quantum State Minimum Time Problem + +```@docs; canonical = false +QuantumStateMinimumTimeProblem +``` +=# + +# _create the minimum time problem_ +min_state_prob = QuantumStateMinimumTimeProblem(state_prob, ψ_goal); + +# _check the previous duration_ +println("Duration before: ", get_duration(state_prob.trajectory)) + +# _solve the minimum time problem_ +solve!(min_state_prob, max_iter=100, verbose=true, print_level=1); + +# _check the new duration_ +println("Duration after: ", get_duration(min_state_prob.trajectory)) + +# _the fidelity is preserved by a constraint_ +println("Fidelity after: ", rollout_fidelity(min_state_prob.trajectory, system)) + +# ----- + +#= + +## Quantum State Sampling Problem + +```@docs; canonical = false +QuantumStateSamplingProblem +``` +=# + +# _create a sampling problem_ +driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) +sampling_state_prob = QuantumStateSamplingProblem([system, driftless_system], ψ_init, ψ_goal, T, Δt); + +# _new keys are added to the trajectory for the new states_ +println(sampling_state_prob.trajectory.state_names) + +# _solve the sampling problem for a few iterations_ +solve!(sampling_state_prob, max_iter=25, verbose=true, print_level=1); + +# _check the fidelity of the sampling problem (use the updated key to get the initial and goal)_ +println("After (original system): ", rollout_fidelity(sampling_state_prob.trajectory, system, state_name=:ψ̃1_system_1)) +println("After (new system): ", rollout_fidelity(sampling_state_prob.trajectory, driftless_system, state_name=:ψ̃1_system_1)) + +# _compare this to using the original problem on the new system_ +println("After (new system, original `prob`): ", rollout_fidelity(state_prob.trajectory, driftless_system)) + +# ----- \ No newline at end of file diff --git a/docs/literate/man/problem_templates.jl b/docs/literate/man/problem_templates.jl deleted file mode 100644 index 75107fd0..00000000 --- a/docs/literate/man/problem_templates.jl +++ /dev/null @@ -1,197 +0,0 @@ -# ```@meta -# CollapsedDocStrings = true -# ``` -using NamedTrajectories -using PiccoloQuantumObjects -using QuantumCollocation - -#= -# Problem Templates - -We provide a number of problem templates for making it simple and easy to set up and solve -certain types of quantum optimal control problems. These templates all construct a -`DirectTrajOptProblem` object, which stores all the parts of the optimal control problem. - -This page provides a brief overview of each problem template, broken down by the state of -the problem being solved. - -Ket Problem Templates: -- [Quantum State Smooth Pulse Problem](#Quantum-State-Smooth-Pulse-Problem) -- [Quantum State Minimum Time Problem](#Quantum-State-Minimum-Time-Problem) -- [Quantum State Sampling Problem](#Quantum-State-Sampling-Problem) - -Unitary Problem Templates: -- [Unitary Smooth Pulse Problem](#Unitary-Smooth-Pulse-Problem) -- [Unitary Minimum Time Problem](#Unitary-Minimum-Time-Problem) -- [Unitary Sampling Problem](#Unitary-Sampling-Problem) -=# - -# ----- -# ## Ket Problem Templates - -#= -#### Quantum State Smooth Pulse Problem -```@docs; canonical = false -QuantumStateSmoothPulseProblem -``` - -Each problem starts with a `QuantumSystem` object, which is used to define the system's -Hamiltonian and control operators. The goal is to find a control pulse that drives the -intial state, `ψ_init`, to a target state, `ψ_goal`. -=# - -# _define the quantum system_ -system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) -ψ_init = Vector{ComplexF64}([1.0, 0.0]) -ψ_goal = Vector{ComplexF64}([0.0, 1.0]) -T = 51 -Δt = 0.2 - -# _create the smooth pulse problem_ -state_prob = QuantumStateSmoothPulseProblem(system, ψ_init, ψ_goal, T, Δt); - -# _check the fidelity before solving_ -println("Before: ", rollout_fidelity(state_prob.trajectory, system)) - -# _solve the problem_ -solve!(state_prob, max_iter=100, verbose=true, print_level=1); - -# _check the fidelity after solving_ -println("After: ", rollout_fidelity(state_prob.trajectory, system)) - -# _extract the control pulses_ -state_prob.trajectory.a |> size - -#= -#### Quantum State Minimum Time Problem -```@docs; canonical = false -QuantumStateMinimumTimeProblem -``` -=# - -# _create the minimum time problem_ -min_state_prob = QuantumStateMinimumTimeProblem(state_prob, ψ_goal); - -# _check the previous duration_ -println("Duration before: ", get_duration(state_prob.trajectory)) - -# _solve the minimum time problem_ -solve!(min_state_prob, max_iter=100, verbose=true, print_level=1); - -# _check the new duration_ -println("Duration after: ", get_duration(min_state_prob.trajectory)) - -# _the fidelity is preserved by a constraint_ -println("Fidelity after: ", rollout_fidelity(min_state_prob.trajectory, system)) - -#= -#### Quantum State Sampling Problem -```@docs; canonical = false -QuantumStateSamplingProblem -``` -=# - -# _create a sampling problem_ -driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) -sampling_state_prob = QuantumStateSamplingProblem([system, driftless_system], ψ_init, ψ_goal, T, Δt); - -# _new keys are added to the trajectory for the new states_ -println(sampling_state_prob.trajectory.state_names) - -# _solve the sampling problem for a few iterations_ -solve!(sampling_state_prob, max_iter=25, verbose=true, print_level=1); - -# _check the fidelity of the sampling problem (use the updated key to get the initial and goal)_ -println("After (original system): ", rollout_fidelity(sampling_state_prob.trajectory, system, state_name=:ψ̃1_system_1)) -println("After (new system): ", rollout_fidelity(sampling_state_prob.trajectory, driftless_system, state_name=:ψ̃1_system_1)) - -# _compare this to using the original problem on the new system_ -println("After (new system, original `prob`): ", rollout_fidelity(state_prob.trajectory, driftless_system)) - - -# ----- -# ## Unitary Problem Templates - -#= -#### Unitary Smooth Pulse Problem - -```@docs; canonical = false -UnitarySmoothPulseProblem -``` - -The `UnitarySmoothPulseProblem` is similar to the `QuantumStateSmoothPulseProblem`, but -instead of driving the system to a target state, the goal is to drive the system to a -target unitary operator, `U_goal`. - -=# - -system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) -U_goal = GATES.H -T = 51 -Δt = 0.2 - -prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt); - -# _check the fidelity before solving_ -println("Before: ", unitary_rollout_fidelity(prob.trajectory, system)) - -# _finding an optimal control is as simple as calling `solve!`_ -solve!(prob, max_iter=100, verbose=true, print_level=1); - -# _check the fidelity after solving_ -println("After: ", unitary_rollout_fidelity(prob.trajectory, system)) - -#= -The `NamedTrajectory` object stores the control pulse, state variables, and the time grid. -=# - -# _extract the control pulses_ -prob.trajectory.a |> size - -#= -#### Unitary Minimum Time Problem - -```@docs; canonical = false -UnitaryMinimumTimeProblem -``` - -The goal of this problem is to find the shortest time it takes to drive the system to a -target unitary operator, `U_goal`. The problem is solved by minimizing the sum of all of -the time steps. It is constructed from `prob` in the previous example. -=# - -min_prob = UnitaryMinimumTimeProblem(prob, U_goal); - -# _check the previous duration_ -println("Duration before: ", get_duration(prob.trajectory)) - -# _solve the minimum time problem_ -solve!(min_prob, max_iter=100, verbose=true, print_level=1); - -# _check the new duration_ -println("Duration after: ", get_duration(min_prob.trajectory)) - -# _the fidelity is preserved by a constraint_ -println("Fidelity after: ", unitary_rollout_fidelity(min_prob.trajectory, system)) - -#= -#### Unitary Sampling Problem - -```@docs; canonical = false -UnitarySamplingProblem -``` - -A sampling problem is used to solve over multiple quantum systems with the same control. -This can be useful for exploring robustness, for example. -=# - -# _create a sampling problem_ -driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) -sampling_prob = UnitarySamplingProblem([system, driftless_system], U_goal, T, Δt); - -# _new keys are addded to the trajectory for the new states_ -println(sampling_prob.trajectory.state_names) - -# _the `solve!` proceeds as in the [Quantum State Sampling Problem](#Quantum-State-Sampling-Problem)]_ - - diff --git a/docs/literate/man/unitary_problem_templates.jl b/docs/literate/man/unitary_problem_templates.jl new file mode 100644 index 00000000..2121f883 --- /dev/null +++ b/docs/literate/man/unitary_problem_templates.jl @@ -0,0 +1,118 @@ +# ```@meta +# CollapsedDocStrings = true +# ``` +using NamedTrajectories +using PiccoloQuantumObjects +using QuantumCollocation + +# ----- + +#= +## Unitary Smooth Pulse Problem + +```@docs; canonical = false +UnitarySmoothPulseProblem +``` + +The `UnitarySmoothPulseProblem` is similar to the `QuantumStateSmoothPulseProblem`, but +instead of driving the system to a target state, the goal is to drive the system to a +target unitary operator, `U_goal`. + +=# + +system = QuantumSystem(0.1 * PAULIS.Z, [PAULIS.X, PAULIS.Y]) +U_goal = GATES.H +T = 51 +Δt = 0.2 + +prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt); + +# _check the fidelity before solving_ +println("Before: ", unitary_rollout_fidelity(prob.trajectory, system)) + +# _finding an optimal control is as simple as calling `solve!`_ +solve!(prob, max_iter=100, verbose=true, print_level=1); + +# _check the fidelity after solving_ +println("After: ", unitary_rollout_fidelity(prob.trajectory, system)) + +#= +The `NamedTrajectory` object stores the control pulse, state variables, and the time grid. +=# + +# _extract the control pulses_ +prob.trajectory.a |> size + +# ----- + +#= +## Unitary Minimum Time Problem + +```@docs; canonical = false +UnitaryMinimumTimeProblem +``` + +The goal of this problem is to find the shortest time it takes to drive the system to a +target unitary operator, `U_goal`. The problem is solved by minimizing the sum of all of +the time steps. It is constructed from `prob` in the previous example. +=# + +min_prob = UnitaryMinimumTimeProblem(prob, U_goal); + +# _check the previous duration_ +println("Duration before: ", get_duration(prob.trajectory)) + +# _solve the minimum time problem_ +solve!(min_prob, max_iter=100, verbose=true, print_level=1); + +# _check the new duration_ +println("Duration after: ", get_duration(min_prob.trajectory)) + +# _the fidelity is preserved by a constraint_ +println("Fidelity after: ", unitary_rollout_fidelity(min_prob.trajectory, system)) + +# ----- + +#= +## Unitary Sampling Problem + +```@docs; canonical = false +UnitarySamplingProblem +``` + +A sampling problem is used to solve over multiple quantum systems with the same control. +This can be useful for exploring robustness, for example. +=# + +# _create a sampling problem_ +driftless_system = QuantumSystem([PAULIS.X, PAULIS.Y]) +sampling_prob = UnitarySamplingProblem([system, driftless_system], U_goal, T, Δt); + +# _new keys are addded to the trajectory for the new states_ +println(sampling_prob.trajectory.state_names) + +# _the `solve!` proceeds as in the [Quantum State Sampling Problem](#Quantum-State-Sampling-Problem)]_ + +# ----- + +#= +## Unitary Variational Problem + +```@docs; canonical = false +UnitaryVariationalProblem +``` + +The `UnitaryVariationalProblem` uses a `VariationalQuantumSystem` to find a control that is +sensitive or robust to terms in the Hamiltonian. See the documentation for the +`VariationalQuantumSystem` in [`PiccoloQuantumObjects.jl`](https://github.com/harmoniqs/PiccoloQuantumObjects.jl) +for more details. +=# + +# _create a variational system, with a variational Hamiltonian, `PAULIS.X`_ +H_var = PAULIS.X +varsys = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [H_var]); + +# _create a variational problem that is robust to `PAULIS.X` at the end_ +robprob = UnitaryVariationalProblem(varsys, U_goal, T, Δt, robust_times=[[T]]); + +# ----- \ No newline at end of file diff --git a/docs/make.jl b/docs/make.jl index 7323c6e8..6d3e26a4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -5,24 +5,13 @@ using Literate push!(LOAD_PATH, joinpath(@__DIR__, "..", "src")) @info "Building Documenter site for QuantumCollocation.jl" -open(joinpath(@__DIR__, "src", "index.md"), write = true) do io - for line in eachline(joinpath(@__DIR__, "..", "README.md")) - if occursin("", line) - comment_content = match(r"", line).captures[1] - write(io, comment_content * "\n") - else - write(io, line * "\n") - end - end -end -# TODO: Callbacks are currently broken pages = [ "Home" => "index.md", - "Manual" => [ - "Problem Templates" => "generated/man/problem_templates.md", + "Library" => [ + "Ket Problem Templates" => "generated/man/ket_problem_templates.md", + "Unitary Problem Templates" => "generated/man/unitary_problem_templates.md", ], - "Library" => "lib.md", ] format = Documenter.HTML(; diff --git a/docs/src/assets/animation.gif b/docs/src/assets/animation.gif deleted file mode 100644 index 69e1adf7..00000000 Binary files a/docs/src/assets/animation.gif and /dev/null differ diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg deleted file mode 100644 index 9ee565d5..00000000 --- a/docs/src/assets/logo.svg +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/docs/src/assets/piccolo_logo.svg b/docs/src/assets/piccolo_logo.svg deleted file mode 100644 index d91bee46..00000000 --- a/docs/src/assets/piccolo_logo.svg +++ /dev/null @@ -1 +0,0 @@ - \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 83d0ef82..b548075a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,60 +1,3 @@ -```@raw html -
- - Piccolo.jl - -
- -
- - - - - - - -
- Documentation -
- - Stable - - - Dev - - - arXiv - - 		- Build Status -
- - Build Status - - - Coverage - - 		- License -
- - MIT License - - 		- Support -
- - Unitary Fund - -
-
- -
- Quickly set up and solve problem templates for quantum optimal control -
-
-``` - # QuantumCollocation.jl **QuantumCollocation.jl** sets up and solves *quantum control problems* as nonlinear programs (NLPs). In this context, a generic quantum control problem looks like @@ -67,7 +10,29 @@ ``` where $\mathbf{Z}$ is a trajectory containing states and controls, from [NamedTrajectories.jl](https://github.com/harmoniqs/NamedTrajectories.jl). -### Problem Templates +----- + +We provide a number of **problem templates** for making it simple and easy to set up and solve +certain types of quantum optimal control problems. These templates all construct a +`DirectTrajOptProblem` object from [DirectTrajOpt.jl](https://github.com/harmoniqs/DirectTrajOpt.jl), which stores all the parts of the optimal control problem. + +----- + +### Get started + +The problem templates are broken down by the state variable of the problem being solved. + +Ket Problem Templates: +- [Quantum State Smooth Pulse Problem](@ref) +- [Quantum State Minimum Time Problem](@ref) +- [Quantum State Sampling Problem](@ref) + +Unitary Problem Templates: +- [Unitary Smooth Pulse Problem](@ref) +- [Unitary Minimum Time Problem](@ref) +- [Unitary Sampling Problem](@ref) + +### Background *Problem Templates* are reusable design patterns for setting up and solving common quantum control problems. @@ -88,29 +53,10 @@ while a *UnitaryMinimumTimeProblem* minimizes time and constrains fidelity, \end{aligned} ``` -In each case, the dynamics between *knot points* $(U_t, a_t)$ and $(U_{t+1}, a_{t+1})$ are enforced as constraints on the states, which are free variables in the solver; this optimization framework is called *direct collocation*. For details of our implementation please see our award-winning IEEE QCE 2023 paper, [Direct Collocation for Quantum Optimal Control](https://arxiv.org/abs/2305.03261). If you use QuantumCollocation.jl in your work, please cite :raised_hands:! - -Problem templates give the user the ability to add other constraints and objective functions to this problem and solve it efficiently using [Ipopt.jl](https://github.com/jump-dev/Ipopt.jl) and [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) under the hood. +----- -## Installation - -This package is registered! To install, enter the Julia REPL, type `]` to enter pkg mode, and then run: -```julia -pkg> add QuantumCollocation -``` +In each case, the dynamics between *knot points* $(U_t, a_t)$ and $(U_{t+1}, a_{t+1})$ are enforced as constraints on the states, which are free variables in the solver; this optimization framework is called *direct trajectory optimization*. -## Example +----- -### Single Qubit Hadamard Gate -```Julia -using QuantumCollocation - -T = 50 -Δt = 0.2 -system = QuantumSystem([PAULIS[:X], PAULIS[:Y]]) -U_goal = GATES.H - -# Hadamard Gate -prob = UnitarySmoothPulseProblem(system, U_goal, T, Δt) -solve!(prob, max_iter=100) -``` +Problem templates give the user the ability to add other constraints and objective functions to this problem and solve it efficiently using [Ipopt.jl](https://github.com/jump-dev/Ipopt.jl) and [MathOptInterface.jl](https://github.com/jump-dev/MathOptInterface.jl) under the hood. diff --git a/docs/src/lib.md b/docs/src/lib.md index dc0310d9..11d9e133 100644 --- a/docs/src/lib.md +++ b/docs/src/lib.md @@ -18,9 +18,4 @@ Modules = [QuantumCollocation.TrajectoryInitialization] ## Trajectory Interpolations ```@autodocs Modules = [QuantumCollocation.TrajectoryInterpolations] -``` - -## Quantum System Templates -```@autodocs -Modules = [QuantumCollocation.QuantumSystemTemplates] ``` \ No newline at end of file diff --git a/src/problem_templates/_problem_templates.jl b/src/problem_templates/_problem_templates.jl index fd6e512d..89d91eab 100644 --- a/src/problem_templates/_problem_templates.jl +++ b/src/problem_templates/_problem_templates.jl @@ -17,6 +17,7 @@ using JLD2 using TestItems include("unitary_smooth_pulse_problem.jl") +include("unitary_variational_problem.jl") include("unitary_minimum_time_problem.jl") include("unitary_sampling_problem.jl") diff --git a/src/problem_templates/quantum_state_sampling_problem.jl b/src/problem_templates/quantum_state_sampling_problem.jl index 0859e3d1..9cc59a2c 100644 --- a/src/problem_templates/quantum_state_sampling_problem.jl +++ b/src/problem_templates/quantum_state_sampling_problem.jl @@ -73,7 +73,7 @@ function QuantumStateSamplingProblem( ) end - traj = merge(trajs, merge_names=(; a=1, da=1, dda=1, Δt=1),free_time=true) + traj = merge(trajs, merge_names=(; a=1, da=1, dda=1, Δt=1), free_time=true) end control_names = [ @@ -88,7 +88,7 @@ function QuantumStateSamplingProblem( for (weight, names) in zip(system_weights, state_names) for name in names - J += KetInfidelityLoss(name, traj, Q=weight * Q) + J += KetInfidelityObjective(name, traj, Q=weight * Q) end end diff --git a/src/problem_templates/quantum_state_smooth_pulse_problem.jl b/src/problem_templates/quantum_state_smooth_pulse_problem.jl index f5d9798a..6cd18912 100644 --- a/src/problem_templates/quantum_state_smooth_pulse_problem.jl +++ b/src/problem_templates/quantum_state_smooth_pulse_problem.jl @@ -128,7 +128,7 @@ function QuantumStateSmoothPulseProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) for name ∈ state_names - J += KetInfidelityLoss(name, traj; Q=Q) + J += KetInfidelityObjective(name, traj; Q=Q) end # Optional Piccolo constraints and objectives diff --git a/src/problem_templates/unitary_sampling_problem.jl b/src/problem_templates/unitary_sampling_problem.jl index 7e334079..f89451da 100644 --- a/src/problem_templates/unitary_sampling_problem.jl +++ b/src/problem_templates/unitary_sampling_problem.jl @@ -120,7 +120,7 @@ function UnitarySamplingProblem( J += QuadraticRegularizer(control_names[3], traj, R_dda) for (weight, op, name) in zip(system_weights, operators, state_names) - J += UnitaryInfidelityLoss(op, name, traj; Q=weight * Q) + J += UnitaryInfidelityObjective(op, name, traj; Q=weight * Q) end # Optional Piccolo constraints and objectives diff --git a/src/problem_templates/unitary_smooth_pulse_problem.jl b/src/problem_templates/unitary_smooth_pulse_problem.jl index 107d0de4..2a560f44 100644 --- a/src/problem_templates/unitary_smooth_pulse_problem.jl +++ b/src/problem_templates/unitary_smooth_pulse_problem.jl @@ -122,7 +122,7 @@ function UnitarySmoothPulseProblem( end # Objective - J = UnitaryInfidelityLoss(goal, state_name, traj; Q=Q) + J = UnitaryInfidelityObjective(goal, state_name, traj; Q=Q) control_names = [ name for name ∈ traj.names diff --git a/src/problem_templates/unitary_variational_problem.jl b/src/problem_templates/unitary_variational_problem.jl new file mode 100644 index 00000000..8bc013b7 --- /dev/null +++ b/src/problem_templates/unitary_variational_problem.jl @@ -0,0 +1,217 @@ +export UnitaryVariationalProblem + + +""" + UnitaryVariationalProblem( + system::VariationalQuantumSystem, + goal::AbstractPiccoloOperator, + T::Int, + Δt::Union{Float64, <:AbstractVector{Float64}}; + robust_times::AbstractVector{<:Union{AbstractVector{Int}, Nothing}}=[nothing for s ∈ system.G_vars], + sensitive_times::AbstractVector{<:Union{AbstractVector{Int}, Nothing}}=[nothing for s ∈ system.G_vars], + kwargs... + ) + +Constructs a unitary variational problem for optimizing quantum control trajectories. + +# Arguments +- `system::VariationalQuantumSystem`: The quantum system to be controlled, containing variational parameters. +- `goal::AbstractPiccoloOperator`: The target operator or state to achieve at the end of the trajectory. +- `T::Int`: The total number of timesteps in the trajectory. +- `Δt::Union{Float64, <:AbstractVector{Float64}}`: The timestep duration or a vector of timestep durations. +- `robust_times::AbstractVector`: Times at which robustness to variations in the trajectory is enforced. +- `sensitive_times::AbstractVector`: Times at which sensitivity to variations in the trajectory is enhanced. +- `unitary_integrator`: The integrator used for unitary evolution (default: `VariationalUnitaryIntegrator`). +- `state_name::Symbol`: The name of the state variable in the trajectory (default: `:Ũ⃗`). +- `variational_state_name::Symbol`: The name of the variational state variable (default: `:Ũ⃗ₐ`). +- `variational_scales::AbstractVector`: Scaling factors for the variational state variables (default: `1.0`). +- `control_name::Symbol`: The name of the control variable (default: `:a`). +- `timestep_name::Symbol`: The name of the timestep variable (default: `:Δt`). +- `init_trajectory::Union{NamedTrajectory, Nothing}`: An optional initial trajectory to start optimization. +- `a_bound::Float64`: The bound for the control variable `a` (default: `1.0`). +- `a_bounds::Vector`: Bounds for each control variable (default: filled with `a_bound`). +- `da_bound::Float64`: The bound for the derivative of the control variable (default: `Inf`). +- `da_bounds::Vector`: Bounds for each derivative of the control variable. +- `dda_bound::Float64`: The bound for the second derivative of the control variable (default: `1.0`). +- `dda_bounds::Vector`: Bounds for each second derivative of the control variable. +- `Δt_min::Float64`: Minimum allowed timestep duration. +- `Δt_max::Float64`: Maximum allowed timestep duration. +- `Q::Float64`: Weight for the unitary infidelity objective (default: `100.0`). +- `Q_v::Float64`: Weight for sensitivity objectives (default: `1.0`). +- `R`: Regularization weight for control variables (default: `1e-2`). +- `R_a`, `R_da`, `R_dda`: Regularization weights for control, its derivative, and second derivative. +- `constraints::Vector`: Additional constraints for the optimization problem. +- `piccolo_options::PiccoloOptions`: Options for configuring the Piccolo optimization framework. + +# Returns +A `DirectTrajOptProblem` object representing the optimization problem, including the +trajectory, objective, integrators, and constraints. + +# Notes +This function constructs a trajectory optimization problem for quantum control using +variational principles. It supports robust and sensitive trajectory design, regularization, +and optional constraints. The problem is solved using the Piccolo optimization framework. + +""" +function UnitaryVariationalProblem end + +function UnitaryVariationalProblem( + system::VariationalQuantumSystem, + goal::AbstractPiccoloOperator, + T::Int, + Δt::Union{Float64, <:AbstractVector{Float64}}; + robust_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars], + sensitive_times::AbstractVector{<:AbstractVector{Int}}=[Int[] for s ∈ system.G_vars], + variational_integrator=VariationalUnitaryIntegrator, + variational_scales::AbstractVector{<:Float64}=fill(1.0, length(system.G_vars)), + state_name::Symbol = :Ũ⃗, + variational_state_name::Symbol = :Ũ⃗ᵥ, + control_name::Symbol = :a, + timestep_name::Symbol = :Δt, + init_trajectory::Union{NamedTrajectory, Nothing}=nothing, + a_bound::Float64=1.0, + a_bounds=fill(a_bound, system.n_drives), + da_bound::Float64=Inf, + da_bounds::Vector{Float64}=fill(da_bound, system.n_drives), + dda_bound::Float64=1.0, + dda_bounds::Vector{Float64}=fill(dda_bound, system.n_drives), + Δt_min::Float64=Δt isa Float64 ? 0.5 * Δt : 0.5 * mean(Δt), + Δt_max::Float64=Δt isa Float64 ? 1.5 * Δt : 1.5 * mean(Δt), + Q::Float64=100.0, + Q_s::Float64=1e-2, + Q_r::Float64=100.0, + R=1e-2, + R_a::Union{Float64, Vector{Float64}}=R, + R_da::Union{Float64, Vector{Float64}}=R, + R_dda::Union{Float64, Vector{Float64}}=R, + constraints::Vector{<:AbstractConstraint}=AbstractConstraint[], + piccolo_options::PiccoloOptions=PiccoloOptions(), +) + if piccolo_options.verbose + println(" constructing UnitaryVariationalProblem...") + println("\tusing integrator: $(typeof(variational_integrator))") + println("\ttotal variational parameters: $(length(system.G_vars))") + end + + # Trajectory + if !isnothing(init_trajectory) + traj = deepcopy(init_trajectory) + else + traj = initialize_trajectory( + goal, + T, + Δt, + system.n_drives, + (a_bounds, da_bounds, dda_bounds); + state_name=state_name, + control_name=control_name, + timestep_name=timestep_name, + Δt_bounds=(Δt_min, Δt_max), + zero_initial_and_final_derivative=piccolo_options.zero_initial_and_final_derivative, + geodesic=piccolo_options.geodesic, + bound_state=piccolo_options.bound_state, + rollout_integrator=piccolo_options.rollout_integrator, + verbose=piccolo_options.verbose + ) + end + + _, Ũ⃗_vars = variational_unitary_rollout( + traj, + system; + unitary_name=state_name, + drive_name=control_name + ) + + variational_state_names = [ + Symbol(string(variational_state_name) * "$(i)") for i in eachindex(system.G_vars) + ] + + for (name, scale, Ũ⃗_v) in zip(variational_state_names, variational_scales, Ũ⃗_vars) + add_component!(traj, name, Ũ⃗_v / scale; type=:state) + traj.initial = merge(traj.initial, (name => Ũ⃗_v[:, 1] / scale,)) + end + + control_names = [ + name for name ∈ traj.names + if endswith(string(name), string(control_name)) + ] + + # objective + J = UnitaryInfidelityObjective(goal, state_name, traj; Q=Q) + J += QuadraticRegularizer(control_names[1], traj, R_a) + J += QuadraticRegularizer(control_names[2], traj, R_da) + J += QuadraticRegularizer(control_names[3], traj, R_dda) + + # sensitivity + for (name, scale, s, r) ∈ zip( + variational_state_names, + variational_scales, + sensitive_times, + robust_times + ) + @assert isdisjoint(s, r) + J += UnitarySensitivityObjective( + name, + traj, + [s; r]; + Qs=[fill(-Q_s, length(s)); fill(Q_r, length(r))], + scale=scale + ) + end + + # Optional Piccolo constraints and objectives + apply_piccolo_options!( + J, constraints, piccolo_options, traj, state_name, timestep_name; + state_leakage_indices=goal isa EmbeddedOperator ? get_leakage_indices(goal) : nothing + ) + + integrators = [ + variational_integrator( + system, traj, state_name, variational_state_names, control_name, + scales=variational_scales + ), + DerivativeIntegrator(traj, control_name, control_names[2]), + DerivativeIntegrator(traj, control_names[2], control_names[3]), + ] + + return DirectTrajOptProblem( + traj, + J, + integrators; + constraints=constraints + ) +end + +# *************************************************************************** # + +@testitem "Sensitive and robust" begin + using LinearAlgebra + using PiccoloQuantumObjects + + system = QuantumSystem([PAULIS.X, PAULIS.Y]) + varsys = VariationalQuantumSystem([PAULIS.X, PAULIS.Y], [PAULIS.X] ) + T = 50 + Δt = 0.2 + + sense_scale = 8.0 + sense_prob = UnitaryVariationalProblem( + varsys, GATES.X, T, Δt, + variational_scales=[sense_scale], + sensitive_times=[[T]], + piccolo_options=PiccoloOptions(verbose=false) + ) + solve!(sense_prob, max_iter=20, print_level=1, verbose=false) + + rob_scale = 1 / 8.0 + rob_prob = UnitaryVariationalProblem( + varsys, GATES.X, T, Δt, + variational_scales=[rob_scale], + robust_times=[[T]], + piccolo_options=PiccoloOptions(verbose=false) + ) + solve!(rob_prob, max_iter=20, print_level=1, verbose=false) + + sense_n = norm(sense_scale * sense_prob.trajectory.Ũ⃗ᵥ1) + rob_n = norm(rob_scale * rob_prob.trajectory.Ũ⃗ᵥ1) + @assert sense_n > rob_n +end \ No newline at end of file diff --git a/src/quantum_constraints.jl b/src/quantum_constraints.jl index 7e3a988b..2e77dafd 100644 --- a/src/quantum_constraints.jl +++ b/src/quantum_constraints.jl @@ -21,7 +21,7 @@ function FinalKetFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = ψ̃ -> [ - ket_infidelity_loss(ψ̃, ψ_goal) - abs(1 - final_fidelity) + abs(QuantumObjectives.ket_fidelity_loss(ψ̃, ψ_goal) - final_fidelity) ] return NonlinearKnotPointConstraint( @@ -44,7 +44,7 @@ function FinalUnitaryFidelityConstraint( traj::NamedTrajectory ) terminal_constraint = Ũ⃗ -> [ - unitary_infidelity_loss(Ũ⃗, U_goal) - abs(1 - final_fidelity) + abs(QuantumObjectives.unitary_fidelity_loss(Ũ⃗, U_goal) - final_fidelity) ] return NonlinearKnotPointConstraint( diff --git a/src/quantum_integrators.jl b/src/quantum_integrators.jl index 4c8e144a..8108a1c8 100644 --- a/src/quantum_integrators.jl +++ b/src/quantum_integrators.jl @@ -3,14 +3,20 @@ module QuantumIntegrators export KetIntegrator export UnitaryIntegrator export DensityMatrixIntegrator +export VariationalUnitaryIntegrator using LinearAlgebra using NamedTrajectories using DirectTrajOpt using PiccoloQuantumObjects +using SparseArrays const ⊗ = kron +# ----------------------------------------------------------------------------- # +# Default Integrators +# ----------------------------------------------------------------------------- # + function KetIntegrator( sys::QuantumSystem, traj::NamedTrajectory, @@ -39,4 +45,40 @@ function DensityMatrixIntegrator( return BilinearIntegrator(sys.𝒢, traj, ρ̃, a) end +# ----------------------------------------------------------------------------- # +# Variational Integrators +# ----------------------------------------------------------------------------- # + +function VariationalKetIntegrator( + sys::VariationalQuantumSystem, + traj::NamedTrajectory, + ψ̃::Symbol, + ψ̃_variations::AbstractVector{Symbol}, + a::Symbol; + scale::Float64=1.0, +) + var_ψ̃ = vcat(ψ̃, ψ̃_variations...) + G = a -> Isomorphisms.var_G(sys.G(a), [G(a) / scale for G in sys.G_vars]) + return BilinearIntegrator(G, traj, var_ψ̃, a) +end + +function VariationalUnitaryIntegrator( + sys::VariationalQuantumSystem, + traj::NamedTrajectory, + Ũ⃗::Symbol, + Ũ⃗_variations::AbstractVector{Symbol}, + a::Symbol; + scales::AbstractVector{<:Float64}=fill(1.0, length(sys.G_vars)), +) + var_Ũ⃗ = vcat(Ũ⃗, Ũ⃗_variations...) + + function Ĝ(a) + G0 = sys.G(a) + Gs = typeof(G0)[I(sys.levels) ⊗ G(a) / scale for (scale, G) in zip(scales, sys.G_vars)] + return Isomorphisms.var_G(I(sys.levels) ⊗ G0, Gs) + end + return BilinearIntegrator(Ĝ, traj, var_Ũ⃗, a) +end + + end \ No newline at end of file diff --git a/src/quantum_objectives.jl b/src/quantum_objectives.jl index 7c085bdd..21ed1cc8 100644 --- a/src/quantum_objectives.jl +++ b/src/quantum_objectives.jl @@ -1,11 +1,9 @@ module QuantumObjectives -export ket_infidelity_loss -export unitary_infidelity_loss -export density_matrix_pure_state_infidelity_loss -export KetInfidelityLoss -export UnitaryInfidelityLoss -export DensityMatrixPureStateInfidelityLoss +export KetInfidelityObjective +export UnitaryInfidelityObjective +export DensityMatrixPureStateInfidelityObjective +export UnitarySensitivityObjective using LinearAlgebra using NamedTrajectories @@ -16,23 +14,22 @@ using DirectTrajOpt # Kets # --------------------------------------------------------- -function ket_infidelity_loss( +function ket_fidelity_loss( ψ̃::AbstractVector, ψ_goal::AbstractVector{<:Complex{Float64}} ) ψ = iso_to_ket(ψ̃) - ℱ = abs2(ψ_goal' * ψ) - return abs(1 - ℱ) + return abs2(ψ_goal' * ψ) end -function KetInfidelityLoss( +function KetInfidelityObjective( ψ̃_name::Symbol, traj::NamedTrajectory; Q=100.0 ) ψ_goal = iso_to_ket(traj.goal[ψ̃_name]) - ℓ = ψ̃ -> ket_infidelity_loss(ψ̃, ψ_goal) - return TerminalLoss(ℓ, ψ̃_name, traj; Q=Q) + ℓ = ψ̃ -> abs(1 - ket_fidelity_loss(ψ̃, ψ_goal)) + return TerminalObjective(ℓ, ψ̃_name, traj; Q=Q) end @@ -40,36 +37,34 @@ end # Unitaries # --------------------------------------------------------- -function unitary_infidelity_loss( - Ũ⃗::AbstractVector, +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real}, U_goal::AbstractMatrix{<:Complex{Float64}} ) U = iso_vec_to_operator(Ũ⃗) n = size(U, 1) - ℱ = abs2(tr(U_goal' * U)) / n^2 - return abs(1 - ℱ) + return abs2(tr(U_goal' * U)) / n^2 end -function unitary_infidelity_loss( - Ũ⃗::AbstractVector, +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real}, op::EmbeddedOperator ) U_goal = unembed(op) U = iso_vec_to_operator(Ũ⃗)[op.subspace, op.subspace] n = length(op.subspace) M = U_goal'U - ℱ = 1 / (n * (n + 1)) * (abs(tr(M'M)) + abs2(tr(M))) - return abs(1 - ℱ) + return 1 / (n * (n + 1)) * (abs(tr(M'M)) + abs2(tr(M))) end -function UnitaryInfidelityLoss( +function UnitaryInfidelityObjective( U_goal::AbstractPiccoloOperator, Ũ⃗_name::Symbol, traj::NamedTrajectory; Q=100.0 ) - ℓ = Ũ⃗ -> unitary_infidelity_loss(Ũ⃗, U_goal) - return TerminalLoss(ℓ, Ũ⃗_name, traj; Q=Q) + ℓ = Ũ⃗ -> abs(1 - unitary_fidelity_loss(Ũ⃗, U_goal)) + return TerminalObjective(ℓ, Ũ⃗_name, traj; Q=Q) end # --------------------------------------------------------- @@ -85,14 +80,45 @@ function density_matrix_pure_state_infidelity_loss( return abs(1 - ℱ) end -function DensityMatrixPureStateInfidelityLoss( +function DensityMatrixPureStateInfidelityObjective( ρ̃_name::Symbol, ψ_goal::AbstractVector{<:Complex{Float64}}, traj::NamedTrajectory; Q=100.0 ) ℓ = ρ̃ -> density_matrix_pure_state_infidelity_loss(ρ̃, ψ_goal) - return TerminalLoss(ℓ, ρ̃_name, traj; Q=Q) + return TerminalObjective(ℓ, ρ̃_name, traj; Q=Q) +end + +# --------------------------------------------------------- +# Sensitivity +# --------------------------------------------------------- + +function unitary_fidelity_loss( + Ũ⃗::AbstractVector{<:Real} +) + U = iso_vec_to_operator(Ũ⃗) + n = size(U, 1) + return abs2(tr(U' * U)) / n^2 end +function UnitarySensitivityObjective( + name::Symbol, + traj::NamedTrajectory, + times::AbstractVector{Int}; + Qs::AbstractVector{<:Float64}=fill(1.0, length(times)), + scale::Float64=1.0, +) + ℓ = Ũ⃗ -> scale^4 * unitary_fidelity_loss(Ũ⃗) + + return KnotPointObjective( + ℓ, + name, + traj; + Qs=Qs, + times=times + ) +end + + end \ No newline at end of file diff --git a/src/trajectory_initialization.jl b/src/trajectory_initialization.jl index 3c329801..7ab0cdd3 100644 --- a/src/trajectory_initialization.jl +++ b/src/trajectory_initialization.jl @@ -424,6 +424,7 @@ function initialize_trajectory( # Construct phase data phase_data = isnothing(phase_operators) ? nothing : π * randn(length(phase_operators)) + return initialize_trajectory( [Ũ⃗_traj], [Ũ⃗_init], @@ -438,6 +439,8 @@ function initialize_trajectory( ) end + + """ initialize_trajectory