Alexander Pitchford (agp1@aber.ac.uk)
Example to demonstrate using the control library to determine control pulses using the ctrlpulseoptim.optimize_pulse function. The (default) L-BFGS-B algorithm is used to optimise the pulse to minimise the fidelity error, which in this case is given by the 'Trace difference' norm.
This in a Symplectic quantum system example, with two coupled oscillators
The user can experiment with the timeslicing, by means of changing the number of timeslots and/or total time for the evolution. Different initial (starting) pulse types can be tried. The initial and final pulses are displayed in a plot
This example assumes that the example-control-pulseoptim-Hadamard has already been tried, and hence explanations in that notebook are not repeated here.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import datetime
from qutip import Qobj, identity, sigmax, sigmay, sigmaz, tensor
from qutip.qip import hadamard_transform
import qutip.logging_utils as logging
logger = logging.get_logger()
#Set this to None or logging.WARN for 'quiet' execution
log_level = logging.INFO
#QuTiP control modules
import qutip.control.pulseoptim as cpo
import qutip.control.symplectic as sympl
example_name = 'Symplectic'
#Drift
w1 = 1
w2 = 1
g1 = 0.5
A0 = Qobj(np.array([[w1, 0, g1, 0],
[0, w1, 0, g1],
[g1, 0, w2, 0],
[0, g1, 0, w2]]))
#Control
Ac = Qobj(np.array([[1, 0, 0, 0,], \
[0, 1, 0, 0], \
[0, 0, 0, 0], \
[0, 0, 0, 0]]))
ctrls = [Ac]
n_ctrls = len(ctrls)
initial = identity(4)
# Target
a = 1
Ag = np.array([[0, 0, a, 0],
[0, 0, 0, a],
[a, 0, 0, 0],
[0, a, 0, 0]])
Sg = Qobj(sympl.calc_omega(2).dot(Ag)).expm()
# Number of time slots
n_ts = 1000
# Time allowed for the evolution
evo_time = 10
# Fidelity error target
fid_err_targ = 1e-10
# Maximum iterations for the optisation algorithm
max_iter = 500
# Maximum (elapsed) time allowed in seconds
max_wall_time = 30
# Minimum gradient (sum of gradients squared)
# as this tends to 0 -> local minima has been found
min_grad = 1e-20
# pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|
p_type = 'ZERO'
#Set to None to suppress output files
f_ext = "{}_n_ts{}_ptype{}.txt".format(example_name, n_ts, p_type)
# Note that this call uses
# dyn_type='SYMPL'
# This means that matrices that describe the dynamics are assumed to be
# Symplectic, i.e. the propagator can be calculated using
# expm(combined_dynamics.omega*dt)
# This has defaults for:
# prop_type='FRECHET'
# therefore the propagators and their gradients will be calculated using the
# Frechet method, i.e. an exact gradient
# fid_type='TRACEDIFF'
# so that the fidelity error, i.e. distance from the target, is give
# by the trace of the difference between the target and evolved operators
result = cpo.optimize_pulse(A0, ctrls, initial, Sg, n_ts, evo_time,
fid_err_targ=fid_err_targ, min_grad=min_grad,
max_iter=max_iter, max_wall_time=max_wall_time,
dyn_type='SYMPL',
out_file_ext=f_ext, init_pulse_type=p_type,
log_level=log_level, gen_stats=True)
INFO:qutip.control.dynamics:Setting memory optimisations for level 0 INFO:qutip.control.dynamics:Internal operator data type choosen to be <class 'numpy.ndarray'> INFO:qutip.control.dynamics:phased dynamics generator caching True INFO:qutip.control.dynamics:propagator gradient caching True INFO:qutip.control.dynamics:eigenvector adjoint caching True INFO:qutip.control.dynamics:use sparse eigen decomp False INFO:qutip.control.pulseoptim:System configuration: Drift dynamics generator: Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = True Qobj data = [[ 1. 0. 0.5 0. ] [ 0. 1. 0. 0.5] [ 0.5 0. 1. 0. ] [ 0. 0.5 0. 1. ]] Control 1 dynamics generator: Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = True Qobj data = [[ 1. 0. 0. 0.] [ 0. 1. 0. 0.] [ 0. 0. 0. 0.] [ 0. 0. 0. 0.]] Initial state / operator: Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = True Qobj data = [[ 1. 0. 0. 0.] [ 0. 1. 0. 0.] [ 0. 0. 1. 0.] [ 0. 0. 0. 1.]] Target state / operator: Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False Qobj data = [[ 0.54030231 0. 0. 0.84147098] [ 0. 0.54030231 -0.84147098 0. ] [ 0. 0.84147098 0.54030231 0. ] [-0.84147098 0. 0. 0.54030231]] INFO:qutip.control.pulseoptim:Initial amplitudes output to file: ctrl_amps_initial_Symplectic_n_ts1000_ptypeZERO.txt INFO:qutip.control.optimizer:Optimising pulse(s) using GRAPE with 'fmin_l_bfgs_b' method INFO:qutip.control.pulseoptim:Final amplitudes output to file: ctrl_amps_final_Symplectic_n_ts1000_ptypeZERO.txt
result.stats.report()
print("Final evolution\n{}\n".format(result.evo_full_final))
print("********* Summary *****************")
print("Final fidelity error {}".format(result.fid_err))
print("Final gradient normal {}".format(result.grad_norm_final))
print("Terminated due to {}".format(result.termination_reason))
print("Number of iterations {}".format(result.num_iter))
print("Completed in {} HH:MM:SS.US".format(datetime.timedelta(seconds=result.wall_time)))
------------------------------------ ---- Control optimisation stats ---- **** Timings (HH:MM:SS.US) **** Total wall time elapsed during optimisation: 0:00:02.998303 Wall time computing Hamiltonians: 0:00:00.112098 (3.74%) Wall time computing propagators: 0:00:02.615654 (87.24%) Wall time computing forward propagation: 0:00:00.022306 (0.74%) Wall time computing onward propagation: 0:00:00.019775 (0.66%) Wall time computing gradient: 0:00:00.222838 (7.43%) **** Iterations and function calls **** Number of iterations: 8 Number of fidelity function calls: 14 Number of times fidelity is computed: 14 Number of gradient function calls: 13 Number of times gradients are computed: 13 Number of times timeslot evolution is recomputed: 14 **** Control amplitudes **** Number of control amplitude updates: 13 Mean number of updates per iteration: 1.625 Number of timeslot values changed: 13000 Mean number of timeslot changes per update: 1000.0 Number of amplitude values changed: 13000 Mean number of amplitude changes per update: 1000.0 ------------------------------------ Final evolution Quantum object: dims = [[4], [4]], shape = (4, 4), type = oper, isherm = False Qobj data = [[ 5.40298994e-01 -1.56965600e-06 -7.01444121e-06 8.41473111e-01] [ 1.56965600e-06 5.40298994e-01 -8.41473111e-01 -7.01444122e-06] [ -7.01444163e-06 8.41473111e-01 5.40298994e-01 1.05774201e-05] [ -8.41473111e-01 -7.01444163e-06 -1.05774201e-05 5.40298994e-01]] ********* Summary ***************** Final fidelity error 6.093073931117445e-11 Final gradient normal 8.705152868582485e-06 Terminated due to Goal achieved Number of iterations 8 Completed in 0:00:02.998303 HH:MM:SS.US
fig1 = plt.figure()
ax1 = fig1.add_subplot(2, 1, 1)
ax1.set_title("Initial Control amps")
ax1.set_xlabel("Time")
ax1.set_ylabel("Control amplitude")
for j in range(n_ctrls):
ax1.step(result.time,
np.hstack((result.initial_amps[:, j], result.initial_amps[-1, j])),
where='post')
ax2 = fig1.add_subplot(2, 1, 2)
ax2.set_title("Optimised Control Amplitudes")
ax2.set_xlabel("Time")
ax2.set_ylabel("Control amplitude")
for j in range(n_ctrls):
ax2.step(result.time,
np.hstack((result.final_amps[:, j], result.final_amps[-1, j])),
where='post')
plt.tight_layout()
plt.show()
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
QuTiP | 4.1.0 |
Numpy | 1.11.3 |
SciPy | 0.18.1 |
matplotlib | 2.0.0 |
Cython | 0.25.2 |
Number of CPUs | 4 |
BLAS Info | INTEL MKL |
IPython | 5.1.0 |
Python | 3.6.0 |Anaconda 4.3.1 (64-bit)| (default, Dec 23 2016, 12:22:00) [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] |
OS | posix [linux] |
Fri Jul 14 17:15:34 2017 BST |