Sylvie Chambon

Antonio Moreno

Anand P Santhanam

Jannick P Rolland
email: jannick@odalab.ucf.edu

Elsa Angelini

Isabelle Bloch
email: isabelle.bloch@enst.fr

CT-PET Landmark-based Lung Registration Using a Dynamic Breathing Model

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Registration of multimodal medical images is a widely addressed topic and is important in many different domains, in particular for oncology and radiotherapy applications. We consider Computed Tomography (CT) and Positron Emission Tomography (PET) of thoracic regions, which provide complementary information about the anatomy and the metabolism of the human body (see Figure 1). Their registration has a significant impact on improving medical decisions for diagnosis and therapy [START_REF] Lavely | Automatic registration of PET and CT studies for clinical use in thoracic and abdominal conformal radiotherapy[END_REF][START_REF] Rizzo | Automatic registration of PET and CT studies for clinical use in thoracic and abdominal conformal radiotherapy[END_REF][START_REF] Vogel | Correction of an image size difference between positron emission tomography (PET) and computed tomography (CT) improves image fusion of dedicated PET and CT[END_REF]. Linear registration is not sufficient to cope with the local deformations produced by cardiac and respiratory motions. Therefore, non-linear registration methods are required to register multimodality images of thoracic and abdominal regions, even with combined PET/CT scanners [START_REF] Shekhar | Automated 3-Dimensional Elastic Registration of Whole-Body PET and CT from Separate or Combined Scanners[END_REF].

Most of the existing non-linear registration methods use intensity information or features in order to calculate the transformation between the images [START_REF] Maintz | A Survey of Medical Image Registration[END_REF][START_REF] Pluim | Image Registration[END_REF][START_REF] Zitovà | Image Registration Methods: A Survey[END_REF]. Thus they have to either find the transformation that maximizes the similarity between the registered image and the target image (iconic methods) or compute a transformation that matches some particular features (landmarks) in both images (geometrical methods). In the case of landmark-based methods, the selection of these particular features is an important task. In many methods, the curvature of the surfaces to reg-ister is used. However, few papers included a detailed study of the selection of the landmarks [START_REF] Betke | Landmark detection in the chest and registration of lung surfaces with an application to nodule registration[END_REF][START_REF] Rohr | Landmark-based elastic registration using approximating thin-plate splines[END_REF]. Moreover, most registration methods are based on image information, and do not take into account the physiology of the human body. However, physiological information can be useful in order to ensure realistic deformations and to guide the registration process. While several papers present breathing models built for medical visualization [START_REF] Zordan | Breathe Easy: Model and Control of Human Respiration for Computer Animation[END_REF], no paper exploits such a model in a registration process. Consequently, in this paper, we propose an approach in which we integrate a physiologically driven breathing model in a non-linear registration procedure based on lung surface landmark points in order to guarantee physiologically plausible deformations.

In Section 2, we summarize existing works which use breathing models combined (or not) with registration algorithms and then we provide an overview of the selected model. The proposed model-based non-linear registration algorithm is detailed in Section 3. Then, the application of a landmark-based registration method adapted to pathological cases combined with the breathing model is described in Section 4. Section 5 discusses some results.

Breathing Models

Thoracic Imaging Registration

A recent study highlighted the effects of breathing during a non-rigid registration process and the importance of taking these into account [START_REF] Sundaram | Towards a Model of Lung Biomechanics: Pulmonary Kinematics Via Registration of Serial Lung Images[END_REF]. Three techniques for respiration-gated radiotherapy are being developed to improve the efficiency of lung tumors radiations:

(1) active: controling the patient's breathing via airflow blockage;

(2) passive or empirical: using external measurements in order to adapt radiation protocols to the tumor's motion [START_REF] Mcclelland | A Continuous 4D Motion Mmodel from Multiple Respiratory Cycles for Use in Lung Radiotherapy[END_REF][START_REF] Nehmeh | Four-dimensional (4D) PET/CT imaging of the thorax[END_REF][START_REF] Wolthaus | Fusion of respiration-correlated pet and ct scans: correlated lung tumour motion in anatomical and functional scans[END_REF]; (3) model-based: employing a breathing model to predict lungs deformations during the breathing cycle [START_REF] Sarrut | Non-rigid registration method to assess reproducibility of breath-holding with ABC in lung cancer[END_REF].

We focus on thoracic volume registration and propose a "patient-specific" registration through the use of a breathing model (technique (3)). Different bio-mathematical representations of the human respiratory mechanics have been developed [START_REF] Mead | Measurement of Inertia of the Lungs at Increased Ambient Pressure[END_REF]. Mathematical tools can be employed and the most popular technique, for medical visualization, is called NCAT (NURBS-based cardiac-torso). It is based on Non-Uniform Rational B-Spline (NURBS) to correct for respiratory artifacts of SPECT images [START_REF] Segars | Study of the Efficacy of Respiratory Gating in Myocardial SPECT Using the New 4-D NCAT Phantom[END_REF]. A multiresolution registration for 4D Magnetic Resonance Imaging (MRI) was proposed in [START_REF] Rohlfing | Modeling Liver Motion and Deformation During the Respiratory Cycle Using Intensity-Based Free-Form Registration of Gated MR Images[END_REF]. In [START_REF] Guerrero | Elastic image mapping for 4-D dose estimation in thoracic radiotherapy[END_REF], a 4D NCAT phantom and a 3D CT image were used to generate 4D CT and to compute an elastic registration. Physically-based models describe the important role of airflow inside the lungs and can be based on Active Breathing Coordinator (ABC) [START_REF] Sarrut | Non-rigid registration method to assess reproducibility of breath-holding with ABC in lung cancer[END_REF] or volume preservation relation [START_REF] Santhanam | Modeling, Simulation, and Visualization of 3D Lung Dynamics[END_REF][START_REF] Zordan | Breathe Easy: Model and Control of Human Respiration for Computer Animation[END_REF]. Authors of [START_REF] Rohlfing | Modeling Liver Motion and Deformation During the Respiratory Cycle Using Intensity-Based Free-Form Registration of Gated MR Images[END_REF][START_REF] Sundaram | Towards a Model of Lung Biomechanics: Pulmonary Kinematics Via Registration of Serial Lung Images[END_REF] used pre-register MRI to estimate the breathing model. In [START_REF] Guerrero | Elastic image mapping for 4-D dose estimation in thoracic radiotherapy[END_REF], the NCAT phantom was used, but, from a modeling and simulation point of view, physicallybased deformation methods are better adapted for simulating lung dynamics as they allow precise generation of intermediate 3D lung shapes. They are easy to adapt to patients, without the need for physical external adaptations for each radiotherapy treatment.

Physics-Based Breathing Model

The modeling approach used in this work was previously discussed in [START_REF] Santhanam | Modeling, Simulation, and Visualization of 3D Lung Dynamics[END_REF] and the two major components involved in the modeling and visualization efforts include:

(1) the parameterization of Pressure-Volume (PV) data of a human subject which act as an ABC; (2) the estimation of the deformation operator from either 4D (3D+t) CT lung data or two 3D CT lung data set.

In step (1) a parameterized PV curve, obtained from normal human subjects, is used as a driver for simulating the 3D lung shapes at different lung volumes. In step (2), the computation takes as inputs the nodal displacements of the 3D lung model and the estimated amount of force applied on the nodes of the meshes (which are on the surface). Displacements are obtained from 4D CT of a normal human subject. The direction and magnitude of the lung surface point's displacement are computed using the fact that the expansion of lung tissues is linearly related to the increase in lung volume. The estimated amount of applied force on each node (that represents the air-flow inside lungs) is estimated based on a PV curve and the lungs's orientation with respect to the gravity, which controls the air flow. Given these inputs, a physics-based deformation approach based on Green's function (GF) formulation is estimated to deform the 3D lung surface models. Specifically the GF is defined in terms of a physiological factor, the regional alveolar expandability (elastic properties), and a structural factor, the inter-nodal distance of the 3D surface lung model.

To compute the coefficients of these two factors, an iterative approach is employed and, at each step, the force applied on a node is shared with its neighboring nodes based on a local normalization of the alveolar expandability coupled with inter-nodal distance. The process stops when sharing of the applied force reaches equilibrium [START_REF] Santhanam | Modeling Simulation and Visualization of Real-Time 3D Lung Dynamics[END_REF].

Using the Breathing Model in Registration

We have conceived an original algorithm in order to introduce the advantages of using the breathing model described above in a registration procedure. Figure 2 shows the computational workflow of the complete algorithm. The inputs consist of one PET volume and two CT volumes of the same patient, corresponding to two different instants of the breathing cycle (intermediate expirations). The first step of the algorithm consists in segmenting the lungs and the tumors on the PET data and on the two CT data sets, using a robust mathematical-morphology-based approach [START_REF] Camara | Explicit Incorporation of Prior Anatomical Information into a Non-Rigid Registration of Thoracic and Abdominal CT and 18-FDG Whole-Body Emision PET Images[END_REF]. The meshes (called CT mesh and PET mesh) corresponding to the different segmentation results are computed. The subsequent steps are detailed next.

Patient-Specific Breathing Model

We first estimate the intermediate 3D lung shapes between the two segmented CT lung datasets. Displacements of lung surface points are computed as follows:

(1) Directions are given by the model (computed from a 4D CT normal data set of reference). (2) Magnitudes are "patient-specific" and are computed from the given 3D CT lung datasets.

In other words, for known directions of displacement the magnitude of the displacement is computed from the two 3D CT lung datasets. With known estimations of applied force and "subject-specific" displacements the coefficients of the GF can be estimated (cf. Section 2.2). Then, the GF operator is used to compute the 3D lung shapes at different intermediate lung volumes. This methodology is further detailed in [START_REF] Santhanam | Modeling, Simulation, and Visualization of 3D Lung Dynamics[END_REF].

CT Selection

By applying the continuous breathing model, we can obtain different instants ("snapshots") of the breathing cycle, generating simulated CT meshes. By comparing each CT mesh with the PET mesh, we select the "closest" one. Let us denote the CT simulated meshes by M 1 , M 2 ,. . . , M N . The mesh M N corresponds to the CT in maximum inhalation and M 1 to maximum exhalation. By using the breathing model, the transformation φ i j between any two instants i and j of the breathing cycle can be computed as: M j = φ i j (M i ). We can compare these CT meshes with the PET mesh (M P ET ). We define a measure of similarity between meshes (or their corresponding volumes) and the mesh that minimizes the criterion (C) is denoted as M C :

M C = arg min i C(M i , M P ET ). (1)
The Root Mean Square (RMS) distance has been chosen as the criterion C, as a first approach:

D RMS (M, A) = 1 2 [d RMS (M, A) 2 + d RMS (A, M ) 2 ] with d RMS (M, A) = 1 |M| p∈M D(p, A)

Deformation of the PET

Registration of M P ET and the original CT mesh M N can be perform with two approaches : (1) A direct registration (dashed line in Figure 3):

M Rd P ET (N ) = f Rd (M P ET , M N ), (2)
where f Rd denotes the transformation that registers directly M P ET and M N , and M Rd P ET the result of registering the PET directly to the CT mesh M N . The transformation f Rd may be computed by any registration method adapted to the problem. As an illustrative example, we choose the original CT to correspond to the end-inspiration, M N , but a similar process could be applied for any CT image. In this approach the deformation itself is not guided by any anatomical knowledge. In addition, if the PET and the original CT are very different (end-inspiration CT), it is likely that this registration will provide physically unrealistic results.

(2) To avoid such potential problems, we propose here an alternative approach: once the appropriate CT (M C ) is selected, we compute the registration between the M P ET mesh and the M C mesh as:

M r P ET (C) = f r (M P ET , M C ), (3)
where f r is the registration transformation and M r P ET denotes the registered mesh. Then, the transformation due to the breathing is used to register the M r P ET to the original CT (continuous line in Figure 3). The transformation due to the breathing between M C and M N can be computed as the following composition:

Φ C,N = φ N -1,N • . . . • φ C+1,C+2 • φ C,C+1 . (4)
We apply to M r P ET the same transformation Φ C,N in order to compute the registration with M N :

M Rbm P ET (N ) = Φ C,N (M r P ET (C)) = Φ C,N (f r (M P ET , M C )), (5)
where M Rbm P ET denotes the PET registered mesh using the breathing model. We can register MP ET to MN following one of the two paths (notations are defined in Section 3.3).

Registration Adapted to Pathologies

The algorithm described in Section 3 can be applied with any type of registration method. We show how the proposed approach can be applied for landmark-based registration of multimodality images in pathological cases, in particular for diagnosis, follow-up and radiotherapy treatments.

Influence of Selected Landmark Points

Features selection is an important task in registration. In this section, we focus on voxel selection but more complex features can be detected [START_REF] Beil | Investigation of approaches for the localization of anatomical landmarks in 3D medical images[END_REF]. The selection can be manual (as in most methods) [START_REF] West | Hybrid point-and-intensity-based deformable registration for abdominal ct images[END_REF], semi-automated [START_REF] Rohr | Landmark-based elastic registration using approximating thin-plate splines[END_REF], or automated [START_REF] Rohr | Landmark-based elastic registration using approximating thin-plate splines[END_REF]. Manual selection is tedious and timeconsuming. The authors in [START_REF] Hartkens | Using Points and Surfaces to Improve Voxel-Based Non-Rigid Registration[END_REF] suggest that semi-automated selection is interesting, integrating experts knowledge in an automatic process. Automatic selection permits reduced execution time with high accuracy. Most of these automatic

Registration with Rigidity Constraints and Breathing Model

Once the different CT meshes are computed and the closest CT mesh, M C , is selected, the PET and the original CT (in our example M N ), are registered as follows:

(1) Select landmarks on the CT mesh M C (with Gaussian or/and mean curvatures); the PET and the original CT using the estimated correspondences.

In step (1), the four variants presented in Section 4.1 can be used. The breathing model is used in step ( 4). The landmark points selected on M C are tracked on the meshes estimated with the breathing model. Consequently, we can assume that the corresponding landmarks selected on the original CT are correct (and actually they represent the same anatomical point) and follow the deformations of the lungs during the respiratory cycle.

Results and Discussion

In this section we present some results we have obtained using the general methodology described in Section 3 and the registration method summarized in Section 4.

We have applied our algorithm on a pathological case exhibiting one tumor and on a normal case. We have one PET and two CT images for each case. First, the breathing model is computed using the meshes of the lungs segmented on CT data. Then, we compare 10 (regularly distributed) instants of the generated model with the lung surface segmented from the PET (meshes have more than 40 000 nodes). The results obtained with the proposed algorithm are physically-based and more realistic than results obtained by registering the PET directly with the original CT. First results confirm this statement as shown in Figures 6 and7. It can be observed that the result of the registration by a direct method (Figure 6(d)) produces unrealistic deformations in the region between the lungs. With the proposed algorithm (Figure 6(h)), the result is visually more accurate. The RMS distance between the registered PET lungs and the original CT lungs is reduced to 11.8 mm. In Figure 7 the improvement of the results is clearly illustrated for the normal case in the region of the right lung close to the liver. All these results have been obtained by using landmarks determined by the combination of mean and Gaussian curvatures plus a uniform selection. This variant provided visually better results although further validation is necessary.

(a) (b) (c) (d) (e) (f) (g) (h)

Conclusion and Future Work

We have developed a CT-PET landmark-based registration method that uses a breathing model to guarantee phys-iologically plausible deformations. The method consists in computing a deformation guided by a breathing model. We also proposed and studied four variants, based on curvature, to select landmarks. Initial results on a pathological case and a normal case are very promising and show the improvement brought by the breathing model. Our algorithm avoids undesired tumor misregistrations and preserves tumor geometry and intensity. Moreover, as the tumor in CT and PET has not necessarily the same size and shape, the registration of these two modalities is very useful because all the information of the PET image is preserved. This is very important in order to know the true extension of the pathology for diagnosis and for the treatment of the tumor with radiotherapy, for example. Moreover, we highlight the best variant to detect landmarks: to uniformly select landmarks by combining mean and Gaussian curvatures.

We are currently performing a deeper evaluation on a larger database, in collaboration with clinicians. Future work includes a refined "snapshot" selection, using further subdivisions of time intervals, a more precise characterization of the tumor movement and its influence on the breathing, and a comparison of the the proposed method with other approaches which include rigid objects in the registration.

Figure 1 .

 1 Figure 1. Corresponding views of the same patient in CT (a)-(b) for two instants of the breathing cycle, and PET (c).

Figure 2 .

 2 Figure 2. Computational workflow of the registration of CT and PET images using a breathing model (cf. Section 3.3).

2

 2 and where D(p, A) = [min q∈A d(p, q)] with d the Euclidean distance.

Figure 3 .

 3 Figure 3. The mesh MC is the closest to the mesh MP ET .We can register MP ET to MN following one of the two paths (notations are defined in Section 3.3).

( 2 )

 2 Estimate corresponding landmarks on the PET (using the Iterative Closest Point algorithm [2]); (3) i = C; (4) Track landmarks from M i to the next CT mesh M i+1 ; (5) If M i+1 = M N , go to step (6) else go to step (4) with i = i + 1; (6) Register (with the method summarized in Section 4.2)

  Figure 5 shows the results of surface comparison between the PET surface, for two instants from the CT data set: the closest and the end-inspiration.

2 Figure 5 .

 25 (a) DRMS = 12.1 (b) DRMS = 24.Superimposition of the contours of the PET (black) and the CT lungs (grey) at two instants of the breathing cycle: (a) closest (MC) and (b) end-inspiration (MN ).

Figure 6 .Figure 7 .

 67 Figure 6. Original PET (a) and CT (b) images with tumor (surrounded by a white circle). Correspondences between the selected points in the PET image and in the CT image are shown in (c), (e) and (g). Registered PET data is shown in (d) for the direct method and in (f) and (h) for the method with the breathing model with the variant combining mean and Gaussian curvatures in (e)-(f) and adding uniformly distributed landmarks in (g)-(h). Illustrations are provided in 2D for the sake of readability.

Acknowledgments: This work has been partially funded by an ANR grant: ANR-05-BLANC-0081.

In the large rectangle, there is no landmark with GAU method whereas there are four landmarks with the MEA method.

In the fusion method (MEA-GAU), these landmarks are selected. In the small rectangle, no landmark is selected with the mean and/or the Gaussian curvatures. However, a landmark is added in this area with the MEA-GAU-UNI method.

methods exploit curvature [START_REF] Rohr | Landmark-based elastic registration using approximating thin-plate splines[END_REF]. In [START_REF] Betke | Landmark detection in the chest and registration of lung surfaces with an application to nodule registration[END_REF], an auto-correlation method is also combined with curvature.

In the present work, landmark selection is automatic and based on Gaussian and mean curvatures, according to the following steps: (1) compute curvature for each voxel of the lung surface;

(2) sort voxels in decreasing order of absolute value of curvature; (3) select voxels based on curvature and distance criteria (detailed in the following paragraph); (4) if a uniform selection is needed then add voxels with zero-curvature in the area where no voxels have been considered as landmarks. This algorithm is proposed to select particular voxels that provide relevant information. Moreover, we intend to obtain an approximately uniform selection to take into account the entire surface of the lungs for computing the deformation. In step 3, we consider V = {v i } i=0..NS , the set of voxels in decreasing order of absolute value of curvature, where N S is the number of voxels of the surface and V L = {v Li } i=0..NL , the set of landmarks, where N L is the number of landmarks. For each voxel v i ∈ V (for i = 0 to N S ) with non-zero-curvature, we add

where d g is the geodesic distance on the lung surface and T is a threshold to be chosen. With this selection process, some regions (the flattest) may contain no landmark, hence the addition of step 4: for each voxel on the surface of the lung v i ∈ V with zero-curvature, if there is no voxel v j ∈ V L with d g (v i , v j ) < T , we add v i in V L .

Four variants are tested: (1) MEA -Mean curvature without step 4;

(2) GAU -Gaussian curvature without step 4;

(3) MEA-GAU -Using mean and Gaussian curvature without step 4; (4) MEA-GAU-UNI -Using mean and Gaussian curvature with step 4. When mean and Gaussian curvatures are employed (methods MEA-GAU and MEA-GAU-UNI), the set V merges the set of voxels in decreasing order of mean curvature and the set of voxels in decreasing order of Gaussian curvature, by taking alternatively a value in each set. These strategies for landmark point selection are compared in Figure 4. Results given by the MEA and GAU methods are different, and it is interesting to combine them (see the results obtained with the MEA-GAU method). The MEA-GAU-UNI method permits to add some points in locally flat regions (see Figure 4).

Rigidity Constraints in Pathological Cases

We have developed a registration algorithm for the thoracic region in the presence of pathologies [START_REF] Moreno | Non-linear Registration Between 3D Images Including Rigid Objects: Application to CT and PET Lung Images With Tumors[END_REF]. The advantage of our approach is that it takes into account the tumors, while preserving continuous smooth deformations. We assume that the tumor is rigid and thus a linear transformation is sufficient to cope with its movements between CT and PET images. This hypothesis is relevant and in accordance with the clinicians' point of view, since tumors are often compact masses of pathological tissue. The algorithm relies on previously segmented structures (lungs and tumors). Landmarks corresponding to homologous points are defined in both images, and will guide the deformation of the PET image towards the CT image. The deformation at each point is computed using an interpolation procedure based on the landmarks, the specific type of deformation of each landmark (depending on the structure it belongs to), and weighted by a distance function, which guarantees that the transformation is continuous. We have shown that a consistent and robust transformation is obtained [START_REF] Moreno | Non-linear Registration Between 3D Images Including Rigid Objects: Application to CT and PET Lung Images With Tumors[END_REF].


Antonio Moreno

Sylvie Chambon

Anand P Santhanam

Roberta Brocardo

Patrick Kupelian

Jannick P Rolland

Elsa Angelini

Isabelle Bloch

Thoracic CT-PET Registration Using a 3D Breathing Model

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Registration of multimodal medical images is a widely addressed topic in many different domains, in particular for oncology and radiotherapy applications. We consider Computed Tomography (CT) and Positron Emission Tomography (PET) in thoracic regions, which provide complementary information about the anatomy and the metabolism of the human body (Fig. 1). Their registration has a significant impact on improving medical decisions for diagnosis and therapy [START_REF] Lavely | Phantom validation of coregistration of PET and CT for imageguided radiotherapy[END_REF][START_REF] Rizzo | Automatic registration of PET and CT studies for clinical use in thoracic and abdominal conformal radiotherapy[END_REF][START_REF] Vogel | Correction of an image size difference between positron emission tomography (PET) and computed tomography (CT) improves image fusion of dedicated PET and CT[END_REF]. Linear registration is not sufficient to cope with local deformations produced by respiration. Even with combined PET/CT scanners which avoid differences in patient orientation and provide linearly registered images, non-linear registration remains necessary to compensate for cardiac and respiratory motions [START_REF] Shekhar | Automated 3-Dimensional Elastic Registration of Whole-Body PET and CT from Separate or Combined Scanners[END_REF].

