Zhe Yuan
email: zhe.yuan@centralesupelec.fr

Yale Gong
email: gong@em-lyon.com

Improving the Speed Delivery for Robotic Warehouses

Keywords: Mobile Fulfilling System, Speed Delivery, Order Picking, Robotic warehouse, Warehouse Management, Optimization of facility layout, Open Queue theory

This paper studied a new order fulfilling system using mobile robots instead of manual picking at online retailers' distribution centers. This new generation of sustainable green warehouse systems can improve productivity and flexibility. We measured the performance of system, and provided design rules for velocity of robots. We builded open queue models for the new order fulfilling system and calculated throughput time of this system given the number of robots. This is one of earliest papers to introduce this new material handling system.

INTRODUCTION

Recently online retailer Zappos had installed a new Kiva Mobile Fulfilling System (MFS) at its distribution center in Shepherdsville, Kentucky in 2008. This new system significantly makes contribution to developing a supply chain to speed delivery [START_REF] Marks | Zappos. com: Developing a supply chain to deliver WOW[END_REF]. Robotic warehouse techniques are researched before, but until 2010, Kiva MFS has been applied to warehouses like Zappos, Walgreens, Staples, Diapers and Gap. Business Week report "Kiva Systems are revolutionizing companies' supply-chain and distribution networks" [START_REF] Scanlon | How kiva robots help zappos and walgreens[END_REF]. In 2012, Amazon bought Kiva system for 775 million USD and implemented a number of Kiva MFS. In August 2015, the company officially changed its name from Kiva Systems to Amazon Robotics [START_REF] Nussbaum | Changing the tooth-to-tail ratio using robotics and automation to beat sequestration[END_REF]. Some other material handling companies and robot companies are developing different robotic picking systems. For example, RMT Robotics, a high-tech company in Toronto, Canada is developing robotic picking systems and robotic case picking systems (see RMTrobotics.com). ADAM (autonomous delivery and manipulation) robots of RMT companies are similar to Kiva's robots. Jingdong (see jd.com), one of the largest online retailers with revenue of 18.51 billion USD and sales of 41.9 billion USD in 2014, has built a MFS and implementing robotic picking in a warehouse of 5,000 m 2 in Shanghai.

Kiva MFS is a "parts-to-picker" warehouse systems. Compared with other order fulfillment systems, Kiva MFS can boost flexibility since it can adapt to changes of product types which is not easy in traditional "AS/DS" equipments. Another advantage of Kiva MFS is in scalability: While traditional materials handling technology need large initial investment on equipments like Mini-load AS/DS (about 750,000 USD per aisle, 3-10 aisles are common in a warehouse) and carousels (about 125,000 USD per This research is supported by China Scholarship Council. carousel, 3-10 carousels are common in a warehouse), Kiva MFS only needs to buy appropriate number of shelves and robots (about 5,000 USD per robot) that fit to their current demand volume [START_REF] Kopytoff | In Warehouses, Kiva Robots Do the Heavy Lifting[END_REF]. When demand scale increases, Kiva MFS simply needs to add robots and shelves incrementally. Kiva MFS is regarded as a new generation of sustainable green warehouse systems since it can cut carbon footprint (traditional systems need to keep on in most time whether they are needed or not. Kiva MFS robots only run when they are needed), has less noise (silent robots compared with noisy fork trucks and conveyor belts ), reduce injuries to warehouse employees since it is more economically designed (in a traditional warehouse, workers need to pick, lift,carry by hands or fork trucks, which bring hazards. Now workers only need to monitor stations), reduce lighting since robots can work with almost no lighting.

Our research is related to two research topics of AGV (Automated guided vehicle system): the fixed guided paths design and vehicle application in warehousing. In guided paths design with a fixed route, an AGV simply moves along the path determined in advance [START_REF] Ventura | A study of the tandem loop with multiple vehicles configuration for automated guided vehicle systems[END_REF], with the objective to minimize transport time including AGV's travel time and waiting time. However, other AGVs use a non-oriented path guided system [START_REF] Singh | Intelligent agent framework to determine the optimal conflict-free path for an automated guided vehicles system[END_REF], prefer softwares to hardwares in path management, and can relatively easily change their paths [START_REF] Tompkins | The vehicle routing problem[END_REF]. Our research is closer to unfixed guided paths. The selection of vehicle routing and scheduling influences the performance of AGV systems (Toth and Vigo, 2001). One of our objectives of the routing is to minimize throughput time of the warehouse. [START_REF] Hu | A framework for the selection of idle vehicle home locations in an automated guided vehicle system[END_REF] studied different scheduling strategies for idle vehicles discussing existing idle AGVs when service requests arrive. [START_REF] Lee | A realtime traffic control scheme of multiple agv systems for collision free minimum time motion: a routing table approach[END_REF] proposed the algorithm of sequential path generation to heuristic can generate collision-free shortest path for multiple AGVs. Qiu et al.

MODEL AND ANALYSIS

We use the following notations: R: the number of robots.

N : the number of pickers.

v: the velocity of robots.

w: the width of aisles.

W : the width of the warehouse.

