Application of Optimal Control Theory to Fourier Transform Ion Cyclotron Resonance.

We study the application of Optimal Control Theory to Ion Cyclotron Resonance. We test the validity and the efficiency of this approach for the robust excitation of an ensemble of ions with a wide range of cyclotron frequencies. Optimal analytical solutions are derived in the case without any pulse constraint. A gradient-based numerical optimization algorithm is proposed to take into account limitation in the control intensity. The efficiency of optimal pulses is investigated as a function of control time, maximum amplitude and range of excited frequencies. A comparison with adiabatic and SWIFT pulses is done. On the basis of recent results in Nuclear Magnetic Resonance, this study highlights the potential usefulness of optimal control in Ion Cyclotron Resonance.

1. Introduction

Performing efficient and robust state control by means of external time-dependent system parameter is a fundamental challenge in many technological developments at macroscopic or microscopic scale [1,2,3,4]. In this latter case, open-loop control protocol, i.e., without any real time feedback from the experiment during the control process, is generally used for practical and technical reasons. The controls are only designed from a modeling of the system dynamics and the efficiency of the control scenario may suffer from the accuracy of the theoretical description. The robustness of a control process with respect to experimental imperfections is therefore a key parameter in view of experimental implementation. Different techniques extending from adiabatic pulses to optimal control theory (OCT) have been developed in this open-loop framework to find the pulse parameters [1,3,5,6,7,8]. Optimal control tackles the question of bringing a dynamical system from one state to another with minimum expenditure of time and resources [1,2,3,4]. The modern version of OCT was born in the 1960s with the Pontryagin Maximum Principle (PMP), which provides a general and rigorous mathematical framework for optimal control techniques [9,10,11,12,13]. OCT has become nowadays a key tool in many different domains extending from space dynamics to robotics or quantum mechanics [1,3,11]. Optimal process is defined from a cost functional (to minimize) which can depend on the state of the system and the control field. For systems with complex dynamics and optimization targets which are difficult to reach, it is necessary to use optimal control algorithms converging iteratively towards the optimal solution. The flexibility of this approach makes it possible to adapt this tool to any experimental situation. Generally, it is possible to include constraints in the algorithms to account for requirements related to a specific material or device [1,3]. The only relative limitation concerns the accuracy of the modeling, even if robustness can be improved by controlling simultaneously an ensemble of systems which differ by the values of one or several constant parameters [14,15,16,17]. OCT has been applied to quantum systems first in the context of physical chemistry to steer chemical reactions or control specific degrees of freedom [5,18], followed by control of spin dynamics [19,20] for applications in Nuclear Magnetic Resonance (NMR) [15,17,21,22,23,24,25] and Magnetic Resonance Imaging [26,27,28,29]. It has become a key tool in this domain to improve the efficiency and the sensitivity of standard experimental setups [3]. In NMR and quantum physics, a well-known optimization method is GRAPE [30], which is a gradient-based algorithm [1]. This approach has been applied with success in many different contexts.

Fourier-Transform Ion Cyclotron Resonance (ICR) Mass Spectrometry [31,32] is a type of mass spectrometry based on cyclotron frequency of ions in a fixed magnetic field [33,34,35]. Ions are excited at their resonant cyclotron frequencies to larger cyclotron radii by an electric field orthogonal to the magnetic field. After the excitation pulse, the ions rotate freely with a frequency characteristic of their mass. The image current induced by the ions on a pair of electrodes is detected. The Fourier transform of the resulting transient signal leads to the mass spectrum after a proper calibration. In a homogeneous magnetic field, ICR allows accessing the highest resolution available in mass spectrometry, while leading to extreme sensitivities. This spectrometry has experienced a recent renewal based on several methodological improvements and the search for very high resolutions, which are required to study complex biological or environmental mixtures. Several techniques developed and used in ICR have been inspired by equivalent approaches in NMR. An example is given by two-dimensional ICR [36,37,38,39,40,41] which was proposed in analogy to two-dimensional NMR spectroscopy [19,20]. Following this fruitful approach and given the success and efficiency of optimal control techniques in NMR, a question which naturally arises is the application of this method in ICR. This paper aims at taking a step toward the answer to this open issue. ICR Mass Spectrometry can provide very high resolution mass spectra over a large range of mass to charge ratio. In the ICR experiment, ions are initially at rest in the center of the trap, and they have to be excited to generate a resonant signal which can cover, in broad band experiments, frequencies from a few 10 s kHz for high up to 1 MHz or higher for the fastest species. However, this implies that all ions have to be excited over this frequency range in an even and controlled manner.

