Skip to content

Commit 2be5a07

Browse files
authored
Migrate and fix Optimal control tutorials from qutip-notebooks to qutip-tutorials (#122)
* migrate and fix opt-control tutorials * removed dependence on qutip-qip * change names * notebooks to markdown * change links
1 parent a8b6a13 commit 2be5a07

8 files changed

+1607
-11
lines changed

tutorials-v5/optimal-control/01-optimal-control-overview.md

+175
Large diffs are not rendered by default.
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,210 @@
1+
---
2+
jupyter:
3+
jupytext:
4+
text_representation:
5+
extension: .md
6+
format_name: markdown
7+
format_version: '1.3'
8+
jupytext_version: 1.16.7
9+
kernelspec:
10+
display_name: Python 3
11+
language: python
12+
name: python3
13+
---
14+
15+
# Calculation of control fields for Hadamard gate on single qubit using L-BFGS-B algorithm
16+
17+
18+
Alexander Pitchford ([email protected])
19+
20+
21+
Example to demonstrate using the control library to determine control
22+
pulses using the ctrlpulseoptim.optimize_pulse_unitary function.
23+
The (default) L-BFGS-B algorithm is used to optimise the pulse to
24+
minimise the fidelity error, which is equivalent maximising the fidelity
25+
to an optimal value of 1.
26+
27+
The system in this example is a single qubit in a constant field in z
28+
with a variable control field in x
29+
The target evolution is the Hadamard gate irrespective of global phase
30+
31+
The user can experiment with the timeslicing, by means of changing the
32+
number of timeslots and/or total time for the evolution.
33+
Different initial (starting) pulse types can be tried.
34+
The initial and final pulses are displayed in a plot
35+
36+
An in depth discussion of using methods of this type can be found in [1]
37+
38+
```python
39+
import datetime
40+
41+
import matplotlib.pyplot as plt
42+
import numpy as np
43+
44+
import qutip_qtrl.pulseoptim as cpo
45+
from qutip import gates, identity, sigmax, sigmaz, about
46+
47+
example_name = "Hadamard"
48+
49+
%matplotlib inline
50+
```
51+
52+
### Defining the physics
53+
54+
55+
The dynamics of the system are governed by the combined Hamiltonian:
56+
H(t) = H_d + sum(u1(t)*Hc1 + u2(t)*Hc2 + ....)
57+
That is the time-dependent Hamiltonian has a constant part (called here the drift) and time vary parts, which are the control Hamiltonians scaled by some functions u_j(t) known as control amplitudes
58+
In this case the drift is simply a rotation about z and the (time-varying) control is a rotation about x
59+
In theory this system is fully controllable (irrespective of global phase) and so any unitary target could be chosen; we have chosen the Hadamard gate.
60+
61+
```python
62+
# Drift Hamiltonian
63+
H_d = sigmaz()
64+
# The (single) control Hamiltonian
65+
H_c = [sigmax()]
66+
# start point for the gate evolution
67+
U_0 = identity(2)
68+
# Target for the gate evolution Hadamard gate
69+
U_targ = gates.hadamard_transform(1)
70+
```
71+
72+
### Defining the time evolution parameters
73+
74+
75+
To solve the evolution the control amplitudes are considered constant within piecewise timeslots, hence the evolution during the timeslot can be calculated using U(t_k) = expm(-i*H(t_k)*dt). Combining these for all the timeslots gives the approximation to the evolution from the identity at t=0 to U(T) at the t=evo_time
76+
The number of timeslots and evo_time have to be chosen such that the timeslot durations (dt) are small compared with the dynamics of the system.
77+
78+
```python
79+
# Number of time slots
80+
n_ts = 10
81+
# Time allowed for the evolution
82+
evo_time = 10
83+
```
84+
85+
### Set the conditions which will cause the pulse optimisation to terminate
86+
87+
88+
At each iteration the fidelity of the evolution is tested by comparaing the calculated evolution U(T) with the target U_targ. For unitary systems such as this one this is typically:
89+
f = normalise(overlap(U(T), U_targ))
90+
For details of the normalisation see [1] or the source code.
91+
The maximum fidelity (for a unitary system) calculated this way would be 1, and hence the error is calculated as fid_err = 1 - fidelity. As such the optimisation is considered completed when the fid_err falls below such a target value.
92+
93+
In some cases the optimisation either gets stuck in some local minima, or the fid_err_targ is just not achievable, therefore some limits are set to the time/effort allowed to find a solution.
94+
95+
The algorithm uses gradients to direct its search for the minimum fidelity error. If the sum of all the gradients falls below the min_grad, then it is assumed some local minima has been found.
96+
97+
```python
98+
# Fidelity error target
99+
fid_err_targ = 1e-10
100+
# Maximum iterations for the optisation algorithm
101+
max_iter = 200
102+
# Maximum (elapsed) time allowed in seconds
103+
max_wall_time = 120
104+
# Minimum gradient (sum of gradients squared)
105+
# as this tends to 0 -> local minima has been found
106+
min_grad = 1e-20
107+
```
108+
109+
### Set the initial pulse type
110+
111+
112+
The control amplitudes must be set to some initial values. Typically these are just random values for each control in each timeslot. These do however result in erratic optimised pulses. For this example, a solution will be found for any initial pulse, and so it can be interesting to look at the other initial pulse alternatives.
113+
114+
```python
115+
# pulse type alternatives: RND|ZERO|LIN|SINE|SQUARE|SAW|TRIANGLE|
116+
p_type = "RND"
117+
```
118+
119+
### Give an extension for output files
120+
121+
```python
122+
# Set to None to suppress output files
123+
f_ext = "{}_n_ts{}_ptype{}.txt".format(example_name, n_ts, p_type)
124+
```
125+
126+
### Run the optimisation
127+
128+
129+
In this step the L-BFGS-B algorithm is invoked. At each iteration the gradient of the fidelity error w.r.t. each control amplitude in each timeslot is calculated using an exact gradient method (see [1]). Using the gradients the algorithm will determine a set of piecewise control amplitudes that reduce the fidelity error. With repeated iterations an approximation of the Hessian matrix (the 2nd order differentials) is calculated, which enables a quasi 2nd order Newton method for finding a minima. The algorithm continues until one of the termination conditions defined above has been reached.
130+
131+
```python
132+
result = cpo.optimize_pulse_unitary(
133+
H_d,
134+
H_c,
135+
U_0,
136+
U_targ,
137+
n_ts,
138+
evo_time,
139+
fid_err_targ=fid_err_targ,
140+
min_grad=min_grad,
141+
max_iter=max_iter,
142+
max_wall_time=max_wall_time,
143+
out_file_ext=f_ext,
144+
init_pulse_type=p_type,
145+
gen_stats=True,
146+
)
147+
```
148+
149+
### Report the results
150+
151+
152+
Firstly the performace statistics are reported, which gives a breadown of the processing times. The times given are those that are associated with calculating the fidelity and the gradients. Any remaining processing time can be assumed to be used by the optimisation algorithm (L-BFGS-B) itself. In this example it can be seen that the majority of time is spent calculating the propagators, i.e. exponentiating the combined Hamiltonian.
153+
154+
The optimised U(T) is reported as the 'final evolution', which is essentially the string representation of the Qobj that holds the full time evolution at the point when the optimisation is terminated.
155+
156+
The key information is in the summary (given) last. Here the final fidelity is reported and the reasonn for termination of the algorithm.
157+
158+
```python
159+
result.stats.report()
160+
print("Final evolution\n{}\n".format(result.evo_full_final))
161+
print("********* Summary *****************")
162+
print("Final fidelity error {}".format(result.fid_err))
163+
print("Final gradient normal {}".format(result.grad_norm_final))
164+
print("Terminated due to {}".format(result.termination_reason))
165+
print("Number of iterations {}".format(result.num_iter))
166+
print(
167+
"Completed in {} HH:MM:SS.US".format(datetime.timedelta(seconds=result.wall_time))
168+
)
169+
```
170+
171+
### Plot the initial and final amplitudes
172+
173+
174+
Here the (random) starting pulse is plotted along with the pulse (control amplitudes) that was found to produce the target gate evolution to within the specified error.
175+
176+
```python
177+
fig1 = plt.figure()
178+
ax1 = fig1.add_subplot(2, 1, 1)
179+
ax1.set_title("Initial control amps")
180+
# ax1.set_xlabel("Time")
181+
ax1.set_ylabel("Control amplitude")
182+
ax1.step(
183+
result.time,
184+
np.hstack((result.initial_amps[:, 0], result.initial_amps[-1, 0])),
185+
where="post",
186+
)
187+
188+
ax2 = fig1.add_subplot(2, 1, 2)
189+
ax2.set_title("Optimised Control Sequences")
190+
ax2.set_xlabel("Time")
191+
ax2.set_ylabel("Control amplitude")
192+
ax2.step(
193+
result.time,
194+
np.hstack((result.final_amps[:, 0], result.final_amps[-1, 0])),
195+
where="post",
196+
)
197+
plt.tight_layout()
198+
plt.show()
199+
```
200+
201+
### Versions
202+
203+
```python
204+
about()
205+
```
206+
207+
### References
208+
209+
210+
[1] Machnes et.al., DYNAMO - Dynamic Framework for Quantum Optimal Control. arXiv.1011.4874

0 commit comments

Comments
 (0)