# Flexible optimization

Highly-configurable non-linear optimization framework for quantum control

The Q-CTRL Python package exposes a highly-flexible optimization engine for general-purpose gradient-based optimization. This feature may be used for systems that do not fit into the structure supported by the basic Optimization feature.

The optimization engine provides a large, modular collection of configuration options in order to support a wide variety of systems and constraints. These options can be combined in different ways to produce optimization solutions tailored to specific systems. In this guide we demonstrate optimization using several different types of systems and constraints, to give a taste of the available flexibility. All the techniques presented here can be cross-combined to yield an even wider class of supported configurations.

## Imports and initialization

# Essential imports
import numpy as np
from qctrl import Qctrl

# Plotting imports
from qctrlvisualizer import plot_controls
import matplotlib.pyplot as plt

# Starting a session with the API
qctrl = Qctrl()


## Worked example: robust control of a single qubit

We first present a detailed worked example of robust optimization in a single-qubit system. Specifically, we consider a single-qubit system represented by the following Hamiltonian:

\begin{align*} H(t) &= \frac{\nu}{2} \sigma_{z} + \frac{1}{2}\left(\gamma(t)\sigma_{-} + \gamma^*(t)\sigma_{+}\right) + \frac{\alpha(t)}{2} \sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\nu$ is the qubit detuning, $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, and $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process and $\sigma_{\{+,-,z\}}$ are the Pauli matrices.

The functions of time $\gamma(t)$ and $\alpha(t)$ are not predetermined, and instead are optimized by the Q-CTRL optimization engine in order to achieve some target operation.

### Defining the graph

The flexible Q-CTRL optimization engine expresses all optimization problems as data flow graphs, which describe how optimization variables (variables that can be tuned by the optimizer) are transformed into the cost function (the objective that the optimizer attempts to minimize).

Below we show how to create a data flow graph for optimizing the single-qubit system described above. Comments in the code explain the details of each step.

# Define standard matrices
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 0], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])

# Define physical constants
nu = 2*np.pi * 0.5 * 1e6 #Hz
gamma_max = 2*np.pi * 0.5e6 #Hz
alpha_max = 2*np.pi * 0.5e6 #Hz
segment_count = 50
duration = 10e-6 #s

# Define the data flow graph describing the system
with qctrl.create_graph():

# Create a constant piecewise-constant (PWC) operator representing the
# detuning term
detuning = qctrl.operations.constant_pwc_operator(
duration=duration,
operator=nu * sigma_z / 2,
)

# Create a complex PWC signal, with optimizable modulus and phase,
# representing gamma(t)
gamma = qctrl.operations.complex_pwc_signal(
moduli=qctrl.operations.bounded_optimization_variable(
count=segment_count, lower_bound=0, upper_bound=gamma_max,
),
phases=qctrl.operations.unbounded_optimization_variable(
count=segment_count,
initial_lower_bound=0,
initial_upper_bound=2*np.pi,
),
duration=duration,
name="gamma",
)
# Create a PWC operator representing the drive term
drive = qctrl.operations.pwc_operator_hermitian_part(
qctrl.operations.pwc_operator(signal=gamma, operator=sigma_m)
)

# Create a real PWC signal, with optimizable amplitude, representing
# alpha(t)
alpha = qctrl.operations.pwc_signal(
values=qctrl.operations.bounded_optimization_variable(
count=segment_count, lower_bound=-alpha_max, upper_bound=alpha_max
),
duration=duration,
name="alpha",
)
# Create a PWC operator representing the clock shift term
shift = qctrl.operations.pwc_operator(alpha, sigma_z / 2)

# Create a constant PWC operator representing the dephasing noise
# (note that we scale by 1/duration to ensure consistent units between
# the noise Hamiltonian and the control Hamiltonian)
dephasing = qctrl.operations.constant_pwc_operator(duration, sigma_z / duration)

# Create the target operator
target_operator = qctrl.operations.target(operator=Y)

# Create infidelity
infidelity = qctrl.operations.infidelity_pwc(
hamiltonian=qctrl.operations.pwc_sum([detuning, drive, shift]),
noise_operators=[dephasing],
target_operator=target_operator,
)


