Aurélien Dumenil
email: aurelien.dumenil@gmail.com

Adrien Kaladji
email: kaladrien@hotmail.fr

Miguel Castro
email: miguel.castro@univ-rennes1.fr

Simon Esneault
email: simon.esneault@therenva.com

Antoine Lucas
email: antoine.lucas@chu-rennes.fr

Michel Rochette
email: michel.rochette@ansys.com

Cemil Goksu
email: cemil.goksu@therenva.com

Pascal Haigron
email: pascal.haigron@univ-rennes1.fr

Finite-element-based matching of pre-and intraoperative data for image-guided endovascular aneurysm repair

Keywords: FEM simulation, tool/tissue interactions, EVAR, endovascular navigation

come

I. INTRODUCTION

BDOMINAL aortic aneurysm (AAA) is the most common type of aneurysm. The endovascular repair of these abdominal aortic aneurysms (EVAR) is a wellestablished technique throughout the medical and surgical communities. A real alternative to open surgery, it in volves deploying one or several endovascularly-inserted stent grafts at the anchoring site in order to repair the aneurysm. This procedure has the advantage of reducing blood loss, intraoperative morbidity, and hospital stay duration [START_REF]Endovascular versus Open Repair of Abdominal Aortic Aneurysm[END_REF], [START_REF]Evaluation des endoprothèses dans le traitement des anévrismes et des dissections de l'aorte thoracique[END_REF]. Its minimally-invasive nature has led to a s ignificant decrease in short-term postoperative mortality. Although increasingly indicated, this technique has some intrinsic limitations, notably due to the devices, imaging techniques, and protocols used, which may result in failures, conversion procedures or reinterventions (repeat procedures).

In order for the stent graft to reach its anchoring site, several guidewires are first inserted into the vascular structure. An extra-stiff guidewire is inserted into the iliac arteries to straighten them. These deformations are difficult to anticipate, as they depend on several factors, including vascular morphology, arterial wall rigidity, calcification status, etc.

Because these interventions are commonly performed under fluoroscopic control, two-dimensional (2D) visualization of the aneurysm can only be achieved via the injection of a contrast agent. In addition to its potential for iatrogenic effects, this projective imaging modality inevitably leads to topographic errors and gives no information on arterial wall quality. Arterial wall quality, however, is a determining factor in stent graft impermeability and stability at the time of deployment [START_REF] Becquemin | Chirurgie Endovasculaire des Anévrysmes de l'Aorte Abdominale[END_REF]. A specially-adapted intraoperative navigation interface could increase deployment accuracy and reveal such information, which preoperative three-dimensional (3D) imaging might otherwise provide.

This study addresses the issue of pre-and intra-operative imaging data matching in the EVAR context, and takes anatomical deformations into account in order to propose a reliable computer-aided endovascular navigation solution. The proposed approach involves the estimation of vascular deformations using the matching of a numerical simulation, Finite element-based matching of pre-and intraoperative data for image-guided endovascular aneurysm repair

Aurélien Duménil, Adrien Kaladji, Miguel Castro, Simon Esneault, Antoine Lucas, Michel Rochette, Cemil Göksu, and Pascal Haigron A derived from preoperatively acquired patient-specific data, and intra-operative data. Few studies on computer-aided endovascular navigation are reported in the literature. Augmented reality in the field of vascular surgery can be achieved by transferring preoperative CT volume data into the intraoperative environment using a registration process [START_REF] Markelj | A review of 3D/2D registration methods for image-guided interventions[END_REF]. Penney et al. proposed 3D/2D intensity-based rigid registration using the vertebra as regions of interest, but registration accuracy was deemed insufficient [START_REF] Penney | Validation of a two-to threedimensional registration algorithm for aligning preoperative CT images and intraoperative fluoroscopy images[END_REF]. This method was then improved and they reported the results from use of their image-guided surgery system during 23 procedures [START_REF] Penney | An Image-Guided Surgery System to Aid Endovascular Treatment of Complex Aortic Aneurysms: Description and Initial Clinical Experience[END_REF]. The method was within the target accuracy of 3mm in 78% of cases. Göksu et al. also proposed featurebased rigid registration [START_REF] Goksu | Endovascular navigation based on real/virtual environments cooperation for computer-assisted TEAM procedures[END_REF]. The features considered were 3D pre-operative vessel centerlines, extracted from pre-operative CTA, and 2D vessel centerlines, obtained from 2D intraoperative angiography. The navigation system provided good accuracy when the insertion of endovascular devices caused no deformity. Penney et al. also proposed a non-rigid registration method capable of deforming the preoperative model locally at the ostia for branched stent graft implantation [START_REF] Raheem | Nonrigid 2D-3D image registration for use in Endovascular repair of Abdominal Aortic Aneurysms[END_REF]. This study, however, focused on the aortic segment, located above the aneurysm (celiacomesenteric segment), rather than on the iliac arteries. The proposed method was based on geometrical considerations, with no attention paid to mechanical behavior law.

Finite element-based calculation makes it possible to assess deformations by taking into account mechanical considerations. Several studies have addressed the digital simulation of aneurismal wall behavior. These studies were aimed at analyzing wall stress distribution and risk of aneurysm rupture in models reconstructed from CT data [START_REF] Doyle | A comparison of modelling techniques for computing wall stress in abdominal aortic aneurysms[END_REF], [START_REF] Doyle | A finite element analysis rupture index (FEARI) as an additional tool for abdominal aortic aneurysm rupture prediction[END_REF]. By taking into consideration patient-specific data, some of these studies explored the impact of geometry [START_REF] Xenos | The effect of angulation in abdominal aortic aneurysms: fluid-structure interaction simulations of idealized geometries[END_REF], [START_REF] Doyle | Vessel asymmetry as an additional diagnostic tool in the assessment of abdominal aortic aneurysms[END_REF] and aneurismal wall calcification on wall stress distribution [START_REF] Speelman | Effects of Wall Calcifications in Patient-Specific Wall Stress Analyses of Abdominal Aortic Aneurysms[END_REF]. Blood flow characteristics before and after stent graft implantation were compared using fluid mechanics simulations [START_REF] Li | Blood flow and structure interactions in a stented abdominal aortic aneurysm model[END_REF], [START_REF] Molony | Fluid-structure interaction of a patient-specific abdominal aortic aneurysm treated with an endovascular stent-graft[END_REF], [START_REF] Howell | Computational Fluid Dynamics Within Bifurcated Abdominal Aortic Stent-Grafts[END_REF]. To our knowledge, only one study focused on the deployment of a complete aortic stent graft, with geometrical characteristics taken from a realistic phantom [START_REF] Bock | Virtual evaluation of stent graft deployment: A validated modeling and simulation study[END_REF]. The vascular structure deformations caused by the progress of the stiff guidewire have also been assessed via finite element simulation [START_REF] Gupta | Investigation of interaction between guidewire and native vessel using finite element analysis[END_REF]. Although the geometry has been defined based on patient preoperative data, no biomechanical model has been tuned according to patientspecific data.

Several studies have also been conducted in order to propose virtual reality-based training systems for endovascular surgery. Solutions simulating interactions between the vascular structure and flexible [START_REF] Tang | A realistic elastic rod model for real-time simulation of minimally invasive vascular interventions[END_REF], [START_REF] Luboz | Endovascular Guidewire Flexibility Simulation[END_REF], [START_REF] Luboz | Real-time Seldinger technique simulation in complex vascular models[END_REF] or stiff guidewires [START_REF] Lenoir | Interactive physically-based simulation of catheter and guidewire[END_REF] have been proposed. A " shortest path" approach has been suggested to estimate guidewire position within a rigid vascular structure [START_REF] Schafer | Real-time endovascular guidewire position simulation using shortest path algorithms[END_REF]. These studies, however, put forth no solutions for the assessment of vascular structure deformations due to real-time constraints.

In the following study, we present first the FEM-based matching methodology, including the workflow of the proposed approach, endovascular device features, preoperative patient features, the tool/tissue interaction model, and feature-based registration. We then present the results obtained on a dataset of 18 patients.

II. METHODS

A. Proposed approach

Fig. 1 illustrates the workflow for the proposed approach. Matching preoperative and intraoperative data involves 3D/2D rigid registration, as well as the mechanical simulation of the interactions between anatomical structures and endovascular material. This includes, specifically, the characterization and description of patient data from preoperative CT imaging and intraoperative angiographic/fluoroscopic imaging, the construction and tuning of the biomechanical model, the simulation of tool/tissue interactions via implicit finite-element analysis, and the projection of deformed model on intraoperative imaging.

B. Endovascular device features

The tool used in this analysis was a Lunderquist® (Extra Stiff Wire Guide, Cook®) stiff guidewire, which surgeon use to facilitate the insertion of aortic stent grafts during EVAR procedures. Among endovascular tools (i.e., flexible guidewires and deployment devices), it is known for most seriously deforming vascular structure. A circular beam with a diameter of 0.889mm represented the guidewire in the model. A linear elastic material model defined its mechanical behavior and its properties were set at standard values for an extra-stiff guidewire (Young's modulus of 200GPa and Poisson's coefficient of 0.3) [START_REF] Gupta | Investigation of interaction between guidewire and native vessel using finite element analysis[END_REF], [START_REF] Burns | Optimization of Medical Guidewires Through Controllability Matrix[END_REF].