L: the length of aisles.

T r: throughput time.

We consider two protocols: (1) R robots only serve for one order picker, (2) serve for multiple order pickers.

Pickers do not share robots

In this section, we consider the situation where R robots serve for 1 order picker. After receiving an order request, a robot follows the protocol:(1) leave order picker station, (2) go to storage grid to pick a shelf which includes requested item , (3) return to a picker station it is serving, (4) after the worker pick the item the robot will return the shelf to its original storage position and the robot will also return the picker station.

The system is a centralized system and the order arrival process to robots are Poisson processes. We assume the storage position follows uniform distribution.

The serve protocol of the robot is as follows: After orders arrive, each robot, originally staying with the order picker, will leave the picker and go to storage positions to pick one shelf which can provide items to that order. Then the robot will return to the picker it is serving.

Lemma 1 If robot r serves the leftmost or rightmost picker, a single robot serve system can be described as a M/T r/1 queue on condition that W = L.

Proof The arrival process is Poisson. We mainly need to show the inter-serve time follows uniform distribution. We assume the beginning position of robot r is the leftmost picker's position. Note a robot receives their orders information wirelessly and use cameras to read navigational barcode stickers on the warehouse floor, and robots will travel in a straight lines. Supposed a robot receives an order with a storage position X r ∼ U [0, W ] and Y r ∼ U [0, L], let Z = X + Y and its density function is h(z). R robots will carry shelves to the picker. R departure processes merge into an arrival process to the picker. The inter-service time of the picker follows exponential distribution. So the picker service system can be described as a G/M/1 system. We describe the system by Figure 1. 2(1-ρ) is the delay time and 1 µ is the service time. For robot r as a server in a M/G/1 queue, we have

E(S 2 ) = 7(W +L) 2 6v 2 , µ = v W +L , ρ = λr(W +L) v
. Substituting E(S 2 ), µ, ρ into EW , we obtain the expected waiting time at robot r is,

EW r = 7λ r (W + L) 2 12v 2 (1 -λr(W +L) v ) + W + L v . (1)
The equation 1 can be used to design the velocity of robots, we have the following design rule, Design rule 1 There exists a minimal acceptable velocity of robots such that expected waiting time at robots is less than ω R .

Proof Since EW r (v) is decreasing function of v. Let the upper bound of EW r (a measure for service level)at robots be ω R , the minimal acceptable velocity of robots v is given by

7λ r (W + L) 2 12v 2 (1 -λr(W +L) v ) + W + L v = ω R . (2)
Or if the velocity of robots is given by mechanical reasons, the maximal throughput to guarantee a service level ω R is Σ R 1 λ r , and λ r can obtain from equation 2. We have a bound of throughput

T H R = Σ R 1 12v(ω R -(W + L)) (W + L)(12vω R -5(W + L))
.

(3)

We substitute the squared coefficient of variation of an inter-arrival time of robots c 2 aR = 1,ρ = λr(W +L) v ,the squared coefficient of variation of an inter-service time of robots 37,Page 2799),and obtain the squared coefficient of variation of an inter-departure time,

