Mohamed N Boushaki

Chao Liu
email: liu@lirmm.fr

Benoît Herman
email: benoit.herman@uclouvain.be

Vincent Trevillot
email: v-trevillot@chu-montpellier.fr@mo-

Mohamed Akkari
email: hamed.akkari.orl@gmail.com

Philippe Poignet
email: poignet@lirmm.fr

Optimization of Concentric-Tube Robot Design for Deep Anterior Brain Tumor Surgery

Most of existing works on the tubes design optimization of concentric-tube robot (CTR) do not include the elastic stability in the optimization criteria. The only work which formulates the elastic stability in the objective function is based on scalarization method which is used in existing multi-objective design optimization. The objective function is formed by a set of weighted objective functions. The selection of the weights is crucial as the optimization results are greatly affected by them and could be misleading if these weights are improperly chosen. As an alternative optimization technique, we use Pareto grid-searching method to avoid this problem and allow a straightforward interpretation of the results following the selection criteria for the parameters to be optimized. This paper shows a three-tube CTR design based on Pareto gridsearching method in order to optimize the reachability and elastic stability of the CTR within a specific curvature range dedicated to the deep anterior brain tumor removal surgery.

I. INTRODUCTION

Typical surgical robots for Minimally Invasive Surgery (MIS) use rigid and straight instruments with an articulated tip. Small sized flexible robots are more suitable for MIS as they can reach surgical sites of high access difficulty [START_REF] Burgner-Kahrs | Continuum Robots for Medical Applications: A Survey[END_REF].

Concentric Tube Robots (CTR) belong to this category of flexible robots with the advantages of good steerability, controllability of the shape and small size. During last decade, CTR has drawn considerable attention [START_REF] Sears | A Steerable Needle Technology Using Curved Concentric Tubes[END_REF] [START_REF] Webster | Toward active cannulas: Miniature snake-like surgical robots[END_REF]. CTR is a kind of flexible robot whose concept is to insert concentrically pre-shaped flexible Nitinol tubes such that each of them is endowed with an axial rotation and linear translation. Previous works show the potential of using this technology in many medical applications, such as neurosurgery [START_REF] Anor | Algorithms for Design of Continuum Robots Using the Concentric Tubes Approach: A Neurosurgical Example[END_REF], intracardiac surgery [START_REF] Bedell | Dupont Design Optimization of Concentric Tube Robots Based on Task and Anatomical Constraints[END_REF], endonasal skull base surgery [START_REF] Burgner | A Bimanual Teleoperated System for Endonasal Skull Base Surgery[END_REF] and lung biopsy [START_REF] Torres | Task-oriented Design of Concentric Tube Robots using Mechanics-based Models[END_REF].

It is obvious that the design of the CTR tubes differs from one medical application to another. The tubes dimensions depend mainly on the specific task and the involved anatomical structure. To meet various performances indices, the robot tubes design should be optimized accordingly. The evolution of CTR tube design optimization can be traced from three aspects. Firstly, the design optimization follows the advances in kinematics modeling. In the early works on tubes design optimization, only basic kinematic model was used [START_REF] Bedell | Dupont Design Optimization of Concentric Tube Robots Based on Task and Anatomical Constraints[END_REF] [START_REF] Burgner | A Bimanual Teleoperated System for Endonasal Skull Base Surgery[END_REF]. Later on, more accurate kinematics model has been involved [START_REF] Burgner | On the Computational Design of Concentric Tube Robots: Incorporating Volume-Based Objectives[END_REF]. With more accurate kinematics, better robot design optimization could be obtained since the kinematic model is involved in the optimization calculation. Secondly, the design optimization evolution could also be analyzed from the view of clinical applications tackled seeing that many specific design optimization algorithms have been proposed for different clinical procedures [START_REF] Anor | Algorithms for Design of Continuum Robots Using the Concentric Tubes Approach: A Neurosurgical Example[END_REF] [START_REF] Bedell | Dupont Design Optimization of Concentric Tube Robots Based on Task and Anatomical Constraints[END_REF] [6] [START_REF] Torres | Task-oriented Design of Concentric Tube Robots using Mechanics-based Models[END_REF]. And another axis of this evolution is the transition from considering a discrete set of point in the objective function to a volumebased objective functions [START_REF] Burgner | On the Computational Design of Concentric Tube Robots: Incorporating Volume-Based Objectives[END_REF].

Elastic instability is a big challenge for CTR, which may occur when the difference of tube angles in the base of the robot is at the neighborhood of π. The rotation angle difference of the tubes at the distal part has no unique solution in the case of elastic instability, which generates a jump at this distal proximity (snapping) [START_REF] Dupont | Design and Control of Concentric-Tube Robots[END_REF]. The elastic stability decreases once the tubes curvatures and/or lengths increase, and when it drops below the stability limit the snapping phenomenon appears. Elastic instability could lead to serious safety problem during the robotized intervention if it has not been considered in the robot design or in the path planning.

The problem of elastic instability was not taken into account in CTR tube design optimization in literature until the work of [START_REF] Ha | Achieving Elastic Stability of Concentric Tube Robots Through Optimization of Tube Precurvature[END_REF], but it is limited to a particular case of two planar tubes with variable pre-curvature of the tubes. Bergeles et al. in [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF] have shown a way to analyze the stability when the number of tubes exceed two, but it is not applicable for arbitrary number of tubes. In [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF], the formulation of the optimization problem is based on scalarization method, and the objective function is a concatenation of weighted subobjective functions which makes the results highly dependent on the a-priori weights. The selection of the weights is a nottrivial and crucial task.

The targeted surgical application of our work presented in this paper is deep frontal lobe brain tumors where elastic instability is of high risk to the operation safety. Therefore, in this work the elastic stability is explicitly considered as an optimization criterion. Instead of scalarization method, a Pareto grid-searching method has been used which incorporates all performance criteria within the optimization process and address them simultaneously to find a set of nondominated designs in the objective space [START_REF] Unal | Multi-criteria Design Optimization of Parallel Robots[END_REF] and therefore avoid the weight choice problem. The objective is to find the tube curvatures which suit best the specifications of the addressed surgical application. Other parameters such as robot structure, lengths of the tubes, and variation intervals of the curvatures to be optimized are pre-selected based on CT scan images in order to respect the anatomical constraints during the insertion and to be as close as possible to the optimal parameters with the help of surgeons.