We explore in this study how optimal control can be used to design efficient and robust excitation pulses in ICR. To the best of our knowledge, this has never been studied. Due to the wide bandwidth of ICR signal, excitation pulses are usually simple chirped adiabatic pulses with a frequency sweep. Some variations have been proposed such as off-resonance monochromatic pulses for selective excitation of given ions. Based on the linearity of ion dynamics, it has also been proposed to generate pulses by Fourier synthesis from a given excitation profile, in an approach called SWIFT (for Stored-waveform Inverse Fourier Transform) [42,43,44]. Optimal control is expected to allow a much wider range of possibilities such as the control of trajectory for given initial and final positions of the ion packet and for a given range of frequencies. To evaluate the contribution of OCT in ICR, we consider in this study the simplest modeling which is experimentally relevant. The experiment is considered in a simplified environment, with a constant magnetic field and a time-dependent homogeneous electric field oriented along a single axis orthogonal to the magnetic field and with no static component. This geometry is unrealistic, as there is no trapping potential, but allows to consider the dynamics of the ions restricted to a plane with a pure cyclotron trajectory and a zero magnetron component. The time-dependent electric field aims at exciting in a robust manner an ensemble of different ions from the center of the cell to a final position which depends in a controlled way of the ion frequency. The linearity of ion dynamics simplifies drastically the derivation of the optimal control law [10,45,46]. If there is no constraint on the intensity of the electric field, linear quadratic optimal control theory (LQOCT) can be applied. Many mathematical results have been established in this case [10,47] and the optimal solution can be derived analytically. When constraints are accounted for, a numerical algorithm has to be used to solve the optimal equations. Note that very few studies have solved optimal control problems of linear systems at the microscopic scale [47,48,49,50,51]. ICR is an interesting example, relevant experimentally, to stimulate further work in this direction.

The remainder of this paper is organized as follows. The formulation of the control problem and the description of the model system are outlined in Section 2. After a brief introduction to the principles of OCT, we apply OCT to ICR in Section 3. We describe the optimal control algorithm which allows taking into account experimental constraints on the control field. Numerical results in different experimental situations are given in Section 4. A comparison is made with the adiabatic and the SWIFT approaches. We conclude in Section 5 with an outlook and future perspectives. The Rotating Wave Approximation (RWA) is discussed in Appendix A. Technical details about the adiabatic and SWIFT techniques are presented, respectively, in Appendix B and Appendix C. The application of LQOCT is described in Appendix D.

2. Formulation of the Control Problem

2.1. The Model System

We consider the simplest modeling of ion trajectories in ICR. The different ions in the experimental cell are confined in the -plane and are subjected to a constant magnetic field and a time-dependent electric field , respectively, along the z- and x-axes of the laboratory frame. Note that optimal control techniques can also be used if two control fields along the x- and y-directions are available. The dynamics are governed by the Lorentz’s equation: where , and are the mass, charge and speed of the ion k. denotes the time derivative of . Equation (1) can be expressed as: with the cyclotron frequency and . The coordinates and describe, respectively, the position and the speed of the ion k in the -plane. We assume that the frequency belongs to the interval , which is defined by the ion packet under study. As described below, the aim of the control process is to excite the different ions in a robust way with respect to the parameter .

The control problem can be defined as follows. Starting from the center of the cell with a zero speed , the goal is to reach at a fixed control time a given radius and phase . As an illustrative example, we force the phase to vary linearly with , contrary to the standard result obtained with chirp pulses, where a quadratic phase dependence is observed (see Appendix B for details). We denote by and , respectively, the radius and the phase of ion k at time t. We assume in a first step that there is no constraint on the electric field. A limitation on the maximum pulse intensity is accounted for in Section 3.2.

To simplify the notations, we omit below the index k. Using Equation (2), it is straightforward to show that is a constant of motion. At , since and , we deduce that so . One of the two coordinates or can be eliminated. This also means that we cannot control simultaneously the position and the speed of the ion with only one control. We arrive at: where . We introduce the vector whose dynamics are governed by: with

The dynamics of this linear system can be explicitly integrated as follows. The eigenvalues of A are and the corresponding eigenvectors can be written as:

At time , the state of the system is given by:

We have: where and

We deduce that: and

2.2. The Rotating Wave Approximation

The oscillating excitation field applied only along the x- axis can be expressed as the sum of two rotating fields, one in the same direction as the ions and the other in the opposite direction. We introduce the Rotating Wave Approximation (RWA) which assumes that the field rotating in opposite direction to the ions has a negligible effect on their trajectories. This approximation is verified if the range of frequencies around the central frequency is not too large, as discussed in Appendix A. Note that RWA is a standard tool in NMR [20,52,53] where it is derived in a similar but different way due to the non-linearity of the system [54]. In particular for ICR, this approximation does not depend on the amplitude of the excitation. Using RWA, we show below that the control of ions is equivalent to the control of an ensemble of springs of different frequencies [47,51].