### Creating the graph object

Once the graph has been defined inside the qctrl.create_graph context manager, an object must be created that ties together the objective function (in this case the infidelity) and the desired outputs (identified by names, in this case for the Rabi rate and clock shift).

graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha", "gamma"]
)


### Running an optimization

With the graph object created, an optimization can be run using the flexible_optimize service. The service returns a new graph object containing the results of the optimization.

graph = qctrl.services.flexible_optimize.run(graph)

100%|██████████| 100/100 [00:11<00:00,  1.67it/s]

print("Optimized cost:\t", graph.result['cost'])

Optimized cost:	 2.303114180068945e-13


### Visualizing the pulses

We may use the qctrlvisualizer package to plot the optimized pulses, which are available in the graph object.

plot_controls(plt.figure(), graph.result['output'])
plt.show() ## Example: time-symmetrized pulses

Next we show how time symmetry can be incorporated into the optimization problem. In some cases, forcing control pulses to exhibit time symmetry can yield natural noise robustness while simplifying the optimization problem (by halving the dimensionality of the search space). In this example, we perform a robust optimization of a single qubit using symmetric pulses. The single-qubit system is represented by the following Hamiltonian:

\begin{align*} H(t) &= \frac{\nu}{2} \sigma_{z} + \frac{1}{2}\left(\gamma(t)\sigma_{-} + \gamma^*(t)\sigma_{+}\right) + \frac{\alpha(t)}{2} \sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\nu$ is the qubit detuning, $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, and $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process.

To enforce time symmetry, we create the control signals in two steps: we first create a standard signal for the first half of the gate duration, and then we create a copy (depending on the same underlying control parameters) that gets reflected and concatenated with the original signal.

# Define standard matrices
sigma_z = np.array([[1, 0], [0, -1]])
sigma_m = np.array([[0, 1], [0, 0]])
Y = np.array([[0, -1j], [1j, 0]])

# Define physical constraints
gamma_max = 2 * np.pi * 8.5e6 #Hz
alpha_max = 2 * np.pi * 8.5e6 #Hz
nu = 2 * np.pi * 6e6 #Hz
segment_count = 50
duration = 154e-9 #s

with qctrl.create_graph():

# Create detuning term
detuning = qctrl.operations.constant_pwc_operator(duration, nu * sigma_z / 2)

# Create a complex PWC signal describing the first half of gamma(t)
half_gamma = qctrl.operations.complex_pwc_signal(
moduli=qctrl.operations.bounded_optimization_variable(
count=segment_count, lower_bound=0, upper_bound=gamma_max
),
phases=qctrl.operations.unbounded_optimization_variable(
count=segment_count,
initial_lower_bound=0,
initial_upper_bound=2 * np.pi,
),
duration=duration / 2,
)
# Define gamma(t) by symmetrizing half_gamma
gamma = qctrl.operations.symmetrize_pwc(
half_gamma,
name="gamma",
)
# Create drive term
drive = qctrl.operations.pwc_operator_hermitian_part(
qctrl.operations.pwc_operator(gamma, sigma_m)
)

# Create alpha(t) similarly
alpha = qctrl.operations.symmetrize_pwc(
qctrl.operations.pwc_signal(
values=qctrl.operations.bounded_optimization_variable(
count=segment_count, lower_bound=-alpha_max, upper_bound=alpha_max
),
duration=duration / 2,
),
name="alpha",
)
# Create clock shift term
shift = qctrl.operations.pwc_operator(alpha, sigma_z / 2)

# Create dephasing noise term
dephasing = qctrl.operations.constant_pwc_operator(duration, sigma_z / duration)

