Andony Arrieula

Hubert Cochet

Pierre Jaïs

Michel Haïssaguerre

Nejib Zemzemi

Mark Potse

In-Silico Data Based Machine Learning Technique Predicts Premature Ventricular Contraction Origin Coordinates

Premature ventricular contraction (PVC) can induce ventricular tachycardia or ventricular fibrillation. Drugresistant PVCs can be cured by catheter ablation, but the accurate localization that this requires can be difficult and time-consuming. An accurate pre-procedural estimate of the origin could make the procedure more efficient.

We propose a machine-learning method for accurate pre-procedural origin estimation. It uses a database of paced 12-lead ECGs with known pacing locations and presents its results on an imaging-based model of the patient. The method was tested using 7 realistic heart-torso models with hundreds of PVCs everywhere in the ventricles.

We found that increasing the number of patients in the training database increased the accuracy of the predictions. The optimal number of pacing sites per patient in the training dataset was about 25, resulting in a prediction error around 15 mm.

We conclude that our method gives a good indication to clinicians to efficiently start a pace-mapping during a catheter ablation procedure. It can be complemented with an intra-procedural method that uses the patient's own paced beats to refine the prediction.

Introduction

Premature Ventricular Contraction (PVC) is a relatively common event where the heartbeat is initiated in the ventricles rather than by the sinoatrial node. Incidental PVCs are usually harmless, but multiple PVCs can induce a ventricular tachycardia or fibrillation.

To cure very frequent and drug-resistant PVCs, the site of origin can be ablated with an endocardial catheter. We previously proposed a patient-specific method to localize PVC origin based on iterative pace mapping, using only a 12-lead ECG and no cardiac imaging [START_REF] Arrieula | An improved iterative pace-mapping algorithm to detect the origin of premature ventricular contractions[END_REF]. However, this method can only work during the intervention, because it needs to be trained with stimulated beats. Cardiologists are often interested in a pre-interventional estimate, and often dispose of cardiac imaging performed before the procedure. To provide such an estimate we propose a machine learning method, based on cardiac imaging and a database of multiple patients.

Recent studies have used machine learning to classify PVC origin by regions, but some work with very large regions [START_REF] Nakamura | Prediction of premature ventricular complex origins using artificial intelligence-enabled algorithms[END_REF], are only applicable for endocardial origins [START_REF] He | Localization of origins of premature ventricular contraction in the whole ventricle based on machine learning and automatic beat recognition from 12-lead ECG[END_REF], or only for left ventricular origins [START_REF] Sapp | Real-time localization of ventricular tachycardia origin from the 12-Lead electrocardiogram[END_REF]. Our proposed method, in contrast, can predict origins in both ventricles and at any depth, and can also predict 3D coordinates of origin, using only the QRS complex of the 12-lead electrocardiogram (ECG).

The method works as follows: cardiac imaging data are used to create a 3D mesh of the ventricles of the patient. We register this mesh to a reference bi-ventricular mesh. tients, also registered to the reference mesh. Thus, using features of the 12-lead ECGs and the corresponding origin coordinates, we can train a model to predict the 3D coordinates of a PVC origin on the reference mesh. The last step is to translate the prediction to the patient's own mesh.

In this study, we are interested in determining the impact of the number of patients and pacing sites in the database on the results.

Methods

Database

To construct the database, we used detailed anatomical models of the heart and torso of 7 cardiac patients. These models were previously created from computed tomography data. Simulations were performed on 3-dimensional computational meshes of normal ventricles with a uniform edge length of 200 µm. The activation of the ventricles was simulated with a monodomain reaction-diffusion equation and the Ten Tusscher -Noble -Noble -Panfilov ionic model [START_REF] Ten Tusscher | A model for human ventricular tissue[END_REF]. The 12-lead ECG was computed using the leadfields method [START_REF] Potse | Scalable and accurate ECG simulation for reactiondiffusion models of the human heart[END_REF]. The simulations were performed with a recent version of the propag-5 software [START_REF] Potse | Scalable and accurate ECG simulation for reactiondiffusion models of the human heart[END_REF][START_REF] Krause | Hybrid parallelization of a large-scale heart model[END_REF].

For each model, a few hundreds of paced beats were simulated. Their origins were chosen randomly but well distributed to cover all ventricular regions. The details of the database are shown in Table 1.

Reference Mesh

For each patient, a 3D mesh of the ventricles was created. The mesh of patient 1 was chosen as the reference mesh. The other meshes and the pacing sites defined on them were scaled to the reference mesh using the iterative closest point method.

Machine Learning Algorithm

As an input of our machine learning algorithm, we used the features proposed by Giffard-Roisin et al. [START_REF] Giffard-Roisin | Noninvasive personalization of a cardiac electrophysiology model from body surface potential mapping[END_REF]. For each QRS complex, the following features were extracted: 1) QRS integral, 2) absolute potential of the global extremum, 3) timing of the global extremum, 4) sign of the potential of the global extremum, 5) number of zero crossing, 6) number of local extrema, 7) sign of the potential of the first extremum. Thus, for each paced beat we had 12 × 7 = 84 features.

We chose the Gradient Boosting method, using the implementation in the scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] package. We performed 500 boosting stages, with a maximum depth of 5 for the trees. At each step, we considered only the square root of the number of features when looking for the best split, which was performed only if at least 1 % of the samples would be in the new nodes. The learning rate to shrink the contribution of each tree was set to 0.05. A random subsample of 80 % of the training dataset was selected at each tree, so the method becomes Stochastic Gradient Boosting.

To assess the effect of the mesh variability on the database, we gradually increased the database by including a new patient geometry at each step. We then had 6 cases: In the first case we learned on simulations performed on a single geometry, in the second case we used two geometries, and so on. For each case, we evaluated prediction errors. For the i-th case, we had C i 6 possible combinations to construct the learning database, built using i geometries. For a given combination, we randomly chose 50 pacing sites per patient, so the size of the database is i × 50 simulations. We then predicted the stimulation sites for the remaining geometries. In order to avoid bias in the choice of the stimulation sites, this step was repeated 100 times and the error between the predicted sites and the real location of the stimulation sites was computed.

In order to reduce the number of stimulation sites used in the training phase, in the 6th case, we gradually decreased the number of selected pacing sites per patient from 50 to 4, and we computed the corresponding prediction error.

Results