c 2 sR = 1 3 , µ = v W +L ,EW r from equation 1, into the equation c 2 dR = c 2 aR + 2ρ 2 c 2 sR -2ρ(1 -ρ)µEW from Whitt (1983)(Equation
c 2 dR = 1+ 2λ 2 r (W +L) 2 3v 2 -2λ r (1-λr(W +L) v )[ 7λr(W +L) 2 12v 2 (1- λr (W +L) v ) + W +L v ].
(2)Superposition of departure streams from robots

The departure processes are merged into one process, which is also arrival process of the picker. We use methods from [START_REF] Whitt | The queueing network analyzer[END_REF] (see page .

λ aP = Σ R 1 λ r (4)
With the asymptotic method, the squared coefficient of variation of an inter-arrival time of the picker is

c 2 aP = Σ R r=1 λr Σ R r=1 λr {1+ 2λ 2 r (W +L) 2 3v 2 -2λ r (1-λr(W +L) v )[ 7λr(W +L) 2 12v 2 (1- λr (W +L) v ) + W +L v ]}. (3)Analysis of pickers
The expected delay time in a picker

ED P = ρ(c 2 aP + c 2 sP ) 2µ(1 -ρ) g(ρ, c 2 aP , c 2 sP ) (5)
is given by [START_REF] Whitt | The queueing network analyzer[END_REF]

(see page 2802). g(ρ, c 2 aP , c 2 sP ) is defined as exp[- 2(1-ρ)(1-c 2 aP ) 2 ) 3ρ(c 2 aP +c 2 sP )
] when c2 aP < 1, and 1 when c 2 aP ≥ 1. We substitute c 2 sP = 1, and λ aP = Σ R 1 λ r into ED P , and obtain

ED P = ρ(c 2 aP +1) 2µ P (1-ρ) g(ρ, c 2 aP , c 2 sP ) = Σ R 1 λr(c 2 aP +1) 2µ P (µ P -Σ R 1 λr) g(ρ, c 2 aP , c 2 sP ).
We consider the mean number of jobs in queue L q . This is also the number of robots who are carrying shelves and waiting before the picker.

L q = (Σ R 1 λ r ) 2 (c 2 aP + 1) 2µ P (µ P -Σ R 1 λ r ) g(ρ, c 2 aP , c 2 sP ) (6)
In this system, there is an upper bound ι on L q . We show we can find a design rule for velocity of robot.

Design rule 2 There exists a minimal acceptable velocity of robots such that the number of robots who are carrying shelves and waiting before the picker L q is less than space limitation ι.

Proof c 2 aP is a decreasing function of v. g(ρ, c 2 aR , c 2 sR ) is decreasing function of v when c 2 aR < 1, and 1 when c 2 aR ≥ 1. So L q is a decreasing function of v. By solving the following function,

L q (v 2 ) = (Σ R 1 λ r ) 2 (c 2 aP (v 2 ) + 1) 2µ P (µ P -Σ R 1 λ r ) g(ρ, c 2 aP (v 2 ), c 2 sP ) = ι. (7)
We can get a minimal acceptable velocity of robots, v 2 .

(4)Analysis of stability of systems For robots, λ r < µ R = v W +L , r = 1, ..., R. For pickers' workstation, Σ R 1 λ r < µ P . (5)Throughput time We estimate the average customer throughput time of orders, which is a sum of components: the delay time that an order spends in the queue at storage grids waiting for being carried by robots ED r , the average service time of a robot1 µ R , the delay time that order spends in the queue at storage grids waiting for being handled by pickers' workstation ED P , and expected picking time by order pickers 1 µ P . From T r = ED r + 1 µ R + ED P + 1 µ P , we have

T r = 7λr(W +L) 2 12v 2 (1- λr (W +L) v ) + W +L v + 1 µ R + Σ R 1 λr(c 2 aP +1) 2µ P (µ P -Σ R 1 µ P . (6)Optimal robot number
We make two assumptions: 1) We assume the work is equally allocated to R robots, and the arrival rate of the Rth a robot is λ r = λ/R.

2) We assume λ r are equal, then Σ R 1 λ r = λ. We obtain the function of T r(R), T r(R) =

7 λ R (W +L) 2 12v 2 (1- λ(W +L) Rv ) + 1 µ R + λ(c 2 aP +1) 2µ P (µ P -λ) g(ρ, c 2 aP , c 2 sP ) + 1 µ P . T r(R) = 7 λ R (W + L) 2 12v 2 (1 -λ(W +L) Rv ) + 1 µ R + λ(c 2 aP + 1) 2µ P (µ P -λ) g(ρ, c 2 aP , c 2 sP ) + 1 µ P . (8) Lemma 2 ED r (R) is a decreasing function of R.
Proof We can observe from

ED r (R) = 7 λ R (W + L) 2 12v 2 (1 -λ(W +L) Rv ) = 7λ(W + L) 2 12v 2 (R -λ(W +L) v ) . ( 9
) Lemma 3 ED P (R) is a nonlinear function of R. When R ≥ 3λ(W +L) 2v
, ED P (R) is an increasing function.

Proof The expected delay time of the pickers' workstation

ED P (R) = λ(c 2 aP + 1) 2µ P (µ P -λ) g(ρ, c 2 aP , c 2 sP ). ( 10
)
The monotony of ED P (R) is only related to c 2 aP (R). We substitute Σ R r=1 λ r = λ into c 2 aP (R) and obtain,

c 2 aP (R) = 1 + 3λ 2 (W +L) 2 2R 2 v 2 -2λ(W +L) Rv = 1 3 + 3 2 [ λ(W +L) Rv -2 3 ] 2 .
As R increases, the function ED P (R) is nonlinear. It first decreases and then increases after reaching the critical point. When 0 < R ≤ 3λ(W +L) 2v , the function ED P (R) decreases, when R ≥ 3λ(W +L) 2v , ED P (R) increases. By solve the function,we can get the optimal robot number.

Lemma 4 When 0 < R ≤ 3λ(W +L) 2v , T r(R) is a decreasing function, we have a R limited by other constraints. When R ≥ 3λ(W +L) 2v , T r(R) is a convex function. Proof Let ∂T r(R) ∂R = 0, we can obtain R 0 from the func- tion R 2 0 = 6λ 2 µ P (µ P -λ) g(ρ, c 2 aR , c 2 sR )[R 0 -λ(W +L) v ] 2 [ v 2λ(W +L) - 1 R0 ]. If R > 3λ(W +L) 2v
, the optimal number of robots

R * = R 0 . If R ≤ 3λ(W +L) 2v
, the optimal number of robots R * = 3λ(W +L) 2v

.

Pickers share robots

In this section, we consider the situation where one picker and R robots.

Assumptions:

Orders arrive following a Poisson process with the rate of λ. Orders are assigned to available robots following F CF S.

The service rate of the picker processing an order is µ p and the service time follows a general distribution. The picker handles robots in a F CF S. A robot puts back a shelf to the closed available spot after the picker finishes the shelf. A robot does not go back to the picker after returning a shelf, but stay there under the shelf until getting the instruction for the next assigned order. Two queuing systems that are connected in the queue network. M/G/R/F CF S for the first queueing. The service time µ 1 depends on the waiting time and processing time of the second queuing system. M/G/1/F ICF S for the second queuing. The arrival rate depends on the first queue but it should be λ.