The derivation starts with the control of speeds which fulfill:

In complex coordinates, we have: where . We consider that where is the carrier frequency of the electric field, , and is small compared to . We also assume that the amplitude and the phase vary slowly in time with respect to the frequency . We express the speed as: , where is the complex speed in the frame rotating at frequency . We deduce that: where is the detuning term. In the RWA, we neglect the rapidly oscillating term and we arrive at:

It is worth noting here that, in the rotating frame, the dynamics are driven by two control parameters, and . Note that we recover the control of an ensemble of springs. An additional step can be done for the position of the ion k, . We set . It is then straightforward to show that:

Since varies slowly with respect to , we can neglect the time derivative , which gives:

If the RWA is valid, we deduce that the speed control leads also to the control of the position of ions. In this study, the validity of RWA is verified in the different examples by a numerical integration of Equation (4).

3. Optimal Control Theory

3.1. A Short Introduction to Optimal Control Theory

We briefly introduce in this section the tools of optimal control theory used in this paper. To keep the introduction as accessible as possible, some mathematical details are not specified. We refer the interested reader to the specialized literature on the subject [1,2,10,13]. We consider a control system described by the following differential equation: where is the state of the system at time t, f a smooth function and the control law. We assume here that there is no constraint on the control amplitude. The goal of a control problem is to bring the state of the system from the initial state as close as possible to a target state in a time while minimizing a given cost functional . For a distance to the target state defined by , a standard functional is: where is a positive constant which expresses the relative weight between the distance to the target state and the second term. This latter can be interpreted as the energy of the control. We formulate the optimal control from the Pontryagin Maximum Principle (PMP) which gives necessary conditions for a control solution to be optimal [2,4,9,11]. We introduce the Pontryagin Hamiltonian (the index P corresponds to Pontryagin) as: where is the adjoint state. This state plays qualitatively the role of a Lagrange multiplier for the optimization problem [10,13]. The state and the adjoint state of the dynamics fulfill the Hamilton’s equation: with the following initial and final conditions and , while the optimal control is given by , i.e.,

In the non-linear case, these conditions can be solved only for simple low-dimensional systems [4,11,25] and numerical algorithms are used for more complex dynamics [1,30,55]. For linear systems, the optimal solutions can be derived explicitly if there is no additional constraint on the control field. This approach is known in the literature as Linear Quadratic Optimal Control [2,10,47] and is applied to ICR in Appendix D. When experimental limitations such as maximum pulse intensity are accounted for in the numerical optimization process, the optimal control law is derived numerically from iterative algorithms, which are described in Section 3.2.

3.2. Optimal Gradient-Based Algorithm

The goal of this section is to develop a first-order gradient-based algorithm suited to this control problem [1]. We use a numerical optimization algorithm to take into account field amplitude constraint of the form . Note that this algorithm can be seen as the counterpart of the GRAPE algorithm in NMR [30] and that other limitations such as spectral constraints or bandwidth limitations could be added [56,57,58,59,60,61,62,63]. In the numerical simulations, the control field is described as a piece-wise constant function. Rapid time variations leading to high frequencies may appear in the optimization process. For question of numerical stability and precision, we apply the algorithm in the system with the RWA and then we use the derived control law in the original dynamical system.

We start from the differential system (6) written in the rotating frame for the ion k as: where and . The two controls satisfy the limitation with . The corresponding target state is . We consider a cost functional with no penalty on the control field defined as:

The Pontryagin Hamiltonian can be expressed as:

The adjoint states fulfill the following relations:

The gradients are given by:

The correction to the control fields and at each step of the algorithm is proportional to these gradients [1]. The final adjoint states can be expressed as: and Equation (8) can be directly integrated backward in time. We thus consider the following gradient-based algorithm.

Choose guess fields and .

Propagate forward the state of every ion k and compute .

Propagate backward the adjoint state of the system from Equation (8).

Compute the corrections and to the control fields, , where is a small positive constant.

Define the new control fields .

Truncate the new control fields and to satisfy the constraint :

Go to Step 2 until a given accuracy is reached.

Similar algorithms are used in NMR for taking into account pulse constraints [15,16,17]. Note that the use of a gradient causes this type of algorithm to converge towards a local maximum of the optimization problem. Numerical simulations with different guess fields allow partly overcoming this limitation, even if the global maximum is not reached with certainty. The efficiency of this algorithm in ICR is illustrated numerically in Section 4.

