Convex optimisation of gradient and shim coil winding patterns

Journal of Magnetic Resonance 244 (2014) 36–45

Contents lists available at ScienceDirect

Journal of Magnetic Resonance journal homepage: www.elsevier.com/locate/jmr

Convex optimisation of gradient and shim coil winding patterns Michael S. Poole a,⇑, N. Jon Shah a,b a b

Institute of Neuroscience and Medicine - 4, Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Straße, 52425 Jülich, Germany Department of Neurology, Faculty of Medicine, RWTH Aachen University, JARA, Aachen, Germany

a r t i c l e

i n f o

Article history: Received 7 January 2014 Revised 20 March 2014 Available online 4 May 2014 Keywords: Gradient coil Shim coil Convex optimisation Field synthesis Boundary element method Wire spacing Sparsity Infinity norm Manhatten norm ‘p -Norm

a b s t r a c t Gradient and shim coils were designed using boundary element methods with convex optimisation. The convex optimisation framework permits the prototyping of many different cost functions and constraints, for example ‘p -norms of the current density. Several examples of gradients and shims were designed and simulated to demonstrate this, as well as to investigate the behaviour of new cost functions. A mixture of ‘1 - and ‘1 -norms of the current density, when used as a regularisation term in the field synthesis problem, was found to produce coils with bunches of equally spaced windings that do not take up all of the available surface. This is thought to be beneficial in the design of coils that will be manufactured from wire with a fixed cross-section. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Magnetic resonance imaging requires 3 different types of magnetic field produced by the main magnet, gradient and radiofrequency (RF) coils. The main magnet generates an exquisitely homogeneous, stable and intense magnetic field usually with superconducting wires and is designed with a great deal of engineering knowledge and experience [1]. Radio-frequency coils are usually simple rectangular or circular loops combined with capacitors in circuits resonant at the Larmor frequency. Gradient coils must produce their magnetic fields to within a few % in the imaging region and are made from coils of copper wire at room temperature. Gradient coils must operate safely in the audio frequency range up to 20 kHz and support a maximum current of more than 800 A with more than 2000 V potential difference across each terminal. Shim coils are similar to gradient coils but need not support such high current and voltages. They are used to correct the magnetic field of the main magnet. The design of gradient and shim coils [2] can be divided into two types of methods, discrete wire and current density method, each with arguments for and against their use. Discrete wire methods model the thin wires that constitute a coil and continuous ⇑ Corresponding author. Fax: +49 2461 61 1919. E-mail address: [email protected] (M.S. Poole). http://dx.doi.org/10.1016/j.jmr.2014.04.015 1090-7807/Ó 2014 Elsevier Inc. All rights reserved.

current density methods model a coil with a surface current density. In discrete methods, the line-integral form of the Biot–Savart law is used to calculate the magnetic field produced by a coil. The position of the wires is parameterised and a cost function that defines the ‘‘quality’’ of the coil [3,4] is usually non-linear. Therefore, it is common and entirely appropriate to use stochastic optimisation techniques for global optimisation of the cost function with discrete wire coil parameterisation [5]. It is conceptually simple to add various kinds of penalty to the cost function to yield the most desirable design. However, a new parameterisation is required when the coil surface is changed requiring considerable mathematical effort for complicated surface shapes [6]. Continuous current density coil design methods [7,8] provide a good approximation to the manufactured coil if enough wires are used. A thin surface with a current density flowing on it can be modelled as incompressible flow with zero divergence. The surface integral form of the Biot–Savart law is used to calculate the magnetic field produced by a coil. Various parameterisations of continuous surface current density are possible that provide a simple linear relationship with the magnetic induction field. The most versatile of these is a parameterisation in which the current carrying surface is modelled as a connected ensemble of flat polygons, known as a mesh. The current density in the surface is then a vector field that is piecewise uniform [9]. The number of free

37

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

parameters over which to optimise can be greatly reduced with assumptions of symmetry [10,11]. Inversion of the Biot–Savart law to obtain coil designs is ill-posed, but easily solved by weighted Tikhonov regularisation with the power dissipation [7,12] or stored energy [13] of the coil with one-step matrix inversion methods. Since the current density has zero divergence and is confined to a surface, a potential function, the stream-function, is defined over the surface [8]. A scalar function whose curl (about the surface normal) yields the current density, the stream function is piecewise linear over a polygonal mesh [9,14]. The vertices of the mesh possess a stream function value that is linearly interpolated over its immediate neighbourhood of triangles. It is these stream function values that act as free parameters of continuous current density coil design methods. A final step in continuous current density methods is the conversion into discrete wires by taking level sets of the stream-function over the surface. Cost functions involving non-quadratic functions of the freeparameters of the coil design have been studied previously and solved efficiently using custom gradient descent based deterministic optimisation algorithms [4,15–18]. In the present paper we introduce a convenient general formulation of the problem where all components of the optimisation are convex. Several prototyping tools are available for the solution of convex optimisation problems. In the present study we used cvx [19,20] and demonstrate how new cost functions are easily prototyped. In particular, the magnitude of the current density is a convex function of the nodal stream function values and therefore we can minimise a convex norm of the absolute current density. Previously the infinity norm of the gradient magnitude was included in the optimisation to provide coils with minimised maximum current density [15]. These coils exhibit a lower maximum temperature for a given field gradient strength [21]. It should be noted that while it has always been possible to add such terms to the cost function, it is only when formulated in this way that they can be solved efficiently and deterministically and associated in any combination. This paper is arranged as follows: the mathematical formulation of the convex optimisation problem is given, followed by its representation in the parsed cvx code as written in Matlab (The Mathworks, Natick, MA, USA). The ‘1 -, weighted ‘2 - and ‘1 -norms of the current density were used for regularisation using triangular and sinusoidal boundary element methods (BEMs) on arbitrary and cylindrical surfaces, respectively. Fundamental studies (similar to those presented in Refs. [22,21]) of the behaviour of this new optimisation were made by trading the field accuracy for different norms of the current density and also by mixing them in varying amounts with fixed field accuracy. Coil patterns are shown for 0th, 1st and 2nd order solid harmonic target fields with resistance and mixed norm minimised solutions and shielded whole-body X-gradient coils were designed for comparison. Simulations of the resultant novel coil designs were performed and are presented herein with discussions of their relevance to coil construction.

