Numerics.Optimal

Optimal control theory

This Wiki page deals with the numerical solution of the optimal control problem for a bilinear control system, formulated in terms of system matrices A, B, N and C, D. Here we use a Dirac-like bra-ket notation with the state vector |x⟩ denoting either the deviation from the ground state, |x⟩= |ψ⟩−|ψ₀⟩ for closed quantum systems, TDSE or the deviation from the equilibrium state |x⟩=ρ−ρ_e for open quantum systems, LvNE within the Lindblad approximation, using suitable vactorization of the (reduced) density matrix ρ.

Control functionals

The task in optimal control is to find field(s) u_k(t) that drive the system from an initial state |x(t=0)⟩=|x₀⟩ to a final state |x(t=T)⟩ such that an observable (a "target") is maximized at time T. Here we assume a single observable only. We also want the energy associated with the field to be bounded for practicality, and the time evolution to be physical. Formally, these are three requirements that can be expressed in terms of three different functionals:

• An output/target functional that should be maximal at (given/fixed) final time T
  
  which allows for a linear and/or quadratic output equation and where constant terms ⟨c_j|x_e⟩ and/or ⟨x_e|D_j|x_e⟩ coming from the shift |x⟩ → |x⟩−|x_e⟩ have been omitted. (Remember that the equilibrium is defined as A|x_e⟩=0, see below) The index j labels the single observable to maximized; multiple observables (with individual weights) could be maximized as well, but this would complicate the analysis at least. Often, the observables are projectors on target states, but other choices are also possible, as long as they can be written in the linear and/or quadratic form given above.

• A cost/energy functional that should be minimal
  
  where the penalty factors α_k > 0 balance the importance of the cost functional w. r. t. the target functional and/or balance the importance among the various field components. The shape functions s_k(t) can be used to enforce certain time-dependent shape(s) of the field pulse(s), e. g., to model typical switch on/off behavior of the control fields [3]. When controlling the dynamics of quantum systems by laser fields, the integrals in the above equation correspond to laser fluences.

• Deviation from exact evolution that should be minimal
  
  where the Lagrange multiplier |z(t)⟩ has been introduced to make sure that the state vector |x(t)⟩ exactly satisfies its evolution equation ∂_t|x(t)⟩ = L(t) |x(t)+x_e⟩. Since we require the evolution to be fulfilled by the Hermitian conjugate ⟨x(t)| as well, we consider the real part of the functional.

The inclusion of further constraints is also possible, see for example [7] and [9].

Control equations

The combined functional J≡J₁−J₂−J₃ becomes extremal if the following three conditions are satisfied simultaneously. (In principle, the derivation is straight-forward by standard variational calculus. For a detailed derivation we recommend the Ph. D. tutorial by Werschnik and Gross [9].)

1. State vector: is propagated forward from |x(t=0)⟩ = |x₀⟩−|x_e⟩ via
   
   which is the standard equation of motion, again with |b_k⟩ ≡ N_k |x_e⟩; for a derivation and its connection to quantum dynamics see the page on bilinear control systems.

2. Lagrange multiplier: is propagated backward from |z(t=T)⟩ = |c_j⟩+D_j|x(T)⟩+D_j|x_e⟩−|x_e⟩ via
   
   For the special case of anti-Hermitian evolution L (i.e. anti-Hermitian A and real symmetric N, for example with closed quantum systems, TDSE) we have the same propagators for state vector |x(t)⟩ and Lagrangian multiplier |z(t)⟩. This is why the Lagrange multiplier is often also interpreted as a quantum state vector in the TDSE case.

3. Optimal control field: is calculated as
   
   which obviously requires the knowledge of the state vector and the Lagrange multiplier.

Note that these equations are somewhat impractical. On the one hand, we need the optimal control field when propagating the state vector |x(t)⟩ and the Lagrange multiplier |z(t)⟩. On the other hand, we need both |x(t)⟩ and |z(t)⟩ to calculate the control field in the first place.