C. Preoperative patient features

CT-scan data was analyzed with the sizing software Endosize® (Therenva, Rennes, France), used in routine practice, and which also allowed us to extract from the CTA the arterial lumen region and its centerlines, as well as vessel wall density (HU) reported to the edge of the vessel lumen (Fig. 2).

For each patient, a co mplete sizing was performed according to the standards of the International Society for Vascular Surgery (ISVS) [START_REF] Chaikof | Identifying and grading factors that modify the outcome of endovascular aortic aneurysm repair[END_REF]. For each patient, the tortuosity index (ratio of centerline to straight line distance) of the iliac arteries was determined. Each element of the mesh describing the arterial lumen was labeled with the corresponding density value obtained from the CTA. The mean of these values was calculated per arterial segment (aortic neck, aneurysm, and common iliac, external iliac, and femoral arteries).

In order to obtain a geometrical representation of the vascular structure, vessels of interest were described using Bspline curves in evenly-spaced planes orthogonal to the vessel centerline. The coordinates of the planes and spline curves were imported into the Ansys DesignModeler geometric modeling software. The vascular wall was reconstituted by interconnecting the spline curves with a surface interpolation tool. A Boolean operator was used to make the junctions between vessels. The vascular structure was meshed using three-noded triangular shell elements. Vascular wall thickness was set at 1.5mm for the aorta [START_REF] Raghavan | Wall stress distribution on three-dimensionally reconstructed models of human abdominal aortic aneurysm[END_REF] and 1mm for the iliac arteries.

The materials properties of the vascular structure can be modeled using behavior laws of varying precision. According to the literature and preliminary test simulations, a linear elastic material model appeared appropriate to describe vascular structure behavior. Indeed, a more precise law did not appear necessary, since the study was not primarily aimed at assessing wall stress distribution, and preliminary simulations tended to show that the model behavior was more sensitive to boundary condition changes than to material properties.

By analyzing tissue density values, it is possible to extract several tissue categories (calcified, healthy, thrombus, etc.) for which specific behavior laws may be used. The mesh elements were classified into two tissue categories. Healthy tissue was composed of elements with a d ensity of less than 240HU, whereas calcified tissue was composed of elements with a density of more than 240HU. Elastic modulus values reported in the literature vary widely in healthy and calcified arteries. Li et al. used a Young's modulus of 1.2MPa for normal arteries and 4.66MPa for aneurysm wall in order to model the fluid/structure interaction in a s tented abdominal aortic aneurysm [START_REF] Li | Blood flow and structure interactions in a stented abdominal aortic aneurysm model[END_REF]. Loree et al. applied circumferential tensile stress to human intimal atherosclerotic plaque and measured a tangential modulus of 1.466 ± 1.284MPa [START_REF] Loree | Static circumferential tangential modulus of human atherosclerotic tissue[END_REF]. Holzapfel et al. performed cyclic quasistatic uniaxial tension tests in the axial and circumferential directions and showed that calcifications have a linear mechanical response with an average Young's modulus of 12.6 ± 4.7MPa [START_REF] Holzapfel | Anisotropic mechanical properties of tissue components in human atherosclerotic plaques[END_REF]. Maier et al. performed uniaxial tension tests to determine the mechanical parameters of calcifications and showed that highly calcified AAA tissue exhibits linear elastic behavior with a Young's modulus varying from 40 to 450MPa [START_REF] Maier | Impact of calcifications on patient-specific wall stress analysis of abdominal aortic aneurysms[END_REF].

In our study, a linear elastic material was assigned to each of the two tissue categories (healthy and calcified) according to the values reported in the literature. Young's modulus was set at a value of 2MPa for elements of the "healthy tissue" category and 10MPa for elements of the "calcified tissue" category.

D. Tool-tissue interaction model D.1. Boundary conditions

Our study focused on macroscopic modeling of the aortoiliac structure. The model was developed by defining the boundary conditions with regard to anatomical-mechanical considerations emerging from the literature and surgeons' knowledge. Parts of the vascular structure with limited mobility could thus be identified. The suprarenal aorta is in contact with the spine, and is enveloped by a fibrous layer that maintains it f irmly against the spine. Furthermore, clinical expertise allowed us to identify the suprarenal aorta, aortic bifurcation, and origin of the internal iliac arteries as structures with limited mobility and highly constrained by surrounding structures.

Based on these considerations and observations, the boundary conditions of the model were defined as follows (Fig. 3):

-The celiac part of the abdominal aorta and the guidewire insertion area in the femoral artery region were considered as fixed.

-Elastic supports were used to model the behavior resulting from the interaction between vascular structures and surrounding tissues. These supports are boundary conditions incorporated within the simulation software which specify spring stiffness per unit area acting exclusively in the direction normal to the face of elements. Elastic supports were used to model the areas of anatomical constraint between the posterior face of the aorta and the anterior face of the vertebral column. Additional elastic supports were considered on the common iliac arteries with stiffness values to be tuned.

D.2. Pre-stress

The aorto-iliac structure is in constant interaction with blood flow and neighboring anatomical structures (muscles, fat, etc.). The configuration of the vascular structure visible on the preoperative CT-scan thus does not reflect its resting state (isolated structure). These neighboring structures are not explicitly represented, but their leverage must be taken into account in the simulation.

Clinical observations have suggested that arteries free of constraint were shorter and straighter than the same arteries in vivo. It was therefore decided to modify the resting state of the iliac arteries by taking these observations into account. Based on Holzapfel's experiments [START_REF] Holzapfel | Biomechanics of soft tissues with application to arterial walls[END_REF], S. Scherer et al. applied an axial pre-stress of 10% to arteries [START_REF] Scherer | Surgical stent planning: simulation parameter study for models based on DICOM standards[END_REF]. In our study, the same standard pre-stress value P was used. It was adjusted depending on the degree of iliac artery tortuosity T.

The shape of the artery at rest was determined based on a modified centerline intermediate between the initial centerline and a straight line between the aortic and iliac bifurcations. This centerline coordinates Cp were calculated using the following relationship: Cp = r * Cs + (1-r) * Ci where Cs corresponds to the straight line coordinates, Ci to the initial centerline coordinates, and r to a coefficient used to obtain the length of the desired artery.

The initial mesh of the iliac artery was then transformed around the new centerline in order to obtain the estimated shape of the artery at rest (Fig. 4). At the start of simulation, the calculated displacements were applied to the iliac artery. Parietal stresses were then cancelled, and the artery was placed back in its initial position, resulting in a pre-stress state.

D.3. Simulation

The simulation was designed to study the steady state of the model at the end of stiff guidewire insertion. The intermediate states corresponding to the progressive insertion of the guidewire were not studied.

At the beginning of the simulation, pre-stress was applied to the iliac artery in which the guidewire was to be inserted. The guidewire was then positioned on a line that constrains it within the vascular structure, minimizing its bending energy. This line was estimated by finding the path with minimal bending energy, using Dijkstra's algorithm and Hooke's law [START_REF] Schafer | Real-time endovascular guidewire position simulation using shortest path algorithms[END_REF]. The guidewire was then placed on this line, and its proximal extremity was fixed at the insertion site.

Tool-tissue interactions were modeled by contact without friction based on a penalty method. Contact elements were defined between the stiff guidewire and the internal surface of the vascular structure. The contact between the stiff guidewire and the aorto-iliac structure was then activated, and the imposed displacements were progressively cancelled until the steady state of the model was achieved (Fig. 5).

E. Feature based rigid registration

In order to match pre-and intra-operative data, the geometrical 3D/2D transformation, between 3D pre-operative coordinate system (associated to CTA and deformation simulation) and 2D intra-operative coordinate system (associated with angiography/fluoroscopy), was estimated through a f eature based registration process (Fig. 6). 3D-2D registration has been largely reported in the literature [START_REF] Markelj | A review of 3D/2D registration methods for image-guided interventions[END_REF]. We just present the outline of the implemented method which consists in finding six extrinsic rigid-body parameters [START_REF] Goksu | Endovascular navigation based on real/virtual environments cooperation for computer-assisted TEAM procedures[END_REF], the four intrinsic perspective projection parameters being known (from intra-imaging device data or calibration).

The image distortion produced by X-ray image amplifier was firstly corrected [START_REF] Maier | Impact of calcifications on patient-specific wall stress analysis of abdominal aortic aneurysms[END_REF]. The 2D data was derived from standard angiography, with a contrast medium injected so as to visualize the vascular structure. Successive subtracted angiography images were combined to visualize the entire vascular structure of interest. Thresholding was performed to obtain a binary image. The centerlines of the projected arteries were subsequently extracted by morphological skeletonization.

The registration distance criteria used to determine the 3D-2D transformation was the Euclidean 2D distance between the transformed (projected) 3D pre-operative vessel centerlines and the 2D vessel centerlines actually observed in the 2D intra-operative image. A Chamfer operator was used to precompute a d istance map from the 2D vessel centerlines. A Powell optimization algorithm was used to iteratively minimize the distance criteria, and thus to estimate the 3D-2D transformation.

The registration error was defined by the mean distance between the transformed 3D centerlines and the 2D intraoperative centerlines. In the following the simulation error was defined in the same way by considering 3D simulated guidewire and 2D observed guidewire.

F. Adaptive tuning