Another contribution of this work is to use the forward kinematic model in the objective function instead of the inverse kinematics which has been used in [START_REF] Bedell | Dupont Design Optimization of Concentric Tube Robots Based on Task and Anatomical Constraints[END_REF] [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF] to avoid the anatomical obstacles. The inverse kinematic model is complex and time consuming to calculate, and it requires a numerical optimization to be solved. The forward kinematic model used in this work is the torsionally compliant model developed in [START_REF] Dupont | Design and Control of Concentric-Tube Robots[END_REF]. Two mappings are required to obtain the pose G of the robot tip. The first one is to calculate the curvatures and the rotation angles of the tubes along the whole robot shaft in terms of the joint variables (base rotations and translations). This mapping is modeled by a function f 1 which is a nonlinear differential equation system with problem of boundary conditions. The second mapping f 2 is to integrate geometrically the obtained curvature u(s) in order to define the robot shaft including the tip pose. Based on the optimal parameters, the tubes are to be fabricated to integrate to our CTR system.

II. MULTICRITERIA OPTIMIZATION USING PARETO

METHODS

The mathematical formulation of a multi-objective optimization problem can be described as follows:

m x in{ϕ 1 (x), ϕ 2 (x), ....ϕ n (x)} subject to x ∈ (1)
where n > 1, x is the decision variables vector, ϕ i are the objective functions constituting the objective space, and is the set of constraints [START_REF] Caramia | Multi-objective Management in Freight Logistics[END_REF]. A decision vector x * ∈ is the Pareto optimal vector for a multi-objective problem if all other vectors x ∈ cause higher values for at least one objective function ϕ i or have the same values for all ϕ i . The collection of all the efficient solutions is called Pareto curve. The trade-off between objective functions is concluded from this Pareto curve.

A. Scalarization technique

Among all methods to solve Pareto optimal solutions for multi-criteria optimization, the scalarization method is the most commonly used one. In scalarisation method, the objective functions to be minimized are transformed to a single-objective scalar function by affecting weighting and summing them [START_REF] Kasprzak | An approach to facilitate decision tradeoffs in Pareto solution sets[END_REF]. The new objective function is defined as:

m x in n i=1 ω i • ϕ i (x), n i=1 ω i = 1, (ω i > 0) subject to x ∈ (2)
The minimum of this new optimization objective belongs to the Pareto curve. This method is easy to understand and implement. However, as the weights are pre-selected, they may not correspond to the relative importance of the objective functions. It is also difficult to change the weights or to find an heuristic iterative variation low making them change in order to reach a part of the Pareto curve. Moreover, selecting different weights vectors to perform several optimizations (in order to reach more parts of the Pareto curve) causes a large computational burden. Two other limitations of this scalarization technique are: a uniform spread of the weights does not produce a uniform propagation of the Pareto curve points i.e. some points of the former have not been considered; non-convex Pareto set parts are not reachable with the minimization of convex combinations of the objective functions.

B. Grid-searching technique

This technique allows to avoid the drawbacks of scalarization previously introduced due to the preselection of the weights. Three steps are to be followed in order to define the best trade-off between objective functions: • First, the Pareto set is generated by taking a sampling grid of the decision space x and then to calculate their corresponding objective functions. The obtained points are plotted in the objective space as illustrated in Fig. 1, where the Pareto curve can be approximated by a polynomial.

• The next issue after tracing the Pareto set is to find the suitable solution on it. Many approaches have been proposed in literature to solve this problem. The approach used in this work is the No-preference approach assuming "equal importance" of each criterion where the mostly used way to find the best solution is to calculate the L 2 norm for all the points as follows:

D L2 (x) = n i=1 (ϕ i (x) -ϕ min i ) 2 , (3)
where ϕ min i is the minimum of ϕ i , and then to select the point corresponding to the minimum of this distance [START_REF] Kasprzak | An approach to facilitate decision tradeoffs in Pareto solution sets[END_REF]. In case of maximization of the objective functions, the maximum of the L 2 norm is used. Note that if the equal importance of the objective functions is not considered, then one of the plot axes should be scaled before calculating the L 2 norm, and this scale corresponds to the weights.

• The last step is to map from the Pareto set back to the decision space. There are many methods in literature which ensure this mapping such as intersection of constant objective contours, etc. However, in this work we select the decision variables that correspond to the optimal objective solution obtained in the previous step. In this work, the specific surgical application is considered; removal of deep frontal lobe brain tumors cf. Fig. 2.a. It is noted that in year 2012, 70 % of the over 250 000 patients with brain tumor died in the world [10]. Usually for anterior skull base surgery, the procedure is transcranial, the access is created by severing a portion of the cranium. The surgeon needs to gently lift up the brain to access to the deep tumor located under the brain in order to accomplish the resection procedure cf. Fig. 2.b. Following this procedure, the brain will swell up after the operation, which creates a brain pressure and causes a large trauma for the patient. To alleviate this problem, an endonasal approach is proposed by using CTR which is flexible and adaptable with the anatomical complexity. Note that, similar as the case in [START_REF] Burgner | A Bimanual Teleoperated System for Endonasal Skull Base Surgery[END_REF] which addresses pituitary gland tumor removal, this operation requires bone drilling in order to access to the surgical site and reach the tumor. But the insertion path in our case from the nostril to the anterior wall of the frontal sinus (location of the drilling point) is surrounded with structures of high sensitivity cf. at the median between the two nostrils as shown in Fig. 3. The Z axis is oriented from the nose to the head back, the Y axis is oriented from the nose to the head top. The insertion point of the robot is at P point in Fig. 3.a and Fig. 3.b, with P defined as the origin of the robot frame F r . The orientation of F r is decided by two rotations: φ which orients the robot from the nostril to the frontal sinus in the Y -Z plan of the frame F b and allows to avoid the nasal turbinates; and The ψ which orients the robot to the left or to the right of the cerebral falx depending on the tumor location. Based on CT scan images and with more clarification and orientation of the surgeon, the space of the passage used to introduce the CTR from the nostril to the frontal sinus is identified roughly as a straight tunnel with variable diameter (from ≈ 10mm for smaller ellipse diameter to ≈ 30mm for bigger ellipse diameter) and a length of ≈ 72mm for the patient sample studied. These allowed workspace borders are illustrated by six green ellipses in Fig. 3. These anatomical constraints (orbits, optic canal, nerves...) should be considered in an eligibility criterion. In order to reduce the friction during the robot insertion, a large passive tube is pre-inserted before introducing the robot. This facilitates the insertion of the robot in the first part of its path which is characterized by high friction according to the surgeons. Once the navigation path is traversed, the distal section of the robot accesses to the anterior skull base to manipulate at the tumor location.