To get out of this dilemma, we can calculate the fields u_k(t) from |x⟩ and |z⟩ at the previous time step, t-Δt. The error of such an approach can be made small by reducing the time step, Δt, accordingly. However, in the second order scheme suggested in [1] and [2], the fields are also calculated at intermediate (half) time steps which is accomplished by a linearization to approximate the derivative of the field (given here only for one-dimensional control and b=0 to keep things simple)

which should allow for larger time steps than the first order scheme.

Numerical schemes

• The system of the three coupled control equations (1-3) can be solved iteratively in the following way:
  We first guess an initial field u⁽⁰⁾, for example constant in time and strong enough to induce notable dynamics, and use it to propagate the initial state |x(t=0)⟩=|x₀⟩−|x_e⟩ forward in time, see Eq. (1), yielding |x⁽⁰⁾(T)⟩.
  Then we iteratively (for step n):

  • Propagate backward the Lagrange multiplier |z⁽ⁿ⁾(t)⟩starting from |z⁽ⁿ⁾(t=T)⟩=|c_j⟩+D_j|x^(n-1)(T)⟩+D_j|x_e⟩−|x_e⟩ using Eq. (2) with the field
    
  • Propagate forward the state vector |x⁽ⁿ⁾(t)⟩ starting from |x(t=0)⟩=|x₀⟩−|x_e⟩ using Eq. (1) with the field
    

  Such algorithms are known to converge monotonically (and rapidly!), for proofs see the impressive work by Rabitz (and Ohtsuki) for the TDSE case [1,2] and for the LvNE case [4]. The condition for the total functional J to increase monotonically is that the observable (matrix D_j) has to be positively definite.

• These algorithms were generalized by the introduction of two coefficients η and ζ describing the mixing of fields obtained for the recent and the previous iteration steps [5,6]
  
  where again v(t) and u(t) are the fields used in backward and forward propagations, respectively. Monotonic convergence is found iff the two parameters η and ζ are between 0 and 2. For the special case of η=ζ=1, we retrieve the Zhu/Rabitz scheme given above, while for η=0 and ζ=1 we retrieve the Krotov method. For an example LvNE system, the convergence behavior for different η and ζ is compared numerically in [7], showing large variations. (But in how far can these results be generalized?)

• Instead of optimizing just the real part of a linear target operator (associated with vector c above), one may consider to optimize the square of that operator, as proposed in [1] for the TDSE case and in [4] for the LvNE case and which is common practice when optimizing populations of quantum states. Assuming here D_j=0, the resulting control equations are the same, except for an additional factor in the optimal control field
  
  Note that ⟨x(t)+x_e|z(t)+x_e⟩ is independent of t, at least in theory (but not in numerical practice!); hence it can be evaluated for arbitrary times. While some authors, e.g., [9] evaluate this at every time t, others recommend choosing t=0 for forward and t=T for backward propagations, see e.g. the Appendix of [1].

References

1. W. Zhu, J. Botina, H. Rabitz, J. Chem. Phys. 108, 1953 (1998)
2. W. Zhu, H. Rabitz, J. Chem. Phys. 109, 385 (1998)
3. K. Sundermann, R. de Vivie-Riedle, J. Chem. Phys. 110, 1896 (1999)
4. Y. Ohtsuki, W. Zhu, H. Rabitz, J. Chem. Phys. 110, 9825 (1999)
5. Y. Maday, G. Turinici, J. Chem. Phys. 118, 8191 (2003)
6. C. LeBris, Y. Maday, G. Turinici, in: Quantum Control: Mathematical and Numerical Challenges (2003)
7. Y. Ohtsuki, G. Turinici, H. Rabitz, J. Chem. Phys. 120, 5509 (2004)
8. J. Werschnik, E. K. U. Gross, J. Opt. B 7, S300 (2005)
9. J. Werschnik, E. K. U. Gross, J. Phys. B 40, R175 (2007)