Jean Philippe Praene1
email: praene@univ-reunion.fr

Fiona Bénard-Sora1

The ongoing energy transition: A comparative analysis of SIDS and European ultraperipheral islands

Keywords:

The vulnerability of island spaces has been recognized since the Earth Summit in Rio in 1992. Global energy pressure and the depletion of resources exacerbate their vulnerability. These territories are supposed to be real laboratories to experiment energy planning strategies, then allowed to help the international community. Nevertheless, the great diversity of island situations makes their case particular and difficult to treat. Some island territories are autonomous, the Small Island Developing States, while others depend on developed metropolitan states. Is there a difference between these two types of territories to succeed energy transition? We adopt electricity generation as a metric for sustainable transition quality. This paper contributes to the understanding of ongoing dynamic through a multidimensional analysis of data allowing to classify these territories: SIDS and outermost regions of the European Union. The combination of Principal Component Analysis and Hierarchical Classification will provide island profiles that will highlight differences between outermost regions and SIDS. Are the outermost regions of the European Union more advanced in the energy transition? In addition to the classification, a ranking will be proposed through the implementation of an energy vulnerability indicator. The classification and ranking will make it possible to draw typical portraits. These will identify the elements of resilience and the elements of vulnerabilities. This analysis, therefore, highlights the strengths and weaknesses will help to guide the energy policies of these territories. The comparison of these two types of territories has not been made so far. This work would highlight the situation of SIDS face island developed territories: outermost regions of the European Union.

INTRODUCTION

Since the 1970s and 1980s and the first major ecological disasters, energy and the exhaustion of energy resources are increasingly present in politico-environmental concerns. Oil spills, oil shocks, and nuclear accidents have made the energy transition an actual social event. The energy transition refers to the transition from an economy based essentially on fossil resources to a lowcarbon energy mix. All the territories of the world are concerned, but the unequal distribution of resources, wealth and the considerable heterogeneity of physical constraints and locations weaken some territories more than others. In these fragile territories, island territories are particularly sensitive to the various consequences of global warming. Although international opinion generally accepts energy transition as essential, its dynamics seems challenging to begin on a global scale.

The island territories must, therefore, take even closer ownership of the issue and the effective implementation of the energy transition. There are more than 50,000 islands spread around the world, for 740 million inhabitants [START_REF] Kuang | A review of renewable energy utilization in islands[END_REF], of which 2,000 are considered small islands [START_REF] Blechinger | Assessment of the Global Potential for Renewable Energy Storage Systems on Small Islands[END_REF]. These territories are very diverse, with different areas, population sizes, geographical features and political status [START_REF] Taglioni | Les petits espaces insulaires face à la variabilité de leur insularité et de leur statut politique[END_REF]. Among these territories, the international community has noted since the 1970s with formalization at the Earth Summit in Rio in 1992, the great particularity of Small Island Developing States [START_REF] Angeon | Les petites économies insulaires : un monde à part entière? Nouveaux regards conceptuels et méthodologiques[END_REF]. These isolated territories with fragile economies [START_REF] Angeon | Les petites économies insulaires et le développement durable : des réalités locales résilientes ? Etudes Caribéennes[END_REF][START_REF] Bouchard | Les petits États et territoires insulaires face aux changements climatiques : vulnérabilité, adaptation et développement[END_REF][START_REF] Briguglio | Small Island Developing States and Their Economic Vulnerabilities[END_REF][START_REF] Encontre | The vulnerability and resilience of small island developing states in the context of globalization[END_REF]have indeed been recognized as individual cases for the environment and development [START_REF]Small Island Developing States[END_REF]. The Summits of Barbados (1994) and Samoa (2014) subsequently confirmed the interest of the international community in SIDS ( UN, 1994). The interest of the international community is based not only on the vulnerability of these territories but also on the issue of the energy transition and the pressure it generates. SIDS is at the forefront of international climate action [START_REF] Gay | Small island developing states and their suitability for electric vehicles and vehicle-to-grid services[END_REF]. Indeed, because of their physical and socio-economic characteristics, SIDS can be real laboratories of experimentations to apprehend the energy transition. The challenge of understanding these territories is therefore global.

Other island territories have also been highlighted for their vulnerability and "usefulness" to the international community. The outermost European regions (ORs) are considered as isolated territories, with a development below that of the continental European zone. As SIDS, their physical and socioeconomic characteristics make them real laboratory experiments to address this critical issue of the XXI century: the energy transition. These territories are French, Portuguese and Spanish.

Less fragile than SIDS, these ultramarine territories are lagging behind their metropolises [START_REF] Sudrie | Quel niveau de développement des départements et collectivités d'outre-mer? Une approche par l'indice de développement humain[END_REF]. A working paper by the Agence Française de Développement shows a delay in the 1990-2010 period around -13% for the French ORs, -4.88% on average for the Portuguese ORs and -1, 73% for the Spanish ORs [START_REF] Sudrie | Quel niveau de développement des départements et collectivités d'outre-mer? Une approche par l'indice de développement humain[END_REF]. All the outermost regions are therefore behind their metropolises, but their situations are very diverse and unequal. A ranking according to the 2010 Human Development Index (HDI) of the ORs and SIDS easily shows that even if a majority of them are at the top of the table, the ORs do not have a clear advantage regarding the standard of living (as defined by the HDI) (Figure 1). The hypothesis of insularity as a discriminating factor seems relevant here. The remoteness and isolation of these territories could be aggravating factors for their development.

Fig. 1 The HDI ranking of ORs and SIDS

The notion of development is today more and more linked to the notions of sustainability and energy transition. Indeed, with the exhaustion of resources, global warming, the increase in the world population and the spread of cities, the energy transition is one of the critical issues of our contemporary societies. All territories around the world are concerned. But not all territories are prepared in the same way. Besides, barriers to the penetration of renewable energies in a territory can be numerous [START_REF] Painuly | Barriers to renewable energy penetration; a framework for analysis[END_REF]. Island territories are among the most vulnerable. However, again, these territories do not all have the same characteristics. For example, accessibility to energy is not always the same, and this is a particular development issue for some territories [START_REF] Surroop | Energy access in Small Island Developing States: Status, barriers and policy measures[END_REF]. In order to study these differences, we propose to focus on two types of island territories: ORs and SIDS. The main difference between the ORs and SIDs is that the ORs are not islands states and depend on a sovereign state (France, Spain or Portugal). These ORs are often not well known because most of the data and analysis is done at the national level. These overseas regions rarely appear in international studies, unlike SIDS.

Considering the hypothesis of a difference between the ORs and the SIDS in order to succeed in the energy transition, electricity generation is adopted as a metric for sustainable transition quality. This paper will look at first to contribute to the understanding the sustainability situation by drawing profiles and then to highlight development differences related to the energy transition by ranking the studied islands. In order to reach the first objective, a multidimensional data analysis will be practiced: the combination of Principal Component Analysis and Hierarchical Classification will allow classifying these territories. This two-step methodology will allow segmenting the studied islands into several clusters and will provide island profiles for the clusters. These will identify the elements of resilience and the elements of vulnerabilities. This analysis, therefore, highlights the strengths and weaknesses will help to guide the energy policies of these territories.

Thanks to the clustering, a ranking will be proposed in a third step through the implementation of a sustainability index. Such an indicator is essential for understanding the energy transition of a territory. It makes it possible to quantify the trajectory of transition and to highlight the most advanced and least developed territories. This ranking will show if the ORs are more advanced in the energy transition. This ranking provides an overview of the current situation of the different islands. We can thus observe the position of the ORs in relation to SIDS, which makes it possible to discuss the influence of a territory's status on its classification.

MATERIALS AND METHOD

Data collection

The implementation of this work, the first step is to identify the islands to study. At the OR level, the European Union has identified nine of them in the world: five French overseas departments (Martinique, Mayotte, Guadeloupe, French Guiana and Reunion), a French overseas community (Saint-Martin), two Portuguese autonomous regions (Madeira and the Azores) and a Spanish autonomous community (the Canary Island) [START_REF] Azevedo | Fact Sheets on the European Union[END_REF]. In all these very diverse territories, we choose to study only seven of them, eliminating the following regions:

-French Guiana because of its particular geographical location compared to other ORs as it is attached to the South American continent.

-Saint Martin because of the lack of data available for this region.

The United Nations (UN) provides a list of 58 SIDS around the world, 38 of which are members of the United Nations. Among the SIDS UN-members, nine are located in Atlantic, Indian Ocean, Mediterranean and South China Sea (AIMS), 16 in the Caribbean and 13 in the Pacific (UN, n.d.).

There are therefore very diverse territories among which Seychelles, Singapore, Comoros, the Dominican Republic, Haiti, Jamaica or Fiji Islands, Papua New Guinea, and the Solomon Islands to name a few. For reasons of availability of data for the selected study variables, 23 SIDS and 7

