Adaptive robust control of longitudinal and transverse electron beam profiles

Feedback control of the longitudinal and transverse electron beam profiles are considered to be critical for beam control in accelerators. In the feedback scheme, the longitudinal or transverse beam profile is measured and compared to a desired profile to give an error estimate. The error is then used to act on the appropriate actuators to correct the profile. The role of the transverse feedback is to steer the beam in a particular trajectory, known as the “ orbit. ” The common approach for orbit correction is based on approximately inverting the response matrix, and in the best case, involves regulating or filtering the singular values. In the current contribution, a more systematic and structured way of handling orbit correction is introduced giving robustness against uncertainties in the response matrix. Moreover, the input bounds are treated to avoid violating the limits of the corrector currents. The concept of the robust orbit correction has been successfully tested at the SwissFEL injector test facility. In the SwissFEL machine, a photo-injector laser system extracts electrons from a cathode and a similar robust control method is developed for the longitudinal feedback control of the current profile of the electron bunch. The method manipulates the angles of the crystals in the laser system to produce a desired charge distribution over the electron bunch length. This approach paves the way towards automation of laser pulse stacking.


I. INTRODUCTION
Many orbit correction algorithms have been devised for synchrotron light sources. The task of orbit correction is to determine a suitable set of currents in the corrector magnets in order to keep the beam on a specified trajectory. The beam position is measured at several locations of the trajectory through "beam position monitors" (BPMs). A static response matrix maps the changes in corrector currents to the changes in the equilibrium position of the beam at each BPM. In principle, orbit correction attempts to invert the response matrix to translate the beam positions into the actuation on each of the corrector magnets. This procedure is usually based on the singular value decomposition (SVD) method which is already in use at many accelerator laboratories over the world [1][2][3][4][5][6][7][8]. In some cases, an adaptive approach was proposed to estimate the response matrix [9], and some also took the corrector's limits into account in the SVD algorithm [10].
Both in theory and experiment one can think of two types of orbit correction algorithms: static orbit correction and dynamic orbit correction [11]. The objective of the static orbit correction is to minimize the steady-state orbit deviation from the desired trajectory. The process runs slowly and the feedback rate is normally approximately 1 Hz. In dynamic orbit correction the objective is to minimize the beam motion on some frequency range. The process rate is of the order of several kHz which may take the corrector magnets' dynamics into account. In this work, we concentrate only on the global steady-state orbit correction with a static response matrix.
The current contribution develops an adaptive robust control scheme that captures uncertainties of the response matrix, and treats the input limits of the system, appropriately. In addition, the system response matrix is estimated on-the-fly, and therefore it can be applied to time-varying systems ( [12], Chap. 6). This approach has been successfully tested at the SwissFEL [13] injector test facility on 10 sets of BPM and corrector magnets. This method may be applied to different systems and in this paper, we develop an automated procedure for the laser pulse stacking system, using the same mathematical approach, to achieve a desired charge distribution over the bunch length. The proposed automation facilitates the laser tuning during the operation.
The paper is organized as follows. Section II defines the problem formulation. In Sec. III, an adaptive robust control scheme is introduced for a multiple-input multiple-output (MIMO) system with a static response matrix. Sections IV and V demonstrate two illustrative applications of the proposed method: The orbit correction (transverse feedback) and the laser pulse stacking system (longitudinal feedback).