2. Methods 2.1. Problem Formulation Two types of BEM were used in this study. One with flat triangular elements and linear shape functions to model the stream-function of the current density [9,14,23]. The values of the stream-function at the nodes of the mesh, w, define the coil. Meshes were made in Blender (Blender Foundation, Amsterdam, Netherlands) and exported to Matlab. The second method is restricted to a finite-length cylindrical surface on which the stream-function is a weighted sum of truncated sinusoidal functions [10,11].

Matrix equations that transform w to the various coil properties were constructed according to Eqs. (9)–(13) in Ref. [15] and are summarised below.  The z-component of the magnetic field at a series of points,

b ¼ Bw;

ð1Þ

 The Cartesian or cylindrical components of the current density in each triangle,

jx ¼ J x w;

jy ¼ J y w;

jz ¼ J z w;

ja/ ¼ J a/ w;

ð2Þ

 The stored energy in the coil;

W ¼ wT Lc w;

ð3Þ

 The resistive power dissipation of the coil,

P ¼ wT Rc w:

ð4Þ

Previously, a weighted linear combination optimisation problem was used with equality constraints. This is simple to solve when each term is quadratic in w using matrix inversion (with Lagrange multipliers for the equality constraints) [14,23]. A non-quadratic term was introduced to minimise the maximum absolute current density [15]. We generalise this by stating the optimisation in the form of a convex optimisation as defined in Ref. [24]:

minimise f 0 ðwÞ subject to f i ðwÞ 6 bi ;

i ¼ 1; . . . ; m;

ð5Þ

where the functions f0 ; . . . ; fm : Rn ! R are convex. Eqs. (1) and (2) are affine. Eqs. (3) and (4) are convex if Lc and Rc are positive semidefinite, which is the case. Various convex cost functions and constraints were tested in the present study using cvx; a toolbox for Matlab that parses the optimisation problem and passes it in standard form to sdtp3 [25] to find the solution. The simplest optimisation for the coil design problem is the Tikhonov regularised minimisation of the root-mean-squared (RMS) residual field.

minimise kBw  bt k2 þ akwk2 ;

ð6Þ

where bt is the target field and a is the user-defined regularisation parameter. The ‘2 -norm regularisation term, kwk2 , has no practical relevance. Substitution of w with an affine transformation Fw retains convexivity. If F is the matrix obtained from Cholesky decomposition of Lc or Rc then the solution of Eq. (6) will yield, respectively, the inductance (stored energy) or resistance (power dissipation) minimised coil. The minimum resistance coil design problem is written in cvx code as. cvx_begin variable x(length(Rc)); minimize (norm(B  x-bt,2) + alph  quad_form (x,Rc)); cvx_end

It should be noted that it is far simpler and quicker to solve the specific problem described by the code above using direct matrix inversion [9,14,23] than this convex optimisation method. In all cases, optimisations were performed on a Mac Pro (Apple Inc., Cupertino, CA, USA) equipped with 2.8 GHz Intel Xeon processors and 4 GB of RAM, using Matlab R2012a (The Mathworks, Natick, MA, USA), CVX [24] 2.0 (beta) build 945 and sdpt3 [25] version 4.0.

38

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

2.2. Normalisation

track thickness, t, of 3 mm and minimum gap between tracks, g _ , of 1 mm were used.

The matrices, B; J x ; J y ; J z ; J a/ ; Lc and Rc , used in the optimisation are derived from physical properties of the coil and therefore have a large range of scales. These matrices were normalised by their Frobenius norm to ensure that they all have approximately the same scale and therefore improved convergence behaviour. With reference to Eq. (6) for example, B ! B=kBkF , where kBkF is typically on the order of 107 . 2.3. Basic regularisation behaviour It is possible to substitute the power dissipation regularisation term, quad_form (x,Rc), by the maximum absolute current density [15], kjk1 , or the sum of absolute current densities, kjk1 , qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 2 2 where j ¼ jx þ jy þ jz or j ¼ ja/ þ jz . See Appendix A for more

2.4. Mixed regularisation The solutions produced by Eq. (7) represent the very extreme solutions of their regularisation type for a given maximum field error. It was shown previously that, depending on the coil design requirements, it can be beneficial to regularise with a combination of terms [21]. In this section, the same coil geometry and target field (X-gradient) as the previous section were used and Dmax was fixed to 5%. The optimisation problem was

minimise a

kjðwÞkp  K pp kjðwÞkq  K qq þ ð1  aÞ K qp  K pp K pq  K qq

ð8Þ

subject to kBw  bt k1 =kbt k1 6 5%;

detail on the relevance of ‘p -norms of the current density magnitude in coil design. It is also common to replace the ‘2 -norm of the field errors with the maximum field error, which is easily changed in cvx by altering the power of the norm to infinity (see Eq. (A5)). In this section the trade-off between maximum field error and the regularisation term was investigated [22,21]. The optimisation is written:

where a was varied from 0 to 1. If wp is the p-norm regularised solution to Eq. (7), then K pq ¼ kjðwp Þkq . The cases where ½p; q ¼ ½1; 2; ½1; 1 and ½2; 1 were studied. It is also simple to recast Eq. (8) into a constrained optimisation problem by converting one of the terms of the minimisation into an inequality constraint and varying the value of that constraint. Coils were simulated in the same manner as described in the previous section.

minimise kjðwÞkp

2.5. Wire patterns for mixed p ¼ 1 and p ¼ 1 regularised coils

subject to kBw  bt k1 =kbt k1 6 Dmax :

ð7Þ

In Eq. (7), p = 1, 2 or 1, and Dmax was varied over the range 0.5% to 20%. Since all areas of all the triangles are the same, the p ¼ 2 case is equivalent to the resistance (see Eq. (A4)). Eq. (7) for the p ¼ 1 case is written in cvx code as cvx_begin variable x(length(Rc)); minimize (max(norms([Jx  x Jy  x Jz  x], 2, 2))); subject to norm (B  x  bt, inf)./max(bt) <= Dmax; cvx_end

The coil surface was a 1000 mm long cylinder with a diameter of 500 mm that was modelled with 2048 flat triangles (1056 mesh nodes). A spherical target field region with 300 mm diameter was specified by 386 surface points. A total of 404 X-gradient coils were designed; 101 for p ¼ 1, 101 for p ¼ 2 (proportional to resistance), 101 for p ¼ 2 weighted to give minimum inductance (using Eq. (6) T with kwk2 replaced with kL1=2 c wk2 or equivalently w Lc w) and 101 for p ¼ 1. The coils were optimised with normalised target fields (i.e. designed to produce 1 Tm1). The continuous current density and related parameters used in the optimisation can be used to predict coil performance. In order to more accurately predict these parameters, contours of w were obtained to yield the wires of the coil. The number of contour levels, N c , used was maximised whilst ensuring that the minimum wire spacing remained above 5 mm. Biot–Savart integration of the thin wire coil model conducting 1 Amp was used to predict the coil efficiency, g (lT m1 A1). The spacing between the centres of the wires, si , was calculated from the current density magnitude, ji , in the mesh triangle and the minimum wire spacing in the coil, s_ , was found: s_ ¼ ðmaxðwÞ minðwÞÞ= maxðji ÞN c (see Appendix A). The total length of the conductor, c, was simply found by summation of the wire segment lengths of the whole coil. Finally, the inductances and resistances of the coils were simulated using the Fasthenry software at 0 Hz [26]. For the X-gradient coils, a maximum track width of 30 mm,

Here we are interested in the shape of typical wire patterns for different target fields. A basic cylindrical surface was assumed in this section and the stream-function was modelled as a weighted sum of sinusoidal functions in z and /. The surface cylinder diameter was 500 mm and its length was 1400 mm. A cylindrical target field region with two diameters, 150 mm and 300 mm, was defined to be 300 mm long and solid harmonic fields up to 2nd order were targeted, e.g. bt / x2  y2 . Of particular interest here is the design of coils with sparse, but equally spaced windings; i.e. mixed ‘1 - and ‘1 -norms. Such coils are expected to be similar to those produced by existing methods used commonly to design shim coils [27] but with fewer constraints on the wire patterns. The optimisation function used was

Fig. 1. Schematic geometry of the whole-body X-gradient coils. Two cylindrical surfaces are modelled using triangles (grey lines) forming the primary and shield coil layers, which are connected at the ends by some radial surfaces. Half of the triangles making up the cryostat surface are shown so as to display the gradient coil surface. The target field region (pink) is partially obscured by the coil surface. Blue lines represent the mid-lines in x and z. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

minimise kjðwÞk1 þ bkjðwÞk1 subject to kBw  bt k1 =kbt k1 6 Dmax ;

ð9Þ

with b ¼ 200 and Dmax ¼ 10%. Using b ¼ 200 was chosen such that kjk1 and bkjk1 were approximately equal for all coils. In a real design situation, b should be chosen for each coil depending on manufacturing considerations. For the zonal coils, 129 sinusoidal functions were used with z variation only. For coils of first degree, 435 sinusoids were used with combinations of 29 sinusoids in z and 15 in /. For the XY shim (second degree coil) 261 sinusoidal functions were used with just 9

39

/ variations. These sinusoidal functions were chosen to be a set of symmetric or antisymmetric functions depending on the symmetry of the target field. The number of contour levels, N c , was fixed to 10.

2.6. Shielded Whole-body X-gradient coils The purpose of these simulations was to assess the performance of shielded, whole-body X-gradient coils. It was assumed that the coil will be made from wire of fixed cross-section, so it is expected

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Fig. 2. Performance metrics of 404 cylindrical unshielded X-gradient coils against maximum field error. (a) and (b) show how well the coils are optimised for R, (c) and (d) for kjk1 and s_ , (e) and (f) for kjk1 and c and (g) and (h) for L. All data are shown on logarithmic axes and marker colour corresponds to the type of optimisation used.

40

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

that the continuous current density power minimisation will not yield the optimal coil performance [28]. The total length of the coil surface was 1200 mm with an inner and outer diameter of 690 mm and 880 mm, respectively. The coil surface was modelled with 4288 flat triangular elements with 2208 mesh nodes. Wires may connect between inner and outer layers via 16 joining surfaces at each end. A target field region was used with 485 points distributed over a 500 mm  500 mm  400 mm diameter ellipsoid. Active shielding was enforced by modelling the instantaneous induced eddy currents in a cylindrical surface with length 1600 mm and 1000 mm diameter made of 2048 flat triangular elements. The coil surface geometry is depicted in Fig. 1. The optimisation problem was

