Fatemeh Chegini
email: chegini@zib.de

Algiane Froehly
email: algiane.froehly@inria.fr

Ngoc Mai

Monica Huynh
email: ngocmaimonica.huynh@unipv.it

Luca Pavarino
email: luca.pavarino@unipv.it

Mark Potse
email: mark.potse@inria.fr

Simone Scacchi
email: simone.scacchi@unimi.it

Martin Weiser
email: weiser@zib.de

N M M Huynh

EFFICIENT NUMERICAL METHODS FOR SIMULATING CARDIAC ELECTROPHYSIOLOGY WITH CELLULAR RESOLUTION

Keywords: EMI model, high-order time integration, spectral deferred correction, algebraic adaptivity, multirate integration scheme, ladder methods, BDDC preconditioner

The cardiac extracellular-membrane-intracellular (EMI) model enables the precise geometrical representation and resolution of aggregates of individual myocytes. As a result, it not only yields more accurate simulations of cardiac excitation compared to homogenized models but also presents the challenge of solving much larger problems. In this paper, we introduce recent advancements in three key areas: (i) the creation of artificial, yet realistic grids, (ii) efficient higher-order time stepping achieved by combining low-overhead spatial adaptivity on the algebraic level with progressive spectral deferred correction methods, and (iii) substructuring domain decomposition preconditioners tailored to address the complexities of heterogeneous problem structures. The efficiency gains of these proposed methods are demonstrated through numerical results on cardiac meshes of different sizes.

INTRODUCTION

Cardiac arrhythmias are significant causes of mortality, and numerical modeling plays a crucial role in understanding their mechanisms and designing effective treatments. While homogenized models like the bidomain and monodomain models represent the action potential at the tissue level, they lack details about arrhythmogenic cellular-level factors such as fibrosis or ion channel distribution. To address this, we consider a recent cell-by-cell model that explicitly represents myocytes. Due to the size of such models, they pose new simulation challenges related to mesh generation, time stepping, and preconditioning.

In this paper, we present a geometry and mesh generation approach for defining large-scale computational meshes and introduce efficient adaptive time stepping and domain decomposition preconditioning techniques. We employ tetrahedral meshes based on the implicit surface meshing of a random network of cells, which mimics real cardiac tissue [START_REF] Potse | A practical algorithm to build geometric models of cardiac muscle structure[END_REF]. This method has been extended to represent 400 cells in 15 million tetrahedra.

Although the cell-by-cell model exhibits localized solution features, the overhead induced by mesh adaptivity during simulation has proven to outweigh efficiency gains. To overcome this, we combine the spectral deferred correction (SDC) method as a higher-order time stepping approach with computationally inexpensive algebraic adaptivity [START_REF] Chegini | Adaptive multirate integration of cardiac electrophysiology with spectral deferred correction methods[END_REF]. This combination, which can be interpreted as multi-rate integration, significantly speeds up simulations in both homogenized and cell-by-cell models [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF]. The approach utilizes nested subset selection on the algebraic level and explicit embedded error estimates, incurring only negligible overhead.

Since the dominant diffusion and elliptic constraints in the cell-by-cell models necessitate implicit time-stepping methods, we also propose a novel domain decomposition preconditioner called balancing domain decomposition by constraints (BDDC) [START_REF] Huynh | Convergence analysis of BDDC preconditioners for hybrid DG discretizations of the cardiac cell-by-cell model[END_REF]. This preconditioner is specifically tailored to the unique structure of cell-by-cell models, and we investigate its convergence properties numerically.

To showcase the overall performance and convergence rate of our proposed approach, we provide demonstrations on 3D cardiac simulations using a combination of SDC, algebraic adaptivity, and BDDC preconditioning.

THE EMI MODEL

The extracellular-membrane-intracellular (EMI) model of electrophysiology [START_REF] Jaeger | Properties of cardiac conduction in a cell-based computational model[END_REF] describes the myocardium as a collection of pairwise disjoint myocytes (Ω i ) i=1,...,N which, together with the extracellular space Ω 0 , cover the whole domain Ω ⊂ R d , d ∈ {2, 3}, occupied by the myocardium, i.e. Ω = N i=0 Ω i . Ions can diffuse freely within each myocyte and in the extracellular domain, subject to conductivities σ i , which leads to electric intra-and extracellular potentials u i ∈ H 1 (Ω i ). Gap junctions between myocytes allow ion currents through passive channels, while active channels regulate ion exchange between myocytes and the extracellular space. The transmembrane current n T σ i ∇u i includes the ion current I ion ij and the capacitive current C m vij . The ion current depends on the transmembrane voltage v ij = u i -u j defined on the membrane F ij = ∂Ω i ∩ ∂Ω j , and w ij comprising ion concentrations and channel states governed by a nonlinear dynamic given by R. This setting results in the partial differential algebraic system

               -div(σ i ∇u i ) = 0 in Ω i -n T i σ i ∇u i = C m vij + I ion ij (v ij , w ij ) on F ij , i ̸ = j n T σ i ∇u i + ϵu i = 0 on ∂Ω i ∩ ∂Ω ẇij = R(v ij , w ij ) on F ij , ij = 0 w ij = 0 ij > 0. (1)
Here, n i denotes the unit outer normal of Ω i for i = 0, . . . , N . We require

I ion ij (v ij , w ij ) = -I ion
ji (v ji , w ji ) and w ij = w ji . The small value ϵ > 0 in the Robin boundary condition on ∂Ω makes the solution unique, which for a pure Neumann problem would be defined only up to a constant. An alternative normalization would be to require Ω ϕu dx = 0 for some nonvanishing function ϕ. The parameter values used for the EMI model are given in Tab. 1.

MESH GENERATION

Computational meshes were generated as described previously [START_REF] Potse | A practical algorithm to build geometric models of cardiac muscle structure[END_REF], with modifications that allowed us to produce larger samples. First we generated a random network of cell centers and links, which were not allowed to deviate more than 15 degrees from a defined "fiber orientation". Each half-link was then assigned the part of the tissue volume that was nearest to it.

Parameter Value

σ >0 0.4 S m -1 σ 0 2.0 S m -1 C m 1 × 10 -4 F m -2
Table 1: Parameters used in the EMI model, based on [START_REF] Jaeger | Properties of cardiac conduction in a cell-based computational model[END_REF]. The boundaries between these domains were discretized inside a tetrahedral mesh. Using these boundaries and the links themselves we determined a level-set function that defined the cell membranes. The level set was discretized by the Mmg remeshing software [START_REF] Dapogny | Three-dimensional adaptive domain remeshing, implicit domain meshing, and applications to free and moving boundary problems[END_REF]. The halflink domains were preserved through level-set discretization and then merged for each center, resulting in cell domains that determined the locations of the intercalated disks, the structures that connect the cells electrically and mechanically. The resulting poor-quality mesh was then improved with a second application of Mmg. In contrast to our previous work [START_REF] Potse | A practical algorithm to build geometric models of cardiac muscle structure[END_REF], we merged the extracellular parts of the cell domains before remeshing. This made the remesher's task much easier, because the junctions of multiple domain boundaries, often joining at small angles, were hard to preserve. This modification, together with recent improvements in the Mmg software, allowed us to generate samples that were orders of magnitude larger than the few dozen myocytes reported before.

A domain decomposition for parallel computations was made by combining the intracellular domains of several cells and expanding them until the extracellular volume was completely partitioned as well. Thus, cells were never split across partition boundaries. Necessarily, the partition boundaries could coincide with cell membranes and intercalated disks.

Fig. 1 shows an overview of the mesh generation methods, Fig. 4 the four meshes used in this study, and Fig. 2 the intra-and extracellular domains of a single cell.

ALGEBRAIC ADAPTIVITY IN HIGH ORDER TIME STEPPING

The depolarization and repolarization fronts travelling through the domain are highly localized, calling for high spatial resolution at the front position but allowing for much coarser meshes in regions far from the front, reducing the problem size considerably. This can be achieved by adaptive mesh refinement and coarsening [START_REF] Colli Franzone | Adaptivity in space and time for reaction-diffusion systems in electrocardiology[END_REF]. Unfortunately, the overhead of error estimation, mesh modification, and re-assembly often compensates the reduction of computational work due to the smaller problems to be solved. Here we consider a recent approach to spatial adaptivity with spectral deferred correction time stepping [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF][START_REF] Chegini | Adaptive multirate integration of cardiac electrophysiology with spectral deferred correction methods[END_REF] that works purely on the algebraic level without incurring significant overhead and thus can achieve a substantial speedup. Combining this algebraic adaptivity with ladder methods [START_REF] Minion | Semi-implicit spectral deferred correction methods for ordinary differential equations[END_REF] improves the efficiency gain further.

Algebraic adaptivity with spectral deferred corrections

Spatial semi-discretization of (1) with finite elements results in a large scale differentialalgebraic system

B ż = Az + F (z) (2)
of index one for z = [u, w] ∈ R N h +n h . Here, A contains the infinitely stiff diffusion term, and F represents the nonlinear ion channel dynamics. Note that w is restricted to the membranes, but may be vectorial, whereas u is defined on the whole domain, but is scalar. Consequently, N h ̸ = n h . Since time derivatives only appear on the membranes, the matrix B is singular, such that implicit L-stable integration methods such as an implicit-explicit (IMEX) Euler scheme

(B -τ A)z k+1 = z k + τ F (z k ) (3)
for a step size τ have to be used [START_REF] Hairer | Solving ordinary differential equations. II: Stiff and differentialalgebraic problems[END_REF]. Spectral deferred correction (SDC) methods turn a low order basic scheme such as (3) into a higher order scheme by solving the defect equation

B δz = Aδz + F (z + δz) -F (z) + [Az + F (z) -ż]
for a given approximate trajectory z on a certain time horizon, using the basic scheme [START_REF] Dutt | Spectral deferred correction methods for ordinary differential equations[END_REF][START_REF] Weiser | Faster SDC convergence on non-equidistant grids by DIRK sweeps[END_REF]. The term in brackets is independent of δz and can be treated explicitly by high-order quadrature on a collocation time grid. Updating the initial guess z ← z + δz and repeating yields a stationary iteration that converges to a high order collocation solution. Choosing a Radau quadrature with n collocation points and taking m iterations yields an L-stable IMEX Runge-Kutta scheme of order min{2n -1, m} with nm stages.