The parameters of the boundary conditions (stiffness of elastic supports) and pre-stress level of the iliac artery were adjusted, taking into account the projection of the simulated guidewire on the intraoperative images. After identifying the error between the projected simulated guidewire and the real guidewire, appropriate model parameters were interactively adjusted to reduce the distance between the two (Fig. 7).

Parameter tuning was performed on a training dataset. It was used to exhibit an analytical law with respect to geometrical and density characteristics of the artery wall. The behavior laws between parameter values and patient data (calcification degree, tortuosity) thus resulted in a unique adaptive model.

G. Experimental data analysis

Preoperative CT-scan data and intraoperative fluoroscopic images were acquired for 18 patients with an AAA treated by EVAR. This data was matched for each patient using rigid registration based on the centerlines of the vascular structures before stiff guidewire insertion.

The study was divided into two distinct phases. In an initial training phase conducted on a first patient group, the biomechanical properties of the model (material properties, boundary conditions, contact properties, etc.) were adjusted based on comparison with intraoperative images. Based on the results obtained, laws establishing the relationship between patient-specific "imaging" data and "biomechanical" data were assessed using polynomial regression so as to obtain a unique patient-adaptive model. In a second phase, the model obtained was investigated with a second patient group.

III. RESULTS

A. Data sets

CT-scan data were obtained for 18 consecutive patients who underwent EVAR in the Department of vascular surgery of the University Hospital of Rennes, France. Interventions were performed when aneurysms exceeded 50mm in diameter. The mean age was 70.9±10.1 years. This study was approved by the hospital's ethics committee.

CT-angiograms were carried out using a 6 4-sclice CT scan (General Electric Medical Systems, Milwaukee, Wisconsin, Fig. 7. Projection of the simulated guidewire on an intra-operative image. Observed guidewire appears in red and simulated guidewire in yellow. On the left, simulation error is not acceptable. On the right, simulation error is considered as acceptable after manual tuning. Fig. 6. Feature based 3D/2D registration. The 3D-2D transformation is estimated using the preoperative (red) and intraoperative (blue) centerlines and then used to evaluate the deviation between simulated (red) and observed (blue) guidewires.

USA, LightSpeed16). Data acquisition parameters were as follows: slice thickness of 1.25mm; tube of 215-260mA and 120kVp. The patients were injected with 120ml of a non-ionic iodinated contrast medium (Hexabrix, iodine concentration of 320mg/mL). CT windowing was performed at 400 HU in the periphery and 40HU in the center.

The first group A, composed of 10 patients, was considered as the training group. The second group B, composed of the remaining eight patients, was considered as the test group.

B. Rigid registration

The registration error of each patient is presented in Fig. 8. The mean registration error for the entire patient population was 1.55±0.45mm, with no significant difference between the training and the test groups (p = 0.756). The registration error for all 18 patients was correlated with the angulation of the aortic neck (p= 0.014), the tortuosity index of the aneurysm (p=0.003), and the aneurysm angle (p=0.010), which seems to be in line with the results obtained by Carrel et al. [START_REF] Carrell | Feasibility and Limitations of an Automated 2D-3D Rigid Image Registration System for Complex Endovascular Aortic Procedures[END_REF].

C. Adaptive setting of simulation parameters

Simulations were carried out using the Ansys Mechanical finite-element solver on a Hewlett Packard Z800 workstation (HP Development Company, California, USA) equipped with a 6-core Xeon processor (3.46GHz).

A first manual tuning of the parameters was performed on the training group (A). Model parameters were adjusted independently for each patient, taking into account simulation results projection on the intraoperative images. Simulation errors (including registration error) for the training group are presented in Table 1. The mean simulation-related error was 2.1±0.8mm.

Based on these parameters, two behavior laws were established between parameter values and patient data. Prestress level was correlated (R² = 0.86) with tortuosity of the iliac artery (Fig. 9): P = -0.6 * T² + 0.9652 * T + 0.6377.

In the same way, stiffness of the elastic support S was correlated (R² = 0.79) with mean density D of the common iliac arterial wall (Fig. 10): S = 0.00002 * D -0.000028 N/mm 3 .

Model parameters for each patient were then modified during a second phase in order to check the previously established laws. For this second phase, simulation errors (including registration error) are presented in Table 2. Mean simulation-related error (including registration error) was 2.3±0.6mm. One patient's intra-operative data was used to match 3D and 2D data on a s econd incidence. Simulation results could thus be compared with intraoperative imaging not used for the tuning of the model parameters. The simulated guidewire was thus projected on two images with a d ifferent angle of incidence (left anterior oblique [LAO] 20° + caudal 13-15° / LAO 30° + cau dal 13-15°). The simulation-related error for this patient was 3.5±2.5mm for the first incidence and 2.0±1.3mm for the second incidence.

D. Evaluation of pre-/ intra-operative matching

The FEM-based matching was applied to the test group (B) using the adaptive model. Simulation errors (including registration error) for this second group are presented in Table 3. Mean simulation-related error was 2.8±0.5mm. Mean simulation calculation time was 301±168 seconds. Typically, 10 minutes were required for data analysis and extraction, 10 minutes for preparation of the simulation, 5 minutes for simulation, and 2 minutes for registration. Simulation results can thus be used to integrate the deformation of the aorto-iliac structure in 3D/2D fusion (Fig. 11 and12). The matching between the vascular and the observed guidewire is largely improved when the deformed model is used.

IV. DISCUSSION

With the extension of EVAR to more and more complex cases (difficult anatomies, and branched and fenestrated stent grafts), a navigation system able to integrate the deformation of vascular structure would be of main interest. The use of such a system could ultimately impact intervention duration, the patient x-ray dose, contrast injection, and stent graft placement.

In this study, we presented a new method for the matching of 3D preoperative data and 2D intraoperative images taking into account anatomical deformations caused by endovascular devices during EVAR procedure.

Vessel deformations were computed during the preoperative phase using a patient-specific finite element model. Patientspecific simulation makes it difficult to take into account the biomechanical properties of each surrounding tissue. In order to focus on evaluating the macroscopic behavior of the vascular structure, we modeled the surrounding tissues with boundary conditions resulting from surgical expertise. This model did not integrate such details as layer-specific material properties or the non-uniformity of wall thickness. The impact of boundary conditions on the behavior of the model was subject to specific study. This work concentrated on the image-based tuning of boundary conditions, though very few data are reported in the literature.

When considering undeformed vascular structures, registration based on pre-operative and intra-operative centerlines appears to be accurate. Registration error, which averages 1.5mm, allows us to correctly assess the finite element simulation results.

In the context of deformation simulations, the interactive tuning of model parameters yielded acceptable results, with simulation errors close to 2mm, and highlighted a unique adaptive model. Indeed, the first tuning led to behavior laws with a good correlation. The mean simulation error was of the same order of magnitude for both the manual and adaptive tuning models, which suggests the correct parameters were tuned. Despite the straightforward relationships between assumed [START_REF] Penney | An Image-Guided Surgery System to Aid Endovascular Treatment of Complex Aortic Aneurysms: Description and Initial Clinical Experience[END_REF]. Simulation calculation time was lower than 10 minutes for each patient. We consider this fast enough to envisage, after further developments, a role for finite element simulation in the clinical workflow (provided it is run in preoperative phase). In order to make this possible, simulations must be sufficiently predictable; this appears to be the case for the 18 patients.

Patient-specific simulation is a difficult question when it comes to the detailed representation of the reality of a patient's non-linear material properties, boundary conditions, and prestress. We deliberately adopted a s implifying approach designed to implement the whole processing chain and address the issue of matching pre-and intra-operative data. Some inevitable limitations can be observed, especially when the different steps of the process are considered individually.

First, the construction of the geometrical model induces some smoothing which may influence simulation results. A possible improvement would be to work directly on a mesh extracted from the segmentation of the vascular structure.

Moreover, we assumed that the aorta was uniform in wall thickness. This assumption can be criticized because it is known that, at least in aneurysms, wall thickness varies [START_REF] Raghavan | Regional distribution of wall thickness and failure properties of human abdominal aortic aneurysm[END_REF]. However, as wall thickness is not easily identifiable from CT scans, most studies assume a uniform thickness, even for the analysis of wall stress.

Linear material properties were chosen to describe vascular structure behavior. This could be criticized, considering that nonlinear models are often reported in the literature, especially in wall stress analysis. Nevertheless, we assumed linear law was sufficient for a first approach given the global behavior of the vascular structure. Moreover, some studies showed that arterial tissue behavior tended to be linear at a blood pressure of 80mmHg [START_REF] Doyle | A comparison of modelling techniques for computing wall stress in abdominal aortic aneurysms[END_REF], [START_REF] Loree | Static circumferential tangential modulus of human atherosclerotic tissue[END_REF], [START_REF] Thubrikar | Wall Stress Studies of Abdominal Aortic Aneurysm in a Clinical Model[END_REF]. We nevertheless plan to use a nonlinear law, as described by Raghavan et al. [START_REF] Raghavan | Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: identification of a finite strain constitutive model and evaluation of its applicability[END_REF], in a future work in order to evaluate the impact of this assumption on the simulation accuracy.

Another point for improvement is the evaluation of iliac geometry at rest; in our study, the geometrical method used was straightforward. Even if our method yielded good results, we intend to improve this point by using an inverse design analysis so as to calculate a stress-free reference configuration [START_REF] Gee | A computational strategy for prestressing patient-specific biomechanical problems under finite deformation[END_REF], [START_REF] Putter | Patient-specific initial wall stress in abdominal aortic aneurysms with a backward incremental method[END_REF].

