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 an open quantum system example, with a single qubit subject to an amplitude damping channel. The target evolution is the Hadamard gate. The user can experiment with the strength of the amplitude damping by changing the gamma variable value
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 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
example_name = 'Lindblad'
Sx = sigmax()
Sy = sigmay()
Sz = sigmaz()
Si = identity(2)
Sd = Qobj(np.array([[0, 1],
[0, 0]]))
Sm = Qobj(np.array([[0, 0],
[1, 0]]))
Sd_m = Qobj(np.array([[1, 0],
[0, 0]]))
Sm_d = Qobj(np.array([[0, 0],
[0, 1]]))
#Amplitude damping#
#Damping rate:
gamma = 0.1
L0_Ad = gamma*(2*tensor(Sm, Sd.trans()) -
(tensor(Sd_m, Si) + tensor(Si, Sd_m.trans())))
#sigma X control
LC_x = -1j*(tensor(Sx, Si) - tensor(Si, Sx))
#sigma Y control
LC_y = -1j*(tensor(Sy, Si) - tensor(Si, Sy.trans()))
#sigma Z control
LC_z = -1j*(tensor(Sz, Si) - tensor(Si, Sz))
#Drift
drift = L0_Ad
#Controls
ctrls = [LC_z, LC_x]
# Number of ctrls
n_ctrls = len(ctrls)
initial = identity(4)
#Target
#Hadamard gate
had_gate = hadamard_transform(1)
target_DP = tensor(had_gate, had_gate)
# Number of time slots
n_ts = 200
# Time allowed for the evolution
evo_time = 2
# Fidelity error target
fid_err_targ = 1e-10
# Maximum iterations for the optisation algorithm
max_iter = 200
# 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 = 'RND'
#Set to None to suppress output files
f_ext = "{}_n_ts{}_ptype{}.txt".format(example_name, n_ts, p_type)
# Note that this call will take the defaults
# dyn_type='GEN_MAT'
# This means that matrices that describe the dynamics are assumed to be
# general, i.e. the propagator can be calculated using:
# expm(combined_dynamics*dt)
# prop_type='FRECHET'
# and the propagators and their gradients will be calculated using the
# Frechet method, i.e. an exact gradent
# fid_type='TRACEDIFF'
# and 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(drift, ctrls, initial, target_DP, n_ts, evo_time,
fid_err_targ=fid_err_targ, min_grad=min_grad,
max_iter=max_iter, max_wall_time=max_wall_time,
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: [[-0.2+0.j 0.0+0.j 0.0+0.j 0.0+0.j] [ 0.0+0.j -0.1+0.j 0.0+0.j 0.0+0.j] [ 0.0+0.j 0.0+0.j -0.1+0.j 0.0+0.j] [ 0.2+0.j 0.0+0.j 0.0+0.j 0.0+0.j]] Control 1 dynamics generator: [[ 0.+0.j 0.+0.j 0.+0.j 0.+0.j] [ 0.+0.j 0.-2.j 0.+0.j 0.+0.j] [ 0.+0.j 0.+0.j 0.+2.j 0.+0.j] [ 0.+0.j 0.+0.j 0.+0.j 0.+0.j]] Control 2 dynamics generator: [[ 0.+0.j 0.+1.j 0.-1.j 0.+0.j] [ 0.+1.j 0.+0.j 0.+0.j 0.-1.j] [ 0.-1.j 0.+0.j 0.+0.j 0.+1.j] [ 0.+0.j 0.-1.j 0.+1.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.5+0.j 0.5+0.j 0.5+0.j 0.5+0.j] [ 0.5+0.j -0.5+0.j 0.5+0.j -0.5+0.j] [ 0.5+0.j 0.5+0.j -0.5+0.j -0.5+0.j] [ 0.5+0.j -0.5+0.j -0.5+0.j 0.5+0.j]] INFO:qutip.control.pulseoptim:Initial amplitudes output to file: ctrl_amps_initial_Lindblad_n_ts200_ptypeRND.txt INFO:qutip.control.optimizer:Optimising pulse using L-BFGS-B INFO:qutip.control.pulseoptim:Final amplitudes output to file: ctrl_amps_final_Lindblad_n_ts200_ptypeRND.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:30.118382 Wall time computing Hamiltonians: 0:00:00.626220 (2.08%) Wall time computing propagators: 0:00:26.829948 (89.08%) Wall time computing forward propagation: 0:00:00.098895 (0.33%) Wall time computing onward propagation: 0:00:00.096445 (0.32%) Wall time computing gradient: 0:00:02.407970 (8.00%) **** Iterations and function calls **** Number of iterations: 128 Number of fidelity function calls: 168 Number of times fidelity is computed: 168 Number of gradient function calls: 168 Number of times gradients are computed: 168 Number of times timeslot evolution is recomputed: 168 **** Control amplitudes **** Number of control amplitude updates: 167 Mean number of updates per iteration: 1.3046875 Number of timeslot values changed: 33400 Mean number of timeslot changes per update: 200.0 Number of amplitude values changed: 66800 Mean number of amplitude changes per update: 400.0 ------------------------------------ Final evolution Quantum object: dims = [[4], [4]], shape = [4, 4], type = oper, isherm = False Qobj data = [[ 0.49311703 -1.03216105e-16j 0.37516863 -3.77560294e-03j 0.37516863 +3.77560294e-03j 0.50421905 -5.03070375e-17j] [ 0.38301167 -8.20510888e-03j -0.39243479 -1.28423923e-02j 0.38198219 +5.64382365e-03j -0.39095475 +1.11005674e-02j] [ 0.38301167 +8.20510888e-03j 0.38198219 -5.64382365e-03j -0.39243479 +1.28423923e-02j -0.39095475 -1.11005674e-02j] [ 0.50688297 +1.03215974e-16j -0.37516863 +3.77560294e-03j -0.37516863 -3.77560294e-03j 0.49578095 +3.64291220e-17j]] ********* Summary ***************** Final fidelity error 0.02068059184705504 Final gradient normal 6.0807607767135144e-05 Terminated due to Max wall time exceeded Number of iterations 128 Completed in 0:00:30.118382 HH:MM:SS.US
t = result.time[:n_ts]
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):
amps = result.initial_amps[:, j]
ax1.plot(t, amps)
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):
amps = result.final_amps[:, j]
ax2.plot(t, amps)
plt.show()
from qutip.ipynbtools import version_table
version_table()
Software | Version |
---|---|
QuTiP | 3.1.0 |
Numpy | 1.9.1 |
matplotlib | 1.4.2 |
Python | 3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2] |
IPython | 2.3.1 |
Cython | 0.21.2 |
OS | posix [linux] |
SciPy | 0.14.1 |
Tue Jan 13 13:32:19 2015 JST |