ORs has been selected for this study, see Table 1. In order to analyze the progress of territories in the energy transition, eight variables were selected for the year 2014 (Table 2). For SIDS, the majority of data come from Data World Bank1, except for CO2 emissions per capita2 and the insularity index that we had to compute according to the method [START_REF] Guerassimoff | Iles et énergies : un paysage de contrastes(Presse des)[END_REF]. For the outermost regions, the data come from Eurostat (Eurostat, 2018) and national and regional sources of statistical data on energy and economy (EDF, 2018a[START_REF] Edf | EDF -La Réunion[END_REF][START_REF] Edf | EDF -La Réunion[END_REF](EDF, , 2018b;;[START_REF] Istac | Canarias en Cifras[END_REF]Pordata, 2018). We selected these variables because of the availability of data for both types of islands. Indeed, all types of data on the ORs are challenging to access because they are generally included in national reports. Besides, the data were chosen to facilitate the reproducibility of the method over time on future databases. This data evaluates the extent of electricity needs, that is the central feature of our analysis.

The share of renewable energies in the electricity mix (ELC.RNEW) % This data makes it possible to evaluate the level of penetration of renewable energies in the electricity mix. It assesses the current choices and opportunities of territory regarding renewable energy.

CO2emissions per capita (CO2.PC) ton metrics/capita

This data is related to the economic development of a territory. The report to the population evaluates lifestyles.

Insular Index (INS.I) -

The insular situation is measured through the remoteness of the territory but also its political system. The calculation method comes from [START_REF] Guerassimoff | Iles et énergies : un paysage de contrastes(Presse des)[END_REF] 2.2 Methodology

PCA and Hierarchical clustering

Principal Component Analysis is an exploratory method of analyzing data from an n*p table of quantitative data, with n the individuals studied and p the variables studied. The goal is to reduce p-variables to new q-dimensions created in order to synthesize the information (q<p). With the creation of new q-dimensions, PCA also offers simple data and variables visualization method:

-n-individuals are optimally represented by minimizing the deformations of the initial point cloud, in a subspace Eq of dimension q (q<p) (individuals factor map);

-the variables are also represented graphically in a subspace Fq by explaining the original links between the variables (variables map).

The part of inertia explained (makes it possible to measure the quality of representations:

𝑟 𝑞 = ∑ 𝜆 𝑖 𝑞 𝑖=1 ∑ 𝜆 𝑖 𝑝 𝑖=1
(1)

Where 𝜆 𝑖 is the eigenvalues. The value of 𝑞 is chosen so that 𝑟 𝑞 is greater than a fixed threshold 𝜁 . We have: 𝑟 𝑞 > 𝜁 .

In our study, PCA has two objectives: firstly, PCA is used to highlight groups of individuals with similar characteristics; secondly, PCA allows to identify the variables that have a significant influence on the variance of the most influential principal component. It means that in our study, the PCA method allows us to see which are the most influential variables in the question of the energy vulnerability of island territory.

The hierarchical classification on principal components makes it possible to identify groups (clusters) from the data of the ACP. Thus, the method is based on the q-dimensions and the coordinates of the n individuals on the q-dimensions to carry out the classification. The graphical representation is advantageous and straightforward because it takes up the factorial axes of the ACP (the q-dimensions) to visualize the clusters. The classification methodology chosen is Ward's method. It consists in grouping the clusters so that the increase of the inertia inter-classes is maximum, or so that the increase of the intraclass inertia is minimum.

Creating an index and Ranking

In addition to the classification, our study proposes the ranking of the territories studied thanks to the creation of a sustainability index (SI). The purpose of this index is to show how far a territory is advanced in the energy transition. It represents three key challenges: developing its territory, conducting a proactive policy by promoting renewable energies and acting on consumption patterns (Benard-Sora & Praene, 2017). Our index will try to translate these difficulties. This index can be understood in two ways: a first way would be to see the index as a vulnerability index (Bénard-Sora & Praene, 2017), a second way would be to see it as an index of sustainability [START_REF] Praene | Sustainable transition in small island developing states: Assessing the current situation[END_REF]. In this work, we choose to speak of sustainibility index (SI), combining both points of view. When this index is reliable, it reflects the sustainability, and conversely, when it is weak, it reflects the vulnerability.

Although the concept of sustainability has been around for a long time, it continues to evolve.

Indeed, in the literature, new technical-economic or institutional techniques are now being discussed. In our article, given the availability of data, we have based ourselves solely on the three fundamental pillars of sustainable development. We have thus assumed that the selected variables succeed in providing a sufficient representation of the concept of sustainability

To create such an index, we consider the results of the PCA, and in particular, the dimension (PCi) created. The index (SI) is defined as follows:

𝑆𝐼 = ∑ 𝜆 𝑖 𝑃𝐶 𝑖 𝑛 𝑖=1 ∑ 𝜆 𝑖 𝑛 𝑖=1
(2)

Each dimension ( 𝐷𝑖𝑚 𝑖 ) must then be decomposed according to the variables 𝑝 𝑗 that participated in its construction during the PCA: it is a linear combination of the variable set. Thus, each variable is weighted by a weight 𝛼 𝑖𝑗 calculated through the PCA. 𝑃𝐶 𝑖 is defined as follows:

𝑃𝐶 𝑖 = ∑ 𝛼 𝑖𝑗 𝑝 𝑗 𝑗 (3)
The variables in the input are of very different natures and order of magnitude. Therefore, it is necessary to proceed before calculating the index to a normalization of the data of the table n*p.

For this normalization, we choose to perform a z-score normalization. Our standardized data so have a mean of 0 and a standard deviation of 1. The particularity of this standardization is the fact it does not affect the shape of the original distribution.

RESULTS AND DISCUSSION

Exploratory data analysis

The objectives of our paper are threefold: i) Assess the determinants of energy vulnerability under insular context. ii) Identify disparities or similarities between different type of island using hierarchical clustering. iii) Create an indicator to help inform the energy transition issue.

In the first step, a principal component analysis (PCA) has been investigated. Most of the time PCA is conducted to reduce the size of the original set of data. In this study, the purpose was not so much the reduction question, but the ability to define some new latent variable. This latent variable is the principal component which is defined by a linear combination of original data. These combinations are useful to define the main characteristics of each principal components.

Our approach was conducted a PCA iteratively. Indeed, the initial full dataset consists of 36 islands. We first explore two PCA calculation. The analysis of the projected plan has highlighted some outliers. These correspond to particular islands cases which actively contribute to the construction of the plane. Six islands (DOM, HTI, PNG, SGP, SUR, TTO) has been identified, their cumulative contribution to the contribution of the first two principal components equals 89.1%. Subsequently, we have removed these cases for the rest of the clustering and PCA. The studied dataset is finally 30 islands. A commonly used method to define the significant numbers of the principal component to consider is to fix the cut-off criteria of an eigenvalue ƛ>1 as defined by Kaiser's rule. As suggested by the results in Table 3, the variations the first three PCs is taken into account. The -LND.TOTL.K2 and ELC.RNEW are positively related to dimension 2. The share of renewable energies in the electricity mix seems to play a close link with the area of the territory.

Many of the territories considered in the study favored solar energy and hydropower. Both of these technologies require a lot of available land area.

-CO2E.PC and POP.DNST are both negatively related to dimension 2. It indicates that these two variables are related: the denser the territory is, the higher its CO2 emissions are. However CO2 emissions have a better projection on the second plan.

-The share of renewable energies seems to be linked to CO2 emissions: the higher this share, the lower the CO2 emissions. Renewable energies are indeed an alternative to fossil fuels, which are highly carbonaceous. The results of the PCA will therefore clearly distinguish islands where renewable energies are already deployed and whose carbon emissions are low. islands and, to a lesser extent, Reunion (RUN) and Bahamas (BHS).

From Figure 3b, we can see that the island is not situated in a specific part of the individuals map.

The third principal component opposes island such as Comoros (COM) and Maldives (MDV) to Acores (ACO), Madeira (MAD), Reunion (RUN). Islands located in the upper part of the plan are heavy electricity consumer. This group also has high values for the share of renewable energy in electricity generation.

This first step of exploratory analysis helps us to understand which variables have significant influence in the definition of the principal components. These parameters have to be taken into account to assess the energy vulnerability of the islands.

Hierarchical clustering on principal components (HCPC)

The second step of the study consists of a hierarchical clustering on the main principal components (Dim 1 and Dim 2). These dimensions correspond to the first two principal components, which represent the best projection plane according to the PCA. HCPC analysis is a widely-used method which investigates the organization of individuals in clusters and among clusters depicting a hierarchy. The selected partitioning approach selected for this study is the K-means method. This method is one of the most commonly used clustering algorithms. In this approach, the center of gravity also characterizes a cluster. The number of clusters is defined thanks to a statistical optimization algorithm called the Gap statistic method. This method proposed an optimum value of three clusters, and the results are shown in Figure 4. Groups are characterized in two ways: the variables and the islands. Each cluster is defined by some categorical variables which have a significant p-value < 0.05. The variables which significantly described the cluster could also be identified by test-value v-test > 1.96 which is equivalent to a pvalue < 0.05. The individual who corresponds in our case to the islands is mainly observed through the paragon. The paragon represents the islands which are the nearest from the center of gravity of cluster. Thus, that means that this island is the one that can qualify for representing the mean behavior of the group.

The description of the three clusters is detailed above:

• Cluster 1: Four variables most characterize the first cluster. Insularity and total land area have some high values for these clusters compared to the overall mean of the data. Islands in these groups also have some low value for two other variables density population and CO2 emission per capita. Indeed, this presents groups has a share of renewable energy in the electricity of 21.6% compared to 12.9% for the overall dataset.

• Cluster 2: The variables characterizing this cluster are not so much different to those characterizing cluster 1. However, their intensity is opposed. Cluster 2 is depicting islands which have a higher CO2 emission per capita due to low penetration of renewable sources in electricity production (10.7%). This cluster represents small islands with low insularity index and the lowest mean total population. The high values of standard deviation for all variables constitute another particular feature of this value. Standard deviation explains, in fact, dispersion

