Symmetry Exploitation in Orbit Feedback Systems of Synchrotron Storage Rings
Idris Kempf, Paul J. Goulart, Stephen R. Duncan, Guenther Rehm
aa r X i v : . [ phy s i c s . acc - ph ] A ug Symmetry Exploitation in Orbit FeedbackSystems of Synchrotron Storage Rings
Idris Kempf, Paul J. Goulart, Stephen R. Duncan, and Guenther Rehm
Abstract
Structural symmetries in the storage ring of synchrotrons are intentionally created duringthe design phase of the magnetic lattices, but they are not considered in the design of controlalgorithms that stabilize the beam of accelerated particles. The choice of control algorithm,however, is limited by the speed requirements of the synchrotron. Standard control algorithmsfor synchrotrons are based on a singular value decomposition (SVD) of the orbit responsematrix. SVD controllers neither exploit the structural symmetries nor exhibit any speedadvantages. Based on the periodicity and the reflection properties of the betatron function,we show that these structural symmetries are inherited by the orbit response matrix. Weshow that the resulting block-circulant and centrosymmetric properties of the matrix canbe used for different computationally efficient decompositions of the controller. We alsoaddress the case of broken symmetry due to odd placements of magnets and monitors.Our efficient decomposition could enable the use of more advanced control techniques forsynchrotrons, such as control algorithms that require real-time optimization. These advancedcontrol techniques could in turn increase the quality of research in synchrotron light sources.
Index Terms
Orbit Feedback, Synchrotron, Symmetries
I. I
NTRODUCTION I N most synchrotrons the magnetic lattices, the beam position monitors and the correctormagnets are placed in repeated patterns around the storage ring [1, Ch. 10.2.4, p. 329].These repeated sections are usually referred to as superperiods or cells , and their patterninvokes a circulant symmetry. Often an additional symmetry is introduced by mirror-reflectingthe pattern in the middle of the storage ring, which invokes a centrosymmetry. The circulantpattern considerably simplifies the design of the synchrotron, while the mirror-reflectioncancels out non-linear effects introduced by quadrupole and sextupole magnets. Althoughthis symmetry is intentionally created in the design phase of the synchrotron, it is most oftennot considered during the synthesis of the orbit feedback system. The feedback system isdesigned around the orbit response matrix R ∈ R SN B × SN C , where S is the number of cellsand N B and N C the number of monitors and corrector magnets per cell, respectively, that The research leading to these results is supported by Diamond Light Source and the Engineering and PhysicalSciences Research Council (EPSRC) with a CASE studentship.I. Kempf, P. Goulart and S. Duncan are with the Department of Engineering Science, University of Oxford,Parks Road, Oxford, OX1 3PJ, UK; (e-mail: { idris.kempf,paul.goulart,stephen.duncan } @eng.ox.ac.uk).G. Rehm was with Diamond Light Source, Didcot, OX11 0DE, UK; He is now with BESSY, Berlin, 12489,Germany (e-mail: [email protected]). relates the corrector magnet inputs u ( z − ) to the horizontal (or vertical) trajectory error ofthe electron beam y ( z − ) as y ( z − ) = R g ( z − ) u ( z − ) + d ( z − ) , (1)where the scalar transfer function g ( z − ) represents the temporal dynamics of the correctormagnets, which are assumed to be identical for all magnets, d ( z − ) the disturbances acting onthe electron beam, such as vibrations transmitted through the girders, and z the Z -transformvariable, respectively. The element on row m and column n , R m,n , of R is characterized bythe betatron function β ∶ R → R + , where L denotes the circumference of the orbit, and givenby [2, eq. (2)] R m,n = √ β Bm β Cn ( πQ ) cos ( πQ − ∣ φ Bm − φ Cn ∣) , (2)where β Xk ∶= β ( s Xk ) , φ Xk ∶= φ ( s Xk ) with s representing distance along the storage ring (startingfrom an arbitrary reference point) and B and C refer to beam position monitors and correctormagnets, respectively. The phase advance φ ∶ R → R + is defined as φ ( s x ) ∶= ∫ s x β − ( s ) ds, (3)and, for a stable electron beam, the betatron tune Q ∶= φ ( L )/ π is always a fractionalnumber [3]. The matrix R typically has a few hundred rows (monitors) and a few hundredcolumns (corrector magnets) and the feedback loop is operated at a frequency of - kHz. The magnetic lattices guide and confine the electron beam around the ring, while theorbit feedback system reduces the trajectory error to a few micrometers. The trajectory errormust be minimized in order to retain certain properties of the synchrotron light that is usedin the beamlines, which are end-stations that use the light for various kinds of experiments,particularly in the X-ray region of the electromagnetic spectrum [4].Most orbit feedback systems use a singular value decomposition (SVD) of R to synthesizean orbit controller, such as in [5] for the European Synchrotron Radiation Facility, [6] forDiamond Light Source, [7] for the MAX IV synchrotron or [8] for the Advanced PhotonSource. In these case, the control input is calculated from u ( z − ) = − K c ( z − ) y ( z − ) , where c ( z − ) is a scalar transfer function that accounts for the dynamics g ( z − ) , and the gain matrix K is calculated from [9, eq. (5)] K = V ( Σ + µ I ) − Σ U ∗ , with R = U [ Σ 0 ] [ V ∗ V ∗ ] , (4)where the zeros in the SVD account for the case that there are usually more actuators thanmonitors, i.e. SN B < SN C . The gain matrix K is a pseudo-inverse of R with a regular-ization parameter µ ∈ R + that accounts for the large difference between the minimum andmaximum non-zero singular value of R , which is common to orbit response matrices ofsynchrotrons [10]. The SVD approach is applicable to any kind of system but does not exploitthe advantages provided by the symmetric structure. Symmetry is always accompanied bycertain redundancies in the mathematical representation of the system [11], and considering thesymmetry speeds-up the controller computations and reduces the memory requirements [12]. Italso benefits the parameter identification [13] and the modeling of parameter uncertainty [10].There may be different reasons that feedback systems have not generally exploited theseexisting symmetries (one exception is [14]). Firstly, the matrices involved in the feedbacksystems have around 100,000 elements and recognizing these symmetries by by inspection is not straightforward. Secondly, the symmetry is often broken by space constraints or by singularcomponents in the storage ring, such as the injection device. Although for this case an SVDseems to be advantageous because it is not reliant on symmetry requirements, it introducesadditional difficulties when the controller robustness is verified with respect to parameteruncertainty [10]. Finally, symmetric decompositions are more difficult to find and in contrastto SVDs, which are supported by most numerical software, no widespread implementation ofan algorithm for symmetric decompositions exists; symmetric decompositions are often foundthrough prior knowledge of the system structure [11].This paper aims to close the gap between the design of the synchrotron, which intentionallyintroduces structural symmetries, and the orbit feedback system, which uses the SVD andtherefore ignores the structural symmetries. We show that the orbit response matrix inher-its the circulant symmetry and/or the centrosymmetry. We present the block-circulant, thecentrosymmetric and two different combined decompositions that are possible when bothstructural symmetries are present. We illustrate the decompositions using the Diamond-II orbitresponse matrix. We also address the unavoidable case of broken symmetry. For each of theblock-circulant, centrosymmetric and the combined-symmetry cases, we derive formulae forapproximating the orbit response matrix using a matrix that has the symmetric properties andminimizes the Frobenius norm error. We show how the asymmetry of the orbit response matrixcan be concentrated in certain elements of the symmetrical decomposition. We conclude ouranalysis by demonstrating the main advantage of exploiting structural symmetries: increasedcomputation-speed. Using a C-language implementation on the device used for the real-time orbit feedback of the Diamond-II and the APS synchrotrons, we compare the SVDapproach with the different symmetrical decompositions and demonstrate the improvement incomputation speed. We show that accelerating the controller computations reduces the time-delay and allows for faster sampling rates or for the deployment of more advanced controlalgorithms, such as algorithms that use real-time optimization [12].This paper is structured as follows. In Section II, block-circulant and centrosymmetric matricesare briefly outlined. More details are included in Appendix A. Section III presents our resultson the symmetric structure of the orbit response matrix and its decompositions and Section IVaddresses the case of broken symmetry. In Section V, our results are summarized in a casestudy of the Diamond-II synchrotron, in which the controller is simulated using symmetricapproximations, the nominal stability is verified and the speed advantages of a controller thatexploits the symmetric structures are demonstrated.II. P RELIMINARY T ECHNICAL M ATERIAL
A. Notation
The set of real, strictly positive real and complex numbers are denoted by R , R + and C ,respectively, with √ − = i . For matrices A and B , let A ⊗ B denote the Kronecker product,diag ( A, B ) the block diagonal concatenation and A ○ B the Hadamard (elementwise) productof two matrices. Let I n represent the identity matrix in R n × n . For a scalar, vector or matrix a ,let a ∗ denote its Hermitian transpose and ¯ a its element-wise complex conjugate. Let Re { a } and Im { a } denote its real and imaginary part, respectively. For A ∈ R m × n , let ∥ A ∥ F denoteits Frobenius norm that is defined as ∥ A ∥ F = √ ∑ mi = ∑ nj = a ij , where a ij denotes the (scalar)entry of A at row i and column j . Let sgn represent the signum function and mod the modulooperator. B. Matrices with Symmetric Structures
In an orbit feedback system, the calculation of optimal set-points for the corrector magnetsrequires at each sampling instant a matrix-vector multiplication. It will be shown that thematrix inherits certain symmetry properties from the storage ring and that these propertiescan be used to simplify the computations needed for the orbit feedback system.
Definition 1.
Let BC ( n, p, m ) ⊂ R np × nm denote the set of block-circulant matrices of order n that have the form B = ⎡⎢⎢⎢⎢⎢⎢⎢⎣ b b . . . b n − b n − b . . . b n − ⋮ ⋮ ⋱ ⋮ b b . . . b ⎤⎥⎥⎥⎥⎥⎥⎥⎦ , b i ∈ R p × m . (5)Consider the Fourier matrix F n ∈ C n × n , defined as F n = √ n ⎡⎢⎢⎢⎢⎢⎢⎢⎣ . . . w . . . w n − ⋮ ⋮ ⋮ w n − . . . w ( n − )( n − ) ⎤⎥⎥⎥⎥⎥⎥⎥⎦ , (6)with w = e i πn and F ∗ F = I n . Every block-circulant ( BC ) matrix is block-diagonalized by theFourier matrix [15], ˆ B = ( F ∗ n ⊗ I p ) B ( F n ⊗ I m ) = diag ( ν , . . . , ν n − ) , (7)with ν i ∈ C p × m and ( F ∗ n ⊗ I p ) ( F n ⊗ I p ) = I np . Equivalently, the block ν j can also be obtainedfrom ν j = n − ∑ k = b k e − i πjkn . (8)The product F n x yields the coefficients of the discrete Fourier transformation of the vector x . Because the Fourier matrix appears in (7), the computation speed of a matrix-vectormultiplication Bx can be increased significantly by transforming it to the Fourier domain, i.e.by computing Bx = ( F n ⊗ I p ) ˆ B ( F ∗ n ⊗ I m ) x . The computational efficiency arises from thepossibility to employ m parallel Fast Fourier Transformations for computing products like ( F ∗ n ⊗ I m ) x and the fact that ˆ B is block-diagonal. For the case that all elements of B arenon-zero, the computation time is reduced by a factor ( npm + ( p + m ) n log n ) / ( n pm ) . (9) Definition 2.
Let CS ( q, t ) ⊂ R q × t denote the set of centrosymmetric matrices of the form R = ⎡⎢⎢⎢⎢⎣ r r J q r J t J q r J t ⎤⎥⎥⎥⎥⎦ , r i ∈ R p × t , (10) where J k = ⎡⎢⎢⎢⎢⎣ ⋰ ⎤⎥⎥⎥⎥⎦ ∈ R n × n and with RJ t = J q R . This formula serves as a rough estimate of the reduction in computation time and does not consider the detailsof the implementation, such as the complex arithmetic and the structure of the block-diagonalized matrix.
A centrosymmetric ( CS ) matrix is block-diagonalized by [16] ˜ R = T T q RT t = diag ( r − r J t , r + r J t ) , (11)where the centrosymmetric transformation is defined as T k = √ [ I k I k − J k J k ] ∈ R k × k , (12)with T T k T k = I k . As for BC matrices, the computation speed of a matrix-vector multiplication Rx can be increased significantly by transforming it to the centrosymmetric domain and onecan show that (9) holds for n = , p = q and m = t (see Appendix A.2). Definition 3.
Let
SCS ( q, t ) ⊂ R q × t denote the set of skew-centrosymmetric matrices of theform D = ⎡⎢⎢⎢⎢⎣ d d − J q d J t − J q d J t ⎤⎥⎥⎥⎥⎦ , d i ∈ R p × t , (13) with DJ t = − J q D . In contrast to CS matrices, a skew-centrosymmetric ( SCS ) matrix is block anti-diagonalizableby (11), i.e. ˜ D = T T q DT t = [ d − d J t d + d J t ] . (14)III. M AIN R ESULTS
In the following, it will be assumed that the storage ring is divided into S sections of equallength L / S . For demonstrating the BC and the CS properties of the orbit response matrix, itwill be assumed that the β -function is periodic with period L / S and centrosymmetric withrespect to L / , respectively, i.e. β can be mirror-reflected about the middle of the storagering. For our convenience, it will be assumed that S , the number of corrector magnets persection and the number of monitors per section are all even. A. Properties of the Orbit Response Matrix
Proposition 1 (Block-Circulant R ) . Suppose that β ( s ) = β ( s + L / S ) and that each one ofthe S storage ring sections contains N B beam position monitors and N C corrector mag-nets, placed at ring locations s B , . . . , s BN B and s C , . . . , s CN C , respectively. Suppose that thisarrangement is repeated for the ( S − ) following sections. Then R ∈ BC ( S, N B , N C ) .Proof: Partition the first N B rows into N B × N C blocks. Let σ ( ⋅ ) denote the module operatorthat is formulated as σ ( n + kN C ) = ( n + kN C − SN C ) + . According to the BC structure (5), we must show that R m + kN B ,σ ( n + kN C ) = R m,n for k = , . . . , S − , m = , . . . , N B and n = , . . . , N C . From the definition of R m,n in (2): R m + kN B ,σ ( n + kN C ) = √ β Bm + kL / S β Cn + kL / S ( πQ ) cos ( πQ − ∣ φ Bm + kL / S − φ Cn + kL ∣) = √ β Bm β Cn ( πQ ) cos ( πQ − ∣ φ Bm + k πQS / L − ( φ Cn + k πQS / L )∣) = √ β Bm β Cn ( πQ ) cos ( πQ − ∣ φ Bm − φ Cn ∣) = R m,n , where we used the fact that φ ( s + kL ) = φ ( s ) + k πQL / S for a periodic β , which can be verifiedfrom (3). Proposition 2 (Centrosymmetric R ) . Suppose that β ( L / + s ) = β ( L / − s ) . In addition,suppose that the position of the monitors and magnets is reflection-symmetric as well, i.e. foreach s Xk there is a s Xp s.t. s Xp = L − s Xk for X = { B, C } . Then R ∈ CS ( SN B / , SN C / ) .Proof: The top-half of the matrix must be a vertically and horizontally reflected versionof the bottom-half of the matrix. For the top-left and bottom-right sub-blocks of the ma-trix, R SN B / − n,SN C / − m must equal R SN B / + n + ,SN C / + m + for all combinations of n = , . . . , ± SN B / − and m = , . . . , SN C / − . After setting ˜ s Ak = L / − s AX + k and notingthat φ ( L / ± s ) = φ ( L / ) ± φ ( s ) , we obtain: R SN B / − n,SN C / − m = √ β ( L / − ˜ s Bn ) β ( L / − ˜ s Cm ) ( πQ ) × cos ( πQ − ∣ φ ( L / − ˜ s Bn ) − φ ( L / − ˜ s Cm )∣) , = √ β ( L / + ˜ s Bn ) β ( L / + ˜ s Cm ) ( πQ ) cos ( πQ − ∣ φ ( ˜ s Cm ) − φ ( ˜ s Bn )∣)´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ = cos ( πQ −∣ φ ( L / + ˜ s Cm )− φ ( L / + ˜ s Bn )∣) , = R SN B / + n + ,SN C / + m + , and analogously for the top-right and bottom-left sub-blocks.Propositions 1 and 2 are intuitive results. If the magnetic lattices, beam position monitorsand corrector magnets are arranged in a symmetric pattern, then the orbit response matrixinherits the same symmetric pattern. The BC property means that a circulant shift of N B and N C elements can be applied to the beam displacement y and magnet inputs u in (1) withoutchanging the system behavior, while the CS property means that each vector can be mirror-reflected about its middle. If Propositions 1 and 2 simultaneously hold, the orbit responsematrix inherits additional properties, which are summarized in the following proposition.Note that if β is periodic with period L / S and CS with respect to L / , then β is CS withrespect to the middle of each of the S sections. Proposition 3 (Centrosymmetric and Block-Circulant R ) . Suppose that the conditions inPropositions 1 and 2 all hold, so that R ∈ BC ( S, N B , N C ) ∩ CS ( SN B / , SN C / ) , and let r k ∈ R N B × N C , k = , . . . , S , denote the sub-blocks of R . Then: r , r S / ∈ CS ( N B / , N C / ) r S / + k J N C = J N B r S / − k Proof:
See Appendix A.3.
B. Decompositions of the Orbit Response Matrix
The BC property ensures that R can be block-diagonalized by pre- and post-multiplication withthe discrete Fourier matrix F S defined in (6). By defining ˆy ∶ = ( F ∗ S ⊗ I N B ) y , ˆu ∶ = ( F ∗ S ⊗ I N C ) u and ˆd ∶ = ( F ∗ S ⊗ I N B ) d , the dynamics (1) can be mapped into the discrete (spatial) Fourierdomain as ˆy ( z − ) = ˆR g ( z − ) ˆu ( z − ) + ˆd ( z − ) , (15)where ˆR ∶ = diag { ˆR , . . . , ˆR S − } . When a vector is mapped into the Fourier domain as in ˆy = ( F ∗ S ⊗ I N B ) y , the Kronecker product between F ∗ S and I N B means that the k th displacementsof each superperiod are grouped. The Fourier transform is then applied to equidistant samplesat s Bm + kL , k = , . . . , S − . This yields the Fourier coefficients for the spatial frequencies ω k = πk / S . The block-diagonal structure of ˆR means that the spatial Fourier coefficientsof the displacements at frequency ω k are not modified by magnetic inputs at frequency ω j for n ≠ j . The Fourier coefficients are, however, influenced by other Fourier coefficients ofthe same spatial frequency that have a different starting point s Cn for the equidistant samples s Cn + kL . Note that one could apply the spatial Fourier transform to any matrix, but the resulting ˆR is block-diagonal if and only if R is BC [15].Analogously to the BC case, the dynamics (1) can be mapped to the CS domain by defining ˜y ∶ = T T SN B / y , ˜u ∶ = T T SN C / u and ˜d ∶ = T T SN C / d . The resulting ˜R ∶ = T T SN B / R T SN C / isblock-diagonal if and only if R is CS . The transformation ˜y = T T SN B / y groups elements k and k + SN C / of y and assigns their sum and differences to ˜y . The block-diagonalized ˜R reflects the fact that the sum (difference) of the displacements, is solely modified by the sum(difference) of the magnetic kicks.When R is both BC and CS , the matrices ˆR as well as ˜R can be further decomposed. For thedecomposition of ˆR , we start by rewriting the complex-valued blocks of the BC decomposition ˆR using (8) as ˆr n = r + ( − ) n r S / + S / − ∑ k = ( r k e − i πnkS + J N B r k J N C e i πnkS ) , where the second part of Proposition 3 was used after reformulating it as r S − k = J N B r k J N C .Separating the real and imaginary parts of ˆr n yieldsRe { ˆr n } = r + ( − ) n r S / + S / − ∑ k = cos ( πnkS ) ( r k + J N B r k J N C ) , Im { ˆr n } = S / − ∑ k = sin ( πnkS ) ( J N B r k J N C − r k ) . Common to matrices with symmetric structures, such as BC or CS matrices, is that they forman algebra (see Appendix A). Because r , r S / and r k + J N B r k J N C are CS , the real part of This can be shown by pre- and post-multiplication with J N B and J N C , respectively. ˆR n is CS , while the imaginary part is SCS , because J N B r k J N C − r k is SCS . Each of theFourier blocks ˆR n can therefore be pre- and post-multiplied by T T N B / and T N C / , which willseparate the real and imaginary part, because CS matrices are block-diagonalized while SCS matrices are block anti-diagonalized by the transformation (12).The decomposition of the CS decomposition ˜R can be found in Appendix A.4, in whichit is shown that, if R is BC as well, each of the blocks of ˜R is CS and can thereforebe decomposed using (11). Note that the BC structure is a more stringent requirement thanneeded, i.e. the doubly CS decomposition only requires that the β -function is CS with respectto L / . Table I summarizes the results from this section and characterizes the BC , CS and theirfurther decompositions by showing the formulae for block-diagonalization and the resultingshapes of the block-diagonalized matrices. The table also addresses symmetric approximationsof R , which are treated in the following section.IV. B ROKEN S YMMETRY
In practice, the regular arrangement of magnetic lattices, beam position monitors and correctormagnets is compromised by space constraints, e.g. there will be one section where the injectiondevice – the entry point for the electrons – will need to be fitted. This will lead to anasymmetric placement of one or more of the aforementioned components. It can also bethat some of the magnets in the magnetic lattices are not perfectly aligned to the vertical orhorizontal plane. This, in turn, leads to an asymmetry of the β -function. Because the symmetryof the orbit response matrix is solely based on the symmetry of the betatron function as wellas on the placement of monitors and magnets, the orbit response matrix will inherit anyasymmetry. If we are to exploit the symmetric structure of the synchrotron for the controller,then the matrix K in (4) must have the symmetry properties. If R is BC and/or CS , the gainmatrix K , which can also be computed from K = ( R T R + µ I SN B ) − R T , will necessarily havethe symmetry properties because each of our symmetric structures form an algebra. A. Symmetric Approximations
Given a perturbed orbit response matrix R p that does not satisfy the symmetry conditions,a matrix R ⋆ is sought that approximates R p and features the BC and/or CS properties. Thisproblem can be formulated as an optimization problem, R ⋆ = arg min X ∈ S ∥ X − R p ∥ F , (16)where S ∈ { BC , CS , BC ∩ CS } and the Frobenius norm was used because this norm leads toclosed-form solutions. If temporal dynamics were involved, a more appropriate choice wouldbe the ∞ -norm, for which R ⋆ would also reflect the stability properties of R p . In [17], asolution is derived for the case that R p is a Toeplitz matrix and in Appendix B, the proofs areextended for a general R p and S ∈ { BC , CS , BC ∩ CS } . The results obtained are summarizedin Table I. They essentially consist of averaging over the sub-blocks of R p according tothe corresponding structure of S , e.g. when S = BC the diagonal block r ⋆ is obtained fromaveraging over all sub-blocks of R p that are lying on the diagonal. The resulting problem would be a linear program.
TABLE I: Symmetric Decompositions for S = Decomposition Centrosymmetric ( CS ) Decomposition of the CS decomposition ( CS − BC )Diagonalization( ˆR ) T T SN B / R T SN C / ( I ⊗ T T SN B / ) T T SN B / R T SN C / ( I ⊗ T SN C / ) Approximation † ( R ⋆ ) ( R p + J SN B R p J SN C ) S ∑ S − k = ( Ω kS ⊗ I N B ) T ( R p + J SN B R p J SN C ) ( Ω kS ⊗ I N C ) Shape of ˆR Shape of ˆ∆ Decomposition Block-Circulant ( BC ) Decomposition of the BC decomposition ( BC − CS )Diagonalization( ˆR ) ( F ∗ S ⊗ I N B ) R ( F S ⊗ I N C ) ( I S ⊗ T T N B / ) ( F ∗ S ⊗ I N B ) R ( F S ⊗ I N C ) ( I S ⊗ T N C / ) Approximation † ( R ⋆ ) S ∑ S − k = ( Ω kS ⊗ I N B ) T R p ( Ω kS ⊗ I N C ) S ∑ S − k = ( Ω kS ⊗ I N B ) T ( R p + J SN B R p J SN C ) ( Ω kS ⊗ I N C ) Shape of ˆR Shape of ˆ∆ Blue, red and gray blocks refer to real, purely imaginary and complex-valued numbers, respectively. † The matrix Ω S is defined in (21). B. Approximation Error
When the gain matrix K is computed using an approximation R ⋆ , the stability propertiesof the resulting closed-loop system might be affected, i.e. the system might be stable if K is computed using R p but unstable when computed using R ⋆ . To measure the amount ofasymmetry, an approximation error is defined as ∆ ∶ = R p − R ⋆ . For a symmetric approxi-mation problem, such as (16), the structure of ∆ can be determined from transforming theoptimization (16) into the symmetric domain S , i.e. by rewriting the norm in (16) using thesymmetric transformations T S ,l , T S ,r as ∥ T ∗ S ,l ( X − R p ) T S ,r ∥ F = ∥ ˜ X − ˜R p − ˜R p – ∥ F , (17)where ˜ X = T ∗ S ,l X T S ,r , ˜R p + ˜R p – = T ∗ S ,r R p T S ,r and where ˜R p has the same structure as ˜ X such that Re ( ˜R p – ) ○ Re ( ˜R p ) = , Im ( ˜R p – ) ○ Im ( ˜R p ) = . (18)Note that the norm is invariant w.r.t. multiplication with an orthonormal matrix [18, Ch. 2.3.5,p. 75]. From (17), it becomes clear that ˜R ⋆ = ˜R p and ˜∆ = ˜R p – and the solution to (16) couldbe found by setting R ⋆ = T S ,l ˜R ⋆ T ∗ S ,r . The shapes of ˜∆ for S ∈ { BC , CS , BC ∩ CS } aredepicted in Table I. Note that for the doubly CS decomposition (column CS − BC in Table I)we use the approximation that yields R ⋆ ∈ BC ∩ CS and (18) therefore does not hold.V. C
ASE S TUDY : D
IAMOND -IIDiamond Light Source is the UKs national synchrotron facility, which produces synchrotronlight for research. The Diamond-I synchrotron is a 3rd-generation light source in whichelectrons circulate around the m storage ring at an energy of GeV. The upcomingDiamond-II conversion will upgrade the synchrotron to a 4th-generation light source thatoperates at . GeV and introduce various changes, such as new lattice technologies and anew orbit feedback system [19]. At Diamond-I, the orbit feedback system uses monitorsto operate × corrector magnets – one set for the vertical and one for the horizontalplane – at a frequency of kHz. The feedback reduces the trajectory error of the electronsto 5 µ m in the horizontal and 600 nm in the vertical plane [20]. Diamond-II will use monitors and × corrector magnets that will be operated at a frequency of kHz. -10-50510 Fig. 1: Diamond-II orbit response matrix with S = sections. − K c ( z ) R p g ( z ) + y ( z ) d ( z ) u ( z ) Fig. 2: Feedback system.
A. Symmetric Approximations
The Diamond-II storage ring will be arranged in S = superperiods and for this case studywe are using a preliminary version of the orbit response matrix, which is shown in Fig. 1for the vertical plane. The BC , CS and BC ∩ CS approximations were applied to the orbitresponse matrix of Fig. 1 and the 2-norm ( ∥ X ∥ ), the average absolute ( avg ij (∣ X ij ∣) ) andthe maximum absolute magnitude ( max ij (∣ X ij ∣) ), respectively, of the resulting approximationerrors are compared to R p in Table II. The small errors for the BC approximation show thatthe BC symmetry property is an accurate assumption for R p , whereas the CS and BC ∩ CS approximations yield significantly larger errors. Their maximum singular values, however, aresubstantially smaller than the maximum singular value of R p .TABLE II: Approximation Error X ∥ X ∥ avg ij (∣ X ij ∣) max ij (∣ X ij ∣) R p
644 2.7112 11.8176 ∆ - BC ∆ - C S
50 0.1378 1.4834 ∆ - BC ∩ CS
50 0.1378 1.4835
B. Orbit Feedback Controller
As a proof of concept for Diamond-II, we are interested in how the standard controller u ( z − ) = K c ( z − ) y ( z − ) performs when K is obtained using a symmetric approximation,i.e. K = ( R ⋆ T R ⋆ + µ I SN B ) − R ⋆ T , while the process model is given by the asymmetricDiamond-II orbit response matrix R p . We will focus on the orbit correction for the verticalplane, but the results are comparable for the horizontal plane. Fig. 2 shows the control systemin its standard configuration. It is assumed that all actuators have the same dynamics g ( z − ) , g ( z − ) = z − d b + b z − − az − , where d = is the delay in terms of time steps and the parameters b , b and a can be foundin [6, p. 207]. An internal model controller is used to form the dynamic part of the controller c ( z − ) [6, eq. (20)].The simulation requires the disturbances d ( t ) as an input. Because no such measurements areavailable for Diamond-II yet, the measurements from Diamond-I are used. The disturbancevector is augmented to fit the dimensions of Diamond-II. First, the = − monitor out-puts are copied and appended to the measurements. Second, the augmented disturbances Frequency (Hz) I B M ( m ) Fig. 3: Simulation of the controller. Also shows the augmented Diamond-I measurements forenabled and disabled feedback.
Frequency (Hz) I B M ( m ) Fig. 4: Close-up of Fig. 3 with the BC and BC ∩ CS approximations.are transformed into mode-space using an SVD of the Diamond-II orbit response matrix andthe power spectrum of the modes plotted, such as in [20, Ch. 3.5, pp. 68–72]. Third, thedisturbance spectrum is scaled to obtain a power spectrum comparable to [20, Fig. 3.11],where the modes associated to large-magnitude singular values show a larger amplitude. Theresulting disturbance profile is depicted in Fig. 3 (labeled by Measured, off).The performance of the controller is measured using the integrated beam motion (IBM), whichis defined as the square root of ∑ Ff = F ∣ y ( f )∣ , where y ( f ) is the discrete Fourier transformof a monitor output and F the frequency in Hz. Fig. 3 shows the average IBM across all beamposition monitors of the storage ring for K computed using R p and the CS approximation of R p . For clarity, the simulation results for the BC and BC ∩ CS approximations are omitted inFig. 3, but are shown in the close-up in Fig. 4. The results show that the controller performsonly slightly worse when a symmetric approximation is used. In Fig. 4, it can be seen thatthe CS approximation yields a slightly larger average trajectory error, which is related to thelarge approximation error.The nominal stability of the controller, i.e. when no uncertainty is present, can be verified bycalculating the poles of the closed-loop transfer functions. The closed-loop transfer functions of the system in Fig. 2 are given by y ( z − ) = ( I N B + R p K L ( z − )) − d ( z − ) , (19a) u ( z − ) = K c ( z − ) ( I N C + R p K L ( z − )) − d ( z − ) , (19b)where L ( z − ) = g ( z − ) c ( z − ) . When the controller uses a symmetric approximation, thecontrol system cannot be diagonalized by an SVD of R p , such as in [6]. One can, however,use an eigendecomposition of R p K to diagonalize (19a) and (19b). E.g. for (19b), one cansubstitute R p K = V DV − , where D is the diagonal matrix of eigenvalues and the columnsof V are the corresponding eigenvectors, which yields u = K V ( I N C + DL ( z − )) − c ( z − ) V − d ( z − ) , and the stability can be verified by computing the poles of the diagonal matrix ( I N C + DL ( z − )) − c ( z − ) . The feedback system is stable if all poles of (19a) and (19b) have magnitude smaller than 1,which is satisfied for K computed using R p as well as using the symmetric approximations. C. Benchmarks on Hardware
The Diamond-I feedback system is implemented on processors, which are distributedaround the storage ring and require a complex distributed network topology, while for Diamond-II the computations will be centralized and the feedback system will be implemented on onemulticore processor [21]. The new setup simplifies the network topology but increases theperformance requirements on the multicore processor. The clock-frequency of the processor is . GHz and the targeted operating frequency of kHz therefore allows for , proces-sor cycles. Without decomposition, the controller computations require × ≈ , multiply-accumulate operations. For demonstrating the speed advantages of the decomposition,the matrix-vector multiplication required by the controller has been implemented on theprocessor. When a decomposition is used, the control input is computed as u ( z − ) = T D ,l ˜K T T D ,r c ( z − ) y ( z − ) , (20)where D ∈ { CS , BC , CS BC , BC CS } refers to the decompositions in Table I and T S ,l , T S ,r to the corresponding transformation, and ˜K is the decomposed gain matrix.Fig. 5 shows the results that were obtained for the implementation of (20) on a single core ofthe processor. The performance is measured as / t , where t is the time required to executeone matrix-vector multiplication and the horizontal gray bars refer to the theoretical speed-upthat was calculated using (9) and the computation frequency of kHz, which correspondsto the frequency of the matrix-vector multiplication without decomposition. Using the BC CS decomposition, the computation frequency is more than eleven times faster than withoutdecomposition and beyond the targeted operating frequency of kHz. The results also showthat for the BC CS decomposition the speed-up of the implementation is significantly largerthan the theoretical prediction. The reason is that, in addition to the number of operationsrequired, the performance of the processor is also limited by memory operations, e.g. thetime needed to transport the matrix data from the memory to the core. The reduced memoryrequirements of the decomposition therefore indirectly benefit the computation time.
11 23 41 59 123
Decomposition k H z Fig. 5: Performance measurement of the controller. The leftmost timing does not use anydecomposition. The gray bars refer to the theoretical speed-up.VI. C
ONCLUSIONS
In this paper, we have shown that the orbit response matrix of a synchrotron inherits themirror-reflective and periodic properties of the betratron function. The mirror-reflective andperiodic properties of the betatron function manifest themselves in a centrosymmetric andblock-circulant structure of the orbit response matrix, which results in a controller gain matrixthat has the same structure. These structural symmetries can be used to decompose the gainmatrix and perform the matrix-vector multiplication, which is required for computing thecontrol inputs, in the symmetric domain. For the symmetries discussed in this paper, thetransformation of a vector into the symmetric domain is computation-efficient and the matrix-vector multiplication in the symmetric domain requires far less multiply-accumulate operationsthan in the original domain.In practice, the mirror-reflective and periodic properties of the betatron function are onlyapproximate, e.g. a symmetry might be broken at some point around the storage ring due oddplacements of monitors or magnets. The asymmetry of the betatron function is inherited bythe orbit response matrix, which, in turn, results in an asymmetric controller gain matrix thatcannot be block-diagonalized using a symmetric decomposition. To recover the symmetricstructure, an optimization problem was formulated in which a matrix was sought that has thecorresponding symmetry properties and approximates the asymmetric orbit response matrix.We presented closed-form solutions for the centrosymmetric, block-circulant and combinedsymmetry cases and an alternative approach that finds the solution by transforming theproblem into the symmetric domain. This alternative approach showed that the approximationerror-matrix always has a particular structure, which complements the approximation in thesymmetric domain.We concluded the paper with a case study of the Diamond-II synchrotron. We approximatedthe Diamond-II orbit response matrix for the centrosymmetric, block-circulant and combinedsymmetry cases and the approximation errors showed that the block-circulant approximationyields a small error, while the centrosymmetric and combined symmetry cases yielded largererrors. The controller gain matrix was computed using the different approximations andthe closed-loop simulations showed that there is only a minuscule difference in trajectoryerror correction between the gain matrices computed using the approximations and using theasymmetric orbit response matrix. The nominal stability for all approximations was verified using an eigendecomposition of the closed-loop transfer functions. The robust stability analysisusing the structured approximation error is currently being considered.We completed the case study using a single-core implementation of the controller computationsand demonstrated the significant improvement of the computation-speed when the matrix-vector multiplication is carried out in the symmetric domain. For the combined symmetrycase, the controller computations in the symmetric domain were more than ten times fasterthan in the original domain and already beyond the targeted operating frequency of kHz.It is expected that, after parallelizing the matrix-vector multiplication, the computations willbe fast enough to enable the use of more advanced control algorithms, such as algorithms thatrequire real-time optimization and explicitly consider system constraints, such as amplitudeand slew-rate constraints of the magnets.R EFERENCES [1] H. Wiedemann,
Particle Accelerator Physics , 4th ed. Springer Berlin Heidelberg, 2007.[2] L. Yu, E. Bozoki, J. Galayda, S. Krinsky, and G. Vignola, “Real-time harmonic closed orbit correction,”
Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectorsand Associated Equipment , vol. 284, no. 2, pp. 268–285, December 1989.[3] M. Martini, “An introduction to transverse beam dynamics in accelerators,” CERN, Geneva, Switzerland,Tech. Rep., March 1996.[4] F. R. Elder, A. M. Gurewitsch, R. V. Langmuir, and H. C. Pollock, “Radiation from electrons in a synchrotron,”
Physical Review , vol. 71, pp. 829–830, June 1947.[5] E. Plouviez and F. Uberto, “The orbit correction scheme of the new EBS of the ESRF,” in , Hamburg,Germany, May 2011, pp. 694–697.[6] S. Gayadeen and S. R. Duncan, “Design of an electron beam stabilisation controller for a synchrotron,”
Control Engineering Practice , vol. 26, pp. 201–210, May 2014.[7] P. Leban, E. Janezic, and M. Sjstrm, “Fast orbit feedback application at MAX IV and SOLARIS storagerings,” in , no. 5, Dresden, Germany, July2014, pp. 1748–1750.[8] J. Carwardine, G. Decker, K. Evans, A. Hillman, F. Lenkszus, R. Merl, and A. Pietryla, “Commissioningof the APS real-time orbit feedback system,” in ,vol. 2, Vancouver, British Columbia, Canada, May 1997, pp. 2281–2283.[9] S. Gayadeen, S. R. Duncan, and G. Rehm, “Optimal control of perturbed static systems for synchrotronelectron beam stabilisation,”
IFAC PapersOnLine , vol. 50, no. 1, pp. 9967 – 9972, 2017, 20th IFAC WorldCongress.[10] S. Gayadeen and S. R. Duncan, “Uncertainty modeling and robust stability analysis of a synchrotron electronbeam stabilisation control system,” in , Maui,Hawaii, December 2012, pp. 931–936.[11] R. D’Andrea and G. E. Dullerud, “Distributed control design for spatially interconnected systems,”
IEEETransactions on Automatic Control , vol. 48, no. 9, pp. 1478–1495, September 2003.[12] I. Kempf, P. J. Goulart, and S. R. Duncan, “Alternating direction of multipliers method for block circulantmodel predictive control,” in , Nice, France,December 2019.[13] P. Massioni and M. Verhaegen, “Subspace identification of circulant systems,”
Automatica , vol. 44, no. 11,pp. 2825–2833, November 2008.[14] S. H. Mirza, R. Singh, P. Forck, and H. Klingbeil, “Closed orbit correction at synchrotrons for symmetricand near-symmetric lattices,”
Physical Review Accelerators and Beams , July 2019. [15] P. J. Davis, Circulant Matrices . Wiley-Interscience, 1979.[16] A. Cantoni and P. Butler, “Eigenvalues and eigenvectors of symmetric centrosymmetric matrices,”
LinearAlgebra and its Applications , vol. 13, no. 3, pp. 275 – 288, January 1976.[17] T. F. Chan, “An optimal circulant preconditioner for Toeplitz systems,”
SIAM Journal on Scientific andStatistical Computing , vol. 9, no. 4, pp. 766–771, July 1988.[18] G. H. Golub and C. F. Van Loan,
Matrix Computations , 4th ed. JHU Press, 2013.[19] Diamond Light Source Ltd., “Diamond-II: Conceptual design report,” Diamond Light Source Ltd., Tech.Rep., May 2019.[20] S. Gayadeen, “Synchrotron electron beam control,” DPhil Thesis, University of Oxford, UK, 2014.[21]
Multicore Fixed and Floating-Point Digital Signal Processor (TMS320C6678) , Texas Instruments, 2010,sPR5691E.[22] J. R. Weaver, “Centrosymmetric (cross-symmetric) matrices, their basic properties, eigenvalues, and eigen-vectors,”
The American Mathematical Monthly , vol. 92, no. 10, pp. 711–717, December 1985.[23] G. Li and Z. Feng, “Mirror-symmetric matrices and their application,”
Tsinghua Science and Technology ,vol. 7, no. 6, pp. 602–607, December 2002. A PPENDIX A A.1 Block-Circulant Matrices
For BC matrices it holds that [15] B ( Ω n ⊗ I m ) = ( Ω n ⊗ I p ) B, Ω n = [ I n − ] , (21)where Ω n is the cyclic shift matrix with Ω T n Ω n = I n and Ω nn = I n . Using (21), it can beshown that BC matrices form an algebra [15], i.e. that sums and products of BC matricesyield another BC matrix. Any BC matrix can be represented as [15] B = n − ∑ k = Ω kn ⊗ b k . (22)When the BC matrix B in (7) is real, the blocks ν j possess additional structure that is inheritedfrom the Fourier matrix F n . If n is even, ν and ν n / are real while for i = , . . . , n / it holds that ν i = ¯ ν n − i . If n is odd, the only real-valued block is ν and the latter holdsfor i = , . . . , ( n )/ . The same pattern of complex conjugates is exploited during a FastFourier Transformation and, according to the properties of ν j , only the first n / blocks mustbe considered for a matrix-vector multiplication. A.2 Centrosymmetric and Skew-Centrosymmetric Matrices
As for BC matrices, it can be shown that CS and SCS matrices also form algebras [22]. Byreversing the order of the second half of the rows and columns, a CS matrix can be permutedinto a BC matrix of order 2: Lemma 1.
The permutation matrices P l = diag ( I q , J q ) and P r = diag ( I t , J t ) permute R ∈ CS ( q, t ) into a BC matrix of order n = : P l RP r = ⎡⎢⎢⎢⎢⎣ r r J t r J t r ⎤⎥⎥⎥⎥⎦ ∈ BC ( , q, t ) . Proof:
Evaluating the product P l RP r yields the result.Lemma (1) shows that CS matrices and BC matrices are closely related. For our purpose, itis sufficient to use Lemma (1) to show that (9) holds for n = . A.3 Proof of Proposition 3
To find the algebraic conditions a simultaneously BC and CS matrix satisfies, consider thepermutation matrix J n acting on the cyclic shift matrix Ω n J n Ω n J n = [ I n − ] = Ω T n = Ω n − n , (23)where for the rightmost equality we used the fact that a cyclic downwards-shift ( Ω T n x ) of avector of length n equals a cyclic upwards-shift by n − places ( Ω n − n x ). Using (23), weobtain J Ω kn J = Ω n − n J Ω k − n J = ⋅ ⋅ ⋅ = Ω kn − kn = ( Ω T n ) k = Ω n − kn , where we used Ω − n = Ω T n and Ω nn = I n . Consider a BC matrix X ∈ BC ( n, l, m ) represented asin (22). If X is also to be CS , then it must commute with the reflection matrix, i.e. n − ∑ k = Ω kn ⊗ x k ! = J nl ( n − ∑ k = Ω kn ⊗ x k ) J nm , = ( J n ⊗ J l ) ( n − ∑ k = Ω kn ⊗ x k ) ( J n ⊗ J m ) , = n − ∑ k = J n Ω kn J n ⊗ J l x k J m , = n − ∑ k = Ω n − kn ⊗ J l x k J m = n − ∑ k = Ω n − kn ⊗ J l x k J m . Note that Ω jn has non-zero entries where Ω kn , k ≠ j , has zero entries and vice-versa. Equatingthe terms with the same power of Ω n yields x n − k = J l x k J m , (24)and in particular x = J l x J m and x n / = J l x n / J m . A.4 Decomposition of the Centrosymmetric Decomposition
When R is BC and CS , the CS decomposition ˜R = T T SN B / R T SN C / can be further decom-posed. Consider the partitioning of R into four equal-sized blocks as in Definition 2 suchthat ˜R = diag ( R − R J SN C / , R + R J SN C / ) . From the BC structure (5), R and R areobtained as R = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣ r r . . . r S r S r . . . r S − ⋮ ⋮ ⋱ ⋮ r S + r S + . . . r ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ , R = ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎣ r S r S + . . . r S r S r S . . . r S − ⋮ ⋮ ⋱ ⋮ r r . . . r S ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎦ . The matrices R , R have the CS blocks r and r S / on their diagonals. In addition, the blocksopposite the diagonals are r k and r k − for R and r S / + k and r S / − k for R , i.e. the opposite blocks satisfy the second part of Proposition 3. This entails that R , R ∈ CS ( SN B / , SN C / ) .Note that if R is CS , then so is R J S / N C . Because the sum of two CS matrices is also CS , the blocks of ˜R are CS and each one of the blocks can be further decomposed by pre-and post-multiplication with T T SN B / and T SN C / , respectively. In case SN B or SN C are notdivisible by , the decomposition is still possible but a different transformation matrix mustbe used (see for example [23]). A PPENDIX B B.1 Block-Circulant Approximation
For approximating a matrix R p ∈ R nl × nm with a matrix X ∈ BC ( n, l, m ) , the optimization (16)is reformulated as min { x k } n k = ∥ n − ∑ k = Ω kn ⊗ x k − R p ∥ F , (25)where x k ∈ R l × m , X was partitioned as in (5) and the BC representation (22) was used.Because Ω kn has non-zero elements where Ω jn , j ≠ k , has zero elements, i.e. ∑ S − k = Ω kS = n,n ,where n,n is a matrix of ones, problem (25) can be rewritten as min { x k } n k = n − ∑ k = ∥ Ω kn ⊗ x k − ( Ω kn ⊗ l,m ) ○ R p ∥ F . (26)By partitioning R p as R p = ⎡⎢⎢⎢⎢⎢⎢⎢⎣ r , r , . . . r ,n − r , r , . . . r ,n − ⋮ ⋮ ⋱ ⋮ r n − , r n − , . . . r n − ,n − ⎤⎥⎥⎥⎥⎥⎥⎥⎦ , (27)where r i,j ∈ R l × m , each summand on the right-hand side of (26) can be rewritten as n − ∑ j = ∥ x k − r j, k + j mod n ∥ F , (28)for k = , . . . , n − . Using (26), (28) and the partitioning (27), the minimization (25) can bereformulated as min { x k } n k = n − ∑ k = n − ∑ j = ∥ x k − r k, k + j mod n ∥ F . (29)The minimum of (29) is attained where its derivative is zero. This can be done element-wise foreach element of x k which – after reconstructing blocks x k – yields x ⋆ k = n ∑ n − j = r k, k + j mod n .Using the cyclic shift matrix, the solution is reconstructed as X ⋆ = ∑ n − k = ( Ω kn ⊗ I l ) T R p ( Ω kn ⊗ I m ) / n . B.2 Centrosymmetric Approximation
For approximating a matrix R p ∈ R q × t with a matrix Y ∈ CS ( q, t ) , the optimization (16) isreformulated as min y ,y ∈ R q × t ∥ ⎡⎢⎢⎢⎢⎣ y J q y J t y J q y J t ⎤⎥⎥⎥⎥⎦ − R p ∥ F . (30) If R p is partitioned as R p = ⎡⎢⎢⎢⎢⎣ r r r r ⎤⎥⎥⎥⎥⎦ with r i ∈ R q × t , the minimization (30) can bereformulated as min y ,y ∈ R q × t (∥ y − r ∥ F + ∥ y − J q r J t ∥ F + ∥ y − r ∥ F + ∥ y − J q r J t ∥ F ) (31)where ∥ U XV ∥ F = ∥ X ∥ F for orthonormal U, V was used [18, Ch. 2.3.5, p. 75]. The min-imum of (31) is attained where its derivative is zero. The minimizers y ⋆ and y ⋆ are ob-tained as y ⋆ = ( r + J q r J t ) / and y ⋆ = ( r + J q r J t ) / and Y ⋆ is reconstructed as Y ⋆ = ( R p + J q R p J t ) / . B.3 Block-Circulant and Centrosymmetric Approximation
For approximating a matrix R p ∈ R nl × nm , where n, l, m > are even, with a matrix Z ∈ BC ( n, l, m ) ∩CS ( nl / , nm / ) , the optimization (16) is reformulated as in (25) and the blocks { z k } n k = n / + are substituted using (24), which yields min { x k } n k = ∥ I n ⊗ z + Ω n / n ⊗ z n / + n / − ∑ k = ( Ω kn ⊗ z k + Ω n − kn ⊗ J l z k J m ) − R p ∥ F , = min { x k } n k = ⎛⎝∥ I n ⊗ z − ( I n ⊗ l,m ) ○ R p ∥ F + ∥ Ω n / n ⊗ z n / − ( Ω n / n ⊗ l,m ) ○ R p ∥ F + n / − ∑ k = (∥ Ω kn ⊗ z k − ( Ω kn ⊗ l,m ) ○ R p ∥ F + ∥ Ω n − kn ⊗ J l z k J m − ( Ω n − kn ⊗ l,m ) ○ R p ∥ F )⎞⎠ . (32)As for the BC approximation in Appendix B-A, the Frobenius norms can be separated fordifferent powers of Ω n . The terms for z and z n / can be rewritten as ∥ Ω kn ⊗ z k − ( Ω kn ⊗ ) ○ R p ∥ F = n − ∑ j = ∥ z k − ρ k,j ∥ F , (33)where k = { , n / } , ρ k,j ∶ = r k, k + j mod n with R p partitioned as in (27). According to (24),sub-blocks z and z n / must be CS . Sub-blocks z k and ρ k,j are therefore partitioned as z k = ⎡⎢⎢⎢⎢⎣ z k J l / z k J m / z k J l / z k J m / ⎤⎥⎥⎥⎥⎦ , ρ k,j = ⎡⎢⎢⎢⎢⎣ ρ k,j ρ k,j ρ k,j ρ k,j ⎤⎥⎥⎥⎥⎦ , and the right-hand side of (33) rewritten as n − ∑ j = (∥ z k − ρ k,j ∥ F + ∥ z k − J l / ρ k,j J m / ∥ F + ∥ z k − ρ k,j ∥ F + ∥ z k − J l / ρ k,j J m / ∥ F ) . (34) Note the similarity between (34) and (31). Setting the derivative of (34) to zero, solving for z k , z k and reconstructing z k yields for k = { , n / } z ⋆ k = n n − ∑ j = ( ρ k,j + J l ρ k,j J m ) . (35)The summands in (32) for k =< !1 , . . . , n / are rewritten as n − ∑ j = ∥ z k − ρ k,j ∥ F + ∥ J l z k J m − ρ n − k,j ∥ F ´¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¸¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¹¶ =∥ z k − J l ρ n − k,j J m ∥ F . Setting the derivative to zero yields for k = , . . . , n / z ⋆ k = n n − ∑ j = ( ρ k,j + J l ρ n − k,j J m ) , which is identical to (35). After reconstruction, the matrix Z ⋆ is obtained as Z ⋆ = n n − ∑ k = ( Ω kn ⊗ I l ) T ( R p + J nl R p J nm ) ( Ω kn ⊗ I m ))