(1)Analysis of robots

The expected delay time in a robot

ED r = φ(ρ, c 2 aR , c 2 sR , R)( c 2 aR + c 2 sR 2 )ED r (M/M/R). (11) φ(ρ, c 2 aR , c 2 sR , R) is given by equation(2.25) in Whitt (1993). ED r (M/M/R) = 1 µR(1-ρ) [ (Rρ) R R!(1-ρ) ]
ζ and ζ is given by equation (2.4) in [START_REF] Whitt | Approximations for the gi/g/m queue[END_REF].

The departure process is also arrival process of the picker. We use methods from [START_REF] Whitt | The queueing network analyzer[END_REF] (see page 2799). The squared coefficient of variation of an inter-arrival time of the picker is

c 2 aP = c 2 dR = 1+(1-ρ 2 )(c 2 aR -1)+ ρ 2 √ R (c 2 sR -1) = 1 + ρ 2 √ R (c 2 sR -1).
(3)Analysis of pickers

The expected delay time in a picker

ED P = φ(ρ, c 2 aP , c 2 sP , 1)( c 2 aP + c 2 sP 2 )ED P (M/M/1). ( 12
)
where c 2 aP is given by equation,

c 2 sP = 1, ρ = λ/(µ P ), EW r (M/M/1) = 1 µ(1-ρ) [ ρ (1-ρ) ]
ζ and ζ is given by equation (2.25) in [START_REF] Whitt | Approximations for the gi/g/m queue[END_REF].

We consider the mean number of jobs in queue L q . This is also the number of robots who are carrying shelves and waiting before the picker.

L q = λφ(ρ, c 2 aP , c 2 sP , 1)( c 2 aP + c 2 sP 2
)ED P (M/M/1) (13)

In this system, congestion near pickers' workstation is most serious since aisle width is limited here. There is an upper bound ι on L q . We show we can find a design rule for velocity of robot.

(4)Analysis of stability of systems For robots, λ < Rµ R . For pickers' workstation, λ < µ P .

NUMERICAL EXPERIMENTS

This section shows how the mobile picking systems practices, especially to optimize the number of robots. An important objective for Kiva warehouse automation operations is to determine the optimal number of robots for "Pickers do not share robots" model.

While increasing the number of robots can improve the throughput time, it will also increase the robots fixed cost and robots waiting cost. While increasing the velocity of robots can improve the throughput time, it will also increase the technology cost of the robots. Our objective is to minimize the total operational cost, including the fixed cost of robots, the robots waiting cost and velocity increasing cost. The objective function is the sum of the fixed cost of robots, the waiting cost of robots and velocity increasing cost. The robots waiting cost Σ R r=1 c r λ r EW r, where c R is the marginal cost of a robot, c r is the expected waiting cost per robot, EW r is the total expected waiting time, c r EW r is the total expected waiting cost. The speed increasing cost c v ∆v, where c v is the unit speed increasing cost, and ∆v is the increment of the velocity of robots, we have min c R R + Σ R r=1 c r EW r + c v ∆v. Since v 0 is the initial velocity of robots, and v r is the velocity of r robot,

we have min c R R + Σ R r=1 c r EW r + c v (v r -v 0 ).

Determine the number of robots

We presented λ = 0.4, the width of the warehouse W = 40 and the length of aisles L = 50. Here we ignored other constraints that may affect the optimal number of robots. When v=10 m/min, by solving the three-dimensional equation( 8), R 0 = 2.33, 3.27 + 3.74i or 3.27 -3.74i. As R is an integer, we deleted the imaginary. Since R 0 ≤ 2λ(W +L) v

, we obtained the optimal number of robots R * = 2λ(W +L) v = 7. Similarly, by solving the three-dimensional equation( 8), when v=20 m/min, R 0 = 1.12, 1.43 + 1.93i or 1.43 -1.93i, we obtained the optimal number of robots R * = 4. When v=30 m/min, R 0 = 0.74, 0.91 + 1.29i or 0.91 -1.29i, we obtained the optimal number of robots R * = 2. When v=40 m/min, R 0 = 0.55, 0.67 + 0.96i or 0.67 -0.96i, we obtained the optimal number of robots R * = 1. When v=50 m/min, R 0 = 0.44, 0.53 + 0.78i or 0.53 -0.78i, we obtained the optimal number of robots R * = 1.

Figure 2 shows that with the increase of the number of robots, the throughput time non-linear increases. Since R is an integer, when v=10 m/min, the number of robots is less than seven, the trend of throughput time declines, and then begin to increase. When v=20 m/min,the number of robots is less than four, the trend of throughput time declines,and then begin to increase. When v=30 m/min,the number of robots is less than two, the trend of throughput time declines,and then begin to increase. When v=40 m/min, the number of robots is less than one, the trend of throughput time declines, and then begin to increase. When v=50 m/min, the number of robots is less

