Kevin Bleakley
email: kevbleakley@gmail.com

Jean-Philippe Vert
email: jean-philippe.vert@mines-paristech.fr

Joint segmentation of many aCGH profiles using fast group LARS

Joint segmentation of many aCGH profiles using fast group LARS

Introduction

Array-based Comparative Genomic Hybridization (aCGH) is a technique that aims to detect chromosomal aberrations on a genomic scale in a single experiment. In tumors for example, chromosomal aberrations include tumor suppressor genes being inactivated by deletion, and oncogenes being activated by duplication. Such changes to the underlying "normal" genomic state, known as copy number variations (CNVs) provide information related to the current disease state of the individual.

Many cancers are known to exhibit recurrent CNVs in specific locations in the genome, for example monoploidy of chromosome 3 in uveal melanoma [START_REF] Speicher | Chromosomal gains and losses in uveal melanomas detected by comparative genomic hybridization[END_REF], loss of chromosome 9 in bladder carcinomas [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF], loss of 1p and gain of 17q in neuroblastoma [START_REF] Bown | 17q gain in neuroblastoma predicts adverse clinical outcome. U.K. cancer cytogenetics group and the U.K. children's cancer study group[END_REF][START_REF] Van Roy | Localization of the 17q breakpoint of a constitutional 1;17 translocation in a patient with neuroblastoma within a 25-kb segment located between the accn1 and tlk2 genes and near the distal breakpoints of two microdeletions in neurofibromatosis type 1 patients[END_REF] and amplifications of 1q, 8q24, 11q13, 17q21-q23 and 20q13 in breast cancers [START_REF] Yao | Combined cDNA array comparative genomic hybridization and serial analysis of gene expression analysis of breast tumor progression[END_REF].

aCGH allows rapid mapping of CNVs of a tumor sample at the genomic level [START_REF] Pinkel | High resolution analysis of DNA copy number variation using comparative genomic hybridization to microarrays[END_REF]. Originally, the technique was based on arrays with several thousand large insert clones (e.g., bacterial artificial chromosomes or BACs) with a megabase resolution. More recently, this has been improved using oligonucleotide-based arrays with up to hundreds of thousands of probes, bringing the resolution down to several kilobases [START_REF] Gershon | DNA microarrays: more than gene expression[END_REF]. Each spot on an aCGH array contains amplified or synthesized DNA of a particular region of the genome. The array is hybridized with DNA extracted from a sample of interest, and in most cases with (healthy) reference DNA. Each of the two samples is then labeled with its own fluorochrome and the ratio of fluorescence between the two samples is expected to reveal the ratio of DNA copy number at each position of the genome. Usually the logarithm of this ratio is then taken and the values (profile) plotted linearly along the genome, as shown for example in Fig. 5(b)-(c).

One challenge related to aCGH profiles is to detect regions of constant copy number, separated by breakpoints, in the presence of significant amounts of noise. Many groups have attempted to answer this question when treating a single profile, with various methods for segmentation and smoothing of 1-dimensional profiles [START_REF] Olshen | Circular binary segmentation for the analysis of array-based DNA copy number data[END_REF][START_REF] Hupé | array CGH data: from signal ratio to gain and loss of DNA regions[END_REF][START_REF] Picard | A statistical approach for array CGH data analysis[END_REF][START_REF] Fridlyand | Hidden markov models approach to the analysis of array CGH data[END_REF][START_REF] Wang | A method for calling gains and losses in array CGH data[END_REF][START_REF] Picard | A segmentation-clustering problem for the analysis of array CGH data[END_REF][START_REF] Huang | Robust smooth segmentation approach for array cgh data analysis[END_REF][START_REF] Tibshirani | Spatial smoothing and hot spot detection for CGH data using the fused lasso[END_REF]. Recently, approaches have been suggested for dealing with multiple aCGH profiles. Indeed, we often have experimental data from a series of patients who have the same disease status, or several groups of patients, each with a particular disease status. The goal of jointly treating several profiles is to extract breakpoints and copy number variations that are globally representative of the disease status of those patients. The STAC algorithm [START_REF] Diskin | STAC: a method for testing the significance of DNA copy number aberrations across multiple array-CGH experiments[END_REF] uses a novel search strategy to find regions of copy number gains and losses that are statistically significant across a set of profiles. It assigns p-values to each location on the genome using a multiple testing corrected permutation approach. The approach of [START_REF] Klijn | Identification of cancer genes using a statistical framework for multiexperiment analysis of nondiscretized array CGH data[END_REF] is based on kernel regression. One disadvantage of this is that results are dependent on the width of the kernel window, so a region may exhibit a statistically significant copy number variation for one width, but not for another, leading to difficulties in interpretation. The approach in [START_REF] Rouveirol | Computation of recurrent minimal genomic alterations from array-CGH data[END_REF] is to first discretize profiles independently into gains, losses and normal. Then, they look for regions in which a significant number of the profiles share gains or losses using a Boolean framework and theoretical results for pattern searching. This method potentially loses pertinent information in the discretization step. Finally, the approach in [START_REF] Picard | Joint segmentation of multivariate Gaussian processes using mixed linear models[END_REF] is to simultaneously segment a set of profiles using mixed linear models to account for both covariates and correlations between probes. They proposed an EM-type algorithm that uses dynamic programming during the segmentation step. Though the computational cost was not provided, the algorithm appears computationally intensive.

In this article, we present an algorithm that jointly renders a set of n profiles piecewise-constant, with the intended goal of uncovering a set of breakpoints and regions of constant copy number that are pertinent with respect to the whole set of profiles. This algorithm generalizes a Fused Lasso-type algorithm that was used to find breakpoints (and smooth) a single profile in [START_REF] Harchaoui | Catching change-points with lasso[END_REF]; there, they transformed the problem into a standard Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] and proposed an extremely fast way to find the whole solution path, that is, an ordered set of pertinent breakpoints. When several profiles must be segmented simultaneously, we propose to impose a fused constraint jointly on the set of profiles, leading to a reformulation of the problem as a group LARS or group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. The originality of our approach lies in the fact that the fused constraint imposes that each newly-selected breakpoint be shared by all profiles. In this formulation, we end up with a stupendously large matrix (with p 2 n 2 entries) to work with, but we show that to run a group LARS algorithm it is not necessary to store this matrix.

In order to find k breakpoints in n profiles of p probes, our method has a complexity O(knp) in time and requires O(np) in memory. In speed trials using Matlab on a 2008 Macbook Pro with 4GB of RAM, for n = 20 profiles of p = 2000 probes, 50 breakpoints were found in 0.72 seconds, an average of 0.014 seconds each. Even at the limit of current aCGH technology, with p ∼ 946 000 and, for example n = 20, the first 50 selected breakpoints were obtained in 304 seconds, an average of about 6 seconds each. We present the performance of the algorithm on simulated data under various noise conditions, and demonstrate that it is able to recover the breakpoints shared by some or all of the profiles as the number of profiles increases. Finally we implement the algorithm to segment aCGH profiles of bladder cancer patients.

Methods

Suppose we have n aCGH profiles each of length p. The p probe values are calculated at identical locations for each of the n profiles. For each profile, we want to find a piecewise constant representation, where the jumps between constant segments represent breakpoints, that is, places where the copy number changes.