4. Numerical Results

We present numerical results obtained either with LQOCT (see Appendix D for details) or with the gradient-based algorithm. A comparison with the SWIFT approach described in Appendix C is also done. Different experimental constraints have to be satisfied by the control pulse. The objective is to excite ions in a wide range of frequencies around a central frequency of the order of 500 kHz. The excitation has to be as uniform as possible in radius and phase in the range and close to zero outside. As a benchmark example, we choose in this section to consider the interval [400, 600] kHz. Using the linearity of the dynamics, these results can be transposed to another range of frequencies by a scaling of the excitation pulse duration and of the pulse amplitude. For instance, if the total process time is increased by a factor , then the range of frequencies and the amplitude of the electric field are divided by the same parameter . The description of the optimal control of this infinite dimensional dynamical system is mathematically quite intricate, even if some results can be established [47,50]. For practical and numerical reasons, it is more convenient to consider a finite set of N systems by discretizing the frequency interval. We thus consider the simultaneous control of each element of this set. In this paper, we consider a regular discretization, but other choices could be possible, and the frequency step is chosen small enough to avoid the discretization effect. Note that the same approach is used in NMR to control a spin ensemble [15,16,17]. The required sharp excitation profile is modeled by the following function: where μ is a free parameter allowing to adjust the slope of the excitation gate. In the numerical simulations, μ is fixed to 0.1. We impose that the final radius is cm and a final phase varying linearly with the frequency, with a maximum variation of the order of 1/1 Hz, which corresponds to 17.45 rad/kHz. The final phase of the ion k is expressed as , where a is a parameter characterizing the slope of the angular variation. The magnetic field is set to 7 T and the maximum electric field amplitude that can be generated is of the order of V·m. The control time can be very long, of the order of few hundred ms, but more stable numerical results were achieved for duration of the order of few ms.

We first present in Figure 1a series of simulations without limitation on the field intensity. The excitation pulse duration is chosen to be equal to 1 ms. The electric field is computed from a set of the order of 500 frequencies regularly spaced in the interval under study. Since the derived solutions are very sharp, this duration can be modified to some extent without changing the control pulse. Figure 1 compares the results achieved by LQOCT and by the SWIFT approach. The optimal solutions can be computed by using or not the RWA. Note that the pulse computed in the RWA is then applied without any approximation to the original system. In the case displayed in Figure 1, very similar efficiencies are obtained for the two optimal excitations. The optimal pulses have a shape similar to that of the SWIFT solution, even if their amplitudes and durations are different. We recall that the analytical expressions of the pulses are different, but that, for a continuous set of frequencies, the control field is expected to be unique. This statement can be rigorously shown in the case of an ensemble of springs under some specific mathematical assumptions [47,50]. However, the optimal control method offers greater flexibility since one can play with different parameters such as the cost functional or the number of discrete frequencies to adjust the final result.

Figure 1

Comparison between the optimal and the SWIFT approaches for the robust control of an ensemble of ions in the frequency range kHz. The small insert is a zoom of the profile around the frequency kHz. (a,e) The evolution of the final radius and phase as a function of f. Note that an arbitrary constant has been added to the phase in order to superimpose the curves (the three lines are practically indistinguishable in (e)). The black, blue (dark gray) and red (light gray) solid lines depict, respectively, the optimal solutions computed without and with the RWA and the SWIFT pulse. The SWIFT and optimal control laws are plotted in (b,c) (optimal without RWA) and (d) (optimal with RWA). The number of discretized frequency points is set to 601 in the optimization process in the range kHz.

We study in Figure 2 the role of the phase slope of the excitation profile on the structure of the pulse. Figure 2 shows that this slope changes the position of the peak of the pulse. This position can be deduced from a Fourier transform of the profile. Very good results were obtained for slopes in the range [0.05, 0.95] with a maximum pulse amplitude almost constant. Pulse distortion appears when the slope parameter a is close to 0 or 1. For , it becomes very difficult to control all the ions which have to reach a fixed target state in a space-fixed frame, independently of their own frequency. Note that similar results were achieved in NMR [64,65], which highlights the similarities between the control of the two dynamics.

Figure 2

Same as Figure 1 but for different slopes of the excitation profile. The parameter a is fixed, respectively, to 0.25 and 0.75 in (a,b). The optimal pulses without RWA are represented.

