Optimal Control Toy Problems

Optimizing a control pulse requires a specific target and a restricted set of allowed pulses. Here the target is to transfer $|0\rangle$ to $|1\rangle$ in a fixed time while dephasing acts. A small deterministic search improves a square pulse within an amplitude bound, then checks whether the same pulse also performs well as a gate.

Helpful background. Solving Lindblad equations introduces numerical propagation. Optimal control discusses more general objectives and constraints. The self-contained example uses Python 3 and NumPy.

A bounded square pulse in a dephasing qubit

In a rotating-frame model, take

$$\frac H\hbar=\frac12(\Omega X+\Delta Z), \qquad \dot\rho=-i[H/\hbar,\rho]+\frac{\Gamma_\phi}{2}(Z\rho Z-\rho).$$

Here $\Delta$, $\Omega$, and $\Gamma_\phi$ have units of inverse time; $\Omega$ and $\Delta$ are angular frequencies. The pulse lasts $T=1$ in the chosen unit and has constant real amplitude $0\le\Omega\le4$. There is no amplitude relaxation or leakage in this model.

Writing $\rho=(I+\mathbf r\cdot\boldsymbol\sigma)/2$ gives

$$\dot{\mathbf r}=A\mathbf r, \qquad A=\begin{pmatrix} -\Gamma_\phi&-\Delta&0\\ \Delta&-\Gamma_\phi&-\Omega\\ 0&\Omega&0 \end{pmatrix}.$$

The noise acts during the rotation. In general, evolving with the drive and then applying a separate dephasing channel is not equivalent. With $\Delta=\Gamma_\phi=0$, a square pulse of area $\Omega T=\pi$ performs the desired transfer.

State transfer is different from gate fidelity

If $\mathbf r(T)=M\mathbf r(0)$, the population-transfer objective is

$$F_{\mathrm{st}}=\langle1|\rho(T)|1\rangle=\frac{1-M_{zz}}2.$$

It tests only the input at the north pole. For the target $R_x(\pi)$, the Bloch rotation is $R=\operatorname{diag}(1,-1,-1)$. The average gate fidelity is

$$F_{\mathrm{avg}}=\frac12+\frac{\operatorname{Tr}(R^{\mathsf T}M)}6.$$

For any trace-preserving qubit channel compared with a unitary target, averaging the pure-state overlaps over the six Pauli-axis states gives exactly this average gate fidelity. It is not merely a heuristic score: those six vectors have zero mean and second moment $I_3/3$, just like the uniform sphere. A possible affine translation cancels in the average. Process or entanglement fidelity uses a different normalization; for a qubit, $F_{\mathrm{avg}}=(2F_e+1)/3$.

Search, refine, and test the same pulse

The code propagates the constant Bloch generator by repeated fourth-order Runge–Kutta steps. For a linear constant equation, one step is the polynomial $I+hA+(hA)^2/2+(hA)^3/6+(hA)^4/24$. Raising this step matrix to the number of steps is equivalent to applying those steps successively.

A coarse search locates a good bracket; golden-section refinement searches inside it. This establishes the best result found within this one-parameter family and interval, not global optimality over all possible control functions.

import numpy as np

T = 1.0
R = np.diag([1.0, -1.0, -1.0])

def propagate(omega, gamma=0.3, detuning=0.0, steps=400):
    A = np.array([[-gamma, -detuning, 0],
                  [detuning, -gamma, -omega],
                  [0, omega, 0]], dtype=float)
    B = T*A/steps
    B2 = B@B
    step = np.eye(3)+B+B2/2+B2@B/6+B2@B2/24
    return np.linalg.matrix_power(step, steps)

def scores(M):
    transfer = (1-M[2, 2])/2
    average = 0.5+np.trace(R.T@M)/6
    six = np.vstack([np.eye(3), -np.eye(3)])
    six_average = np.mean([(1+(R@r)@(M@r))/2 for r in six])
    assert np.isclose(average, six_average)
    return transfer, average

assert np.allclose(propagate(np.pi, gamma=0), R, atol=1e-9)
assert np.allclose(propagate(0), np.diag([np.exp(-0.3), np.exp(-0.3), 1]))
grid = np.linspace(0, 4, 81)
values = [scores(propagate(w))[0] for w in grid]
k = int(np.argmax(values))
assert 0 < k < len(grid)-1
lo, hi = grid[k-1], grid[k+1]
ratio = (np.sqrt(5)-1)/2
for _ in range(50):
    a = hi-ratio*(hi-lo)
    b = lo+ratio*(hi-lo)
    if scores(propagate(a))[0] < scores(propagate(b))[0]:
        lo = a
    else:
        hi = b
omega_best = (lo+hi)/2
baseline = scores(propagate(np.pi))
best = scores(propagate(omega_best))
assert 0 <= omega_best <= 4
assert best[0] >= baseline[0]
assert np.linalg.norm(propagate(omega_best, steps=800)
                      -propagate(omega_best, steps=400)) < 1e-8