CONCLUDING

This paper studies a new order fulfilling system featured by a number of robots transporting movables shelves to order pickers. This new generation of sustainable green warehouse systems can improve productivity and flexibility simultaneously. We measure the performance of system, and provide design rules for warehouses. We compare two protocols: To share or no to share robots. The models is this paper can be used to determine the number of robots, determine the velocity of robots, and compute throughput given service level.This is one of earliest papers to introduce this new material handling system.

  Fig. 1. A queue networks model for Kiva system (2002) developed efficient algorithms for AGV's scheduling and routing problems. The focus strategy of our paper is on a random selection strategy and the minimize throughput time selection strategy.

( 1 )

 1 Analysis of robots We use ED to denote the expected delay time, and ES to denote the expected service time. The expected waiting time EW = ED + ES. From Wolff (1989)(Equation 154,Page 280), we have EW = λE(S 2 ) 2(1-ρ) + 1 µ , where ED = λE(S 2 )

Fig. 2 .Fig. 3 .

 23 Fig. 2. Throughput time with the change of Rthan one, the trend of throughput time declines, and then begin to increase. When the number of robots is larger than 30, the throughput time have no great changes.3.2 Determine the velocity of robots

λr) g(ρ, c

aP , c 2 sP )+


Francisco J Blanco-Rodríguez

Stéphane Le Dizès

Elliptic instability of a curved Batchelor vortex -corrigendum

. The correct values are twice those indicated in this paper for all the modes. This modifies the values given in table 2 and formulas (C2n-q), (C3n-q), (C4n-q). For instance, in table 2 the correct value of R AB for the mode (-2, 0, 1) at

)).

 for the Rankine vortex is σ max /ε 2 = [(0.428 log(8/ε) -0.455) 2 -0.113] 1/2 while we obtain for the Lamb-Oseen vortex σ max /ε 2 = 0.5171 log(8/ε) -0.9285. The Lamb-Oseen vortex ring is thus less unstable than the Rankine vortex ring as soon as ε > 0.039 for the same reason as previously indicated.

Figure 8. Maximum elliptic instability growth rate σmax/ε 2 = σ * max 2πR 2 /Γ versus ε = a/R for a vortex ring of radius R (a), and 2 counter-rotating vortices distant of 2R (b), both of core size a. Black lines: W0 = 0 (mode (-1, 1, 1)), red lines: W0 = 0.2 (mode (-2, 0, 2)), green lines: W0 = 0.4 (mode (-2, 0, 1)). Solid lines are for Re = ∞, dashed lines for Re = 5000. For both cases, the perturbation wavenumber is assumed not to be discretized (see figure 9).

0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 1.2 ε σ * max 2πR 2 /Γ 0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 1.2 ε 0 0.1 0.2 0.3 0 0.2 0.4 0.6 0.8 1 1.2 ε (a) (b) (c)
Figure 9. Elliptic instability growth rate σmax/ε 2 = σ * max 2πR 2 /Γ for a vortex ring of radius R = 1/ε for Re = ∞ (solid lines) and Re = 5000 (dashed lines). Thick lines are the maximum growth rate curves plotted in figure 8(a), thin lines are the growth rate curves assuming that the instability wavenumbers are discretized and only take the values kn = 2πn ε, n = 1, 2, 3, • • •. (a): W0 = 0 (mode (-1, 1, 1)), (b): W0 = 0.2 (mode (-2, 0, 2)), (c): W0 = 0.4 (mode (-2, 0, 1)). 1,1)), red lines: W0 = 0.2 (mode (-2, 0, 2)), green lines: W0 = 0.4 (mode (-2, 1)). Left plot: normalization with Γ/(2πL 2 z ). Right plot: normalization with Γ/(2πR 2 ). The vertical lines on the left plot indicate the value of Lz/R where the local strain rate vanishes for each configuration. . Maximum elliptic instability growth rate versus L/R of two helical vortices (solid lines) and of an double-array of straight co-rotating vortices (dashed lines) for a/R = 0.1 and Re = ∞. Black lines: W0 = 0 (mode (-1, 1, 1)), red lines: W0 = 0.2 (mode (-2, 0, 2)), green lines: W0 = 0.4 (mode (-2, 0, 1)). Left plot: normalization with Γ/(2πL 2 z ). Right plot: normalization with Γ/(2πR 2 ). The vertical lines on the left plot indicate the value of L/R where the local strain rate vanishes for each configuration.

F

Figure 10 .

 10 Figure10. Maximum elliptic instability growth rate versus Lz/R of an array of alternate vortex rings (solid lines) and of an array of straight counter-rotating vortex pairs of alternate sign (dashed lines) for ε = a/R = 0.1 and Re = ∞. Black lines: W0 = 0 (mode (-1, 1, 1)), red lines: W0 = 0.2 (mode (-2, 0, 2)), green lines: W0 = 0.4 (mode (-2, 1)). Left plot: normalization with Γ/(2πL 2 z ). Right plot: normalization with Γ/(2πR 2 ). The vertical lines on the left plot indicate the value of Lz/R where the local strain rate vanishes for each configuration.

  Figure11. Maximum elliptic instability growth rate versus L/R of two helical vortices (solid lines) and of an double-array of straight co-rotating vortices (dashed lines) for a/R = 0.1 and Re = ∞. Black lines: W0 = 0 (mode (-1, 1, 1)), red lines: W0 = 0.2 (mode (-2, 0, 2)), green lines: W0 = 0.4 (mode (-2, 0, 1)). Left plot: normalization with Γ/(2πL 2 z ). Right plot: normalization with Γ/(2πR 2 ). The vertical lines on the left plot indicate the value of L/R where the local strain rate vanishes for each configuration.