In electrophysiology simulations, the SDC corrections have a dominantly local support around the fronts. Based on an error estimator based on the linear convergence of the stationary iteration, thresholding corrections for determining the relevant subregion still requiring higher accuracy, and restricting subsequent iterations to that subregion has been shown to reach the desired uniform tolerance [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF]. The selection of the subregions and the formulation of the restricted problems can be done purely on the algebraic level of finite element coefficients, as shown in Fig. 5, which is computationally cheap and thus turns the reduction in problem size directly into faster solution. Wall clock time speedup factors over a conventional SDC iteration are given in Tab. 3 for grids of different size. Since the depolarization front width stays fixed, the speedup increases slightly with the domain size.

Ladder methods

The speedup that can be achieved with the algebraic adaptivity outlined above is limited by the number m of SDC iterations. This holds for the extreme case when after the first iteration, which needs to be performed on the whole domain, the following iterations involve just a negligible number of degrees of freedom. A further improvement is possible if the first iteration does not involve all n collocation points, but just one. So-called ladder methods, refining the collocation grid every second iteration, reduce the computational cost, but also yield a larger discretization error, and have been found not to improve the efficiency of ODE integration because they require more iterations for the same accuracy [START_REF] Minion | Semi-implicit spectral deferred correction methods for ordinary differential equations[END_REF].

This changes if the computational cost for the iterations is not uniform. With algebraic adaptivity, the computational cost of the first iteration is much higher than for the later iterations. Consequently, the upper bound for the speedup is increased from m to nm. As a very coarse model of impact of ladder methods on accuracy and effort, we will assume the following:

1. Standard SDC converges linearly with a contraction factor ρ < 1. This is in general satisfied, with 0.15 ≤ ρ ≤ 0.6 the usual range for reaction-diffusion equations [START_REF] Weiser | Faster SDC convergence on non-equidistant grids by DIRK sweeps[END_REF].

2. Algebraic adaptivity does not affect the reached accuracy and hence convergence rate. With a proper choice of drop tolerance, this is in general the case [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF].

3. The number of degrees of freedom retained in iteration i scales as N r i for some r < 1.

The actual number of retained degrees of freedom is, of course, highly problem dependent, but a factor r = 0.5 appears to be realistic.

4.

The larger basic time step in ladder methods increases the contraction factor to ρ q , with q < 1. Asymptotically, this is a rather pessimistic assumption, since after 2n iterations, the full collocation grid is reached and the ladder method's contraction rate reaches ρ. Up to then, of course, larger contraction factors are observed.

With these assumptions, and performing m = 2n -1 iterations to reach the full time discretization order of Radau collocation, we face a computational effort proportional to for standard SDC with algebraic adaptivity. For ladder methods with algebraic adaptivity, the number of iterations is about m = (2n -1)/q to reach the same truncation error. With a slower convergence rate ρ q , the number of retained degrees of freedom scales as N r iq and thus decreases slower, but there are also just ⌊i/2 + 1⌋ collocation points in each iteration. Thus, with a = r q and m = ⌈(2n -1)/q⌉, the computational work is roughly proportional to

W := m-1 i=0 nr i = n 1 -r 2n-1 1 -r F.
⌈(2n-1)/q⌉ i=0 i 2 + 1 r qi ≈ 2 -a + a m+1 ((m + 2)a -(m + 3)) 2(1 -a) 2 =: W lad . ( 4
)
The expected speedup W/W lad due to using ladder methods with algebraic adaptivity is shown in Fig. 3 for r = 0.5 and different values of n and q. As can be expected, larger speedups are predicted for more collocation points (larger n) and less impact of the ladder approach on the convergence rate (larger q). Also apparent is that for small n and q, an actual slowdown must be expected. The empirical speedup observed in the examples considered here is reported in the last column of Tab. 3.

DOMAIN DECOMPOSITION PRECONDITIONERS

The basic IMEX Euler scheme (3) requires the solution of globally coupled equation systems of the form B -τ A. The matrix B -τ A is positive definite if the nullspace of the pure Neumann problem is eliminated by applying Robin boundary conditions as defined in [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF]. For large grids with many myocytes as shown in Fig. 4, iterative solvers are essential, and the method of choice is preconditioned conjugate gradients (PCG).

Among the numerous available options we consider a BDDC preconditioner [START_REF] Dohrmann | A preconditioner for substructuring based on constrained energy minimization[END_REF] specifically tailored to the geometrical structure of the EMI model [START_REF] Huynh | Convergence analysis of BDDC preconditioners for hybrid DG discretizations of the cardiac cell-by-cell model[END_REF]. BDDC splits the mesh in several non-overlapping subdomains and decomposes the system matrix M = B -τ A into a corresponding block-diagonal system by replicating degrees of freedom associated to the interface between subdomains in each block. Approximate global continuity is preserved by adding few constraints such as continuity at subdomain corners or vanishing average jumps over subdomain interfaces. Since there are very few such constraints, the Schur complement with respect to the block-diagonal part is relatively small, and can therefore be solved efficiently using direct solvers. This yields an approximate solution with artificial jumps across subdomain interfaces, which are subsequently removed by averaging the solution at the interfaces and computing local solutions on subdomains given these Dirichlet boundary data. The resulting approximate solution is usually quite good, and thus BDDC used as a preconditioner exhibits very low condition numbers of the preconditioned system. The EMI geometry suggests a natural decomposition of the computational domain Ω along the cell membranes into subdomains Ω i representing the myocytes and extracellular domain. Since the potentials u i are discontinuous across membranes, there are two different degrees of freedom associated to each interface vertex, and both of them are replicated in both subdomains. This BDDC design has been proved to exhibit the expected almost mesh-independent condition bound and demonstrated to yield a highly effective preconditioner in 2D problems [START_REF] Huynh | Convergence analysis of BDDC preconditioners for hybrid DG discretizations of the cardiac cell-by-cell model[END_REF].

Here, we consider larger 3D problems. This leads to a large extracellular subdomain, which is challenging for direct solvers and causes imbalance with the many smaller myocyte subdomains. We therefore subdivide the extracellular space into several subdomains. Introducing an interface condition of the same mathematical structure as gap junctions between these subdomains allows extending the BDDC design seamlessly to subdivided extracellular space. On these interfaces, we impose an ion current I ion (v) = σ virt v with a large conductivity σ virt , realizing an interior penalty discontinuous Galerking discretization, such that the solution is almost continuous.

NUMERICAL RESULTS

We have applied the algorithmic approaches described in Sec. 4 and 5 to four grids of different size shown in Fig. 4. Their properties are given in Tab. 2. Since the potential values jump across membranes, the number of degrees of freedom in a piecewise linear finite element discretization exceeds the number of mesh vertices by a small fraction.

We simulated the EMI model ( 1) on these grids with SDC on a Radau IIA time grid with three collocation points. The number of degrees of freedom retained in every SDC iteration is shown in Fig. 5 over the course of the whole simulation. The clearly visible significant reduction of system size translates into a reduction of wall clock simulation time, given as speedup over plain SDC in Tab. 3. Since the domains are larger than the one considered before [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF], but the excitation front width remains the same, the observed speedup is slightly larger (4-6 here, 3-4 before). Since the computational cost of the SDC iterations with algebraic adaptivity We tested both unpreconditioned CG and BDDC-preconditioned CG for solving the full systems. As expected, BDDC reduces the iteration count dramatically, and leads to resolutionindependent convergence due to the bounded condition number of the preconditioned system, see Tab. 4. Despite the increased computational effort per iteration introduced by the preconditioner, the computational time is also reduced by a large factor.

Conclusion

We have made significant progress towards achieving efficient large-scale EMI model simulations, focusing on three key aspects: (i) robust creation of artificial but structurally realistic myocardium geometries and grids, (ii) efficient spatial adaptivity for higher-order time stepping, achieved through a combination of algebraic adaptivity and ladder SDC methods, and (iii) development of BDDC preconditioners tailored for complex 3D piecewise continuous elliptic problems.

The numerical results demonstrate that the proposed methods are effective, efficient, and allow scaling EMI computations to larger domain sizes than previously attainable.

Figure 1 :

 1 Figure 1: Overview of the mesh generation methods. Left: a network with 6 centers (black dots), 5 links, and halflink domains. The dashed lines represent cell membranes. Middle and right: transverse and longitudinal section of the level-set function, where red colors code the inside of the cells. The narrow spaces between the cells (about one tenth of the cell diameter in reality) made it necessary to discretize the domain boundaries before computation of the level-set function.

Figure 2 :

 2 Figure2: Each myocyte (orange) is surrounded by an extracellular subdomain (green). As illustrated here for a location near a model boundary, the extracellular domains can be considerably larger than the cells when these are sparse.

Figure 3 :

 3 Figure 3: Expected speedup due to ladder methods applied in addition to algebraic adaptivity for different numbers n of collocation points and larger contraction factors ρ q .

F

  . Chegini, A. Froehly, N.M.M. Huynh, L. Pavarino, M. Potse, S. Scacchi, and M. Weiser

Figure 4 :

 4 Figure 4: Different 3D meshes, detail in Table 2. Dimensions are indicated in micrometers; the largest model (d) is 1 mm long. Individual myocytes are differentiated by color.

Figure 5 :

 5 Figure 5: Number of degrees of freedom of the EMI model on grid (a) from Tab. 2 used in the SDC iterations versus simulated time. Left: pure algebraic adaptivity. Right: algebraic adaptivity combined with ladder methods.

Table 2 :

 2 Properties of the meshes shown in Fig.4.

  is no longer uniform, ladder methods yield a further moderate but consistent increase of the speedup by 20-50%, which is well within the range predicted by the quite coarse efficiency model worked out in Sec. 4.2. The number of retained dofs, however, exhibits a strong oscillation, which suggests that the combination of algebraic adaptivity and ladder methods deserves more investigation.

	mesh AA AA + Ladder ratio
	(a) 4.526	5.342	1.180
	(b) 4.753	5.973	1.257
	(c) 5.432	8.252	1.519
	(c) 5.863	8.847	1.509

Table 3 :

 3 Computational speedup due to algebraic adaptivity (AA) applied to standard SDC and to ladder SDC methods, over non-adaptive SDC. With 3 collocation points used, the additional speedup due to ladder methods is within the range expected from the coarse efficiency model (4).

Table 4 :

 4 Iteration count and condition number(k 2 ) estimate for unpreconditioned and BDDC-preconditioned CG on problems of different size.

Acknowledgement. This work was supported by the European High-Performance Computing Joint Undertaking EuroHPC under grant agreement No 955495 (MICROCARD) co-funded by the Horizon 2020 programme of the European Union (EU), the French National Research Agency ANR, the German Federal Ministry of Education and Research, and the Italian ministry of economic development.