Fig. 1 shows the effect of geometric variability in the database. The i-th bar group corresponds to the prediction error in the i-th case described above. For instance, the first bar group provides the mean and the standard deviation (20.7 ± 9.5) of the prediction error when learning on 50 simulations performed on a single geometry, while the last group provides the mean and the standard deviation (14.8 ± 7.1) of the prediction error when learning on 6 different geometries. The mean error of the prediction decreases as the number of patients in the database increases.

Fig. 2 shows the sensitivity of the method to different numbers of pacing sites per patient in the training database. The figure shows that the prediction decreases with increasing number of pacing sites per patient in the training database. However, the improvement in the prediction accuracy was negligeable when more than 30 pacing sites were used.

Discussion

In this study, we proposed a machine learning method to localize PVC origins, using the 12-lead ECG and cardiac imaging data, before the clinical intervention.

For the preparation of an ablation procedure it is important to know in advance in which ventricle the origin lies. With 15 mm accuracy, our method will in most cases be able to indicate the correct ventricle, and if the prediction is close to the interventricular septum the cardiologist will be warned that access to both ventricles may be needed during the procedure. Our results show that the accuracy of the model can be improved by adding more patients to the training database. It also shows that we don't need more than 30 pacing sites per patient in order to be efficient, which is clinically feasible. However, this method requires cardiac imaging, which is not always available for catheter ablation of PVCs. Since they are generated randomly, some pacing sites in the current database cannot be reached by a catheter. This suggests to construct a model that learns only on accessible pacing sites, which has not been done in the current work. The current accuracy of the model, around 15 mm of error, is not enough for direct ablation, but useful for pre-interventional planning. Moreover, the method can be complemented with an intra-procedural method that uses the patient's own paced beats to refine the prediction, as we proposed in our previous work [START_REF] Arrieula | An improved iterative pace-mapping algorithm to detect the origin of premature ventricular contractions[END_REF]. Also, this test was performed with only 6 patients in the training database, and our data suggest that additional patients can increase the accuracy of the predictions. This method can already be validated with cathlab data after procedure, by comparing our predictions and the ablation sites.

Figure 1 .

 1 Figure 1. Mean prediction error over the pacing sites of the considered patient as a function of the number of patients in the training dataset.

Figure 2 .

 2 Figure 2. Fitted curves showing the mean prediction error over the pacing sites of the considered patient as a function of the number of pacing per patient in the training dataset.

Table 1 .

 1 Number of paced beats by patient and tissue.

	In addition we use a database of 12-lead ECGs previously
	recorded after pacing at known locations in multiple pa-
	Patient	LV endo	LV epi	RV endo	RV epi	All sites
	#1 (M,44)	99	94	99	95	387
	#2 (F,31)	99	88	99	97	383
	#3 (F,20)	99	95	99	99	392
	#4 (M,59)	99	95	99	98	391
	#5 (F,15)	99	95	99	97	390
	#6 (M,35)	99	92	99	96	386
	#7 (M,44)	99	95	99	99	392
	All patients	693	654	693	681	2721

Acknowledgements

This work was supported by the Région Nouvelle-Aquitaine, grant nr. 2017 -1R50109 -00013434; the European Research Council (H2020 grant agreement number 715093, ECSTATIC); and the French National Research Agency, grant references ANR-10-IAHU04-LIRYC and ANR-11-EQPX-0030. This work was granted access to the HPC resources of CINES under the allocation 2021-A0090307379 made by GENCI.


Ihab S Mohamed
email: mohamedi@iu.edu

Guillaume Allibert
email: allibert@i3s.unice.fr

Philippe Martinet

Sampling-Based MPC for Constrained Vision Based Control

Ihab S. Mohamed 1 and Guillaume Allibert 2 and Philippe Martinet 3 Abstract-Visual servoing control schemes, such as Image-Based (IBVS), Pose Based (PBVS) or Hybrid-Based (HBVS) have been extensively developed over the last decades making possible their uses in a large number of applications. It is well-known that the main problems to be handled concern the presence of local minima or singularities, the visibility constraint, the joint limits, etc. Recently, Model Predictive Path Integral (MPPI) control algorithm has been developed for autonomous robot navigation tasks. In this paper, we propose a MPPI-VS framework applied for the control of a 6-DoF robot with 2D point, 3D point, and Pose Based Visual Servoing techniques. We performed intensive simulations under various operating conditions to show the potential advantages of the proposed control framework compared to the classical schemes. The effectiveness, the robustness and the capability in coping easily with the system constraints of the control framework are shown.

I. INTRODUCTION

Vision Based Control, also named Visual Servoing (VS), has been developed using directly 2D information extracted from images (IBVS), indirectly by computing 3D information from it (PBVS), or by combining both 2D and 3D information in the same scheme (HBVS). Nice state of the art have been published in [START_REF] Hutchinson | A tutorial on visual servo control[END_REF]- [START_REF] Chaumette | Visual Servoing[END_REF]. As referred in [START_REF] Chaumette | Visual Servoing[END_REF], it is still necessary to instantiate Visual servoing to dedicated applications, to develop new non linear control strategies (i.e multi sensor based control), and to explore new Direct Visual Servoing approaches by the use of a data-driven methodology. In addition, the study of stability, singularity locus, local minima of Visual Servoing scheme are still open problems to be addressed in deep. Many advanced control frameworks have been applied to Visual servoing, like General Predictive Control [START_REF] Gangloff | 6 dof high speed dynamic visual servoing using gpc controllers[END_REF], Model Predictive Control (MPC) [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF], Linear Quadratic Gaussian [START_REF] Hashimoto | LQ optimal and nonlinear approaches to visual servoing[END_REF], and Computed Torque Control [START_REF] Dahmouche | Dynamic visual servoing from sequential regions of interest acquisition[END_REF], to cite some.

