Optimized Beam Matching Using Extremum Seeking

Proceedings of 2005 Particle Accelerator Conference, Knoxville, Tennessee OPTIMIZED BEAM MATCHING USING EXTREMUM SEEKING Eugenio Schuster ∗ , Christopher K. Allen † , and Miroslav Krstić ‡ Abstract The matching problem for a low energy transport system is approached from a control theoretical viewpoint. The beam dynamics are modeled using the KV (KapchinskijVladimirskij) envelope equations. Multi-Parameter Extremum Seeking, a real-time non-model based optimization technique, is considered in this work for the lens tuning in the beam matching system. Numerical simulations illustrate the effectiveness of this approach. Figure 1: Matching channel. INTRODUCTION In this work we approach the beam matching problem, where the beam must be matched to the acceptance ellipse of an accelerating structure or transport section. Specifically we consider a fixed geometry matching section consisting of four quadrupole lenses. The objective of this system is to take any arbitrary initial beam state and “match” it to the acceptance ellipse of the following section, i.e., any given initial state xini to a prescribed final state xf in , through the control of the lens strengths in the transport (matching) channel. A review of beam transport including matching was recently presented in [1]. We assume the matching channel to be composed of discrete beamline elements, such as lenses, and drifts. These elements are cascade along the beam axis, considered the z axis, to form the matching channel. This configuration is depicted in Figure 1. The input to the lenses, labeled θ1 , θ2 , θ3 , and θ4 in the figure, represent the focusing strength of the lenses and are the parameters of the channel that may be varied. We consider the usage of extremum seeking as optimization technique. The performance of this technique in terms of computational demand, globality, and ability to deal with cost functions full of local minima, is reported. The paper is organized as follows. In Section 2 the optimization problem is defined. Section 3 introduces the fundaments of extremum seeking. The results of the simulation study are presented in Section 4. The paper is closed by a summary in Section 5. Figure 2: Focusing function. envelope equations [2]. Let the z coordinate represent the position along the design trajectory, and thus the xy plane is the transverse plane for the particle beam. At each z location, let X(z) and Y (z) represent the semi-axes of the beam envelope in the x and y planes, respectively. The KV equations then appear as ǫ2 2K − X3 X +Y X 2K ǫ2 Y ′′ + θ(z)Y − − Y3 X +Y Y X ′′ − θ(z)X − Assuming a continuous, elliptically-symmetric particle beam, we model its dynamics using the KV coupled∗ E. Schuster is with the Department of Mechanical Engineering and Mechanics, Lehigh University, 19 Memorial Drive West, Bethlehem, PA 18015-0385, USA, [email protected] † C.K. Allen is with Los Alamos National Laboratory, Los Alamos, NM 87545, USA, [email protected] ‡ M. Krstić is with the Department of Mechanical and Aerospace Engineering, University of California at San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0411, USA, [email protected] 0 (1) = 0 (2) where the prime indicates differentiation with respect to z, that varies from 0 to L. The function θ(z) is the focusing (control) function. K is the beam perveance, ǫX and ǫY are the effective emittances of the beam in the x and y planes, respectively. The focusing function θ(z), shown in Figure 2, can be written as  κθ1 z ∈ [Ld , Ld + Lq ]    κθ2 z ∈ [2Ld + Lq , 2Ld + 2Lq ]  θ(z) = κθ3     κθ4 0 PROBLEM DEFINITION = z ∈ [3Ld + 2Lq , 3Ld + 3Lq ] z ∈ [4Ld + 3Lq , 4Ld + 4Lq ] otherwise (3) where κ is a constant, Ld is the drift length, and Lq is the quadropole lens length. The matching channel parameters θ1 , θ2 , θ3 , and θ4 must satisfy the following constraints: 0 ≤ θ1 , θ3 ≤ 50 and −50 ≤ θ2 , θ4 ≤ 0. We are given initial conditions for the beam at z = 0, the transport system’s entrance location. These conditions characterize the beam coming from the preceding section of the transport or accelerator system.They may be translated into initial conditions for the beam envelopes in the ′ ′ ) and in the y plane (Yini , Yini ). In x plane (Xini , Xini 4269 0-7803-8859-3/05/$20.00 c 2005 IEEE Proceedings of 2005 Particle Accelerator Conference, Knoxville, Tennessee Figure 3: Extremum seeking control scheme. matching systems we are also given desired final conditions, or target conditions, at z = L, the exit location of the matching channel. We denote this target conditions as ′ ′ ) and (Ytar , Ytar ). They are prescribed by the (Xtar , Xtar acceptance requirements of the next section of the transport or accelerator system. T Denoting x = [X X ′ Y Y ′ ] , we define xini = ′ ′ T Yini Yini ] and xf in = x(0) = [Xini Xini T  ′ ′ x(L) = Xf in Xf in Yf in Yf in . In addition, we define a target value for x denoted as xtar = ′ ′ T Ytar Ytar ] , and desired beam profiles [Xtar Xtar (beam trajectories) for X(z) and Y (z) denoted as Xdes (z) and Ydes (z) respectively. Given xini , xtar , Xdes (z) and Ydes (z), we use an extremum seeking procedure to minimize the cost function J given by beam “run.” Thus, we employ the discrete time variant [4] of extremum seeking. The implementation is depicted in Figure 3, where q denotes the variable of the Z-transform. The high-pass filter is designed as 0 < h < 1, and the modulation frequency ω is selected such that ω = απ, 0 < |α| < 1, and α is rational. The static nonlinear block J(θ) corresponds to one run of the KV system. The objective is to minimize J. If J has a global minimum its value is denoted by J ∗ and its argument by θ∗ . In this case we are dealing with a multi-parameter extremum seeking procedure, where the variables are written as       θ̂1 (k) θ1 (k) ξ1 (k)  θ̂2 (k)   ξ2 (k)   θ2 (k)   , θ̂(k) =  θ(k) =    θ̂3 (k)  , ξ(k) =  ξ3 (k)  . θ3 (k) θ4 (k) ξ4 (k) θ̂4 (k) The extremum seeking constants a and b shown in Figure 3 are diagonal matrices of dimension 4. In addition, cos(ωk) and cos(ωk − φ) denote column vectors of dimension 4, where each one of the components has a specific frequency ωi and phase φi , for i = 1 . . . 4. In each iteration of the extremum seeking procedure, θ(k) is used to compute the focusing function θ(z) according to (3), which is in turn fed into the KV equations (1) and (2). Given xini , the KV equations are integrated to obtain X(z), Y (z), and xf in . The output of the nonlinear static map, J(k) = J(θ(k)), is then obtained by evaluating (4) and used to compute θ(k + 1) according to the extremum seeking procedure in Figure 3, or written equivalently as 1 J ={k1 J1 + k2 J2 + k3 J3 } 2 (4) Jf (k) J1 =KX (Xf inal − Xtarget )2 + KY (Yf inal − Ytarget )2 (5) 2 2 ′ ′ J2 =KdX Xf′ inal − Xtarget +KdY Yf′inal − Ytarget (6) L   w(z) KiX (X(z) − Xdes (z))2 J3 = = −hJf (k − 1) + J(k) − J(k − 1) (8) ξ(k) = Jf (k)b cos(ωk − φ) θ̂(k + 1) = θ̂(k) − γξ(k) (10) θ(k + 1) = θ̂(k + 1) + a cos(ω(k + 1)). (11) (9) 0 +KiY (Y (z) − Ydes (z))2 dz, where k1 , k2 , k3 , KX , KY , KdX , KdY , KiX , and KiY are weight constants, and wz is an integral weight function. EXTREMUM SEEKING Extremum seeking control, a popular tool in control applications in the 1940-50’s, has seen a resurgence in popularity as a real time optimization tool in different fields of engineering [3]. Extremum seeking is applicable in situations where there is a nonlinearity in the control problem, and the nonlinearity has a local minimum or a maximum. The parameter space can be multivariable. In this paper we use extremum seeking for iterative optimization of θ to make xf in as close as we can to xtar . For each new value of θ we run the KV equations and obtain xf in . We point out that, since xtar is given arbitrarily, xf in is obtained via solving a system of nonlinear differential equations, and the lens input applied through θ is highly constrained in its waveform, there may not exist θ such that perfect matching is achieved, xf in = xtar , thus we try to obtain the best possible approximate matching. We change θ after each 0-7803-8859-3/05/$20.00 c 2005 IEEE SIMULATION RESULTS (7) The physical parameters used in the simulations presented in this section are K = 2.7932 × 10−6 , ǫX = 6 × 10−6 , ǫY = 6 × 10−6 , κ = 2.6689, Ld = 0.1488, Lq = 0.0610, and L = 0.988. In addition, the extremum seeking parameters are h = 0.4, ω1 = ωbase × π, 2 3 4 ω2 = ωbase × π, ω3 = ωbase × π, ω4 = ωbase × π, where M (ω1 ) (ω1 ) ωbase = 0.95, γ1 = 0.1, γ2 = 0.1 M (ω2 ) , γ3 = 0.1 M M (ω3 ) , (ω1 ) γ4 = 0.1 M M (ω4 ) , and φ1 = −φ(ω1 ), φ2 = −φ(ω2 ), φ3 = −φ(ω3 ), φ4 = −φ(ω4 ), where M (ω) and φ(ω) are respectively the magnitude and phase of the frequency response of the high-pass filter. For all the simulations, the initial condition of the beam at the entrance of the channel is     0.00109372 0.00147377  −0.00786479   −0.00601315  . (12) ,x = xini = 0.00201395  tar  0.00328979  0.01172626 0.00768638 T Denoting θ = [θ1 θ2 θ3 θ4 ] , the initial condition for the extremum seeking parameters in all the simulations is T equal to θ(0) = [25, −25, 25, −25] . 4270 Proceedings of 2005 Particle Accelerator Conference, Knoxville, Tennessee Parameter Evolution Cost Function 5 40 4.5 −15 30 4 20 3.5 Filtered θ J [dB] −20 −25 −30 10 −10 x 10 X Y Xtar Ytar X,Ydes 3 θ1 θ2 θ3 θ4 0 2.5 2 −20 −35 Beam profile for estimated θ −3 50 −10 1.5 −30 1 −40 0.5 −40 −45 50 100 150 200 250 300 350 400 450 500 −50 50 100 150 200 Iteration (a) KX = 2000, KY = 2000, KdX = 1, KdY = 1, KiX = KiY = 50, k1 = 1, k2 = 1, k3 = 1. (13) The integral weight w(z) is equal to 0 for 0 ≤ z < 0.2, 2 for 0.2 ≤ z < 0.8, 50 for 0.8 ≤ z < 0.9, and 100 for 0.9 ≤ z < L. We try not only to match the final section of the beam profile but also to reduce excursions in the middle section. In this case the amplitudes of the sinusoidal excitations are varied according to the following law: = [2.25, 2, 1.75, 1.5] = [0.25, 0.75, 0.5, 0.5] = [0.1, 0.25, 0.1, 0.1] = [0.05, 0.05, 0.05, 0.05] = [0.01, 0.01, 0.025, 0.025] if if if if if −35dB ≤ J −45dB ≤ J < −35dB −46dB ≤ J < −45dB −50dB ≤ J < −46dB J < −50dB Comparing xf in (below) with xtar , we can note that we do have an acceptable matching of the target conditions. In this case we are converging to θ̂f in (below).     0.00109266 34.855771  −30.710796   −0.00734280  . , θ̂ = xf in = 0.00327968  f in  14.736266  0.01063013 300 350 400 450 500 0 0 0.2 (14) −30.669087 The time evolution of θ1 , θ2 , θ3 , θ4 in Figure 4-b shows that a steady value is reached after 150 iterations. This can be also noted from Figure 4-a, where the cost function does reach a steady value after 150 iterations, showing a very fast convergence. Figure 4-c shows the beam profile for θ̂f in . Not only the matching of the target conditions is acceptable, but also the matching of the desired profile. This is explained by how the cost function was defined. When k3 = 0, the global minimum of the functional J corresponds to θ∗ = [ 38 −38 38 −38 ]T , giving a perfect matching for the target conditions. In this case, we achieve very similar final conditions reducing at the same time the excursion of X(z) and Y (z). It is interesting to 0.4 0.6 0.8 z (b) Figure 4: Cost function (a), θ evolution (b), beam profile for θf inal (c). In this case we take Xdes (z) and Ydes (z) as a combination of two linear functions as shown in Figure 4-c (dotted line). The slope of the last section of the desired beam profile coincides with the target conditions for the derivatives in order to facilitate their matching. The use of only one linear function, connecting Xini and Yini , with Xtar and Ytar respectively, would be in conflict with the target conditions for the derivatives. Figures 4 shows the extremum seeking results when the cost function parameters are given by a1,2,3,4 a1,2,3,4 a1,2,3,4 a1,2,3,4 a1,2,3,4 250 Iteration (c) note how different is the value of θ̂f in (above) from the global minimum θ∗ and at the same time how acceptable is the matching of the target conditions. This is a sample of how complex the functional map J(θ) is. An in-depth analysis of the functional map shows that we are dealing with a spiky cost function with numerous local minima. CONCLUSIONS A multi-parameter extremum seeking procedure has been implemented, and successfully tested in simulations, for the tuning of the lens strengths in a 4-lens matching channel. Although the scheme shows a very fast convergence, considering that we are tuning simultaneously four parameters, it is not clear at this point if the scheme can be used for real-time optimization. However, based on the promising results obtained in the simulation study, it is anticipated that the scheme can play an important role in an off-line design process. A rigorous performance comparison with other methods previously used is still pending and will be part of our future work in the short term. We must highlight at this point the capability of the scheme of avoiding getting stuck in local minima with relatively large values of the cost function. The modification of the amplitude of the sinusoidal excitation as a function of the value of the cost function is key in this achievement. This suggests the possibility of designing an extremum seeking scheme that automatically adapts their gains or sinusoidal amplitudes to permanently seek a lower value of the cost function. This potential scheme would be very useful for applications with spiky cost function maps as the one considered in this work. 4271 REFERENCES [1] S.M. Lund and B. Bukh, “Stability properties of the transverse envelope equations describing intense ion beam transport,” Phys. Rev. ST Accelerators and Beams, vol. 7, 024801, 2004. [2] I.M. Kapchinskij and V.V. Vladimirskij, Proc. Int. Conf. on HighEnergy Accelerators and Instrumentation, CERN, 1959, pp. 274288. [3] K. Ariyur and M. Krstic, Real-Time Optimization by Extremum Seeking Feedback, Wiley, 2003. [4] J.-Y. Choi, M. Krstic, K. Ariyur and J.S. Lee, “Extremum seeking control for discrete-time systems,” IEEE Transactions on Automatic Control, vol. 47, no. 2, pp. 318-323, 2002. 0-7803-8859-3/05/$20.00 c 2005 IEEE