minimise ckjðwÞkA;1 þ kjðwÞk1 subject to kBw  bt k1 =kbt k1 wT Ew

OR

wT Rc w

6 Dmax ;

ð10Þ

6 Q max :

In Eq. (10), c defines the importance placed upon sparsifying the current density at the expense of increasing the maximum current density. In this work, c was varied between 1  103 and 1  101 in 50 logarithmic steps. wT Rc w is the traditional power minimisation term. The constraint Q max places a limit on the total amount of power deposited by the induced eddy currents in a surrounding surface that represents the cold shield of a superconducting magnet. Here, Q max ¼ 100 Wm=T was used throughout to provide T 1 good active shielding. The matrix E ¼ M Tec ðL1 e Þ Re Le M ec , where M ec is the mutual inductance matrix between the coil surface and the eddy current surface, Le is the self inductance matrix for the eddy current surface and Re is the resistance matrix of the eddy current surface [29]. The number of contours, N c , of the stream-function was chosen so as to generate a gradient field efficiency, g = 110 lT m1 A1. Coil performance parameters were simulated in the same way as in previous sections. In the impedance simulations, the coil was modelled with a wire thickness of 4 mm and a fixed width equal to ðs_  1Þ mm – i.e. a cross-sectional area of tw ¼ 4ðs_  1Þ mm2.

3. Results The coils produced from the solution to Eq. (6) with direct matrix inversion and cvx were identical, as expected. Furthermore, the coils produced from the solution to Eq. (7) with p ¼ 1 are qualitatively very similar to those produced with the algorithm in Ref. [15] although an exact comparison was not made since the coils are constrained by the maximum field error as opposed to the ‘2 -norm error in the prior work. 3.1. Basic regularisation behaviour In Fig. 2 the appropriate performance metrics are plotted for all coils against the maximum field error produced by the coils. In the first column of the figure (a, c, e and g) the performance metrics relate directly to those figures that were used in the optimisation; i.e. the performance metric calculated from the continu1 ous current density model. These metrics are R1 ; kjðwÞk1 1 ; kjðwÞk1 1 and L and reflect how well the resistance, maximum current density, current density sum and inductance were minimised. In the right-hand column of Fig. 2b, d, f and h the performance metrics are the figures-of-merit (FoMs) as simulated by Biot–Savart integration and using Fasthenry (with variable track width). If the continuous current density is a good approximation to the discrete wire model, and if the optimisation is appropriately formulated, then these results should exhibit similar outcomes. These FoMs are g2 =R; gs_ ; g=c and g2 =L and reflect the coil efficiency normalised by resistance, minimum wire spacing, total conductor length and inductance. The optimisation problem contained 993 free variables (the stream-function values at the mesh nodes) and was solved in 250, 318, 430 and 327 s for the resistance, kjk1 ; kjk1 and inductance minimised solutions, respectively. 3.2. Mixed regularisation The results of 303 X-gradient coils designed to have the same maximum field error with mixed regularisation types (Eq. (8))

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Fig. 3. Mixed regularisation results for X-gradient coils with Dmax ¼ 5%. L-curves of the trade-off between norms p ¼ 2 and p ¼ 1, (d), p ¼ 1 and p ¼ 2, (e), p ¼ 1 and p ¼ 1, (f) are shown. Examples of wire patterns are shown in (a), (b), (c), (g), (h), (i), with 1 quadrant drawn.

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

(a)

41

3.3. Wire patterns for mixed p ¼ 1 and p ¼ 1 regularised coils Six resistance minimised coils were designed with different target fields and another 6 were designed with a mixture of sparsity and maximum current density optimisation. The field accuracy was kept the same for all coils. Their wire patterns are shown in Fig. 5 where parts (a, c, e, g, i and k) have minimum resistance and parts (b, d, f, h, j and l) have minimum mixed p ¼ 1 and p ¼ 1-norm. Also, the fundamental discrete loop coils for each target field obtained from Ref. [3] are overlaid for comparison. 3.4. Shielded whole-body X-gradient coils

(b)

Two types of shielded whole-body X-gradient coils were designed; one coil with minimised resistance and 50 coils with varying mixtures of weighted p ¼ 1 and p ¼ 1 -norm regularisation. After simulating the coil efficiency and impedance, the mixed norm coil with the best g2 =R value was chosen for comparison. Table 1 and Fig. 6 show the performance properties and the wire patterns of the two coils.

4. Discussion

(c)

Fig. 4. Optimisation metrics and discrete wire simulated FoMs compared. Parts (a)– (c) compare the resistance, wire spacing and coil length FoMs respectively. Values obtained from coils in the trade-off between p ¼ 2 and p ¼ 1 are shown as 2 $ 1, etc.

are shown in Fig. 3d–f. The change in the norm values are plotted to give an L-curve, where each norm measure of optimality is normalised by the least optimal in that set. For example, in Fig. 3d kjk1 =K 21 is plotted against kjk2 =K 12 , where kjk1 is the maximum current density of each coil and K 21 ¼ kjðw2 Þk1 is the maximum current density of the coil optimised for minimum resistance. Exemplary coil patterns are shown in Fig. 3 for 6 coils optimised with mixed regularisation optimised principally and secondarily for: resistance and maximum current density, (a); sparsity and resistance, (b); maximum current density and sparsity, (c); maximum current density and resistance, (g); resistance and sparsity, (h); sparsity and maximum current density, (i). In order to assess if the optimisation has the desired effect on the coil performance, the simulated coil performance metrics, FoMs, are plotted against their continuous current density counterparts used in the optimisation in Fig. 4. The FoMs for the resistance, maximum current density and sparsity are plotted in parts a, b and c, respectively. The mean time taken to compute the optimal solutions with mixed p ¼ ½2; 1; p ¼ ½1; 2 and p ¼ ½1; 1 norms was 434, 440 and 389 s, respectively.