A new algorithmic methodology based on Path Integral (PI) optimal control theory has been proposed by Kappen in [START_REF] Kappen | Path integrals and symmetry breaking for optimal control theory[END_REF], for solving the nonlinear Stochastic Optimal Control problem. Traditionally, based on dynamic programming, this problem is defined by a partial differential equation (PDE) known as the Hamilton-Jacobi-Bellman (HJB) equation which can not be solved analytically. In PI control theory, non-linear PDE can be transformed into an expectation over all possible trajectories using the Feynman-Kac (FK) France, {philippe.martinet}@inria.fr lemma. This transformation allows to solve the problem by sampling methods, such as forward-in-time Monte-Carlo approximation, instead of solving the HJB equation backwardin-time. Inspired by the PI control theory, Williams et al. [START_REF] Williams | Model predictive path integral control: From theory to parallel computation[END_REF] proposed a sampling-based model predictive control algorithm known as Model Predictive Path Integral (MPPI) control framework, which has been successfully applied to autonomous driving task in 2D environments. More recently, a generic and elegant MPPI control framework has been presented in [START_REF] Mohamed | Model predictive path integral control framework for partially observable navigation: A quadrotor case study[END_REF], which enables the robot to navigate autonomously in either 2D or 3D environments that are inherently uncertain and partially observable.

Although MPPI and MPC follow the same control strategy, we believe that the MPPI control framework significantly outperforms the conventional MPC strategy because it is a sampling-based and derivative-free optimization method, and it does neither require the computation of gradients (i.e., derivatives), and neither the first-or second-order approximation of the system dynamics and quadratic approximation of the cost functions. Additionally, discontinuous cost functions (i.e., indicator functions) can be easily handled and added to the running cost function. For instance, in the context of autonomous flying tasks, a large-weighted indicator function can be employed as part of the running cost, for penalizing the collision with ground or obstacles.

In this paper, the MPPI-VS framework based on PI control theory is presented. More precisely, we propose a realtime and inversion-free control method for all image-based, 3D point-based (3DVS), and position-based visual servo schemes, which has been validated on a 6-DoF Cartesian robot with an eye-in-hand camera. The paper is organized as follows. In section II, classical visual servoing concepts are summarized for 2D image point, 3D point, and Pose features. Section III described the proposed MPPI-VS framework. Section IV is dedicated to the results and discussion. The last section presents the conclusion and perspectives of the work.

II. CLASSICAL VISUAL SERVOING CONTROL SCHEMES

Classical Visual Servoing techniques are summarized in [START_REF] Chaumette | Visual servo control, part I: Basic approaches[END_REF]. They are based on the establishment of an interaction matrix L s that characterizes the evolution of a sensor feature s(t) regarding to the relative velocity between the sensor and the environment part from which the feature is extracted. The main objective of all vision-based control schemes is to minimize the error e(t) between the current visual features s(t) and the desired features s * , which is typically defined as e(t) = s(t) -s * . By imposing an exponential decrease of the error ( ė(t) = -λ s • e(t)), and considering an eye-in-hand configuration, we obtain a simple proportional control law:

v c = -λ s • L s + • (s(t) -s * ) (1
u; v]     - fu Z 0 u Z uv fv -fu -u 2 fu fu fv v 0 - fv Z v Z fv + v 2 fv -uv fu - fv fu u     [t; θu] c * Rc 03 03 L θu , L θu =I3-θ 2 [u]×+ 1-sinc θ sinc 2 θ 2 [u] 2 ×
For the estimation of the interaction matrix L s , we can consider different cases of L s # : at each iteration (#0), at each iteration with fixed depth (#1), at equilibrium (#2), and finally the half sum of cases #0 and #2 [START_REF] Malis | Improving vision-based control using efficient secondorder minimization techniques[END_REF].

Case

III. MPPI CONTROL STRATEGY FOR VISUAL SERVOING

In this section, we present the control strategy of our proposed sampling-based MPC approach (MPPI) for visual servoing systems; then, we state the mathematical formulation of MPPI in the presence of constraints such as the visibility, three-dimensional, and control constraints.

A. Review of MPPI

The MPPI control strategy is a sampling-based and derivative-free optimization method to MPC that can be easily applied in online to the real system, without requiring the first-or second-order approximation of the system dynamics and quadratic approximation of the objective functions. The control cycle of MPPI is described in Fig. 1. At each time-step ∆t, MPPI samples thousands of control trajectories from the system dynamics using a Graphics Processing Unit (GPU) to ensure a real-time implementation. Afterward, based on parallel nature of sampling, each of these trajectories is individually executed and then evaluated according to its expected cost. In sequential, the optimal control sequence U, over a finite prediction time-horizon t p ∈ {0, 1, 2, . . . , T -1}, is updated based on a weighted average cost over these generated trajectories, where U = (u 0 , u 1 , . . . , u T -1 ) ∈ R m×T , and T ∈ R + refers to the number of timesteps. Finally, the first control u 0 is applied to the system, while the remaining control sequence of length T -1 is slid-down to be used for providing a warm-starting to the optimization at the next time-step. Let δu t ∈ R m be a zero-mean Gaussian noise vector with a variance of Σ u , i.e., δu t ∼ N (0, Σ u ), where δu t represents the random noise associated with the commanded control input u t to the system. Then, the control input is defined by the exploring update v t = u t +δu t . Suppose that the number of the samples (namely, trajectories or rollouts) drawn from the discrete-time dynamics system, x t+1 = f (x t , v t ), is K, where x t ∈ R n denotes the state of the system at time t. Moreover, let S (τ t,k ) ∈ R + be the cost-to-go of the k th trajectory from time t onward. Then, based on the detailed derivation given in [START_REF] Williams | Model predictive path integral control: From theory to parallel computation[END_REF] and well-summarized in [START_REF] Mohamed | Model predictive path integral control framework for partially observable navigation: A quadrotor case study[END_REF], the optimal control sequence {u t } T -1 t=0 can be readily updated using the following iterative update law:

u t ← u t + K k=1 exp -(1/λ) S (τ t,k ) δu t,k K k=1 exp -(1/λ) S (τ t,k ) (2)
where λ ∈ R + is so-called the inverse temperature which determines the level of the selectiveness of the weighted average. We defined the cost-to-go of each trajectory τ over the predefined prediction time-horizon as:

S (τ ) = φ (x T ) + T -1 t=0 q (x t , u t , δu t ) (3)
in which φ (x T ) refers to the terminal cost. Whilst q (x t , u t , δu t ) denotes the instantaneous running cost, which composed of the sum of state-dependent running cost q (x t ) and quadratic control cost, and is defined as follows:

q = q (x t ) State-dep. + 1 -ν -1 2 δu T t Rδu t + u T t Rδu t + 1 2 u T t Ru t Quadratic Control Cost
where R ∈ R m×m is a positive definite control weight matrix, and ν ∈ R + is so-called the exploration noise which determines how aggressively MPPI explores the state-space. The impact of changing ν has been studied in [START_REF] Mohamed | Model predictive path integral control framework for partially observable navigation: A quadrotor case study[END_REF].