A Pierquin

T Henneron
email: thomas.henneron@univ-lille1.fr

S Clénet

S Brisset

Model Order Reduction of Magnetoquasistatic Problems Based on POD and Arnoldi-based Krylov Methods

Keywords: Krylov subspace, Model Order Reduction, Proper Orthogonal Decomposition, Quasi-static problem

The Proper Orthogonal Decomposition method and the Arnoldi-based Krylov projection method are investigated in order to reduce a finite element model of a quasistatic problem. Both methods are compared on an academic example in terms of computation time and precision.

I. INTRODUCTION

T O analyse electromagnetic devices, the Finite Element Method (FEM) associated with a time-stepping scheme is often used to solve Maxwell's equations. The discretisation of the space and time domains can yield a large-scale equation system which computation can be dramatically prohibitive when a fine mesh and a high number of time steps are used. To tackle this issue, model order reduction approaches can be investigated. These approaches enable to reduce the size of the problem to solve. In the literature, the Proper Orthogonal Decomposition (POD) [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] and the Arnoldi-based Krylov projection (AKP) [START_REF] Wittig | Two-step lanczos algorithm for model order reduction[END_REF] approaches have been used to solve many problems in engineering. Both methods consist in performing a projection of the solution onto a reduced basis. For the POD approach, the projector is calculated from a set of solutions of the problem (the snapshots) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, part I: Coherent structures[END_REF]. For the AKP approach, the projector is deduced from a Krylov subspace built from the equation system to solve in the frequency domain. In computational electromagnetics, the POD method has been already applied to study magnetoquasistatic and electroquasistatic problems [START_REF] Henneron | Model order reduction of quasi-static problems based on pod and pgd approaches[END_REF], [START_REF] Schmidthausler | Reduction of linear subdomains for non-linear electro-quasistatic field simulations[END_REF]. The AKP approach has been used to solve magnetoharmonic problem using the partial element equivalent circuit (PEEC) method [START_REF] Nguyen | Equivalent circuit synthesis method for reduced order models of large scale inductive peec circuits[END_REF]. To the author's knowledge, any linear electromagnetic problem hasn't been solved in the time domain using the projector constructed from AKP method in the frequency domain.

In this paper, we propose to apply and compare the POD and AKP methods in the case of quasistatic problems solved by the vector potential formulation in the time domain and the frequency domain. First, the numerical model deduced from the vector formulation is presented. Second, the POD and the AKP methods are described. Finally, an eddy current problem in time and frequency domains is studied by the proposed reduction approaches. The results obtained with reduced models are compared with those deduced from the full model. The two model order reduction methods are also compared in terms of accuracy, computation time and complexity. 1). The problem is solved for x ∈ Ω and for t in the time interval [0, T ]. A magnetoquasistatic problem can be solved by using the vector potential formulation. The potential A is defined such that B(x, t) = curlA (x, t), E(x, t) = -∂A (x, t) ∂t and A (x, t) × n = 0 on Γ B and Γ E where B is the magnetic flux density and E the electric field. These expressions are combined with Ampère's equation in order to obtain the strong formulation

II. LINEAR MAGNETOQUASISTATIC PROBLEM Let us consider a domain