The convex optimisation framework used in this work is a generalisation of previous coil design methods. It relies on the linear relationship between the basis functions of w and the magnetic field and current density. If the desired properties of the coil can be written as convex functions of w, then the coil design problem can be solved deterministically. This is illustrated by incorporating the ‘1 -norm of the current density amplitude, kjk1 , into the optimisation. This norm is most well known for its use in image denoising [30], compressed sensing MRI reconstruction [31] and passive shimming of magnets [32]. The ‘1 -norm has also been used in MRI for B0 shimming [33]. Appendix A gives more detail on the use of ‘p -norms of j. Furthermore the convenience of existing convex optimisation tools was demonstrated by prototyping several different optimisation functions using cvx and sdpt3. Any physical model that provides linearity between the free variables and the optimisation functions can be used, and here we have used two BEMs; one with flat triangular elements and linear shape functions and another with a finite cylindrical surface and truncated sinusoidal shape functions. 4.1. Basic regularisation behaviour Simulations were performed to elucidate the most general behaviour of different optimisation functions by trading the maximum magnetic field error for P; kjk1 ; kjk1 and W. These parameters generally exchange monotonically such that an increase in the maximum field error permits a consequential decrease in P (blue dots in Fig. 2a), kjk1 (red dots in Fig. 2c), kjk1 (green dots in Fig. 2e) and W (yellow dots in Fig. 2g). The same effect was seen when the coils were simulated using discrete wires in that an increase in the maximum field error permits a consequential increase in g2 =R (blue dots in Fig. 2b), gs_ (red dots in Fig. 2d), g=c (green dots in Fig. 2f) and g2 =L (yellow dots in Fig. 2h). As expected, g2 =R for coils designed with minimum resistance are the highest. Similarly, each FoM is highest when optimised in the appropriate manner. Most similar are the coils designed using minimum resistance and inductance. The most different are the coils designed by minimising kjk1 and kjk1 in that the wire spacing FoM, gs_ , is lowest for sparse coils and g=c is lowest for coils designed with minimum maximum current density.

42

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 5. Winding patterns of coils that generate Z0 (a and b), Z (c and d), Y (e and f), Z2 (g and h), ZY (i and j) and XY (k and l) fields with only half of each coil drawn. The lefthand coil has minimum P and the right-hand coil has minimum kjk1 and kjk1 .

4.2. Mixed regularisation Fig. 3d is similar (but inverted) to Fig. 2c in Ref. [21]. For fixed field error the maximum current density can be exchanged for resistance in a smooth manner resulting in a typical L-curve. This shows that large reductions in maximum current density can be achieved for modest increases in resistance, and conversely, that large reductions in resistance can be obtained by permitting some small increase in maximum current density. As was shown previously [21], it is therefore often sensible to chose a coil with a mixture of p ¼ 2 and p ¼ 1 regularisation. Fig. 3e shows that only small decreases in kjk1 can be obtained for a significant increase in resistance, and (f) shows a greater curvature, primarily as a result of large K 11 and K 11 values (i.e. the optimal coils in each case are very different from each other). The wire patterns in Fig. 3 illustrate the typical patterns obtained for X-gradient coils. Coils with primarily minimum resistance (a and h) are smooth, have varying wire spacing and take up most of the available surface. Coils with minimised maximum current density (c and g) exhibit less smooth patterns take up most of

the surface due to the even and highly spread windings. Finally, the sparse coil windings (b and i) are smooth, but have sharp transitions from dense regions of winding and regions where there are no wires at all. Minimum mixed norm coils have favourable characteristics of each type. The evenly spaced concentric windings appear attractive for cheap hand-winding techniques. The inner and outer loops often appear to be slightly more separated than the other loops. This is due to the finite resolution of the basis functions that make up the stream function. The mesh triangles are either of a finite size or the sinusoidal functions have a maximum spatial frequency, meaning that sharp transitions are not possible despite being favoured by the optimisation. This could be remedied by forcing these loops to conform to concentricity with the other loops. Increasing the mesh resolution or number of shape function would help but greatly increase the computational cost as the number of free-variables increased. Alternatively, multigrid approaches could be used in which the size of the triangles in the mesh is reduced in regions where j j j varies greatly, or the by using the mesh Laplacian of w. Another approach would be to change the shape

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45 Table 1 Performance parameters for whole-body X-gradient coils designed with minimum resistance and minimum mixed p ¼ 1 and p ¼ 1-norms. The table is separated into 4 parts; input parameters, optimisation results (for 1 Tm1), discrete wire and Biot– Savart simulated data and Fasthenry simulated results. Regularisation type

Resistance

ðckjkA;1 þ kjk1 Þ

c

– 33

1.265 30

kBw  bt k1 =kbt k1 (%)

12.0

12.0

wT Rc w

1:38  106

Nc

T

w Lc w kjðwÞk1

4

1:63  106 7:9  104

6:9  10

12:5  105

4:89  105 2090

max DBz (%) max Bz;cryo (%) P cryo (Wm/T) s_ (mm) c (m) L (lH) R (mX)

109.9 8.4 2.7 100 6.5 438 1131 259

110.2 10.8 2.5 100 11.3 437 991 131

g2 =L (T2 m2 A2 H1) g2 =R (T2 m2 A2 X1)

1:07  105

1:22  105

4:7  108

9:2  108

Calculation time (s)

g (lT m1 A1)

5

