# Calculation of control fields for symplectic dynamics using L-BFGS-B algorithm¶

Alexander Pitchford ([email protected])

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.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import datetime

In [2]:
from qutip import Qobj, identity, sigmax, sigmay, sigmaz, tensor
import qutip.logging 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'


### Defining the physics¶

In [3]:
#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()


### Defining the time evolution parameters¶

In [4]:
# Number of time slots
n_ts = 1000
# Time allowed for the evolution
evo_time = 10


### Set the conditions which will cause the pulse optimisation to terminate¶

In [5]:
# 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
# as this tends to 0 -> local minima has been found


### Set the initial pulse type¶

In [6]:
# pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|
p_type = 'ZERO'


### Give an extension for output files¶

In [7]:
#Set to None to suppress output files
f_ext = "{}_n_ts{}_ptype{}.txt".format(example_name, n_ts, p_type)


### Run the optimisation¶

In [8]:
# 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,
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.pulseoptim:System configuration:
Drift dynamics generator:
[[ 1.0+0.j  0.0+0.j  0.5+0.j  0.0+0.j]
[ 0.0+0.j  1.0+0.j  0.0+0.j  0.5+0.j]
[ 0.5+0.j  0.0+0.j  1.0+0.j  0.0+0.j]
[ 0.0+0.j  0.5+0.j  0.0+0.j  1.0+0.j]]
Control 1 dynamics generator:
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
[ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
[ 0.+0.j  0.+0.j  0.+0.j  0.+0.j]
[ 0.+0.j  0.+0.j  0.+0.j  0.+0.j]]
Initial operator:
[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]
[ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
[ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]
[ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]
Target operator:
[[ 0.54030231+0.j  0.00000000+0.j  0.00000000+0.j  0.84147098+0.j]
[ 0.00000000+0.j  0.54030231+0.j -0.84147098+0.j  0.00000000+0.j]
[ 0.00000000+0.j  0.84147098+0.j  0.54030231+0.j  0.00000000+0.j]
[-0.84147098+0.j  0.00000000+0.j  0.00000000+0.j  0.54030231+0.j]]
INFO:qutip.control.pulseoptim:Initial amplitudes output to file: ctrl_amps_initial_Symplectic_n_ts1000_ptypeZERO.txt
INFO:qutip.control.optimizer:Optimising pulse using L-BFGS-B
INFO:qutip.control.pulseoptim:Final amplitudes output to file: ctrl_amps_final_Symplectic_n_ts1000_ptypeZERO.txt


### Report the results¶

In [10]:
result.stats.report()
print("Final evolution\n{}\n".format(result.evo_full_final))
print("********* Summary *****************")
print("Final fidelity error {}".format(result.fid_err))
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:08.261914
Wall time computing Hamiltonians: 0:00:00.208757 (2.53%)
Wall time computing propagators: 0:00:07.326453 (88.68%)
Wall time computing forward propagation: 0:00:00.052961 (0.64%)
Wall time computing onward propagation: 0:00:00.051577 (0.62%)
Wall time computing gradient: 0:00:00.611217 (7.40%)

**** Iterations and function calls ****
Number of iterations: 13
Number of fidelity function calls: 18
Number of times fidelity is computed: 18
Number of gradient function calls: 17
Number of times gradients are computed: 17
Number of times timeslot evolution is recomputed: 18

**** Control amplitudes ****
Number of control amplitude updates: 17
Mean number of updates per iteration: 1.3076923076923077
Number of timeslot values changed: 17000
Mean number of timeslot changes per update: 1000.0
Number of amplitude values changed: 17000
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.40302194e-01  -3.23386879e-07   9.59430701e-08   8.41471056e-01]
[  3.23386880e-07   5.40302194e-01  -8.41471056e-01   9.59430701e-08]
[  9.59430675e-08   8.41471056e-01   5.40302194e-01   2.00178239e-07]
[ -8.41471056e-01   9.59430665e-08  -2.00178237e-07   5.40302194e-01]]

********* Summary *****************
Final fidelity error 4.955449758551298e-14
Terminated due to Goal achieved
Number of iterations 13
Completed in 0:00:08.261914 HH:MM:SS.US


### Plot the initial and final amplitudes¶

In [11]:
t = result.time[:n_ts]

fig1 = plt.figure()
ax1.set_title("Initial Control amps")
ax1.set_xlabel("Time")
ax1.set_ylabel("Control amplitude")
for j in range(n_ctrls):
amps = result.initial_amps[:, j]
ax1.plot(t, amps)

ax2.set_title("Optimised Control Amplitudes")
ax2.set_xlabel("Time")
ax2.set_ylabel("Control amplitude")
for j in range(n_ctrls):
amps = result.final_amps[:, j]
ax2.plot(t, amps)

plt.show()


### Versions¶

In [12]:
from qutip.ipynbtools import version_table

version_table()

Out[12]:
SoftwareVersion
OSposix [linux]
IPython2.3.1
matplotlib1.4.2
Cython0.21.2
QuTiP3.1.0
Python3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2]
SciPy0.14.1
Numpy1.9.1
Tue Jan 13 13:35:23 2015 JST