# Create infidelity
infidelity = qctrl.operations.infidelity_pwc(
qctrl.operations.pwc_sum([detuning, drive, shift]),
qctrl.operations.target(Y),
[dephasing],
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["gamma", "alpha"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [00:39<00:00,  2.52it/s]

Optimized cost:	 1.5084890137873168e-10 ## Example: non-linear dependence on control pulses

We next present an example showing how a system with non-linear dependence between controls and Hamiltonian may be optimized. We consider an artificial single-qubit system represented by the following Hamiltonian:

\begin{align*} H(t) &= \alpha_1(t)\sigma_x + \frac{\alpha_1(t)^2}{\alpha_{\text{max}}}\sigma_y + \alpha_2(t)\sigma_z \,, \end{align*}

where $\alpha_1(t)$ and $\alpha_2(t)$ are real time-dependent pulses, and we note the non-linear dependence of the Hamiltonian on $\alpha_1(t)$. While here we have chosen a simple artificial non-linearity for the sake of convenience and brevity, the approach we demonstrate extends trivially to far more complex, physically-relevant situations.

Expressing non-linear dependence between control parameters and Hamiltonians is achieved by performing standard arithmetical operations on the control parameters prior to constructing signals. In this case the operations are very simple—we simply square and re-scale the parameters describing the $\alpha_1(t)$ values—but in general any combination of supported primitive operations may be used.

# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_y = np.array([[0, -1j], [1j, 0]])
sigma_z = np.array([[1, 0], [0, -1]])

# Define physical constraints
alpha_max = 2*np.pi * 6e6 #Hz
segment_count = 40
duration = 200e-9 #s

with qctrl.create_graph():

# Create the values backing the alpha_1(t) signal
alpha_1_values = qctrl.operations.bounded_optimization_variable(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
)

# Create the alpha_1(t) signal
alpha_1 = qctrl.operations.pwc_signal(
values=alpha_1_values, duration=duration, name="alpha_1"
)

# Apply a non-linear transformation to create the alpha_1_squared(t) signal
alpha_1_squared_values = alpha_1_values * (alpha_1_values / alpha_max)
alpha_1_squared = qctrl.operations.pwc_signal(
values=alpha_1_squared_values, duration=duration, name="alpha_1_squared"
)

# Create the alpha_2(t) signal
alpha_2 = qctrl.operations.pwc_signal(
values=qctrl.operations.bounded_optimization_variable(
count=segment_count, lower_bound=-alpha_max, upper_bound=alpha_max
),
duration=duration,
name="alpha_2",
)

# Create Hamiltonian terms
x_term = qctrl.operations.pwc_operator(alpha_1, sigma_x)
y_term = qctrl.operations.pwc_operator(alpha_1_squared, sigma_y)
z_term = qctrl.operations.pwc_operator(alpha_2, sigma_z)

# Create infidelity
infidelity = qctrl.operations.infidelity_pwc(
qctrl.operations.pwc_sum(terms=[x_term, y_term, z_term]),
qctrl.operations.target(operator=sigma_x),
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha_1", "alpha_1_squared", "alpha_2"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [00:31<00:00,  2.21it/s]

Optimized cost:	 0.0 ## Example: band-limited pulses with bounded slew rates

Next we show how to optimize a system in which the rates of change of the controls are limited. Using this constraint can help to ensure that optimized controls can be reliably implemented on physical hardware (which may be subject to bandwidth limits). We consider a standard single-qubit system subject to dephasing noise:

\begin{align*} H(t) &= \frac{1}{2} \alpha_1(t)\sigma_{x} + \frac{1}{2} \alpha_2(t) \sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\alpha_1(t)$ and $\alpha_2(t)$ are real time-dependent pulse and $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process. In this case, we enforce a maximum slew rate constraint on $\alpha_1(t)$ and $\alpha_2(t)$, to cap the variation between adjacent segment values.

To implement a bounded slew rate, we create signal values using a Q-CTRL helper function, anchored_difference_bounded_variables. This function produces values that are constrained to satisfy the slew rate requirement, and in addition are anchored to zero at the start and end of the gate.

# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
Y = np.array([[0, -1j], [1j, 0]])

# Define physical constraints
alpha_max = 2 * np.pi * 8.5e6  # Hz
max_slew_rate = alpha_max / 10
segment_count = 250
duration = 400e-9  # s

with qctrl.create_graph():

# Create alpha_1(t) signal
alpha_1_values = qctrl.operations.anchored_difference_bounded_variables(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
difference_bound=max_slew_rate,
)
alpha_1 = qctrl.operations.pwc_signal(
values=alpha_1_values, duration=duration, name="alpha_1"
)

# Create alpha_2(t) signal
alpha_2_values = qctrl.operations.anchored_difference_bounded_variables(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
difference_bound=max_slew_rate,
)
alpha_2 = qctrl.operations.pwc_signal(
values=alpha_2_values, duration=duration, name="alpha_2"
)

# Create drive term
drive = qctrl.operations.pwc_operator(alpha_1, sigma_x / 2)
# Create clock shift term
shift = qctrl.operations.pwc_operator(alpha_2, sigma_z / 2)

# Create dephasing noise term
dephasing = qctrl.operations.constant_pwc_operator(duration, sigma_z / duration)

# Create infidelity
infidelity = qctrl.operations.infidelity_pwc(
qctrl.operations.pwc_sum([drive, shift]),
qctrl.operations.target(Y),
[dephasing],
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha_1", "alpha_2"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [05:11<00:00,  2.12s/it]

Optimized cost:	 1.587647027965894e-11 ## Example: band-limited pulses with filters

Next we show how to optimize a system using a different technique for enforcing pulse band limits, namely filtering. As above, producing band-limited controls can help to ensure that optimized controls can be reliably implemented on physical hardware, but in some cases one desires tighter control over the frequency domain behavior of the controls than what can be achieved with a slew rate bound. The technique we use here is to pass the controls through a linear time-invariant filter in order to enforce explicit bandwidth limits. We consider, once again, a standard single-qubit system subject to dephasing noise:

\begin{align*} H(t) &= \frac{1}{2} \alpha_1(t)\sigma_{x} + \frac{1}{2} \alpha_2(t) \sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\alpha_1(t)$ and $\alpha_2(t)$ are real time-dependent pulse and $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process. In this case we force $\alpha_1(t)$ and $\alpha_2(t)$ to be signals produced from an ideal sinc filter, which forces the signals to exhibit a hard upper bound in the frequency domain.

To implement this behavior, we create optimizable signals as usual, and then use the convolve_pwc (to apply a filter) and discretize_stf (to obtain a piecewise-constant discretization of the smooth filtered signal) functions to derive band-limited signals. These band-limited signals are then used to construct the Hamiltonian.

# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_z = np.array([[1, 0], [0, -1]])
Y = np.array([[0, -1j], [1j, 0]])

# Define physical constraints
alpha_max = 2 * np.pi * 8.5e6  #Hz
sinc_cutoff_frequency = 2 * np.pi * 48e6  #Hz
nu = 2 * np.pi * 6e6  #Hz
segment_count = 50
duration = 400e-9  #s

with qctrl.create_graph():

# Create initial alpha_1(t) signal
alpha_1_values = qctrl.operations.bounded_optimization_variable(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
)
alpha_1 = qctrl.operations.pwc_signal(
values=alpha_1_values, duration=duration,
)
# Create filtered signal
alpha_1_filtered = qctrl.operations.convolve_pwc(
alpha_1,
qctrl.operations.sinc_integral_function(sinc_cutoff_frequency),
)
# Discretize to obtain smoothed alpha_1(t) signal
alpha_1 = qctrl.operations.discretize_stf(
stf=alpha_1_filtered,
duration=duration,
segments_count=segment_count,
name="alpha_1",
)

# Similarly, create alpha_2(t) signal
alpha_2_values = qctrl.operations.bounded_optimization_variable(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
)
alpha_2 = qctrl.operations.pwc_signal(
values=alpha_2_values, duration=duration,
)
# Create filtered signal
alpha_2_filtered = qctrl.operations.convolve_pwc(
alpha_2,
qctrl.operations.sinc_integral_function(sinc_cutoff_frequency),
)
# Discretize to obtain smoothed alpha_2(t) signal
alpha_2 = qctrl.operations.discretize_stf(
stf=alpha_2_filtered,
duration=duration,
segments_count=segment_count,
name="alpha_2",
)

# Create drive term
drive = qctrl.operations.pwc_operator(alpha_1, sigma_x / 2)
# Create clock shift term
shift = qctrl.operations.pwc_operator(alpha_2, sigma_z / 2)

# Create dephasing noise term
dephasing = qctrl.operations.constant_pwc_operator(duration, sigma_z / duration)

# Create infidelity
infidelity = qctrl.operations.infidelity_pwc(
qctrl.operations.pwc_sum([drive, shift]),
qctrl.operations.target(Y),
[dephasing],
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha_1", "alpha_2"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [00:38<00:00,  1.66it/s]

Optimized cost:	 3.695218195466676e-12 ## Example: incorporating smoothing via linear filters on control signals

We next present an example showing how systems containing linear time-invariant filters can be optimized. Unlike the example above, where we employed such filters in order to produce smooth controls, here we consider the case where such filters are present in the system itself. In this particular example, we consider the common scenario in which electronic filters are present on the control lines. Specifically, we consider a basic single-qubit system:

\begin{align*} H(t) &= \frac{1}{2} L(\alpha_1)(t)\sigma_{x} + \frac{1}{2} L(\alpha_2)(t) \sigma_{z} + \beta(t) \sigma_{z} \,, \end{align*}

where $\alpha_1(t)$ and $\alpha_2(t)$ are real time-dependent pulses, $\beta(t)$ is a small, slowly-varying stochastic dephasing noise process, and $L$ is a sinc filter applied to the pulses. The effect of these filters is to transform the control signals before they reach the quantum system, via convolution with the filter impulse response. Failing to take the filters into account during the optimization can lead to poor results, since in that case the system used for the optimization does not accurately model reality.

As above, we use the convolve_pwc function to apply a filter (obtained via the sinc_integral_function function) to a piecewise-constant signal. The resulting smooth tensor function (STF) may then be manipulated in exactly the same way as piecewise-constant functions, up to and including calculating the infidelity. We stress that even though highly irregular pulses are generated for $\alpha_1(t)$ and $\alpha_2(t)$, these pulses are chosen so that, when smoothed by the filters in the system, they will effect a high-quality gate.

# Define standard matrices
sigma_x = np.array([[0, 1], [1, 0]])
sigma_z = np.array([[1, 0], [0, -1]])

# Define physical constraints
alpha_max = 2*np.pi * 8.5e6 #Hz
nu = 2*np.pi * 6e6 #Hz
sinc_cutoff_frequency = 2*np.pi * 48e6 #Hz
segment_count = 50
duration = 250e-9 #s

with qctrl.create_graph():
# Create alpha_1(t) signal
alpha_1_values = qctrl.operations.bounded_optimization_variable(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
)
alpha_1 = qctrl.operations.pwc_signal(
values=alpha_1_values, duration=duration, name="alpha_1",
)
# Create filtered signal
alpha_1_filtered = qctrl.operations.convolve_pwc(
alpha_1,
qctrl.operations.sinc_integral_function(sinc_cutoff_frequency),
)

# Similarly, create filtered alpha_2(t) signal
alpha_2_values = qctrl.operations.bounded_optimization_variable(
count=segment_count,
lower_bound=-alpha_max,
upper_bound=alpha_max,
)
alpha_2 = qctrl.operations.pwc_signal(
values=alpha_2_values, duration=duration, name="alpha_2",
)
alpha_2_filtered = qctrl.operations.convolve_pwc(
alpha_2,
qctrl.operations.sinc_integral_function(sinc_cutoff_frequency),
)

# Create drive term (note the use of STF functions instead of PWC functions,
# because we are dealing with smooth signals instead of piecewise-constant
# signals).
drive = qctrl.operations.stf_operator(alpha_1_filtered, sigma_x / 2)
# Create clock shift term
shift = qctrl.operations.stf_operator(alpha_2_filtered, sigma_z / 2)

# Create dephasing noise term
dephasing = qctrl.operations.constant_stf_operator(sigma_z / duration)

# Create infidelity (note that we pass an array of sample times, which
# governs the granularity of the integration procedure)
infidelity = qctrl.operations.infidelity_stf(
np.linspace(0, duration, 150),
qctrl.operations.stf_sum([drive, shift]),
qctrl.operations.target(sigma_x),
[dephasing],
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha_1", "alpha_2"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [00:29<00:00,  1.39s/it]

Optimized cost:	 8.070826425341707e-11 ## Example: CRAB optimization

In chopped random basis (CRAB) optimization, pulses are defined via optimizable linear combinations from a set of basis functions, which can greatly reduce the dimensionality of the optimization search space. Traditionally, a randomized Fourier basis is used, although the same technique has also seen success with other bases, for example Slepian functions. In this example, we perform a CRAB optimization (in the Fourier basis) of a qutrit system in which we effect a single-qubit gate while minimizing leakage out of the computational subspace. The system is described by the following Hamiltonian:

\begin{align*} H(t) = & \frac{\chi}{2} (a^\dagger)^2 a^2 + \gamma(t) a + \gamma^*(t) a^\dagger + \frac{\alpha(t)}{2} a^\dagger a \,, \end{align*}

where $\chi$ is the anharmonicity, $\gamma(t)$ and $\alpha(t)$ are, respectively, complex and real time-dependent pulses, and $a = |0 \rangle \langle 1 | + \sqrt{2} |1 \rangle \langle 2 |$.

The Q-CTRL optimization engine provides a convenience function, real_fourier_pwc_signal, for creating optimizable signals in a Fourier basis, suitable for use in a CRAB optimization. Other bases are supported by the framework, but require the user to manually provide operations that compute the appropriate linear combinations.

# Define standard matrices
a = np.array([[0., 1., 0.],
[0., 0., np.sqrt(2)],
[0., 0., 0.]], dtype=np.complex)
[1., -1.,0],
[0, 0, np.sqrt(2)]], dtype=np.complex)/np.sqrt(2)
qubit_projector = np.array([[1, 0, 0],
[0, 1, 0],
[0, 0, 0]], dtype=np.complex)

# Define physical constraints
chi = 2 * np.pi * -300.0 * 1e6 #Hz
gamma_max = 2 * np.pi * 30e6 #Hz
alpha_max = 2 * np.pi * 30e6 #Hz
segment_count = 200
duration = 100e-9 #s

with qctrl.create_graph():
# Create gamma(t) signal in Fourier bases. To demonstrate the full
# flexibility, we show how to use both randomized and optimizable
# basis elements. Elements with fixed frequencies may be chosen too.
gamma_i = qctrl.operations.real_fourier_pwc_signal(
duration=duration,
segments_count=segment_count,
randomized_frequencies_count=10,
)
gamma_q = qctrl.operations.real_fourier_pwc_signal(
duration=duration,
segments_count=segment_count,
optimizable_frequencies_count=10,
)
gamma = qctrl.operations.tensor_pwc(
gamma_i.durations,
qctrl.operations.complex_value(real=gamma_i.values, imag=gamma_q.values) * gamma_max,
name="gamma",
)

# Create alpha(t) signal
alpha = qctrl.operations.real_fourier_pwc_signal(
duration=duration,
segments_count=segment_count,
initial_coefficient_lower_bound=-alpha_max,
initial_coefficient_upper_bound=alpha_max,
optimizable_frequencies_count=10,
name="alpha",
)

# Create anharmonicity term
anharmonicity = qctrl.operations.constant_pwc_operator(
duration, ad2a2 * chi / 2,
)
# Create drive term
drive = qctrl.operations.pwc_operator_hermitian_part(
qctrl.operations.pwc_operator(gamma, 2 * a)
)
# Create clock shift term
shift = qctrl.operations.pwc_operator(alpha, ada / 2,)

# Create infidelity (note the use of a projector, to calculate the infidelity
# in the qubit subspace)
infidelity = qctrl.operations.infidelity_pwc(
qctrl.operations.pwc_sum([anharmonicity, drive, shift]),
)

# Create graph object
graph = qctrl.factories.graphs.create(
definition=infidelity, fetches=["alpha", "gamma"],
)

# Run the optimization
graph = qctrl.services.flexible_optimize.run(graph)

print("Optimized cost:\t", graph.result['cost'])

# Plot the optimized controls
plot_controls(plt.figure(), graph.result['output'])
plt.show()

100%|██████████| 100/100 [00:31<00:00,  1.83it/s]

Optimized cost:	 2.2890578321721478e-12 ## Wiki

Comprehensive knowledge base of quantum control theory