7:8  105

4:95  10 628

kjðwÞkA;1

functions to some basis that could possess higher spatial frequencies such as wavelets. The data presented in Fig. 4 demonstrate that the convex minimisation terms (plotted on the horizontal axes) have the desired effect on the coil properties, for fixed field error. In the case of the sparse coils, the aim is to minimise the total wire length by minimising the absolute sum of current densities, which is shown in Fig. 4c to be appropriate. A small kink appears in these data indicating that for some range of coils, a decrease in kjk1 actually leads to an decrease in coil efficiency per meter of wire. This could be due to either an increase in c or decrease in g over this small range of coils. 4.3. Wire patterns for mixed p ¼ 1 and p ¼ 1 regularised coils A visual comparison between winding patterns of coils designed with minimum resistance and minimum mixed norm is

(a)

43

possible with the results in Fig. 5. Some interesting findings were obtained in that for some coils the winding pattern tends toward fundamental, well-known solutions [3]. The Z0 coil shown in Fig. 5b has become a simple solenoid whose length is as short as possible whilst satisfying the field accuracy constraint and providing the specified weighting between wire spacing and sparsity. Although the data are not shown, if the sparsity weighting is increased, then the wires become closer, the solenoid splits in two and tends towards the Helmholtz coil solution (shown in green). The field error constraint is 10% and if this is reduced, then the solenoid splits into more separate bunches reminiscent of designs for highly homogeneous MRI magnets [1], albeit just in the axial direction rather than the axial and radial directions. It may be profitable to use this technique as a starting point for the design of superconducting magnets. A similar effect is seen for the Z (c and d) and the Z2 (g and h) coils. The wires of the Z coils are centred around the position of the Maxwell pair (shown in green), and the mixed kjk1 and kjk1 coil has evenly spaced windings for a finite portion of the available space. The Z2 coils do not match the positions of the designs from Ref. [3], but it is likely that the field accuracy is different. The Y coil windings (e and f) change from the well known ‘‘fingerprint’’ pattern to a concentric series of windings; evenly spaced with large wire-free areas. The shape of these loops is interesting in that they are not similar to the rectangular type (green) that is a solution from Ref. [3] or the well-known solutions using the method in Ref. [27]. Instead they are in the form of a squashed oval. The ZY coil windings (i and j) become dramatically more sparse taking up only a small portion of the available space. This follows closely the shape of the discrete loop coil from Ref. [3]. The windings are mostly evenly spaced and concentric, but the smaller central loops are more spaced out because they must maintain Bz ðrÞ / zy within 10%. In the limit of a purely sparsity promoting optimisation and poor field accuracy so that just one loop is yielded, this type of optimisation may also be useful in obtaining RF coil patterns. In the XY coil (k and l) the windings are similar to the Y coil in that they are slightly squashed ovals, but in order to maintain the 10% accurate field, 4 loops become separated from the main concentric loops. The gap between is caused by the sparse optimisation target. Although not shown the same effect was seen in the Y coil when the sparse weighting was increased and/or when the field constraint was made more restrictive.

(b)

Fig. 6. Winding patterns of two whole-body shielded X-gradient coils designed with minimum P, (a), and mixed kjkA;1 and kjk1 , (b). One quarter of the winding pattern is shown and the z axis is an unwrapped hybrid ðzÞ  ðqÞ  ðzÞ axis showing the connections between primary and shield surfaces.

44

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

4.4. Shielded whole-body X-gradient coils With these examples, the simulated performance shows that the coil efficiency normalised by the heat dissipation can be increased significantly if the coil is to be made from fixed width wires. In this case the width of wire used to wind the coil is limited by the minimum wire spacing. When simulated with Biot–Savart and Fasthenry this result is corroborated by an increase in g2 =R of 96%. Therefore it is estimated that just half the amount of heat would be generated whilst providing the same gradient strength. This is primarily due to the increase in s_ that can be achieved via minimising the maximum current density, but also through sparsification of the coil windings. Connections between primary and shield layers in X-gradient coils have been previously presented [34]. Here, space is left for incorporating passive shim trays, but also ensures that the connections between primary and shield are purely radial for easy connection. The mixed norm winding pattern of the primary layer is similar to Figs. 3c and 5f in that the wires form concentric loops, with some wire-free areas, while the shield windings are smooth and unevenly spaced like the resistance minimised coil. This is because shield wires must accurately cancel the field from the primary coil at the magnet surface and so must be a spatially smoothed version of the primary wire pattern. Sparsity cannot be promoted and the maximum current density need not be reduced. The large improvement in g2 =R using the mixed kjkA;1 and kjk1 minimisation were obtained in simulation. It is important to fully test these results by manufacturing this type of coil and assess its performance in a real MRI experiment. These results are intended for coils wound from fixed width wire(s) and are not applicable to coils built with variable track widths using, for example, water-jet cutting construction techniques. 5. Conclusion The inverse boundary element method and generalised convex optimisation techniques were combined to provide a more flexible framework to design gradient and shim coils. This new method was used to design coils by minimising different types of norms of the current density. The general performance characteristics as well as the behaviours of the winding patterns of such coils were investigated using cvx and sdpt3. The ‘1 -norm of the current density when minimised promotes sparsity of the windings, whereas minimising the ‘1 -norm causes the wires to become more spread out, and evenly so, as previously reported. A mixture of norms provides coils with favourable characteristics and we show that a mixture of ‘1 - and ‘1 -norm minimisation provides interesting coils whose windings are concentric and compact. This may be of use when constructing coils by hand from fixed cross-section wire because it provides a closer approximation to the real resistance than traditional continuous current density approximations. A whole-body X-gradient coil was designed using this technique that was simulated to produce almost half the heat as a standard resistance minimised coil for the same gradient field. Other forms of optimisation could be proposed using this framework by altering the minimisation and constraint functions to convex functions that reflect other desirable characteristics of gradient (e.g. less vibration [35] or lower peak temperature [18]), shim or other types of coil. Acknowledgments We kindly acknowledge funding by the BMBF (Grant No. 13N9121). The first author is a gratefully recipient of a Marie Curie