Most of the existing non-linear registration methods are based on image information and do not take into account any knowledge of the physiology of the human body. Landmark-based registration techniques do take physiology into account by forcing homologous points to match. In this direction, several breathing models were built for medical visualization, for correcting artefacts in images or for estimating lung motion for radiotherapy applications, but few papers exploit such models in a registration process (Section 2).

In this paper, we propose to integrate a physiologically driven breathing model into a 3D non-linear registration (Sections 3 and 4). The registration problem is defined between two CT volumes and one PET volume (Fig. 1).

Breathing models

Breathing Models and Thoracic Imaging Registration -Currently, respirationgated radiotherapies are being developed to improve the efficiency of radiations of lung or abdominal tumors [START_REF] Sarrut | Deformable registration for image-guided radiation therapy[END_REF]. Three techniques have been proposed so far: (i) active techniques controling the patient's breathing (airflow is blocked); (ii) passive or empirical techniques using external measurements in order to adapt radiation protocols to the tumor's motion [START_REF] Mcclelland | A Continuous 4D Motion Model from Multiple Respiratory Cycles for Use in Lung Radiotherapy[END_REF][START_REF] Nehmeh | Four-dimensional (4D) PET/CT imaging of the thorax[END_REF][START_REF] Wolthaus | Fusion of respiration-correlated PET and CT scans: correlated lung tumour motion in anatomical and functional scans[END_REF]; (iii) model-based techniques employing a breathing model to evaluate lungs deformations during the breathing cycle [START_REF] Sarrut | Non-rigid registration method to assess reproducibility of breathholding with ABC in lung cancer[END_REF]. Different bio-mathematical representations of the human respiratory mechanics have been developed [START_REF] Mead | Measurement of Inertia of the Lungs at Increased Ambient Pressure[END_REF]. Among Mathematical tools, the most popular technique is based on Non-Uniform Rational B-Spline (NURBS), surfaces which are bidirectional parametric representations of an object. In [START_REF] Segars | Study of the Efficacy of Respiratory Gating in Myocardial SPECT Using the New 4-D NCAT Phantom[END_REF], NURBS surfaces were used to correct for respiratory artifacts of SPECT images, building NCAT (NURBS-based cardiac-torso) model. A multi-resolution registration approach for 4D Magnetic Resonance Imaging (MRI) was proposed in [START_REF] Rohlfing | Modeling Liver Motion and Deformation During the Respiratory Cycle Using Intensity-Based Free-Form Registration of Gated MR Images[END_REF] with NCAT. In [START_REF] Guerrero | Elastic image mapping for 4-D dose estimation in thoracic radiotherapy[END_REF], a 4D NCAT phantom and an original CT image were used to generate 4D CT and to compute an elastic registration. Physically-based models, describing the important role of airflow inside the lungs, can be based on Active Breathing Coordinator (ABC allows clinicians to pause the patient's breathing at a precise lung volume) [START_REF] Sarrut | Non-rigid registration method to assess reproducibility of breathholding with ABC in lung cancer[END_REF] or on volume preservation relations [START_REF] Zordan | Breathe Easy: Model and Control of Human Respiration for Computer Animation[END_REF][START_REF] Santhanam | Modeling, Simulation, and Visualization of 3D Lung Dynamics[END_REF]. In [START_REF] Pollari | Evaluation of cardiac PET-MRI registration methods using a numerical breathing phantom[END_REF], segmented MRI data were used to simulate PET volumes at different instants of the breathing cycle. These estimated PET volumes were used to evaluate different PET/MRI registration processes. Authors of [START_REF] Rohlfing | Modeling Liver Motion and Deformation During the Respiratory Cycle Using Intensity-Based Free-Form Registration of Gated MR Images[END_REF][START_REF] Sundaram | Towards a Model of Lung Biomechanics: Pulmonary Kinematics Via Registration of Serial Lung Images[END_REF] used pre-register MRI to estimate a breathing model. CT registration using a breathing model was presented in [START_REF] Sarrut | Non-rigid registration method to assess reproducibility of breathholding with ABC in lung cancer[END_REF] but a specific equipement is needed. From a modeling and simulation point of view, physically-based deformation methods are better adapted for simulating lung dynamics and are easy to adapt to the patient, without the need for physical external adaptations.

Physics-Based Dynamic 3D Surface Lung Model -We employ an approach that was previously discussed in [START_REF] Santhanam | Modeling, Simulation, and Visualization of 3D Lung Dynamics[END_REF] and in which the two major components involved in the modeling efforts include: (1) Parameterization of PV (Pressure Volume) data from a human subject which acts as an ABC; (2) Estimation of the deformation operator from 4D CT lung data sets. In step (1) a parameterized PV curve, obtained from a normal human subject, is used as a driver for simulating the 3D lung shapes at different lung volumes. In step (2), the computation takes as inputs the nodal displacements of the 3D lung models and the estimated amount of force applied on the nodes of the meshes (which are on the surface of the lungs). Displacements are obtained from 4D CT data of a normal human subject. The direction and magnitude of the lung surface point's displacement are computed using the volume linearity constraint, i.e. the fact that the expansion of lung tissues is linearly related to the increase in lung volume. The amount of applied force on each node (that represents the air-flow inside lungs) is estimated based on a PV curve and the lungs's orientation with respect to the gravity, which controls the air flow. Given these inputs, a physics-based deformation approach based on Green's function (GF) formulation is estimated to deform the 3D lung surface models. Specifically the GF is defined in terms of a physiological factor, the regional alveolar expandability (elastic properties), and a structural factor, the inter-nodal distance of the 3D surface lung model. To compute the coefficients of these two factors, an iterative approach is employed and, at each step, the force applied on a node is shared with its neighboring nodes, based on local normalization of the alveolar expandability, coupled with inter-nodal distance.

The process stops when this sharing of the applied force reaches equilibrium.

For validation purposes, a 4D CT dataset of a normal human subject with four instances of deformation was considered [START_REF] Santhanam | Modeling Simulation and Visualization of Real-Time 3D Lung Dynamics[END_REF]. The simulated lung deformations matched the 4D CT dataset with 2 mm average distance error.

Combining Breathing Model and Image Registration

We have conceived an original algorithm in order to incorporate the breathing model described above in our multimodal image registration procedure. Fig. 2 shows the complete computational workflow. The input consists of one PET volume and two CT volumes of the same patient, corresponding to two different instants of the breathing cycle (end-inspiration and end-expiration, for example, collected with breath-hold maneuver). The preliminary step consists in segmenting the lung surfaces (and, eventually, the tumors) on the PET data and on the two CT data sets, using a robust mathematical-morphology-based approach [START_REF] Camara | Explicit Incorporation of Prior Anatomical Information into a Nonrigid Registration of Thoracic and Abdominal CT and 18-FDG Whole-Body Emision PET Images[END_REF], and extracting meshes corresponding to the segmented objects. Computation of a Patient-Specific Breathing Model -For each patient, we only have two segmented CT datasets, therefore we first estimate intermediate 3D lung shapes between these two datasets and then, the displacements of lung surface points. Directions are given by the model (computed from a 4D CT normal data set of reference) while magnitudes are "patient-specific" (computed from the given 3D CT lung datasets). With known estimations of applied force and "subject-specific" displacements the coefficients of the GF can be estimated (Section 2). Then, the GF operator is used to compute the 3D lung shapes at different intermediate lung volumes.

CT Selection -Let us denote the CT simulated meshes M 1 , M 2 ,. . . , M N with M 1 corresponding to the CT in maximum exhalation and M N to maximum inhalation. By using the breathing model, the transformation φ i,j between two instants i and j of the breathing cycle can be computed as: M j = φ i,j (M i ). Our main assumption is that even if the PET volume represents an average volume throughout the respiratory cycle, using a breathing model, we can compute a CT volume that can be closer to the PET volume than the original CT volumes. By applying the continuous breathing model, we generate simulated CT meshes at different instants ("snapshots") of the breathing cycle. By comparing each CT mesh with the PET mesh (M P ET ), we select the "closest" one (i. e. with the most similar shape). The mesh that minimizes a measure of similarity C (here the root mean square distance) is denoted as

M C : M C = arg min i C(M i , M P ET ).
Deformation of the PET -Once the appropriate CT (M C ) is selected, we compute the registration, f r , between the M P ET mesh and the M C mesh as:

M r P ET (C) = f r (M P ET , M C ), (1)
where M r P ET (C) denotes the registered mesh. Then, the transformation due to the breathing is used to register the PET to the original CT (continuous line in Fig. 3) incorporating the known transformation between M C and M N :

Φ C,N = φ N -1,N • . . . • φ C+1,C+2 • φ C,C+1 . (2)
We apply Φ C,N to M r P ET (C) in order to compute the registration with M N :

M Rbm P ET (N ) = Φ C,N (M r P ET ) = Φ C,N (f r (M P ET , M C )), (3)
where M Rbm P ET (N ) denotes the PET registered mesh using the breathing model. A direct registration, denoted f Rd , can also be computed between M P ET and the original CT mesh M N (dashed line in Fig. 3): M Rd P ET (N ) = f Rd (M P ET , M N ), where M Rd P ET (N ) is the result of registering the PET directly to the CT mesh M N (note that this could be done with another instant M i ). In the direct approach the deformation itself is not guided by any anatomical knowledge. In addition, if the PET and the original CT are very different, it is likely that this registration procedure will provide physically unrealistic results.

Registration Method Adapted to Pathologies

The algorithm described in Section 3 can be applied with any type of registration method, to estimate f Rd and f r . These functions may be computed by any registration method adapted to the problem. We show here how the proposed approach can be adapted for registration of multi-modality images in pathological cases.

Registration with Rigidity Constraints -We have previously developed a registration algorithm for the thoracic region taking into account the presence of tumors, while preserving continuous smooth deformations [START_REF] Moreno | Non-linear Registration Between 3D Images Including Rigid Objects: Application to CT and PET Lung Images With Tumors[END_REF]. We assume that the tumor is rigid and that a linear transformation is sufficient to cope with its displacements between CT and PET scanning. This hypothesis is relevant and in accordance with the clinicians' point of view, since tumors are often compact masses of pathological tissue. The registration algorithm relies on segmented structures (lungs and tumors). Landmark points are defined on both datasets to guide the deformation of the PET volume towards the CT volume. The deformation at each point is computed using an interpolation procedure where the specific type of deformation of each landmark point depends on the structure it belongs to, and is weighted by a distance function, which guarantees continuity of the transformation.

Registration with Rigidity Constraints and Breathing Model -Here, the following procedure is used to compute f r (in our example M N is the original CT):

1. Selection of landmark points on the CT mesh M C (based on Gaussian and mean curvatures and uniformly distributed on the lung surface) [START_REF] Chambon | CT-PET Landmark-based Lung Registration Using a Dynamic Breathing Model[END_REF];

2. Estimation of corresponding landmark points on the PET mesh M P ET (using the Iterative Closest Point (ICP) algorithm [START_REF] Besl | A Method for Registration of 3-D Shapes[END_REF]); 3. Tracking of landmark points from M C to the CT mesh M N using the breathing model; 4. Registration of the PET and the original CT using the estimated correspondences with the method summarized in the previous paragraph.

The breathing model used in step [START_REF] Vogel | Correction of an image size difference between positron emission tomography (PET) and computed tomography (CT) improves image fusion of dedicated PET and CT[END_REF] guarantees that the corresponding landmarks selected on the original CT are correct (and actually they represent the same anatomical point) and follow the deformations of the lungs during the respiratory cycle.

Results and Discussion

We The breathing model was initialized using the lung meshes from the segmented CT. Ten meshes (corresponding to regularly distributed instants) are generated and compared with the PET. The computation time can reach two hours for the whole process (a few seconds for segmentation, a few minutes for landmark point selection and about ninety minutes for registration). Although this is not a constraint because we do not deal with an on-line process, this computation time will be optimized in the future. As illustrated in Fig. 4 and 5 (one normal case and one pathological case), the correspondences between landmark points on the original CT and the PET are more realistic in the results obtained with the breathing model (images (e) and (f)) than without (images (b) and (c)). Using the model, it can be observed that the corresponding points represent the same anatomical points and that the uniqueness constraint is respected, leading to visually better looking PET registered images. In particular, the lower part of the two lungs is better registered using the model, the lung contour in the registered PET is closer to the lung contour in the original CT, cf. Fig. 4(g-i). In the illustrated pathological case, the tumor is well registered and not deformed. Moreover, the distance between the registered PET lungs and the original CT lungs is lower than using the direct approach.

In this paper, we consider the impact of the physiology on lung surface deformation, based on reference data of normal human subjects. Therefore the methodology presented in this paper will further benefit upon the inclusion of patho-physiology specific data once established. The use of normal lung physiology serves to demonstrate improvements in CT and PET registration using a physics-based 3D breathing lung model. Current work includes a quantitative comparison and evaluation on a larger database, in collaboration with clinicians. Acknowledgments: This work has been partially funded by an ANR grant: ANR-05-BLANC-0081.

Fig. 1 .

 1 Fig. 1. CT images (a,b) corresponding to two different instants of the breathing cycle and PET image (c) of the same patient (coronal views).

Fig. 2 .

 2 Fig. 2. Registration of CT and PET volumes using a breathing model.

Fig. 3 .

 3 Fig.3. Registration framework on PET (MP ET ) and CT mesh (MN ) -The MC mesh is the closest to the MP ET mesh. We can register MP ET to the MN mesh (original CT) following one of the two paths.

Fig. 4 .Fig. 5 .

 45 Fig. 4. (a) Original PET, (d) CT images in a normal case. Correspondences between selected points in the PET image and in the CT image are shown in (b) for the direct method and (e) for the method with the breathing model (corresponding points are linked). The registration result is shown in (c) for the direct method and in (f) for the method with the breathing model. Details of registration on the bottom part of right lung, (g) CT, (h) PET registered without breathing model, (c) with breathing model. The white crosses correspond to the same coordinates.

  have applied our algorithm on a normal case and on four pathological cases, exhibiting one tumor. In all cases, we have one PET (of size 144 × 144 × 230 with resolution of 4×4×4 mm 3 or 168×168×329 with resolution of 4×4×3 mm 3 ) and two CT volumes (of size 256 × 256 × 55 with resolution of 1.42 × 1.42 × 5 mm 3 to 512 × 512 × 138 with resolution of 0.98 × 0.98 × 5 mm 3 ), acquired during breathhold in maximum inspiration and in intermediate inspiration, from individual scanners.


Louis Gouillart
email: louis.gouillart@c2n.upsaclay.fr

Andrea Cattoni

Julie Goffard

Frederique Donsanti

Gilles Patriarche

Marie Jubault

Negar Naghavi

Stéphane Collin

Development of reflective back contacts for high-efficiency ultrathin Cu(In,Ga)Se2 solar cells

Keywords: Solar cells, Copper indium gallium selenide, Ultrathin films, Transparent back contact, Reflective back contact, Nanostructured back mirror

Because of poor light absorption, Cu(In,Ga)Se2-based (CIGS) solar cells with an ultrathin absorber layer (<500 nm) require the development of reflective back contacts. To enhance rear reflectance in CIGS ultrathin devices, we investigate novel back contact architectures based on a silver metallic mirror covered with a thin layer of In2O3:Sn (ITO), which is fully compatible with nanopatterning for further light trapping improvements. First, numerical electromagnetic simulations of complete solar cells have been performed for a 490 nm thick CIGS absorber with various back contacts. We predict a short-circuit current density of JSC = 34.0 mA/cm 2 for a 490 nm thick CIGS absorber with a silver nanostructured mirror. Second, we have fabricated and characterized 490 nm thick CIGS solar cells with transparent back contacts made of ITO, and reflective back contacts made of silver covered with ITO. Solar cells with a transparent ITO back contact exhibit an average efficiency of 10.0 %, compared to 9.3 % for standard molybdenum back contacts. A 5 nm thick Ga2O3 layer is revealed at the ITO/CIGS interface by transmission electron microscopy and energy dispersive X-ray spectroscopy. When silver is added, the reflective back mirror leads to a JSC improvement of 4.6 mA/cm 2 (from 22.4 to 27.0 mA/cm 2 ). These results pave the way for efficient ultrathin CIGS solar cells on reflective back contacts.

Introduction

Cu(In,Ga)Se2-based (CIGS) solar cells are one of the most promising thin-film technologies, with a record efficiency of 22.9 % achieved with a 2-3 µm thick absorber [1], [START_REF] Jackson | Effects of heavy alkali elements in Cu(In,Ga)Se2 solar cells with efficiencies up to 22.6%[END_REF]. However, the scarcity and high cost of Indium are a drawback for industrial production of competitive modules. With a CIGS absorber thickness lower than 500 nm it is possible to reduce deposition time and materials consumption, resulting in a decreased manufacturing cost [START_REF] Edoff | Technological and economical aspects on the influence of reduced Cu(In,Ga)Se2 thickness and Ga grading for coevaporated Cu(In,Ga)Se2 modules[END_REF], [START_REF] Horowitz | Cost and potential of monolithic CIGS photovoltaic modules[END_REF].

Ultrathin CIGS-based solar cells exhibit lower short-circuit current densities (JSC) and efficiencies mainly due to lower light absorption, enhanced recombination [START_REF] Vermang | Introduction of Si PERC Rear Contacting Design to Boost Efficiency of Cu(In,Ga)Se Solar Cells[END_REF] and low reflectivity of the conventional Mo back contact [START_REF] Jehl | Thinning of CIGS solar cells: Part II: Cell characterizations[END_REF], [START_REF] Mollica | Light absorption enhancement in ultra-thin Cu(In,Ga)Se2 solar cells by substituting the back-contact with a transparent conducting oxide based reflector[END_REF]. These loss mechanisms can be overcome by substituting the Mo back contact with passivating and more reflective back contact. Up to now, there are few studies of reflective back contacts for CIGS solar cells, such as: direct use of reflective or metallic back contacts [START_REF] Orgassa | Alternative back contact materials for thin film Cu(In,Ga)Se2 solar cells[END_REF]- [START_REF] Schleussner | Reactively sputtered ZrN for application as reflecting back contact in Cu(In,Ga)Se2 solar cells[END_REF], a combination of transparent conducting oxide and metallic reflector [START_REF] Mollica | Light absorption enhancement in ultra-thin Cu(In,Ga)Se2 solar cells by substituting the back-contact with a transparent conducting oxide based reflector[END_REF], passivating reflective back contact [START_REF] Vermang | Highly reflective rear surface passivation design for ultra-thin Cu(In,Ga)Se2 solar cells[END_REF], or nanostructured back contacts [START_REF] Goffard | Light Trapping in Ultrathin CIGS Solar Cells with Nanostructured Back Mirrors[END_REF], [START_REF] Yin | Optoelectronic Enhancement of Ultrathin CuIn1-xGaxSe2 Solar Cells by Nanophotonic Contacts[END_REF]. Nanostructured back reflectors are a promising strategy for enhanced absorption in ultrathin CIGS layers, as they lead to calculated JSC up to 36.3 mA/cm 2 in a 150 nm thick CIGS layer compared to 23.5 mA/cm 2 with a standard Mo back contact [START_REF] Goffard | Light Trapping in Ultrathin CIGS Solar Cells with Nanostructured Back Mirrors[END_REF]. Nanostructured mirrors have also yielded significant JSC improvements in other types of solar cells such as amorphous Si:H [START_REF] Ferry | Optimized Spatial Correlations for Broadband Light Trapping Nanopatterns in High Efficiency Ultrathin Film a-Si:H Solar Cells[END_REF] and ultrathin GaAs [START_REF] Yang | Ultra-thin GaAs single-junction solar cells integrated with a reflective back scattering layer[END_REF], [START_REF] Chen | 200nm-thick GaAs solar cells with a nanostructured silver mirror[END_REF].

In this paper we present the first steps towards the development of a reflective back contact based on a silver metallic mirror encapsulated with a thin layer of In2O3:Sn (ITO). Reflective back contacts for CIGS solar cells need not only to provide enhanced rear reflectance but also to form an ohmic contact with the absorber. Several studies report on the formation of a resistive interfacial layer of Ga2O3 after CIGS deposition on oxide layers such as hydrogendoped In2O3, ITO, MoO3, SnO2:F, ZnO or ZnO:Al [START_REF] Keller | Bifacial Cu(In,Ga)Se2 solar cells using hydrogen-doped In2O3 films as a transparent back contact[END_REF]- [START_REF] Nakada | Microstructural and diffusion properties of CIGS thin film solar cells fabricated using transparent conducting oxide back contacts[END_REF]. In particular, ITO is known to prevent Na diffusion from the glass to the CIGS absorber and to form a detrimental layer of Ga2O3 for CIGS deposition temperatures above 520 °C [START_REF] Nakada | Microstructural and diffusion properties of CIGS thin film solar cells fabricated using transparent conducting oxide back contacts[END_REF], [START_REF] Nakada | Novel device structure for Cu(In,Ga)Se2 thin film solar cells using transparent conducting oxide back and front contacts[END_REF]. As a result, the stability and electrical properties of ultrathin devices on transparent ITO-based substrates were first investigated by co-evaporation of CIGS at low temperature (450°C) with incorporation of Na from a NaF post-deposition treatment.

In this study, the improvement of light absorption in a 490 nm thick CIGS absorber is first simulated using a nanostructured and a flat mirror encapsulated in transparent conducting oxides such as ITO. were extracted from ellipsometric measurements on a UVISEL ellipsometer from Horiba Jobin-Yvon. The ITO resistivity was measured with a four-point probe system before (ρ = 9.0x10 -4 Ω.cm) and after annealing (ρ = 7.0x10 -4 Ω.cm). Finally, a Sentech reflectometer was used to determine the spectral reflectance of four different substrates: SLG/600-nm-Mo, SLG/150-nm-Ag, SLG/150-nm-Ag/30-nm-ITO and SLG/150-nm-Ag/200-nm-ITO substrates.

Current-voltage I(V) measurements of 0.1 cm 2 solar cells were carried out with a Newport class AAA solar simulator under AM1.5G illumination and in the dark. Devices on transparent substrates were placed on a black paper before light I(V) measurements to avoid additional rear-reflection from the substrate holder. External Quantum Efficiency (EQE) measurements were performed with a IQE200 Newport Instrument. A slab of a 490 nm thick CIGS layer deposited on an ITO-based transparent back contact was prepared by a focused ion beam, and analysed with TEM in a FEI Titan Themis XFEG instrument at an acceleration voltage of 200 kV. Composition mappings were determined by EDX.

Results and discussion

Absorption simulation

We first investigate numerically the optical effects of a nanostructured back mirror, a flat back mirror, and a transparent back contact as compared to a conventional Mo back contact. The spectral absorption is simulated with a rigorous coupled wave analysis method, and the optical indices are calculated from ellipsometric measurements for CIGS and annealed ITO, and taken from the literature for other materials. More details can be found in reference [START_REF] Goffard | Light Trapping in Ultrathin CIGS Solar Cells with Nanostructured Back Mirrors[END_REF].

The architecture of flat ultrathin devices is sketched in Fig. 1a), while the optimized geometry of a complete device including a nanostructured back mirror is presented in Fig. 1b) and c).