Joshua Steyer

Fatemeh Chegini

Mark Potse

Axel Loewe

Martin Weiser

Continuity of Microscopic Cardiac Conduction in a Computational Cell-by-Cell Model

Conduction velocity in cardiac tissue is a crucial electrophysiological parameter for arrhythmia vulnerability. Pathologically reduced conduction velocity facilitates arrhythmogenesis because such conduction velocities decrease the wavelength with which re-entry may occur. Computational studies on CV and how it changes regionally in models at spatial scales multiple times larger than actual cardiac cells exist. However, microscopic conduction within cells and between them have been studied less in simulations. In this work, we study the relation of microscopic conduction patterns and clinically observable macroscopic conduction using an extracellularmembrane-intracellular model which represents cardiac tissue with these subdomains at subcellular resolution. By considering cell arrangement and non-uniform gap junction distribution, it yields anisotropic excitation propagation. This novel kind of model can for example be used to understand how discontinuous conduction on the microscopic level affects fractionation of electrograms in healthy and fibrotic tissue. Along the membrane of a cell, we observed a continuously propagating activation wavefront. When transitioning from one cell to the neighbouring one, jumps in local activation times occurred, which led to lower global conduction velocities than locally within each cell.

Introduction

Being one of the main reasons for sudden cardiac death, which causes around 700,000 deaths per year in Europe alone, cardiac arrhythmias remain a major public health issue [START_REF] Chugh | Sudden cardiac death in 2017: Spotlight on prediction and prevention[END_REF]. An important tissue property affecting arrhythmia vulnerability is the conduction velocity (CV) of cardiac excitation waves, which can be derived from local activation times (LATs), defined as the moment at which the excitation wave front reaches a certain point in space. The spatial distribution [START_REF] Verma | Regional conduction velocity calculation from clinical multichannel electrograms in human atria[END_REF] and effect [START_REF] Azzolin | Personalized ablation vs. conventional ablation strategies to terminate atrial fibrillation and prevent recurrence[END_REF] of these two quantities have been investigated in homogenised models in which a model element averages over several hundreds of myocytes. Such models, however, do not directly consider myocyte size and shape, an inhomogeneous arrangement of them, connectivity patterns, and microscopic conduction from cell to cell which has been found to play a role in arrhythmogenesis [START_REF] De Bakker | Ventricular tachyrdia in the infarcted, Langendorff-perfused human heart: Role of the arrangement of surviving cardiac fibers[END_REF]. In this study, we investigate the distribution of LAT and CV values using the novel extracellular-membraneintracellular (EMI) model, which explicitly represents these three domains in the mesh [START_REF] Tveito | Modeling Excitable Tissue: The EMI Framework[END_REF][START_REF] Jaeger | From Millimeters to Micrometers; Re-introducing Myocytes in Models of Cardiac Electrophysiology[END_REF]. To do so, we used a 3D myocyte configuration embedded in a bath for which excitation propagation is initiated for three different cases to study purely longitudinal and transversal wave propagation, as well as a combination of them. We then determined LATs both on the membrane and the intracellular space and derived local CV, which were both analysed in terms of difference between the local (subcellular and cellto-cell) and global behaviour of the excitation propagation.

Methods

We used the EMI model described in [START_REF] Tveito | Modeling Excitable Tissue: The EMI Framework[END_REF], which represents the extracellular medium Ω e , intercalated discs Γ g , the cell-to-extracellular membrane Γ m and the intracellular domain Ω i = j Ω j consisting of individual myocyte domains, Ω j , explicitly in the mesh. Its model equations read

               ∇ • (σ i ∇ϕ i ) = 0 in Ω i , ∇ • (σ e ∇ϕ e ) = 0 in Ω e , V m = ϕ i -ϕ e on Γ m , w g = ϕ j -ϕ k where j ̸ = k on Γ g , C m Vm + I ion = -n T i σ i ∇ϕ i = n e σ e ∇ϕ e on Γ m , C m ẇg + I g = -n T j σ i ∇ϕ j = n k σ i ∇ϕ k on Γ g , (1)
where ϕ denotes the potential in the corresponding domain, while V m describes the membrane voltage between Ω i and Ω e , and w g between two adjacent cell domains, Ω j and Ω k . n denotes the unit outer normal of the respective domain and the conductivities of the intra-and extracellular domains are σ i = 0.3 S/m and σ e = 2.0 S/m, and the membrane capacitance is C m = 10 -4 F/m 2 . To model membrane dynamics, we used the Aliev-Panfilov model [START_REF] Aliev | A simple two-variable model of cardiac excitation[END_REF] with its ion current given by

Computing in

I ion (V m , w g ) = g a V m (V m -a)(V m -1) + V m w g , (2)
and the gating dynamics being governed by

∂w g ∂t = - 1 4 ε(V m , w g )(w g + g s V m (V m -a -1)), (3)
where

ε(V m , w g ) = ε 0 +µ 1 w g /(V m +µ 2 ), g a = g s = 8.0, a = 0.1, ε 0 = 0.01, µ 1 = 0.07 and µ 2 = 0.3.
Being purely phenomenological, the normalised state variables of the dynamics V m and w g are restricted to [0,1] a.u. Similar to [START_REF] Chegini | Efficient adaptivity for simulating cardiac electrophysiology with spectral deferred correction methods[END_REF], we used a linear (Ohmic) relationship between the current through the gap junctions and the respective voltage between the two cells:

I g = w g R g , (4)
where R g = 4.5•10 -4 Ωm 2 . This system was solved using the Kaskade 7 finite element toolbox [START_REF] Weiser | Higher-order time integration using spectral deferred correction method (SDC) in a cell by cell discretization of cardiac excitation[END_REF]. The 3D mesh [START_REF] Potse | A practical algorithm to build geometric models of cardiac muscle structure[END_REF] is shown in The myocytes' largest extension (i.e. preferential orientation) is along the x-axis, which is thus the direction of fastest propagation and may be referred to as longitudinal direction. For transversal direction, we consider the y-axis. The intercalated discs with the gap junctions are located at x = 150, 250 µm along the longitudinal axis, such that we have two planes (y, z) with gap junction properties given by Eq. 4.

In total, we ran three different simulations with different initial conditions leading to different excitation dynamics: a) initial depolarisation of all leftmost myocytes along the x-axis (causing an almost purely longitudinal propagation and referred to as the left case), b) depolarisation of all lowermost cells (purely transversal propagation, bottom case) and c) depolarisation of one cell in the corner of the myocyte configuration (longitudinal and transversal propagation, single case).

To quantify the spatial propagation of the excitation along the cell membrane (quantified via V m (x, t)) and the intracellular potential (ϕ i (x, t)), where x ∈ R 3 denotes the position, we chose V th m = 0.5 a.u and ϕ th i = 0.45 a.u. as thresholds to determine the local activation times (LATs). Using a projection of these LAT values onto the primary axis of propagation, i.e. either x or y, we then defined longitudinal (∥) and transversal (⊥) CVs as:

CV ∥ (x) = ∂LAT(x) ∂x -1 , (5)
CV ⊥ (x) = ∂LAT(x) ∂y

-1 . ( 6
)

Results

As already mentioned, LATs are obtained with the thresholds imposed on the time evolution of V m and ϕ i . An exemplary snapshot of ϕ i (x, t) is given in Fig. 2 for a purely longitudinal propagation. By determining the nearest neighbour nodes, we compared the LAT values obtained via V m and ϕ i . In all of the three simulations, the Pearson correlation was always ρ > 0.999, which is why we here look at both LAT datasets combined. The projections of these LAT values onto these axes are given in Fig. 3. For all cases, we fitted a linear function to a) the data within each "column" of cells (perpendicular to the propagation direction, green lines) and b) the complete data set (red line). The upper left panel shows the projection of LAT values of the central and right myocyte columns onto the x-axis (i.e. the axis of propagation). Both values grow linearly in propagation direction distance with a jump at the intercalated disc (x = 250 µm). The upper right panel illustrates the distribution of LAT values along the transversal axis obtained from the "bottom" case. Note that, in contrast to the fully longitudinal propagation, the separate myocyte planes overlap, since the cells are intertwined both along the y and z direction (Fig. 1). Finally, for the depolarisation of one cell in the corner, we can project the LATs onto both the longitudinal and transversal direction, given in the bottom row of Fig. 3, respectively. While the LATs in transversal direction show linear increase in propagation direction, both within the myocyte planes and globally, the longitudinal data for this initialisation does not show a clear trend and is broadly distributed. We now quantify the time delay of the activation propagation when transitioning from one myocyte to the next via the gap junctions in the intercalated disc. Both the mean and standard deviation of the LAT delays, ∆LAT, are given in Tab. 1.

propagation

∆LAT (ms) fully long. 0.8 fully trans.

1.6 ± 0.9 long. (single) -0.1 ± 0.3 trans. (single)

1.5 ± 0.7

Table 1: Mean and standard deviation of LAT delays on the membrane during transition from one cell layer to the neighbouring one.

The local and global CVs defined as the inverse slopes of the single myocyte column regressions and the inverse of the global slope, according to Eq. 6, are given in Tab. 2. Overall, the global CVs are smaller than the locally determined ones.

Discussion