A scenario of a sphere tumor defined by its center position coordinates c and its radius a = 10 mm is shown by purple color in Fig. 3. The frequent scenario in reality is that the tumor is located at the left of at the right of the frontal lobe.

V. CTR CONFIGURATION SELECTION

By considering the aforementioned anatomical description, a robot configuration is proposed in this section following the surgeons advices.

A. Kinematics description

To fulfill accessibility and dexterity requirement, the candidate robot to be optimized is a three-tube CTR. The outermost tube T 1 is straight and the other two tubes T 2 and T 3 have two parts (straight and curved) as shown in Fig. 4. The tubes deployment (see Fig. 4) decides the number of tubes existing in one section, the number of sections, their order, lengths and curvatures (constant or variable curvature). In terms of reachability, the proposed configuration allows to benefit from both advantages of constant curvature and variable curvature of the robot distal section. The distal section can be composed by two pre-curved tubes which is a variable curvature section (S B4 ) or by only one pre-curved tube which means a constant curvature section (S C5 ) as seen in Fig. 4. The translation limits, and the curved and straight tube lengths define this transition constant/variable curvature of the distal section. The distal section manipulates the tool to accomplish the requested task (grasping, aspiration, electrocoagulation etc). For the proposed robot structure, the constant curvature section allows to reach the regions far

B. Actuation description

Concerning the actuation unit, it is composed of three blocks with one for each tube. Note that for the proposed actuation unit and for each tube, the rotation and translation are decoupled. Thanks to this decoupling, each rotation is controlled by one dedicated motor, which improves the control performance, facilitates the calibration procedure and reduces the modeling complexity. The tubes translations are produced by the blocks translations, which are ensured by the rotations of a nut-screw systems and guidance rails (see Fig. 5). Once a block translates, it pulls the blocks behind it. This mechanism allows to reduce the tubes lengths and hence the tubes torsions. The Block1 (front block) is dedicated to the outer-most tube T 1 , and the Block3 (arrear block) is dedicated to the inner-most one T 3 . The insertion of the robot from the nostril till the frontal sinus is ensured by the translation of the Block1 which pulls the two other blocks, producing a translation of the three tubes simultaneously. The distal proximities of the tubes are at the same limit. The shape of the robot shaft is held during this insertion. After this navigation insertion, the tube T 1 is held, then if only Block2 is translated, the distal section of the robot is a variable curvature section (see Fig. 4b), otherwise a part of the inner-most tube T 3 extends out of tube T 2 and forms a distal constant curvature section. Once the distal extremity of the robot accesses to the anterior cranial fossa through the posterior table of the frontal sinus, tube T 2 and/or T 3 are/is actuated in rotation and/or translation to accomplish the manipulation task.

VI. FORMULATION OF GRID-SEARCHING OPTIMIZATION ALGORITHM FOR CTR A. Decision space

For our specific surgical application, the manipulation part of the robot extends from the drilling point to the posterior part of the tumor. For a such distance, the lengths of the curved parts of the two inner tubes do not have significant effect on the optimization results. Moreover, reducing the number of decision parameters allows to decrease the computation time significantly. The lengths of tubes are then pre-selected suitable for the task respecting the workspace and anatomical constraints. The curvatures k 2 and k 3 of the tubes T 2 and T 3 respectively are considered as the decision variables. The localization of their ranges is firstly estimated from the geometry of the anatomical constraints, then obtained by executing the Algorithm.1 as introduce later.

B. Objective space

The objective space is defined by two criteria (objective functions):

• Reachability criterion: In order to evaluate the reachability of the surgical workspace (tumor), this workspace is sampled homogeneously to a set of points W (i). Once the euclidean distance between the workspace center and the tip position is less or equal than the work-space radius a, we check which point from W (i) is closer to this tip position, and we increase the score of reachability of this point W (i). For a given configuration (k 2 , k 3 ), a point from W (i) can be reached more than once. Consider m is the number of the point reached from W (i), and n is the number of the work-space samples points, the reachability is evaluated as:

C reach (k 2 , k 3 ) = m n × 100 (4)
• Stability criterion: For the elastic stability analysis, previous works rely on the evolution of the difference of distal tubes angles in terms of the base difference tubes angles [START_REF] Dupont | Design and Control of Concentric-Tube Robots[END_REF].

There is no systematic and general approach so far which allows to evaluate the elastic stability once the tube number exceed three. The case of three tubes has been treated by Bergeles et al in [START_REF] Bergeles | Concentric Tube Robot Design and Optimization Based on Task and Anatomical Constraints[END_REF].

In this work, a new way of stability evaluation is proposed. It relies on numerical evaluation of the tip trajectory smoothness. In case of elastic instability, jumps appear in the robot tip trajectory even if the tubes rotations variations are continuous and smooth. The elastic stability could be quantified by evaluating the euclidean distance χ 2 (1) + χ 2 (2) covered by the robot tip, since the robot shaft rotation is around the Z axis. In order to find out the jumps (case of instability) on the distal angle difference in ascending and descending, the bvpinit Matlab function is used to select the symmetric solution which shows a distal angle jump leading to a tip position jump. The new proposed approach is to quantify these jumps or discontinuities of the robot tip trajectory. Firstly, the numerical derivative of the euclidean distance in term of the joint variation (rotation sampling step) is calculated for each configuration (k 2 , k 3 ), then the standard deviations of these derivatives are calculated. The stability of a given configuration (k 2 , k 3 ) is inversely proportional to the derivative standard deviation.

C stab (k 2 , k 3 ) = normalization(1/std( ∆ χ 2 (1) + χ 2 (2) q )) (5)

C. Constraints