Although the time required for the process seems to fit the requirement of clinical workflow, a cer tain amount of user interaction is needed at different steps, especially during preparation of the simulation. We plan to address this issue and hope to automatize certain steps (contact zones, boundary conditions, mesh parameters, and simulation options).

We took into account the effect of spine on the vascular structure through elastic supports, but we neglected other surrounding tissues, including muscle, fat, and thrombus. This simplifying hypothesis could contribute to simulation error. However, full consideration of the surrounding tissue would require the fine segmentation of each tissue; this would be very time-consuming and difficult to do in clinical practice. In order to remain consistent with the clinical workflow and improve simulation accuracy, we would have to find a good balance between a detailed biomechanical simulation and a simplified simulation jointly operating with an intraoperative imagebased correction step.

V. CONCLUSION

We presented a n ew method for the matching of 3D preoperative data (images, models) and intraoperative 2D observations affected by anatomical deformations during EVAR procedure. Patient-specific simulation of tool/tissue interactions is especially complex because the intraoperative step-for which observations showing the endovascular device interacting with tissues are available-is generally not exploited. The proposed solution involves a finite element model whose mechanical parameters (boundary conditions and pre-stress level) were adjusted based on comparisons with intra-operative images. The image-based adaptive tuning seems an effective tool that yields acceptable simulation results, especially when boundary conditions are not well known. These first results show that the proposed adaptive model is patient-specific at a macroscopic level, and can thus properly estimate the deformation of the vascular structure. This approach constitutes the basis of an augmented-reality endovascular navigation system. With further development, the proposed approach could be used to restore to practitioners an estimation of tissue deformation at different steps of the EVAR procedure.

Fig. 1 .

 1 Fig. 1. Workflow of the finite element-based matching of pre-and intraoperative data.

Fig. 2 .

 2 Fig. 2. Pre-operative sizing. (A) corresponds to the measurement point on the centerline immediately below the renal arteries, (B) to the end of the aneurysm neck, (C) to the aortic bifurcation, (D) and (E) to the left and right iliac bifurcations respectively. Calcification distribution appears in white on the vascular structure.

Fig. 5 .

 5 Fig. 5. Simulation of tool/tissue interactions. The pre-stress is applied to the iliac artery (left). The guidewire is positioned within the vascular structure (middle). The guidewire is relaxed and deforms the aorto-iliac structure until steady state of the model is reached (right).

Fig. 4 .

 4 Fig. 4. Geometrical modification of the iliac artery mesh for pre-stressing. The straight line Cs is represented by the dashed line, the initial centerline Ci in light gray, the new centerline Cp in black, the initial mesh in light orange, and the new mesh in dark orange.

Fig. 3 .

 3 Fig. 3. Boundary conditions of the finite element model. The proximal (A) and distal (D) extremities of the model are fixed. Elastic supports (in blue) are added to the suprarenal aorta (B) and the common iliac arteries (C).

Fig. 10 .

 10 Fig. 10. Elastic support stiffness versus mean density of the common iliac arterial wall (for group A).

Fig. 9 . 8 Fig. 8 .

 988 Fig. 9. Pre-stress level versus iliac artery tortuosity (for group A).

Fig. 11 .

 11 Fig. 11. Roadmap with preoperative aorto-iliac structure. The inserted stiff guidewire has deformed the vascular structure. It is not superimposed with the preoperative aorta model.

Fig. 12 .

 12 Fig. 12. Roadmap with simulated aorto-iliac structure deformed by the stiff guidewire. The observed guidewire is superimposed with the deformed aorta model.

TABLE 3 SIMULATION

 3 ERROR FOR THE ADAPTIVE MODEL (GROUP B)patient data and model parameters, results from the test group were almost as good as those from the training group. Simulation error averaged less than 3mm, which is within the range of acceptable clinical accuracy that Penney et al.

	Patient	Mean	Min -Max
	11	2.6 ± 1.6	0.0 -6.6
	12	3.6 ± 2.5	0.0 -9.1
	13	4.3 ± 1.9	0.0 -7.8
	14	1.7 ± 0.8	0.0 -3.0
	15	1.2 ± 0.9	0.0 -4.2
	16	3.1 ± 1.9	0.0 -6.2
	17	3.6 ± 1.6	0.0 -5.7
	18	2.1 ± 1.5	0.0 -4.8

ACKNOWLEDGEMENTS

This work has been partially conducted in the experimental platform TherA-Image (Rennes, France) supported by Europe FEDER.

This work was supported in part by the French national research agency (ANR) through the TecSan program (project ANGIOVISION n°ANR-09-TECS-003).


Romain Couillet
email: romain.couillet@supelec.fr

Samir M Perlaza
email: samir.medina-perlaza@supelec.fr

Hamidou Tembine
email: hami-dou.tembine@supelec.fr

Mérouane Debbah
email: merouane.debbah@supelec.fr

Electrical Vehicles in the Smart Grid: A Mean Field Game Analysis

In this article, we investigate the competitive interaction between electrical vehicles or hybrid oilelectricity vehicles in a Cournot market consisting of electricity transactions to or from an underlying electricity distribution network. We provide a mean field game formulation for this competition, and introduce the set of fundamental differential equations ruling the behavior of the vehicles at the system equilibrium, namely the mean field equilibrium. This framework allows for a consistent analysis of the evolution of the sale-and-purchase price of electricity as well as of the instantaneous total demand. Simulations precisely quantify those parameters and suggest that following the charge and discharge policy at the equilibrium allows for a significant reduction of the daily electricity peak demand. I. INTRODUCTION Electrical vehicles (EV) and plug-in hybrid electrical vehicles (PHEV) have been recognized as natural components of future electricity distribution networks, known as smart grids [1], [2], [3]. As opposed to classical vehicles, EV and PHEV are equipped with batteries which can be charged or discharged by using a simple plug-in connector compatible with the local electricity distribution grid. Thus, EV and PHEV can be conceived as both energy consuming devices and mobile energy sources [4], [5], [6], [7].

In the former case, EV and PHEV can be seen as devices straining the energy demand of energy suppliers

and, thus, adding a new constraint to reliably distribute the electricity. In the latter case, EV and PHEV can be used to store or even to transport the energy from one geographical area to another and then to increase the reliability of the energy supply in certain zones or time intervals. That is, either areas where the energy supply relies only non-easy predictable sources, e.g. wind turbines, solar panels, etc., or at time intervals where the demand naturally suddenly increases (rush hours).

Within this framework, it is therefore an important economical and social challenge to enforce charge and discharge policies to EV and PHEV in an optimal manner. Here, optimality must be interpreted in the sense of individual revenue obtained by the EV and PHEV owners when participating in the energy trades and also in terms of reliability of the energy supply process to the fixed consumers. In this paper, we consider that a way to improve reliability is to allow EV and PHEV to buy and sell energy to or from the smart grid, as in a classical Cournot competition [START_REF] Cournot | Researches into the mathematical principles of the theory of wealth[END_REF]. Clearly, the price at which the energy is sold and bought depends on the existing demand in the grid and also on the demand and offer resulting from all the vehicles connected to the network. This competitive interaction resulting from the energy trade where each vehicle owner decides the amount of energy to be sold or bought, given a global price, can be analyzed using tools from dynamic game theory [START_REF] Bas | Dynamic Noncooperative Game Theory[END_REF]. For instance, in [START_REF] Saad | A noncooperative game for double auction-based energy trading between PHEVs and distribution grids[END_REF], the coexistence of a number of PHEV groups aiming to sell part of their stored energy to the smart grid is studied using non-cooperative game theory [START_REF] Fudenberg | Game theory[END_REF]. An algorithm based on best response dynamics [START_REF] Fudenberg | Game theory[END_REF] is proposed to allow PHEV groups to reach a Nash equilibrium [START_REF] Nash | Equilibrium points in n-person games[END_REF].

Nonetheless, in practical scenarios, the number of vehicles might be drastically large and, thus, elements from classical game theory might not necessarily bring enough insight about the global behavior of the market. To overcome this problem, in this paper, we study the energy trade when the number of vehicles tends to infinity and all vehicles are considered alike, following the paradigm of [START_REF] Aumann | Markets with a continuum of trackers[END_REF]. See also [START_REF] Carmona | Nash equilibria of games with a continuum of players[END_REF] for recent results. More precisely, we shall model this interaction as a mean field game [START_REF] Lasry | Jeux à champ moyen. i -le cas stationnaire[END_REF], [START_REF]Jeux à champ moyen. ii -horizon fini et contrôle optimal[END_REF]. In contrast to games with a finite number of players, where each player follows the evolution of the state of the game and the actions taken by all the other players in order to maximize a given individual benefit, in the mean field game formulation, every player's action is driven not by the individual actions of each other player but by the collective (or mean) behavior of all players. Typically, the notion of (Nash) equilibrium in the context of mean field games is known as mean field equilibrium (MFE) [START_REF] Lasry | Jeux à champ moyen. i -le cas stationnaire[END_REF], [START_REF]Jeux à champ moyen. ii -horizon fini et contrôle optimal[END_REF]. The MFE is found as the solution of a coupled system of stochastic partial differential equations (SPDE) which includes a (backward) Hamilton-Jacobi-Bellman (HJB) equation and a (forward) Fokker-Planck-Kolmogorov (FPK) equation.