• Cluster 3: This last cluster is composed of Jamaica, Mauritius, and Canary. This group is a set of islands with high values for electricity consumption, GDP and population. These islands have a low insularity and are heavily dependent on fossil fuels for electricity generation. Except for Mauritius which reach to 20%, the other two islands have a share of renewable energy around 10%.

This hierarchical clustering reveals no specificity regarding organization between SIDS and ORs.

Indeed, one expected result from the first of this study were to identify any specific disparities between SIDS and ORs. Even if certain ORs have a better projection on the third principal component, all the ORs are part the cluster 2. This second group is the biggest, and both contain SIDS and ORs. Thus, we could conclude that the status of the island from our dataset does not influence the definition of a specific cluster for the ORs. Indeed, these regions are provided with the recognition of necessary derogations in European law. Thus, a first hypothesis was to observe the effect, thereby differentiating these regions by a specific group mainly composed by the ORs.

The clustering shows that ORs behavior and characteristics are the same as SIDS with some medium GDP.PPA, POP.TOTL, INS.I, but some high value for C02 emission due to economic development and low renewables sources integration in the electricity mix. A particular point that can be noticed is that 43% of the SIDS are in cluster 1. Moreover, this cluster is only composed of SIDS.

Sustainability indicator

The composite indicator considered in this section is based on the previous PCA analysis. PCA is a well-established method of dimensionality reduction. The eigenvectors (modes) which represent the maximum variance directions have been defining previously. The first three components explain 71.2% of the total variance of the overall data set. Thus, the first three modes are fixed as a consistent cut-off for the definition of the composite indicator. The diagram presents in Figure 5 shows the way our composite indicator (SI) is defined. Based on the PCA, SI expects to capture the complexity of sustainable energy dimensions for the islands. As shown, the indicator summarizes several aspects of sustainability: environmental, economic, energetic, social and geographical dimensions. As shown in Figure 5, the original weighting scheme for each variable is not an equal weight approach. The contribution to a principal component (PC) construction gives the weight of each variable in the definition of a PC. It can be noticed that some variables are part of two different principal components such as the share of renewable sources or insularity index. Thus, the PCA approach implicitly gives more importance to the variables and consequently to the themes (environment, social, economy, energy) which significantly explain the variance of the overall data set. We built a ranking based on a PCA weighting approach with linear aggregation and Z-scores normalization.

The score of sustainability proposed by our indicator SI is presented in Figure 6. Small islands are generally assumed to be vulnerable territories mainly driven by common vulnerabilities such as high carbon energy production, isolation. A high score depicts the level of sustainability. The islands ratings are sorted in descending order. A remarkable point is the fact that the ORS is spread across the entire graphic.

Regarding the indicator, Mauritius appears as the most sustainable territory. This result is mostly due to high values on economic, energy and environment aspects. The three islands from cluster 3 (Mauritius, Jamaica, Canary) have the top score for SI. The graphic shows two levels: the first one between MAD and ATG, and the second from DMA to STP. These two levels indicate that most islands are very similar according to their principal characteristics. The two steps of sustainability (blue and orange) are suggested that most of the islands have a low level of sustainability. Their vulnerability is explained renewable energy sources in electricity mix that is very low, less than 6% except for the case of Dominica. Indeed, hydropower currently accounts for 21% of Dominica's total electricity mix.

Focusing on the second step shows in orange, a disparity should be highlighted between ORS and SIDS. Even if these islands are on the same stage, some difference is visible while exploring the database. Indeed, the of renewable energy is relatively high (more than 20%). The SIDS in this step has a mean value of around 10%. The particular case is Samoa (WSM) is also explained by 33% of hydropower in their electricity mix. Another remarkable point is that ORS insularity index is less than half that SIDS (40.9 versus 148.3). As the indicator, does not include any governance, policy parameters or variables that could qualify the specific status of SIDS or ORS, no disparity could be noticed between the two types of studied islands.

It is observed that the organization of the islands described by SI is in the same trend of hierarchical clustering. Thus, it could be supposed that the indicator summarizes correctly the information contained by the dataset. The index has demonstrated its usefulness to simplify the analysis of the data set.

CONCLUSION

It is widely agreed that small islands are among territories of opportunities to experiment innovative policy or solutions to face their sustainable transition. In the same time, these islands are generally considered as the most vulnerable territories. With this study, we have provided a statistical approach for the construction of indicator that can be useful to policy-makers that want to create their index of sustainability. The use of principal components has helped us to select variables which explain a significant part of the total variance. This step has been used to select variables. Then a hierarchical clustering based on the principal components shows that the islands are projected into three groups. This trend has also followed the index SI. The proposed index summarizes the principal aspects of sustainability which is given by available online data. This work shows that at present, the island nature is decisive and no significant differences have been identified between the two groups. As the ORs belong to so-called rich European countries, the comparison with islands states that do not belong to this category is interesting. Indeed, the current situation does not show any specific disparity between the types of islands. Their mode of governance, therefore, does not seem to be a discriminating factor concerning sustainability issues. However, the implementation of the transition will require substantial financial engineering that allows for proper investment planning. It is in this sense that it will be interesting to observe the evolution of the values of our indicator. Indeed, given the European Union's support for the ORs, it is likely that the transition dynamic may be stronger than for many SIDS. This gap is likely to be very significant if some SIDS do not have external support for the structural transformation of energy production in their territory. Our work has therefore made it possible to evaluate in a first step a current state, which can be taken as a reference. This situation could be used in future work to qualify the dynamics of the transition to a more sustainable territory in the future.

However, the ranking respects the overall trend define in the hierarchical sectionnevertheless, the lack of data which characterize political status, governance, and energy policy. A point of interest is to investigate alternative weighting by introducing the negative and positive effect of variables based on the correlation and contribution of principal components construction. This will be the aim of further research investigation.

  Variable Factor Map allows to analyze the relation of the variables with each other and their link with the dimensions created, see Figure 2a. The representation of the variables on the two first projection plan highlights four main conclusions: First plan (Dim1, Dim2) -Three variables are mainly related to dimension 1: GDP.PPA, POP.TOTL and ELC.CONS.GWH is positively related to dimension 1 while INS.I is negatively related. This finding clearly indicates the opposite relationship between the INS.I variable and the other variable remainder in the data set. Insularity seems to play a particular role in the characterization of the studied territorieshowever, the quality of projection of the INS.I on the first plan is not good. Thus, it seems important to observe and discuss the insularity on the second projection plane.

Figure 2 .

 2 Figure 2. Islands scores on the first and second principal components.

Figure 4 .

 4 Figure 4. Hierarchical clustering projection of the first and second principal components.

Figure 5 .

 5 Figure 5. Diagram of the composite indicator construction.

Figure 6 .

 6 Figure 6. Sustainability indicator of SIDS and ORS.

Table 1 List of studied islands

 1

		Islands	ID	Regions
		Antigua and Barbuda	ATG	Caribbean
		Bahamas	BHS	Caribbean
		Barbados	BRB	Caribbean
		Belize	BLZ	Caribbean
		Cape Verde	CPV	AIMS
		Comoros	COM	AIMS
		Dominica	DMA	Caribbean
	SIDS	Fiji Grenada Jamaica	FJI GRD JAM	Pacific Caribbean Caribbean
		Kiribati	KIR	Pacific
		Maldives	MDV	AIMS
		Mauritius	MUS	AIMS
		Micronesia	FSM	Pacific
		Samoa	WSM	Pacific
		Sao Tome and Principe	STP	AIMS
		Seychelles	SYC	AIMS
		Azores	AZO	Atlantic
		Canary island CAN	Atlantic
	ORs	Guadeloupe Madeira Martinique	GUA MAD MAR	Caribbean Atlantic Caribbean
		Mayotte	MAY	Indian Ocean
		Reunion	RUN	Indian Ocean

Table 2 List of variables Variable Unit Special features

 2

	Gross Domestic Product (GDP.PPA)	$ Purchasing Power Parity (PPP)
	Population Density (POP.DNST)	Capita/km2	This indicator considers the pressure of urbanization for the energy transition.
	Electricity Consumption	Gigawatt hour
	(ELC.CONS.GWH)	(GWh)

Area (LND.TOTL.K2) km2

This data is a good indicator of the development potential of an island (economies of scale, renewable energies).

Population (POP.TOTL) Capita

This data makes it possible to quantify the need on the territory.

Energy consumption is closely linked to the economic development of a territory. PPP GDP is used to compare the GDP of different territories with different currencies. It allows expressing in a standard unit the purchasing power of different currencies.

Table 3 Eigenvalues of correlation matrix. Principal Component Eigenvalue ƛ Percentage of variance Cumulative percentage of variance

 3

	PC1	2,4	30,4	30,4
	PC2	1,9	24,1	54,5
	PC3	1,3	16,7	71,2

https://data.worldbank.org

http://hdr.undp.org/en/data


Chadi El Boustani
email: chadi.el-boustani@strains.fr

Jérémy Bleyer

Mathieu Arquier
email: mathieu.arquier@strains.fr

Mohammed Khalil Ferradi

Karam Sab
email: karam.sab@enpc.fr

Jeremy Bleyer
email: jeremy.bleyer@enpc.fr

Mohammed Khalil Ferradi
email: mohammed-khalil.ferradi@strains.fr

An efficient and

Keywords: Connection modelling, Dual approaches, Error estimators, Second-order cone programming, Interior point method