One profile, one chromosome

Our starting point is the following framework for finding breakpoints in a onedimensional piecewise-constant signal with white noise, as introduced by [START_REF] Harchaoui | Catching change-points with lasso[END_REF]. Let Y = (y 1 , . . . , y p ) be the observed signal. Consider the following constrained optimization problem:

min β∈R p β -Y 2 2 subject to p i=1 |β i -β i-1 | < µ, (1)
where µ is a fixed non-negative constant and by convention β 0 = 0. The constraint p i=1 |β iβ i-1 | < µ can be seen as a convex relaxation of bounding the number of jumps in β, an idea also implemented in the fused Lasso [START_REF] Tibshirani | Spatial smoothing and hot spot detection for CGH data using the fused lasso[END_REF]. When µ is small enough, this constraint causes the solution β to be made up of runs of equally-valued β i separated by an occasional jump from one constant to another, i.e., a piecewise constant function. For µ large enough the constraint is no longer effective and the solution is merely β = Y . This problem is convex, meaning that any standard convex optimization package can solve it for a given µ or a sequence {µ j } J j=1 , though this remains computationally intensive. Making the change of variable

u 1 = β 1 , u 2 = β 2 -β 1 , . . . , u p = β p -β p-1 , [7] rewrite (1) as min u∈R p Au -Y 2 2 subject to p i=1 |u p | < µ,
where A is a p × p matrix whose lower-diagonal and diagonal are 1 and upper-diagonal is 0. This is exactly a Lasso regression problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], whose complete solution path for µ can be obtained for example by the lars package [START_REF] Efron | Least Angle Regression[END_REF] in a matter of seconds if p is of the order of several hundred. However, once p gets into the thousands, implementations like lars, which require storing a p × p matrix and that have complexity O(k 3 + pk 2 ) for calculating the first k breakpoints, greatly slow down or do not run. It was shown in [START_REF] Harchaoui | Catching change-points with lasso[END_REF] that the algorithm can be reformulated to take advantage of the structure of the matrix A without storing it, resulting in a computation time in O(pk) with only the need to store vectors of length p.

Many profiles, one chromosome

Let us now consider a n × p matrix Y containing n profiles of length p. Our aim is to apply a similar procedure to the one profile case, but jointly to the whole set of profiles. We propose the following constrained optimization problem:

min β∈R n×p β -Y 2 2 subject to p i=1 β i -β i-1 2 < µ, (2)
where µ is a fixed non-negative constant, β i is the vector of length n containing the value of the i th probe for each of the n individuals and by convention, β 0 = 0. The choice of L 2 norm in the constraint has the effect of creating long runs of consecutive vectors β i that are equal, interspersed with occasional jumps from one "constant" vector to the next [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. This mirrors the one profile case with its piecewise constant approximation, but here it is n profiles jointly. Intuitively, places where the set of profiles jointly "break" tend to be places where a large number of individual profiles exhibit a break. We now introduce a practical framework for solving such a problem. As for the one profile case [START_REF] Harchaoui | Catching change-points with lasso[END_REF], we make the change of variable

u 1 = β 1 , u 2 = β 2 -β 1 , . . . , u p = β p -β p-1 ,
where all these objects are now n-dimensional vectors. This gives us the representation:

min u∈R np Au -Y 2 2 subject to p i=1 u i 2 < µ,
where u is the np-dimensional vector of the p vectors u i of length n stacked on top of each other, by momentary abuse of notation Y is the np-dimensional vector of the columns of Y stacked on top of each other and A is now as follows: for each 1 in the matrix A of the one profile case, replace it with an n × n identity matrix and for each 0, replace it with an n × n matrix of zeros. The matrix A thus becomes an np × np matrix.

In fact, this can be rewritten as a group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], i.e., min

u∈R np p i=1 A i u i -Y 2 2 subject to p i=1 u i 2 < µ, (3)
where A i is the matrix of size np × n of the columns n(i -1) + 1 up to ni of A, and u i the i th column of u. Here, each group i is the set of n variables, one from each profile, found in position i on the genome. The group Lasso and group LARS algorithms [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] select variables "in groups" rather than one by one as for Lasso and LARS. This is useful when we have prior knowledge as to relationships between subsets of variables. In our case, the relationship is that each group represents the changes in value of n profiles between two adjacent locations on the genome. Selection of a group by the algorithm corresponds to choosing a location on the genome where a significant amount of change happens to a significant number of profiles. Whereas the standard algorithms for solving Lasso and LARS are almost identical, generalizations to group Lasso and group LARS are less so. The group LARS algorithm we describe here is therefore not the solution path to (3), but to a similar problem. We chose to work with group LARS as similar ideas to those in [START_REF] Harchaoui | Catching change-points with lasso[END_REF] could be used to implement an extremely fast algorithm.

"Many profiles, one chromosome" fast group LARS

The group LARS algorithm we propose explicitly follows the steps given in [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], but avoids the suggested matrix formulation by generalizing the methodology of [START_REF] Harchaoui | Catching change-points with lasso[END_REF]. To find breakpoints, we follow the following steps, which have computational complexity O(np), resulting in a complexity O(npk) to find the first k breakpoints.

• Step 1: start with u n×p = 0, k = 1 and r n×p = Y .

• Step 2: compute the currently most correlated set:

A 1 = arg max i A ′ i r 2 2 /p i .
By abuse of notation, r has been momentarily written as an np-dimensional vector of the columns of r stacked on top of each other. As ∀i, j, p i = p j = n here, we can ignore the division by p i . Due to the structure of A, to calculate the most correlated set it suffices to calculate the cumulative sums by row of r n×p , starting in the p th column and moving back along the rows to the 1 st column. We call this resulting n × p matrix C. The column of C with the largest L 2 norm corresponds to the first selected group, i.e., breakpoint. This step takes O(np) operations, and stores a n × p matrix.

• Step 3: compute the current direction γ. We represent γ as an n × p matrix. Each column of γ whose index is not in the active set A k is identically zero. Calculating the rest of γ usually requires having to calculate

G nk A k = (A ′ A k A A k ) -1
, the inverse of an nk × nk matrix. In fact, due to the analogous structure of A with A from the one profile case in [START_REF] Harchaoui | Catching change-points with lasso[END_REF], it suffices to calculate the inverse

G A k = G k
A k of the relevant k × k matrix in the one profile case, which is a tridiagonal matrix with a closed form expression [START_REF] Harchaoui | Catching change-points with lasso[END_REF]. Then, with the generic notation B[:, i] meaning the i th column(s) of any matrix B, the remaining non-zero columns of γ are obtained in O(nk) by the matrix multiplication:

γ[:, A k ] = C[:, A k ] × G ′ A k .
• Step 4: for all i / ∈ A k , compute how far the group LARS algorithm will progress in direction γ before A i enters the most correlated set. This can be measured by an α i ∈ [0, 1] such that

A ′ i (r -α i Aγ) 2 2 = A ′ i ′ (r -α i Aγ) 2 2
, where i ′ is arbitrarily chosen from A k and r and γ are again momentarily given as np-dimensional vectors for notational simplicity. In [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], it was shown that α i is always well-defined. This equation is quadratic in α i . To solve it simultaneously for all i in O(np):