B. MPPI-VS Control Strategy

From the relation ṡ(t) = L s • v c , a discrete-time model can be approximated using the Newton-Euler method as

s(t + 1) = s(t) + L s (t)v c (t)∆t (4)
where s(t) defines the state vector of the model, composed of a set of n p 2D image points in MPPI-IBVS, a set of n p 3D-points in MPPI-3DVS, or by a pose vector in MPPI-PBVS. Similarly, the approximated interaction matrix is given either by L s = L PIm ∈ R 2np×6 , L P ∈ R 3np×6 or L P os ∈ R 6×6 for MPPI-{IBVS,3DVS,PBVS} control schemes. MPPI-VS control strategy directly uses the estimation of the interaction matrix. In CASE #0 and #3 of the approximate L PIm where the depth Z i of each point-like feature needs to be estimated, we use the discrete-time model of the 3D-point features given in (4) for predicting the future evolution of the object depths Z i (t), assuming that, at each time-step ∆t, the initial depth Z i (0) of a set of n p point-like features is known.

1) Handling Visibility and 3D Constraints: One of the most attractive features of MPPI, compared to the classical MPC, is its capability of coping easily with the hard and soft constraints. More precisely, a large-weighted indicator function can be employed as part of the state-dependent running cost function q(x) for handling the constraints. In MPPI-VS control scheme, the instantaneous state-dependent running cost function, q(x) ≡ q(s), is defined as

q(s t ) = q 1 (s t ) + q 2 (s t ) (5)
where q 1 (s t ) = (s t -s * ) T Q(s t -s * ) denotes the statedependent cost which is a simple quadratic cost for enforcing the current state s t to reach its desired value s * . q 2 (s t refers to an indicator function used for handling the visibility constraints (expressed as s min ≤ s t ≤ s max , where s min = [u min ; v min ] and s max = [u max ; v max ] refer to the lower and upper ranges, in pixels, of the image point coordinates) and the three-dimensional constraints such as workspace or joint limits. The definition of q 2 (s) depends on the control scheme.

In MPPI-IBVS control scheme, q 2 is formulated as q 2 (s k ) = 10 7 C 1 +10 5 C 2 , where C 1 and C 2 are Boolean variables that are used to heavily penalize each trajectory that violates the visibility and 3D constraints, respectively. For instance, C 1 is turned on (i.e., C 1 = 1) if, at any time, s exceeds its given bounds, i.e., C 1 = 1 ⇔ ∃s| (s > s max ) ∨ (s < s min ) . We can also add the second term, if needed, for constraining the future evolution of the 3D-point features, namely,

C 2 = 1 ⇔ ∃P i | (P i > P max ) ∨ (P i < P min ) ∀i = 1, . . . , n p .
To name just a few, this term can be added to ensure that the camera remains in the workspace, simply by limiting its position in the workspace. The special case of camera retreat, present in large rotation around the optical axis, can easily be taken into account. C 2 will be active if at any time the Z position of the camera exceeds a threshold Z max , i.e.,

C 2 = 1 ⇔ ∃Z i | (Z i > Z max ).
In MPPI-PBVS control scheme, the controller might produce input that would ultimately lead to that some visual features may leave the camera's Field of View (FoV). This is mainly because the control law is explicitly expressed in Cartesian space and there is no direct control to the visual features on the image plane. Therefore, we propose two methods to guarantee that the object always remains within the camera's FoV. Concerning the first proposed method, the indicator function q 2 is defined as a large-weighted exponential penalty function, as follows:

q 2 = β np i=1 e -α xmax-|xi| + e -α ymax-|yi| (6)
so that the visibility constraints are constantly satisfied which can be written as |p i | ≤ p max ∀i = 1, . . . , n p , where

p i = [x i ; y i ] = [ Xi Zi ; Yi Zi ] ∈ R 2
denotes the normalized Cartesian coordinates of an image point, p max = [x max ; y max ], p min = -p max are the minimum and maximum bounds, in meters, of the point feature projected on the normalized image plane. β and α are positive scalar variables. In the second method, we proposed an alternative solution based on augmenting the 3D-point features within the nominal state vector. In further words, we combine both MPPI-PBVS and MPPI-3DVS control schemes. Thus, in this case, we have s = [P os; ×6 . In the present scheme, the penalty function q 2 (s) is replaced by a large-weighted quadratic state-dependent cost function which modulates how fast the current state vector of the 3D-point features P converges to its desired one P * . Therefore, the running cost function q(s) is reformulated as

{P i } i ] = c * t c ; θu; {X i ; Y i ; Z i } i ∈ R (6+3np) , whilst L s = L P os ; L P ∈ R (6+3np)
q(s) = w 1 5 i=0 (s i -s * i ) 2 :=q1(s), si∈ P os, s * i ∈ P os * + w 2 n-1 i=6 (s i -s * i ) 2 :=q2(s), si∈ P , s * i ∈ P * (7)
where s * = [P os * ; P * ], n = 3n p + 6 is the length of the state vector s, whereas w 1 , w 2 ∈ R + refer to the cost weighting of q 1 and q 2 , respectively. It should be noted that assigning very high weights, particularly for w 2 , is sufficient for enforcing the visual features to stay within the image plane, without imposing additional constraints.

2) Handling Control Constraints: Most robotic systems have constraints on the actuators (range or velocity limitations). These constraints are known as control constraints which are considered to be hard constraints that should not be violated. In MPPI control strategy, there exists two ways for handling the control constraints. First, they can be implemented as a natural part of the running cost by adding an appropriate term to penalize all the trajectories that violate the constraints. However, the common issue associated with this way is that the control constraints acting as soft constraints, and, accordingly, it is difficult to ensure that the control input always remains within its allowed bounds. Therefore, we propose an alternative method based on pushing the control constraints into the dynamics system [START_REF] Williams | Information theoretic MPC for model-based reinforcement learning[END_REF], meaning that x t+1 = f x t , g(v t ) , where g(v t ) is an element-wise clamping function that is used to restrict the control input v t ∼ N (u t , Σ u ) to remain within a given range, for all samples drawn from the dynamics system. Thus, g(v) can be defined as g(v) = max v min , min(v, v max ) , where v min and v max are the lower and upper bounds of the control input. The major advantage is that the problem formulation of MPPI, which takes into account control constraints, is converted into an unconstrained one, without violating the constraints and affecting the convergence of the MPPI algorithm as g(v) has an impact only on the dynamics system.