print('amplitude:', omega_best, 'baseline scores:', baseline, 'best scores:', best)

# Evaluate this fixed pulse at other parameters; do not reoptimize it.
for gamma, delta in [(0, 0), (0.3, 0.2), (0.6, 0.5)]:
    M = propagate(omega_best, gamma, delta, steps=800)
    assert np.linalg.svd(M, compute_uv=False).max() <= 1+1e-10
    print('gamma, detuning:', gamma, delta, 'scores:', scores(M))

The noiseless limit checks the factor $1/2$ in the Hamiltonian, and the no-drive limit checks the dephasing rate. Step refinement checks numerical propagation error. Comparing the fixed pulse at additional noise strengths and detunings tests parameter sensitivity; optimizing a new pulse at every point would answer a different question.

The search gives $\Omega\approx3.19145$, increasing transfer fidelity from about $0.930278$ to $0.930830$. At the same time, average gate fidelity decreases from about $0.910371$ to $0.910064$. The improvement over a calibrated $\pi$ pulse can be small. Compare its size with the propagation error before interpreting it physically. State transfer and average gate fidelity need not improve together, because they reward different behavior. The six-state identity assumes a linear trace-preserving channel; applying it to a normalized postselected map would require additional care.

Useful extensions with clear constraints

Try two or four constant pulse segments, each satisfying $\Omega_x^2+\Omega_y^2\le\Omega_{\max}^2$, and propagate them in time order. Integrate each segment separately so a numerical step does not straddle a discontinuity. A smooth raised-cosine pulse $\Omega_x(t)=\Omega_{\max}\sin^2(\pi t/T)$ has area $\Omega_{\max}T/2$; its noiseless $\pi$ pulse needs peak amplitude $2\pi/T$, which violates the bound $4$ used here when $T=1$. Comparisons must respect the same constraints.

If adding power or smoothness penalties, express them in dimensionless time $s=t/T$ and normalized controls $u=\Omega/\Omega_{\max}$, for example $\lambda\int_0^1u(s)^2ds$. Report physical fidelity separately from the penalized objective. A derivative penalty on a discontinuous piecewise-constant pulse is not a finite ordinary integral; use a stated discrete difference penalty or smooth the pulse family.

Markovian dephasing cannot generally be refocused away as if it were a constant unknown detuning. Pulse optimization can change how long a state occupies vulnerable superpositions, but the answer remains tied to the assumed noise model and amplitude limit.

Exercises

A perfect transfer that is a poor gate. Consider the channel that discards its input and prepares $|1\rangle$. What are its transfer and average fidelities relative to $R_x(\pi)$?

Solution

For input $|0\rangle$, $F_{\mathrm{st}}=1$. The channel has $M=0$ and a translation toward the south pole, so $F_{\mathrm{avg}}=1/2$. Perfect performance on one input does not certify a gate.

Why six inputs suffice for this average. Derive the average gate-fidelity formula for $\mathbf r'=M\mathbf r+\mathbf c$.

Solution

For a pure input, the overlap with the ideal output is $[1+(R\mathbf r)\cdot(M\mathbf r+\mathbf c)]/2$. Both the uniform sphere and the six Pauli vectors have $\langle\mathbf r\rangle=0$ and $\langle r_ir_j\rangle=\delta_{ij}/3$. Averaging therefore gives $1/2+\operatorname{Tr}(R^{\mathsf T}M)/6$.

Pulse area under a bound. Find the minimum duration of a noiseless resonant raised-cosine $\pi$ pulse with peak amplitude at most $4$.

Solution

Its area is $\Omega_{\max}T/2$. Setting this to $\pi$ gives $T\ge\pi/2$, whereas the square pulse can reach area $\pi$ in $T=\pi/4$ under the same peak bound.

A reliable improvement. What additional check is needed if the optimized fidelity exceeds the baseline by only $10^{-7}$?

Solution

Propagate both pulses with tighter numerical accuracy and verify that the change in their fidelity difference is much smaller than $10^{-7}$. Also refine the optimization tolerance. More displayed decimal places do not establish that the improvement is resolved.

References

• C. Brif, R. Chakrabarti, and H. Rabitz, “Control of quantum phenomena: Past, present and future,” New Journal of Physics 12, 075008 (2010). doi:10.1088/1367-2630/12/7/075008.
• D. D’Alessandro, Introduction to Quantum Control and Dynamics, Chapman & Hall/CRC (2007).
• S. J. Glaser et al., “Training Schrödinger’s cat: Quantum optimal control,” European Physical Journal D 69, 279 (2015). doi:10.1140/epjd/e2015-60464-1.
• M. A. Nielsen, “A simple formula for the average gate fidelity of a quantum dynamical operation,” Physics Letters A 303, 249–252 (2002). doi:10.1016/S0375-9601(02)01272-0; author manuscript.