The closest contribution to our specific problem setting is [START_REF] Guéant | Mean field games and oil production[END_REF]. Therein, a mean field game approach to the study of oil production in infinite-time horizon is developed. In [START_REF] Guéant | Mean field games and oil production[END_REF], the selfish players are oil producers and the mean field variable is the oil selling price. In this article, we develop a similar framework to [START_REF] Guéant | Mean field games and oil production[END_REF] but on a finite time horizon, applied to both EV and PHEV, with vehicle owners as the selfish players and electricity price as the mean field variable. In general, it is difficult to prove that the mean field solutions are well-defined limits of the finite game, see e.g. the discussions in [START_REF]Mean field games and applications[END_REF]. Similarly, it is in general difficult to prove the existence and uniqueness of solutions of one equilibrium and, if so, to derive a numerical method that is provably able to converge to an equilibrium, see e.g. [START_REF] Achdou | Mean field games: Numerical methods[END_REF]. We will not try to prove any of these aspects here, our main target being rather to informally explore the potentials of the mean field game setting to the economical problem of optimal policies for EV and PHEV penetrations in the smart grid. Similar to [START_REF] Guéant | Mean field games and oil production[END_REF], only the numerical results will convey a justification of the correct behavior of our method.

The reminder of this article unfolds as follows. In Section II, we describe the problem formulation in the case where only electrical vehicles interact with the energy market. Therein, the problem is formulated as a continuous time differential game with finite horizon. This formulation is then written under the form of a mean field game and the differential equations describing the mean field equilibrium are presented.

In Section III, the same analysis presented for EV is carried out for the case of PHEV. In Section IV, we provide numerical simulations and derive conclusions for both scenarios. Finally, in Section V, we conclude this work.

II. PURELY ELECTRICAL VEHICLES

A. System Model

Consider a finite set K = {1, . . . , K} of EV connected to an electricity distribution network. The average consumption rate of vehicle k ∈ K at time t ∈ [0, T ] is denoted by g

(k)
t . This consumption rate is measured in units of electricity per time. We assume the consumption function g the energy provisioning rate of vehicle k at time t, that is, the rate at which vehicle k buys or sells its energy. We relate the variable x

(k) = {g (k) t , 0 ≤ t ≤ T } of EV k to
d dt x (k) t = α (k) t -g (k) t . (1)
In the following, we denote