IV. SIMULATION-BASED EVALUATION

We conduct extensive simulations considering a 6-DoF Cartesian robot with an eye-in-hand camera configuration to evaluate and demonstrate the potential advantages of our proposed control strategy with a comparison to the classical schemes (C-IBVS, C-3DVS, and C-PBVS).

A. Simulation Settings and Performance Metrics

We validated the MPPI-VS control strategy as well as the classical schemes using four points (i.e., n p = 4) in the image plane. These four dots, forming a square in the xy-plane, have the coordinates in the reference frame . During all the simulations, it is assumed that the MPPI algorithm runs with a time horizon t p of 3.5 s, a control frequency of 50 Hz (i.e., T = 175), and generates 2000 samples each timestep ∆t with an exploration variance ν of 1000. Moreover, the control weighting matrix R is set to 1 2 λΣ -1 u , assuming that the random noise associated with the control input has a variance of Σ u = Diag (0.02, 0.01, 0.01, 0.02, 0.02, 0.01). The remaining parameters are defined in Table I. The inverse temperature λ is set to a high value in MPPI-IBVS scheme, and low values in both MPPI-3DVS and MPPI-PBVS control schemes. It is important to bear in mind that fine-tuning λ is mainly based on the statedependent running cost function q (s). The real-time MPPI-VS algorithm is executed on an NVIDIA GeForce GTX 1080 Ti desktop GPU. All the control schemes were developed using Visual Servoing Platform (ViSP) [START_REF] Marchand | Visp for visual servoing: a generic software platform with a wide class of robot control skills[END_REF] integrated with the Robot Operating System framework.

1) Initial Camera Configurations: 120 initial camera configurations were randomly extracted from a uniform distribution within the robot's workspace, with a guarantee that the robot kinematics has initially the capability of reaching these generated poses, and the visual features are initially located within the camera's FoV.

2) Desired Camera Configurations: We considered two different desired poses of the camera, all expressed in the reference frame F o attached to the object with z-axis pointing downward. The first desired camera pose is given by o P os * c1 = [ o t c ; θu] = [0, 0, -0.75, 0, 0, 0] T in ([m], [deg]). We have also defined the second camera desired configuration o P os * c2 = [0.076, 0.202, -0.727, 10, -10, -15] T in ([m], [deg]) in such a way that the object lies in the top-right corner of the image, to demonstrate how robust the proposed control strategy is against violating the system constraints particularly visibility constraints.

3) Performance Metrics: The visual servoing task is considered to be successful if the following metrics are satisfied: R LM : The camera does not reach a local minimum (i.e., R LM = 0) during the simulations. Such configuration corresponding to a local minimum occurs when v c = 0 and o P os c = o P os * c . Mathematically, we used the meansquared error (MSE) as a metric that measures the positioning error , which can be formulated as MSE = 1 n e T e in which e = o P os c -o P os * c and n = 6. Thence, we considered that the local minimum is completely avoided if and only if

MSE 1 = 2 n e T 1 e 1 < 1 and MSE 2 = 2 n e T 2 e 2 < 2
, where e 1 (in [m]) and e 2 (in [rad]) refer to the translational and rotational errors between o P os c and o P os * c , respectively. During the simulations, we set the two thresholds as 1 = 10 -5 and 2 = 10 -4 . R JL : The robot joint limits are avoided (i.e., R JL = 0). If the robot reaches one of the joint limits, the velocity of each joint will be set to zero by the low level robot controller. P out : The considered visual features always remain within the camera's FoV (i.e., P out = 0). Finally, the number of successful servoing tasks, with respect to the total initial camera configurations, is denoted as N success , while S rate indicates the success rate.

B. Intensive Simulation Results

Six different tables (II:VII) summarize the overall performance of each control scheme individually, given the 120 initial camera configurations. In total, we conduct a set of 24 different tests to evaluate and compare the performance of MPPI-VS control scheme with classical visual servoing. The robustness, the parameters influence and the capacity to handle different constraints are also studied. 1) Prediction Process Influence and Comparison with C-IBVS: In the four first tests (see Table II), we analyzed the impact of the image prediction process on the performance of the MPPI-IBVS control strategy by taking into account the four cases of the approximate interaction matrix L PIm . It can be clearly seen that the prediction process exerts a high influence as MPPI provides better performance in CASE #0 and #3, in which the 2D visual features and their depth information are assumed to be estimated each iteration. As expected, the worst performance is achieved when L PIm is assumed to be constant, as depicted in Test 3.

More accurately, for Test 1, we conducted 3 trials. For all the trials, our proposed control scheme succeeded in completing 119 out of 120 tasks/configurations, while there existed only one initial camera configuration led the robot to reach one of its joint limits. This failure case can be easily tackled by increasing the prediction horizon as illustrated in Test 5 where we adopted t p = 5 s and K = 1000 instead of t p = 3.5 s and K = 2000. In order to illustrate the superiority of MPPI-IBVS, intensive simulations of the C-IBVS control scheme are carried out in Test 6 and 7 in which we set λ s to 0.5 and 0.2, respectively. As anticipated, the success rate of visual servoing tasks was reduced to 70%, with 36-R JL failure cases. In fact, these failure cases occurred due to the fact that their corresponding initial camera configurations have large rotations around the camera optical axis as well as their positions along the z-axis are much closer to the minimum allowable limit of the third joint.

2) MPPI Parameters Influence and Robustness: In Table III is presented the influence of MPPI parameters as well as the robustness to modeling errors. From Test First, a white noise generated from a uniform distribution is added to the visual features and their 3D information (i.e., depth estimation) with a maximum error of ±1 pixel for 2D features and ±0.5 cm for 3D points. In addition, errors on the intrinsic parameters of the camera are also added: +30 % in f , ±20 % in ρ u and ρ v (with f u = f ρu and f v = f ρv ), ±15 % in u 0 and v 0 . As can be seen, our proposed control scheme behaves perfectly as there exists a maximum of 4-R JL failure cases that occurred, demonstrating its robustness.