This simulated nanostructured back mirror consists of square-shaped dielectric nanogrid on top of a flat ITO-coated Ag stack. The dielectric nanogrid has a refractive index of n = 1.5, a height of 200 nm, a period of 600 nm and hole width of 350 nm. The experimental fabrication of such a dielectric nanogrid is under development, using a scalable process based on direct nanoimprint of a sol-gel derived metal oxide [START_REF] Bottein | Environment-controlled sol-gel soft-NIL processing for optimized titania, alumina, silica and yttria-zirconia imprinting at sub-micron dimensions[END_REF], [START_REF] Dalstein | Nanoimprinted, Submicrometric, MOF-Based 2D Photonic Structures: Toward Easy Selective Vapors Sensing by a Smartphone Camera[END_REF]. The simulated absorptions of CIGS devices on Mo, ITO, as well as flat and nanostructured reflective back contacts are plotted in Fig. 2.

As shown in Fig. 2a, using a standard Mo back contact leads to substantial absorption losses for long wavelengths (> 600 nm), because of the low reflectivity at CIGS/Mo interface.

Moreover, optical losses are also observed at short wavelengths (< 500 nm) due to the low bandgap of CdS (2.4 eV) buffer layer. Assuming a perfect collection of photogenerated carriers, a theoretical short-circuit current density JSC can be deduced from the simulated CIGS absorption [START_REF] Goffard | Light Trapping in Ultrathin CIGS Solar Cells with Nanostructured Back Mirrors[END_REF]. It is limited to 28.8 mA/cm 2 with a Mo back contact. When Mo is replaced by ITO (Fig. 2b) a similar spectral absorption with a JSC of 28.5 mA/cm 2 is calculated. The introduction of an Ag layer results in enhanced back reflectance and significant improvement of CIGS light absorption at long wavelengths (> 800 nm), (Fig. 2c), which leads to an increased JSC of 32.2 mA/cm 2 . Finally, if the flat mirror is replaced by an optimised nanostructured back mirror made of a nanogrid, additional absorption resonances occur close to the bandgap (Fig. 2d). Multi-resonant absorption in the thin CIGS layer leads to a JSC of 34.0 mA/cm 2 , which represents an additional 1.8 mA/cm 2 gain as compared to the flat back mirror.