INTRODUCTION

The study of large scale problems including material non-linearities such as plasticity and boundary non-linearities such as contact using finite element method remains a problematic matter for inexperienced users. Non-linear finite element calculations therefore require some expertise when aiming at obtaining accurate estimates of quantities of interest such as pressure, openings or penetration, while trying to keep computational costs at a moderate level. In the present paper, we propose an alternative approach for computing structures in contact such as 3D steel assemblies by using a solution algorithm which exhibits very robust convergence properties and does not require any fine-tuning algorithmic parameters while still being very competitive in terms of computational costs compared to standard approaches such as penalisation and augmented Lagrangian techniques. Inspired from the classical approach for limit analysis [START_REF] Lyamin | Lower bound limit analysis using non-linear programming[END_REF][START_REF] Lyamin | Upper bound limit analysis using linear finite elements and nonlinear programming[END_REF][START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF], the dual approach provides upper and lower bounds for the system's potential or complementary energy and other quantities on interest for the engineer. The dual principles can be formulated as convex minimization problems where non-linearities can be easily introduced if expressed as a complementarity condition on self-dual cones. The solving algorithm used is a state-of-the-art interior point method (IPM) which relies on the latest advances of the mathematical programming community.

The formulation of elastic problems with contact is discussed in section 2 along with the main idea for the IPM in section 3. The extensions for error estimation, remeshing and material non-linearities are presented in section 4. The whole process is illustrated in a typical steel assembly example presented in section 5 using Strains's (https://strains.fr) DS-Steel software in which these methods have been implemented. The results are then compared with respect to computations made using Abaqus to ensure the validity of results and to compare computational cost and the quality of the results. Detailed analysis and results can be found in upcoming publication [START_REF] Boustani | Dual finite-element analysis using second-order cone programming for structures including contact[END_REF].

DUAL APPROACH TO MECHANICAL PROBLEMS

Contact constitutive equations

Contact conditions are enforced using a combination of Signorini's unilateral contact constraints and Coulomb's frictional constraints. When associated contact is assumed, the unilateral and frictional laws can be expressed over the contact surface Γ 𝑐𝑐 between the deformable bodies as follows:

� if ‖𝝈𝝈 𝑻𝑻 ‖ < -𝜇𝜇𝜎𝜎 𝑁𝑁 then 𝒈𝒈 𝑻𝑻 = 𝜆𝜆 𝝈𝝈 𝑻𝑻 ‖𝝈𝝈 𝑻𝑻 ‖ with 0 ≤ 𝜆𝜆 ≤ 𝑔𝑔 𝑁𝑁 𝜇𝜇 if ‖𝝈𝝈 𝑻𝑻 ‖ = -𝜇𝜇𝜎𝜎 𝑁𝑁 then 𝒈𝒈 𝑻𝑻 = 𝜆𝜆 𝝈𝝈 𝑻𝑻 ‖𝝈𝝈 𝑻𝑻 ‖ with 𝜆𝜆 = 𝑔𝑔 𝑁𝑁 𝜇𝜇 (1)
Which ensures that 𝜎𝜎 𝑁𝑁 ≤ 0, 𝑔𝑔 𝑁𝑁 ≥ 0 but also couples the normal and tangential gaps. The previous conditions can be rewritten as a pair of second-order cones using a set of kinematic (𝑔𝑔 𝑁𝑁 , 𝒈𝒈 𝑻𝑻 ) and static variables (𝜎𝜎 𝑁𝑁 , 𝝈𝝈 𝑻𝑻 ):

(𝜎𝜎 𝑁𝑁 , 𝝈𝝈 𝑻𝑻 ) ∈ 𝐾𝐾 𝜇𝜇 and (𝑔𝑔 𝑁𝑁 , 𝒈𝒈 𝑻𝑻 ) ∈ 𝐾𝐾 -1/ 𝜇𝜇 where

𝐾𝐾 𝛼𝛼 = {(𝑥𝑥, 𝒚𝒚) ∈ ℝ × ℝ 2 | ‖𝒚𝒚‖ + 𝛼𝛼𝑥𝑥 ≤ 0} (2)
Since 𝐾𝐾 -1/𝜇𝜇 is the polar cone of 𝐾𝐾 𝜇𝜇 , the pair of static and kinematic variables satisfies a complementarity condition over dual cones.

Energy principles as conic programs

Classical energy principles can be used to formulate finite elements problems as minimisation programs. These principles can be adapted to consider contact conditions and/or plasticity.

While the total complementary energy is written as follows:

𝛱𝛱 𝑐𝑐 (𝝈𝝈) = � 1
Where 𝕀𝕀 𝐾𝐾 -1/𝜇𝜇 and 𝕀𝕀 𝐾𝐾 𝜇𝜇 are the indicator functions for the plastic cones expressed in equation ( 2).

The two energy principles can be then expressed as a pair of dual conic programs:

1. A kinematic program (TPE) in which a kinematically admissible displacements field is calculated, equilibrium equations are thus verified in their weak form; 2. A static approach (TCPE) in which a statically admissible stress field is obtained, with kinematic constraints being verified in their weak form.