We have analysed local CV values in the EMI model [START_REF] Tveito | Modeling Excitable Tissue: The EMI Framework[END_REF][START_REF] Jaeger | From Millimeters to Micrometers; Re-introducing Myocytes in Models of Cardiac Electrophysiology[END_REF], which explicitly represents intra-and extracellular media, as well as the membrane separately in the mesh. As opposed to what would be possible with homogenised models, such as the frequently used monodomain or bidomain models [START_REF] Plank | The openCARP simulation environment for cardiac electrophysiology[END_REF], this allowed studying the the local behaviour of the potential propagation. We have shown that LATs undergo jumps when the wave transitions from one myocyte to the next. CVs, calculated by the inverse slope of the linearly fitted LAT values reveal that the potential and excitation propagation within a cell is faster than globally, which in turn can be explained by the time delay in the order of a few milliseconds caused at each cell-to-cell transition. This is in accordance with what has been theorised e.g. in [START_REF] Keener | Mathematical Physiology: I: Cellular Physiology[END_REF] and shown in 2D simulations using the EMI model [START_REF] Jaeger | Efficient, cell-based simulations of cardiac electrophysiology; The Kirchhoff Network Model (KNM)[END_REF]. The absolute value of global and local CV is a result of the choice of model parameters and gap junction distribution [START_REF] De Souza | Effect of Gap Junction Distribution, Size, and Shape on the Conduction Velocity in a Cell-by-Cell Model for Electrophysiology[END_REF]. The setup used here motivated by previous work led to CV values one to two orders of magnitude smaller than typically observed in healthy human tissue [START_REF] Verma | Regional conduction velocity calculation from clinical multichannel electrograms in human atria[END_REF][START_REF] Lubrecht | Automatic reconstruction of the left atrium activation from sparse intracardiac contact recordings by inverse estimate of fibre structure and anisotropic conduction in a patient-specific model[END_REF]. There are several limitations to this analysis. First, the method to calculate CVs used here is not suitable for the longitudinal propagation in the single cell initialisation case due to multimodal LATs, as can be seen in Fig. 3, and will have to be done node-wise. Furthermore, we assumed the intracellular medium to be homogeneous without considering potential conduction anisotropy due, to e.g., the nucleus or other organelles in the cell. The findings presented here provide a first insight into electrophysiological excitation propagation behaviour (CV) on a microscale quantitatively. If the model parameters are adjusted such that physiological CV values are obtained, it may form a basis for future work analyzing for example the genesis of intracardiac electrograms in such microstructural models [START_REF] Steyer | Simulated Excitation Patterns in the Atria and Their Corresponding Electrograms[END_REF] or extending the setup to pathological scenarios allowing for mechanistical studies of the role of microscopic cell arrangement and distribution of, e.g., fibrosis.

  Cardiology 2023; Vol 50 Page 1 ISSN: 2325-887X DOI: 10.22489/CinC.2023.385

Fig. 1 .

 1 Its dimensions are 400 µm × 240 µm × 240 µm and it consists of 507648 tetrahedral mesh elements. The intracellular domain, Ω i comprises n x × n y × n z = 7 × 3 × 7 = 147 myocytes.

Figure 1 :

 1 Figure 1: Myocyte arrangement used for the simulations embedded in Ω e (semi-transparent grey box) and an exemplary single myocyte (top) with edges indicated in white. Note that colours are used to discriminate between individual myocytes, not varying electrophysiological properties.

Figure 2 :

 2 Figure 2: Exemplary snapshot of a purely longitudinal propagation of ϕ i (t) initialised by depolarising all leftmost myocytes.

Figure 3 :

 3 Figure 3: Spatial distribution of LAT values on the membranes as projections onto the axis of propagation. Green lines indicate the linear regression of LAT values within one cell layer (see Fig. 2) and the red line the linear regression of all LATs along this axis. Dark grey background indicates intercalated disc regions, which are wider in the right column due to intertwining of myocytes along the yaxis. Bottom row: Diagonal excitation starting from a single myocyte in the bottom left corner of the mesh.

Table 2 :

 2 Mean and standard deviation of the single myocyte layer LAT slopes and the global one along propagation direction x i .

Acknowledgments

This work was supported by the European High-Performance Computing Joint Undertaking EuroHPC under grant agreement No 955495 (MICROCARD) cofunded by the Horizon 2020 programme of the European Union (EU), the German Federal Ministry of Education and Research and the French National Research Agency ANR.


MIP-based heuristics for combinatorial design of reconfigurable rotary transfer machines for production of multiple parts

Keywords: Machining economics, preliminary machine design, combinatorial optimization, batch machining, line balancing, reconfigurable manufacturing systems

This paper deals with a problem of the optimal configuration of a rotary transfer machine with turrets for machining multiple parts. This is a hard combinatorial optimization problem appearing at the preliminary machine design stage. Such machines are multi-positional, i.e. parts are sequentially machined on several working positions. At each working position, several machining modules (spindle heads) can be installed to process the tasks assigned to this position. Machining modules are activated sequentially or simultaneously. Sequential activation is realized by the use of turrets. Simultaneous activation is possible if machining modules are related to the different sides of the part, and if they can work in parallel. There are horizontal and vertical spindle heads, and turrets to access different sides of the parts on a working position. At the preliminary design stage, the following decisions must be made: the choice of orientations of the parts on the rotary table; the partitioning of the given set of tasks into positions and their assignment to machining modules (selection or design of machining modules to use), and the choice of cutting modes for each spindle head and turret. The objective is to minimize the total cost of equipment used. The number of possible solutions for this combinatorial design problem increases exponentially with the number of part types to be produced, and this represents a computational burden for decision-makers (usually process engineers). In this paper, in order to help decision makers deal efficiently with the manufacturing of multiple batches of parts, we develop a powerful heuristic framework which can be used in real life industrial cases. We test the developed methodology on the real-life cases provided by one of our industrial partners and demonstrate its efficiency. The proposed model and algorithms allow to minimize the cost of designed machines.

Introduction

The use of decision aid tools to help process engineers design efficient machining systems is an important issue in production economics [START_REF] Hitomi | Integrated component scheduling models for chips hooter machines[END_REF]. The first work on such problems was initiated for mass production systems only machining one type of part. The solution techniques developed were based on graph theory, mathematical programming and the application of other solution approaches primarily employed for assembly line balancing [START_REF] Boysen | Assembly line balancing: Which model to use when?[END_REF][START_REF] Battaïa | A taxonomy of line balancing problems and their solution approaches[END_REF][START_REF] Battaïa | Hybridizations in line balancing problems: A comprehensive review on new trends and formulations[END_REF]. Examples of such models can be found in [START_REF] Szadkowski | An approach to machining process optimization[END_REF][START_REF] Szadkowski | Critical path concept for multi-tool cutting processes optimization[END_REF][START_REF] Dolgui | On problem of optimal design of transfer lines with parallel and sequential operations[END_REF][START_REF] Dolgui | Graph approach for optimal design of transfer machine with rotary table[END_REF].

In the current manufacturing context, in order to respond to volatile markets, manufacturing systems need to increase their flexibility and reconfigurability. These new circumstances make the design problem even more challenging and complex than in the case of mass production systems, but provide new opportunities for the development of competitive advantages if designed appropriately. A comparison of the profitability of a reconfigurable manufacturing system and several dedicated manufacturing systems, each of them used for a given product, was carried out in [START_REF] Dolgui | Profitability of a multi-model manufacturing line versus multiple dedicated lines[END_REF]. We consider in this paper the case of a reconfigurable machine with a rotary table. According to [START_REF] Cesar | Achieving plant responsiveness from reconfigurable technology: Intervening role of SCM[END_REF], Reconfigurable Manufacturing Systems (RMS) represent the link between reconfigurability, strategy and technology in order to support frequent changes and evolutions in markets. In fact, the reconfigurability of manufacturing systems plays a primordial role in the responsiveness of the Supply Chain [START_REF] Koren | Value creation through design for scalability of reconfigurable manufacturing systems[END_REF]. Since its introduction by [START_REF] Koren | Reconfigurable Manufacturing Systems[END_REF], the concept of RMS has gained a lot of attention both in academia and in practice [START_REF] Bortolini | Reconfigurable manufacturing systems: Literature review and research trend[END_REF][START_REF] Yelles-Chaouche | Reconfigurable manufacturing systems from an optimisation perspective: a focused review of literature[END_REF]Pansare et al., 2021aPansare et al., , 2021b;;[START_REF] Gu | Mass-Individualisationthe twenty first century manufacturing paradigm[END_REF].

Different families of optimization approaches have been developed in the literature in order to support decision makers in design, planning and operation of RMS. These include multicriteria decision making [START_REF] Goyal | Optimal configuration selection for reconfigurable manufacturing system using NSGA II and TOPSIS[END_REF][START_REF] Chaube | An adapted NSGA-2 algorithm based dynamic process plan generation for a reconfigurable manufacturing system[END_REF][START_REF] Khezri | Towards a sustainable reconfigurable manufacturing system (SRMS): multi-objective based approaches for process plan generation problem[END_REF], fuzzy techniques [START_REF] Singh | Justification for the selection of a reconfigurable manufacturing system: a fuzzy analytical hierarchy based approach[END_REF][START_REF] Abdi | Fuzzy multi-criteria decision model for evaluating reconfigurable machines[END_REF], artificial intelligence algorithms [START_REF] Li | Degradation-aware decision making in reconfigurable manufacturing systems[END_REF], optimization methods [START_REF] Youssef | Modelling and optimization of multiple-aspect RMS configurations[END_REF][START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF] including metaheuristics [START_REF] Dou | A multi-objective particle swarm optimisation for integrated configuration design and scheduling in reconfigurable manufacturing system[END_REF][START_REF] Khettabi | Sustainable reconfigurable manufacturing system design using adapted multi-objective evolutionary-based approaches[END_REF].