Optical characterization of back contacts

Before deposition of CIGS on transparent and reflective back contacts, the spectral reflectance of 4 different back contacts deposited on SLG were first measured and compared: a 600 nm thick Mo, a 150 nm thick Ag layer, and 150 nm thick Ag layers with 30 nm and 200 nm ITO coatings (Fig. 3a). As expected, Mo has a low reflectivity and cannot act as a back mirror, resulting in low CIGS absorption for ultrathin CIGS solar cells (Fig. 2a). On the contrary, the silver layer provides a high reflectivity (R>90 % for λ>500 nm), but it cannot be used as a back contact due to diffusion of Ag in CIGS during the growth. When the silver layer is covered with ITO coatings, the reflectivity is significantly reduced due to parasitic absorption.

Reflectivity dips are attributed to Fabry-Perot resonances in the ITO layer. After a 10-minuteannealing in air at 540 °C, the reflectivity is strongly enhanced and exceeds 90 % for wavelengths above 550 nm. This annealing experiment performed in air proves that the ITO layer efficiently encapsulates the Ag mirror, preventing its oxidation and morphological changes that negatively impact its reflectivity; the stability of such reflective back contact under CIGS deposition was only tested after the complete solar cell fabrication. To determine the origin of the enhanced reflectivity of Ag/ITO stacks after annealing, ellipsometric measurements were carried out on a 200 nm thick ITO layer sputtered on SLG before and after annealing in air at 540 °C during 10 minutes. The refractive index n and extinction coefficient k are displayed in Fig. 3b). The annealing results in an increase of n and a decrease of k in the visible and infrared domains (λ>500 nm). The origin of these effects was not investigated in detail, though it is commonly accepted that annealing of amorphous ITO generally results in an increased layer density, an improved carrier mobility and a reduced free-carrier absorption. It leads to a significant improvement of the Ag/ITO reflectivity and a spectral shift of the reflectivity dips toward shorter wavelengths (Fig. 3a).

TEM-EDX back contact characterization

A TEM-EDX analysis was performed in order to investigate the ITO/CIGS interface. Fig. 4a) and b) show a TEM image of an ultrathin 490-nm-CIGS layer deposited on an ITO back contact along with the average composition profile of absorber elements deduced from EDX

analysis. An average GGI of 0.28 was calculated from the average EDX signal of the CIGS layer, which is close to the GGI value of 0.32 determined by XRF. It is worth mentioning that the EDX signal of Cu is overestimated as Cu from the substrate holder is also detected. Hence CGI cannot be deduced from this analysis.

This EDX analysis confirms that the CIGS composition is homogeneous from the front to the back interface, as expected from a one-stage CIGS deposition process. Moreover it reveals that Ga segregates at the back contact.

Closer views are shown in Fig. 4. with (c) an additional dark field TEM image and its corresponding compositional EDX mappings of O, Ga, In, and Se elements. They reveal a segregation of both Ga and O elements at the ITO/CIGS interface, which suggests that a thin Ga2O3 layer is formed during CIGS deposition. This layer should be highly resistive and possibly n-doped [START_REF] Keller | Bifacial Cu(In,Ga)Se2 solar cells using hydrogen-doped In2O3 films as a transparent back contact[END_REF]. However, in this study the low substrate temperature (450 °C) during CIGS deposition results in a very thin Ga2O3 layer of only 5 nm that does not seem to deteriorate solar cell performances.

Solar cell performances

The performances of ultrathin CIGS solar cells are summarized in Table 1 for Mo and ITO According to simulation results (Fig. 2), the short-circuit current could also be increased by adding a flat silver back mirror between the ITO layer and the glass substrate, with an expected gain of 4.6 mA/cm 2 over the Mo reference. We have performed preliminary experiments of 490 nm thick CIGS solar cells deposited on reflective multilayer stacks (MLS) based on silver covered by ITO. They led to poor electrical performances due to process issues that should be overcome by further optimizations. However, the reflective MLS resulted in a strong EQE enhancement. The EQE of ultrathin solar cells with a Mo back contact, an ITO back contact, and a reflective MLS are plotted and compared in Fig. 5b. The reflective back contact leads to an improvement for wavelengths above 500 nm and results in a promising JSC = 27.0 mA/cm 2 . The same EQE measurements are also plotted and compared to numerical simulation in Fig. 2. A very good agreement between experiments and simulations is obtained.

Conclusion

Ultrathin

List of figures

EQE

Wavelength (nm)

Mo ITO MLS

Fig. 5 ,

 5 Fig. 5, together with the experimental EQE of a 490 nm thick CIGS cell on a reflective back contact (Ag). The reference solar cells on a SLG/Al2O3/Mo substrate exhibit an average efficiency of 9.3 %.

Fig. 1

 1 Fig. 1 Sketch of an ultrathin CIGS solar cell a) on planar back contacts such as molybdenum, TCO, or reflective multi-layer stack (MLS) and b) on a nanostructured back mirror consisting of a periodically patterned dielectric on top of reflective MLS. c) Three dimensional view of the ultrathin device on a nanostructured back mirror with optimized geometry. CIGS is deposited on top of a square-shaped dielectric nanogrid with a height of 200 nm, a period of 600 nm and hole width of 350 nm.

Fig. 2

 2 Fig. 2 Simulated absorption within each layer of 490 nm thick CIGS solar cells under AM1.5G illumination, with back contacts made of a) molybdenum, b) ITO, c) reflective multi-layer stack (MLS) and d) reflective MLS with nanopatterns. Respective experimental EQEs are given in a), b) and c) for comparison.

Fig. 3 .

 3 Fig. 3.

Fig. 4 a

 4 Fig. 4 a) TEM dark field image of ultrathin CIGS layer deposited on ITO back contact with b) corresponding composition profile from average EDX signal, c) TEM dark field image of ultrathin CIGS layer close to the ITO back contact with associated EDX mappings of O, Ga, In, and Se elements.

Fig. 5

 5 Fig. 5 Experimental I(V) and External Quantum Efficiency (EQE) curves of 490 nm thick CIGS solar cells on molybdenum (black), ITO (red) and reflective multilayer stack (blue) back contacts.

  CIGS solar cells with transparent and reflective back contacts were first numerically investigated. A JSC of 32.2 mA/cm 2 is calculated for 490 nm thick CIGS solar cells with a flat back mirror, corresponding to a 3.4 mA/cm 2 gain as compared to standard molybdenum back contacts. An additional optical gain is predicted with a nanostructured reflective back contact, leading to a JSC of 34.0 mA/cm 2 . Flat reflective stacks of Ag/ITO were fabricated and are stable under a 10-minute-annealing in air at 540 °C. Their spectral reflectance is over 90 % for wavelengths above 550 nm, which should substantially increase absorption in ultrathin CIGS absorbers. Complete ultrathin CIGS solar cells were fabricated on molybdenum and

	transparent ITO back contacts, with respective average efficiencies of 9.3 % and 10.0 %. ITO
	forms an ohmic back contact with CIGS, as confirmed by an average FF of 70.3 %. A TEM-
	EDX investigation shows that the thickness of the Ga2O3 formed at ITO/CIGS interface
	during CIGS deposition is limited to 5 nm and has no impact on solar cell performances.
	Preliminary experiments on reflective back contacts made of Ag and ITO led to a JSC = 27.0
	mA/cm

2

. These results provide a promising route toward efficient, ultrathin CIGS solar cells deposited on reflective back contacts.

Table captions Table 1 .

 captions1 Summary of I(V) parameters for 490 nm thick CIGS solar cells with various back contacts. A two-diode model is used to fit dark I(V) curves (J01 and J02: saturation currents for respective ideality factors of 1 and 2, RSH: shunt resistance, RS: series resistance).

Acknowledgements

This work is supported by the ARCIGS-M project within the European Union's Horizon 2020 research and innovation program under grant agreement No. 720887.


Hong Cai

Loïc Quéval

Martin Hennebel

Cyril Gisbert

Marie Petitet

" Long

Long-term Unit Commitment Problem with Optimal Zone Configuration: a Case Study in Martinique

I. Introduction

Island regions are not interconnected with the mainland and the isolation of these regions has brought challenges to the electric power systems. For instance, the power grids in island regions are more fragile to the rapid variance of the electricity consumption or production than the continental grid. Meanwhile, more intermittent renewable energy, such as wind and solar, has been connected to the power grid. Some of the overseas territories of France are island territories in the Atlantic, Pacific and Indian Oceans.

In this work, we propose a long-term unit commitment (UC) model with simplified network (zone) constraints to apply in these areas.

One of the main purposes of this work is to examine whether a better zonal configuration could help to increase the performance of the long-term UC model. Also, by applying this long-term UC model, we can create a robust generation plan for these territories and would be able to test, for instance, whether a larger amount of renewable energy could be introduced in these areas and whether the performance of the system could be improving by constructing a new transmission line.

II.

Method

The method is made of four steps (Fig. 1). The first step is to simulate short-term market outcomes by using an economic dispatch model with a detailed nodal network. We consider different levels of load and renewable generation. This nodal network model considers both the network and economic constraints, therefore, the solution provides price signals to optimally utilize the available resources. The second step is to find the optimal zonal configuration. The price signals given by the first step are used to decide the zonal configurations. The similarity of prices for interconnected nodes implies a lower degree of congestion.