International Incoming Fellowship of the European Community’s Seventh Framework Programme under the project number FP7PEOPLE-2012-IIF-331751: MCODE. Appendix A. ‘p -Norms of the current density in coil design Consider the ith triangular element, PQR, of a surface shown in Fig. A1a with its base, PQ, defined by the vector ei and height di . The stream-function, wi ðrÞ, is a linear function in the triangle defined by the nodal values, wiP ; wiQ ; wiR . In the following we define wiQ ¼ wiP for simple visualisation of the problem, but since current density in each triangular element is a linear sum of the three current density basis functions [9], the calculation below is general. The current density in this triangle is given by ji ¼ ei Dwi = j ei jj di j¼ ei Dwi =2Ai , where Dwi ¼ wiR  wiP , the current density magnitude, ji ¼j ji j and the area is Ai . Converting current density into current carrying conductors (wires or tracks) with rectangular cross-section, we use N i wires, such that each wire carries I Amps of current where I ¼ Dwi =N i . Finally, we end up with N i wires in the triangular element with thickness t, width wi , spacing si and gap between conductors, g i , where each wire has a length, lin . The power dissipated in this element is P i ¼ I2 Ri , and Ri is given by

Ri ¼

Ni q X

twi

lin :

Combining Eq. (A1) with above paragraph,

Ri ¼

ðA1Þ

n¼1

PNi

n lin

¼j ei j N i =2 and equations in the

q Ai jj : twi I

ðA2Þ

A.1. Case 1: Coils built from CNC cut copper sheet Thin cuts are made in copper sheets to ensure the correct current flow in the coil and so we assume g i ! 0 and wi ! si (see Ref. [28]). We know in this case wi ¼j di j =N i ¼ I=ji so the resistance of this element becomes the familiar

Ri ¼

qAi tI2

2

ji :

ðA3Þ

Equation (A3) justifies the minimisation of the ‘2 -norm of the current density in order to obtain coils of minimum resistance and therefore minimum power dissipation;

P ¼ I2

X qX 2 q 2 qT Ri ¼ Ai ji ¼ kjkA;2 ¼ j diagðAÞj; t t t i i

ðA4Þ

where j is a vector of the current density magnitudes, ji , in each element. A is a vector of the areas of each element, Ai , and the weighted ‘p -norm is

(a)

(b)

Fig. A1. Illustration of the conversion from current density in a triangle to current carrying wires. In (a), an example triangular element, PQR, of a surface has a current density, ji , flowing in it. The physical implementation, (b), of this is a set of current carrying wires or copper tracks.

M.S. Poole, N. Jon Shah / Journal of Magnetic Resonance 244 (2014) 36–45

kxka;p :¼

M X

!1=p am j xm jp

:

ðA5Þ

m¼1

In practice, our optimisation variables are w where the transform from w ! j is convex. Eq. (A4) can also be written in quadratic form with w (Eq. (4)): P ¼ wT Rc w, where Rc is sparse, square, symmetric, positive semidefinite but not diagonal. Similarly, the stored energy (Eq. (3)), W ¼ wT Lc w where Lc is full, square, symmetric and positive semidefinite, and therefore diagonalisable. Both P and W can thus be expressed as a weighted ‘2 -norm after basis transformation. A.2. Case 2: Coils built from fixed width wire The coil of minimum resistance made from fixed width wire has w ! mini si , where mini si ¼ I=maxi ji . Then, Eq. (A2) becomes

Ri ¼

q Ai tI2

ji max ji : i

ðA6Þ

Now the power dissipation of the whole coil becomes

P ¼ I2

X X q q Ri ¼ max ji Ai ji ¼ kjk1 kjkA;1 i t t i i

ðA7Þ

The product kjkA;1 kjk1 is not convex, but if we minimise

ckjkA;1 þ kjk1 then we will minimise kjkA;1 kjk1 for some value of c, and therefore P will be minimal for coils built with wire whose width is constant and equal to the minimum wires spacing in the coil design. When all triangle areas are equal, which is the case in Sections 2.3, 2.4 and 2.5, diagðAÞ is proportional to an identity matrix and kjkA;p / kjkp . References [1] J. Bird, D. Houlden, N. Kerley, D. Rayner, D. Simkin, S. Pittard, Design considerations for ultra high field mri magnet systems, in: Ultra High Field Magnetic Resonance Imaging, Biological Magnetic Resonance, vol. 26, Springer, US, 2006, pp. 19–43. [2] R. Turner, Gradient coil design: a review of methods, Magn. Reson. Imaging 11 (7) (1993) 903–920. [3] F. Roméo, D.I. Hoult, Magnet field profiling: Analysis and correcting coil design, Magn. Reson. Med. 1 (1) (1984) 44–65. [4] S.E. Ungersma, H. Xu, B.A. Chronik, G.C. Scott, A. Macovski, S.M. Conolly, Shim design using a linear programming algorithm, Magn. Reson. Med. 52 (3) (2004) 619–627. [5] S. Crozier, L.K. Forbes, D.M. Doddrell, The design of transverse gradient coils of restricted length by simulated annealing, J. Magn. Reson. Ser. A 107 (1) (1994) 126–128. [6] S. Shvartsman, M.C. Steckner, Discrete design method of transverse gradient coils for mri, Concepts Magn. Reson. Part B: Magn. Reson. Eng. 31B (2) (2007) 95–115. [7] M.J.E. Golay, Field homogenizing coils for nuclear spin resonance instrumentation, Rev. Sci. Instrum. 29 (4) (1958) 313–315. [8] W. Edelstein, F. Schenck, Current Streamline Method for Coil Construction, US Patent 4,840,700, 1989. [9] S. Pissanetzky, Minimum energy mri gradient coils of general geometry, Meas. Sci. Technol. 3 (7) (1992) 667–673.