To achieve efficiency in the implementation of RMS, several optimization problems must be properly addressed starting with the preliminary design of the system [START_REF] Battaïa | Variety-oriented design of rotary production systems[END_REF][START_REF] Battaïa | Optimal cost design of flow lines with reconfigurable machines for batch production[END_REF][START_REF] Battaïa | Design of reconfigurable machining lines: A novel comprehensive optimisation method[END_REF][START_REF] Goyal | Design of reconfigurable flow lines using MOPSO and maximum deviation theory[END_REF], the choice of its layout (Haddou [START_REF] Benderbal | Machine layout design problem under product family evolution in reconfigurable manufacturing environment: a two phase-based AMOSA approach[END_REF], configuration management [START_REF] Moghaddam | Configuration design in scalable reconfigurable manufacturing systems (RMS): a case of single-product flow line (SPFL)[END_REF][START_REF] Ashraf | Configuration selection for a reconfigurable manufacturing flow line involving part production with operation constraints[END_REF][START_REF] Wang | A tree-based decision method for the configuration design of reconfigurable machine tools[END_REF][START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF][START_REF] Singh | Composite performance metric for product flow configuration selection of reconfigurable manufacturing system (RMS)[END_REF], process planning (Battaïa et al., 2017a;[START_REF] Xia | Reconfigurable machining process planning for part variety in new manufacturing paradigms: Definitions, models and framework[END_REF][START_REF] Khezri | Towards a sustainable reconfigurable manufacturing system (SRMS): multi-objective based approaches for process plan generation problem[END_REF], setup planning and scheduling [START_REF] Dou | Bi-objective optimization of integrating configuration generation and scheduling for reconfigurable flow lines using NSGA-II[END_REF]. The more integrated is the tool developed, considering many problems and system parameters at the same time, the more efficient are the results obtained, from a resource use point of view. [START_REF] Bensmaine | A new heuristic for integrated process planning and scheduling in reconfigurable manufacturing systems[END_REF] considered integrated process planning and scheduling in RMS. [START_REF] Dou | Bi-objective optimization of integrating configuration generation and scheduling for reconfigurable flow lines using NSGA-II[END_REF]Dou et al. ( , 2012) ) dealt with integrated configuration generation and scheduling for reconfigurable flow lines. Battaïa et al. (2017a) studied an integrated process planning and system configuration problem.

In recent studies [START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF][START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF], more attention is paid to the impact of the use of RMS on the sustainability of production systems. [START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF] developed a heuristic-based, non-linear mixed integer approach for minimizing energy consumption by selecting the most suitable modular machines from a set of candidate machines in a sustainable reconfigurable manufacturing environment. [START_REF] Massimi | A heuristic-based non-linear mixed integer approach for optimizing modularity and integrability in a sustainable reconfigurable manufacturing environment[END_REF] studied a multiobjective integer linear programming with three criteria: the total production time, cost and the sustainability-metric value related to liquid hazardous waste and greenhouse gas emissions. Moreover, the last objective was spit in two, namely: the amount of the greenhouse gases emitted by machines and the hazardous liquid wastes. The first two objectives were also used by [START_REF] Khettabi | Sustainable reconfigurable manufacturing system design using adapted multi-objective evolutionary-based approaches[END_REF] who developed a non-linear multi-objective integer program for which several evolutionary approaches have been developed. The reconfigurability of the manufacturing systems significantly increases their complexity and the number of possible manufacturing system configurations [START_REF] Farid | Measures of reconfigurability and its key characteristics in intelligent manufacturing systems[END_REF]. Therefore, the optimisation of such systems raises new challenges for system designers, process planners and production managers [START_REF] Battaïa | Design of reconfigurable machining lines: A novel comprehensive optimisation method[END_REF]. This paper yields a methodological support for decision makers in order to provide them with fast and efficient solutions for the preliminary design of reconfigurable rotary transfer machines with turrets for machining multiple parts (see Fig. 1). The number of possible configurations of such a machine increases exponentially with the number of different parts to be produced. The study is the first attempt in the literature to model the considered problem mathematically as a combinatorial optimization problem. Due to the solution methods developed in this article, the designers can deal efficiently with the computational burden.

On rotary machines, parts are sequentially machined on m (1, 2, …, m) working positions. The zero position of the machine is exclusively used for loading and unloading operations. At other positions, several machining modules (spindle heads) can be installed, with a sequential or parallel intervention on the parts. Turrets are used for sequential activation of spindle heads.

Parallel intervention is possible on the part, since the horizontal and vertical spindle heads (separate or mounted in turrets) can access its different sides (Battaïa et al., 2017ab).

In our study, motivated by a long-time collaboration with a company designing machining equipment, designers can use only one vertical turret or vertical spindle head mounted at a single position but capable of covering all working positions. Several horizontal spindle heads or turrets can be used, but only one per position. Fig. 1 shows an example of design solution with 6 working positions, where a common vertical spindle head is used on positions 1, 3, 4, 5, two horizontal turrets are installed on position 1 and 3, and one horizontal spindle head is used on position 4. At the preliminary design stage, the designers have to decide:

-how to place each part on the rotary table, -which tasks to execute for each part at each working position, -the spindle heads and turrets to be installed and the cutting modes for them. The quality of these decisions impacts the cost and productivity of the machine. In this paper, we develop a methodological support to help designers to take design decisions at this step which are as close as possible to optimum.

In this paper, integrated process (equipment) design and task assignment is considered for the preliminary design of a reconfigurable rotary machine. Such machines are widely used in automotive, mechanical equipment and airspace industries. The machine is configured for batch production of different parts. During our long-term collaboration with our industrial partner designing rotary machines, we have worked on several related optimization problems in the past. Initially, we studied the design of mass production rotary machines for single product manufacturing [START_REF] Dolgui | Graph approach for optimal design of transfer machine with rotary table[END_REF]. Then, the variety-oriented design of rotary production systems was investigated in [START_REF] Battaïa | Variety-oriented design of rotary production systems[END_REF]. Later on, we considered the integrated process planning and system configuration for such machines (Battaïa et al., 2017ab) for the particular case of mixed-model machining. Contrary to all these previous studies, in this paper, we deal Vertical spindle head is used on positions 1, 3, 4, 5 with batch production on reconfigurable rotary machines. To the best of our knowledge, this is the first attempt in the literature to address this problem.

Previously, the design of multi-part machining lines was addressed in [START_REF] Dolgui | Minimizing setup costs in a transfer line design problem with sequential operation processing[END_REF], but for only for sequentially executed tasks and with need of setup between different parts. In this case, the main optimization objective was to minimize the number of stations and the second objective was to minimize the setup cost. As a consequence, the developed model is not suitable for the setting considered in this paper. We can also mention studies on assembly line design for Printed Circuit Boards (PCB) which have some similarities with the machining lines considered in this article [START_REF] Hitomi | Integrated component scheduling models for chips hooter machines[END_REF][START_REF] Sun | Component allocation and feeder arrangement for a dual-gantry multihead surface mounting placement tool[END_REF]. Nevertheless, the differences in technologies, the induced technical constraints and parameters make those models inapplicable for the combinatorial design of reconfigurable rotary transfer machines for production of multiple parts.

In Section 2, we define the general optimization framework. In Section 3, we develop the methodological support based on mathematical models and heuristic algorithms. In Section 4, the results of case studies are reported and discussed. Concluding remarks and research perspectives are given in Section 5.

Problem statement

We consider the problem of combinatorial design and task assignment of a rotary transfer machine with m working positions for machining d0 types of parts with required output Od, d =1, 2, …, d0.

In practice, batch production can be organized on rotary machines in the following three modes:

-different parts of a batch are machined sequentially (problem A1); -different parts of a batch are machined simultaneously (problem A2); -different batches of parts are machined sequentially (problem A3). In the first case, the rotary transfer machine is reconfigured once the processing of Od parts of type d is finished, i.e. the fixtures of parts are changed and some spindles are mounted or dismounted if necessary.

In the second case, the parts are loaded in sequence =(1, 2, …, 0 ), where i{0, 1, 2, …, d0}, i=1, 2, …, 0, 0 is multiple to m0+1 (i=0 means that no part is loaded). We define m0 as a maximum possible value of m which depends on the type of the rotary machine to be designed. In this case, no reconfiguration of the rotary transfer machine is allowed during the processing of a batch, but can be done when batch mix changes. We assume O d =  d , d=1, 2, …, d0, where  d is the occurrence number of d in . Using sequence , we can define in a oneto-one manner function (i,k) of part number on the k-th working position.

In the third case, the parts are grouped in  batches with required output O, =1,2,…, , which are processed sequentially. Parts of -th batch are loaded in sequence =(1, 2, …,  ), where j{0, 1, 2, …, d0}, j=1, 2, …, ,  is multiple of m0+1, (i=0 means that no part is loaded). Using sequences , we can define in a one-to-one manner function (i,k), i=1,…,O+m0-1, of part number on the k-th working position after i turns of the rotary table.

In this case,

O d =    =    d O 1 , d=1, 2, …, d0
, where  d is the occurrence number of d in .

After the end of a current batch, the rotary transfer machine is reconfigured, i.e. the fixtures of parts are changed and some spindles are mounted or dismounted if necessary.

Let N d be the set of machining operations (tasks) needed for machining of elements of the d-th part d=1, 2, …, d0, located on nd sides. )) ( (

= =
, where hrs(d) equals j, j=1,2, if the elements of the s-th of the part d can be machined by spindle head or turret of type j. The choice of the part orientation defines which sides of the part will be accessible for horizontal spindle heads, and which only for vertical ones.

There is no solution if  d s N p p H   = ) (
for some d{1, 2, …, d0} and s{1, 2, …, nd}.

Let subset Nk, k=1,...,m contain the operations from set N assigned to the k-th working position. Let sets Nk1 and Nk2 be the sets of operations assigned to working position k that are concerned by vertical and horizontal machining, respectively. Finally, let bkj be the number of machining modules (not more than b0) of type j installed at the k-th working position and, respectively, subsets Nkjl, l=1,...,bkj contain the operations from set Nkj assigned to the same machining module. This assignment has to respect the technological constraints of the machining process required. These can be grouped in the three following families.

A number of known technological factors (such as fixed sequences of operations for machining part elements, the presence of roughing, semi-finishing and finishing operations, etc.) determines an order relation on the set N, which defines possible sequences of operations.

These precedence constraints can be specified by a directed graph G OR =(N,D OR ), where an arc (p,p)D OR if and only if the operation p has to be executed before the operation p. It should be noted that if such operations p and p belong to different sides of the part then they cannot be executed at the same position without violating the precedence constraint.

The required precision (tolerance) of mutual disposition of machined part elements as well as a number of additional factors imply the necessity of performing some pairs of operations from N at the same working position, by the same turret, by the same spindle head or even by the same spindle (for different parts) for each pair. Such inclusion constraints can be given by undirected graphs

G SP =(N,E SP ), G ST =(N,E ST ), G SM =(N,E SM ) and G SS =(N,E SS ),
' ' d s N , s{1,…,nd}, q " " d s N , s{1,…,nd} such that (p,q)  E SS  E SM  E ST and hrs(d)hrs(d).
t b (Pdkjl)=L(Ndkjl)/Гdkjl+ a , (1)
where L(Ndkjl)=max{(p)|pNdkjl}, and  a is the time for advance and disengagement of tools.

We assume that if the turret of type j is installed at the k-th position, then the execution time of operations from Ndkjl is equal to (2):

t h (Pdkj)= g bkj +  = kj b l 1 t b (Pdkjl), |j=1, 2, (2)
where  g is the additional time for one rotation of turret. If the spindle head is installed, then

t h (Pdkj)= t b (Pdkjl), |j=1,2. If all Ndkjl are empty, then t h (Pdkj)=0. If bkj =1 then t h (Pdkj)=t b (Pdkj1).
The execution time t p (Pdk) is defined as in (3):

t p (Pdk)= r +max{t h (Pdkj)|j=1,2}, (3)
where  r is the additional time for table rotation.

When different parts of a batch are machined sequentially, then the time td for machining all the elements of the d-th part is equal to (4):

t d (P)=max{t p (Pdk)| k=1,…,m}. ( 4
)
Taking into account that, at the beginning and the end of machining of O d parts, not all the working positions are occupied, the total time T(P) for machining O d parts is equal to (5):

T(P) =  = - + 0 1 ) 1 )( ( d d d d m O P t . (5)
When different parts of a batch are machined simultaneously, then time T(P) of execution of all corresponding operations after 0 turns of rotary table is defined as in (6):

T(P)=   = 0 1 i max{ ) ( ) , ( k k i p P t  |k=1,,m0}. (6)
When different batches of parts are machined sequentially, then the time T(P) of execution of all corresponding operations after  turns of rotary table is defined as in (7):

T(P)=  - +  =   1 1 0 m O i max{ ) ( ) , ( k k i p P t   |k=1,,m0}, (7)
and the time T(P) for machining all the batches is equal to (8):

T(P)=   =  1 T(P). ( 8
)
We assume that the given throughput is provided, if the total time T(P) does not exceed the available time T0. It is easy to see that the constraint on the throughput is provided if and only if it is satisfied for Гdkjl=min{γ2(p)|pNdkjl}, d=1,…,d0, k=1,…,m, j=1, 2, l=1,…,bkj.

Let C1, C2, C3, and C4 be the relative costs for one position, one turret, one machining module of a turret, and one spindle head, respectively. Since the vertical spindle head (if it is installed) is common for several positions, its size (and therefore the cost) depends on the number of positions to be covered. Let h k min and h k max be the minimal and maximal positions of the common vertical spindle head. Then, its cost can be estimated as C4+( h k max -h k min )C5 , where C5 is the relative cost for covering one additional position by a vertical spindle head. If the vertical spindle turret is installed, its cost can be estimated by C2+C3bk1. In a similar way, the cost C(bk2) for performing a set of operations Nk2 by associated bk2 machining modules can be assessed as in (9):

C(bk2) =       + = = . 1 if , 1 if , 0 if 0 2 2 3 2 2 4 2 k k k k b b C C b C b (9)
The machine cost Q(P) is calculated as the total cost of all equipment used, i.e. ( 10):

+ + + - + =   = = ) )( | (| )) | (| 1 |)( (| ) ( 1 3 2 1 12 1 12 1 4 1 k m k k m k k b C C N sign N sign sign C m C P Q N  = + - m k k h h b C k k C 1 2 min max 5 ) ( ) ( , ( 10
)
where sign(a) = 1 if a > 0, and sign(a) = 0 if a  0.

Since m=m0 for problems A2 and A3, we let C1=0 in these cases. If the vertical turret is installed, then the second and forth summands are equal to 0, since Nk12≠ for some ,d0, k=1,...,m, l=1,...,bkj; e) the feed per minute Гdkjl for each subset Ndkjl,d=1,2,…,d0,k=1,...,m,j=1,2,l=1,,bkj in such a way that the machine cost is as small as possible, and none of the constraints are violated.

In the next section, we present mixed-integer programming (MIP) formulations of the considered design problems.

MIP formulations

We can build set N based on graph For A2: F i is an auxiliary variable which determines the time of execution of operations from N at all the positions when machining of part i is finished.

G SSM =(N,E SSM =E SS  E SM ). Let ) , ( SSM i SSM i SSM i E N G = , i=1,…,n SSM ,
For A3: F I is an auxiliary variable which determines the time of execution of operations from N at all the positions when processing of part i of -th family is finished; i s F  is an auxiliary variable for determining the time of execution of operations from N at the first i positions after the i-th turn of the rotary table for processing -th family; i f F  is an auxiliary variable for determining the time of execution of operations from N at the last i positions after the O+m0-i-th turn of the rotary table for processing -th family.

We calculate in advance parameters tpp = max((p), (p))/min(2(p),2(p))+ a . They represent the minimal time necessary for the execution of operations p and p in the same machining module. It is assumed that (p,p)E DM if min(γ2(p),γ2(p)) < max(γ1(p),γ1(p)). For each operation pN, we calculate a set B(p) of block indices from {1,2,…, 2m0b0} and a set K(p) of position indices from {1,2,…,m0}, where operation pN can be potentially assigned.

Let I(k)=[2(k-1)b0+1,2kb0], I(k,j)=[2(k-1)b0+(j-1)b0+1,2(k-1)b0+jb0],
and I(k,j,l)=[2(k-1)b0+(j-1)b0+l,2(k-1)b0+ (j-1)b0+l], respectively.

Cost calculation

The objective can be represented as in (11):

Min     = = = = + - + + + 0 0 0 1 2 1 2 4 3 2 1 21 4 1 1 ) 2 ( m k j kj m k k m k k Y C C C Y C Z C ) ( min 1 max 1 5 1 4 1 2 1 3 3 0 0 Y Y C Y C Y C m k j b l kjl - + +    = = = (11) 4.2 Assignment constraints 1 ) ( =   pq p B q X ; pN (12)   +        - = = q p p B j k I q q p p B j k I q qX qX k k j ) ( ) ( ) , ( ) ( ) ( ) ' , ' ( 1 1
' 2 1 '   , ) 1 ( ) ' ( ) ' ( ) ,
(

q p p B j k I q X q    -  (p,p)D OR ; p, pN; kK(p); j=1,2 (13) q p p B k I q X ) ( ) ( ) (     = ' ) ' ( ) ' ( ) ( ' q p p B k I q X     ; (p,p)E SP ; p, pN; kK(p)K(p) (14) q p p B j k I q X ) ( ) ( ) , (     = ' ) ' ( ) ' ( ) , ( '
q p p B j k I q X    
; (p,p)E ST ; p,pN; kK(p)K(p); j=1,2 (15) …,d0; s=1,…,nd; kK(p);

q p p B k I q X ) ( ) ( ) (     + ' ) ' ( ) ' ( ) ( ' q p p B k I q X     1, (p,p)E DP ; p, pN; kK(p)K(p) (16) q p p B j k I q X ) ( ) ( ) , (     + ' ) ' ( ) ' ( ) , ( ' q p p B j k I q X     +Ykj22, (p,p)E DT ; p, pN; kK(p)K(p); j=1,2 (17) q p X ) (  + q p X ) ' (  1 (p,p)E DB ; p, pN; qB(p)B(p) (18) q p X ) (  =  p d s N ; d=1,
{hrs(d)= j| r=1,…,rd}=; qI(k,j)B(p) (19) q p j k I p B q X p K k ) ( ) , ( ) ( ) (       = ' ) ' ( ) , ( ) ' ( ' ) ' ( q p j k I p B q X p K k       ; p, p d s N ; j=1,2; d=1,…,d0; s=1,…,nd (20) d kjl Y       ) ( ) , , ( ,
) ( p B l j k I q p q p d X N ; d=1,…,d0; k=1,…,m0; j=1,2; l=1,…,b0 (21)
     ) ( ) , , ( , )