We cluster nodes with similar prices into one zone. The third step is to decide the simplified network constraints for the long-term unit commitment. There are two different sets of simplified network constraints which could be applied in the longterm unit commitment problem, i.e., net transmission capacity (NTC) and flowbased market coupling (FBMC). In this study, we use the NTC network constraints for our study. The simplified network constraints are mainly to limit the power exchange between two zones and reduce the network congestion in real time. The fourth step is to apply the simplified network constraint in the long-term unit commitment problem. The nodal economic dispatch model uses an optimal power flow approach, which takes into account the physical and technical constraints [START_REF] Schweppe | Spot pricing of electricity[END_REF], to give the price of power for each node. Most power systems are alternating current (AC) and the problem that we ideally would like to solve is an alternating current optimal power flow (AC-OPF) problem. This is a difficult problem to solve, due to nonlinearities and non-convexities. Thus, in practice, most of system operators us a lossless direct current (DC) approximation is used. The resulting problem is called direct current optimal power flow (DC-OPF). This approximation has been proved to efficiently approximate the real network [START_REF] Overbye | A comparison of the AC and DC power flow models for LMP calculations[END_REF]. The nodal economic dispatch model is an optimization problem given by,

𝑀𝑀𝑀𝑀𝑀𝑀 𝐺𝐺 ∑ 𝑐𝑐𝑐𝑐 𝑐𝑐 × 𝐺𝐺 𝑐𝑐 𝑐𝑐 (1)
Subject to:

𝑁𝑁𝑁𝑁 𝑛𝑛 = ∑ 𝐺𝐺 𝑐𝑐, + 𝑐𝑐𝑔𝑔𝑀𝑀𝑀𝑀𝑔𝑔 𝑛𝑛 + 𝑐𝑐𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔 𝑛𝑛 + 𝑐𝑐𝑔𝑔𝑀𝑀𝑔𝑔 𝑛𝑛 -𝑄𝑄 𝑛𝑛 , ∀𝑀𝑀 𝑐𝑐∈𝐶𝐶 𝑛𝑛 (2)
𝐺𝐺 𝑐𝑐 ≤ 𝑐𝑐𝑔𝑔𝑔𝑔𝑔𝑔 𝑐𝑐 , ∀𝑐𝑐 ∈ 𝐶𝐶 (3)

𝐿𝐿𝐿𝐿 𝑛𝑛,𝑛𝑛𝑛𝑛 = 𝑔𝑔𝑏𝑏𝑏𝑏𝑐𝑐𝑏𝑏𝑔𝑔𝑔𝑔 𝑛𝑛,𝑛𝑛𝑛𝑛 (∆ 𝑛𝑛 -∆ 𝑛𝑛𝑛𝑛 ), ∀ 𝑀𝑀, 𝑀𝑀𝑀𝑀 (4)
𝑁𝑁𝑁𝑁 𝑛𝑛 = ∑ 𝐿𝐿𝐿𝐿 𝑛𝑛,𝑛𝑛𝑛𝑛 𝑛𝑛𝑛𝑛 -∑ 𝐿𝐿𝐿𝐿 𝑛𝑛𝑛𝑛,𝑛𝑛 𝑛𝑛𝑛𝑛 , ∀𝑀𝑀 (5)
-𝑝𝑝𝑔𝑔𝑔𝑔𝑔𝑔 𝑛𝑛,𝑛𝑛𝑛𝑛 × 𝑆𝑆𝑏𝑏 𝑛𝑛,𝑛𝑛𝑛𝑛 ≤ 𝐿𝐿𝐿𝐿 𝑛𝑛,𝑛𝑛𝑛𝑛 ≤ 𝑝𝑝𝑔𝑔𝑔𝑔𝑔𝑔 𝑛𝑛,𝑛𝑛𝑛𝑛 × 𝑆𝑆𝑏𝑏 𝑛𝑛,𝑛𝑛𝑛𝑛 , ∀𝑀𝑀, 𝑀𝑀𝑀𝑀

∆ 𝑛𝑛 ′ = 0, 𝑀𝑀 ′ = 𝑔𝑔𝑔𝑔𝑔𝑔𝑐𝑐𝑠𝑠𝑔𝑔𝑠𝑠𝑔𝑔 (6)
𝐺𝐺 𝑐𝑐 ≥ 0, 𝑄𝑄 𝑛𝑛 ≥ 0 (7)
Where 𝐺𝐺 𝑐𝑐 is the generation of power plant c, 𝑁𝑁𝑁𝑁 𝑛𝑛 is net injection of node n, and 𝑆𝑆𝑏𝑏 𝑛𝑛,𝑛𝑛𝑛𝑛 is the binary variable, 0 indicating the transmission line on an off statue and 1 indicating an on statue.

The objective of the nodal pricing model is to minimize the operation cost (Eq.

(1)), considering the various constraints. The areas applying nodal pricing are constrained by the energy balance (Eq. ( 2)), the maximum generation capacity of thermal power plants (Eq. ( 3)), and the restrictions on power transmission (Eq. ( 4)-( 6)). The energy balance (Eq. ( 2)) ensures that at node 𝑀𝑀 , net input or withdrawal 𝑁𝑁𝑁𝑁 𝑛𝑛 is equal to the difference between power generation (including all the conventional power generation 𝐺𝐺 𝑐𝑐 at node n, wind generation 𝑐𝑐𝑔𝑔𝑀𝑀𝑀𝑀𝑔𝑔 𝑛𝑛 , solar generation 𝑐𝑐𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔 𝑛𝑛 and biomass generation 𝑐𝑐𝑔𝑔𝑀𝑀𝑔𝑔 𝑛𝑛 ) and nodal demand 𝑄𝑄 𝑛𝑛 .

The generation of renewable is given exogenously. The dual variables of Eq. ( 2), which are the marginal costs/benefits of increasing injections in the nodes by one unit, are the nodal prices. Thermal power generation is restricted by the maximum generation requirement gmax c (Eq. ( 3)). The DC approximation (see [START_REF] Breuer | Determination of alternative bidding areas based on a full nodal pricing approach[END_REF] and [START_REF] Sauma | Proactive planning and valuation of transmission investments in restructured electricity markets[END_REF]) is used to determine the load flows 𝐿𝐿𝐿𝐿 𝑛𝑛,𝑛𝑛𝑛𝑛 (Eq. ( 4)) in each line and the resulting injection or withdrawal 𝑁𝑁𝑁𝑁 𝑛𝑛 (Eq. ( 5)). Eq. ( 6) is to ensure that the load flow over a line should not exceed its thermal capacity limit. 𝑆𝑆𝑏𝑏 𝑛𝑛,𝑛𝑛𝑛𝑛 describes the status of a line (on or off). It is used to enforce the N-1 security constraint. Under this constraint, the failure of a single power line should not cause power outages.

Enforcing the N-1 constraint improves security, however, might at a cost. Eq. ( 7)

is to specify the slack bus. Eq. ( 8) requires the generation and the demand to be positive.

Step 2: Zonal clustering

The zonal clustering algorithm groups nodes into a pre-defined number of zones.

Generally, it is required that there is not much intra-zonal congestion and that nodes within a zone should have a similar impact on inter-zonal lines. Different clustering algorithms have been proposed in the literature. The two most used are the nodal pricing method [START_REF] Burstedde | From nodal to zonal pricing: a bottom-up approach to the second-best[END_REF] - [START_REF] Breuer | Determination of alternative bidding areas based on a full nodal pricing approach[END_REF] and nodal Power Transfer Distribution Factors (nodal PTDFs) method [START_REF] Klos | The scheme of a novel methodology for zonal division based on power transfer distribution factors[END_REF].

The capacity limitation of the transmission network prevents full uses the cheapest energy in the network, i.e., there is no enough capacity to transfer the power from the cheapest places to other places. The price at each node is also decided by the network topology and energy location. Nodal prices differences indicate the level of grid congestion. Less prices differences between nodes mean less congestion.

Therefore, the clustering of nodes with similar nodal prices will result in a low intra-zonal congestion zone configuration.

The nodal PTDFs quantifies how the power injection at a given node impacts a given line. The nodal PTDF-values based method takes the congested lines as inter-zonal line, and clusters nodes with similar PTDF-values for the congested lines within the same zone. In this way, a zonal configuration is obtained with low intra-zonal congestion and nodes within the same zone have a similar impact on the interzonal links.

In this work, we apply the nodal pricing method to group the nodes into a zone.

However, we also consider the bottlenecks (i.e., the most congested lines) in the systems. That is, we identify the bottlenecks and set them as the inter-zonal lines before we group the nodes with similar prices. The zonal clustering is an optimization problem given by, Minimize:

𝑀𝑀𝑀𝑀𝑀𝑀 ∑ ∑ (𝑝𝑝 𝑛𝑛 𝑔𝑔 𝑛𝑛,𝑧𝑧 -𝑐𝑐 𝑘𝑘 𝑔𝑔 𝑛𝑛,𝑧𝑧 ) 2 z n (9)
Subject to:

∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 𝑧𝑧 = 1∀𝑀𝑀 (10)
∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 𝑛𝑛 ≥ 𝑍𝑍𝑁𝑁 𝑧𝑧 ∀𝑧𝑧 (11)
𝑐𝑐 𝑧𝑧 ∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 𝑛𝑛 = ∑ (𝑝𝑝 𝑛𝑛 𝑛𝑛 𝑔𝑔 𝑛𝑛,𝑧𝑧 ) ∀𝑧𝑧 (12)
𝑔𝑔𝑔𝑔 𝑛𝑛,𝑧𝑧 (∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 -1

𝑛𝑛 ) ≤ 10000𝑔𝑔 𝑛𝑛,𝑧𝑧 (∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 -1 𝑛𝑛 ) ∀𝑀𝑀, 𝑧𝑧 (13)
𝑔𝑔𝑔𝑔 𝑛𝑛,𝑧𝑧 �∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 -

1 𝑛𝑛 � ≥ �1 -10000(1 -𝑔𝑔 𝑛𝑛,𝑧𝑧 �](∑ 𝑔𝑔 𝑛𝑛,𝑧𝑧 -1 𝑛𝑛 ) ∀𝑀𝑀, 𝑧𝑧 (14)
𝑔𝑔𝑔𝑔 𝑛𝑛,𝑧𝑧 = ∑ 𝑔𝑔𝑀𝑀𝑀𝑀𝑠𝑠 𝑛𝑛,𝑛𝑛𝑛𝑛 × 𝑔𝑔 𝑛𝑛𝑛𝑛,𝑧𝑧 × 𝑔𝑔 𝑛𝑛,𝑧𝑧, 𝑛𝑛𝑛𝑛 ∀𝑀𝑀 ≠ 𝑀𝑀𝑀𝑀

Where 𝑔𝑔 𝑛𝑛,𝑧𝑧 is a binary variable. 𝑔𝑔 𝑛𝑛,𝑧𝑧 = 1 indicates that the node n is within zone z while 𝑔𝑔 𝑛𝑛,𝑧𝑧 = 0 indicates node n is not in zone z. 𝑔𝑔𝑔𝑔 𝑛𝑛,𝑧𝑧 is the total number of connections with the other nodes within the same zone z for node n. 𝑐𝑐 𝑧𝑧 is the average price for zone z.

The objective of the clustering technique is to minimize the price differences within the same zone (Eq. ( 9)). Eq. ( 10) is to force each node to only belong to one zone. Eq. ( 11) is to enforce that each zone has at least 𝑍𝑍𝑁𝑁 𝑧𝑧 node (e.g., 1 or 2 nodes). Eq. ( 12) defines the average price 𝑐𝑐 𝑘𝑘 in zone z. Eq. ( 13) and ( 14) are to enforce that any node should have at least one connection with another node within the same zone if the number of nodes in a zone is greater than 1. Eq. ( 15)