3) Handling Constraints: In Table IV are presented the results when handling constraints. In Test 12, we repeated the simulations of Test 1 under the assumption that the maximum control input of the camera velocity screw v max is limited to 0.5 m/s for the translational speed and 0.3 rad/s for the rotational speed, while the minimum control input v min equals to -v max . The obtained results of the 105 successful tasks demonstrate that our proposed control scheme performs perfectly with a high capability of handling the control constraints. Note that the 15-R JL failure cases can be easily avoided and tackled by involving the 3D constraints, especially joint limits constraints, as a part of the MPPI-IBVS optimization problem. In all tests where MPPI-IBVS is utilized, the simulations are carried out taking into account the visibility constraints, which are defined by the following inequalities:

u min = 0 v min = 0 ≤ s t ≤ u max = 640 v max = 480 , in [pixels]. (8)
Additionally, to better evaluate the capability of handling the constraints, the simulations of Test Tests 19 to 21, simulations are conducted to verify the performance of our proposed PBVS control scheme. We studied the behavior of MPPI-PBVS without applying the visibility constraints, together with the performance of the classical scheme. The obtained results illustrate that both control strategies produce control input ultimately leading to that some visual features leave the camera's FoV (i.e., P out = 0). This result was expected since no position-based controllers (classical or MPPI-based) provide any guarantee that the projections of the visual features remain within the camera's FoV. This is especially true in the case of MPPI due to the stochastic nature of the control sequence. [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF] provides better motion in the image space while the 3D camera trajectory obtained by applying ( 6) is considerably shorter. It is noteworthy that the trajectories obtained by [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF], particularly in the image plane, and satisfying the visibility constraints are mainly based on the assigned value to w 2 . The results show that higher this value is, the better is the performance.

V. CONCLUSIONS AND FUTURE WORK

In this work, we have proposed a real-time samplingbased MPC approach for predicting the future behavior of the VS system, without solving the online optimization problem which usually exceeds the real system-sampling time and suffers from the computational burden. There is no need for estimating the interaction matrix inversion or performing the pseudo-inversion in real-time. The proposed approach directly uses the approximate interaction matrix, i.e., it is an inversion-free control method. The visibility, three-dimensional (i.e., 3D), and control constraints as well as parametric uncertainties can be easily handled. Our next step will be the validation on a real Cartesian robot.

)

  In this expression, v c = [v c ; ω c ] is the kinematic screw applied to the camera, with the translational velocity v c = [v x ; v y ; v z ] and the rotational velocity ω c = [ω x ; ω y ; ω z ]. Let's define a 3D point P = [X; Y ; Z] whose image projection is P Im = [u; v]. Let's also define the rotation matrix between the desired camera frame and the current camera frame c * R c and the pose (position, orientation) between this two frames as P os = [t; θu]. The following