( p B l j k I q p q p d X N |N d | d kjl Y ; d=1,…,d0; k=1,…,m0; j=1,2; l=1,…,b0 (22)
Ykjl   = 0 1 d d d kjl Y ; k=1,…,m0; j=1,2; l=1,…,b0 (23)
Yk11 + Yk212Zk; k=1,…,m0 (30)
(m0-k+1)Yk11+Y1min  m0+1; k=1,…,m0 (31)

Y1max  kYk11; k=1,…,m0 (32) Y1max  m0Y1 (33) Y1min  m0Y1 (34) ds kj Y       ) ( ) , ( ,
) ( p B j k I q N p q p d s X
; d=1,…,d0; s=1,…,nd; k=1,…,m0; j=1,2 (35) …,d0; s=1,…,nd; k=1,…,m0; j=1,2 (36)

     ) ( ) , ( , ) ( p B j k I q N p q p d s X | d s N | ds kj Y ; d=1,
 = d n s ds k Y 1 1  1; d=1,…,d0; k=1,…,m0 (37)
-  1 d r h ds kj h j j n s Y rs d    = = 2 , 1 1 ; d=1,…,d0; r=1,…,nd (38)  = d r r d r h 1 =1; d=1,…,d0 (39)
 e d r ) , ( d r h  |e|-1, eE DH , k=1,…,m0 (40)
Equations ( 12) provide the assignment of each operation from N exactly to one machining module. Equations ( 13) -( 17) are used to model precedence constraints, inclusion constraints for positions and turrets, as well as exclusion constraints for positions, turrets, and machining modules, respectively. Equations ( 18) prohibit the assignment of operations from

d s
N to machining modules of type j if there is no feasible orientation of part d for such an execution. Equations ( 19) guarantee the assignment of operations from d s N to the same type of spindle head or turret. Constraints ( 20) -( 25) define the existence of machining module l of type j at position k. Equations ( 26) -( 34) are used to calculate Zk, k=1,…,m0, Y1, as well as Y1min and Y1max. Constraints ( 35) -( 40) provide the choice of feasible orientation of each part d. j=1,2; d=1,…,d0; k=1,…,m0; l=1,…,b0; qI(k,j,l) pN d ; j=1,2; d=1,…,d0; k=1,…,m0; l=1,…,b0; qI(k,j,l)

Time calculation

d kjl F  tpp q p X ) (  ; pN d ;
B(p) (41) d kjl F  tpp( q p X ) (  + q p X ) ' (  -1); p,

B(p)B(p)

(42) …,d0; k=1,…,m0; j=1,2 (44) Expressions ( 41) -( 42), ( 43) and ( 44) are used to estimate the execution times of operations from N d by the l-th machining module, vertical spindle head and at the k-th position, respectively.

d k F 11  ((p)/γ2(p)+ a )( q p X ) (  + q p X ) ' (  -1); pN d ; pN; d=1,…,d0; k, k=1,…,m0; k k or p N d ; qI(k,1,1)B(p); qI(k,1,1)B(p) (43) d k F  ) 1 ( 2 0 1 3 2 0 0 -  +  +  +   = = d kj g b l b l kjl g kj g d kjl Y b Y Y F ; d=1,
d F  d k F + r ; d=1,…,d0; k=1,…,m0 (45) d k T  Fd-T0(1-Zk); d=1,…,d0; k=1,…,m0 (46) ) ( 0 0 1 1 d d d m k d k d d F T O F - +   = = T0 (47) i F  ) , ( k i k F  + r ; i=1,…,0; k=1,…,m0 (48) 0 1 0 T F i i    = . ( 49
) i F   ) , ( k i k F   + r ; =1,…,; i=1,…,; k=1,…,m0 (50) i s F   ) , ( k i k F   + r ; =1,…,; i=1,…,m0-1; k=1,…,i (51) i f F   ) , ( 0 k i m O k F - +      + r ; =1,…,; i=2,…,m0; k= i,,m0 (52)     =  =   1 1 ( i O i F  +  - = 1 1 0 m i i s F  +  = 0 2 m i i f F  -     + -  = 2 0 m i i F  )  T0 (53)
Expressions ( 45) -( 47), ( 48) -( 49), and ( 50) -( 53) provide the required productivity for problems A1, A2 and A3, respectively.

Bound constraints

d k F  r d k t  - k=1,…,m0; d=1,…,d0 (62) Fd  d t  d t  d=1,…,d0 (63) i F  i T  i T  i=1,…,0 (64) i F   i t   i t   = i=1,…, (65) i s F   } ,..., 1 | max{ ) , ( i k t k i =    } ,..., 1 | max{ ) , ( i k t k i k =    = i=1,…,m0-1 (66) i f F   } ,..., | max{ 0 ) , ( m i k t k i O = -       } ,..., | max{ 0 ) , ( m i k t k i O k = -       = i=2,…,m0 (67)
where, for problem A1 : d t =max{(p)/2(p)+ a + r |pN d };

d d d d d d d d k O O t T t / ) ( 0 ' , 1 ' ' ' 0  - =  = , for problem A2: } ,..., 1 | max{ 0 ) , ( m k t T k i i = =  ,    = - = 0 ' , 1 ' ' 0 i i i i i T T T and } ) , ( , ,..., 1 | max{ 0 d k i i T t i d k =   = = ,
and for problem A3:

} ,..., 1 | max{ 0 ) , 2 ( 0 m k t t k i m i = = + -    ,       =   =        =    - - = ' 1 ' ' , 1 ' ' ' ' ' , 1 ' ' 0 / ) ( i i i i i i i t O t O T t and } ) , 2 ( , ,..., 1 , ,..., 1 | max{ 0 d k i m i t t i d k = + -   =  =  =    .

Calculations of B(p) and K(p)

Based on matrices H(d) and the set E DH , we can build a matrix H of possible orientations of all parts. Each row of H defines in a one-to-one manner the orientation H and the corresponding partition of N to N1 and N2. For such a partition, we can calculate B(p,H) and K(p,H) by using algorithms from [START_REF] Dolgui | Graph approach for optimal design of transfer machine with rotary table[END_REF][START_REF] Dolgui | MIP approach to balancing transfer lines with blocks of parallel operations[END_REF]. Then, B(p) and K(p) are defined as

 H  = H H p B p B ) , ( ) (
and

 H  = H H p K p K ) , ( ) (
, respectively.

As has been mentioned, the modelled optimization problem is difficult to tackle as the number of parts to be machined increases. As a consequence, designers need heuristic approaches in order to deal with such problem instances. In the next section, we develop a new heuristic framework in order to provide designers with efficient solution methods for largescale problems.

Heuristics

In this section, we present a heuristic algorithm based on the idea of a heuristic approach developed for mass production machining lines by [START_REF] Guschinskaya | A heuristic multi-start decomposition approach for optimal design of serial machining lines[END_REF]. This approach is adapted to the case of combinatorial design of reconfigurable rotary transfer machines for production of multiple parts. The main approach used in this heuristic framework is to fix random sets of decision variables and optimize the values of remaining decision variables. This process is repeated during available solution time.

Algorithm 1:

Step 0. Generate the initial population using a greedy heuristic or truncated resolution with a MIP solver.

Step 1. Let Cmin = , TRtot = 0, TRnimp = 0.

Step 2. Choose two solutions X and X from the population.

Step 3. Let N f ={pN|q * (p,X)=q * (p,X)} where 1) , ( * = X p pq X .

Step

4. If |N\N f | < Nmin a random set is removed from N f . If |N\N f | > Nmax a random set is added to N f .
Step 5. Let B(p)={q * (p,X)} for pN f and B(p) be unchanged for pN\N f .

Step 6. Solve the obtained problem ( 1) -( 57) with a MIP solver within time limit TMIP.

Step 7. Compute the value C of the objective function (1) for the obtained solution.

Step 8. Replace the worst solution in the population with the new solution.

Step 9. With probability Pm let N f =N\N f and execute Steps 4 -8.

Step 10. If Cmin > C , then set Cmin = C, TRnimp = 0 and keep the current solution as the best, set TRnimp = TRnimp + 1, otherwise.

Step 11. Set TRtot = TRtot + 1.

Step 12. Stop if one of the following conditions holds:

• a given solution time is exceeded;

• TRtot is greater than the maximum number of iterations authorized;

• TRnimp is greater than a given value.

Go to

Step 2, otherwise.

The five versions HEUR1, HEUR2, HEUR3, HEUR4 and HEUR5 of the general heuristic framework were developed as follows:

• In HEUR1, model ( 10) -( 67) developed in this paper is used for all sub-problems where only the calculation of the bounds ( 54) -( 60) is adapted to each sub-problem.

• In HEUR2, model ( 10) -( 67) is rebuilt for each sub-problem.

• In HEUR3, all known solutions are grouped in clusters according to the same orientations of all parts to be machined. In this case, only solutions from the same cluster are chosen at Step 2 of Algorithm 1, B(p) are modified according to the orientations of parts and model ( 10) -( 67) is rebuilt.

• HEUR4 and HEUR5 are modified versions of HEUR2 and HEUR3, respectively, where at Step 5 of Algorithm 1, the assignment of any operation pN f can be attributed to any machining module at the corresponding working position.

In order to compare the performances of the developed heuristic methods, we use a number of case studies provided by our industrial partner. The results of this comparison are presented in the next section.

Experimental study

Small industrial case studies

The first dataset includes 25 industrial instances provided by our industrial partner. The characteristics of 25 instances are presented in Table 1. In this table, |N| is the number of operations, OSP is the order strength of precedence constraints, DM, DT, DP, SS, and SM are the densities of graphs G DM , G DT , G DP , G SS , and G SM respectively, LS is the sum of lengths of loading sequences, NB is the number of batches, and NO is the number of rows of matrix H. Experiments were carried out on an ASUS notebook (1.86 Ghz, 4 Gb RAM) with an academic version of CPLEX 12.2. In Tables 2 and3, we compare results for CPLEX12.2 (maximal solution time 3600 sec) and the proposed heuristics HEUR1, HEUR2, HEUR3, HEUR4, HEUR5 with TRtot =200, TRnimp =80, TMIP=10 sec, Pm = 0.5 and maximal solution time 600 sec. The obtained results show that for this dataset all heuristics found optimal solutions for at least 23 problem instances out of 25, i.e. for 92%, HEUR1 found 24 optimal solutions, i.e. 96%, and HEUR4 found all optimal solutions as CPLEX. The solution time for finding optimal solutions was relatively short for all 25 instances, the longest one is less than 1 minute. For 23 problem instances, the shortest time was attained by CPLEX, for 2 problem instances by HEUR3 and for one problem instances by HEUR5. These results show that the heuristic methods developed are capable of finding optimal solutions for small size problem instances. Further, we evaluate their performances on the problem instances of larger size randomly generated on the basis of the data provided by our industrial partner. Random generation is used in order to increase the number of different problem instances tested.

Results of tests for randomly generated problems of larger size

The second dataset includes 1200 randomly generated problems: 100 for each combination of optimization problem (A1, A2 and A3) and number of types of parts (4, 6, 8 and 10) to be machined. Their characteristics are presented in Table 4. Constraints were generated using the data provided by our industrial partner. Table 5 reports the obtained results for all heuristic methods, HEUR1-HEUR5, CPLEX and CPLEXO (with a heuristic solution used for a warm start). The methods are compared on the following indicators: NS is the number of found feasible solutions, OS is the number of found optimal solutions, BS is the number of found best known solutions (for instances which have not been solved to optimality), AD and MD are respectively average and maximal deviations (in percentages) of the found values of the objective function from the best known values, calculated on the basis of 100 instances for each optimization problem and each number of parts to be machined. The deviations were not calculated if the number of feasible solutions found by a method was less than 100. Table 6 presents a synthesis of the parameters for all problems and all instances. For the tested instances, neither CPLEX nor any heuristic could find all optimal solutions for all problems within limited computational time.

Overall, HEUR5 found more optimal solutions and more solutions with the best value of the objective function than other heuristic methods. Moreover, HEUR5 found more solutions with the best value than CPLEX within limited computational time. In general, CPLEX and CPLEXO provided similar results. However, CPLEXO found more feasible, optimal and best value solutions than CPLEX in total. CPLEXO also provided more optimal and best value solutions than any heuristic method, but has larger average deviation and maximal deviation than heuristic HEUR5 which provided the best results for these indicators. From this analysis, we can conclude that the best methods are HEUR5 and CPLEXO. In Figures 2 to 4, we present comparative diagrams of HEUR5 and CPLEXO for all problems and all numbers of types of parts to be machined to compare their performances in more detail. It can be observed that for A1, both methods found a comparable number of best solutions and the average deviation is very close for 4 and 6 types of parts, and then becomes much smaller for 8 and 10 types of parts for HEUR5. For A2, CPLEXO found more best solutions for 6 and 8 parts and fewer for 10 types of parts. The average deviation grows with the number of parts for CPLEXO and remains smaller for HEUR5. The same can be observed for A3.

In conclusion, we can note that HEUR5 could be the best option for large-scale problems, while for medium-scale problems CPLEXO remains the most efficient method. The analysis of the realized tests shows the applicability of the developed approaches in reallife conditions, where the design of complex reconfigurable rotary transfer machines for production of multiple parts can be substantially improved both in terms of the time reduction and the increase in quality of the obtained solutions. The lines appropriately designed can be used longer and for a more diversified portfolio of machine parts, both advantages improve the competitivity of the machining systems in the context of volatile markets and frequent customer demand changes.

Conclusion

The preliminary design of reconfigurable machines with rotary tables was considered. The machines are designed for batch production of different products. The mixed integer programming approach is used for modelling the selection of the best configuration and equipment. The selection deals with the simultaneous choice of: process plan; machining modules; machining conditions (cutting speed); line balancing and task assignment. Three classes of problems were considered, each of them for a specific mode of processing different parts. Thus, the models and algorithms proposed in this paper cover a very large range of industrial settings. The obtained problems are very hard optimization problems with mixed decision variables (discrete and continue). For these new decision problems, mixed integer programming models (MIP) were proposed. The solver Cplex was used to solve real case studies of small size.

For larger instances, five powerful new heuristics were developed and tested, first, on 25 real life problems of relatively moderate sizes, then on a set of randomly generated instances of larger sizes. The results of tests demonstrated the excellent performances of the suggested algorithms. The MIP model allows small instances to be solved quickly. The heuristics obtain similar results in term of quality, but with larger computational time. For large instances, the heuristics outperform the MIP model with Cplex, and give excellent results after a relatively short calculation time.

We can conclude that the proposed models and algorithms can be used for solving very complex real-life problems, exactly or approximately, depending on the size of problem. They can be useful in real life combinatorial design of reconfigurable machines with rotary tables at the preliminary design stage and help to reduce the computational burden for decision makers.

This study demonstrates how mathematical modelling of complex machining systems can be combined with advanced decision aid techniques to obtain sophisticated, but very efficient models and algorithms applicable in real life manufacturing practice for a range of systems. These are able to improve the decision process drastically as well as the quality of the obtained solutions. Despite the complexity of real-life systems, a large number of specific constraints and dependences among parameters, the article shows that it is possible to develop explainable and powerful mathematical models and algorithms for decision aid, which provides users with efficient solutions in a reasonable time. This is a pioneer work for the type of reconfigurable machines considered with rotary tables, turrets and machining modules for different modes of processing multiple parts. A profound knowledge of the process and production practices, combined with advanced mathematical tools for decision aid and optimization in the domain of line balancing and scheduling were necessary to reach the objectives of this research.

The article demonstrates that due to the use of sophisticated modelling techniques and operations research approaches, the efficient solution of real-life complex manufacturing problems becomes possible, in spite of their combinatorial complexity and a large number of constraints and interdependent parameters.

Our conclusions are of course limited to the considered classes of reconfigurable rotary machines and for considered modes of batch processing. Further studies will be focused on the introduction of new production modes and technologies as well as on extension of the proposed algorithms for other classes of machining systems. Another possible direction for further research is the introduction of sustainability criteria in the decision process and models (energy consumption, gas emissions, etc.)

Figure 1 .

 1 Figure 1. A rotary transfer machine with turrets.

  s=1, 2, …, nd, be a subset of operations for machining of elements of the s-th side of part d. Part d can be located at the zero position in different orientations H(d) that can be

  All operations pN are characterized by the following parameters: -the length (p) of the working stroke for operation pN, i.e. the distance run by the tool in order to complete operation p; -range [γ1(p), γ2(p)] of feasible values of feed rate which characterizes the machining speed; -set H(p) of feasible orientations of the part (indexes r{1, 2, …, rd} of rows of matrix H(d)) for execution of operation p d s N by spindle head or turret of type j (vertical, if hrs(d)=1, and horizontal, if hrs(d)=2).

  where the edges (p,p)E SS ((p,p)E SM , (p,p)E ST , (p,p)E SP ) if, and only if, operations p and p must be executed by the same spindle, at the same machining module (turret, position). At the same time, the possibility of performing operations from N at the same working position, by the same turret or by the same spindle head, is also defined by a number of constructional and technological constraints, for instance, mutual influence of combining operations, possibility of tool location in spindle head, turret, etc. These exclusion constraints can also be defined by undirected graphs G DM =(N,E DM ), G DT =(N,E DT ), and G DP =(N,E DP ) where the edge (p,p)E DM ((p,p)E DT ), (p,p)E DP )) if, and only if, operations p and p cannot be executed by the same machining module (turret, position). Infeasible combinations of part orientations are given by a set E DH . Each e = {(d1,r1),(d2,r2),…,(dk,rk)} represents a collection of pairs (part number d and row number of H(d)). These simultaneously prohibit orientation r1 for part d1, orientation r2 for part d2, and orientation rk for part dk. The set E DH includes {(r,d), (r,d)}, if there exist p