is to calculate the total number of connections with the other nodes within the same zone, where 𝑔𝑔𝑀𝑀𝑀𝑀𝑠𝑠 𝑛𝑛,𝑛𝑛𝑛𝑛 is the connecting matrix, which tells if there is a direct connecting between two nodes n and nn.

Step 3: Optimal zonal parameters

This section gives the procedure for implementing the NTC constraints. They are given by,

NEX 𝑧𝑧 = ∑ (𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧,𝑧𝑧𝑧𝑧 -𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧𝑧𝑧,𝑧𝑧 𝑧𝑧𝑧𝑧 ), ∀𝑧𝑧 ∈ ⌈1, 𝑍𝑍⌉ (16)
0 ≤ 𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧,𝑧𝑧𝑧𝑧 ≤ 𝑀𝑀𝑏𝑏𝑐𝑐 𝑧𝑧,𝑧𝑧𝑧𝑧 , ∀𝑧𝑧, 𝑧𝑧𝑧𝑧 (17)
Where NEX 𝑧𝑧 is the net position of a zone and 𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧,𝑧𝑧𝑧𝑧 is total export from zone z to zone zz.

NEX 𝑧𝑧 is defined as the difference of its total export 𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧,𝑧𝑧𝑧𝑧 and import 𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧𝑧𝑧,𝑧𝑧 (Eq. ( 20)). The NTC model does not have specific limitations on specified lines.

However, it restricts the total transfer 𝐵𝐵𝐵𝐵𝐵𝐵 𝑧𝑧,𝑧𝑧𝑧𝑧 between two pricing areas to a predefined maximum trading volume 𝑀𝑀𝑏𝑏𝑐𝑐 𝑧𝑧,𝑧𝑧𝑧𝑧 , (Eq. ( 21)).

Step 4: Long-term unit commitment (UC) model

For power system operation, one of the most important issue is to decide for each period which generating units should run to satisfy the demand. In a typical electrical system, there are a variety of resource and power plants available for generating electricity, and each of these units has its own characteristic. For instance, the start-up cost and the variable cost are different. Unit commitment problem is to help to make decisions in such a situation [21][24].

Unit commitment is one of the classic optimization problems in power systems.

The goal is to decide the ON/OFF status of all the generating units, which meet the forecasted loads and reserve requirements and provide a least-cost power generation plan. There are two main ways of defining this least-cost objective function. The first one is used by TSO, which minimize the total cost (i.e., minimizing fuel costs and startup/shutdown cost) or maximize the social welfare while meeting a forecasted hourly load. The second way is used by the generation companies, which maximize the total profit based on their bidding strategy. In this study, we focus on the first type.

Different mathematical formulations of the UC problem exist. One of the most used formulations of UC problems is the mixed integer linear programming (MILP) methods. Due to fact that there could be thousands of generation units and transmission lines existing in a system, the unit commitment problem modeled by MILP could be a computationally challenging problem with many integer variables and constraints. Different optimization techniques including Lagrangian relaxation and branch-and-bound based MILP methods have been used to solve large-scale UC problems [START_REF] Barth | A stochastic unit commitment model for the evaluation of the impacts of integration of large amounts of intermittent wind power[END_REF], [START_REF] Hobbs | The next generation of electric power unit commitment models[END_REF].

The physical conditions of the operating generators are the main constraints for the UC model, including the generation capacity limits, minimum ON time, minimum OFF time, ramp up/down rate, reserve requirements. In the following, we will illustrate the typical constraints in a UC model used in this project. More UC constraints could be found on [START_REF] Huang | Electrical power unit commitment: deterministic and two-stage stochastic programming models and algorithms[END_REF].

𝑔𝑔𝑀𝑀𝑀𝑀 ∑ ∑ (𝑔𝑔𝑏𝑏 𝑐𝑐 × 𝑉𝑉 𝑐𝑐𝑐𝑐 + 𝑔𝑔𝑔𝑔 𝑐𝑐 × 𝑊𝑊 𝑐𝑐𝑐𝑐 ) 𝑐𝑐 𝑐𝑐 + ∑ ∑ (𝑐𝑐𝑐𝑐 𝑐𝑐 × 𝐺𝐺 𝑐𝑐,𝑐𝑐 ) 𝑐𝑐 𝑐𝑐 + 𝑉𝑉𝑉𝑉𝐿𝐿𝐿𝐿 × ∑ ∑ 𝐿𝐿𝑆𝑆 𝑐𝑐,𝑐𝑐 𝑐𝑐 𝑐𝑐 (18) 𝑈𝑈 𝑐𝑐,𝑐𝑐 -𝑈𝑈 𝑐𝑐,𝑐𝑐-1 ≤ 𝑈𝑈 𝑐𝑐,𝜏𝜏 ∀c, t, τ = t, … , min{𝑏𝑏 + minion 𝑐𝑐 -1, |T|} (19)
𝑈𝑈 𝑐𝑐,𝑐𝑐-1 -𝑈𝑈 𝑐𝑐,𝑐𝑐 ≤ 1 -𝑈𝑈 𝑐𝑐,𝜏𝜏 ∀c, t, τ = t, … , min{𝑏𝑏 + minioff 𝑐𝑐 -1, |T|} (20)
𝑉𝑉 𝑐𝑐,𝑐𝑐 ≥ 𝑈𝑈 𝑐𝑐,𝑐𝑐 -𝑈𝑈 𝑐𝑐,𝑐𝑐-1 ∀𝑐𝑐, 𝑏𝑏

𝑊𝑊 𝑐𝑐,𝑐𝑐 ≥ -𝑈𝑈 𝑐𝑐,𝑐𝑐 + 𝑈𝑈 𝑐𝑐,𝑐𝑐-1 ∀𝑐𝑐, 𝑏𝑏

𝑐𝑐𝑔𝑔𝑔𝑔𝑔𝑔 𝑐𝑐 × 𝑈𝑈 𝑐𝑐,𝑐𝑐 ≤ 𝐺𝐺 𝑐𝑐,𝑐𝑐 ≤ 𝑐𝑐𝑔𝑔𝑀𝑀𝑀𝑀 𝑐𝑐 × 𝑈𝑈 𝑐𝑐,𝑐𝑐 , ∀𝑐𝑐, 𝑏𝑏 (22)
𝐺𝐺 𝑐𝑐 ≥ 0

-𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑀𝑀 𝑐𝑐 ≤ 𝐺𝐺 Based on two status of a generator 𝑈𝑈 𝑐𝑐 , we can know that when a generator is online and commits to supply capacity, i.e., 𝑈𝑈 𝑐𝑐 = 1 and when the generator commitment status is in an "off" state, i.e., 𝑈𝑈 𝑐𝑐 = 0, the power dispatch level become zero and has no any operational cost.

Some specific operations are also enforce in the UC constraints, such as the minimum ON time, i.e.,minion 𝑐𝑐 and minimum OFF time, i.e., minioff 𝑐𝑐 , and also specify startup action and shutdown action on each unit at a time period t, respectively. A generator can't be started up or shut down arbitrarily in consecutive hours, Eq [START_REF] Schweppe | Spot pricing of electricity[END_REF] and ( 24) respectively indicate two generator's requirements: the shortest ON duration has to be met before a generator being shut down and the shortest OFF duration is also required before a generator being restarted up.

The startup action 𝑉𝑉 𝑐𝑐,𝑐𝑐 and the shutdown action 𝑊𝑊 𝑐𝑐,𝑐𝑐 are determined by the generator commitment statuses in the previous time period t -1 and the current time period t, as in Eq. ( 25) and [START_REF] Tuohy | Unit commitment for systems with significant wind penetration[END_REF]. Any operational actions can incur start-up or shut-down costs, which are considered in the objective function.

A generator output in an hour is subject to the maximum generation limit 𝑐𝑐𝑔𝑔𝑔𝑔𝑔𝑔 𝑐𝑐 , and the minimum generation limit 𝑐𝑐𝑔𝑔𝑀𝑀𝑀𝑀 𝑐𝑐 , . When a generator is scheduled online (𝑈𝑈 𝑐𝑐,𝑐𝑐 = 1), the generation capacity is active giving bounds on dispatch level, shown in Eq. ( 27); otherwise, a generator output is forced to zero.

The generation must be positive Eq. ( 28).

In addition, a generator output can be adjusted, increasing 𝑔𝑔𝑠𝑠𝑝𝑝 𝑐𝑐 or decreasing 𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑔𝑀𝑀 𝑐𝑐 between two successive time periods. The generation difference between two adjacent time periods is called ramping. A basic constraint to address generation ramping is presented in Eq [START_REF] Yang | Subhourly unit commitment with feasible energy delivery constraints[END_REF].

Operating reserve is one type of ancillary operations to support the power balance on the demand sides. The current operating reserve services being offered in electric energy markets include synchronous or non-synchronous, regulation reserves, spinning reserves, and non-spinning reserves. Generally, the sources of energy provided from different reserve services are different: regulation service mainly supplied from online generators, partial spinning reserve provided from generators already connected to the grid or system resources, and non-spinning reserve provided from quick-start generators, system resources or interruptible loads. The response times of reserve services can vary greatly, depending on the control reserve deployment time. To achieve the optimization of energy and reserve in practice, one can obtain an efficient energy and reserve offering strategy by Heuristic method [START_REF] Ni | Optimal integrated generation bidding and scheduling with risk management under a deregulated power market[END_REF] or consider the reserve determination on precontingency and post-contingency conditions [START_REF] Bouffard | Market-clearing with stochastic securitypart I: formulation[END_REF].

In our models, if we consider the spinning reserve problem, the spinning reserve 𝑆𝑆 𝑐𝑐,𝑐𝑐 is generally accounted for partial online generating capacity or off-line generation resources. Their outputs are constrained by predetermined maximum spin reserve 𝑔𝑔 𝑐𝑐 𝑚𝑚𝑚𝑚𝑚𝑚 , shown as Eq. (30).

Meanwhile, the generators that participate in biding spinning reserve must meet the spin reserve requirements given by TSOs. Eq (31) describes an operating condition that the total spinning reserve at bus i should not less than the fixed reserve requirement 𝑔𝑔𝑔𝑔 𝑖𝑖,𝑐𝑐 .

III. Case study 1. Presentation

We apply the method described in the section II to the power network of Martinique. Martinique is an island in the Caribbean Sea and is an overseas department of France. Martinique has a tropical and humid climate with an average temperature of 23 degrees to 29 degrees Celsius. The average monthly temperatures remain remarkably stable, varying by only a few degrees year-round.

Therefore, the power demand is also very stable during the year.

In Martinique, the main renewable energy resource is solar energy (about 15 %).

They also have little share of biomass and wind energy (less than 5 %) [31].

Currently, the system operator limits that the maximum share of intermittent renewable energies that can be injected into the network at a given time to 35%

[31]. Beyond this threshold, disconnections are possible in order to preserve the stability of the electrical system. However, the system operator is considering relaxing this restriction. Figure 1 shows a typical generation curve in Martinique.

In a normal day, there are two types of peak hours. The first peak hours are in the daytime (10 h-14 h) with solar power. The second peak hours are in the nighttime 18 h-20 h) without solar power. As congestion happens mostly during the peak hours, we mostly focus in these periods in this work, especially in the first twosteps model. [source: opendata-martinique.edf.fr]

Results of step 1

We only consider the HTB grid (63 kV) in Martinique, as shown in Figure 3. In total, there are 14 hubs (nodes) and 25 HTB transmission lines. The thermal power plants are located at the nodes "Bellefontaine", "Hydrobase" and "Galion".

We choose four typical hours in the working days for our first two steps analysis as show in Note that the N-1 constraint might affect the optimal dispatch. For instance, different prices and quantities of power produced would be given even it is the same reference load and generation profiles. Even more, the bottleneck in the system might also be changed when different transmission lines are disable. of the AC transmission line is good indicator to check the bottlenecks in the system. Table 2 gives the average utilization rate under the N-1 constraints for the four different study hours. The transmission line between the "Desbrosses"