II. PROBLEM FORMULATION
In this section and the section that follows, we derive an optimal control algorithm for a general static system with multiple inputs and multiple outputs. The algorithm is then applied on two different systems, namely, the orbit correction and the laser pulse stacking systems. For the orbit correction configuration, the inputs and outputs are, respectively, the input currents to the corrector magnets and the beam position measured at the BPMs. For the laser system, the crystals' angles and the current profile are considered as the system inputs and outputs, respectively. The two systems are described in more detail in the corresponding sections. The control method is iterative; the measured output in one run is used to update the input for the next run.
Consider the following nonlinear MIMO system with u ∈ R m and y ∈ R n denoting the input and output, respectively.
where f is a static nonlinear function, uniquely mapping u to y. We assume at least as many measurements as actuators (n ≥ m).
The system can be linearized around the operating point, y 0 , with the associated input, u 0 : whereR ∈ R n×m is the first order approximation of f. The matrixR is commonly referred to as the system response matrix. We model the linear term with matrix R which might differ from the actual response matrix,R. The difference between the two response matrices can be expressed in a multiplicative manner. That is, where Δ ∈ R m×m is an initially unknown perturbation matrix capturing the uncertainty in our model of the system. Even though the perturbation matrix, Δ, is not fully known, we can restrict the problem to uncertainties up to a certain size. In other words, we assume that Δ is an unknown matrix with the following norm upper-bound, For the sake of generality, we do not make any additional assumptions on the structure of Δ at this stage. Let us define y d as the desired output vector which is assumed constant in time. This can indicate the desired beam position at the BPMs, or the desired current profile in the laser problem. The objective is to find the appropriate actuation, u, to generate the desired output and, as it was already pointed out, this is done in an iterative manner.
Since the perturbation, Δ, is initially not known, we model the initial output response as, whereŷ kþ1 is the expected output at the next iteration, and where subscript k is the iteration counter.
Calculating the least-square solution for y d ¼ŷ kþ1 in (5), gives us: under the assumption that ðR T RÞ −1 exists. Therefore, the first order approximation of output vector at the next iteration would be, whereR represents the actual system response matrix which might be different from R. The error is considered as the discrepancy between the measured output and the desired output. From (7) and (3) the linearized error transmission matrix can be readily calculated: Note that if ‖MðΔÞ‖ < 1 (for any induced norm), the error asymptotically decays to zero as the iteration proceeds, i.e., e k → 0 as k → ∞. We now show that, with the modeled uncertainty, the spectral norm of MðΔÞ can potentially exceed one. For simplicity in the analysis, we assume that R is an invertible (and thus square) matrix. Therefore, the MðΔÞ simplifies to To observe the stability of MðΔÞ, we calculate the spectral radius (largest absolute eigenvalue): ρðMðΔÞÞ ≔ maxfjλ 1 j; …; jλ n jg: In general, the spectral radius provides a lower bound on the 2-norm, and we have where the second inequality in (10) follows from ‖AB‖ 2 ≤ ‖A‖ 2 ‖B‖ 2 . Now, taking the largest spectral radius over all allowed Δ's (i.e., ‖Δ‖ 2 ≤ γ), leads us to where κ ≔ ‖R‖ 2 ‖R −1 ‖ 2 denotes the condition number of the response matrix, R. The condition number of a matrix is the ratio of its largest to smallest singular value which gives an indication how spread out its singular values are. Inequality (11) gives an upper-bound for the worst-case spectral radius of MðΔÞ. Now, we show that the equality can actually happen if Δ is a full matrix without any structure [14]. To show this, define Z ≔ 1 Note that ‖Z‖ 2 ≤ 1. Select Z ¼ VU T , where U and V are, respectively, the left and right matrices of the singular value decomposition of R ¼ UΣV T . Then we have ‖Z‖ 2 ¼ 1 where the one before last equality follows since U T ¼ U −1 and eigenvalues are unchanged under similarity transformation. Therefore, for an invertible response matrix with full unstructured perturbation matrix, we have Equation (13) implies that for a response matrix with a large condition number, such that γκ is greater than one, the error starts to grow. This illustrates the reason why a robust control scheme is needed to handle the uncertainties properly. This will be discussed in the following section.

III. ROBUST OPTIMAL CONTROL SCHEME
The robust control design for systems with uncertainty in the model, has been fully described in [15] and need not be discussed further here. Also in [16], the authors studied the robust optimal control for uncertain static systems, similar to the system described by (2). In the current approach, the optimal control is derived at each iteration and the constraints on the input are respected. These constraints are mathematically expressed as: The constraint (14a) is an element-wise inequality, i.e., each actuator input must lie within a bound to avoid hitting the actuator limit. The second quadratic constraint given in (14b) restricts the input changes at each iteration to be less than a constant scalar denoted by η.
The algorithm is split in three steps. In the first step, referred to as "initialization," the input is derived from an optimization problem which minimizes the norm of the error for the worst case uncertainty in the model. Once this input is applied to the system and the measurements are captured, another optimization is performed to given an estimation of the uncertainty matrix, Δ. This is discussed in the step "uncertainty estimation." Estimating the perturbation matrix, leads us to the third step to find the optimal input for the next iteration. After the initialization, steps 2 and 3 are repeated subsequently until the convergence is achieved. The three steps in the optimal control algorithm are explained as follows.

A. Initialization
Let us assume that we are at iteration k ¼ 0, with the objective to achieving the desired output vector, y d , at the next iteration. The initial input is conservatively determined through the following optimization problem which minimizes the tracking error for the worst case perturbation: The cost function defined in (15), has an upper limit as follows: Therefore, our task reduces to finding the robust optimal input, δu 0 , through the following optimization problem which minimizes the upper bound on the error norm:

B. Perturbation estimation
Solving (17) gives the initial input to the system, i.e., u 1 ¼ u 0 þ δu 0 , which results in y 1 at the output. For iteration k ≥ 1, the perturbation, Δ, can be estimated iteratively using the input-output data. On this basis, another optimization is performed to minimize the difference between the measurement and the model estimate, Remark.-In the uncertainty estimation given by (18), we use only the last input-output data and the previous measurements and inputs are not considered. This is mainly to better estimate the system around the last operating point. However, the solution to (18) may not be unique. To see this, we reformulate the cost function in the following way: where ⊗ is the Kronecker product, and where vecðΔÞ is a column vector obtained by stacking the columns of Δ. For simplicity in the notation, we defineũ k ≔ u k − u 0 . The matrix R ⊗ũ T k is a matrix of size n × m 2 , which is made up of m different scalings of R. This implies that, in the best case, the rank of R ⊗ũ T k is m whereas vecðΔÞ ∈ R m 2 , and therefore, the optimal solution will not be unique. However, a typical solution is given by taking the pseudoinverse of R ⊗ũ T k , as follows, where † denotes the pseudoinverse, and the constraint on the matrix norm is assumed inactive.
To have a unique minimizer, several experiments should be run to collect enough input-output data. Assume that N experiments are performed, and the minimization in (19) is calculated over the complete data set. We define Φ k , a regressor, to be: Therefore, the original optimization problem (19) converts to: where Y k is given as follows, . . .
To have Φ k to be of full column rank, we need to run N ≥ m 2 n experiments. On the other hand, as the iteration proceeds, the input vectors converge, i.e., u k → u ∞ , which eventually reduces the rank of Φ k to m. Hence, in this paper, the uncertainty is estimated in an optimistic way.
▪ The perturbation matrix can be restricted to a special form or structure such as diagonal, Diagonal perturbation matrices indicate that the uncertainties are modeled only on the actuators. For the orbit correction problem, we will use a diagonal Δ to assign the uncertainty only to the currents in the corrector magnets.
The form of the optimization given in (18) implies that the perturbation matrix, Δ, is estimated only at each time using the last input-output data. However, this may not give us complete information about the uncertainty as it is an optimistic estimation. An improved version of the uncertainty estimation procedure is discussed in [12]. From the estimated Δ, the response matrix can be updated in two different ways: adaptive and nonadaptive. Figure 1 schematically illustrates the concept of the two ways of updating the response matrix. In the nonadaptive way [see Fig. 1], the uncertainty is always added to the initial response matrix, i.e., And since the perturbation matrix has a norm upper bound of γ, it lies in a ball with radius γ‖R 0 ‖. Moreover, the center of the ball is always fixed at R 0 . However, in the adaptive manner [see Fig. 1(b)], the response matrix is recursively updated as and therefore, the center of the ball moves. The maximum radius of the ball is given by γ‖R k ‖. In the adaptive case the norm of the perturbation matrix can shrink as iteration proceeds, depending on the applied inputs. More details on this matter are given in ( [12], Chap. 6).

C. Finding the optimal input
Estimating the perturbation Δ k from (18), leads us to the next step, which is to determine the optimal control input for k ≥ 1: The optimal solution of (27) is then used to update the input, and apply it to the system. The steps defined in (III B) and (III C) are repeated iteratively until convergence is reached. In the following two sections, we apply the adaptive robust method to control the transverse and the longitudinal profiles of the electron beam in a linear accelerator.

IV. ROBUST ORBIT CORRECTION
In the global orbit correction problem, with many sensors and actuators, the response matrix dimension can be large. This might result in an ill-conditioned response matrix with a very wide range of singular values. This makes the system response very sensitive to unstructured uncertainties, especially if an inverse-based controller is used to correct the orbits. Figure 2 depicts a schematic of the orbit correction system in a linear accelerator. The focusing and defocusing quadruples are omitted to reduce complexity of the figure. The beam position (in the transverse plane of x and y) is measured by the beam position monitors (BPMs) as the beam passes through. The BPMs are distributed between several corrector magnets. These correctors are used to influence the beam trajectory in x and y axes through the currents, I x and I y in the magnet.
The response matrix of the system is measured by perturbing the currents in the magnets and measuring the effects on the BPMs. For BPM ith's x-position change the model is of the form, where N denotes the number of correctors which we assume is equal to the number of BPMs. Note that we assume a general case that the model includes coupling of the vertical correction and the horizontal orbit. From (29), we can construct the response matrix relating the current changes to the beam position changes for x and y direction, 0   In the nonadaptive case, the estimated response matrix remains within a range from the initial model, R 0 . In the adaptive case, the center of the ball moves iteratively and the radius shrinks as the norm of the perturbation matrix, Δ k , decreases at every iteration. plotted in Fig. 4. For this experiment, 10 BPMs were used. As we can see, the matrix is close to lower triangular which comes from the fact that the downstream correctors have no impact on the upstream BPMs in a linear accelerator. The matrix has a checkered pattern which indicates that x and y directions are decoupled.