Let

  P=<P1,...,Pk,...,Pm> be a design decision with Pk=( The execution time t b (Pdkjl) of operations from Ndkjl with the feed per minute Гdkjl[max{γ1(p)|pNdkjl}, min{γ2(p)|pNdkjl}] is equal to (1):

  vertical spindle head is installed, then the second summand is equal to C4 and the third summand is equal to 0is no vertical machining, then the second, the third and the fourth summands are equal to 0, since N1=, Nk12=, k=1,…,m, studied problem is to determine: a) the number of positions m; b) orientations of parts H(d); c) the number bkj of machining modules of type j (j=1 for vertical and j=2 for horizontal) installed at the k-th position, k=1,...,m; d) subsets Ndkjl of operations from N d assigned to the l-th machining module of type j at the k-th position d=1,2,…

hF

  auxiliary variable, equal to 1 if elements of the d-th part are machined with the r-th orientation; d kjl F an auxiliary variable for determining the time of execution of operations from N d in the l-th machining module of spindle head or turret type j at the k-th position; d k an auxiliary variable for determining the time of execution of operations from N d at the k-th position; Different auxiliary variables are introduced for A1, A2 and A3. For problem A1: Fd is an auxiliary variable which determines the time of execution of all the operations from N d ; d k T is an auxiliary variable which is equal to d F if the k-th position exists and 0 otherwise.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of HEUR5 and CPLEXO for A1

  be connectivity components of G SSM including isolated vertices. Only one vertex Xpq decision variable, equal to 1 if the operation p from N is assigned to the block q=2(k-1)b0+(j-1)b0+l, i.e. l-th machining module of spindle head or turret of type j at the k-th position; auxiliary variable, equal to 1 if at least one task for machining elements of the d-th part is executed in the l-th machining module of spindle head or turret of type j at the k-th position;

	(operation or task) i is chosen from each SSM i N	, let (p)=i for all p SSM i N	and included
	into N. Let us introduce the following notation:
	ds kj Y auxiliary variable, equal to 1 if at least one operation from d s N is assigned to spindle
	head or turret of type j at the k-th position;
	d
	kjl