x t = (x (1)
t , . . . , x battery level and provisioning rates of all EV at time t, respectively. Consider now a predefined period

[0, T ]. We denote x (k) = {x (k) t , 0 ≤ t ≤ T } and α (k) = {α (k)
t , 0 ≤ t ≤ T } the trajectories of the battery level and provisioning rates for EV k, respectively. We also denote x = {x t , 0 ≤ t ≤ T } and α = {α t , 0 ≤ t ≤ T } the joint trajectories of the battery levels and provisioning rates. We finally denote X the set of all possible state trajectories and A k the set of all possible controls α (k) of player k.

The price at which vehicles either sell or buy electricity at time t is evaluated by the function p t : R K → R, α t → p t (α t ). The time dependency of the price p t models a realistic dynamic pricing policy accounting for the energy demand for other needs than EV battery loading. This function can be tuned to create incentives for EV to sell or buy energy at specific time periods. In addition to electricity price, other factors influence the selling and buying behavior of EV owners. We model these, for player k, by the following set of functions. The function h

(k) t : R → R, α → h (k)
t (α) models the (psychological) cost for player k to buy or sell electricity at rate α at time t. Indeed, EV owners are more likely to trade energy at some convenient time intervals, e.g. during nighttime when the EV is plugged to the distribution network. The function f

(k) t : [0, 1] → R, x → f (k)
t (x) models the cost for vehicle k to possess only a fraction x of energy reserves at time t. For instance, during periods of high energy consumption, the interest of EV owners is to have as fully loaded batteries as possible. Finally,

κ (k) [0, 1] → R, x → κ (k) (x)
models the cost for EV k to end the trade period [0, T ] with a fraction x of battery load. This function guarantees that EV owners do not sell all their battery content at the end of the trade. For simplicity, we will assume that the aforementioned functions are sufficiently regular for the following derivations to make sense. Specific realistic choices of functions are considered in Section IV.

The goal of EV k is to determine the consumption rates α (k) that minimizes its total cost over a time window [0, T ], given by

J k x 0 , α (k) , α (-k) = T 0 α (k) t p t (α t ) + h (k) t (α (k) t ) + f (k) t (x (k) t ) dt + κ (k) (x (k) T ) (2)
for a given initial battery level vector x 0 at time t = 0.

B. Classical Game Formulation

We model the energy trades resulting from the interactions among the electrical vehicles and the smart grid by a K-player continuous-time differential game of pre-specified fixed duration T > 0. Let K, the set of EV, be the set of players. Player k controls the its state trajectory x (k) through the control variable

α (k)
t . The state trajectory x of the game is determined by the initial state x 0 and by the transport equation October 1, 2011 DRAFT [START_REF] Simpson | Cost-benefit analysis of plug-in hybrid electric vehicle technology[END_REF]. The cost function of player k is defined by [START_REF] Bradley | Design, demonstrations and sustainability impact assessments for plug-in hybrid electric vehicles[END_REF] and we assume that the information available for player k at time t is the actual price p t (α t ).

The objective for player k is to choose the control functions α (k) ∈ A k that minimize the cost J k given the initial condition x 0 ∈ R K and the instantaneous price p t (α t ). We denote α (-k) ∈ A 1 × . . . × A k-1 × A k+1 × . . . × A K the trajectories of actions of all players but k. Here, the interdependence between players appears through the sale and purchase electricity price; the control α (k) t depends on the price p t (α t ), which depends itself on all the other players' controls α

(-k) t .
The formulation of the game is completed by restricting the deterministic function g t and the set A 1 × . . . × A K to lie in a subset of admissible consumption rates and control profiles α t , in order for the solution of (1) to exist and be unique. This condition is necessary to ensure that the state trajectory

x is unique and, thus, to define unequivocally the cost functions [START_REF] Bradley | Design, demonstrations and sustainability impact assessments for plug-in hybrid electric vehicles[END_REF].

Finally, we interpret the notion of equilibrium of the game in the sense of Nash [START_REF] Nash | Equilibrium points in n-person games[END_REF], [START_REF] Aubin | Mathematical methods of game and economic theory[END_REF] defined as follows,

Definition 1 (Nash Equilibrium (NE) [START_REF] Nash | Equilibrium points in n-person games[END_REF]): The K-tuple of control functions α = α (1) , . . . , α

(K) ∈ A 1 × . . . × A K is a Nash equilibrium if for all k ∈ K and for any admissible control α (k) ∈ A k , it holds that J k x 0 , α (k) , α (-k) J k x 0 , α (k) , α (-k) , (3)
for a given initial state x 0 .

Our interest in the NE lies in the fact that, at a state of NE, all the EV are using a control policy which is optimal with respect to the control policy of all other EV. Otherwise stated, if the system is at equilibrium, the unilateral deviation in control of any player would lead to a higher cost and, thus, none of the players has an interest in unilaterally changing its control function. Nonetheless, analyzing the NE of such a game, where K is a large number, is a difficult problem. In fact, even if a NE exits, it would lead to solutions that are inherently difficult to exploit. In particular, it is clear that, under this formulation, any change in the battery level of a given player impacts all other players which must react as a consequence. We aim at reducing this complexity by adopting some additional, but reasonable, conditions.

C. Mean Field Game Formulation

The game described in the previous section can be remarkably simplified by considering a large population of symmetric players. We assume that the number of players K grows to infinity in a fluid manner, in such a way that the set of players will now be considered continuous and that the individual cost functions become identical. This implies that the set of controls is now the same for all players. We will denote A k = A. In this context, the impact of a particular control adopted by a single individual represents a negligible cost variation on all the other players' costs. Here, players are sensible to the distribution of all other players' states. We model this distribution by a limiting counting measure

F : [0, T ]×B → [0, 1], F (t, B) = lim K→∞ 1 K K k=1 δ {x (k) t ∈B}
where B is the Borel σ-field on [0, 1]. We assume F exists unequivocally and has a density m(t, x)

= d dx F (t, x).
In order to relax the strong constraint g (k) t g t , now unrealistically identical for all users, the individual state of each player is assumed to be a noisy version of the mean state trajectory. That is, we define the evolution of the state x (k) t

x t by the following stochastic differential equation (SDE),

dx t = α t dt -g t (dt + σ t dW t ) + dN t . (4)
Here, W t is a Brownian motion and the term dW t must be interpreted in the sense of Itô [START_REF] Itô | Multiple Wiener integral[END_REF]. As such, the energy consumption becomes a random variable g t [dt + σ t dW t ], at time t. The term dN t is a reflective variable to ensure that x t remains in [0, 1]. Now, under the assumption of independence between the different dimensions of the Brownian motion W t , the analysis of the game reduces to the state trajectory analysis of a single player. The cost function of a given player can be obtained from (2) by dropping the player index and considering the density m(t, •), as follows,

J (x 0 , α, m 0 ) = E T 0 (α t p t (m) + h t (α t ) + f t (x t )) dt + κ(x T ), (5)
where m ∈ M denotes the trajectory of the probability distributions, i.e. m = {m(t, •), 0 ≤ t ≤ T }, and M is the set of all such trajectories. The initial state conditions for the considered player are

x 0 ∈ [0, 1]
and m 0 m(0, •), the initial density function. In this context, the energy sale and purchase price function

writes p t : M → R + , m → p t (m).
Indeed, the price is a function of the total instantaneous demand 1 0 α t (x)m(t, x)dx. However, for computational ease, we will instead consider that prices are fixed not by the total consumption 1 0 α t (x)m(t, x)dx but by the expected consumption g t + d dt ( 1 0 xm(t, x)dx), where both quantities only differ by an additional Brownian motion term when σ t > 0. In practice, this suggests that the energy providers which set the instantaneous prices do not have the information on the instantaneous demand at time t but are able to track the density trajectory m. We therefore define p t as

p t (m) = D(t, •) -1 g t + d dt 1 0 x m(t, x)dx ,
where D(t, p) is the total energy demand function at time t for a given price p, and the inverse is with respect to the operation of function composition. Under the above assumptions, the continuous time differential game discussed in Section II-B becomes a mean field game as introduced in [START_REF] Lasry | Jeux à champ moyen. i -le cas stationnaire[END_REF], [START_REF]Jeux à champ moyen. ii -horizon fini et contrôle optimal[END_REF].

D. Mean Field Equilibrium

Our interest now is to transpose the notion of NE into the corresponding notion of equilibrium in the mean field game, namely the mean field equilibrium (MFE) [START_REF] Lasry | Jeux à champ moyen. i -le cas stationnaire[END_REF], [START_REF]Jeux à champ moyen. ii -horizon fini et contrôle optimal[END_REF]. Based on Definition 1, a MFE is fully described a control α ∈ A solution to the stochastic control problem,

v(0, x 0 ) = inf α∈A [J (x 0 , α, m 0 )] dx t = α t dt -g t [dt + σ t dW t ] + dN t , (6)
with initial state x 0 and m 0 .

Let us define the partial cost function on the time interval [u, T ] as

v(u, x • ) = inf α∈A E T u (α t p t (m) + h t (α t ) + f t (x t )) dt + κ(x T )
where the eligible trajectories x t are such that

x u = x • . We will denote v = {v(t, •), 0 ≤ t ≤ T }
∂ t v(t, x) = -inf α∈R , {α∂ x v(t, x) + C(α, x, m (t, x), t)} + g t ∂ x v(t, x) - 1 2 g 2 t σ 2 t ∂ 2 xx v(t, x) (7)
∂ t m(t, x) = -∂ x [(α t -g t )m(t, x)] + 1 2 g 2 t σ 2 t ∂ 2 xx m(t, x). (8)
Equation ( 7) is a backward HJB equation which determines the optimal trajectory α in the control problem [START_REF] Sojoudi | Optimal charging of plug-in hybrid electric vehicles in smart grids[END_REF], while Equation ( 8) is a forward FPK equation which determines the optimal density measure m .

We assume here that the cost h(α t ) for control is quadratic and reads

h(α t ) = 1 2 h t α 2 t ,
with h t > 0 representing the unwillingness of the car owner to buy or sell energy at time t. This choice is seemingly non-natural as it implies that users are more willing to buy or sell small quantities rather than October 1, 2011 DRAFT large quantities of energy. Nonetheless, under the mean field game formulation, this has to be understood as the fact that, on average, only a limited population of users at time t is willing (or able) to buy energy.

As such, intuitively, making the (psychological) cost of buying or selling energy larger for larger amounts of energy forces only part of the population to buy or sell. As for the particular choice of a quadratic cost rather than any other cost function, it is convenient for calculus mostly.

Under this assumption, solving

inf α∈R , {α∂ x v(t, x) + C(α, x, m (t, x), t)}
for all t, it is immediate to see that the optimal trajectory α is explicitly given by

α t = - 1 h t [∂ x v(t, x) + p t (m )] , (9)
possibly subject to some boundary conditions to ensure that x t ∈ [0, 1] at all times. In the remainder of the article, we will assume this condition always met, so that at no time we will consider EV owners with completely full or completely empty batteries.

The HJB equation now becomes

0 = ∂ t v(t, x) - 1 h t [∂ x v(t, x) + p t (m )] + g t ∂ x v(t, x) - p t (m ) h t [∂ x v(t, x) + p t (m )] + f t (x) + 1 2h t [∂ x v(t, x) + p t (m )] 2 + 1 2 σ 2 t g 2 t ∂ 2 xx v(t, x),
which can be simplified as

∂ t v(t, x) = 1 2h t (∂ x v(t, x) + p t (m )) 2 + g t ∂ x v(t, x) -f (t, x) - 1 2 σ 2 t g 2 t ∂ 2 xx v(t, x)
and the FPK equation is

∂ t m(t, x) = 1 h t [∂ x v (t, x) + p t (m )] + g t ∂ x m(t, x) + 1 h t ∂ 2 xx v (t, x)m(t, x) + 1 2 g 2 t σ 2 t ∂ 2 xx m(t, x).
This defines the two fundamental differential equations to be solved, either explicitly or numerically, for determining the MFE.

In the next section, we improve the EV framework by turning the purely electrical vehicles into PHEV, introducing therefore the possibility for players to select between two alternative sources of energy.

III. PLUG-IN HYBRID VEHICLES

A. System Model

In this section, we consider that vehicles in the set K are PHEV. A PHEV can operate both with an electrical energy source and an alternative energy source, for instance oil. Thus, a PHEV interacts with the electricity distribution grid by buying and selling electricity with an elastic price. In the case of oil, it can be only bought at a fixed price, which is a natural assumption on a daily or even weekly basis.

We describe the energy reserves of PHEV k by the two-dimensional vector z

(k) t = (z (k) 1,t , z (k) 2,t ) T ∈ [0, 1] 2
, where z 1,t is the amount of energy stored in the batteries and z 2,t the level of the oil reserve. We denote the provision rates of electricity and oil of PHEV k by µ

(k) 1,t ∈ R and µ (k)
2,t ∈ R, respectively. In addition, we denote β (k) : R + ×[0, 1] 2 → [0, 1], (t, z) → β (k) (t, z) the function that determines the relative proportion of energy drawn from the batteries of PHEV k at time t. Typically, taking

β (k) (t, z) = z 1 /(z 1 + z 2 ),
with z = (z 1 , z 2 ) T , translates a policy where energy is consumed indistinctly of the energy source.

Note that, depending on the typical distances covered by PHEV owners at time t (e.g. weekdays against weekends), β (k) (t, z t ) may explicitly depend on t. Alternatively, we may have considered β (k) (t, z t ) an additional control variable which can be set optimally by the car owner depending on the status of the energy market. Nonetheless, for simplicity of analysis, we do not consider this scenario here. We relate the variables z t , µ t and β (k) by the following two-dimensional transport equation,

d dt z t =   µ 1,t µ 2,t   -   β(t, z t ) 1 -β(t, z t )   g t . (10)
We define the cost of PHEV k in the time window [0, T ] as,

L k z 0 , µ (k) , µ (-k) = T 0 r t µ (1)
t , . . . , µ

(K) t + q (k) t (µ (k) t ) + s (k) t (z (k) t ) dt + ξ (k) (z (k) T ), (11)
for a given initial state z 0 ∈ [0, 1] 2K , i.e. the energy reserves of all PHEV, where µ = (µ (1) , . . . , µ (K) ),

with µ (k) = {µ (k) t = (µ (k) 1,t , µ (k) 2,t ), 0 ≤ t ≤ T }. Here, r t : R 2K → R 2 , µ → r t (µ (1) , . . . , µ (K) ) = (r 1,t (µ (1) 1 , . . . , µ (K) 1 ), r 2,t (µ (1) 2 , . . . , µ (K)
2 )) evaluates the instantaneous prices r 1,t of electricity and r 2,t of oil, given the controls µ (k) = (µ

(k) 1 , µ (k)
2 ). In particular, we assume here that the price for oil is fixed, given by r 2,t (µ

(1) 2 , . . . , µ (K)
2 ) = r 2 . Note that in this case the trajectory of the state x = {x t = (x (1) t , . . . , x (N ) t ), 0 ≤ t ≤ T } is determined by the initial state x 0 = (x (1) 0 , . . . , x (K) 0 ) and the transport equation [START_REF] Saad | A noncooperative game for double auction-based energy trading between PHEVs and distribution grids[END_REF]. The function q (k) t : R 2 → R, µ → q (k) t (µ) evaluates the psychological cost of trading a quantity µ 1 of electricity and a quantity µ 2 of oil at time t, where µ = (µ 1 , µ 2 ) T . The function s

(k) t : [0, 1] 2 → R, October 1, 2011 DRAFT z → s (k)
t (z) denotes the cost for PHEV k to be in state z at time t. Finally, ξ (k) : [0, 1] 2 → R, z → ξ (k) (z) is the cost for PHEV k to be in state z at time T . These are the analogous to the functions

h (k) t , f (k)
t , and κ (k) in (2), respectively. In the following, we formulate the finite-number of players differential game.

B. Classical Game Formulation

The interaction between all PHEV is modeled by a K-player continuous-time stochastic differential game of pre-specified fixed duration T > 0. As for the case of EV, the set of players is the set K of PHEV and the objective for player k is to choose the control trajectory µ [START_REF] Fudenberg | Game theory[END_REF], given the initial conditions z 0 and the control functions µ (-k) adopted by all the other players. We denote the set of all possible controls µ (k) of player k over the time period [0, T ] by U k . The state trajectory z is determined by z 0 and by the evolution function [START_REF] Saad | A noncooperative game for double auction-based energy trading between PHEVs and distribution grids[END_REF]. The set of state trajectories z = (z (1) , . . . , z (K) ), with z (k) = {z (k) t , 0 ≤ t ≤ T }, is denoted by Z. As in the previous case, the information available for each player k at time t is the initial state z 0 and the actual instantaneous price p t (m). Finally, we impose the necessary conditions for [START_REF] Saad | A noncooperative game for double auction-based energy trading between PHEVs and distribution grids[END_REF] to have a unique solution and thus we ensure that the resulting game is well formulated.

(k) = {µ (k) t , 0 ≤ t ≤ T } that minimizes the cost L k in

C. Mean Field Game Formulation

In this section, we proceed similarly to Section II-C. We use here a finite-game counting measure

F (K) : R + × B × B → [0, 1] of the form F (K) (t, B 1 , B 2 ) = 1 K K k=1 δ {z (k) 1,t ∈B1, z (k) 2,t ∈B2} , (12)
and we assume that it admits a limiting measure F with density m(t, •, •), i.e.

lim K→∞ F (K) (t, B 1 , B 2 ) = F (t, B 1 , B 2 ) (13) = B1 B2 m(t, z 1 , z 2 )dz 1 dz 2 . ( 14
)
The function m(t, z 1 , z 2 ) evaluates the density of players lying in state z = (z 1 , z 2 ) T . As previously, the individual state of each player is assumed to be a noisy version of the deterministic state trajectory in [START_REF] Saad | A noncooperative game for double auction-based energy trading between PHEVs and distribution grids[END_REF] determined by the following SDE,

dz t =   µ 1,t µ 2,t   dt -   β(t, z t ) 1 -β(t, z t )   g t     1 1   + σ t dW t   + dN t October 1, 2011 DRAFT
for a given initial state z 0 . In particular, W t = (W 1,t , W 2,t ) T is a two-dimensional Brownian motion with independent components. Similar to previously, σ t determines the variance of the noise at time t and the term dN t is a reflective term in order to guarantee that z t ∈ [0, 1] 2 . As in the EV case, the analysis of the game reduces to the analysis of the behavior of a single player. The cost function L (k) L is now assumed identical to all players and is given by

L (z 0 , µ, m 0 ) = E T 0 (r t (m) + q t (µ t ) + s t (z t )) dt + ξ(z T ), (15)
where m ∈ M is the trajectory of the probability distributions, that is,

m = {m(t, •, •), 0 ≤ t ≤ T } and
M is the set of all possible trajectories m. The initial state conditions are z 0 ∈ [0, 1] 2 and m(0, •, •). The price for electricity is given by the function r 1,t : M → R + , with

r 1,t (m) = D(t, •) -1 g t 1 0 1 0 β(t, z 1 , z 2 )m(t, z 1 , z 2 )dz 1 dz 2 + d dt 1 0 1 0 z 1 m(t, z 1 , z 2 )dz 1 dz 2 . ( 16
)
The price for oil is still constant, given by r 2,t = r 2 .

The next section is dedicated to determining the MFE for this game.

D. Mean Field Analysis

Under the above game formulation, the optimal control problem which represents the equilibrium of the game formulates as

u(0, z 0 ) = inf µ∈A L (z 0 , µ, m 0 ) dz t =   µ 1,t µ 2,t   dt -   β(t, z t ) 1 -β(t, z t )   g t [dt + σ t dW t ] + dN t . (17)
We introduce the value function

v(u, z • ) = inf µ∈A E T u (r t (m) + q t (µ t ) + s t (z t )) dt + ξ(z T ) (18)
with

z u = z • .
As in the EV case, we consider the cost function as quadratic, that is,

q t (µ t ) = 1 2 q 1,t (µ 1,t ) 2 + 1 2 q 2,t (µ 2,t ) 2 , with (q 1,t , q 2,t ) ∈ R 2 .
October 1, 2011 DRAFT

The HJB equation is here given by

-∂ t v(t, x, y) = inf µt∈R 2 {µ 1,t r 1,t (m ) + µ 2,t r 2 + q t (µ t ) +( µ 1,t -g t β(t, x, y) )∂ x v(t, x, y) + (µ 2,t + g t (β (t, x, y) -1)) ∂ y v(t, x, y)} +f t (x, y) + 1 2 σ 2 t g 2 t [( β(t, x, y) 2 ∂ 2 xx v(t, x, y) +2β(t, x, y)(1 -β(t, x, y))∂ 2
x,y v(t, x, y)

+(1 -β(t, x, y)) 2 ∂ 2 yy v(t, x, y) , ( 19
)
where m is solution to the FPK equation

∂ t m (t, x, y) = -∂ x (µ 1,t -β (t, x, y) g t )m (t, x, y) -∂ y (µ 2,t + (β (t, x, y) -1)g t )m (t, x, y) + 1 2 g 2 t σ 2 t β (t, x, y) 2 ∂ 2 xx m (t, x, y) + (1 -β (t, x, y)) 2 ∂ 2 yy m (t, x, y) +2β (t, x, y) (1 -β (t, x, y))∂ 2 xy m (t, x, y) ,
where µ t = (µ 1,t , µ 2,t ), the cost minimizing control, is

µ 1,t = - 1 q 1,t (r 1,t (m ) + ∂ x v(t, x, y)) (20)
µ 2,t = - 1 q 2,t (r 2 + ∂ y v(t, x, y)) . ( 21
)
Assuming σ t = 0, we obtain more compact forms. In particular, after substitution of the expression of µ t , the HJB equation becomes

∂ t v (t, x, y) = 1 2q 1,t (∂ x v (t, x, y) + r 1,t (m )) 2 + 1 2q 2,t (∂ y v (t, x, y) + r 2 ) 2 +g t β (t, x, y) ∂ x v (t, x, y) +g t (1 -β (t, x, y))∂ y v (t, x, y) -f t (x, y) . ( 22
)
where m is solution to

∂ t m (t, x, y) = 1 q 1,t (∂ xx v ) + 1 q 2,t (∂ yy v ) + g t [∂ x β (t, x, y) -∂ y β (t, x, y)] m (t, x, y) + 1 q 1,t (r 1,t (m) + ∂ x v ) + β (t, x, y) g t ∂ x m (t, x, y) + 1 q 2,t (r 2 + ∂ y v ) + (1 -β (t, x, y))g t ∂ y m (t, x, y) (23)
with v the solution to [START_REF] Itô | Multiple Wiener integral[END_REF].

In particular, for β(t, z t ) = z1,t z1,t+z2,t , which we will use in Section IV, we have that

[∂ x β (t, x, y) -∂ y β (t, x, y)] = 1 z 1,t + z 2,t . (24)
IV. SIMULATIONS

In this section, we provide simulation results for the electrical vehicle schemes developed in Section II and Section III.

A. EV analysis

We first consider the scenario of Section II. We assume a realistic three-day scenario (t = 0 at midnight the first day and t = T = 1 seventy-two hours later) where players have an average consumption rate that depends on specific periods of the days. The scenario is typical of a Friday to Sunday energy consumption, with higher overall electricity consumption on Friday and different patterns of car usage on Friday than on Saturday and Sunday. Since it is difficult to provide a universal system parametrization, we will take arbitrary scalings in the energy consumption functions. Some insights on the impact of different scalings will be provided.

The car electricity consumption function g t is depicted in Figure 1, where we see in particular that consumption is higher on Friday and with a peak around 5pm, while consumption is lower on weekend days with different peak times. The variance σ 2 t on the consumption is taken equal to 0.01 at all time, ensuring a standard deviation of the order of 10%. The demand function D(t, p) is such that the price p is a quadratic function of the total electricity demand from both electrical vehicles and other electricity services. Specifically, we take here

p t = g t + d dt xm(t, x)dx + + d t 2
where d t stands for the demand of electricity in services other than electrical cars, with [x] + = max(x, 0).

We therefore assume that this demand is deterministic and is not altered by price evolution, which is unrealistic to some extent but helps appreciating the impact of EV consumption on prices. The function d t is depicted in dashed line in Figure 4, up to a constant corresponding to the total average EV consumption;

that is, the dashed line represents the total electricity consumption if EV consumption were distributed equally in time. For simplicity of understanding, we assume h t = 30 constant; that is, we do not consider that the car owners have any particular incentive to charge or discharge at some specific time periods. 1We take f (t, x) = (1 -x) 2 to impose consumers to keep a certain level of electricity in their batteries, and the boundary condition κ(x) = (1 -x) 2 in order to avoid large sales at the last minute. The initial condition on m(0, •) is a triangle distribution m 0 centered at 0.5 and with support [0.3, 0.7]. The boundary conditions on m and v are such that

∂ x m(0, •) = ∂ x m(1, •) = ∂ x v(0, •) = ∂ x v(1, •) = 0 in order to force the energy content to lie in [0, 1].
To solve the system of equations ( 7), ( 8) in (m, v), we proceed by solving sequentially the HJB and FPK equations using a simple fixed-point algorithm until convergence. We do not ensure here that this algorithm does converge, neither do we ensure that the solution obtained is the solution sought for. Using a finite difference method on a sampling of 144 points in the time axis (every 30min) and of 100 points in the battery level axis, the above scheme leads to the density evolution m depicted in Figure 2. A few observations can be already made from this figure. We easily observe daily sequences of increases and decreases of the average battery levels. We see in particular that during nighttime, the battery levels increase, indicating that energy is purchased at time and consumed during daytime. It is interesting to note that, due to the small variance σ 2 t that was chosen, the overall tendency is for m (t, •) to concentrate into a single mass when t → 1. This is a usual phenomenon which determines the steady state if time were to continue with constant values for all time-dependent system parameters.

From the expression of m , v , and the equations derived in Section II, it is now possible to obtain much information about the system. In particular, it is interesting to follow the electricity bought or sold by electrical vehicles at all time, that is the quantity

g t + d dt xm (t, x)dx
or the overall electricity consumption in the market given by

g t + d dt xm (t, x)dx + d t
and the price p t (m ) defined here as

p t (m ) = g t + d dt xm (t, x)dx + + d t 2 .
This is depicted in Figure 3, Figure 4 and in Figure 5, respectively.

We see first in Figure 3 that the peaks of electricity bought by electrical vehicles take place during the night where the overall demand is low, while they are at their lowest during peak demand periods.

This is a natural outcome of the fact that prices are high during peak demand periods. However, we also see that the difference of amplitude between lowest and highest purchases is not large. This is due to the fact that, while prices are high in peak demand periods, the EV owners still have a strong incentive not to find their batteries empty, driving them to keep buying electricity at peak periods. This behaviour can be hindered by relaxing the constraint f (t, x).

Of more interest is Figure 4, where the differences between electricity consumption with or without incentives on EV behaviour. This figure depicts in dashed line the overall energy consumption if the EV purchases were equally distributed in the three-day period (that is, with no incentive), and in plain line the overall consumption under our current assumptions. It is seen here that the price incentives on electricity purchases produces a much expected peak demand reduction in the critical day periods, and a simultaneous increase of consumption during low consumption periods. Note importantly that our analysis does not consider changes in d t when the price for electricity changes; only the part of electricity reserved for EV drives prices which in turn drive the EV behaviour, which is a natural assumption if different price conditions are applied to EV and other services. The price evolution is depicted in Figure 5, where it is seen in this setting that the price is mostly driven by the function d t .

B. PHEV analysis

In this second section, we wish to analyze the behavior of hybrid vehicles as described in Section III. Since solving three-dimensional differential equations is time-consuming, we only provide results for the time scale discretized in 12 samples and for the "spatial" scales discretized both in 16 samples. For each differential equation, the resolution is performed by iterating the resolution of the two-dimensional differential equations along time and electricity scales for each fixed oil tank level, and time and oil scales for each fixed battery level. Then the system of HJB and FPK differential equations is solved by further iterating a fixed point algorithm as in the previous section. For simplicity of interpretation, we consider here a time-independent scenario where both g t = 0.2 and (q 1,t , q 2,t ) = (125, 125) are constant with time. 2 We take the electricity price to be r 1,t = (D(t, r 1,t )) + + 0.5, where now the demand is solely due to the electricity being bought by PHEV; that is, we do not consider other sources of electricity consumption in order to focus on the oil/electricity interaction solely. The oil price is set to r 2,t = r 2 = 0.7. This is a natural choice as it is expected that an approximate quantity g t = 0.2 will be asked for at any time to cover for the energy consumed, hence a price for electricity r 1,t 0.7. We impose a constraint s t (z) = 20(2 -z 1 -z 2 ) 2 , where z = (z 1 , z 2 ) T . The relative consumption β of oil and electricity is proportional to the total quantity of energy, that is β(t, z) = z 1 /(z 1 + z 2 ) and therefore

1 -β(t, z) = z 2 /(z 1 + z 2 )
. We take σ t = 0 for simplicity. The boundary constraints are identical to those in the previous section. As for the terminal constraint on v, it imposes that v(T

, z) = ξ(z) = (10(2 -(z 1 + z 2 )) 2 .
We consider the scenario where m(0, •) is a (properly truncated and scaled) Gaussian distribution with mean (0.4, 0.6) T and covariance 0.02I 2 , with I 2 the 2×2 identity matrix. That is, we assume that, initially, most vehicles have more oil than electricity. This is depicted in Figure 6. We then let the system evolve freely under the above set of constraints. It is natural to guess that the overall behavior is a decrease of either or both quantities of oil and electricity to zero if the prices are too high, or an increase of either or both quantities to one, if the prices are more reasonable. What is interesting to observe is the trajectory jointly followed by the players. The resulting final distribution m (T, •) is depicted in Figure 7. What we observe in the aforementioned conditions is that the initial distribution has shifted towards an increase of both electricity and oil levels, with a stronger increase of the mean battery level. Another observation is that the distribution tends to stretch along the z 1 = z 2 diagonal in the figure, translating the fact that oil and electricity are seen almost as equivalent goods due to the loosely constraining energy cost policy.

Among the different further analyses, in Figure 8, we consider a section of the distribution of the optimal transaction policy µ 1,t and µ 2,t at time t = 0 + , for z 2,t = 0.5 and z 2,t = 0.9 (we remind that both µ 1,t and µ 2,t are functions of t, z 1,t and z 2,t ). That is, we observe the initial behavior of players with half-filled oil tanks and almost completely filled oil tanks. It is seen that, for users with a very low level of electricity, buying electricity is an appealing choice. This can be interpreted by the fact that, as few players are in strong need for energy, it is possible to acquire a large quantity of electricity at a reasonable price. Those players with low reserves of electricity are the main beneficiaries. For users with already a reasonable level of electricity though, electricity and oil are seen as equivalent goods. As a matter of fact, our results also show that, at time t = 0 + , the price of electricity equals r 1,t = 0.706 r 2 .

That is, the players with low electricity levels draw as much of the electricity overhead (compared to oil) as is needed to reach an equilibrium price with oil. Now, it is also observed that, for users with large quantities of oil, electricity becomes a compelling purchase in order to further increase the total quantity of energy (since f imposes z 1,t + z 2,t to be close to 2), hence a larger incentive for buying electricity when the battery level is not large. When both battery and tank levels are alike, we see that the quantity of electricity purchased is the same as the quantity of oil purchased. Obviously, from the very generic settings of both EV and PHEV problems, many more scenarios can be carried out so to evaluate the actual impact of the EV and PHEV on realistic smart grid scenarios. The simulations above and their interpretations only provide a framework of understanding of fully rational vehicle owner's behavior, which needs be reported to real-life conditions with extreme care. Battery level Rate of energy purchase µ 1,t , z 2,t = 0.5 µ 2,t , z 2,t = 0.5 µ 1,t , z 2,t = 0.9 µ 2,t , z 2,t = 0.9 Fig. 8. Optimal transactions at time t = 0 + for players with different oil and battery levels.

V. CONCLUSION

In this article, we proposed a game theoretical framework to model the behavior of electrical vehicle and hybrid electricity-oil vehicle owners aiming at selfishly minimizing their operating cost. As the number of selfish players is large, and players are assumed alike, we then turned the problem into a mean field game, for which we obtain the fundamental differential equations describing the mean field equilibrium of the game. Using numerical methods, we drew conclusions which give new insights on the way to optimize the electrical vehicle penetration in the future smart grid.

= 1

 1 be deterministic and known by EV k. The amount of energy stored in the battery of vehicle k at time t is denoted by x (k) t ∈ [0, 1] evaluated in energy units. Here, x for a fully charged battery. We denote by α (k) t

  the trajectories of the partial cost function on the horizon [0, T ]. The density measure m obtained at the MFE and the associated optimal trajectories of the partial cost function v at the MFE are solution of the following system of coupled stochastic partial differential equations (SPDE), for a given initial distribution m 0 = m(0, •) and an initial state x 0 ,

Fig. 1 .

 1 Fig. 1. Mean energy consumption gt of EV as a function of time.

Fig. 3 .

 3 Fig. 3. Electricity purchased by EV as a function of the time t.

Fig. 5 .

 5 Fig. 5. Evolution of the price pt(m ) as a function of the time t.

Fig. 6 .Fig. 7 .

 67 Fig.6. Initial distribution m(0, •) at time t = 0, as a function of both levels of battery and oil tank.

  Fig. 2. Density solution m (t, x) as a function of the time t and the battery level x.

				m (t, x)
		0.8
	Battery level x	0.4 0.6
		0.2
		0	0 6 12 18 24 30 36 42 48 54 60 66
			Time t (in hours)
		0.72
		0.71
	EV purchases	0.69 0.7
		0.68
		0.67	g t + d dt	xm (t, x)dx
		0 6 12 18 24 30 36 42 48 54 60 66 72
			Time (in hours)

October 1, 2011 DRAFT

Note that the determination of a correct ht is highly subjective and is better kept constant for the sake of interpretation.

Such a large value for the entries of ht is motivated by faster algorithm convergence reasons, although it inhibits as a counterpart fast variations of m along time.October 1, 2011 DRAFT