Ω of boundary Γ = Γ H ∪ Γ B , Γ H ∩ Γ B = ∅, with a stranded inductor Ω s and a conducting subdomain Ω c of boundary Γ c = Γ E ∪ Γ Jed , with Γ E ⊂ Γ B (Fig.
curl 1 µ curlA (x, t) + σ ∂A ∂t = N (x)i(t), (1)
with µ the magnetic permeability, σ the electric conductivity, N the unit current density vector in Ω s and i the current flowing through the stranded inductor.

To solve this problem, the potential A is discretized using edge elements and N using facet elements. Applying the Galerkin method to (1), a system of ordinary differential equations is obtained:

M dX(t) dt + KX(t) = F i(t). (2)
The vector of unknowns X(t) is composed of the circulations of A along all edges of the mesh, X(t) = [A 1 (t), . . . , A Nun (t)] T . In the case of a sinusoidal supply of frequency f , the problem can be solved in the frequency domain to determine the steady state. Then, the problem is solved for x ∈ Ω and for f ∈ [f min , f max ]. In this case, ( 2) becomes

(M jω + K)X(ω) = F i(ω) (3)
with j the imaginary unit and ω = 2πf .

III. MODEL ORDER REDUCTION

The principle of the model order reduction is to obtain from (2) or (3) an equation system of reduced size. An orthonormal reduced basis Ψ is determined to approximate the solution X of the full problem under the form ΨX r , with X r the solution of the reduced problem:

M r dX r (t) dt + K r X r (t) = F r i(t) (4)
or

(M r jω + K r )X r (ω) = F r i(ω) (5)
with

M r = Ψ T M Ψ, K r = Ψ T KΨ, F r = Ψ T F .
The Krylov subspace projection and the POD approach are two different ways to determine the reduced basis Ψ.

A. Krylov Projection method

For the AKP approach, the solution of the full model is projected onto a Krylov subspace. To determine this subspace, an expansion of the transfer function of the problem is performed. The Laplace's transform is applied to (2):

(M s + K)X(s) = F i(s). ( 6
)
The transfer function has the following expression (see ( 6))

h(s) = (M s + K) -1 F . ( 7
)
h is approximated by a Padé's expansion at the expansion point s exp

h(s) = ∞ j=0 h j (s -s exp ) j . ( 8
)
The vectors h j can be expressed as

h j = -(M s exp + K) -1 M j (M s exp + K) -1 F (9) = G j g, with G = -(M s exp + K) -1 M and g = (M s exp + K) -1 F .
The set of vectors (h j ) j=0,Ns-1 generates the Krylov subspace K Ns (G, g) = span{g, G 1 g, G 2 g, . . . , G Ns-1 g}. The Arnoldi's algorithm is used to obtain an orthonormal basis Ψ of the Krylov subspace K Ns (G, g) [START_REF] Odabasioglu | Prima: passive reducedorder interconnect macromodeling algorithm[END_REF]. Thus, Ψ allows to project the solution X onto K Ns (G, g). Several expansion points s k can be used to express (7):

h(s) = ∞ j k =0 h j k (s -s k ) j k , ( 10
)
with h j k = -(M s k + K) -1 M j k (M s k + K) -1 F . The
h j k vectors associated with each expansion point s k are socalled point moments. The expansion point s k will be equal to 2πf k or j2πf k according to the problem [START_REF] Feldmann | Efficient linear circuit analysis by Pade approximation via the lanczos process[END_REF], [START_REF] Cangellaris | Electromagnetic model order reduction for system-level modeling[END_REF].

B. Proper Orthogonal Decomposition

With the snapshot POD approach [START_REF] Sirovich | Turbulence and the dynamics of coherent structures, part I: Coherent structures[END_REF], the projector Ψ is defined from a set of solutions of the full model (i.e. snapshots). The problem (2) is solved for the N s first time-steps to obtain X(t i ), i = 1, . . . , N s , then these solutions are concatenated in the snapshot matrix S = [X(t 1 ), . . . , X(t Ns )]. In the frequency domain, (3) is solved for N s different frequencies and the snapshot matrix is defined by S = [X(ω 1 ), . . . , X(ω Ns )].

Applying the singular value decomposition to the matrix S gives

S = V ΣW T (11)
with V Nun×Nun and W Ns×Ns orthonormal matrices and Σ Nun×Ns a rectangular diagonal matrix of the singular values.

The Ψ basis corresponds to the normalised matrix V Σ or SW . In practice, it is more efficient to compute W by solving the eigenvalues problem on the correlation matrix C of S because C has a smaller size:

C = 1 N s S T S = W Σ T V T V ΣW T = W ∆W T . ( 12
)
Then the columns of SW = V Σ are normalised to obtain an orthonormal basis Ψ.

IV. APPLICATIONS

The application example is composed of two conducting plates submitted to a magnetic field created by a stranded inductor. Due to its symmetries, only one eighth of the geometry is modeled (Fig. 2). The two model order reduction methods are applied in time-domain and in frequency-domain. In the time-domain, the stranded inductor is supplied by a current i(t) = sin(2π5000t). In the frequency-domain, a current with a magnitude equal to one is imposed for all frequencies. The results obtained with the reduced models are compared with those deduced from the full model.

Stranded inductor

Conducting plate Fig. 2. Geometry and mesh of the device for the stranded inductor and the conducting plate.

A. Quasi static problem

For the quasi static problem in the time domain, the full model ( 2) is solved by using an implicit Euler scheme. The modeling is performed on 4 periods with 20 points per period. For the POD method, the snapshots are the solutions of the N s first time steps of the full model. For the AKP approach, only one expansion point is chosen for the Padé's expansion. The value of the expansion point is equal to the angular velocity of the current: s exp = 2π5000.

In order to evaluate the influence of the size of the reduced basis on a global quantity, the evolutions of the Joule losses versus time, obtained from the POD and AKP models, are compared with the one of the full model in Fig. 3 and Fig. 4. We can observe that the Joule losses converge towards the reference with both approaches when the size of the reduced basis increases. Figures 5(a For a given size of the reduced basis, the error obtained with the AKP method is smaller than the error obtained with the POD model. Nevertheless, both approaches give the same error when the size of the reduced basis is higher than 12. Figure 6(a) presents the distribution of the eddy current density at a given time step in a cross section of the conducting plate obtained from the full model. Figures 6(b) and 6(c) present the norm of the difference between the eddy current density from the full model and those deduced from the POD and AKP models. The error density is smaller with the AKP method, but considering the scale, the error of the POD model remains acceptable. Figure 7 presents the computation time of the reduced models with respect to the reduced basis size. For both approaches, the computation time is smaller than the one required for the full model. The POD model is faster than the AKP model. This difference comes from the worse matrix conditioning obtained with the AKP method. In fact, with the POD method, due to the Euler scheme, the matrix to inverse is (M 1 dt + K), and with the AKP approach, the matrix to inverse is (M s exp + K). The profile of the two matrices is the same, but the values of 1 dt and s exp are not equal in both cases, that changes consequently the matrix conditioning, and then changes the convergence of the biconjugate gradient (BiCG) used to solve the system. The convergence of BiCG is longer in the case of AKP.

For non linear systems, the POD method can be applied in the same way as in linear case. Recently, the POD method has been associated with the Discrete Empirical Interpolation Method to speed up the accounting of the non linearity with the POD. For electromagnetic problem, this kind of approach has been used to study a single phase transformer [START_REF] Henneron | Model order reduction of non-linear magnetostatic problems based on POD and DEI methods[END_REF]. However, the AKP method can not be used directly but the nonlinear behaviour can be taken into account by linearization [START_REF] Bai | Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems[END_REF].

B. Magnetoharmonic problem

For the magnetoharmonic problem, we consider a frequency interval from 100 Hz to 20 kHz with a frequency step equal to 100 Hz. For this case, we choose snapshots and expansion points uniformly distributed on the frequency interval. The reduced basis size is fixed at 6. For the POD approach, six snapshots are uniformly distributed on the frequency spectrum. For the AKP method, 3 expansion points are chosen (at 100 Hz, 10 kHz and 20 kHz) and for each expansion point the two first moments are computed in order to obtain the same size for the two reduced models. Table I presents the error on the magnetic energy and on the Joule losses introduced by the reduced models. For this example, the errors of both reduced models are in a same range. Figure 8(a) presents the real part of the distribution of the eddy current density in a cross section of the conducting plate obtained from the full model. Figures 8(b) and 8(c) present the difference between the real part of the eddy current density from the full model and those deduced from the POD and AKP models. As for the first study in the time domain, the error density is smaller with the AKP method. We can also observe that the eddy current density obtained from both reduced models are close to the one deduced from the full model. The computation time of the reduced models are similar; 199s for the POD approach and 208s for the AKP method. The full model requires 1907s to compute the solution on the frequency interval. Then, the speed up is greater than 9.

V. CONCLUSION

The Proper Orthogonal Decomposition method and the Arnoldi-based Krylov projection method combined with a FEM vector potential formulation have been developed in order to solve 3D magnetodynamic and magnetoharmonic problems. For the studied example, it has been shown that the two approaches enable to reduce the computation time significantly. The errors on the global and local quantities introduced by both reduced models are small compared with the results obtained from the full model. For a given size of the reduced model, the errors associated with the AKP method are smaller than those of the POD approach. Nevertheless, the computation time of the reduced model from the POD method is smaller than this of the AKP approach.

Fig. 1 .

 1 Fig. 1. Spatial domain of the problem.

  ) and 5(b) present in logscale the error on the Joule losses and on the magnetic energy,εP = P ref -P red P ref and ε E = E ref -E red E ref with P ref and P red (resp. E ref and E red ) the evolution of the Joule losses (resp. magnetic energy) obtained from the full and reduced models.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Evolution of the Joule losses obtained from the POD model with 1,2 and 3 snapshots.

  Error on the Magnetic energy.

Fig. 5 .

 5 Fig. 5. Error with respect to the reduced basis size.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Eddy current density of the full model and error density for a reduced basis of size 8, in the conducting plate viewed from the top.

Fig. 8 .

 8 Fig. 8. Real part of the eddy current density at 10 kHz and error density, in the conducting plate viewed from the top.

TABLE I ERROR

 I ON MAGNETIC ENERGY AND JOULE LOSSES FOR REDUCED MODELS OF SIZE 6.

	Reduced model	POD	Krylov-Arnoldi
	E ref -E red E ref	1, 7656.10 -6	1, 7322.10 -6
	P ref -P red P ref	2, 0241.10 -5	2, 3489.10 -5


archive ouverte pluridisciplinaire HAL, est destinée au dépôt

Theoretical predictions for the elliptical instability in a two-vortex flow -corrigendum

Stéphane Le Dizès, Florent Laporte

To cite this version:

Stéphane Le Dizès, Florent Laporte. Theoretical predictions for the elliptical instability in a two-vortex flow -corrigendum. 2017. hal-01430605

L'

CORRIGENDUM

As pointed out by Pantxika Otheguy, a few misprints are present in formulas of section 6. In expressions (6.1a) and (6.2a), b 2 /a 2 1 under the square root should read b 4 /a 4 1 . Similarly, in expressions (6.1b) and (6.2b), b 2 /a 2 2 under the square root should read b 4 /a 4 2 . These misprints do not affect the results which were obtained with the correct expressions.