1. define γ * as the n × p matrix we get by taking the row-wise cumulative sums of γ, starting this time at column 1 and finishing at column p. 2. define δ * as the n × p matrix we get by taking the row-wise cumulative sums of γ * , starting at column p and finishing at column 1. 3. the set of p quadratic equations is symbolically aα 2 + bα + c, with α = (α 1 , . . . , α p ). Let A k (j) mean the j th index in the active set, when sorted into ascending order. Then: 4. to calculate a, first let

a * = δ * . × δ * -(δ * [:, A k (1)]. × δ * [:, A k (1)]) × 1 1×p ,
where '×' means matrix multiplication and '.×' means element by element multiplication. a is then the vector of length p given by the column-wise sums of a * . 5. to calculate b,

let b * = C. × δ * -(C[:, A k (1)]. × δ * [:, A k (1)]) × 1 1×p . Then, b is -2 times the column-wise sums of b * . 6. to calculate c, let c * = C. × C -(C[:, A k (1)]. × C[:, A k (1)]) × 1 1×p
. Then, c is the column-wise sums of c * . 7. the quadratic equation can then be solved.

• Step 5: select the smallest non-negative α j such that j / ∈ A k , then update the active set to include this new member: A k+1 = A k ∪ j.

• Step 6: perform in O(np) the updates: u = u + α j γ, r = Y -Aγ (with γ momentarily written as an np-dimensional vector), k = k + 1 and update C with respect to the new r as was done in Step 2. Then go back to Step 3 and repeat until either p iterations are performed or a chosen number K max < p of breakpoints are found.

Many profiles, many chromosomes

CGH profiles usually span many chromosomes. Constraints linking the last data point on each chromosome to the first data point on the next are meaningless. Removing these constraints means that the matrix A must be changed in the following way: replace 1s with 0s whenever they corresponds to entries of the matrix indexed by two different chromosomes. This has minor follow-on effects in the algorithm, including calculation of cumulative sums and the matrix G A k , and the updating of r. The speed of the algorithm is unaffected.

Experiments

Speed trials

All trials used Matlab on a 2008 Macbook Pro with 4GB of RAM. Fig. 1(a) indicates linearity in n, and 50 breakpoints were found in 0.72 seconds, an average of 0.014 seconds each. Fig. 1(b) shows linearity in p. Fig. 1(c) shows, for n and p fixed, a near-linear relationship in k, i.e., subsequent breakpoints do not take longer to find than earlier ones. This confirms the theoretical O(npk) complexity.