Fig. 1 :

 1 Fig. 1: MPPI control loop.

  F o : P 1 = [-0.1; -0.1; 0], P 2 = [0.1; -0.1; 0], P 3 = [0.1; 0.1; 0], and P 4 = [-0.1; 0.1; 0] [m]

  Table givesthe expressions of the corresponding interaction matrices (L P , L PIm , L P os ). Γ represents the intrinsic parameters (f u , f v , u 0 , v 0 ) of a pinhole camera, and I 3 ∈ R 3×3 is the identity matrix.

	s(t)	Ls (L

P , L P Im , L P os ) [X; Y ; Z] -I3 P , P is the skew matrix of P [

TABLE I :

 I MPPI-VS Parameters

	Scheme	Parameter	Value	Parameter	Value
	MPPI-IBVS	λ	100	Q	2.5I 8
	MPPI-3DVS	λ	10 -2	Q	35I 12
	MPPI-PBVS	λ β, α (6)	10 -3 150, 10 3	Q w 1 , w 2 (7)	35I 6 35, 150

TABLE II :

 II Prediction Influence and Comparison

	Test No. Control Scheme RLM Pout RJL Nsuccess Srate
	Test 1	MPPI (CASE #0) 0	0	1	119 99.16%
	Test 2	MPPI (CASE #1) 0	0 29	91	75.8%
	Test 3	MPPI (CASE #2) 0	0 54	66	55%
	Test 4	MPPI (CASE #3) 0	0	6	114	95%
	Test 5	MPPI (CASE #0) 0	0	0	120	100%
	Test 6 Classic (CASE #0) 0	0 37	83 69.17%
	Test 7 Classic (CASE #0) 0	0 36	84	70%

TABLE III :

 III MPPI Parameters Influence and Robustness , the influence of both control input updates Σ u , prediction horizon t p , and number of sampled trajectories K is individually studied. In Test 8, we adopted Σ u = Diag (0.009, 0.009, 0.009, 0.03, 0.03, 0.02), while in Test 9 t p is set to 1 s. Finally, we set K to 100 in Test 10. In general, the intensive simulations demonstrate that the impact of either having short-time horizons t p or changing the control input updates Σ u is appreciably higher than the influence of decreasing the number of samples K, as the success rate S rate in Tests 8 and 9 is significantly lower than that in Test 10. Concerning the quality of successful servoing tasks, it is interesting to notice that assigning very low values to t p and K affects the convergence rate 1 and quality of the trajectories in both the image and 3D space. The robustness of the MPPI-IBVS is then evaluated with respect to the camera modeling errors and measurement noise. Due to the lack of space, only the test combining all the errors is given (seeTest 11).

	Test No. Control Scheme RLM Pout RJL Nsuccess Srate
	Test 8 MPPI (CASE #0) 0	0 19 101 84.17%
	Test 9 MPPI (CASE #0) 0	0 35	85	70.8%
	Test 10 MPPI (CASE #0) 0	0	7	113 94.17%
	Test 11 MPPI (CASE #0) 0	0	4	116 96.67%
	8 to 10

TABLE IV :

 IV Handling Control and Visual Constraints

	Test No. Control Scheme RLM Pout RJL Nsuccess Srate
	Test 12 MPPI (CASE #0) 0	0 15 105 87.5%
	Test 13 MPPI (CASE #0) 0	0 25	95 79.17%
	Test 14 Classic (CASE #0) 0	0 42	78	65%

  1 are repeated in Test 13, considering o P os * c2 as the desired configuration instead of o P os * c1 . Moreover, in Test 14, we repeated the simulations of C-IBVS given in Test 6 with respect to o P os * c2 . The intensive simulations presented in Table IV demonstrate the efficiency and capability of our proposed control strategy in coping easily with the visibility constraints, as the visual features always remain within the camera's FoV (i.e., P out = 0 for all tests). 4) 3DVS Control Schemes: In Table V, the results of MPPI-3DVS are presented. Four intensive simulations are

TABLE V :

 V 3DVS Performance 3DVS and C-3DVS (with λ s = 0.5), considering the two desired camera configurations and the MPPI-3DVS parameters listed in TableI. As excepted when 3D feature points are considered, the simulations show that both control strategies give very satisfactory results, where the maximum failure cases were 7 (6 for R JL and only one for P out ) occurred in Test 18.5) Unconstrained MPPI-PBVS V.S. C-PBVS: TableVIshows the results of MPPI-PBVS without constraint. In the

	both MPPI-
	Test No. Control Scheme RLM Pout RJL Nsuccess Srate
	Test 15 MPPI ( o P os * c 1 ) 0 Test 16 Classic ( o P os * c 1 ) 0 Test 17 MPPI ( o P os * c 2 ) 0 Test 18 Classic ( o P os * c 2 ) 0	0 0 0 1	0 4 5 6	120 116 96.67% 100% 115 95.83% 113 94.17%
	carried out in Tests 15 to 18 for assessing the performance of

TABLE VI :

 VI PBVS Performance in Unconstrained Case

	Test No. Control Scheme RLM Pout RJL Nsuccess Srate
	Test 19 MPPI ( o P os * c 1 ) 0 Test 20 Classic ( o P os * c 1 ) 0 Test 21 Classic ( o P os * c 2 ) 0	45 0 25 0 21 0	75 62.50% 95 79.17% 99 82.50%

TABLE VII :

 VII 6) Constrained MPPI-PBVS: In Table VII are presented the results of MPPI-PBVS with constraints. To take into PBVS Performance in Constrained Case the difficulties mentioned above, we propose to handle visual constraints in MPPI-PBVS schemes. The two methods proposed in (6) and (7) to cope with the visibility constraints are used for the desired configuration o P os * c1 (see Tests 22 and 24). Results show the validity and effectiveness of our proposed methods in handling visibility constraints, regardless of which desired camera configuration is considered (Test 22 and 23), as there exists neither R JL nor P out failure cases. It is also observed that despite MPPI-PBVS on the basis of

	Test No.	Control Scheme RLM Pout RJL Nsuccess Srate
	Test 22 MPPI ( o P os * c 1 , (6)) 0 Test 23 MPPI ( o P os * c 2 , (6)) 0 Test 24 MPPI ( o P os * c 1 , (7)) 0	0 0 0	0 0 0	120 100% 120 100% 120 100%

account

Note that the simulation time of Test 9 is adopted to be 180 s as the convergence rate of MPPI, in this case, is extremely slow.


M Bacal

P Svarnas

S Béchu

J Pelletier

Study of the Plasma near the Plasma Electrode by Probes and Photodetachment in ECR-driven Negative Ion Source

Keywords: Negative ions, plasma electrode bias, magnetic filter, extracted beams PACS : 32.80.Gc, 52.27.Cm, 29.25.Ni, 52.80.Pi 9 cm -

The effect of the plasma electrode bias on the plasma characteristics near the extraction aperture in a large volume hybrid multicusp negative ion source, driven by 2.45 GHz microwaves, is reported. Spatially resolved negative ion and electron density measurements were performed under various pressures (1-4 mTorr) by means of electrostatic probe and photodetachment technique.

INTRODUCTION

Contemporary negative ion sources operate with a magnetic filter field extending up to the plasma electrode (PE) which contains the extraction aperture 1 . It is generally accepted that a positive PE bias (V PE ) reduces the extracted electron current. In some experiments, an optimum V PE for negative ion extraction was found [2][3][4] . A possible cause of these effects was suggested in 3 : electrons are trapped in the weak magnetic field and lost along the field lines, whereas positive and negative ions are not trapped because of their large Larmor radius. More negative ions arrive from the bulk plasma to ensure plasma neutrality. However, a definite explanation of the H -ion current peak, when the PE bias is varied, is still not available.

The electrostatic particle simulation 1,[START_REF] Sakurabayashi | Conference Proceedings[END_REF][START_REF] Hatayama | [END_REF] without applied PE bias showed that the presence of the weak magnetic field produces important modifications in the positive ion flow and as a result in the structure of the plasma potential. A characteristic peak in the plasma potential (V p ) dependence on the distance from the PE (D) is predicted, where many negative ions are localized.

The purpose of this work is the experimental study of the spatial distribution of the plasma parameters in the vicinity of the PE, specifically on the axis of the extraction aperture, using electrostatic probes and laser photodetachment.

EXPERIMENTAL SETUP and DIAGNOSTICS

The H -ion source Camembert III has been described in 4 and a draft is shown in Fig. 1. The magnetic multipole chamber Camembert III was previously equipped with 16 filaments, installed on the top flange. It was operating in a dc mode, similar to a filamented ion source. Since 2003, a 2-D network of seven elementary independent ECR plasma sources 7 has been installed on the top flange. They are powered by microwaves (2.45 GHz) at 0.9 kW, equally distributed at the seven elementary sources. Each source contains a permanent magnet which provides the magnetic field required for ECR (875 G). This magnetic field also confines the fast electrons. It represents a first magnetic filter and the region around the magnet represents the "driven region". Opposite to the top flange there are the PE with the extraction hole, followed by the extraction electrode equipped with a pair of magnets for electron separation, and the acceleration electrodes. The stray magnetic field from the magnets located in the extraction electrode extends into the plasma and provides a second magnetic filter. This magnetic field attains a maximum of 20 G at a distance of 1.6 cm from the PE. This region can be denoted as the "extraction region". A large volume of cold plasma separates the above two extreme regions.

The various plasma parameters were measured with an electrostatic probe made of tungsten wire of 0.5 mm diameter, mounted on a three-axis manipulator installed on a flange at the PE level. A tip of 15 mm was exposed to the plasma, while the rest was shielded by a ceramic cylinder. This probe was L-shaped so that its longer part was coaxial with the laser beam, used for measuring the negative ion density by photodetachment 8 . The measurements were effected on the axis passing through the center of the extraction opening.

A Nd-YAG laser (repetition rate 10 Hz, photon energy 1.2 eV) was used to illuminate the plasma, through a 6 mm diaphragm, and the photodetached electrons from the H -ions were collected by the L-shaped, positively biased, probe. The probe transient current, generated by the collected photodetached electrons, was recorded on a digital oscilloscope (500 MHz, 2 GSample/s) operating in averaging mode, via capacitive decoupling. Since previous studies showed that in the region near the PE the negative-ion relative density n H -/n e is much higher than in the center of the source, the analysis of the probe signals was carried out according to the theory and recommendations of 9 .

SPATIAL VARIATION of the PLASMA PARAMETERS

The results obtained at a pressure p of 3 mTorr are presented in Fig. 2 to 6. The electron density spatial variation is shown in Figure 2. The electron density goes down linearly with the distance to the PE. The higher the positive bias V PE , the lower the electron density. This shows that the electrons are efficiently swept away from the PE neighborhood due to the applied positive PE bias. The spatial variation of the negative ion density (Fig. 3) clearly indicates its enhancement at a distance close to the plasma electrode where the electron density tends to be reduced. The optimum PE bias is V PE = 5.5 V. V PE = 4 V

V PE = 5 V V PE = 6 V V PE = 7 V
electron density (10

V PE = 4 V V PE = 5 V V PE = 6 V V PE = 7 V H - density (10
Fig. 4 shows the variation of the relative negative ion density n H -/n e versus D. Note that this relative density goes up with the PE bias and attains the value of 1.9 for the highest V PE . An indication of a weak plasma potential peak versus D exists for V PE > 4 V. This "maximum" is located at a distance of 1 cm (Fig. 5). The spatial variation of the electron temperature (Fig. 6) indicates that for V PE < 6 V the electron temperature is low (0.3 -0.5 eV) at all the distances studied, but goes up abruptly near the PE when V PE > 6 V.

V PE = 4 V V PE = 5 V V PE = 6 V V PE = 7 V
V PE = 4 V V PE = 5 V V PE = 6 V V PE = 7 V plasma potential (V)
distance from the PE (mm) FIGURE 6. Variation of the electron temperature with distance from the plasma electrode, for four values of the plasma electrode bias V PE .

A comparison on the pressure dependence of the negative ion densities at 150 mm and 7.5 mm from the PE shows that the latter is enhanced particularly at the lowest pressure (Fig. 7). Even at a distance of 20 mm the effect of the magnetic field and the PE bias starts to go down, as can be seen in Fig. 3.

FIGURE 7.

Comparison between the negative ion density close to the extractor aperture (solid symbols) and at the center of the source (open symbols) vs. pressure. The extracted negative ion current is also compared (dotted line). 4 The extracted negative ion current I H -is correlated better with the plasma parameters near the extraction opening (Fig. 7) than with those in the center of the source (150 mm from the PE). The optimum pressure at the center of the source lies between 4 and 5 mTorr, as opposed to 1.5-3 mTorr at 7.5 mm from the PE. The latter approaches the optimum pressure for the extracted H -ion current. Fig. 8 shows the ratio n H -/n e between the negative-ion density and the electron density found at a distance of 7.5 mm from the PE when V PE is varied. In the same figure the ratio I H -/I e of the extracted currents is as well presented.

V PE = 4 V V PE = 5 V V PE = 6 V V PE = 7 V

CONCLUSION

Compared to the electron density at the center of the source (3x10 10 cm -3 at 3 mTorr), the electron density at 1 cm from the plasma electrode is reduced by a factor 30, when a PE bias V PE = 6 V is applied.

The plasma potential varies abruptly near the PE for extreme bias values, when ion or electron sheaths are formed, but stays approximately constant for the medium value V PE = 5.5 V, which is close to the plasma potential far from the PE. Let us denote this value V PE * . This value of the PE bias -V PE * -is important since the maximum negative ion density occurs at this bias (not shown in Fig. 3). The electron temperature values increase abruptly when V PE > V PE * (see Fig. 6). The highest negative ion densities are found in the neighborhood of the plasma electrode, at distances where the electron density is reduced .

Note that the ratio n H -/n e attains the value 1.9 (Fig. 4), which means that the negative ions constitute the majority of the negative species. We should remind that in the filamented version of Camembert III values of n H -/n e as high as 12 were measured in a discharge of 2.5 kW 9 . It was observed that the maximum n H -/n e ratio went up when the plasma density or discharge power was enhanced. This may explain the lower n H -/n e value (i.e. ~2) found in this, ECR driven, version of CAMEMBERT III, where the plasma density is limited by the available microwave power.

In conclusion the explanation of the H -ion density peak can be summarized as follows 4 : the depletion of the electron population, when the positive bias of the PE is enhanced, disturbs the local plasma neutrality and new negatively charged particles should replace the electrons collected by the PE. In a collisionless situation, the electrons present in the bulk plasma are unable to cross the transverse magnetic field in front of the PE. Therefore, negative ions from the central plasma replace the depleted electrons in the magnetized region, as suggested in 3 . Indeed, the measured H - ion density increases close to the PE and peaks when V PE is slightly greater than V PE * . This results in higher H -ion extracted current. When the PE bias is further increased the H -ions are accelerated toward the PE. The negative-ion residence time decreases

FIGURE 1 .

 1 FIGURE 1. Schematic representation of the ECR-driven hybrid multucusp source Camembert III (a) and the magnets forming the multicusp configuration (b).

FIGURE 2 .

 2 FIGURE 2. Variation of electron density with distance from the plasma electrode, for four values of the plasma electrode bias V PE .

FIGURE 3 .

 3 FIGURE 3. Variation of the negative ion density with distance from the plasma electrode, for four values of the plasma electrode bias V PE .

FIGURE 4 .

 4 FIGURE 4. Variation of the relative negative ion density with distance from the plasma electrode, for four values of the plasma electrode bias V PE .

FIGURE 5 .

 5 FIGURE 5. Variation of the plasma potential with distance from the plasma electrode, for four values of the plasma electrode bias V PE .

FIGURE 9 .

 9 FIGURE 9. Ratios of the negative ion to electron density (dashed line) and the extracted negative ion to electron current (solid line) vs. plasma electrode bias.4

  / n e 7.5 mm from the PEI H _ / I e

due to the larger velocity and the H -density near the PE as well as the extracted negative-ion current go down.