(

Indicator functions can be removed from the objective part and be passed as conic constraints obtaining then the classical from of quadratic programs.

Because of duality, both minimization problems provide respectively an upper and lower bound of the real potential energy of the system (after inverting the sign of the complementary energy):

-𝛱𝛱 𝑐𝑐 (𝝈𝝈 𝒉𝒉 ) ≤ -𝛱𝛱 𝑐𝑐 (𝝈𝝈 * ) = 𝛱𝛱 𝑝𝑝 (𝒖𝒖 * ) ≤ 𝛱𝛱 𝑝𝑝 (𝒖𝒖 𝒉𝒉 ) (6)
These bounds provide the engineer the capacity to judge the quality of convergence since one can evaluate the proximity of approximate solutions to the exact solution. The relative gap between the upper and lower bound therefore represents the above-introduced error estimator and serves as an excellent indicator of convergence which can be defined as:

𝛥𝛥𝛥𝛥 = 𝛱𝛱 𝑝𝑝 (𝒖𝒖 𝒉𝒉 ) + 𝛱𝛱 𝑐𝑐 (𝝈𝝈 𝒉𝒉 ) 𝛱𝛱 𝑐𝑐 (𝝈𝝈 𝒉𝒉 ) (%) (7)

INTERIOR POINT METHOD FOR CONIC PROGRAMS

Quadratic second-order cone programming

Both static and kinematic approaches result in a similar mathematical formulation of secondorder cone programs (SOCP). These energy principles are a special case called quadratic second-order cone program (QP) which consist in minimizing a quadratic objective function under linear and conic constraints such as the positive orthant or the second order cone as shown in equation ( 8).

𝑚𝑚𝑚𝑚𝑚𝑚. 𝒙𝒙, 𝒙𝒙 𝒊𝒊

1 2 𝒙𝒙 𝑻𝑻 𝑷𝑷𝒙𝒙 + 𝒄𝒄 𝑻𝑻 𝒙𝒙 + 𝒄𝒄 𝑰𝑰 𝑻𝑻 𝒙𝒙 𝑰𝑰 𝑠𝑠. 𝑡𝑡. 𝑨𝑨𝒙𝒙 = 𝒃𝒃 𝑨𝑨 𝑰𝑰 𝒙𝒙 -𝒙𝒙 𝑰𝑰 = 𝒃𝒃 𝑰𝑰 𝒙𝒙 𝑰𝑰 ∈ 𝐾𝐾 (8)
Most standard SOCP program formats assume no quadratic objective term, but the previous problem can be always reformulated as a problem involving a linear objective function only by introducing auxiliary variables and additional conic constraints [START_REF] Alizadeh | Second-order cone programming[END_REF][START_REF] Andersen | On implementing a primal-dual interior-point method for conic quadratic optimization[END_REF][START_REF] Boyd | Convex optimization[END_REF]. However, in our implementation, the quadratic term will be kept as such in the objective function since its physical meaning is interesting.

Interior-point method solution procedure

A lot of research has been carried out since the 1960s on numerical algorithms to solve optimisation problems. Today, these algorithms provide the engineer an efficient and reliable technology to solve convex programs especially SOCPs. We will focus on the primal-dual interior-point algorithm, a subclass of the interior-point methods, which is currently considered as one of the most efficient ways to solve QP problems. The idea of an IPM is to find a solution to the Karush-Kuhn-Tucker (KKT) conditions for problem (8) by following the neighbourhood of a curve called the central path and parametrized by a barrier parameter 𝜂𝜂 ≥ 0. The optimal solution is obtained for 𝜂𝜂 → 0. Find (𝒙𝒙, 𝒙𝒙 𝑰𝑰 , 𝝀𝝀, 𝝂𝝂 𝑰𝑰 , 𝒔𝒔 𝑰𝑰 ) such that

𝑷𝑷𝒙𝒙 + 𝒄𝒄 + 𝑨𝑨 𝑻𝑻 𝝀𝝀 + 𝑨𝑨 𝑰𝑰 𝑻𝑻 𝝂𝝂 𝑰𝑰 = 0 𝒄𝒄 𝑰𝑰 -𝝂𝝂 𝑰𝑰 -𝒔𝒔 𝑰𝑰 = 0 𝑨𝑨𝒙𝒙 -𝒃𝒃 = 0 𝑨𝑨 𝑰𝑰 𝒙𝒙 -𝒙𝒙 𝑰𝑰 -𝒃𝒃 𝑰𝑰 = 0 𝒙𝒙 𝑰𝑰 ∘ 𝒔𝒔 𝒊𝒊 = 𝜂𝜂𝒆𝒆 𝑰𝑰 𝒙𝒙 𝑰𝑰 ∈ 𝐾𝐾 , 𝒔𝒔 𝑰𝑰 ∈ 𝐾𝐾 * (9)
Where 𝒆𝒆 = (1, 𝟎𝟎) in the case of Lorentz cones. The IPM relies on an iterative algorithm that calculates descent steps using Newton's method on the residuals equations and introducing the Nesterov-Todd scaling matrices [START_REF] Nesterov | Primal-dual interior-point methods for self-scaled cones[END_REF]. Therefore, at each iteration, a linear Newton-like system is solved:

� 𝑷𝑷 + 𝑨𝑨 𝑰𝑰 𝑻𝑻 𝑭𝑭 𝟐𝟐 𝑨𝑨 𝑰𝑰 -𝑨𝑨 𝑇𝑇 -𝑨𝑨 0 � � ∆𝒙𝒙 ∆𝝀𝝀 � = -𝑹𝑹 (10)
With 𝑹𝑹 being the residuals vector and 𝑭𝑭 a specific scaling matrix depending on the types of cones. After each iteration, the variables are updated to ensure that the conic constraints will always be strictly verified which ensures the verification of the physical phenomena modelled by the cones (such as non-penetration in contact pairs and plasticity cones). Considering 𝑨𝑨 of full rank, the linear system (10) is never singular, thus convergence in ensured in a limited number of iterations (ranging from 15 to 25 iterations in most cases), which is a key advantage of this method.

In comparison, incremental Newton-Raphson algorithm, based on tangent matrix and residual force, often leads to problems of convergence, especially when the initial point is "far" from the solution. Engineers must divide the loads to perform step by step analyses until a satisfactory state is calculated.

Finite-element discretisation

While classical quadratic 10-node displacement tetrahedra are used in the kinematic approach, special improved equilibrium elements are necessary to handle the static approach. Over a mesh made of 4-node linear tetrahedra, a "generalized" stress field, which already includes the equilibrium equation: div(𝝈𝝈) + 𝒃𝒃 = 0, is interpolated. This reduces the number of degrees of freedom from 6 × 4 = 24 per tetrahedra to 24 -3 = 21 per tetrahedra and removes a first set of linear redundancies. Since the tetrahedra are considered disconnected, the second equilibrium equation: ⟦𝝈𝝈. 𝒏𝒏⟧ = 0, ensuring the continuity of the stress vector on element interfaces, should be written on all internal tetrahedron interfaces. This will cause enormous redundancies in the linear constraints matrix thus making the solved-system highly indefinite.

The solution to eliminate these redundancies is to use the super-elements method introduced in [START_REF] De Almeida | A set of hybrid equilibrium finite elements models for the analysis of three-dimensional solids[END_REF][START_REF] Kempeneers | Pure equilibrium tetrahedral finite elements for global error estimation by dual analysis[END_REF] and consists in subdividing each tetrahedron into 4 and then removing all internal degrees of freedom, leaving only the ones from external faces.

ADVANTAGES AND EXTENSION TO MATERIAL NON-LINEARITIES

Error Estimator and remesh procedure

To limit the number of elements used in 3D analysis, a general adaptive remesh scheme is implemented using a constitutive equation type of error based on duality. This type of error qualifies, both on the local and global level, the "closeness" of the kinematic and static solutions. In short, if the gap in constitutive relation in zero, both solutions are correct and present an accurate dual representation of the optimal solution. This versatile type of error, introduced by Ladeveze [START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF] based on the Fenchel inequality, allows us to consider each physical phenomenon such as plasticity or contact, separately. It can be written in the case of elastic continua with contact interfaces as follows:

𝑒𝑒 2 (𝒖𝒖, 𝝈𝝈) = 𝑒𝑒 𝛺𝛺 2 (𝒖𝒖, 𝝈𝝈) + 𝑒𝑒 𝛤𝛤 𝑐𝑐 2 (𝒖𝒖, 𝝈𝝈) ≥ 0 ∀(𝒖𝒖, 𝝈𝝈) (11)
Which is constituted by the elastic volume part of the error and the surface contact part of the error. Thus, this error can also be expressed as: 𝑒𝑒 2 (𝒖𝒖, 𝝈𝝈) = 𝛱𝛱 𝑝𝑝 (𝒖𝒖) + 𝛱𝛱 𝑐𝑐 (𝝈𝝈) [START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF] Therefore, the gap given by equation ( 12) is nothing else than a qualification of the global error and can be used to assess global convergence.

The complete procedure involves the following steps:

1. An initial static and kinematic calculation are performed over an automatically software-generated coarse mesh Ω (0) ; 2. On this mesh, an error calculation is made, and global, elementary and relative errors are computed along with the relative gap between objectives functions; 3. Nodal relative errors are then used to compute a new anisotropic mesh size map; 4. A new mesh Ω (1) is generated and new static and kinematic analyses are performed. The process is repeated until a target objective gap is reached.

An extension to material plasticity

IMP provides a reliable solution for non-linear convex problems that can be written in the form of conic problems. The most interesting part of this method is its ability to provide solutions for problems traditionally unsolvable using a Newton-Raphson procedure such as limit analysis. The use of this method is rather recent [START_REF] Kanno | Nonsmooth mechanics and convex optimization[END_REF][START_REF] Lyamin | Lower bound limit analysis using non-linear programming[END_REF][START_REF] Lyamin | Upper bound limit analysis using linear finite elements and nonlinear programming[END_REF][START_REF] Makrodimopoulos | Lower bound limit analysis of cohesive-frictional materials using second-order cone programming[END_REF][START_REF] Makrodimopoulos | Upper bound limit analysis using simplex strain elements and second-order cone programming[END_REF] and is being generalized in various yield problems in non-smooth mechanics such as viscoplastic fluids [START_REF] Bleyer | Advances in the simulation of viscoplastic fluid flows using interior-point methods[END_REF].

Considering material non-linearities is rather a direct extension of the dual approach exposed in this article provided that the yield condition can be written in the form of a complementarity condition over second order cones and the dissipation derives from a potential. That is the case with most standard materials, specifically steel, where a Von-Mises type plasticity can be used. The Von Mises criterion is formulated 𝑓𝑓(𝝈𝝈) = ‖𝒔𝒔‖ -√2𝑘𝑘 ≤ 0, with 𝒔𝒔 ∈ ℝ 5 the deviatoric stress tensor and 𝑘𝑘, the shear limit which can written as a second-order cone of dimension 6: (√2𝑘𝑘, 𝒔𝒔) ∈ 𝐾𝐾 6 . In the case of limit analysis, the elastic response of the structure is ignored therefore the plastic failure of the structure can be analyzed. In this case, the hessian matrix is removed from the objective thus this problem is reduced to minimizing a linear term under conic constraints. It's impossible to form a classic linear system to solve this problem, even using the incremental Newton-Raphson. Only convex optimization methods, such as interior-point algorithm, can tackle this issue.

EXAMPLE AND RESULTS

A typical beam-column steel assembly is used to illustrate the method. Counter-calculations are made using Abaqus finite elements software which used classical Newton-Raphson method to solve non-linear problems. To eliminate all forms of variations and to focus on the performance of the solvers only, the same series of 6 iteratively refined meshes is used in all the studies. This model is calculated using the IP algorithm and various configurations in Abaqus, corresponding to different modelling choices of contact enforcement: the penalisation (PEN) technique and the augmented Lagrangian (LA) technique. 4 gives a general description of the model. It consists of a HEB200 central column with two IPE360 beams attached over the flanges using welded end plates and bolts. The end plates have a 15𝑚𝑚𝑚𝑚 thickness and 6 M18 bolts are used to attach each beam. Web stiffeners are used to prevent its buckling. For the HEB column, their thickness is 12.7𝑚𝑚𝑚𝑚 to coincide with the IPE's flange thickness, whereas for the IPE beams, 15𝑚𝑚𝑚𝑚 thick stiffeners are used over the supports to correctly channel shear forces. We suppose that the bolt hole is equal to its diameter and, to prevent rigid body motions, one of the bolt heads is glued to the plate. A 45𝑀𝑀𝑀𝑀𝑀𝑀 normal pressure is applied over the top section of the HEB column and the displacement of the nodes belonging to the surface defined by the intersection of the IPE's web stiffeners are blocked in all the 3 directions, the HEB sections remaining free. All details including bolts heads and cores are modeled using 3D quadratic tetrahedrons. Material non-linearity is not considered in this example but can be included easily as discussed in section 4. Since Abaqus offers displacement-based elements only, the calculation times are compared with respect to the kinematic approach. The solved system size is nearly the same as in Abaqus, the system sparsity patterns and conditioning are also comparable. The main difference lies in contact modelling and resolution strategy. To closely compare the numerical performance of the algorithm, all the calculations were made using the same machine. OpenMP technology was used to parallelize over 8 threads. Table 1 shows the CPU times for the different approaches. The IPM is largely comparable to Abaqus penalty approach in terms of CPU times, with even a small speed-up factor when compared to the augmented Lagrangian approach. Compared to the latter, which generally yields more accurate results, the IPM shows a speed-up factor of 2 to 4, the speed-up factor increasing with the mesh size. This is a major advantage of the IPM since it can handle large scale sparse problems efficiently. And while Abaqus number of iterations required for convergence tends to increase (from 20 to 30 for the AL approach), the IPM number of iterations remains fairly constant varying between 18 and 20 iterations for all the remeshes. To assess the global convergence level with mesh refinement, total quantities are considered. The influence of mesh size over the relative gap between the static and kinematic IP approaches has been represented in figures 5 and 6. The quality of convergence of the finite element analysis can therefore be directly evaluated by checking the relative gap between the approaches which reaches acceptable values (less than 5% directly after the first remesh iteration) and 1.68% on the fourth remesh iteration. Let us now consider some local quantities. Figure 7 shows the evolution of the total vertical displacement of the column over the six meshes. We can clearly see the upper and lower bound offered by the dual approach that we used. Abaqus AL approach shows nearly the same results as the kinematic approach whereas with the penalty contact constraints, one can easily overestimate displacements thus leading to misinterpretations. These aspects are generally unknown to engineers who do not have means to estimate displacements from above, contrary to what offers equilibrium-based computations. The maximum gap between the end plate and the column's flange presented in figure 8 shows the consequences of the contact constraints enforcement method. With the penalty approach, the maximum gap is overestimated. An important point is that with the IPM, no contact constraint violation is permitted, thus we have zero penetration, whereas in Abaqus, a penetration tolerance is allowed. The results show that a "negative" gap, equivalent to up to 8% of the maximum gap, is possible with the penalty as seen in figure 9. However, with the AL approach, negative gaps are limited to 0.0001% of the maximum gap, thus yielding better results that are comparable to the IP kinematic approach. It is clear that with the penalty approach in Abaqus and the removal of supplementary contact constraints, one can manage to obtain results in a more reasonable time. Figure 10 shows the local contact pressure over the end plates, one obtained using the IP static approach and the two other from Abaqus. We can clearly see that the stresses postprocessed by Abaqus from a displacement field are underestimated. While the AL approach gives slightly comparable results to the IP static approach (247𝑀𝑀𝑀𝑀𝑀𝑀 compared to 268𝑀𝑀𝑀𝑀𝑀𝑀), the penalty approach largely underestimates these values (89.2𝑀𝑀𝑀𝑀𝑀𝑀). The importance of the dual approach is clear and produces better exploitable results for engineers.

CONCLUSIONS

The study of a deformable continua in contact has been discussed in this article. Using a dual approach, a kinematically admissible displacement field and a statically admissible stress field can be obtained which gives the engineer more reliable results along with an error estimator and a convergence indicator. The IPM handles with no major difficulty the contact nonlinearities with no user-parameter settings. The method can also be extended to consider material non-linearities such as plasticity.

Numerically, the IPMs shows stable and robust results with even better computational times as shown in the steel assembly example. The solution clearly verifies conic constraints, and this can be seen over the contact interfaces where no penetration is permitted. It also provides a better evaluation of quantities of interest for the engineer such as displacements obtained by the kinematic approach and stresses obtained by the static approach.

Figure 1 :

 1 Figure 1 : Local coordinate system and geometry of the friction cones

Figure 2 :

 2 Figure 2 : Reference model In the case of elastic continua with contact interfaces Γ 𝑐𝑐 between solids Ω 1 and Ω 2 , the total potential energy functional can be written as follows: 𝛱𝛱 𝑝𝑝 (𝒖𝒖) = � 1 2 𝜺𝜺: ℂ: 𝜺𝜺 𝑑𝑑𝑑𝑑 𝛺𝛺

Figure 3 :

 3 Figure 3: General idea of an IP algorithm This central path is no other than the unique solution to the following perturbation of the KKT system [1,2,4].

Figure 4 :

 4 Figure 4 : Description of the model

Figure

  Figure 4 gives a general description of the model. It consists of a HEB200 central column with two IPE360 beams attached over the flanges using welded end plates and bolts. The end plates have a 15𝑚𝑚𝑚𝑚 thickness and 6 M18 bolts are used to attach each beam. Web stiffeners are used to prevent its buckling. For the HEB column, their thickness is 12.7𝑚𝑚𝑚𝑚 to coincide with the IPE's flange thickness, whereas for the IPE beams, 15𝑚𝑚𝑚𝑚 thick stiffeners are used over the

Figure 5 :Figure 6 :

 56 Figure 5 : Total potential and complementary energy

Figure 7 :Figure 8 :

 78 Figure 7 : Total vertical displacement of the column

Figure 9 :

 9 Figure 9 : Initial mesh, gap iso-values (in m) over one of the end plates

Figure 10 :

 10 Figure 10 : Fourth remesh iteration, normal pressure iso-values (in MPa) over one of the end-plates.

Table 1 :

 1 CPU times and speed-up factors of the IPM over Abaqus AL and PEN approaches

	Remesh iteration	Mesh size	IP kinematic approach (𝒔𝒔)(𝑵𝑵 𝒊𝒊𝒕𝒕𝒆𝒆𝒊𝒊 )	Abaqus AL approach (𝒔𝒔)(𝑵𝑵 𝒊𝒊𝒕𝒕𝒆𝒆𝒊𝒊 )	Speed-up factor	Abaqus PEN approach (𝒔𝒔)(𝑵𝑵 𝒊𝒊𝒕𝒕𝒆𝒆𝒊𝒊 )	Speed-up factor
	0	𝑁𝑁 𝑒𝑒 = 9150	10.2 (17)	25.0 (21)	2.5	9.0 (8)	0.9
	1	𝑁𝑁 𝑒𝑒 = 33051	28.3 (18)	69.0 (22)	2.4	30.0 (9)	1.1
	2	𝑁𝑁 𝑒𝑒 = 42696	35.0 (17)	97.0 (24)	2.8	38.0 (9)	1.1
	3	𝑁𝑁 𝑒𝑒 = 85187	78.8 (17)	278.0 (25)	3.5	98.0 (9)	1.2
	4	𝑁𝑁 𝑒𝑒 = 222082 320.6 (18) 1368.0 (30)	4.3	435.0 (10)	1.4


Chafik Samir

Ines Adouani

Bézier curves and C 2 interpolation in Riemannian Symmetric Spaces

Keywords: Riemannian Bézier curves, Regression on Riemannian manifolds, Curve fitting, Mean square acceleration, Special orthogonal group

We consider the problem of interpolating a finite set of observations at given time instant. In this paper, we introduce a new method to compute the optimal intermediate control points that define a C 2 interpolating Bézier curve. We prove this concept for interpolating data points belonging to a Riemannian symmetric spaces. The main property of the proposed method is that the control points minimize the mean square acceleration. Moreover, potential applications of fitting smooth paths on Riemannian manifold include applications in robotics, animations, graphics, and medical studies.

Introduction

The problem of constructing smooth interpolating curves in non-linear spaces, or manifolds plays an important role in a wide variety of applications. For instance, interpolation in the rotation group SO(3) has immediate application not only in computer graphics and animation of 3D objects [START_REF] Barr | Smooth interpolation of orientations with angular velocity constraints using quaternions[END_REF][START_REF] Popiel | Bézier curves and c2 interpolation in riemannian manifolds[END_REF], [START_REF] Fang | Real time motion fairing with unit quaternions[END_REF], but also in applications ranging from robot motion planning to machine vision [START_REF] Arnould | Fitting smooth paths on riemannian manifolds: Endometrial surface reconstruction and preoperative mri-based navigation[END_REF][START_REF] Camarinha | Splines of class C k on non-euclidean spaces[END_REF][START_REF] Zefran | On the generation of smooth three-dimensional rigid body motions[END_REF]. Such applications encourage us to further search for some efficient methods to generate smooth interpolating curves on non-linear spaces.

Motivated by potential applications in engineering science and technology, our goal is to develop a new framework for generating C 2 Bézier curves on Riemannian manifolds that interpolates a given ordered set of points at specified time instants. While quite general, we will focus on a special class of Riemannian symmetric spaces. The task of constructing interpolating curve on SO(n) has attracted the attention of several authors. One of the most widely cited approaches is the work of Shoemake [START_REF] Shoemake | Animating rotation with quaternion curves[END_REF] on SO [START_REF] Barr | Smooth interpolation of orientations with angular velocity constraints using quaternions[END_REF], who adopts a re-parametrization of the rotation matrices based on unit quaternion representation. Shoemake's approach can essentially be viewed as a generalization of the de Casteljau's algorithm for Bézier curves to SU (2) in which two elements of SO(3) are interpolated by the geodesic that joins them. Although this algorithm seems computationally efficient, unfortunately the resulting curve depends on the choice of local system coordinates. A few years later, taking into account the Shoemake algorithm, a more careful geometric analysis of unit-quaternion-based method was introduced by Barr et al. [START_REF] Barr | Smooth interpolation of orientations with angular velocity constraints using quaternions[END_REF], Hart et al. [START_REF] Hart | Visualizing quaternion rotation[END_REF], Ge and Ravani [START_REF] Ge | Geometric construction of bézier motions[END_REF], and Nielson et al. [START_REF] Nielson | Animated rotationsusing quaternions and splines on a 4d sphere[END_REF]. Despite the fact of producing an intrinsic curves, these approaches does not generalize to higher-dimensional manifold.

In this paper, we present a novel framework to treat the interpolation problem in the setting of Riemannian geometry and Bézier curve approach. We show that it makes sense to define a C 2 interpolating Bézier curve on Riemannian symmetric spaces as the result of a least squares minimization and a recursive algorithm. In particular, we will focus on a special class of Riemannian symmetric spaces: the special orthogonal group SO(n). Indeed, working in such Riemannian manifold allow nice properties to solve the issues above. The key point to give explicit solution for the interpolation problem and ensures the C 2 differentiablity condition at joint points is the use of global symmetries in these last points. In fact, we will first derive equations for control points of a C 2 Bézier curve on the Euclidean space R m . Then, building upon prior works [START_REF] Arnould | Fitting smooth paths on riemannian manifolds: Endometrial surface reconstruction and preoperative mri-based navigation[END_REF][START_REF] Samir | C 1 interpolating bézier path on riemannian manifolds, with applications to 3d shape space[END_REF], we use these equations to find the control points of a C 1 interpolating Bézier curve on Riemannian manifolds as a generalization of the Bézier based fitting in the Euclidean space and by means of methods of Riemannian geometry. These results are sufficient to give explicit formula for control points of the C 2 interpolating Bézier curve on SO(n). The proposed method will be shown to enjoy a number of nice properties and the solution is unique in many common situations.

The rest of the paper is organized as follows. In section 2, we present our new algorithm to construct a C 2 Bézier curve on the Euclidean space. This will help with the visualization of its main features and motivate its generalization on SO(n). In section 3, the generalization of our approach on the Lie group SO(n) is prescribed. We conclude the paper with numerical examples and a conclusion .

C 2

Interpolating Bézier Curves on R m

In this section, we first describe our approach on the Euclidean space R m . For simplicity we will assume that the time instants are t i = i. In this work, we only use Bézier curves of degree 2 and 3 such that the segment joining p 0 and p 1 , as well as the segment joining p N 1 and p N are Bézier curves of order two, while all the other segments are Bézier curves of order three. Explicitly, the Bézier curve k of degree k 2 {2, 3} are expressed in R m with a number of control points b i , represented as their coefficients in the Bernstein basis polynomials by :

2 (t; b 0 , b 1 , b 2 ) = b 0 (1 t) 2 + 2b 1 (1 t)t + b 2 t 2 , 3 (t; b 0 , b 1 , b 2 , b 3 ) = b 0 (1 t) 3 + 3b 1 t(1 t) 2 + 3b 2 t 2 (1 t) + b 3 t 3 .
Moreover, we assume that there exists two artificial control points ( b b i , b b + i ) on the left and on the right hand side of the interpolation point p i for i = 1, ..., (N 1). Consequently, the Bézier curve on R m is given by:

(t) = 8 > > < > > : 2 (t; p 0 , b b 1 , p 1 ), if t 2 [0, 1] 3 (t (i 1); p i 1 , b b + i 1 , b b i , p i ), if t 2 [i 1, i], i = 2, ..., N 1 2 (t (N 1); p N 1 , b b + N 1 , p N ), if t 2 [N 1, N] Then is C 1 on [t i , t i+1 ], for i = 0, ..N 1.
To ensure that is C 1 at knots p i , for i = 1, ..N 1, we shall make the following assumption:

˙ k i (b i 0 , ..., b i k i ; t i + 1)| t=i = ˙ k i+1 (b i+1 0 , ..., b i+1 k i+1 ; t i)| t=i i = 0, ..., N 2. (1)
This differentiability condition allows us to express b b + i in terms of b b i as:

b b + 1 = 5 3 p 1 2 3 b b 1 , (2)
b b + i = 2p i b b i , i = 2, ..., N 2 (3) b b + N 1 = 5 2 p N 1 3 2 b b N 1 , (4)
We are left with the task of computing the control points b b i , for i = 1, ..., N 1, that generate the C 1 Bézier curve . In [START_REF] Samir | C 1 interpolating bézier path on riemannian manifolds, with applications to 3d shape space[END_REF], we have shown that solutions of the problem of minimization of the mean square acceleration of the Bézier curve are exactly the control points of the curve:

min b b 1 ,..., b b N 1 E( b b 1 , ..., b b N 1 ) := min b b 1 ,..., b b N 1 Z 1 0 k ¨ 0 2 (t; p 0 , b b 1 , p 1 )k 2 + N 2 X i=1 Z 1 0 k ¨ i 3 (t; p i , b b i , b b i+1 , p i+1 )k 2 + Z 1 0 k ¨ N 1 2 (t; p N 1 , b b N 1 , p N )k 2 (5) It turns out that the optimal solution Y = [ b b 1 , ..., b b N 1 ] T 2 R (N 1
)⇥m of ( 5) is the unique solution of a tridiagonal linear system

Y = A 1 CP = DP with j=N +1 X j=0 d ij = 1. ( 6
)
where A is a tridiagonal sparse square matrix of size (N 1) ⇥ (N 1) with a dominant diagonal, C a matrix of size (N 1) ⇥ (N + 1) and P the matrix of p i 's of size (N + 1) ⇥ m given by:

A (1,1:2) = [16 6] (7)
A (2,1:3) = [6 36 9] (8)

A (i,i 1:i+1) = [9 36 9], (9)
A (n 1,n 2:n 1) = [9 36] (10)
C (1,1:2) = [16 6] (11)
C (2,2:3) = [6 36 9] (12)
C (i,i:i+1) = [9 36 9], i = 3, ..., n 2 (13)
C (n 1,n 1:n+1) = [9 36] (14)
Now, let us assume that is C 1 , so that (1) is met and the solution Y given by ( 6) is obtained. The additional C 2 condition for a C 1 curve is the equality of the second derivative at the joint point p i , for i = 1, ..., N 1:

¨ k i (b i 0 , ..., b i k i ; t i + 1)| t=i = ¨ k i+1 (b i+1 0 , ..., b i+1 k i+1 ; t i)| t=i i = 0, ..., N 2.
It is obvious that with this C 2 condition the position of the control points b b i and b b + i that generate the curve will be modified. Therefore, it is more convenient to use another notation. Let us denote by b i and b + i the new control points on the left and on the right hand side of the interpolation point p i , for i = 1, ..., N 1. Computing the acceleration of on respective intervals and taking into account that is C 1 , we shall replace b + 1 by ( 2), b + i by (3), and b + N 1 by (4). We deduce that :

b 2 = 1 3 p 0 1 2 b 1 + 8 3 p 1 , (15a)
b i+1 = b + i 1 + 4p i 4b i , i = 2, ..., N 2 (15b)
p N = 2p N 1 + 2b + N 1 6b N 1 + 3b + N 2 , (15c)
We see at once that points that will be modified by the additional C

C 2

Interpolating Bézier Curves On SO(n)

Our objective in this section is to work out concretely the extension of our approach used to find control points that define a C 2 Bézier curve in the Euclidean space to the Riemannian manifold SO(n). In other words, given R 0 , ..., R N a set of (N + 1) distinct points in SO(n) and 0 = t 0 < t 1 < ... < t N = N an increasing sequence of time instants, we present a conceptually simple framework to construct a C 2 Bézier curve

: [0, N] ! SO(n) such that (t k ) = R k , k = 0, ..., N.
For the most part of Riemannian manifolds, the generalization of our approach is not straightforward. For the case treated here, of the Lie group SO(n), since it is a symmetric space and all the important geometric functions have nice, closed-form expressions, the problem of finding a C 2 Bézier curve that interpolates a given set of points in such space can be completely solved.

Let us start by briefly sketch the differential structure of SO(n). We illustrate this with the geometric toolbox described in table.1. For more details concerning the differential geometry of SO(n), see [START_REF] Helgason | Differential Geometry, Lie Groups, and Symmetric spaces[END_REF], [START_REF] Adams | Lectures on Lie Groups[END_REF].

Set: SO(n) = {R 2 R n⇥n | R T R = In and det(R) = 1} Tangent spaces: TRSO(n) = {H 2 R n⇥n | RH T + HR T = 0} Inner product: < H1, H2 >R= trace(H T 1 H2) Exponential: Exp R (H) = Exp I (R T H) = R exp(R T H). Logarithm: Log R 1 (R2) = R1 log(R T 1 R2)
The shortest geodesic arc joining R 1 to R 2 in SO(n) can be parameterized explicitly by:

↵(t, R 1 , R 2 ) = R 1 exp(t log(R T 1 R 2 )), t 2 [0, 1]. (16)
and we write:

↵(t, R 1 , R 2 ) := @ @u | u=t ↵(t, R 1 , R 2 ).
Furthermore, for each R 1 2 SO(n), there exists a symmetry ) R2 the derivative of the geodesic symmetry ' R1 at R 2 . Then, the following result can be easily proved and will be very important for the derivation of the results presented along this section.