We now focus on ion control with amplitude constraint. The numerical simulations were carried out by assuming the RWA. The same set of discretized frequencies is chosen. We optimize piecewise constant functions with a time step lower than 1 μs to avoid discretization effect. The dynamics are integrated numerically through the formulas given in Section 2. More than 1000 iterations are usually needed to converge to an efficient solution. In a first step, we apply the gradient-based algorithm described in Section 3.2 with only one control field, namely , and the phase of the electric field is set to 0. We consider the same control problem as before and the optimal solutions derived above are used as guess field for the optimization algorithm. Figure 3 displays the best result achieved with this limitation. The maximum field amplitude can be reduced from 130 to 100 V·m while maintaining an almost ideal excitation profile. This reduction was made possible by a distribution of the energy along the control interval. Outside of ms, the amplitude of the optimized field is much larger than the one of the guess pulse. As a comparison, Figure 3 also presents the profile obtained from the optimal pulse of Figure 1 whose amplitude has been arbitrarily limited to 100 V·m, showing the non-trivial transformation made by the algorithm.

Figure 3

(a) The final radius r as a function of the frequency f. The black, blue (dark gray) and red (light gray) curves represent, respectively, the ideal profile, the one obtained with the optimization algorithm and the one corresponding to the optimal pulse of Figure 1 whose amplitude has been abruptly limited (see the text for details). Note that the black and blue lines in (a) are almost superimposed. The amplitude of the optimal fields in the rotating frame with (black curve) and without (red or light gray curve) constraints are depicted in (b).

The optimization algorithm fails to converge towards a very good excitation profile, when the maximum amplitude is much smaller than 100 V·m. Note that more advanced optimal algorithms could be tested in this case to verify these convergence issues [66]. This obstacle can be partly overcome by considering in a second step two control fields (in the rotating frame) denoted and . An example is displayed in Figure 4 for a maximum amplitude of 100 and 50 V·m. An almost perfect excitation profile is achieved in these two cases. Note the different structures of the fields along the x- and y-directions, namely even and odd functions. This observation was also made in some optimal control problems in NMR [17].

Figure 4

Plot of the optimal amplitudes (blue or dark gray) and (red or light gray) for a maximum amplitude of 100 V·m (a)) and 50 V·m (b). (c) The corresponding total amplitude .

A systematic analysis of the efficiency of the optimized control fields with respect to the maximum pulse amplitude and to the control duration is provided in Figure 5. The efficiency of the control process is measured from the cost functional given in Equation (7). As could be expected, we observe that better results are achieved for larger maximum amplitude and control time. However, the final fidelity varies in a quite complex way with the control time. A saturation is observed for times of the order of few milliseconds. It is not clear if this point is due to an intrinsic limitation of the control protocol or to convergence problems of the algorithm. Further investigations are needed to clarify this issue.

Figure 5

Evolution of the logarithm of the cost functional as a function of the control time for different maximum amplitudes (black, 30 V·m; blue or dark gray, 50 V·m; red or light gray, 70 V·m).

5. Conclusions

We applied optimal control techniques to the robust excitation of ions in ICR. We considered the simplified but realistic conditions of a two-dimensional trajectory and of a homogeneous magnetic field. In this model system, we propose different ways to solve the optimal control problems. Such methods are directly inspired from NMR in which OCT is a standard and efficient tool. In the case without pulse limitation, the linearity of the dynamical equations allows using LQOCT, which has the advantage to lead to an analytical formula of the control law. Very good results were obtained both for the final radii and phases of the ions. A specific range of frequencies was considered in this study, but the same approach can be extended to broadband excitation from 100 to 900 kHz. However, this solution is both in shape and in amplitude very similar to the SWIFT pulse. The two solutions are expected to be equal for a continuous range of frequencies. More original control laws are derived when the pulse intensity is limited. Due to this constraint, optimal iterative algorithms have to be used, and we adapt to ICR the standard GRAPE algorithm, well-known in NMR. Even if this algorithm has some limitations, it allows reducing the pulse intensity, by a factor larger than three in the examples under study. On the basis of NMR results, this algorithm is expected to be very efficient in the case of other excitation profiles. The very encouraging and promising results obtained in this investigation must now be confirmed by experimental implementation. Numerical simulations of this study are not fully realistic. Effects such as the magnetron motion, field geometry, field inhomogeneities or ion collisions are neglected. However, the model system we consider describes quite faithfully the main cyclotronic behavior and permits to grasp rapidly the main features of ion trajectories. Numerical codes were developed to account for such experimental details. The relative simplicity of the application of numerical optimal algorithms makes it possible to adapt it straightforwardly to a new class of control problems. They could thus be combined with such codes. We are therefore quite confident about the extension of optimization procedures to these additional experimental constraints and limitations. Work is in progress on these different issues.

All pulse shapes are available upon request to the corresponding author.