The first constraint is to consider only the decision variables satisfying (k 2 < k 3 ) during the grid sampling. This heuristic condition is recommended in [START_REF] Dupont | Design and Control of Concentric-Tube Robots[END_REF] because the proximal sections are constrained by the navigation path and then should be of small curvature. Moreover, the distal sections which are composed of inners tubes dedicated to manipulate the tip should be with high curvatures. The second constraint is defined by an eligibility condition under the notation C eli . The confliction of the robot shaft with the anatomical structure is checked using the ellipses borders. Once the Y coordinate of the robot χ(2) is equal to the Y coordinate as an ellipse, this shaft robot point is checked to see if it is inside the ellipse or not by using the two radii of the ellipse to calculate the distance from the borders. The robot shaft is considered satisfying the eligibility condition if and only if each point of the robot which has the same Y coordinate of an ellipse is inside its corresponding ellipse. This condition is first to be checked, i.e the configurations outside of the borders are rejected and will not be evaluated in the reachability and stability criteria.

D. Grid-searching algorithm execution

After the first execution of the Algorithm.1 with large sampling interval (15mm -1 ) which allows to locate the variation ranges of the decision variables (k 2 , k 3 ) ; where the optimal decision variables intervals corresponding to the highest reachability and stability may lie in. A second execution is conducted with smaller sampling interval (5mm -1 ) to ensure a finer mapping from the decision variables(k 2 , k 3 ) to the objective space (C stab , C reach ). For each configuration, (k 2 ,k 3 ) of the grid satisfying (k 2 < k 3 ), the eligibility condition C eli is checked. If it is respected, then two selection criteria are evaluated, i.e. the stability criterion C stab and the reachability criterion C reach as defined in last subsection. In the end, the criteria C stab and C reach are evaluated along all the configurations satisfying the eligibility condition. The next step is to calculate the L 2 norm of all obtained solutions in the objective space. As the objective of this optimization problem is to maximize both objective functions, the solution with highest L 2 norm is chosen and their corresponding (k 2 , k 3 ) are defined as the optimal decision variables.

Algorithm 1

Initialization:

Reachable(W (i)) = 0 with i = 0 : n for k2 = k2init : k 2f in , k3 = k3init : k 3f in do conf ig + + for d2 = 0 : d 2f in , d3 = 0 : d 3f in , θ2, θ3 = 0 : 2π do u(s) = f1(k2, k3, d2, d3, θ2, θ3) G(s) = f2(u(s)) if χ(2) = Y

VII. RESULTS AND DISCUSSIONS

The grid-searching optimization method as explained in last section was executed twice: the first time with larger parameters interval, (k 2 , k 3 ) and curvatures sampling resolution, then the second time to refine the optimization procedure with smaller parameters interval and curvature sampling. In order to accelerate the computation, the algorithm has been run on a server with the characteristics: Intel(R) Xeon(R), CPU E5-2695-v3 2.30GHz, (27 cores / 396 Go ram).

As the stability evaluation approach is new, the results obtained through the second execution are interpolated and shown in Fig. 6. It can be seen from this illustration that the stability criterion index increases along both X and Y axis, i.e. the curvatures (k 2 , k 3 ) are decreasing (curvature radii are increasing) which is consistent with the evolution stability in terms of curvatures presented in [START_REF] Dupont | Design and Control of Concentric-Tube Robots[END_REF].

Theoretically, the reachability improves by increasing the curvatures. However, in our study this rule can not be well satisfied as the anatomical obstacles are considered along the whole robot shaft, i.e. many points from the workspace are eliminated by the eligibility condition (the reason why the number of points in Fig. 7 is restricted (16 points) is because of the fact that non-suitable solutions have already been eliminated in the first execution of the algorithm).

From the mapping to the objective space shown in Fig. 7, the solution with highest In this work, a new optimization method for CTR design is presented using grid-searching Pareto which avoids the weight pre-selection problem of scalarization optimization methods. A new approach for elastic stability evaluation was proposed. Based on the specific task requirements and the proposed optimization method, optimal tube parameters can be decided following the proposed optimization approach.

The obtained optimized tubes parameters will be used to curve the tubes to be integrated to our CTR platform, and the addressed deep tumor removal surgical procedure will be tested through ex vivo lab test and in vivo cadaver test.

Fig. 1 :

 1 Fig. 1: Mapping from decision space to objective space.

Fig. 2 :

 2 Fig. 2: Transcranial intervention for frontal lobe tumor resection

  Fig. 3: Work-space characterization

TranslationFig. 4 :

 4 Fig. 4: Robot shaft sections

  Fig. 5: Robot

F

  b : base frame, its origin is O, Fr : robot frame, its origin is P φ: X-rotation angle of Rob, ψ: Y -rotation angle of Rob after its X-rotation c: tumor center coordinates f 1 : solves the nonlinear differential equation system of the forward kinematic model f 2 : numerical integration of the curvatures u(s) to define the robot shape G(s): homogeneous matrix of the pose of a point s from the robot shaft k 2 ,k 3 : curvatures of tubes T 2 and T 3 expressed in mm -1 ∆q: sampling step of the rotation joints, ∆s: sampling length χ: cartesian tip position in F b a: radius of the sphere tumor = 10 mm l s1 : length of the straight part of T 1 = 70 mm, l s2 : length of the straight part of T 2 = 30 mm l c2 : length of the straight part of T 3 = 31 mm, l s1 : length of the curved part of T 2 = 41 mm l s2 : length of the curved part of T 3 = 43 mm ∆k: curvature sampling step in the refining = 10mm -1 k 2 = 55 : ∆k : 70mm -1 ,k 3 = 35 : ∆k : 50mm -1

  L 2 norm with respect to the P min of coordinates (min(C reach ), min(C stab )) = (0, 46, 0, 48) is selected (D L2 = 0.5631). The point P opt surrounded by a black circle is the one with highest L 2 norm, and all the other points corresponding to solutions with different decision variables (k 2 , k 3 ) are included in the big blue circle with a radius equal to the D L2 . The coordinates of this solution point corresponding to the reachability and stability respectively are (0.72, 0.99). The decision variables corresponded to this solution point P opt are k 2 = 1 70 mm -1 , k 3 = 1 50 mm -1 .

TABLE I :

 I Notations and Dimensions

  Calculate C stab (conf ig) for all configurations Calculate C reach (conf ig) for all configurations

	Reachable(conf ig, ind) + +
	save(χ) to be used in C stab evaluation
	save(Reachable(conf ig)) to be used in C reach
	evaluation
	end if
	end if
	end if
	end for
	end for