and "Schoelcher" nodes has the highest utilization rate. The utilization rate for this transmission line is more that 80 % except for the off-peak hour without solar power (78 %). Therefore, we consider this transmission line to be the bottleneck in the system. For the other transmission lines, the utilization rates generally are less than 60 %, therefore, they are not considered as the bottlenecks in the system for our analysis.

Finally, the average nodal prices given by the economic dispatch model are summarized in Table 3.

Results of Step 2

Here, we cluster the nodes with similar nodal prices into zones, while excluding the bottleneck from them.

For nodes Lamentin, Petit Bourg, Marin, François, Galion, Trinité, Marigot, the nodal price differences for these nodes are less than 1% (Table 3). This means that there is no congestion between theses nodes. Therefore, it is reasonable that these nodes are grouped into a zone. The other nodes could be grouped into two or three zones, depending on the clustering criteria, i.e., whether a zone can only include a single node. Figure 4 gives two proposed zonal configurations. The 3zone configuration requires that there should be at least 2 nodes within the same zone (Eq. ( 11)), while the 4-zone configuration allows for zones with a single node. The objective value for the clustering model, as showed in Table 4, will be improved if we allow there is only one node in a zone, as show on the right hand side of Figure 4. Larger improvement could be found in the cases with high peak hours. Case II:4-zone

Results of step 3

It is difficult to decide the NTC value accurately for the long-term UC problem.

One possible way is to use the solution provided by the nodal pricing as a benchmark. The NTC constraints are used in the European day-head market, which sets a limitation on the commercial power exchanges among countries. An ideal NTC value in this case should allow the maximum utilization of the network.

However, due to the loop flow effect, the real time dispatching given by the dayahead model might overload the lines and therefore, needs to be rescheduled. This rescheduling is also called re-dispatching or counter trading. One of the reasons for re-rescheduling is that the NTC model does not consider the limitation of the network constraints within a zone. It is important to point out that re-dispatching will lead to an extra cost.

However, the purpose of our cases study is not to increase the power exchange among zones, but to increase the efficiency of the whole system. Many studies show that nodal pricing model is consider to give a better result than NTC model.

This has been proved by a large stream of literatures [13][15]. Therefore, we use the nodal economic dispatch model solution to calculate the NTC value.

Most congestion happens in the peak hours. NTC can limit the exchange in those hours and therefore could decrease the need for re-dispatching. In the off-peak hours, we could have less restriction because most of time, there is no or less congestion, and therefore the re-dispatching cost is not much.

Table 5 summarizes the power exchanges between zones for different configurations given by the second step model. The amount of power exchange for the hours without solar power is larger than that with solar power. In this case, if we use the result with the solar power, this might reduce the power exchange among zones and increase the total cost. Table 6 summarizes the total dispatching cost given different zone configuration and 2 different set of NTC value. We compare the number to the model with full network details. We find that the total dispatching cost in the 3-zone case is the one mostly close to the one with the full network. Furthermore, the cases with a higher NTC value have a higher dispatching cost than the case with a lower NTC value. However, the results here cannot sufficiently conclude that one zone configuration is superior to the other one. To fully evaluate it, we much consider the re-dispatching cost given the generation plan.

Another possibility here to evaluate the performance of different zone configuration and NTC value is to check if the generators are on in the "correct" time. Table 7 summarize the total on hours for each generator in different cases.

We compare these numbers to the cases with full network model, as show in Table 7. Generally, we could see the cases with a higher NTC model will perform better than the ones with lower NTC model. This is mainly because that the higher NTC value could cover more situation in the peak hours.

Figure 1 :

 1 Figure 1: Overview of the method

Figure 2 :

 2 Figure 2: Generation curve from Monday 22nd to Sunday28th, Oct. 2018.[source: opendata-martinique.edf.fr]

Figure 3 :

 3 Figure 3:HTB power network in Martinique [source: opendata-martinique.edf.fr]

Figure 4 :

 4 Figure 4: Proposed zonal configurations. The bottleneck is shown in blue.

Table 4 :

 4 Objective values for Step 2 (a) Case I: 3-zone configuration (b) Case II: 4-zone configuration

  nn), (c,cc), (l,ll), (i,ii) gdxxrw %datadir%%data%.xls MaxDupeErrors=999999999 cmerge=1 squeeze=N @temp.tmp" $gdxin %data% $load l n p map_pn map_ln z zonemember $load lineup nodeup genup $offmulti circuits(l) = lineup(l,"Circuits"); p_max(l) = lineup(l,"ThermalLimit")*circuits(l); voltage(l) = lineup(l,"Voltage"); slack(n) = 1$nodeup(n,'Slack bus'); react(l) = lineup(l,"Reactance")/circuits(l); resis(l) = lineup(l,"Resistance")/circuits(l); incidence(l,n) = 0; incidence(l,n)$(lineup(l,"From") eq n.val) = 1; incidence(l,n)$(lineup(l,"To") eq n.val) = -1; bvector(l) = (react(l) / (react(l)**2 + resis(l)**2))$(react(l) or resis(l)); gvector(l) = (resis(l) / (react(l)**2 + resis(l)**2))$(react(l) or resis(l));h(l,n) = bvector(l) * incidence(l,n); b(n,nn) = sum(l$incidence(l,n), incidence(l,n) * h(l,nn)); b(n,nn)$slack(n)=0; b(n,nn)$slack(nn)=0; g_max(p) = genup(p,'Generation capacity'); g_min(p) = genup(p,'Min Generation'); cg(p) = genup(p,'Cost'); Q(n) = nodeup(n,'Demand share')*210; g_wind(n) = nodeup(n,'Wind share')*0; g_solar(n) = nodeup(n,'Solar share')*0; g_bio(n)= nodeup(n,'Biomass share')*0

  𝑐𝑐,𝑐𝑐 -𝐺𝐺 𝑐𝑐,𝑐𝑐-1 ≤ 𝑔𝑔𝑠𝑠𝑝𝑝 𝑐𝑐 Where 𝑈𝑈 𝑐𝑐,𝑐𝑐 , 𝑉𝑉 𝑐𝑐,𝑐𝑐 , 𝑊𝑊 𝑐𝑐,𝑐𝑐 are binary value. 𝑈𝑈 𝑐𝑐,𝑐𝑐 is the on/off status of the power plant c at time t, 𝑉𝑉 𝑐𝑐,𝑐𝑐 is the startup action and 𝑊𝑊 𝑐𝑐,𝑐𝑐 is the shutdown action. 𝐺𝐺 𝑐𝑐,𝑐𝑐 is the generation of plant c at time t, and 𝑆𝑆 𝑐𝑐,𝑐𝑐 represents the partial online generating capacity or off-line generation resources.

		(25)
	0 ≤ 𝑆𝑆 𝑐𝑐,𝑐𝑐 ≤ 𝑔𝑔 𝑐𝑐 𝑚𝑚𝑚𝑚𝑚𝑚	(26)
	∑ 𝑆𝑆 𝑐𝑐,𝑐𝑐 ≥ 𝑔𝑔𝑔𝑔 𝑖𝑖,𝑐𝑐 𝑐𝑐	(27)

Eq. (

22

) is the UC objective function, which is composed of two component costs.

The first component cost is determined by the startup 𝑉𝑉 𝑐𝑐𝑐𝑐 decision and shutdown decision 𝑊𝑊 𝑐𝑐𝑐𝑐 on each generator (𝑔𝑔𝑏𝑏 𝑐𝑐 is the startup cost and 𝑔𝑔𝑔𝑔 𝑐𝑐 is the shutdown cost). The second component cost comes is primarily made up of fuel cost 𝑐𝑐𝑐𝑐 𝑐𝑐 and possible unserved energy penalty 𝑉𝑉𝑉𝑉𝐿𝐿𝐿𝐿 (i.e., load loss 𝐿𝐿𝑆𝑆 𝑐𝑐,𝑐𝑐 ). This penalty is usually used to avoid the load-shedding, where scheduled generators are not able to satisfy demands.

Table 1

 1

	, which includes two peak hours and two off-peak hours with
	different solar generation.

Table 1 :

 1 Reference load and generation for the case study

Table 2 :

 2 Average utilization rate of the transmission linesData for the four study hours are used as input to run the first steps model. The utilization rate (i.e., the actual power flow divided by the line nominal capacity)

Table 3 :

 3 Average nodal prices (EUR/MWh)1

		Peak hours	Off-Peak hours
	Node	With solar Without Solar With solar Without Solar
	Saint Pierre 148.02	148.01	146.10	154.33
	FL	148.06	148.06	146.10	154.33
	Bellefontaine 148.02	148.01	146.10	154.33
	Schoelcher 146.79	148.21	145.89	155.93
	Desbrosses 160.26	162.27	147.39	157.05
	Dillon	157.54	159.10	147.22	156.56
	Hydrobase 157.54	159.10	147.22	156.56
	Lamentin	156.76	158.21	147.18	156.44
	Petit Bourg 156.76	158.25	147.18	156.44
	Marin	156.76	158.39	147.18	156.44
	François	156.76	158.56	147.18	156.44
	Galion	156.76	158.67	147.18	156.44
	Trinité	156.76	158.67	147.18	156.44
	Marigot	156.76	158.67	147.18	156.44

Table 6 :

 6 Total dispatching cost (1000 EURO)

Table 7 :

 7 Total on hours for the generators

	Generator	3 zone Lower NTC	High NTC	4 zone Lower NTC	High NTC	Full network
	1	2903	3644	1981	3103
	2	8760	8760	8760	8760
	3	7061	273	5747	0	698
	4	8544	5917	8760	8541
	5	8760	8434	8760	7232
	6	8760	8760	8760	8039
	7	6528	7303	8758	6934
	8	127	2334	0	5754
	9	0	4379	0	3972
	10	368	2569	0	1770
	11	0	1465	0	2262
	12	2197	6525	8719	4836
	13	8757	6480	8760	7769
	14	8721	6480	0	6029
	15	8760	8756	8723	5519
	16	28	261	87	89	732
	17	0	5	0	0	77
	18	1	0	0	205	62
	19	0	21	0	0	46
	20	4	14	13	6	53

Table 8 :

 8 On hours compared to the full network modelAt this stage, we still could not conclude which zone configure would be a better choice for the long-term UC model with simplified network constraints. On of the main results is that we must consider the re-dispatching cost given by the generation plans. Another issue is about how to define the NTC values. as we can see from the preliminary results, the NTC values have a large impact on the results.In this study we use the same NTC value for all the hours. However, in the future study, we might need to designate different NTC values to different hours.

	IV. Future study

n Wind generation at node n [MW] gbio n Biomass generation at node n [MW] ntc z,zz Net transfer capacity between zone z and zz [MW] pmax n,nn Thermal transmission limit of transmission line l, n, nn [MW]

II: GAMS Code for step 1 and 3

	$ontext
	Author: Hong Cai
	$offtext
	option solprint = on;

The average prices here exclude the case when load shedding is enforced in the model. The average prices here are lower than the price in realities. This is because that we don't consider the startup cost in the first step model and assume all the power in the network are available.

The research topic and the test case have been proposed by EDF R&D. The study has been carried out independently by the candidate and the results do not engage the responsibility of EDF R&D.

This research benefited from the support of the Fondation Mathématique Jacques Hadamard program PGMO (PGMO-IROE no. P-2019-0014).

Table 5: Power exchanges based on nodal solution. Case I: 3-zonal configuration (Fig. 4a), Case II: 4-zonal configuration (Fig. 4b).

Results of step 4