Y Ykjl auxiliary variable, equal to 1 if the l-th machining module of spindle head or turret of type j is installed at the k-th position;

Y1min auxiliary variable, equal to k if k is the minimal position covered by vertical spindle head or turret; Y1max auxiliary variable, equal to k if k is the maximal position covered by vertical spindle head or turret; Y1 auxiliary variable which is equal to 1 if the vertical spindle head or turret is installed;

Zk auxiliary variable, equal to 1 if at least one operation is assigned to the k-th position;

Table 2

 2 The best value of the objective function

	Examples HEUR1 HEUR2 HEUR3 HEUR4 HEUR5 CPLEX
	1	429.1	299.1	43.6	611.8	39.4	4.4
	2	25.3	59.6	6.6	158.5	7.8	0.3
	3	128.7	203.5	43.9	364.3	45.5	1.4
	4	128.8	116	33.5	263.1	33.2	1.1
	5	139.1	148.9	30.8	337	30.2	1
	6	234	397.3	22.9	600.7	26	2.5
	7	113.3	206	9	357	14.4	1.6
	8	305.2	91.6	69.8	161.3	13	266.1
	9	36.5	64.2	17.1	130.6	16	0.8

Table 3 Solution time (sec)

Table 4

 4

	Parameters of problems

Table 5 .

 5 Results for randomly generated problems

Table 6 .

 6 Synthesis of the results for randomly generated problems

	Method	NS	OS	BS	AD	MD
	HEUR1	1200	573	589	3.63	62.96
	HEUR2	1200	746	780	1.6	29.35
	HEUR3	1200	761	827	1.3	27.78
	HEUR4	1200	695	712	2.3	21.21
	HEUR5	1200	908	1065	0.37	17.65
	CPLEX	1100	928	1034	-	27.52
	CPLEXO	1200	1023	1068	0.92

Table 7

 7 presents a synthesis of the parameters for all instances per type of problem. We can see that for A1, the conclusions are the same as previously and for A2 and A3 they are different for HEUR5. For these problems, HEUR5 found more best value solutions than CPLEXO and still had smaller average and maximal deviations.

	A1	A2	A3
	METH
	NS OS BS AD MD NS OS BS AD	MD NS OS BS AD	MD
	HEUR1 400 170 173 6.2 62.96 400 219 224 2.1 16.52 400 184 192 2.5	16.09
	HEUR2 400 247 255 2.2 29.35 400 255 265 1.3 15.6 400 244 260 1.31 16.09
	HEUR3 400 250 262 1.9 27.78 400 260 288 0.95 12.39 400 251 277 1.02 16.09
	HEUR4 400 230 236 2.7 21.21 400 243 248 1.9 16.52 400 222 228 2.3	15.93
	HEUR5 400 322 361 0.5 17.65 400 295 354 0.3 8.2 400 291 350 0.36 16.09
	CPLEX 379 355 368 -	5.43 365 328 337 -	27.52 356 245 329 -	0
	CPLEXO 400 355 368 0.8 44.33 400 335 352 1.025 23.89 400 333 348 0.9425 15.6

Table 7 .

 7 Synthesis of the results per type of optimization problem

				4						6
		NS	OS	BS	AD	MD	NS	OS	BS	AD	MD
	HEUR1	300	236 236 1.08	12.5	300 157 161 5.33	62.96
	HEUR2	300	292 292 0.07	5.08	300 222 228 0.83	17.78
	HEUR3	300	293 293 0.12	10.61 300 200 209 1.36	27.78
	HEUR4	300	294 294 0.07	6.85	300 206 209 1.45	13.93
	HEUR5	300	295 295 0.05	5.08	300 257 275 0.25	10.61
	CPLEX	300	300 300 0	0	299 281 283 -	16.81
	CPLEXO	300	300 300 0	0	300 281 284 0.39	13.89
				8						10
		NS	OS BS AD	MD	NS OS BS AD	MD
	HEUR1	300	107 112 3.86	27.94 300 73 80 4.24	21.21
	HEUR2	300	132 139 2.77	29.35 300 100 121 2.73	28.79
	HEUR3	300	158 178 1.66	17.65 300 110 147 2.06	21.21
	HEUR4	300	125 131 3.35	17.78 300 70 78 4.28	21.21
	HEUR5	300	210 258 0.59	17.65 300 146 237 0.6	13.64
	CPLEX	276	192 242 -	27.52 225 155 209 -	0
	CPLEXO	300	242 261 1.4	44.33 300 200 223 1.89	31.65

Table 8 .

 8 Synthesis of the results per number of parts to be machined

Table 8

 8 presents a synthesis of the parameters for all instances according to the number of types of parts to be machined. We can see that for 4, 6, and 8 types the conclusions are the same as for Table6, but for 10 types of parts, HEUR5 found more best value solutions than CPLEXO and still with smaller average and maximal deviations.

Acknowledgements

This research was funded by the region Pays de la Loire, France (Project RMS) and by the National Agency for Scientific Research, France (project ReconfiDurable, ANR-21-CE10-0019-06).

Results of tests for the small industrial case studies