borders then if the shaft is inside the ellipses borders then Test Reachability: if χ -c ≤ a then ind = index(min(χ -W (i))

This work was partially supported by the French Investissements d'Avenir program (Labex CAMI) under reference ANR-11-LABX-0004.


Martin Delacourt
email: martin.delacourt@univ-orleans.fr

Nicolas Ollinger
email: nicolas.ollinger@univ-orleans.fr

Permutive one-way cellular automata and the niteness problem for automaton groups

Keywords: reset Mealy automata, one-sided cellular automata, permutive cellular automata, periodicity problem, reversible computation

The decidability of the niteness problem for automaton groups is a well-studied open question on Mealy automata. We connect this question of algebraic nature to the periodicity problem of one-way cellular automata, a dynamical question known to be undecidable in the general case. We provide a rst undecidability result on the dynamics of one-way permutive cellular automata, arguing in favor of the undecidability of the niteness problem for reset Mealy automata.

Introduction

Finite-state automata provide a convenient nite description for dierent kinds of behavior generated by their computations. As such, Mealy automata [START_REF] Glushkov | The abstract theory of automata[END_REF] provide a nite description for the family of automaton (semi)groups that has proven its usefulness to generate interesting counter examples in the eld of group theory [START_REF] Bartholdi | Groups dened by automata[END_REF]. Several decision problems inspired by algebraic questions have been studied on automaton groups: the word problem is decidable whereas the conjugacy problem is undecidable [START_REF] Uni¢ | The conjugacy problem in automaton groups is not solvable[END_REF]. The general case of the niteness problem remains open [START_REF] Akhavi | On the niteness problem for automaton (semi) groups[END_REF] although special cases have been solved: Gillibert [START_REF] Gillibert | The Finiteness Problem for Automaton Semigroups is Undecidable[END_REF] proved that the problem is undecidable for semigroups and Klimann [START_REF] Klimann | Automaton semigroups: The two-state case[END_REF] that it is decidable for reversible Mealy automata with two states. The status of the niteness problem remains open for the class of reset Mealy automata.

Cellular automata [START_REF] Kari | Theory of cellular automata: A survey[END_REF] provide a nite description for a family of discrete dynamical systems, the endomorphisms of the shift dynamical system [START_REF] Hedlund | Endomorphisms and Automorphisms of the Shift Dynamical Systems[END_REF]. Decision problems inspired by dynamical questions have been investigated on cellular automata since the work of Amoroso and Patt [START_REF] Amoroso | Decision procedures for surjectivity and injectivity of parallel maps for tessellation structures[END_REF]. The computation nature of cellular automata lead to sophisticated construction techniques to establish the undecidability of various decision problems like the nilpotency problem [START_REF] Kari | The nilpotency problem of one-dimensional cellular automata[END_REF] or more recently the periodicity problem [START_REF] Kari | Periodicity and immortality in reversible computing[END_REF]. The status of the periodicity problem remains open for one-way cellular automata.

Without much surprise, cellular automata can be a valuable tool to establish undecidability results on Mealy automata. Indeed, Gillibert's result is inspired by Kari's proof of the undecidability of the nilpotency problem.

In this paper, we study the computational power of reversible permutive oneway cellular automata [START_REF] Boyle | Expansive invertible onesided cellular automata[END_REF][START_REF] Dartnell | Combinatorial constructions associated to the dynamics of one-sided cellular automata[END_REF]. Our rst contribution is a precise formalization of the connection between both open problems: the niteness problem for reset Mealy automata is decidable if and only if the periodicity problem for oneway cellular automata is decidable. Our second contribution is a technique to embed computation inside reversible one-way cellular automata using permutive automata. The technique is applied to prove a rst undecidability result on these objects.

Denitions

For a detailed introduction on Mealy automata, the reader is referred to Bartholdi and Silva [START_REF] Bartholdi | Groups dened by automata[END_REF] and to Kari [START_REF] Kari | Theory of cellular automata: A survey[END_REF] for cellular automata.

Mealy automata

A Mealy automaton is a deterministic complete 1-to-1 transducer (A, Σ, δ, ρ), where A is a nite set of states, Σ a nite alphabet, δ = (δ i : A → A) i∈Σ is the set of transition functions and ρ = (ρ x : Σ → Σ) x∈A the set of production functions.

The transition

x i|ρx(i) ----→ δ i (x) is depicted by x δ i (x) i ρ x (i)
The production functions naturally extend to functions on the set of nite words :

ρ = (ρ x : Σ * → Σ * ) x∈A with ρ x (au) = ρ x (a)ρ δa(x) (u).
The semigroup generated by the automaton is the set of all compositions of the production functions H = ρ x : x ∈ A . An automaton semigroup is a semigroup generated by a Mealy automaton.

A Mealy automaton is invertible if ρ x is a permutation of Σ for every x ∈ A. Note that it implies that every ρ x is also a permutation on Σ k for every k ∈ N. The group generated by a invertible Mealy automaton is G = ρ x , ρ -1

x : x ∈ A . An automaton group is a group generated by a Mealy automaton. An invertible Mealy automaton generates a nite group if and only if it generates a nite semigroup [START_REF] Akhavi | On the niteness problem for automaton (semi) groups[END_REF].

Finiteness problem Given an invertible Mealy automaton, decide if the generated group is nite.

A Mealy automaton is reset if, for each transition x i|ρx(i)

----→ δ i (x), the output state δ i (x) depends only on the input letter i and not on the input state x, that is δ i (x) = f (i) for some letter-to-state map f : Σ → A. To simplify notations, such an automaton will be denoted as (A, Σ, f, ρ).

When studying the decidability of the niteness problem restricted to reset automata, one can focus on the case Σ = A and f = Id as stated below.

Lemma 1. The group generated by a reset Mealy automaton (A, Σ, f, ρ) is nite if and only if it is the case for the automaton (Σ, Σ, 1, ρ ) with ρ x (i) = ρ f (x) (i).

Proof. Let G be the group generated by (A, Σ, f, ρ) and H be the group generated by (Σ, Σ, 1, ρ ). As every generator ρ x of H is a generator ρ f (x) of G then H is a subgroup of G. A generator ρ y of G with y ∈ A \ f (Σ) is not a generator of H, however as y can only be an initial state of a transition, it only impacts the size of G by a factor n! where n is the size of Σ.

Cellular automata

A one-way cellular automaton (OCA) F is a triple (X, r, δ) where X is the nite set of states, r is the radius and δ : X r+1 → X is the local rule of the OCA. A conguration c ∈ X Z is a biinnite word on X. The global function F : X Z → X Z synchronously applies the local rule: F(c) i = δ(c i , . . . , c i+r ) for every c ∈ X Z and i ∈ Z. The spacetime diagram ∆ : Z × N → X generated by an initial conguration c is obtained by iterating the global function: ∆(k, n) = F n (c) k for every k ∈ Z and n ∈ N. Following Hedlund [START_REF] Hedlund | Endomorphisms and Automorphisms of the Shift Dynamical Systems[END_REF] characterization of cellular automata as endomorphisms of the shift, we assimilate an OCA, with minimal radius, and its global function.

Notice that OCA are the restriction of classical cellular automata (CA) where a cell only depends on other cells on the right side. The identity function Id, the left shift map σ l and the XOR rule are OCA, respectively encoded as (X, 0, 1), (X, 1, (x, y) → y) and ({0, 1}, 1, ⊕), whereas the right shift map σ r is not.

A state x ∈ X is quiescent if δ(x, . . . , x) = x. A conguration is nite if it contains the same quiescent state x everywhere but on nitely many positions.

An OCA is (left

) permutive if the map x → δ(x, x 1 , . . . , x r ) is a permutation of X for every (x 1 , . . . , x r ) ∈ X r . An OCA is periodic of period T > 0 if F T = Id.
Periodicity problem Given a cellular automaton, decide if it is periodic.

An OCA is reversible if its global function is bijective with an inverse that is also an OCA. A periodic OCA F of period T is reversible of inverse F T -1 . Following Hedlund [START_REF] Hedlund | Endomorphisms and Automorphisms of the Shift Dynamical Systems[END_REF], the inverse of a bijective OCA is always a CA, however usually not one-sided. The following lemmas assert that this technical issue disappears by considering only permutive OCA. Lemma 2. Every reversible OCA is permutive.

Proof. Let F be a reversible OCA. As both F and F -1 are OCA, F is bijective on

X N too. Let (x, x , x 1 , . . . , x r ) ∈ X r+2 and c = x 1 • • • x r x ω r . If x = x , as F(xc) = F(x c), we have δ(x, x 1 , . . . , x r ) = δ(x , x 1 , . . . , x r ).
Lemma 3. Every bijective permutive OCA is reversible.

Proof. Let F be a bijective permutive OCA. By permutivity

f : X × X N → X dened by f (y, c) = x where F(xc) = yF(c) is well dened and permutive in its rst argument. Let u, u ∈ X -N and v ∈ X N . Let w ∈ X -N and v ∈ X N be such that F -1 (u v ) = wv. Let w ∈ X -N be dened recursively by w i = f (u i , w i-1 • • • w 0 v). By construction F -1 (uv ) = w v. As F -1 (uv ) and F -1 (u v ) are equal on N for all u, u , v the CA F -1 is an OCA.
Note that the inverse of a bijective permutive OCA can have a larger radius. However, as bijectivity is decidable for cellular automata [START_REF] Amoroso | Decision procedures for surjectivity and injectivity of parallel maps for tessellation structures[END_REF] and as bijectivity is preserved by grouping cells [START_REF] Delorme | Bulking ii: Classications of cellular automata[END_REF], when studying the decidability of the periodicity problem restricted to permutive OCA, one can focus on the case of reversible permutive OCA with radius 1 and inverse radius 1 syntactically characterized by the following lemma (already proven in [START_REF] Boyle | Expansive invertible onesided cellular automata[END_REF]). Lemma 4. An OCA (X, 1, δ) is reversible with inverse radius 1 if and only if it is permutive and for all x, y, x , y ∈ X, if δ(x, y) = δ(x , y ) then π x = π x where π y maps x to δ(x, y).

From now on, we only consider these OCA.

3 Linking niteness and periodicity Reset Mealy automata with Σ = A and f = Id and permutive OCA of radius 1 are essentially deterministic complete letter-to-letter transducers of a same kind, as depicted on Figure 1. The following proposition formalizes this link. Proof. First note that the following equations hold for all u 0 , u 1 , . . . , u k ∈ Σ:

ρ u k (u k-1 , u k-2 , . . . , u 0 ) = ρ u k (u k-1 )ρ u k-1 (u k-2 ) . . . ρ u1 (u 0 ) δ(u 0 , u 1 , . . . , u k ) = ρ u1 (u 0 )ρ u2 (u 1 ) . . . ρ u k (u k-1 )
By extension, for all k > t > 0 and for all words u ∈ Σ t and v ∈ Σ k , the following equation holds:

ρ u (v) k = δ t (v k , v k-1 , . . . , v k-t ).
Suppose now the group generated by the reset Mealy automaton is nite. Let a ∈ Σ be any letter and let n be the order of ρ a , then ρ a n = (ρ a ) n = Id thus δ n = Id, the OCA is periodic.

Conversely, let n be the period of the OCA. By previous remarks, for all k > 0 and words

u, v ∈ Σ n , w ∈ Σ k , the image ρ u (vw) is v w for some v ∈ Σ n .
The set of ρ u generates a subgroup of permutations of Σ n , which is nite, when u takes all possible values in Σ n . The automaton group is nite.

Corollary 1. The niteness problem restricted to reset Mealy automata is decidable if and only if the periodicity problem restricted to OCA is decidable.

Notice that the situation described by this corollary is optimal: if the problem is decidable, Mealy automata is the right setting to prove this result and the decidability of the periodicity for OCA will be a consequence; if the problem is undecidable, cellular automata is the right setting to prove this result and the undecidability of the niteness problem will be a consequence. The remainder of this paper is dedicated to prove that computational phenomena do appear inside the dynamics of permutive OCA, advocating for the undecidability of the problem. Conjecture 1. The niteness problem is undecidable. Given a permutive OCA, we show how to build a reversible OCA that can simulate every spacetime diagram of the original. The idea is to slow down the computation by delaying each state using a xed number of distinct copies per state and perform a transition only after going through every copy. Adjacent columns of states are then desynchronized to obtain reversibility as a consequence of permutivity. Denition 1. Let F be an OCA (X, 1, δ) and let

1 ≤ k ≤ n -1. The (n, k)- embedding F of F is the OCA (X , 1, δ ) where X = 1≤i≤n {x (i) :
x ∈ X} and such that: α) , y (β) = δ(x, y) (1) if α = n and β = k x (1+(α mod n)) otherwise Fig. 2 illustrates the embedding.