To test the speed of the algorithm at the current limit of aCGH technology, we performed speed trials using sample data freely provided by Affymetrix (http://www.affymetrix.com) for the Affymetrix Genome-Wide Human SNP Array 6.0, which includes around 946 000 copy number probes. The dataset includes 5 sets of 5 replicates. Three of the five individuals have abnormal copies of the X chromosome, with 3, 4 and 5 copies respectively. We calculated the log-ratio of 20 profiles (from 5 replicates on 4 individuals) against the first profile of the first (normal) individual after removing Y chromosome probes, leaving 937 223 probes per profile.

Again, the algorithm ran linearly in k. For 10 profiles, the group LARS algorithm took 3.4 seconds to find each subsequent joint breakpoint, and for 20 profiles, 6.1 seconds. Thus, the algorithm is computationally practical, extremely fast even, at the upper limits of current technology. We remark that for these 20 profiles, the fast group LARS algorithm correctly selected as the most important joint breakpoint (out of the 937 223 possible choices) the jump from 0 to the first probe on the X chromosome.

Performance on simulated data

We performed a series of simulations in order to verify that the algorithm behaved well with respect to our stated goals, namely, recover breakpoints shared by several of a set of profiles. We designed four experiments to move gradually from an artificial to a more realistic setting.

1. All profiles share the same breakpoints. We simulated profiles of length 1000, each divided into ten constant-valued segments of length 100. Hence, there are 10 breakpoints, located before probes 1, 101, 201, . . . , 901. The constant value of each of the 10 segments was randomly drawn from a uniform distribution on [-1, 1]. White noise from a N (0, σ2 ) was then added independently to each of the 1000 probe values. As shown in Fig. 2, we simulated with varying levels of noise: σ 2 ∈ {0.01, 0.1, 0.2, 0.5}. In particular, we see that with σ 2 = 0.5, the noise is often significantly larger than the distance between subsequent underlying constant segments. For a given value of σ 2 , we randomly generated one profile in this way.

We then asked the fast group LARS algorithm to find 10 breakpoints only. If these ten breakpoints corresponded exactly to the ten real breakpoints, we stopped. Otherwise, we added a second randomly generated profile and ran the group LARS jointly on these two, and so on. For each of 1000 such trials, we calculated the number of profiles needed to find the 10 real breakpoints as its first 10 predictions. A good algorithm should correctly select these breakpoints, given enough profiles.

located at exactly the same probe on each profile. This is the same experimental condition as (1) except that each breakpoint that was fixed at location i can now be at any of {i-2, i-1, i, i+1, i+2}, chosen uniformly. For simplicity, the breakpoint before probe 1 was kept fixed. For a given σ 2 , we then run the algorithm as in [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF], though in each loop we initially find the first 50 ordered breakpoints. We stop adding profiles one by one when the following event happens: every one of the first m ≥ 10 ordered predicted breakpoints is included in one of the 10 breakpoint zones and there is at least one predicted breakpoint in each of the 10 zones. i.e., we find all of the 10 zones before adding a non-existant breakpoint.

3. All profiles have a subset of a predefined set of breakpoints. Potential breakpoints are located in the 10 locations described in [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF]. First, the breakpoint before the first probe is automatically generated, i.e., the value of the constant segment from probe 1 to 100. Then, the potential breakpoint 2 before probe 101 is randomly included with probability 0.7. If it is included, we randomly select the value of the segment from 101 to 200 as in [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF]. Otherwise, the value of the segment from 1 to 100 is continued from 101 to 200. We iterate this method up to the 10 th possible breakpoint. Thus, on average there are 7.3 breakpoints per profile (the first is always chosen and the 9 others chosen independently with probability 0.7). This is closer to what we see in reality, with aCGH profiles of patients with the same disease state sharing certain key breakpoints, yet but not all.

We then proceed as in (1) for 1000 trials. One minor detail: if the first few randomly generated profiles only exhibit a subset of size s < 10 of the 10 possible breakpoints and the group LARS algorithm finds these s breakpoints before any others, we stop the trial at this point, as we cannot expect all 10 breakpoints to be found if some of them have not yet been randomly exhibited.

4. All profiles have a subset of a predefined set of breakpoints though the exact location of each breakpoint can vary slightly between profiles. This experiment is a direct combination of ( 2) and (3). Again, 1000 trials were performed for each level of noise. The minor detail mentioned in (3) is treated in the same way here.

Simulation results from experiments 1-4 are shown in Figures 34. The main result is that, given enough profiles, the algorithm correctly selected the 10 breakpoint locations/regions for every experimental condition and every noise level. Specifically, in lower noise conditions (σ 2 = 0.01, 0.1 and 0.2), rarely are more than fifty profiles needed to correctly select the true breakpoints. In more realistic conditions (σ 2 = 0.5), with high probability, up to 75-200 profiles were necessary to correctly select the whole set of true breakpoints, though often much fewer were required.

Application to bladder tumor CGH profiles

We considered a publicly available aCGH data set of 57 bladder tumor samples [START_REF] Stransky | Regional copy number-independent deregulation of transcription in cancer[END_REF]. Each aCGH profile gave the relative quantity of DNA for 2215 probes. We removed the probes corresponding to sexual chromosomes, because the sex mismatch between some patients and the reference used made the computation of copy number less reliable, giving us a final list of 2143 probes. remark that even though (c) was forced to have the same breakpoints as (b), this does not translate into a poor smoothed version of (c), rather, the forced breakpoints are tiny jumps that can be ignored by biologists.

Fig. 5(a) confirms nearly all of the duplications and deletions associated with bladder cancer found in [START_REF] Kallioniemi | Comparative genomic hybridization for molecular cytogenetic analysis of solid tumors[END_REF][START_REF] Kallioniemi | Identification of gains and losses of DNA sequences in primary bladder cancer by comparative genomic hybridization[END_REF][START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF]: frequent duplication of 8q22-24, 17q21 and 20q is observed, and frequent deletion of 8p22-23, 13q, 17p, 11p and all of chromosome 9. The two known duplications that could not be confirmed here were 12q14-15 and 11q13. Fig. 5(a) suggests other potentially important CNVs, including frequent duplication of 1q, 5p and deletion of 4q and 10q.

Discussion

Segmentation of a single aCGH profile into regions of constant copy number, separated by breakpoints, is a well-studied problem [START_REF] Olshen | Circular binary segmentation for the analysis of array-based DNA copy number data[END_REF][START_REF] Hupé | array CGH data: from signal ratio to gain and loss of DNA regions[END_REF][START_REF] Picard | A statistical approach for array CGH data analysis[END_REF][START_REF] Fridlyand | Hidden markov models approach to the analysis of array CGH data[END_REF][START_REF] Wang | A method for calling gains and losses in array CGH data[END_REF][START_REF] Picard | A segmentation-clustering problem for the analysis of array CGH data[END_REF][START_REF] Huang | Robust smooth segmentation approach for array cgh data analysis[END_REF][START_REF] Tibshirani | Spatial smoothing and hot spot detection for CGH data using the fused lasso[END_REF]. Recently, attempts have been made to deal simultaneously with many profiles [START_REF] Diskin | STAC: a method for testing the significance of DNA copy number aberrations across multiple array-CGH experiments[END_REF][START_REF] Rouveirol | Computation of recurrent minimal genomic alterations from array-CGH data[END_REF][START_REF] Klijn | Identification of cancer genes using a statistical framework for multiexperiment analysis of nondiscretized array CGH data[END_REF][START_REF] Picard | Joint segmentation of multivariate Gaussian processes using mixed linear models[END_REF][START_REF] Robin | Simultaneous occurrences of runs in independent Markov chains[END_REF]. In these methods, the search for important shared CNV regions tends to occur as the final step, either by choice of a level of significance [START_REF] Diskin | STAC: a method for testing the significance of DNA copy number aberrations across multiple array-CGH experiments[END_REF] or in post-processing [START_REF] Rouveirol | Computation of recurrent minimal genomic alterations from array-CGH data[END_REF]. These methods do not use the biological prior information that different profiles are likely to share at least some breakpoints. Also, they may require the choice of pertinent kernel windows [START_REF] Klijn | Identification of cancer genes using a statistical framework for multiexperiment analysis of nondiscretized array CGH data[END_REF], lose information through discretization [START_REF] Diskin | STAC: a method for testing the significance of DNA copy number aberrations across multiple array-CGH experiments[END_REF][START_REF] Rouveirol | Computation of recurrent minimal genomic alterations from array-CGH data[END_REF] or involve prohibitive computational complexity for large p [START_REF] Picard | Joint segmentation of multivariate Gaussian processes using mixed linear models[END_REF].

To our knowledge, we have introduced for the first time a way to explicitly code the prior biological information of expecting patients with the same disease to share certain CNVs. Our method forces breakpoints to be located in the same places for all profiles. This has the effect of selecting breakpoint locations where many, but not necessarily all, profiles exhibit a breakpoint. This corresponds exactly to one of the underlying biological goals in CNV studies. As shown in Fig. 5(c), it is important to note that a profile forced to have breakpoints where it clearly does not, still ends up with a good quality smoothed representation. We also showed that superimposing all smoothed versions on one graph allows intuitive visual interpretation of the data in a lower-dimensional form. On a real bladder cancer data set [START_REF] Stransky | Regional copy number-independent deregulation of transcription in cancer[END_REF], the smoothed versions we obtained (Fig. 5(a)) confirmed nearly all of the known CNVs described in the articles [START_REF] Kallioniemi | Comparative genomic hybridization for molecular cytogenetic analysis of solid tumors[END_REF][START_REF] Kallioniemi | Identification of gains and losses of DNA sequences in primary bladder cancer by comparative genomic hybridization[END_REF][START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF]. Furthermore, our proposed algorithm is extremely fast. Even at the limit of current aCGH technology, it is practical, taking a few minutes on a single laptop computer.

The piecewise-constant versions of the original profiles can be seen as extracted low-dimensional features. These can potentially be used to implement classification algorithms to discriminate between two or more classes of aCGH profile, e.g., different disease states. For example, each piecewise-constant profile could simply be treated as a vector of the constant values. As the breakpoints are forced to be in the same place on all profiles, the vector representation is the same size for each profile, directly opening the way for the use of many well-known classification methods. This is a promising research direction.

The question of how many breakpoints to choose, i.e., when to 'stop' the algorithm, remains open. There are at least two possible solutions. First, if the algorithm were to be associated with a classification algorithm, a stopping criteria using internal crossvalidation on the learning set can be defined. Second, it might be useful to calculate how much each subsequently selected breakpoint closes the distance between the set of smoothed and original profiles, and define a stopping criteria based on this.

Figure 1 :

 1 Figure 1: Speed trials. (a) CPU time for finding 50 breakpoints when there are 2000 probes and the number of profiles varies from 1 to 20. (b) CPU time when finding 50 breakpoints with the number of profiles fixed at 20 and the number of probes varying from 1000 to 10000 in intervals of 1000. (c) CPU time for 20 profiles and 2000 probes when selecting from 1 to 50 breakpoints.

Figure 2 :

 2 Figure 2: Simulating an aCGH profile with different levels of white noise. (a) is the underlying piecewise-constant profile. White noise from a N (0, σ 2 ) is added to the probes of (a) with: (b) σ 2 = 0.01, (c) σ 2 = 0.1, (d) σ 2 = 0.2 and (e) σ 2 = 0.5.

Fig. 5 (

 5 a) shows the result of superimposing the smoothed versions of the 57 bladder tumor aCGH profiles, when the algorithm has selected 80 ranked common breakpoints. Figs 5(b) and (c) show 2 of the original 57 profiles and their associated smoothed version, where (b) was a profile exhibiting much instability, and (c) only on chromosome 9. We

Figure 3 :

 3 Figure 3: Simulation conditions 1 (left) and 2 (right). Left: histograms of the number of profiles required to correctly predict 10 real breakpoints with no mistakes in the presence of white noise. Right: histograms of the number of profiles required to correctly predict all real breakpoints when each profile exhibits a breakpoint in each of 10 tightly defined regions, in the presence of white noise. The noise is N (0, σ 2 ) with (a) σ 2 = 0.01, (b) σ 2 = 0.1, (c) σ 2 = 0.2 and (d) σ 2 = 0.5. Each experiment was performed 1000 times.

Figure 4 :

 4 Figure 4: Simulation conditions 3 (left) and 4 (right). Left: histograms of the number of profiles required to correctly predict all real breakpoints when each profile exhibits a subset of a predefined set of 10 breakpoints, in the presence of white noise. Right: histograms of the number of profiles required to correctly predict all real breakpoints when each profile exhibits a breakpoint in a subset of a predefined set of 10 tightly defined regions, in the presence of white noise. The noise is N (0, σ 2 ) with (a) σ 2 = 0.01, (b) σ 2 = 0.1, (c) σ 2 = 0.2 and (d) σ 2 = 0.5. Each experiment was performed 1000 times.

Figure 5 :

 5 Figure 5: Graphical representation. (a) superimposition of the smoothed versions of 57 bladder tumor aCGH profiles [21] with 80 breakpoints. Vertical lines divide chromosomes 1-22. (b) a profile exhibiting many CNVs, and its smoothed version. (c) a profile only showing a deletion on chromosome 9, and its smoothed version. Smoothed profiles are obtained by replacing the set of probe values between consecutive breakpoints with their mean value.

All profiles share the same breakpoint 'regions', though the breakpoints are not all


W Frelin

L Berthet
email: luc.berthet@edf.fr

Y Brument
email: yves.brument@edf.fr

M Petit
email: marc.petit@supelec.fr

G Perujo

J C Vannier

THERMAL BEHAVIOR OF LV CABLES IN PRESENCE OF HARMONIC CURRENTS

The proliferation of non-linear loads on public and industrial networks causes increases of harmonic current distortion. This paper deals with the thermal impact of harmonic distortion on LV cables, especially LV cables with neutral conductor. The thermal impact of harmonics is shown by measurements and FEM simulations on two different LV cable structures.

Introduction

Nonlinear loads (computers, compact fluorescent lamps, variable speed drives...) produce harmonic currents and therefore increase the harmonic distortion of the currents flowing through power distribution systems. The RMS value of the currents in phase or neutral conductors increases, and thus the I 2 R losses rise. This increase of losses creates additional heating in the cable [START_REF] Desmet | Thermal analysis of the influence of nonlinear, unbalanced and asymmetric loads on current conducting capacity of LV-cables[END_REF]- [START_REF] Desmet | Modelling and sensitivity analysis of the Thermal behaviour of LV-cables for different current conditions[END_REF], which may have a critical impact on the cable life expectancy [START_REF] Pachecho | Power Quality Impact on Thermal Behaviour and Life Expectancy of Insulated Cables[END_REF]. In power distribution systems, LV cables supplying nonlinear loads are impacted by harmonics, in particular the neutral conductor [START_REF] Thomas | A Survey of Neutral Currents in Three-Phase Computer Power Systems[END_REF]. When the harmonic content of phase currents includes triplen harmonics, the neutral current is non-zero. This current circulating in the neutral conductor can rise above phase current and bring to light a dimensioning problem of the neutral conductor. Previous works on a 50 sqmm LV-cable, computed by FEM methods, presented in [START_REF] Frelin | Evaluation of thermal heating in LV-cables in presence of harmonics[END_REF] show the cable temperature increase in presence of harmonics. This paper focuses on the thermal behavior of LV cables with neutral conductor, as shown in Fig. 1, when supplying non linear loads. These cables are composed of three 150 sqmm phase conductors and one neutral conductor with different sections. The two structures studied in this paper are made as follows: (i) a 90 sqmm solid neutral conductor (left scheme) and (ii) a 70 sqmm stranded neutral conductor with a lead sheath (right scheme). The first configuration is the one that is chosen for the new cable generation. After introducing the configuration of the experimental system, FEM modeling of LV cables is presented. Next, the thermal behavior of the two LV cables given by FEM simulation is presented and compared to measurement results.

Configuration of experimental system

The experimental system, shown in Fig. 2 and Fig. 3, has three parts: (1) a 230V 3 phase power supply is feeding a 40kVA transformer that supplies a five meter length LV cable and three 1-phase rectifiers with resistive loads or resistive and capacitive loads. Current harmonics are produced by the 1-phase rectifiers with RC loads, and in this configuration triplen harmonic currents circulate in the neutral conductor. (2) A thermal acquisition system with twelve thermocouples on each LV cable, which gives sheath and core temperatures. Fig. 1 gives position of sheath and core thermocouples, sheath thermocouples are pressed against conductors using special glue which aids heat transfer and core thermocouples are positioned on neutral and phase conductor cores without glue. (3) An acquisition system which measures phase and neutral currents. This acquisition system computes current amplitudes for each harmonic order (up to harmonic order 50) and records each current waveform. Before measurement, thermocouples and current acquisition system were gauged in laboratory, with accuracy materials.

FEM Modeling

Finite Element Method (FEM) is used to solve the heat transfer equation coupled to an electrical model of the conductors. Stranded conductors (the phases of the two LV cables considered and the neutral of one of them) are considered for simulation as solid conductor with the same section. The electrical current for each conductor in linear and non-linear conditions is described in table I. Let us note that in a symmetric and balanced network, the neutral current without harmonics is equal to zero.

Table I Electrical equations

Linear conditions [START_REF] Desmet | Thermal analysis of the influence of nonlinear, unbalanced and asymmetric loads on current conducting capacity of LV-cables[END_REF] Non linear conditions (2) 3), which results from a simplification of Maxwell's equations, allows the calculation of the current density repartition in each conductor. By using this equation, skin and proximity effects are taken into account. The parameter V is adjusted for each phase conductor to have the same current. For neutral conductor, V is adjusted to fix neutral current to its value.

( ) ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = = 0 3 2 cos 2 3 2 cos 2 cos 2 1 1 1 N C B A I t I I t I I t I I π ω π ω ω ( ) ( ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ + ⋅ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = + = ∑ ∑ ∑ ∑ ∞ = ∞ = ∞ = ∞ = 1 3 3 1 1 1 3 cos 2 3 3 2 cos 2 3 2 cos 2 cos 2 h h h N h h h C h h h B h h h A t h I I t h I I t h I I t h I I ϕ ω ϕ π ω ϕ π ω ϕ ω Equation (
( )

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ - - = - = ∧ ∇ ∧ ∇ + z r z z r A i dz dV j L V A A i ωσ σ σ μ μ μ ωσμ 0 0 (3)
With μ 0 the permeability of free space, μ r the relative permeability, ω the pulsation, σ the electrical conductivity, A z the vector potential, V the electric potential difference along the cable, L the cable length, j the current density and i the complex number. The aim of electromagnetic calculation is to estimate Joule losses in cable, which are in fact the heat quantity Q in thermal calculation. Equation (4) explains Joule losses in cable in linear and non linear conditions respectively.

⎪ ⎩ ⎪ ⎨ ⎧ + + + = + + = ∑ ∞ = 2 1 2 3 3 2 2 2 2 1 1 2 3 3 2 2 2 2 1 1 Nh Nh h h h h h h h j j I R I R I R I R P I R I R I R P (4)
Where R 1 , R 2 , R 3 are the resistance (phase 1, phase 2, phase 3 respectively) at fundamental frequency. I 1 , I 2 , I 3 are the fundamental current (phase 1, phase 2, phase 3 respectively). R 1h , R 2h , R 3h, R Nh are the resistance for harmonic h (phase 1, phase 2, phase 3 and neutral respectively). I 1h , I 2h , I 3h , I Nh are the current for harmonic h (phase 1, phase 2, phase 3 and neutral respectively). Heat transfer in the inner of LV cable is dominated by conduction, thus the temperature repartition in each domain may be calculated from the Poisson equation ( 5).

( )

Q T = ∇ ⋅ ∇ - λ (5)
With λ the thermal conductivity, T the temperature and Q the heat quantity.

On the sheath surface the heat is transferred to the environment by natural convection and heat radiation. The heat flux φ through the sheath surface S is given by the Newton law:

( ) T S h h r c Δ + = ϕ (6)
With h c the convection coefficient for natural convection and h r an equivalent convection coefficient for heat radiation. For natural convection, the convection coefficient h c can be expressed from

D Nu h c λ = (7)
With λ the air conduction coefficient, D the LV cable diameter and Nu the Nusselt number which is calculated for natural convection by [START_REF] Bianchi | Transfert Thermique[END_REF] ( ) With Ra the Rayleigh number, g the gravitational acceleration, β the expansion coefficient of air, v the kinematic viscosity of air and a the thermal diffusivity. The heat flux φ r dissipated by heat radiation between two surfaces S 1 and S 2 , respectively at the temperature T 1 and T 2, can be calculated by [START_REF] Jannot | Transfert thermique[END_REF]. Considering that S 2 >>S 1 , the heat flux expression given by [START_REF] Jannot | Transfert thermique[END_REF] can be simplified and expressed by equation (9), with S 1 the sheath surface (S 1 =S).

( )

4 2 4 1 T T S b r - = ε σ ϕ (9)
With σ b the Boltzmann constant, ε the thermal emissivity, T 1 the sheath temperature and T 2 the ambient temperature (T 1 and T 2 in Kelvin). Introducing the equivalent convection coefficient h r , equation (9) can be rewritten as follows

( )( ) 2 1 2 2 2 1 T T T T h with T S h b r r r + + = Δ = ε σ ϕ (10)
The global convection coefficient is the sum of the convection coefficient given by natural convection [START_REF] Jannot | Transfert thermique[END_REF] and the equivalent convection coefficient given by heat radiation (10). For an ambient temperature of 23°C and a thermal emissivity of 0.9, the global convection coefficient value is spread from 12.7 to 14. Electrical conductivity is temperature dependent according to:

( ) (
)

20 1 20 - + = T T α σ σ (11)
With σ 20 the electric conductivity at 20°C, α the temperature coefficient of the conductor material, T the conductor temperature, σ (T) the electrical conductivity at the temperature T. Fig. 4 shows the simulation process to obtain the LV cable temperature. The electromagnetic calculation gives the Joule losses in the cable which represent the heat quantity in thermal module. In a first time, the electromagnetic calculation is done with an initial temperature in the conductor core of 20°C and a corresponding electrical conductivity. Then thermal calculation is launched and a new temperature profile is known. As electrical conductivity depends on temperature, a new electromagnetic calculation is required with an update of the electrical conductivity at the temperature calculated previously. Simulations are computed until a constant temperature is obtained.

Fig. 4 Simulation process

Table II gives the electrical and thermal properties of the LV cable materials.

Table II Electrical and thermal properties of LV cable

Core Insulator Filler Screen Sheath λ (W/K.m) 160 0.25 0.16 26 0.2 σ (S.m -1 ) 3.10 7 10 -7 10 -15 1.8.10 6 10 -15

Results

In order to investigate thermal behaviour of the cable under linear and non linear supply conditions, different simulations and measurements were carried out with a constant RMS current, fixed to 210A for each case. Table III summarizes the harmonic spectrum for the three different cases, the first case corresponding to linear conditions and the other two cases to non-linear conditions with different harmonic contents. Fig. 5 and Fig. 6 summarize results given by simulations and measurements for cases described in table III and the two different LV cable structures. Fig. 5 gives core temperatures on phase 2 and neutral conductors obtained by simulations and measurements for the different cases. Fig. 6 gives sheath temperatures on phase 2 and neutral conductors obtained by simulations and measurements for the different cases. Each measurement result is composed of two or three values (due to the thermocouple number) and presents a small gap due to the difficulty to lay properly thermocouples. These figures show a good agreement between FEM simulations and measurements. It can be noticed that without harmonics, the thermal behaviour of the two LV cables is the same because there is no current in neutral conductor and phase conductors are the same. Then these figures show that neutral and phase conductor temperatures are higher with harmonic currents because the neutral conductor acts as a fourth thermal source. The simulated temperature distribution given by case 3, for the two cable structures, is presented in Fig. 7, cable with solid neutral conductor on the left and cable with stranded neutral conductor on the right. This figure shows that LV cable with solid neutral conductor is less sensitive to harmonics than LV cable with stranded neutral conductor. As the solid neutral conductor cross section is higher, compared to the stranded neutral conductor cross section, Joule losses and temperatures are lower with the same harmonic content. Fig. 7 Cable temperature of LV cable with solid neutral conductor (left) and LV cable with stranded neutral conductor (right) for case 3.

The first conclusion of the different simulations and measurements presented previously is that harmonics increase temperature in LV cable. Another point of view is to consider that harmonics reduce the transport capacity, i.e. the maximum current supported by the cable.

In order to show the reduced level of current that can circulate in LV cables in the presence of harmonics, the current which gives a 90°C temperature in one of the four conductors is calculated. Both cable structures are considered and computed in underground conditions because this is the main use of these LV cables. Fig. 8 gives the maximum fundamental current value which can circulate in both LV cables for each case described in table III and shows the difference of transport capacity with the two neutral conductor structures. It can be seen that harmonics reduce this fundamental current, in other words if harmonic content increases, the amount of loads connected to the LV cable has to be reduced.

Conclusion

Simulations and measurements presented in this paper show temperature increases on a 150 sqmm LV cable in presence of harmonics. The circulation of triplen harmonic currents in the neutral conductor, which in this case acts as a fourth heating source, contributes to the thermal heating of the LV cable. It can be seen that the LV cable with solid neutral conductor is less sensitive to harmonics due to a smaller neutral conductor resistance. This paper shows a good agreement between FEM simulations and measurements. At last, it is shown in simulations that the transport capacity is reduced when harmonic distortion is higher.

Fig. 1 .Fig. 2 .

 12 Fig. 1 . Cable geometry and thermocouple position

  ISEF 2009 -XIV International Symposium on Electromagnetic Fields in Mechatronics, Electrical and Electronic Engineering Arras, France, September 10-12, 2009

Fig. 5 .

 5 Fig.5. Phase 2 and neutral core temperature for the two different LV cable structures given by simulations and measurements (LV cable with solid neutral conductor on the left and LV cable with stranded neutral conductor on the right)

Fig. 6 .

 6 Fig.6. Phase 2 and neutral sheath temperature for the two different LV cable structures given by simulations and measurements (LV cable with solid neutral conductor on the left and LV cable with stranded neutral conductor on the right)

ISEF 2009 -

 2009 XIV International Symposium on Electromagnetic Fields in Mechatronics, Electrical and Electronic Engineering Arras, France, September 10-12, 2009

  Fig. 8 Influence of harmonics on the transport capacity

Arras, France, September 10-12, 2009

Acknowledgment

The authors wish to thank Department LME of EDF R&D, especially Y.Brument and R. Fonder, for their valuable material and scientific support during the test. The authors also wish to acknowledge Y. Jannot for his valuable discussions and suggestions.


W Frelin
email: frelin@crans.org

L Berthet
email: luc.berthet@edf.fr

M Petit
email: marc.petit@supelec.fr

J C Vannier

Transformer Winding Losses Evaluation when Supplying Non Linear Load

Keywords: Eddy-current, harmonics, Transformer, winding losses

Nonlinear loads produce harmonic currents, which induce additional losses in transformers and cause temperature rise, especially in the windings. The estimation of winding eddycurrent loss, in the presence of harmonics, is based on the knowledge of winding eddy-current loss at power frequency. This paper deals with the transformer winding loss estimation and gives a simplified winding resistance expression that takes winding eddy-current into account and allows to estimate associated losses. The impact of transformer winding characteristics, as conductor size and layer number, on the resistance variation with frequency is shown.

I. INTRODUCTION

The growing use of power electronics devices in industry, office and home equipment has led to study the impact of harmonic distortion on power networks materials. Current and voltage harmonics may cause a large number of problems for electrical equipment:

-Malfunction of circuit breakers and electronic equipment [START_REF] Subjak | harmonics-causes, effects, measurement, and analysis: An update[END_REF]. -Additional losses in motors, transformers and conductors [START_REF] Tofoli | Analysis of losses in cables and transformers under quality related issues[END_REF], [START_REF] Lin | Harmonic losses of single phase induction motors under nonsinusoidal voltages[END_REF]. -Additional heating due to these additional losses. IEEE Std C57.110-1998 [START_REF]IEEE Recommended Practice for Establishing Transformer Capability When Supplying Nonsinusoidal Load Currents[END_REF] was created to establish uniform methods for determining the capability of transformers to supply non linear loads. Transformer losses are separated into no load loss (core loss) and load loss [START_REF]IEEE Recommended Practice for Establishing Transformer Capability When Supplying Nonsinusoidal Load Currents[END_REF]. Load loss (P ll ) is subdivided into I²R loss and stray loss due to stray electromagnetic flux in the windings, core, core clamps, magnetic shields, enclosure or tank walls, etc. The stray loss is subdivided into winding eddy-current loss (P ec ) and stray loss in components other than the windings (P osl ) as explained in [START_REF] Subjak | harmonics-causes, effects, measurement, and analysis: An update[END_REF]. Where P ll is load loss. P j =R dc .I ² is loss due to load current and dc resistance of the windings.

P ec is winding eddy-current loss.

P osl is stray loss in components other than the windings. R ac is AC winding resistance. R dc is DC winding resistance.

I is the RMS current

When a transformer supplies non linear loads, I 2 R increases with the RMS current variation and winding eddy-current loss increases with the square of the load current and the square of the frequency. Equation ( 2) allows calculating winding eddycurrent loss through a proportionality factor applied to the winding eddy-current loss at fundamental frequency [START_REF]IEEE Recommended Practice for Establishing Transformer Capability When Supplying Nonsinusoidal Load Currents[END_REF]. The reliance on the square of the frequency is recognized in this expression with the h² term. This expression shows the necessity to know fundamental eddy-current loss to determine the harmonic impact of eddy-current loss on the transformer.

hl o EC EC F P P ⋅ = - where 2 1 1 2 2 1 1 max max ∑ ∑ = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = h h h h h h hl I I h I I F (2)
Where P ec is winding eddy-current loss. P ec-o is winding eddy-current loss at the measured current and the power frequency.

F hl is harmonic loss factor for winding eddy-current. h is harmonic order. I h is the RMS current at harmonic frequency of order h. I 1 is the RMS current at fundamental frequency.

The aim of this paper is to give a winding resistance expression, which takes winding eddy-current into account and allows estimating P EC-o and then winding loss in presence of harmonics. This analytical expression shows the impact of transformer winding characteristics, as conductor size and transformer winding geometry, on the resistance frequency dependence. After theoretical results on a simple conducting plate, the method to obtain the analytical expression of winding resistance is developed. Finally, some numerical results on two different transformers are presented to show the impact of winding characteristics on the AC winding resistance.

Fig. 1 shows a simple conducting plate, assumed to be infinite in the y and z directions and carrying a sinusoidal current density J in the y direction, with sinusoidal magnetic fields of constant amplitudes h 1 and h 2 outside the plate in the z direction. The distribution of magnetic field H in metals at relatively low frequencies is described by the diffusion equation ( 3), which is calculated from Maxwell's equations.

t H H

∂ ∂ μσ = ∇ 2 (3)
Where μ is the magnetic permeability and σ is the electrical conductivity.

Considering that the magnetic field is only with z component and as there are no variations in the y and z directions, the magnetic field equation can be written:

t H x H z z ∂ ∂ μσ ∂ ∂ = 2 2 (4) A solution of this equation is of the form H z =Re(h(x)e jωt ), thus equation (4) becomes: 0 ) ( ) ( 2 2 2 = + x h k dx x h d (5) With ωμσ j k - = 2
Where ω is the pulsation and j the complex number defined by j² = -1.

The equation ( 5) solution is

x k j x k j e a e a x h - + = 2 1
) ( [START_REF] Vandelac | A novel Approach for Minimising High-Frequency Transformer Copper Losses[END_REF] The integration constants a 1 , a 2 are determined from the known boundary conditions.

⎪ ⎩ ⎪ ⎨ ⎧ + = = - + = = - - jka jka jka jka e a e a h a h e a e a h a h 2 1 1 2 1 2 ) ( ) (
Hence, the magnetic field expression is

( ) [ ] ( ) [ ] ( ) ka j a x jk h a x jk h x h 2 sinh sinh . sinh . ) ( 2 1 + + - - = (7)
The distribution of current density J is calculated from Maxwell's equations.

x

H J z y ∂ ∂ - = (8)
A solution of this equation is of the form J y =Re(j(x)e jωt ), thus equation ( 8) becomes :

dx x dh x j ) ( ) ( - = (9)
Hence, the current density expression is

( ) [ ] ( ) [ ] ( ) ⎭ ⎬ ⎫ ⎩ ⎨ ⎧ + - - = ka j a x jk h a x jk h jk x j 2 sinh cosh . cosh . ) ( 2 1 (10)
The loss per square meter is obtained by integrating ρ.J 2 .

∫ =

+ - a a dx x j P ds 2 ) ( 2 1 ρ (11)
Consequently, the Joule losses in the plate may be determined by the integral (11) and the current density expression (10) as follows

⎭ ⎬ ⎫ ⎩ ⎨ ⎧ - + = ) ( 2 2 ) ( 2 2 2 2 1 1 2 2 2 1 ϕ ϕ σ ωμσ f h h f h h P ds (12) With δ ϕ e = and ωμσ δ 2 = ) 2 cos( ) 2 cosh( ) 2 sin( ) 2 sinh( ) ( 1 ϕ ϕ ϕ ϕ ϕ - + = f (13) ) 2 cos( ) 2 cosh( ) cos( ). sinh( ) sin( ). cosh( ) ( 2 ϕ ϕ ϕ ϕ ϕ ϕ ϕ - + = f (14)
Where e is the thickness of the plate and δ is the depth of penetration.

This expression can be rewritten as shown in equation (15).

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ - + = ) ( . 4 ) ( . ) 1 ( 2 2 2 1 1 2 2 2 1 2 2 ϕ ϕ σδ f h h f h h h P ds (15)
For conductors with small thickness compared to the depth of penetration (φ<<1), functions f 1 (φ) and f 2 (φ) can be approximated by g 1 (φ) and g 2 (φ) given in ( 16) and ( 17) respectively.

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⋅ + = ≈ ) ( 45 4 1 1 ) ( ) ( 8 4 1 1 ϕ ε ϕ ϕ ϕ ϕ g f (16) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + ⋅ - = ≈ ) ( 90 7 1 2 1 ) ( ) ( 8 4 2 2 ϕ ε ϕ ϕ ϕ ϕ g f (17)
Fig. 2 shows comparison, between functions f 1 (φ), f 2 (φ) and g 1 (φ), g 2 (φ) for different values of the ratio e/δ. Fig. 3 shows the relative error in per cent between the functions and their approximations.

Fig. 3. Error in per cent between functions f 1 (φ)& g 1 (φ) and f 2 (φ)& g 2 (φ)

versus ratio e/δ.

It can be seen that the two functions g 1 (φ) and g 2 (φ) give good results when the ratio between thickness and depth of penetration is less than 1.4 (errors less than 5%). Thus, the Joule losses in the plate can be expressed according to (18) substituting f 1 (φ) and f 2 (φ) in equation ( 15) by g 1 (φ) and g 2 (φ). For analysis, the structure of transformer windings can be approximated by a superposition of several conducting plates [START_REF] Dowell | Effects of eddy currents in transformer windings[END_REF]. Fig. 4 shows the geometry of MV/LV transformer with the LV and MV windings respectively in the inner and in the outer. The curvature of the conductors, edge and end effects are neglected.

⎪ ⎭ ⎪ ⎬ ⎫ ⎪ ⎩ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - = 45
Using ampere's law, magnetic fields on both the left-hand side h 1 and right-hand side h 2 are respectively given by ( 19) and ( 20) for layer m.

d i n N m h ) 1 ( 1 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ - = (19) d i n N m h 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ = (20)
With m the layer number from the inner to the outer. N the total number of turns in the windings. n the total number of layers in the windings. î the peak current in a turn. d the height of windings. Thus, the power per unit length, according to the y direction, dissipated in layer m (21) can be calculated using Joule losses expression in a plate (18) multiplied by the height of windings and magnetic fields expressions given by ( 19) and (20).

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⋅ ⋅ ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = 45 4 3 3 1 2 4 2 2 m m e d e I n N P m δ σ (21)
Total loss per unit length in transformer windings is then given by the sum of power dissipated in each layer (22). Introducing the DC resistance per unit length of the conductor (23), the AC resistance per unit length can be identified in the expression of total losses per unit length in transformer winding (22). Eddy-current losses are represented by the second term into brackets. Transformer geometry impact can be seen by the presence of layer number -n -and conductor thickness -e -in the expression. IV. NUMERICAL APPLICATION Generally, LV windings are made with rectangular or sheet conductors, while MV windings are made with small round conductors. Table I gives the winding characteristics of two typical 100 kVA transformers which are used for numerical application presented next. The ratio e/δ, given in the last column, is calculated for secondary conductors at 2.5kHz.

∑ = ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⋅ ⋅ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = =
⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⋅ = n N de n R dc σ (23) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = ⋅ =
As primary winding characteristics are almost the same, winding ratio evolution is similar. The main difference between the transformers is the secondary winding geometry (Table I).

Then we present in Fig. 5 and Fig. 6 the R ac /R dc ratio evolution as a function of frequency for the secondary windings of transformer n°1 and n°2 respectively. The ratio is either calculated with (24) using g 1 and g 2 , or with f 1 and f 2 . Theses figures show that winding resistance with sheet conductors is less sensitive to frequency than winding resistance with rectangular conductors. Thus when a transformer is supplying non linear loads, winding eddycurrent loss increase is much higher for rectangular conductors. Moreover (24) calculated with the simplified functions g 1 and g 2 gives a conservative result for transformer n°1 as soon as e/δ is greater than 1.4 (i.e. frequency greater than 800 Hz in the present case).

In presence of harmonics, eddy-current loss P EC can be calculated using resistance expression (24).

∑ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =
⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = - δ (26)
Where δ h is the depth of penetration for harmonic h. δ R is the depth of penetration at the power frequency. Substituting (26) in (25), eddy-current loss expression can be rewritten and gives the same result that the IEEE Std C57.110-1998 (2). Rac/Rdc ratio Transformer n°2 with with g 1 and g 2 functions Fig. 6. Secondary winding R ac /R dc ratio evolution as a function of frequency for transformers n° 2 presented in table I, and using functions {g 1 , g 2 }equation ( 24) .

The assumption that eddy-current losses in transformer windings are proportional to the square of the frequency (25) gives accurate results for small conductors or at low frequency (e/δ < 1.4). For large conductors and at high frequency (e/δ > 1.4), it leads to conservative results [START_REF] Sergey | Corrected Harmonic Loss Factor for Transformers supplying Nonsinusoidal Load currents[END_REF][START_REF] Elmoudi | Corrected Winding Eddycurrent Harmonic Loss Factor for transformers Subject to Nonsinusoidal[END_REF] as shown in Fig. 5.

V. CONCLUSIONS

The F hl coefficient [START_REF]IEEE Recommended Practice for Establishing Transformer Capability When Supplying Nonsinusoidal Load Currents[END_REF] is generally used to determine additional losses and capability of transformer when non linear loads are present, but it is necessary to know fundamental eddy-current loss to determine the harmonic impact on loss increase in transformer.

Fig. 1 .

 1 Fig. 1. Studied conducting plate.

Fig. 2 .

 2 Fig.2. Functions f 1 (φ), f 2 (φ), g 1 (φ) and g 2 (φ ) versus ratio e/δ.

Fig. 4 .

 4 Fig. 4. Geometry of MV/LV transformer.

2

 2

  loss at the measured current I and power frequency, P EC-O

Fig. 5 .

 5 Fig.5. Secondary winding R ac /R dc ratio evolution as a function of frequency for transformers n°1 presented in table I, and using functions {g 1 , g 2 }equation (24) -or {f 1 , f 2 }.

This paper gives a winding resistance expression which takes winding eddy-current into account and permits to estimate winding losses at power frequency, P EC-0 , and also in presence of harmonics. This expression shows the impact of the following winding parameters: conductor thickness and layer number. This expression also underline that transformers composed of thin foil windings have winding eddy-current losses smaller than transformer with rectangular conductors. Thus, transformers with thin foil windings are less sensitive to harmonics. This work also shows that the F hl coefficient gives conservative results for large conductors.