' R1 : SO(n) ! SO(n), R 2 ! R 1 R T 2 R
Lemma 1. Let R 1 2 SO(n). i) (d' R1 ) 1 R2 = (d' R1 ) 'R 1 (R2) , for all R 2 2 SO(n) ii) (dExp R1 ) 1 H = (dExp R1 ) H (d' R1 ) Exp R 1 (H)
, for all H 2 T R1 SO(n) Let us now denote by k (t, V 0 , ..., V k ) the Bézier curve of order k 2 {2, 3} on SO(n) with a number of control points V i for i = 0, ..., k. Furthermore, similar to the Euclidean case, we will suppose that there exists two artificial control points ( b

Z i , b Z + i )
on the left and on the right hand side of the interpolation point R i for i = 1, ..., (N 1). Hence, the Bézier Curve : [0, N] ! SO(n) is defined by:

(t) = 8 > > < > > : 2 (t; R 0 , b Z 1 , R 1 ), if t 2 [0, 1] 3 (t (i 1); R i 1 , b Z + i 1 , b Z i , R i ), if t 2 [i 1, i], i = 2, ..., N 1 2 (t (N 1); R N 1 , b Z + N 1 , R N ), if t 2 [N 1, N]
In order to obtain equations that govern the control points of the C 2 Bézier curve on SO(n), one should begin to compute ( b Z i , b Z + i ), for i = 1, ..., N 1, control points of the Bézier curve that ensure the C 1 diffirentiablity condition of at knots R i on SO(n). To do this, our main idea is to treat the fitting problem on the tangent space T Ri SO(n) at a point R i 2 SO(n) as for the Euclidean case. Consequently, for each i = 1, ..., N 1, we would like to transfer the data R 0 , ..., R N in each tangent space T Ri SO(n) using Riemannian logarithmic map. The mapped data are then given by Q