A. Experimental results at the SwissFEL Injector Test Facility
This experiment demonstrates robust orbit correction using the algorithm introduced in the previous section. In the robust approach, the control inputs are determined in a structured way, along with the handling of the actuator limits. In this experiment the optimization schemes, defined in (18) and (27), are used to generate the optimal currents into the corrector magnets. For the optimal control input determination, the cost function is slightly modified to avoid getting close to the actuators limits. The optimization problems in (17) and (27) where the pðu k þ δu k Þ is a penalty, added to keep a safe distance from the actuator limits. In this experiment, we choose where α is a constant. Similar expressions can be employed to make pðuÞ more sensitive to values close to the limits and less sensitive for actuation far from the limits. Since the orbit system response is in a triangular form, any small perturbation in the corrector currents affect all downstream BPMs. Therefore, we model the uncertainty with perturbations on the inputs, i.e., on the currents in the corrector magnets. As a result, the perturbation matrix, Δ, is diagonal. The algorithms are implemented in Matlab scripts using the CVX convex optimization solver [17,18]. The communication to the machine has a control protocol with a roundtrip communication period of approximately 1 second. Figure 5 illustrates the experimental results of robust orbit correction at the SwissFEL injector test facility using 10 BPMs. The beam is brought to the origin after approximately 20 seconds. The beam position is initially off-axis over the whole BPM chain. In this experiment all BPMs are equally weighted, however one can assign different weights on different locations of the beam line. Figure 6    The saturating limit of the corrector magnets is u max ¼ 10A. The penalty function pðuÞ, introduced in the objective function in (31), prevents the inputs to get close to the saturation. Correctors 9 and 10 (close to the downstream end of the beamline) came close to the saturation limit.

V. LASER PULSE STACKING AUTOMATION
This section demonstrates another application of the proposed adaptive robust algorithm in accelerator physics. In this section, the focus is on the longitudinal profile of the beam. The system to be investigated is the "laser pulse stacker" used in the SwissFEL test machine. The production of high-brightness electron beams is challenging as it introduces high demands such as a flat-topped pulse shape on the driving laser pulse applied to a copper photocathode. The following section describes the pulse stacking system in more detail.

A. Laser pulse stacking
A technique for flat-topped pulse generation relies on stacking several laser pulses in time [19]. In this scheme, up to six crystals are used in sequence, which produce up to 64 individual pulses, each with a different delay, depending on the thickness and the group birefringence of the crystals. The combination of those pulses leads to a long quasi flattopped pulse in the temporal domain [20]. This laser pulse is then used to extract the electrons from the cathode.
In particular, a laser pulse stacker is composed of N birefringent crystals that transform a single Gaussian pulse into a stack of 2 N Gaussian pulses that overlap and produce a long flat-topped pulse [19]. Figure 7 illustrates the concept of laser pulse stacking for 2 crystals. The initial single pulse, which is linearly polarized in vertical direction, is passed through the first crystal with the optic axis tilted by some degrees relative to the vertical axis. The crystal splits the input pulse into two pulses separated in time by Δt 1 . The delay between the two pulses is constant and depends on the characteristics of the crystal such as length and refractive indices. The relative intensity of the two laser pulses can be adjusted by rotating the optic axis. The two resulting pulses from crystal 1 are oriented at different angles to the vertical. The two intermediate pulses then pass through the second crystal to be each divided into 2 pulses separated by Δt 2 . The whole process is repeated by adding more crystals so that with N crystals the initial single pulse is divided into 2 N overlapped pulses.