45

[10] Z.J.J. Stekly, Continuous, transverse gradient coils with high gradient uniformity, in: Proceedings of the Society for Magnetic Resonance in Medicine, vol. 4, 1985, p. 1121. [11] J.W. Carlson, K.A. Derby, K.C. Hawryszko, M. Weideman, Design and evaluation of shielded gradient coils, Magn. Reson. Med. 26 (2) (1992) 191–206. [12] D.I. Hoult, R. Deslauriers, Accurate shim-coil design and magnet-field profiling by a power-minimization-matrix method, J. Magn. Reson. Ser. A 108 (1) (1994) 9–20. [13] R. Turner, Minimum inductance coils, J. Phys. E: Sci. Instrum. 21 (10) (1988) 948–952. [14] R.A. Lemdiasov, R. Ludwig, A stream function method for gradient coil design, Concepts Magn. Reson. Part B: Magn. Reson. Eng. 26B (1) (2005) 67–80. [15] M. Poole, P. Weiss, H.S. Lopez, M. Ng, S. Crozier, Minimax current density coil design, J. Phys. D: Appl. Phys. 43 (9) (2010) 095001. 13pp. [16] P.T. While, L.K. Forbes, S. Crozier, Designing gradient coils with reduced hot spot temperatures, J. Magn. Reson. 203 (1) (2010) 91–99. [17] C.T. Harris, W.B. Handler, B.A. Chronik, Electromagnet design allowing explicit and simultaneous control of minimum wire spacing and field uniformity, Concepts Magn. Reson. Part B: Magn. Reson. Eng. 41B (4) (2012) 120–129. [18] P.T. While, M.S. Poole, L.K. Forbes, S. Crozier, Minimum maximum temperature gradient coil design, Magn. Reson. Med. 70 (2013) 584–594. [19] CVX Research, Inc., CVX: Matlab Software for Disciplined Convex Programming, version 2.0, (August 2012). [20] M. Grant, S. Boyd, Graph implementations for nonsmooth convex programs, recent advances in learning and control (a tribute to m. vidyasagar), in: S.B.V. Blondel, H. Kimura (Eds.), Lecture Notes in Control and Information Sciences, Springer, 2008, pp. 95–110. [21] M.S. Poole, P.T. While, H. Sanchez Lopez, S. Crozier, Minimax current density gradient coils: analysis of coil performance and heating, Magn. Reson. Med. 68 (2012) 639–648. [22] C. Cobos Sanchez, M. Fernandez Pantoja, M. Poole, A. Rubio Bretones, Gradientcoil design: a multi-objective problem, IEEE Trans. Magn. 48 (6) (2012) 1967– 1975. [23] M. Poole, R. Bowtell, Novel gradient coils designed using a boundary element method, Concepts Magn. Reson. Part B: Magn. Reson. Eng. 31B (2007) 162– 175. [24] S. Boyd, L. Vandenberghe, Convex Optimisation, Cambridge University Press, 2004. [25] K.C. Toh, M.J. Todd, R.H. Tutuncu, Sdpt3 — a matlab software package for semidefinite programming, Optim. Methods Software 11 (1999) 545–581. [26] M. Kamon, M.J. Tsuk, J.K. White, Fasthenry: a multipole-accelerated 3-d inductance extraction program, IEEE Trans. Microwave Theory Tech. 42 (1994) 1750–1758. [27] H. Siebold, Gradient field-coils for mr imaging with high spectral purity, IEEE Trans. Magn. 26 (2) (1990) 897–900. [28] P.T. While, J.G. Korvink, N.J. Shah, M.S. Poole, Theoretical design of gradient coils with minimum power dissipation: Accounting for the discretization of current density into coil windings, J. Magn. Reson. 235 (2013) 85–94. [29] G.N. Peeren, Stream function approach for determining optimal surface currents, J. Comput. Phys. 191 (1) (2003) 305–321. [30] L.I. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Phys. D: Nonlinear Phenom. 60 (14) (1992) 259–268. [31] M. Lustig, D. Donoho, J.M. Pauly, Sparse mri: the application of compressed sensing for rapid mr imaging, Magn. Reson. Med. 58 (6) (2007) 1182–1195. [32] B. Dorri, Method of Optimizing Passive Shim Placement in Magnetic Resonance Magnets, US Patent 5,045,794, 1991. [33] J. Lee, M. Lustig, D. hyun Kim, J.M. Pauly, Improved shim method based on the minimization of the maximum off-resonance frequency for balanced steadystate free precession (bssfp), Magn. Reson. Med. 61 (6) (2009) 1500–1506. [34] S. Shvartsman, M. Morich, G. Demeester, Z. Zhai, Ultrashort shielded gradient coil design with 3d geometry, Concepts Magn. Reson. Part B: Magn. Reson. Eng. 26B (1) (2005) 1–15. [35] G. Li, C. Mechefske, Structural-acoustic modal analysis of cylindrical shells: application to mri scanner systems, Magn. Reson. Mater. Phys. Biology Med. 22 (6) (2009) 353–364.