= (Q i 0 , ..., Q i N ) with Q i k = Log Ri (R k ) for k = 0, ..., N.
Applying our approach used to define a C 2 Bézier curve on the Euclidean space R m in each tangent space T Ri SO(n), for i = 1, ..., N 1, provides a natural and intrinsic method to compute control points ( b

Z i , b Z + i ) of the desired C 1 Bézier curve on SO(n). Theorem 1. Let R 0 , ..., R N be a finite sequence of distinct points in the special orthogonal group SO(n) with R T i R k , i 6 = k, sufficiently close to I n . For each i = 1, ..., N 1, Q = (Q i 0 , ..., Q i N ) are the corresponding mapped data in the tangent space T Ri SO(n) at R i defined by Q i k = Log Ri (R k ) for k = 0, ..., N. Set t 0 = 0 < ... < t N = N a sequence of time instants.
Then, there exists a unique matrix

X i = [(B 1 1 ) , ..., (B 1 N 1 ) ] T 2 R n(N 1
)⇥n containing the (N 1) control points that generate the C 2 Bézier curve i , in each tangent space T Ri SO(n) and a matrix Q = [ Qi 0 , ..., Qi N ] T of size n(N + 1) ⇥ n containing the new (N + 1) interpolation points in each tangent space T Ri SO(n).