∀x (α) , y (β) ∈ X δ x (
Lemma 5. The (n, k)-embedding of a permutive OCA F is reversible.

Proof. The local rule of the inverse OCA τ can be dened by τ (z (1) , y (k+1) ) = x (n) for all x, y ∈ X such that δ(x, y) = z; τ (x (i) , * ) = x ((i-1 mod n)+1) otherwise. The idea of the embedding is to desynchronize adjacent columns by shifting them vertically of some constant k between 1 and n -1. When two consecutive columns are not correctly arranged, they do not interact. Lemma 6. There exists an injective transformation of spacetime diagrams of F (in X Z×N ) into spacetime diagrams of F (in X Z×Z ): for every c ∈ X Z , there exists a unique conguration c ∈ X Z for F with

∀m ∈ Z, p ∈ N, ∀1 ≤ i ≤ n, F m(n-k)+pn+i-1 (c ) m = (F p (c) m ) (i) x y z -→ x (1)
x (2) x (3) z (1) z (2) z (3) y (1) y ( 2) Remark 1. If the OCA has a particular state 0 such that δ(x, 0) = x for every

y (3)
x ∈ X, we can keep a unique version of state 0 by identifying the 0 (i) as a unique state 0 where δ (0, y (k) ) = (δ(0, y)) (1) δ (0, * ) = 0 where * can be anything

δ (x (i) , 0) = x (1+(i mod n))
This lemma allows to transfer results from F to F , in the sequel we prove the following result for permutive OCA and obtain it for reversible ones:

Reachability problem Given a reversible OCA with a quiescent state 0 ∈ X and two states x, y ∈ X, decide if y appears in the spacetime diagram generated by the initial conguration ω 0.x0 ω . Theorem 1. The reachability problem is undecidable for reversible OCA.

Main construction

To prove the theorem, we need to embed some Turing complete computation into permutive OCA. This goal is achieved by simulating multi-head walking automata. This section describes the simulation of these automata by permutive OCA.

Multi-head walking automata

A multi-head walking automaton consists of a nite number of heads on the discrete line, each one of them provided with a state out of a nite set. At each step, they can only interact (read the state) with the heads that share the same cell, update their state and move. Initially, all the heads are in position 0.

Denition 2. A k-head walking automaton is dened by

(Σ, I, F, (f i , g i ) 1≤i≤k )
where Σ is a nite alphabet that does not contain ⊥, I ∈ Σ is the initial state, F ⊆ Σ is a set of nal states, ∀i, f i : (Σ ∪ {⊥}) k → Σ and ∀i, g i : (Σ ∪ {⊥}) k → {-1, 0, 1} are the update functions for the state and the position.

A conguration of this automaton is a k-tuple a = (a i , s i ) 1≤i≤k ∈ (Z × Σ) k and its image conguration is

(a i + g i (b i ), f i (b i )) 1≤i≤k where ∀1 ≤ i ≤ k, b i (j) = s j if a i = a j and ⊥ otherwise.
Starting from conguration (0, I) k , the automaton computes the successive images and stops only if the position of every head is 0 and their states are nal.

A multi-head walking automaton can mimic a counter with 2 heads. Turing completeness is achieved by simulating 2-counter Minsky machines.

Proposition 2. The halting problem of 4-head walking automata is undecidable.

Finite congurations

Every nite conguration can be written c = ω 0.u 0 0 ω for some u 0 ∈ X n , n ∈ N, hence its successive images can only be

F k (c) = ω 0v k .w k 0 ω for some v k ∈ X k , w k ∈ X n ,
that is: the nite non-quiescent word extends to the left, one cell at each step. We also have the following result.

Lemma 7 (Coven et al [START_REF] Coven | Prevalence of odometers in cellular automata[END_REF]). In the spacetime diagram associated to a nite conguration, every column is periodic.

P-signals

Most cellular automata constructions use signals as elementary geometrical building blocks. Unfortunately, it is not possible to embed a classical signal in a permutive OCA nor is it possible to directly simulate multi-head automata. We provide P-signals as a technique to replace signals of speed k by an (n, k)-embedding of the front line of the spacetime diagram of the XOR with only one non-zero cell. This will eectively allow to send a signal through space in a reversible permutive OCA.

For our construction, we x n = 7, hence we can use speeds from 1 to 6. Due to remark 1, the states of these signals will be denoted Z 8 = {0, 1, 2, 3, 4, 5, 6, 7}. Next lemma states that dierent speeds induce similar spacetime diagrams up to a vertical shift on each column. Lemma 8. Denote F k and F k the OCA corresponding to signals of speeds k and k and let c be the conguration ω 0.10 ω then we have

∀m ≤ 0, ∀p ≥ 0, F p k (c) m = F p-(k-k )m k (c) m .

2-recognizability

Our construction relies on the ability to identify some specic sets of spacetime positions. We achieve this goal by considering products of independent P-signals.

To avoid boring calculations, we rely on the following lemma to assert some regularity of P-signals and rely on the theory of p-recognizable sets of tuples of integers [START_REF] Bruyere | Logic and p-recognizable sets of integers[END_REF] to characterize these positions.

Lemma 9. The spacetime sequence ∆ : N 2 → Z n+1 of the (n, k)-embedding F of the XOR, where ∆(x, y) = F y ( ω 0.10 ω ) -x , is 2-recognizable.

Proof. The spacetime sequence ∆ of the (n, k)-embedding F is generated by the

2-substitution s : Z n+1 → Z 2
n+1 uniquely dened, as per Fig. ??, by

s(0) = 0 0 0 0 s    n . . . 1    = n • • • (k + 1) k • • • 1 n • • • (k + 1) • • • 1 0 • • • 0 n • • • • • • • • • 1 0 • • • 0
Indeed, the substitution rule is compatible with the local rule of the OCA.

Computation windows

Let F denote the OCA we are constructing. We use 4 P-signals of speeds 2, 3, 5 and 6 as foundations of F they allow us to build computation windows. Let c be the initial conguration which is null everywhere except for c(0) = (

We use Lemma 9 twice for each of the following lemmas, that is, we have a characterization of the set of positions where some specic state pairs appear. Both lemmas are illustrated in Fig. 3, 4 and 5.

Lemma 10.

∀m ≤ 0, p ≥ 0,

F p (c) m = (1, 0, _, _) F p (c) m+1 = (0, 3, _, _) ⇔ ∃h ∈ N * , m = -8 h p ≡ -3m -1[-14m]
Lemma 11.

∀m ≤ 0, p ≥ 0,

F p (c) m = (_, _, 7, 6) F p (c) m+1 = (_, _, 5, 0) ⇔ ∃h ∈ N * , m = -8 h p ≡ -12m -1[-14m]
Hence we add a fth layer using alphabet {0, 1} with the rule: δ 5 ((1, 0, _,_,x),(0,3,_,_,_)) = 1x δ 5 ((_, _,7,6,x),(_,_,5,0,_)) = 1x δ 5 ((_, _,_,_,x),(_,_,_,_,_)) = x otherwise. Actually, the same arguments work for columns -2•8 h , h ∈ N * and -4•8 h , h ∈ N * , hence we use them all. We therefore call computation windows the vertical segments in the spacetime diagram where the fth layer contains 1:

(m, p) in a computation window In light gray, speeds 2, 3, 5 and 6 Psignals together allow to draw the whole computation windows.

⇔ ∃h ∈ N * , m = -2 h , p ≡ p[14 • 2 h ] and 3 • 2 h ≤ p < 12 • 2 h .
Fig. 6. The heads are supported by speed 4 P-signals that pass through the computation windows. The windows are large enough so that this property remains even after they are shifted up or down. Fig. 7. When the head arrives in a computation window, it is here shifted upward, then upward again and downward in the third window. The state also changes then (red, later blue, orange and blue again). In light gray, another head whose position does not change.

Computing heads

Every head of the walking automaton is simulated on a new layer. It is composed of a support that is mainly a speed 4 P-signal, and an internal state that belongs to the set of states Σ of the walking automaton. We need 4 heads, hence there are 4 additional layers (6 to 9).

The idea is that the heads move like speed 4 P-signal carrying internal states until they reach a computation window (as illustrated in Fig. 6). Then, depending on the context, they update their internal state and eventually shift upward or downward. Each shift corresponds to a move of the head of the walking automaton. Hence its position is encoded by the global shift applied to it, call it the height of the head.

More formally, each layer representing a head uses the alphabet ([1..7] × Σ) ∪ (0, ⊥). Outside computation windows, the rst component follows the rule of speed 4 signals, and the second component is maintained if possible or otherwise taken from the right neighbor.

Suppose now that we compute the new state of cell m at time p with (m, p) inside a computation window. Denote by x, y and z the states of cells (m, p), (m + 1, p) and (m, p + 1) in the spacetime diagram. The following rules apply to layer 6 (the same applies for other layers):

if y 6 = (0, ⊥), then apply standard rules. if y 6 = (0, ⊥), look at every other layer where the support state is the same and apply f 1 and g 1 , the update functions of the walking automaton, with these states. If the result of g 1 is -1 (resp. 0, 1), apply the rule of a speed 3 (resp. 4, 5) signal to the rst component. If the support of z 6 is not 0, then the internal state is given by f 1 .

Simulation

Finally, the OCA we build has 9 layers : c(m) = (0, 0, 0, 0, 0, (0, ⊥), (0, ⊥), (0, ⊥), (0, ⊥)) elsewhere. First note that the height can vary of at most 1 each time the head crosses a column -2 h , h ∈ N, hence: Lemma 12. Given any simulating head, while its height is between -h and h in column m = -2 h , h ∈ N, every non (0, ⊥) value of the layer in this column is inside a computation window. This means that each head has to apply the update rule when arriving in column -2 h , h ∈ N. Now check that every time it does, one step of the computation of the walking automaton is simulated. The main point is to ensure that two heads of the walking automaton are on the same cell at step h if and only if the support of their simulating heads coincide on column 2 h -1. This is true since, for any head whose height is s, its support takes value 1 exactly in cells (2 h -1, p) with p ≡ 4 • (2 h -1) + s [7 • 2 h ]. This property is due to the XOR rule, with the (7, 4)-embedding. Fig. 7 illustrates this behaviour with two heads (one does not move to simplify the representation). This completes the simulation and the proof of theorem 1.

One step further

Orbit periodicity problem Given a reversible OCA with a quiescent state 0 ∈ X and a state x ∈ X, decide if the spacetime diagram generated by the initial conguration ω 0.x0 ω is periodic. Theorem 2. The orbit periodicity problem is undecidable for reversible OCA.

If the computation halts, the windows and the 4 P-signals that help determine the computation heads keep progressing eternally. To reach periodicity, it is necessary (and enough thanks to Lemma 7) to kill them all. It is possible to do so by giving killing orders to the heads. At the halting step h 0 , each head is given the responsibility to kill one of the P-signals while sacricing itself. They have to slow down or speed up to meet the corresponding P-signal. Again, we use Lemma 9 to prove that the meeting can happen on column 2 h0+1 with a local context that does not happen elsewhere.

Fig. 1 .

 1 Fig. 1. linking reset Mealy automata to permutive OCA

4

  Computing with permutive OCA As shown in Fig 1, the time goes up in every representation of this paper.

Fig. 2 .

 2 Fig. 2. The (3, 1)-embedding is given by the transformation of the local rule (to the left). The space-time diagram of the XOR OCA with a unitary conguration (at the center) is transformed into a space-time diagram of a reversible OCA (to the right).

Fig. 3 .

 3 Fig. 3. Speeds 2 and 3 P-signals determine the positions of the lower points of the computation windows as stated in Lemma 10.

Fig. 4 .

 4 Fig. 4. Speeds 5 and 6 P-signals determine the positions of the upper points of the computation windows as stated in Lemma 11.

Fig. 5 .

 5 Fig. 5.

  4 P-signals to determine the computation windows on the fth layer, and 4 to simulate the 4 heads of the Turing universal walking automaton. The initial conguration is the nite conguration c with c(0) = (1, 1, 1, 1, 0, (1, I), (1, I), (1, I), (1, I))