B. Longitudinal feedback on the current profile
The objective of this section is to generate a desired beam current profile. To measure the current pulse, the beam is deflected by a transverse deflecting cavity and a camera takes an image of the longitudinal profile which is then processed to extract the charge distribution over the bunch length [21]. Figure 8 illustrates an example of the beam current profile. The pulse shape is influenced mainly by the radio frequency (rf) accelerating cavities and the laser stacker system. In this control scheme, we modify the pulse shape by manipulating the crystals' angles which, in turn, changes the relative intensities of the stacked pulses in the laser profile. In order to generate a flat-topped current pulse, the flat-topped region of the profile is divided into several cells. The beam current is averaged at each cell to The arrows represent the polarization axis. The initial single pulse is passed through 2 birefringent crystals to produce 4 stacked pulses. The first crystal with length L 1 splits the input pulse into two intermediate pulses, separated by Δt 1 . The resulting pulses are then passed through the second crystal with length L 2 to be further split and delayed by Δt 2 . The process is repeated to finally produce 2 N overlapped pulses by passing the initial pulse through N crystals.
give an estimate of flatness error which is then used to steer the crystals (see Fig. 8). In the SwissFEL test injector machine, four crystals in the laser stacker system are motorized so the angles can be set remotely. We can define a response matrix that relates the changes in crystal angles to the changes in the beam intensity at each cell. Dividing the flat-topped region into more cells gives more information about the pulse shape, however since the number of actuators is limited (in this case 4 crystals), this may lead to an under actuation situation. Therefore, the number of cells in the flat-topped region is determined empirically.

C. Experimental results at the SwissFEL injector test facility
The same algorithm described in Sec. III is applied to the laser stacker system to tune the crystal angles by directly measuring the current profile of the beam. In this experiment, the flat-topped region is split into 16 cells and 4 crystals are used to modify the pulse shape. Therefore, the response matrix is a tall matrix of size 16 by 4, which is initially estimated by perturbing the angles and measuring all 16 cells. This response matrix is used to correct the pulse shape based on the adaptive robust algorithm introduced in Sec. III. In this case, no special form is assumed for the perturbation matrix, i.e., Δ is a full (or unstructured) matrix. Figure 8 illustrates the final pulse shape after 50 iterations compared to the initial one. The processing time of each iteration is approximately 2 seconds which is determined by the computation time and speed of the crystal motors. Figure 9 shows the standard deviation of the flat-topped region as a criterion of flatness. The relative standard deviation is improved by a factor of 3, starting from an arbitrary beam current profile.
The response matrix, R, is updated iteratively in an adaptive manner as previously described in (26): where Δ k is the perturbation matrix estimated in (18). Figure 9 also plots the 2-norm of Δ k vs iteration index. The norm, ‖Δ k ‖ 2 , decreases to the level of the residual noise. In this experiment, γ is set to 0.3.
The concept of cell-splitting of the current profile is illustrated more clearly in Fig. 10. The current at each cell is averaged and subtracted from the mean of the flat-topped region. This gives an output vector of length 16 which represents the error in the flatness. The error vector is then used to steer the 4 crystals through the optimization defined in Sec. III. Figure 11 compares the longitudinal profile of the electron beam, captured on the transverse profile monitor of the deflecting cavity, before and after applying the algorithm. The final image is taken after 50 iterations which shows that the charge is more uniformly distributed   The relative standard deviation of the flat-topped region (in blue) and norm of the perturbation update matrix, Δ k , (in green) vs iteration index. The response matrix is updated iteratively according to the following update law: over the bunch length. The maximum rotation angle of crystals (with respect to their initial state) is approximately 4 degrees.

VI. CONCLUSION
This research demonstrates two applications of an adaptive robust control method that can be beneficial for both transverse and longitudinal profile correction of the electron beam. Each of the two systems studied in this work can be modeled as a static response matrix which relates the actuation changes to the output changes. The control objective is to generate a desired output vector or profile by manipulating the input actuators without violating the input limits. The control approach is based on an adaptive robust method which seeks to minimize the 2-norm of the error under the assumption that the model is not exact and may include some bounded uncertainty. The algorithm estimates the uncertainty iteratively and finds the optimal control input to achieve the desired trajectory at the output such that the actuator constraints are respected. In the first system, i.e., orbit correction, the objective is to have the beam at zero (or any other desired) position in all BPMs. The corrector magnets are used as the actuators with predefined limits on the input currents. The experimental results at the SwissFEL injector test facility show that the beam position is brought to the origin after about 20 iterations without violating the input constraints. A similar method was applied to the laser pulse stacker system to produce a flat-topped beam current profile. This application facilitates the laser pulse tuning in a sense that it gives the possibility to the operators to automatically adjust the crystal angles within a limited amount of time. The proposed control approach can be extended to other systems that deal with response matrices with a high condition number.  11. The longitudinal profile of the deflected beam, measured at the transverse profile monitor. The automatic laser pulse stacking procedure, generates more uniform charge distribution by controlling the beam current pulse shape.