Algorithm 1 Construction of the C 1 interpolating Bézier curve on SO(n).

Input: N 3, R = [R0, ..., RN ] T a matrix of size n(N +1)⇥n containing the (N +1) interpolation points on SO(n).

Output: b Z and R.

1: for i = 1 : N 1 do 2: Calculate Q = [Q i 0 , ..., Q i N ]
T a matrix of size n(N + 1) ⇥ n containing the (N + 1) interpolation points on TR i SO(n): 3:

for k = 0 : N do 4:

Q i k = Log R i (R k ) = Ri log(R T i R k ) 5: Calculate Xi = [(B i 1 ) , ..., (B i N 1 )
] T a matrix of size n(N 1) ⇥ n containing the (N 1) control points of the C 2 Bézier curve i on TR i SO(n), and Q = [ Qi 0 , ..., Qi N ] T a matrix of size n(N + 1) ⇥ n containing the new interpolation points on TR i SO(n) using the prescribed method on section 2. 6:

Calculate control point b Z i with b Z i = Exp R i ((B i i ) ) 7:
Calculate the new interpolation points Rk = Exp R i ( Qi k ).

8:

end for 9: end for 10: return b Z and R, Proposition 1. Under the same hypotheses of Theorem 1, there exists a unique matrix

Z = [ b Z 1 , ..., b Z N 1 ] T 2 R n(N 1
)⇥n , containing the (N 1) control points that generate the Bézier curve interpolating the points R i at t i on SO(n), for i = 0, ..., N. The rows of b Z are given by: b

Z i = Exp Ri (x i ), i = 1, ..., N 1. ( 17
)
where xi , represent the row i of X i in T Ri SO(n), for i = 1, .., N 1. Moreover, the new (N + 1) interpolation points in SO(n) are given by:

Rk = Exp Ri ( Qi k ), k = 0, ..., N; i = 1, ..., N 1. (18)
Algorithm 1 provides a detailed exposition of the steps of the proof of Theorem 1 and Proposition 1.

Corollary 1. The Bézier path : [0, 1] ! SO(n) is C 1 on SO(n).
Proof. The following result may be proved in much the same way as Corollary 3.3. in [START_REF] Samir | C 1 interpolating bézier path on riemannian manifolds, with applications to 3d shape space[END_REF].

We are now in a position to formulate the main theorem of this section, which contains the counterpart of the equations derived in the last section that generate control points of a C 2 Bézier curve on R m . Let us assume that is C 1 , so that the solution b Z is obtained. Let us denote by Z i and Z + i the new control points on the left and on the right side of the interpolation point Ri that generate the C 2 Bézier curve on SO(n). The key point to find the control points Z i , for i = 1, ..., N 1 is similar to the Euclidean case. That is, we might know Z 1 (and therefore Z + 1 by the C 1 differentiability condition on SO(n)) and wish to define iteratively Z i for i = 2, ..., N 1 (and obviously Z + i in much the same way as Z + 1 ).

Algorithm 2 Construction of the C 2 interpolating Bézier curve on SO(n).

Input: N 3, R = [ R0, ..., RN ] T a matrix of size n(N +1)⇥n containing the (N +1) interpolation points on SO(n).

Output: Z. 1: Calculate b Z = [ b Z 1 , ..., b Z N 1 ] T using Algorithm 1. 2: Set Z 1 = b Z 1 . 3: Calculate control point Z + 1 : 4: Z + 1 = Exp R1 ( 2 3 Exp 1 R1 (Z 1 )) 5: Calculate control point Z 2 : 6: Z 2 = Exp Z + 1 ⇣ 1 3 ⇣ (d' R1 ) Z 1 ⇣ ↵(1, R0, Z 1 ) ⌘ 4 ↵(0, Z 1 , R1) ⌘⌘ 7: for i = 2 : N 2 do do 8: Z + i = Exp Ri ( Exp 1 R1 (Z i )) 9: Z i+1 = Exp Z + i ⇣⇣ (d' Ri ) Z i ↵(1, Z + i 1 , Z i ) 2 ↵(0, Z i , Ri) ⌘⌘ 10: end for 11: Calculate control point Z + N 1 : 12: Z + N 1 = Exp RN 1 ( 2 3 Exp 1 RN 1 (Z N 1 )) 13: return Z,
Theorem 2. Let R0 , ..., RN be a set of distinct points in the special orthogonal group SO(n) given by Eqn. (18) and ↵(t) the shortest geodesic arc joining control points of the curve on SO(n) given by Eqn. (16). Let X 1 = [(B 1 1 ) , ..., (B 1 N 1 ) ] T be the matrix of size n(N 1) ⇥ n containing the control points of the C 2 Bézier curve 1 in T R1 SO(n). Then, there exists a unique matrix Z = [Z 1 , ..., Z

N 1 ] T 2 R n(N 1)⇥n , containing the (N 1) control points that generate the C 2 Bézier curve interpolating the points Ri at t i on SO(n), for i = 0, ..., N. The rows of Z are given by: i) Z 1 = Exp R1 ((B 1 1 ) ). ii) Z 2 = Exp iii) Z i+1 = Exp

Z + i ⇣⇣ (d' Ri ) Z i ↵(1, Z + i 1 , Z i ) 2 ↵(0, Z i , Ri )
⌘⌘ , i = 2, ..., N 2.

We illustrate the proposed method to construct a smooth interpolating path on SO(3) from four rotation matrices R 1 , R 2 , R 3 , and R 4 . We display the result in Figure 1 where rotations are applied to rotate a 12 sided dice and the given time instants are displayed in a box. We can easily check that the resulting curve path is smooth including at the interpolation points.

Conclusion

In this paper, we have introduced a new framework and algorithms to study the fitting problem of C 2 Bézier curves to a finite set of time-indexed data points on the special orthogonal group SO(n). The proposed method takes into account the global symmetries defined in the joint points. Therefore, the presented approach is valid on any locally symmetric space and other Riemannian symmetric spaces. In the future, we intend to extend the theory and then apply it to more general nonlinear manifolds.

2

 2 condition are b b i and hence b b + i , for i = 2, ..., N 1. The point b b 1 remains invariant and consequently it will be the case for b b + 1 . We thus get b 1 = b b 1 , with b b 1 is the first row of the matrix Y obtained as a solution of the optimization problem (5). However, the endpoint p N is affected as we can deduce from Eqn. (15c). Nevertheless, it follows that giving the control point b 1 allows us to find all the other control points including b 2 with Eqn. (15a) and hence b + 2 with (3), then b i+1 for i = 2, ..., N 2 with (15b) and therefore b + i , for i = 3, ..., N 2 with (3) and b + N 1 with (4).

Fig. 1 :

 1 Fig. 1: Example of an interpolating path on SO(3) applied to rotate a 12 sided dice at given time instants (1, 5, 9, 13).

Table 1 :

 1 Geometric toolbox for the Riemannian manifold SO(n)

  1 that reverses geodesics through R 1 . It is easy to check that ' R1 is an isometry and thus SO(n) turns into a Riemannian symmetric space. For R 1 , R 2 2 SO(n), let us denote by (dExp R1 ) H the derivative of Exp R1 at H 2 T R1 SO(n) and by (d' R1