Emre

Nabil

Jérome Wissem

Keywords: Electrical machines, Meta-models, Adjoint variable method, Uncertainty propagation, Finite element method

This dissertation deals with the approaches for optimization and reliability analysis of electrical machines modelled by the nite element method. The nite element method is the most sophisticated tool to model the electromagnetic phenomenon.

However, it is computationally expensive. Thus, its usage for optimization and reliability analysis (iterative processes) should be made with caution since only a limited number of evaluations of the model can be tolerated. Furthermore, the impact of the manufacturing process on the electrical machines is scarcely studied in the literature.

The integration of this aspect in the design phase is one of the contributions of this thesis alongside the main contribution, which is the development and comparison of optimization approaches for electrical machines.

We present the approaches adapted to the subject and develop new ones. On the one hand, the nite element model can be seen as a "black-box" for which we develop a non-intrusive approach based on Kriging meta-models. On the other hand, we consider an intrusive approach as we look inside the "black-box," we upgrade the model to provide the derivatives of the quantities of interest. The derivatives are essential to some optimization and reliability analysis tools. They are computed eciently using the adjoint variable method. Finally, the methods are compared to give insight into the advantages and the shortcomings of each of them.

Lastly, a real case study is considered; it consists of studying the impact of the manufacturing process on a claw-pole machine manufactured by Valeo. From the production line, machines are withdrawn to measure their dimensions and characterize their deviation from the nominal one. Then a statistical analysis is conducted to assess the reliability and impact on the performances. vii Résumé Ce mémoire aborde les approches d'optimisation et d'analyse de la abilité des machines électriques modélisées par la méthode des éléments nis. La méthode des éléments nis est l'outil le plus sophistiqué pour modéliser le phénomène électromagnétique. Cependant, elle est coûteuse en temps de calcul. Ainsi, son utilisation pour l'optimisation et l'analyse de abilité (processus itératifs) doit être faite avec prudence car seul un nombre limité d'évaluations du modèle peut être toléré.

De plus, l'impact du processus de fabrication sur les machines électriques est peu étudié dans la littérature. L'intégration de cet aspect dans la phase de conception est l'un des apports de cette thèse aux côtés de la principale contribution, qui est le développement et la comparaison des approches d'optimisation pour les machines électriques.

Nous exposons les approches adaptées au sujet de la thèse et en développons de nouvelles. D'une part, le modèle d'éléments nis peut être considéré comme une "boîte noire" pour laquelle nous développons une approche non intrusive basée sur les méta-modèles de Krigeage. D'autre part, nous considérons une approche intrusive car nous regardons à l'intérieur de la "boîte noire", nous améliorons le modèle pour fournir les dérivées des quantités d'intérêt. Les dérivées sont essentielles à certains outils d'optimisation et d'analyse de abilité. Elles sont calculées ecacement en utilisant la méthode de la variable adjointe. Enn, les méthodes sont comparées pour donner un aperçu des avantages et des inconvénients de chacune d'entre elles.

Enn, une étude de cas réel est abordée; elle consiste à étudier l'impact du procédé de fabrication sur la machine à gries fabriquée par Valeo. Sur la chaîne de production, les machines sont prélevées pour mesurer leurs dimensions et caractériser leur écart par rapport aux dimensions nominales. Ensuite, une analyse statistique est menée pour évaluer la abilité et l'impact sur les performances. ix

Résumé substantiel

En résumé, les objectifs de cette thèse ont été atteints. Une comparaison des approches appropriées pour l'optimisation des machines électriques a été eectuée.

Les avantages et les inconvénients de chaque méthode ont été mis en évidence. L'intégration de l'analyse de abilité dans la phase de conception a été évaluée sur un cas test réel (une machine à gries). L'impact du processus de fabrication sur la géométrie de la machine a été étudié, et une analyse quantitative a été réalisée.

Introduction

La méthode des éléments nis est l'outil le plus sophistiqué pour modéliser le phénomène électromagnétique. Cependant, elle est coûteuse en temps de calcul.

Ainsi, son utilisation pour l'optimisation et l'analyse de abilité (processus itératifs) doit être faite avec prudence car seul un nombre limité d'évaluations du modèle peut être toléré. De plus, l'impact du processus de fabrication sur les machines électriques est peu étudié dans la littérature. L'intégration de cet aspect dans la phase de conception est l'un des apports de cette thèse aux côtés de la principale contribution, qui est le développement et la comparaison des approches d'optimisation pour les machines électriques.

Etat de l'art

Ce chapitre examine les algorithmes d'optimisation et la méthode des éléments nis pour mettre en perspective les contributions de cette thèse. Tous les algorithmes d'optimisation ne sont pas adaptés aux modèles à base d'éléments nis en raison du coût de calcul de ces modèles. Nous présentons les principales catégories d'algorithmes d'optimisation et leurs extensions au contexte de l'optimisation sous incertitudes. Une brève présentation de la modélisation électromagnétique des machines électriques utilisant la méthode des éléments nis est également présentée.

Deux approches adaptées au sujet de cette thèse sont sélectionnées pour être investiguées.

Méthode non-intrusive

Ce chapitre présente une approche non-intrusive dans laquelle le modèle éléments nis est considéré comme une boîte noire, et un modèle de substitution est créé en évaluant le modèle coûteux en quelques points, puis en utilisant le modèle de substitution peu coûteux pour l'optimisation et la quantication des incertitudes. Nous présentons une introduction approfondie de l'optimisation à l'aide de méta-modèles xi et des dés posés par la mise en place de "bons" méta-modèles et la manière dont on peut évaluer leur qualité, puis nous proposons une nouvelle stratégie d'exploitation des modèles de krigeage.

Méthode intrusive

Ce chapitre examine une approche intrusive pour calculer les dérivés d'un modèle à éléments nis. Dans cette approche, la majeure partie du travail est eectuée sur le modèle plutôt que sur l'algorithme. Les dérivées sont utilisées par les algorithmes d'optimisation et les méthodes de quantication de l'incertitude les plus ecaces.

Nous présentons la méthode de la variable d'état adjointe et la façon dont elle peut être obtenue à partir d'un code éléments nis. La méthode de la variable adjointe permet de calculer ecacement les dérivées à partir du code d'éléments nis.

Benchmark et comparaisons

Ce chapitre est consacré aux tests numériques des approches développées dans le cadre de l'optimisation. Pour ce faire, nous abordons deux benchmarks bien connus traités par les chercheurs de la communauté électromagnétique. Il s'agit des cas tests issus des workshops TEAM (Testing Electromagnetic Analysis Methods) [1] [2].

Nous présentons une comparaison entre les approches étudiées et d'autres approches pour l'optimisation des dispositifs électromagnétiques.

Machine à gries

Ce chapitre traite l'impact des imperfections du processus de fabrication sur les performances de la machine à gries. Nous étudions principalement les imperfections géométriques. Dans ce chapitre, nous présentons la machine à gries et son mode de fonctionnement puis nous introduisons la procédure de métrologie adoptée pour les mesures métrologiques de la machine électrique. Ces mesures sont analysées, et la variabilité de la géométrie de la machine électrique due au processus de fabrication est modélisée. Ensuite, un modèle éléments nis paramétré de la machine est présenté. Puis, la propagation de l'incertitude utilisant la variabilité caractérisée et le modèle de la machine est réalisée en utilisant la simulation de Monte Carlo avec le modèle de Krige.

Conclusion

Cette thèse traitait des approches d'optimisation et d'analyse de la abilité des machines électriques modélisées par la méthode des éléments nis. La méthode des éléments nis est l'outil le plus sophistiqué pour modéliser les phénomènes électromagnétiques dans les machines électriques. Elle permet de modéliser nement les champs électromagnétiques dans le domaine étudié et de manipuler des géométries complexes, comme la machine à gries. Cependant, elle est coûteuse en temps de calcul à cause des matériaux non linéaires, des géométries 3D, de la dépendance au temps. Par conséquent, son utilisation dans le cadre de l'optimisation et de l'analyse de abilité (processus itératifs) doit être faite avec précaution, car seul un nombre limité d'évaluations de l'outil de simulation peut être toléré.

Les algorithmes d'optimisation générale ne peuvent pas être appliqués tels quels ; certaines modications indispensables sont nécessaires pour les rendre utilisables pour l'optimisation des machines électriques. Dans le premier chapitre, nous avons présenté la littérature sur les algorithmes d'optimisation, et nous avons discuté de leur adaptabilité au sujet de cette thèse. À partir des diérentes catégories, deux approches ont été envisagées, qui dépendent des capacités du modèle utilisé.

1. Une approche non-intrusive basée sur le méta-modèle de Krigeage 2. Une approche intrusive basée sur la méthode de la variable adjointe L'approche non-intrusive utilisant des modèles de substitution est largement utilisée dans le contexte de l'optimisation de modèles coûteux. De nombreux chercheurs ont souligné les dicultés que pose l'utilisation d'une telle stratégie. Dans cette thèse, nous avons mis l'accent sur certains des problèmes très connus et proposé une nouvelle méthodologie d'utilisation des méta-modèles et la manière d'accélérer les temps d'optimisation. Nous avons souligné les inconvénients de l'approche conventionnelle, qui consiste à faire tenir un seul méta-modèle sur l'ensemble de l'espace de conception et à l'enrichir séquentiellement à l'aide de critères de remplissage. Nous avons proposé une nouvelle stratégie qui consiste à construire de nombreux métamodèles sur des régions spéciques de l'espace de recherche Ensuite, de manière itérative, on élague les régions qui ne sont pas prometteuses. Ce processus explore tout l'espace de recherche ; il nous permet donc de produire des solutions ables. D'autre part, pour l'approche intrusive, nous proposons comment calculer les dérivés d'un modèle d'éléments nis. La plupart des travaux ont été eectués sur le modèle plutôt que sur l'algorithme (nous avons utilisé l'algorithme SQP pour cette thèse). Nous avons présenté la méthode de la variable adjointe et la manière dont le gradient peut être dérivé d'un code FEM. Nous avons développé une méthode ecace pour calculer les dérivées des sensibilités de forme pour les paramètres géométriques, qui sont essentiels pour le calcul du gradient. La méthode de la variable adjointe a xiii été comparée à la méthode des diérences nies pour valider et mettre en évidence son ecacité en termes de précision et de temps de calcul.

Pour la comparaison entre des approches intrusives et non-intrusives pour l'optimisation des dispositifs électromagnétiques à l'aide de la méthode des éléments nis, nous avons traité deux problèmes de référence bien connus de la littérature connus sous le nom de TEAM Workshop : les problèmes 22 et 25 [1] [2]. Nous avons utilisé deux métriques pour la comparaison, la première est la qualité de la solution, et la seconde est le coût de calcul.

En termes de performances, l'approche intrusive surpasse les autres stratégies pour les deux cas de test ; cela a été possible grâce au calcul du gradient en utilisant la méthode de la variable adjointe. Cela améliore considérablement la qualité des solutions mais s'accompagne d'un coût de manipulation intrusive du code de la méthode des éléments nis. L'approche non-intrusive demeure une excellente alternative en termes de mise en ÷uvre et de qualité des solutions. L'approche développée a permis de surmonter certains des problèmes très connus qui se posent lors de l'utilisation de méta-modèles pour l'optimisation.

Dans le dernier chapitre, nous avons présenté une démarche complète pour effectuer une analyse de abilité de la machine électrique ; nous avons commencé par la présentation de la machine électrique. Ensuite, nous avons présenté les diérentes mesures métrologiques qui ont été eectuées pour caractériser la variabilité de sa forme. Puis, nous avons présenté un modèle éléments nis paramétrique de la machine et nous avons exposé comment la mesure peut être prise en compte dans le modèle numérique. L'évaluation numérique de l'impact des variables géométriques sur les performances de la machine a été réalisée en utilisant une approche de méta-modèle pour réduire le coût de calcul. La variabilité des performances a ensuite été caractérisée, et les variables qui inuencent cette variabilité ont également été identiées.

Introduction

The societal context Electrication is gaining momentum. If the early twentieth century saw electricity revolutionized our daily lives, especially in static applications, the beginning of the twenty-rst century sees electricity also profoundly modify all embedded applications. The Electricity Fairy is a very convenient way of transporting and converting energy. The tremendous advances in power electronics has led to new solutions for the vast majority of energy conversion systems. Not to mention, for example, the signicant steps are taken with all-electric cars, e.g. Tesla, or battery-powered aircrafts, e.g. Airbus E-FAN, making them accessible for mass usage. Nevertheless, they are still currently limited by the storage and charging capacities of much less energy than oil. Electricity spreads actively inside systems throughout the energy conversion chain of embedded applications.

The car is becoming hybrid; an electric traction chain assists the internal combustion engine. Naval propulsion is increasingly done by pods containing electrical machines powered by diesel generators. Hydraulic and pneumatic energies are being progressively replaced in aircraft by electric power, causing high expansion of the power of the alternators driven by the engines of the plane (1.4MVA for the Boeing 787 against 300kVA for A330).

The gains are signicant in terms of simplicity of installation, modularity, maintenance, and source of failure (an electric cable never leaks). Electric actuators multiply for the sole purpose of comfort: adjusting the electric seats, the electric window of a car, or motorization of rear-view mirrors.

In all these more electrical systems, the ends of the electromechanical chain are often electrical machines; either an actuator that generates electricity while driven by an internal combustion engine or an actuator to generate mechanical energy that is used to achieve the nal function. These electrical machines are mostly rotating machines. In this context, their improvement is a vital issue in terms of weight, bulk, eciency, manufacturing cost, reliability and recyclability. Furthermore, air pollution has become an urgent problem of environmental sanitation, aecting both developed and developing countries. Increasing amounts of potentially harmful gases and particles are released into the atmosphere, damaging INTRODUCTION human health and the environment. They also undermine the long-term resources needed for the sustainable development of our planet.

Studies in France and around the world have shown that the share of transport in total emissions of air pollutants has increased steadily over the last few years. Some countries have taken initiatives to promote the use of clean or low-emission vehicles.

The reduction of emissions requires a signicant eort in improving the eciency of all energy consumers in vehicles and especially auxiliaries such as air conditioning and alternator. The improvement of the performances of electrical devices in cars, which this thesis examines, is in line with this spirit.

The industrial context

The automotive market is very competitive, which drives automakers to improve their products continually. As a result, electricity and electronics are of great importance in today's automobiles. The electrication of the passenger compartment enhances the comfort and the pleasure of driving and riding the vehicles. Besides, cars are becoming safer thanks to new electrical devices that improve active safety (anti-lock braking system) and passive safety (airbags). Also, the replacement of mechanical or hydraulic devices, such as power steering, by electric actuators can reduce costs and be more comfortable to control.

The increasing demand for electricity in cars is due to the use of these electronic devices. Thus requiring a more signicant on-board production of this energy from the alternator on which the constraints necessary for the functioning of the electrical systems, mainly constraints of eciency, compactness and low-cost manufacturing, will be exercised.

In this context, the company Valeo has performed several studies on claw-pole alternators (Figure 1). The overall goal of these studies was to design and optimize the machine to improve specic performances (Power, noise, . . . ) while ensuring a low cost and compactness. For the designer to reach these goals, it is necessary to have models that address most of the machine's aspects (electrical, thermal, magnetic, mechanical). Nevertheless, these models are generally expensive in computational time, requiring the solution of large system of algebraic equations. Thus, optimizing the machine with these models is cumbersome and takes a tremendous amount of time.

Additionally, most of these studies did not include uncertainties introduced by the manufacturing process. In the design phase, the variability of the electrical machine geometry and material properties has not been taken into account. This variability can impact electrical machine performances and can lead to non-compliant products. Thus, optimization and uncertainty quantication methodologies are of

The complexity of the design

With the ever-growing complexity of products, modelling and simulation allow us to create an all-digital prototype, to understand and optimize the critical performances in order to ensure that the product will fulll the specications correctly during its life cycle. Therefore, the usage of high-delity models has become mandatory.

A model is regarded as representative of a system or a phenomenon. It is a ctitious system, which assembles equations associated with particular physical hypotheses to draw specic conclusions. A model is oriented; from input variables, it provides a result. To estimate the eect of the causes, in the same way, the inverse model is the one that reverses this causality.

A model is ideal for design when its inversion is unique; if a single cause produces the desired eect. Indeed, the reversal makes the development faster and less cumbersome since the work is done only once. Nevertheless, the models are generally not invertible, and there are many degrees of freedom that allow meeting the same specications. Optimization tools can be used to improve a design whose model is not invertible.

Using optimization tools early in the design phase of the development cycle dramatically increases the device potential for improvement. If the designer uses the optimization tools downstream in the development cycle, the main choices will already be made, and the contribution of the optimization will remain moderate. This approach is still not widespread at Valeo, where optimization is used later in the development cycle. Thus, there is a real need for design methodologies to deal with the complexity of the models.

THE COMPLEXITY OF THE DESIGN INTRODUCTION

The computational cost To be able to model the physics related to the studied phenomena, numerical methods are used; one of the famous and the most robust ones is the Finite Element Method, i.e. FEM. It is a method that enables us to determine an approximate solution on a spatial domain by calculating a eld (of scalars, vectors) corresponding to the solution of the given equation.

The method consists of splitting the space domain into small elements, also called meshing, and looking for a simplied formulation of the problem on each element, transforming the system of any equations into a system of equations. The latter is represented by a matrix. The matrices for all the elements are assembled, forming a large matrix; the resolution of this global system gives the approximate solution to the problem. However, the time needed for a FEM simulation depends on the number of small elements. Figure 2 shows the evolution of the computational time versus the number of elements that compose a studied domain. It is worth noting that the data are retrieved by running a nite element code that we have implemented, simulating an electromagnetic device for dierent element sizes. We can see that the computational time evolves polynomially to the number of elements. The claw-pole machine has a very complex geometry; consequently, a signicant number of elements are needed (around two million) to be able to model its behaviour nely. Therefore, the time needed for the simulations is considerable. The existence of nonlinearity or a dynamic behaviour in the physical phenomena can multiply the computational time considerably. The designer should take this feature into account since it will signicantly impact the duration of the design phase by using methodologies that are well-suited.

INTRODUCTION

The manufacturing process

The manufacturing processes, because of their imperfections, lead to dispersions on the characteristics of the products. In the case of electrical machines, these imperfections can impact its geometry and properties. These dispersions induce "variabilities" on quantities of interest and the performances for the product. For example, it becomes crucial for an automobile alternator to study the dispersion inuence on the output current and eciency. This kind of study aims to ensure that the product fullls the specications contracted with the customer.

The classical approach based on a deterministic model (in the sense that the model entries are entirely known) to address this problem is no longer sucient. It is necessary to consider the quantities of inputs as uncertain quantities and to place oneself within the conceptual framework of uncertainties.

The uncertainty in the manufacturing process can be reduced by improving the production lines; however, this will induce a high cost. One approach consists of knowing the uncertainty on the geometry on the material process and predict the impact of the uncertainty on the performances of the product. This approach will enable us to control the inuence of the variability of the geometry or the material characteristics to see if the design variables respect the specications and if their variability is not too high. If it is the case, this approach will also enable to determine which parameters on the geometry or the material behaviour have the most impact on the variability of performances, therefore to act on these parameters in order to have a more robust product or to act on the manufacturing process in order to reduce the variability of these inuential parameters.

Thesis organization

The remainder of this dissertation is organized as follows.

Chapter 1 reviews optimization algorithms and the nite element method to put the contributions of this thesis into perspective. Not all optimization algorithms are adapted for nite element models because of the computational cost of these models.

We present the main categories of optimization algorithms and their extensions to the context of optimization under uncertainty. A brief presentation electromagnetic modelling of the electrical machines using the nite element method is also presented.

Two approaches adapted to the subject of the thesis are selected to be investigated.

Chapter 2 deals with a non-intrusive approach where the nite element model is considered as a black-box, and a surrogate model is t by evaluating the expensive model at few points, then, use the cheap surrogate model for optimization and uncertainty quantication. We present a thorough introduction of optimization INTRODUCTION using meta-models and some of the challenges of tting "good" meta-models and how to assess their quality, then we propose a novel strategy for exploiting Kriging models.

Chapter 3 examines an intrusive approach to compute the derivatives from a nite element model. In this approach, most of the work is done on the model rather than the algorithm. The derivatives are used by the most ecient optimization algorithms and uncertainty quantication methods. We present the adjoint state variable method and how it can be derived from a nite element code. The adjoint variable method enables to compute the derivatives eciently from the nite element code.

Chapter 4 is dedicated to numerical tests of the developed approaches in the context of optimization. For this task, we address two well-known benchmarks treated by electromagnetic community researchers. They are ones of the TEAM (Testing Electromagnetic Analysis Methods) workshops [1] [2]. We present a comparison between the investigated approaches and other approaches for optimizing electromagnetic devices.

Chapter 5 presents the impact of manufacturing process imperfections on the performance of the claw-pole machine. We mainly study geometry imperfections.

In this chapter, we present the claw-pole machine and its operating mode then introduce the metrology procedure adopted for the measurement of the electrical machine. These measurements are analyzed, and the variability of the geometry of the electrical machine induced by the manufacturing process is modelled. After, a parametrized nite element model of the machine is presented. Next, the uncertainty propagation using the variability characterized and the model of the machine is conducted employing Monte Carlo simulation assisted by Kriging meta-model.

The nal conclusions and considerations for future work are given in the last chapter.

Chapter 1

Literature Review

As described in the introduction, one can use optimization tools and models to improve electrical machines' design. One can categorize their use by the type and how they are integrated into the design phase. This chapter begins by briey describing several deterministic optimization algorithms. The second section reviews how to incorporate uncertainty in an optimization process. The third section focuses on the modelling of electrical machines modelled by nite element method. The review of optimization tools and nite element method is presented so that the contributions of this thesis are put into perspective.

Deterministic Optimization

Optimization and especially numerical optimization have had a signicant boost in the last fty years with computers' advent. Optimization is a design tool that helps designers identify the optimal solutions from several possible choices or an innite set of choices. Optimization is increasingly applied in industry since it provides engineers with a reasonable and exible means to identify optimal designs before physical deployment. It is often the last step of computational analysis. First, one should start by studying the physical phenomenon, describing it by equations, studying these equations, and having shown that they can be solved with a computer.

Next, identifying the objective, the variables, and the constraints that should be satised. Then, optimize the system by adjusting the variables to change the solution in the desired direction.

An objective is a measure of the performance of the system that we want to minimize or to maximize. In electrical machines, we may want to maximize the performances or minimize production cost, or both.

The variables are the parameters of the system for which we want to nd values. The variables may be the parameters that dene the geometry or material properties of the electrical machine.

The constraints are the functions that dene the allowable values for the variables and the requirements that should be satised for the system to be acceptable.

Mathematically, one can express any optimization problem as

min p f (p) s.t. g i (p) ≤ 0, i = 1 . . . n c (1.1)
where p = (p 1 , p 2 , . . . , p nv ) are the variables to be determined, f (p) is the objective function to be minimized and g i (p) are n c constraints to be satised. By convention, the standard form denes a minimization problem. A maximization problem can be treated by negating the objective function. Likewise for the constraints, in the case of g(p) ≥ 0 then, it is transformed to -g(p) ≤ 0 g(p) ≤ r where r is a real number, then, it becomes g(p) -r ≤ 0 g(p) = 0 can be transformed to two constraints g(p) ≤ 0 and g(p) ≥ 0.

However, some algorithms handle equality constraints naturally.

the problem can be written in the standard form. The optimization problem should be carefully dened based on the performances required and the corresponding constraints that are given by a model. the constraints at this variable. Then, the optimization algorithm uses this information to approach the optimal solution p opt .

Model

In this section, we present a brief overview of the formulation of optimization problems and the algorithms that enable us to solve them with insight into the pros and cons of each of them.

Optimization problems formulation

As noted before, an essential step in the optimization process is classifying the optimization problem, since algorithms for solving optimization problems are dedicated to a particular type of problems. Indeed, most algorithms developed are designed to solve specic types of problems and are not very ecient for all dierent types.

Here, we provide insights on how optimization problems are classied.

A One vs many objectives

A single objective function denes single-objective problems. Multi-objective problems exist when a compromise is to be sought between several conicting objectives.

For example, optimizing an electrical machine might involve minimizing bulk while maximizing performances. In practice, problems with multiple objectives may be reformulated as single-objective problems by either applying a weighting coecient on the dierent objectives or by replacing some of the objectives with constraints.

Nevertheless, there exist some dedicated algorithms that solve multi-objective problems.

B Unconstrained vs constrained problems

It is vital to identify problems where constraints exist on the variables. Constrained optimization problems result from applications in which there are constraints on the variables. The constraints on the variables can vary from simple bounds or systems of equalities and inequalities to a more complex nonlinear relationship among the variables. Constrained optimization problems can be reformulated and expressed as unconstrained ones by adding a penalty term in the objective function. Naturally, problems with constraints are more complicated to solve and use dedicated algorithms.

C Continuous vs discrete problems

In some cases, the variables are discrete, most often in the form of integers or binaries. The optimization problem is said to be discrete. Thus, several algorithms designed to solve this type of problems have been developed that do not require any hypothesis on the function. The use of black-box models, as well as the application of suited algorithms, made it possible to solve challenging optimization problems. Nevertheless, there still some challenging aspects that make black-box optimization a thriving research area.

F Intuitions on the classications of optimization problems

It is worth noting that there are no problems that are easier to solve than others;

for example, there are some continuous linear problems that are more dicult to solve than nonlinear problems.

Nevertheless, there are some problems for which we can guarantee the optimality of the solution, e.g. quadratic problems; that is to say, the solution found is a global one. Generally, we approach hard problems with "many" of these guaranteed optimality problems, such as approaching nonlinear problems with a series of quadratic problems. These approximations generally lead to acceptable solutions.

Optimization algorithms

The objective of this section is to present dierent types of optimization algorithms. We will present ve classes of algorithms: derivative-based, heuristic-based, population-based algorithms, exact algorithms and approximation-based algorithms.

In every class, we can nd many algorithms that exploit dierent ideas.

To simplify the presentation of algorithms, we focus on the presentation of single

H(p) = [H i,j (p)] = ∂ 2 f ∂p i ∂p j ij (1.4)
It should be noted that whenever the derivatives are available and are easy to compute, a derivative-based method should be preferred. A typical algorithm gradually improves the estimate of the solution step by step. It needs to decide at each iteration the location of the new estimate. However, it has only very local information about the landscape of the objective function the current function value and the derivatives. The derivatives express the slope of the function at a point. Hence, it seems natural that the derivative information is used to dene a search direction, and the algorithm searches along this direction for a new (better) point. Derivative-based algorithms are summarized in the algorithm below.

The designer chooses a starting point p 0 at iteration k = 0, the algorithm exploits the derivatives to nd the "best" direction d 0 and the step α 0 . Consequently, the next iterate p 1 is found, and the algorithm repeats the same process until con- end What are the convergence criteria that enable to stop the algorithm?

The direction d k has to be a descent direction, that is to say, the function f is decreasing in the direction d k , mathematically, it is expressed as

f (p k + εd k ) < f (p k )
for a suciently small ε or equivalently d k T ∇f (p k ) < 0.

The step length α k is also a key ingredient for the success of the algorithm. In the literature, it is referred as line search. The line search consists of determining the step α k to be taken along the descent direction d k . It is often chosen in order to verify the following two conditions:

1. the function f must decrease suciently along the descent direction; 2. the step α k should not be too small in order to avoid slowing down the convergence.

To nd the optimal step length (neither too small nor too big) an exact line search can be performed by minimizing the following problem.

min α≥0 f (p k + αd k ) (1.5)
However, this type of line search generally requires a lot of computation time without signicantly improving the convergence of the algorithm. The inexact line search provides an ecient way of computing an acceptable step length α k that reduces the objective function suciently rather than minimizing the objective function exactly [3]. One of the most famous inexact line search methods is the Wolfe conditions [4].

A step length α k is said to satisfy the Armijo and Wolfe conditions if the following two inequalities hold:

1. f (p k + α k d k ) ≤ f (p k ) + c 1 α k d k T ∇f (p k ) CHAPTER 1. LITERATURE REVIEW 2. -d k T ∇f (p k + α k d k ) ≤ -c 2 d k T ∇f (p k ) with 0 < c 1 < c 2 < 1. c 1 is
usually chosen to be quite small, e.g. 10 -4 , while c 2 is much larger, e.g. 0.9.

Both conditions can be interpreted as providing an upper and lower bound on the admissible step length values. Iterative algorithms are used to solve the inequalities; usually, the rst value of α k that veries both inequalities is chosen.

All line search methods are carried out along descent directions. Dierent ways to determine these directions are discussed: the gradient descent, the conjugate gradient method, Newton's method and its generalization, quasi-Newton methods.

Gradient descent

This method uses the opposite of the function gradient as a descent direction

d k = -∇f (p k ) (1.6)
This direction is obviously a descent direction because we have

d k T ∇f (p k ) = -||∇f (p k )|| 2 < 0
More importantly, it is the steepest descent direction at position p k . The gradient descent method can take many iterations to compute a local minimum, and the line search can also be very long. Other methods, such as conjugate gradient, can be good alternatives. Generally, such methods converge in fewer iterations, but the cost of each iteration is higher [5].

Conjugate gradient

Conjugate gradient construct the descent direction by adding the opposite of the gradient to the previous descent direction multiplied by a real scalar β k

d k = -∇f (p k ) + β k d k-1 (1.7)
At the rst iteration, there is no previous direction then the algorithm takes simply the opposite of the gradient [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. Several methods exist for calculating the term β k+1 , among them the Polack-Ribière method emerges [START_REF] Fletcher | Function minimization by conjugate gradients[END_REF].

β k = ∇f (p k ) T (∇f (p k ) -∇f (p k-1 )) ∇f (p k-1 ) T ∇f (p k-1 ) (1.8)
An additional issue arises for this method to ensure that d k is a descent direction, i.e.

d k T ∇f (p k ) < 0 d k T ∇f (p k ) = -||∇f (p k )|| 2 + βd k-1 T ∇f (p k ) (1.9) Equation (1.9) is negative if d k-1 T ∇f (p k ) = 0.

Newton method

Newton method constructs the descent direction based on a second-order Taylor approximation of around the iterates [START_REF] Avriel | Nonlinear programming: analysis and methods[END_REF]. The second Taylor expansion of f around p k

f (p k + d k ) ≈ f (p k ) + d k T ∇f (p k ) + 1 2 d k T H(p k )d k (1.10)
The minimum of the approximation can be found by setting the rst derivative to zero.

0 ≈ ∇f (p k ) + H(p k )d k (1.11)
Then, the Newton descent direction can be deduced

d k = -H(p k ) -1 ∇f (p k ) (1.12)
This direction is a descent direction if H(p k ) is a denite positive matrix.

Newton's method can guarantee convergence to local optima under mild assumptions, e.g. convexity, positive deniteness, . . . . For the case where the hessian H of the function is not positive denite, the algorithm may diverge. Thus, another category of algorithms was developed to cope with this drawback, i.e. quasi-Newton methods.

Quasi-Newton method

The quasi-newton method has the same iterates as Newton's method except for the hessian, which is not the actual hessian of the objective function but an approximation of it [START_REF] Dennis | Quasi-newton methods, motivation and theory[END_REF]. This method is ecient when the actual hessian is not denite positive, or it is not very easy to compute.

d k = -B k -1 ∇f (p k ) (1.13)
where B k is an approximation of the hessian. There are dierent methods to approximate the hessian the most famous ones is BFGS for Broyden -Fletcher -Goldfarb -Shanno [START_REF] Avriel | Nonlinear programming: analysis and methods[END_REF]. The approximation is written as

B k+1 = B k + y k y k T y k T s k - B k s k s k T B k s k T B k s k (1.14)
where y k = ∇f (p k+1 ) -∇f (p k ) and s k = p k+1 -p k . It has been proved that B k+1 is denite positive if B k is of the same nature and y k T s k > 0. Generally, B 0 is set equal to the identity matrix.

Extension to constrained problems

In constrained optimization, the general aim is to transform the problem into an easier subproblem that can then be solved and used as the basis of an iterative process. The general idea relies on the solution of the Karush-Kuhn-Tucker (KKT) conditions [10][11]. The KKT enables the transformation of the constrained problem to an unconstrained one by including KKT multipliers (equivalent to the Lagrange multipliers for equality constraints).

To recall, the constrained problem is written as

min p f (p) s.t.g i (p) ≤ 0, i = 1 . . . n c (1.15)
Using KKT, it can be transformed

min p,λ L(p, λ) = f (p) + nc i=1 λ i g i (p) (1.16)
where L is called the Lagrangian function and λ are the KKT multipliers.

The KKT equations are necessary conditions for the optimality of constrained optimization problems and stated as follows.

1. Stationarity ∇L(p, λ) = 0

2. Feasibility λ i ≥ 0 and g i (p) ≤ 0 for all i = 1, . . . , n c 3. Complementary slackness λ i g i (p) = 0

The optimal solution should satisfy these conditions. These conditions are necessary but not always sucient, and additional information is required. The necessary conditions are sucient if the objective function and the constraints are convex functions. Thus, for solving complex nonlinear problems, many algorithms solve successive "simple" problems for which the KKT condition are sucient to prove optimality. From which, we will detail sequential quadratic programming (SQP), one of the most robust derivative-based algorithms for nonlinearly constrained optimization problems.

SQP The SQP method can be viewed as a generalization of quasi-Newton's method for constrained optimization [12][13]. It nds a step away from the current point by minimizing a quadratic model of the problem.

min d k f (p k ) + d k T ∇f (p k ) + 1 2 d k T B k d k s.t.g i (p k ) + d k T ∇g i (p k ) ≤ 0, i = 1 . . . n c (1.17)
Dierent algorithms can solve this problem; one of the most used is the active-set.

The active-set algorithm solves the KKT conditions using "guess and check" to nd critical points. Guessing that each inequality constraint is inactive (g i (p k ) + d k T ∇g i (p k ) < 0) is conventionally the rst step. After solving the remaining system for d k , feasibility can be checked. If any constraints are violated, they should be considered active (g i (p k ) + d k T ∇g i (p k ) = 0) in the next iteration. If any multipliers λ i are found to be negative, their constraints should be considered inactive in the next iteration.

The solution to the problem gives the descent direction for the iterate.

d k = -B k -1 [∇f (p k ) + nc i=1 λ i ∇g i (p)] (1.18)
The multipliers λ i for inactive constraints are taken to be zero.

The choice of distance α k to move along the direction generated by the subproblem is not as evident as in the unconstrained case, where we simply choose a step length that approximately minimizes f along the search direction. For constrained problems, we would like to have the next iterate not only to decrease f but also to satisfy the constraints. Often these two aims conict, so it is necessary to weigh their relative importance and dene a merit or penalty function, which we can use to determine whether or not one point is better than another.

M (p k ) = f (p k ) + nc i=1 r k i max(0, g i (p k )) (1.19)
where r k i are suciently positive parameters.

This merit function is used to nd the step length α k to take in the descent direction d k found in (1.18) and, thus, nding the next iterate.

p k+1 = p k + α k d k (1.20)
There exist many algorithms for nonlinearly constrained problems such as interiorpoint and trust-region algorithms. There is not one general algorithm that will work eectively for every kind of nonlinear problems. We chose to detail SQP because it proved to be robust for the cases that we treated.

Extension to multi-objective problems

The ideas swing around transforming a multi-objective problem into many single objective ones. We present two methods: weighted-sum and -constraint. weighted-sum is the most used method; it transforms the objective functions into a weighted sum of them [START_REF] Marler | The weighted sum method for multi-objective optimization: new insights[END_REF]. Then, to use a single objective optimization algorithm. Changing the weights can lead to dierent solutions (but not necessarily) that form the Pareto front. This approach is straightforward to conduct, but it cannot lead to a complete description of the Pareto front if the latter is not convex.

-constraint method only considers one objective function while subjecting the others to constraints [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identication and system optimization[END_REF]. A Pareto front can be obtained by varying the parameters of the constraints.

-constraint method tackles the drawback of the weighted-sum approach.

Extension to black-box problems

The gradient is a fundamental ingredient for derivative-based algorithms; nevertheless, this information is not always available when using black-box models. Finite dierence (FD) can be used to estimate the derivatives using the formula below

∂f ∂p i (p) ≈ f (p 1 , . . . , p i + ε, . . . , p nv ) -f (p 1 , . . . , p i , . . . , p nv ) ε (1.21)
where ε is the nite dierence parameter. This equation is called the forward nite dierence; other approximations also exist, such as backward FD and centred FD.

They exploit the same idea as forward FD nevertheless centred FD may lead to better results than the others.

The choice of ε is essential for the success of the approximation. A very small value will lead to signicant numerical roundo errors, and high value will not lead to the correct approximation of the gradient. Furthermore, the objective function and/or the constraints can be noisy, which renders the derivatives' computation using nite dierence not reliable.

B Derivative-free algorithms

A heuristic algorithm is one that is intended to solve a problem in a faster and more ecient fashion than traditional methods by sacricing optimality, accuracy or precision for speed. From these, we can cite Nelder-Mead simplex, Pattern search, DIRECT (DIviding RECTangle), . . . Engineers and scientists widely use heuristic algorithms. They are a class of optimization methods that are easy to implement, do not require derivatives and are often claimed to be robust to deal with problems with discontinuities or where the functions are noisy.

Nelder-Mead

The Nelder-Mead (NM) is a heuristic numerical method that exploits the concept of a simplex [START_REF] Nelder | A simplex method for function minimization[END_REF]. A simplex can be a line segment on a line, a triangle on a plane, a tetrahedron in three-dimensional space, etc. Starting initially from such a simplex, the simplex undergoes simple transformations during iterations: it deforms, moves and progressively reduces until its vertices approach a point where the function is locally minimal. The full description is shown in the algorithm below.

Algorithm 2: Nelder-Mead algorithm set the initial simplex ;

while Not terminated do order according to the values at the vertices: The Nelder-Mead method frequently gives notable improvements in the rst few iterations and quickly produces satisfying results. Also, it typically requires only one or two objective function evaluations per iteration. This feature is needed in applications where each function evaluation is time-consuming. Aside from being simple to understand and use, this is the main reason for its reputation in practice.

f (p 1 ) ≤ f (p 2 ) ≤ • • • ≤ f (p nv+1
Contrarily, the lack of convergence theory is often considered in practice as a numerical breakdown of the algorithm, even for smooth and well-behaved functions [START_REF] Lagarias | Convergence properties of the neldermead simplex method in low dimensions[END_REF]. The constraints handling is also an issue as will be discussed in in the following.

DIRECT

The modication of Lipschitzian optimization motivated the DIRECT optimization algorithm development [START_REF] Jones | Lipschitzian optimization without the lipschitz constant[END_REF]. It was created to solve challenging global optimization problems with bound constraints. DIRECT is a sampling algorithm; it requires no knowledge of the objective function derivatives. Alternatively, the algorithm samples points in the search domain and uses the information it has obtained to decide where to search next. The name DIRECT originates from the shortening of the expression "DIviding RECTangles," which outlines the way the algorithm moves towards the optimal solution.

In the rst step, DIRECT transforms the search space to be the unit hypercube.

The objective function is then sampled at the center-point of this cube. The latter is afterwards divided into smaller hyper-rectangles whose center-points are also sampled. DIRECT identies a set of potentially optimal rectangles at each iteration.

All potentially optimal rectangles are further divided into smaller rectangles. The procedure described above is repeated for a predened number of iterations. When there is no limitation for the number of iteration, DIRECT will exhaustively sample the domain, it will cluster sample points around local and global optima [19] [20].

Other algorithms

There exist other well-known algorithms that exploit dierent ideas.

The pattern search (PS) method samples a set of points around the current point, looking for the one where the value of the objective function is lower than the value at the current point. The set of points is called a mesh. The mesh is formed by adding the current point to a scalar multiplied by a set of vectors called a pattern. If the pattern search algorithm nds a point in the mesh that improves the objective function at the current point, the new point becomes the current point at the next iteration of the algorithm [START_REF] Torczon | On the convergence of pattern search algorithms[END_REF].

Additionally to his signicant contributions to the development of quasi-newton algorithms, Michael JD Powell developed many derivative-free algorithms from which we cite COBYLA, BOBYQA, ... [START_REF] Powell | A view of algorithms for optimization without derivatives[END_REF].

Extension to constrained problems

Dierent strategies can treat the constrained problems depending on the algo-rithm; one of the most used is the penalty method. The penalty method replaces a constrained optimization problem by one or a sequence of unconstrained problems whose solutions ideally converge to the solution of the original constrained problem.

To recall, the constrained problem is written as

min p f (p) s.t.g i (p) ≤ 0, i = 1 . . . n c (1.22)
Using the penalty method, it can be transformed

min p f (p) + P (p) (1. 23
)
where P is the penalty function. In the literature, we can nd dierent penalty methods such as the augmented Lagrangian function, barrier function, interior point, etc.

Each of them exploits dierent ideas with more or less complicated implementation [START_REF] Bazaraa | Nonlinear programming: theory and algorithms[END_REF]. The exterior penalty method is commonly used because of its simple implementation. The exterior penalty function is written as

P (p) = M nc i=1 max(0, g i (p)) (1.24)
M is the penalty parameter. The idea is to start with a relatively small value of M .

Subsequently, we will solve a series of unconstrained problems with monotonically increasing values of M chosen so that the solution of each new problem is close to the previous one. This will prevent any notable diculties in nding the minimum from one iteration to the next. Sometimes, a well-chosen value of M can lead to the optimum in one iteration.

Extension to multi-objective problems

As for the derivative-based algorithms, heuristic-based algorithms can be extended to multi-objective problems through weighted-sum or -constraint methods.

Extension to black-box problems

Heuristic-based are well suited for black-box problems since no additional information is needed except the response of the model to a set of variables.

C Population-based algorithms

The previous optimization algorithms maintain a single current solution. Populationbased algorithms maintain multiple solutions. They are usually motivated by biological evolution, such as the evolution of species, behaviour of reies, etc. [24] [25].

One of the most popular is the genetic algorithm (GA), which is inspired from natural selection [START_REF] Goldberg | Genetic algorithms and machine learning[END_REF]. GA starts by creating a random initial population, then creates a sequence of new populations. At each iteration, the algorithm uses the individuals in the current generation to create the ones of the next population. For this goal, the algorithm performs dierent operations, i.e. elitism, selection, crossover and mutation. The algorithm is detailed below. As for the heuristic-based algorithms, the penalty method can be used to handle constraints in population-based algorithms. Nevertheless, some other strategies were developed, such as the one developed by K. DEB [START_REF] Deb | An ecient constraint handling method for genetic algorithms[END_REF].

Extension to multi-objective problems

Population-based algorithms are well suited to multi-objective problems; they are sometimes dedicated to only multi-objective optimization since population-based algorithms dominate multi-objective optimization literature. This is explained by the fact that there are many solutions in each population. These solutions ideally converge the desired Pareto front. NSGA-II is one of the most known optimization algorithms to deal with multi-objective problems [START_REF] Deb | An ecient constraint handling method for genetic algorithms[END_REF].

Extension to black-box problems

No additional modications are needed for the population-based algorithms to handle black-box problems, except that the fact that most black-box problems are costly and special measures should be undertaken to make the optimization tractable in a limited time.

D Exact algorithms

This category of algorithms is less addressed in the literature because of the implementation complexity and the computational cost needed for solving specic problems. However, exact methods are the only ones that can guarantee global optimality of the solution. This feature is highly needed but generally sacriced for the sake of computation time by using other algorithms. Branch and bound (B&B)

is the most known of exact algorithms [START_REF] Lawler | Branch-and-bound methods: A survey[END_REF]. Initially, it was used to solve integerlinear problems. The process of resolution is simple. First of all, the problem is solved as if the variables are continuous (relaxed problem); if the solution variables are integers, then it is the optimal solution. Otherwise, one variable is chosen to split the initial problem into two sub-problems, then solve the corresponding relaxed problems to compute the bounds. The process continues until a satisfactory solution is found.

An extension to optimizations of the general case (nonlinear and continuous) was introduced, it relies on the estimation of bounds using interval arithmetic [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF].

f i , f i of f on V i ; if f ≥ f i then insert (f i , V i ) in L; if f ≥ f i then update f = f i ; eliminate from L all the domains (z, Z) for which z > f ; end end end until f -min (z,Z)∈L z < ε;
B&B is a systematic enumeration of candidate solutions by splitting the search space into smaller spaces. The algorithm explores each smaller space then estimates the upper and lower bounds on each one of them. A space is pruned if it cannot produce a better solution than the best one found so far by the algorithm.

Extension to constrained problems

Constraints can be handled relatively easy by modifying the previous algorithm slightly. The constraints bounds need to be calculated, then the elimination will be based on the objective function bounds and the constraints bounds. The subdomains with at least one of the constraints that have a lower bound higher than 0 are eliminated.

Extension to multi-objective problems

Multi-objective B&B is less addressed in the literature because of the double complexity, rst, for the multi-objective aspect and secondly the cost of the exact search [START_REF] Brisset | Pareto-based branch and bound algorithm for multiobjective optimization of a safety transformer[END_REF].

Extension to black-box problems

Minimal cases are presented in the literature where exact algorithms were able to deal with black-box problems [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF] [31] [START_REF] Fontchastagner | Résolution du problème inverse de conception d'actionneurs électromagnétiques par association de méthodes déterministes d'optimisation globale avec des modèles analytiques et numériques[END_REF]. The estimation of the bounds of the functions is dicult or even impossible. Thus, the usage of exact algorithms is not suitable since these bounds are key-ingredient for their success.

E Approximation-based algorithms

The basic idea is that if the function is noisy and likely expensive to evaluate, then that function can be sampled at some points, and a t of it is created. Then, the optimization is not performed on the original function, but on the cheap to evaluate and smooth t [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF]. These ts are referred to by many names such as approximation models, predictors, surrogate models, meta-models, and response surfaces. These terms will be used synonymously throughout this dissertation. The goal of using a surrogate model is to provide a smooth functional relationship of satisfactory delity to the real function with the added advantage of quick computational speed.

The approximation could be used in combination with a derivative-based algorithm or other algorithms. The details on how to build and to exploit approximations eectively are still keeping black-box optimization a thriving research area.

It is worth noting that the approximation itself is not an optimization algorithm; however, its combination with another or many optimization algorithms aforementioned is an optimization algorithm.

Many approximation methods exist, such as polynomial approximation, radial basis function, neural networks, etc. Nevertheless, the underlying usage procedure is almost the same with minor discrepancies, and the abstract algorithm is presented below.

Because the process of tting the models and locating the new sample points can be viewed as optimization problems themselves, the burden associated with approximation-based algorithms can be signicant. On the other hand, other methods, such as genetic algorithms or derivative-based algorithms, require minimal computational eort in determining where to evaluate the functions next. However, they

Response surface methodology

The most simple surrogate model is a second-order polynomial function; it is used by the response surface methodology to have a local approximation of the real function [START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF]. The surrogate model is written as

f (p) = a + nv i=1 b i p i + nv i=1 nv j=i c ij p i p j (1.25)
where f is the approximation or the predictor of the real function f . The hat notation will be commonly used throughout this dissertation to denote the predictors. a, b i and c ij are the parameters of the predictor. They are identied using a regression strategy based on n evals evaluations of the real function. The regression aims to reduce the error between the real function and its predictor in the evaluated points.

It can be written as an optimization problem.

min a,b i ,c ij n evals k=1 (f (p k ) -(a + nv i=1 b i p k i + nv i=1 nv j=i c ij p k i p k j )) 2 (1.26)
The number of evaluations n evals should be at least (nv+1)(nv+2)

2 to be able to capture the behaviour of the real function. These evaluations should be chosen carefully to prevent numerical problems and also to enable the predictor to represent as best as possible the real function in a local region [START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF].

Once the parameters are found, the predictor f can be used to optimize the real function. Generally, the second-order polynomial cannot capture the behaviour of the real function on the whole domain; thus, multiple response surfaces can be t in dierent locations to capture the maximum of information about the real function.

Other surrogate models enable to t only one response surface on the whole domain and capture the global behaviour such as radial basis function and its generalization Kriging.

Kriging

Kriging is a data interpolation scheme that has its roots in the eld of geostatistics [START_REF] Cressie | The origins of kriging[END_REF]. The technique was adapted for data coming from deterministic computer simulations. This form of data collection and approximation is known as Design and Analysis of Computer Experiments (DACE), Kriging is based on a regression term and a stochastic term [START_REF] Martin | Use of Kriging Models to Approximate Deterministic Computer Models[END_REF]. The regression part is a polynomial function, while the stochastic term aims to compensate for the error due to regression and is constructed based on the location of the sampled points. While the idea behind Kriging is simply put here, the details are left for the next chapter.

Extensions

In recent years, numerous advances have been made in approximation-based algorithms. Constraint handling strategies and extension to multi-objective optimization have been discussed [START_REF] Schonlau | Computer experiments and global optimization[END_REF][38] [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF][40] [START_REF] Martin | Use of Kriging Models to Approximate Deterministic Computer Models[END_REF].

A signicant area where little research has been done to date is parallelization.

El Bechari et al. discussed the limitations of some of the strategies [START_REF] Bechari | Enhanced Meta-Model-Based Optimization Under Constraints Using Parallel Computations[END_REF]. The parallelization strategies represent one of the main contributions of this dissertation and are discussed in the next chapter.

Discussion

The models of electrical machines are usually considered as computationally expensive black-box; the time needed for one simulation can vary from minutes to hours or even days for more complex phenomena. Thus, not all the algorithms are suited for this kind of model. The choice of the algorithm is critical for the success of optimization, mainly in terms of the quality of the solution and the computational eort needed for obtaining a reliable solution. Additionally, these models can be noisy due to the discrete aspect which furthermore limits the choice of the algorithms.

The ideal algorithm is the one that has the least cost leads to the global solution handles the constraints eciently is capable of dealing with numerical noise Nevertheless, such an algorithm does not exist. The No free lunch theorem demonstrates that if an algorithm performs well on a particular class of problems [START_REF] Ho | Simple explanation of the no-free-lunch theorem and its implications[END_REF], then it necessarily pays for that with degraded performance on the set of all remaining problems. Thus an inevitable compromise should be found. We consider that the cost is a higher priority criterion than the other criteria and to choose the algorithms that oer a margin of improvement to the studied problem.

The derivative-based algorithms oer the best performances in term of cost and constraints handling, thanks to the rigorous foundations of the algorithms nevertheless they perform poorly for problems with many local minima due to their local properties, and also they perform inadequately in terms the noise handling since the derivatives are approximated by nite dierence. At the same time, this approximation is highly sensitive to numerical noise.

Exact and population-based algorithms are not adapted to our problem because of the excessive computational cost needed for performing an optimization.

Heuristic-based algorithms are good at handling the numerical noise since no derivatives are needed. However, they perform very bad at constraints handling because the choice of the penalty function highly conditions the results. This is caused when the constraints do not have the same magnitude. If some constraints are dominant, the algorithm will steer towards a solution that satises those constraints at the expense of searching for the minimum. Similarly, the value of the penalty parameter should be xed so that the penalty term's magnitude is not much smaller than the magnitude of the objective function. If an imbalance exists, the objective function's inuence could direct the algorithm to head towards an unbounded minimum, even in the presence of unsatised constraints. In either case, convergence may be exceedingly slow.

As noted before, approximation-based algorithms are very exible since the surrogate model can be combined with one or many of the algorithms mentioned above.

Their development aims to reduce the computational cost needed for optimization.

Dierent ideas were studied in the literature. One particular optimization algorithm is known as Ecient Global Optimization (EGO), developed by Jones et al. [START_REF] Jones | Ecient Global Optimization of Expensive Black-Box Functions[END_REF]. It employs Kriging surrogate models as the approximation method. Other ideas, such as genetic algorithm assisted by surrogate models, are also of interest.

Uncertainty

The optimization problem types described in the previous section implicitly assume that the variables for the given problem are known precisely. For many actual problems, however, the problem variables cannot be known deterministically for various reasons such as measurement error or tolerances on the manufacturing process. In most cases, deterministic optimization without considering the uncertainties of vari- Experimental uncertainty This comes from the variability of experimental measures. The experimental uncertainty is unavoidable and can be noticed by repeating a measurement many times using the same settings for all inputs/variables.

Variable uncertainty This arises from the variability of the input variables of the model. For example, the dimensions of a workpiece in the manufacturing process may not be exactly as designed, which would induce variability in its performance.

This dissertation deals with the last category of uncertainties that enables one to get relationships between the model input variables and their impact on its responses. Thus, the uncertainty on the variables induces uncertainty on the response of the model. This means that the uncertainty propagates through the model to generate uncertainty on the responses. This process is called uncertainty propagation, and it is performed using simulation and numerical tools.

The models are generally deterministic in the sense that model inputs, i.e. variables, are assumed to be known. However, we have seen that this is not always true. Deterministic models are then no longer directly adapted to handle this type of input data. It is necessary to use ecient methods to propagate the uncertainties on the input data to the model's output quantities.

One of the challenges to using models is to capture the eect of uncertainty on the model output in an ecient manner. The problem can be stated as: given the probability density function (PDF) of the random variables (e.g. those of the design variables and parameters), what are the statistics of a system output? The issue is how to propagate the eect of uncertainty through the model. This process is often referred to as uncertainty quantication. In this section, we will start by detailing the uncertainty propagation methods and, afterwards, the integration of these methods in an optimization procedure is presented.

Uncertainty propagation

As stated before, uncertainty on the variables leads to uncertainty on performances.

Thus, the rst stage of any uncertainty propagation procedure is to characterize the input of the model. The input variables are considered as random variables characterized by their PDF. Other statistics such as mean, variance, range, skewness and covariance can also be characterized. The pdf of random variables are not always needed but generally turn out to be useful for the variables characterization.

Once the random input variables are characterized, numerical tools are needed for the propagation of uncertainty, that is to say, compute statistics on the indicators of performance. In the literature, many algorithms enable to propagate uncertainty.

We cite ve groups.

Monte-Carlo simulation and variants [START_REF] Hammersley | Monte carlo methods[END_REF].

Approximations methods [START_REF] Liu | Prise en compte des incertitudes dimensionnelles introduites par les procédés de fabrication dans les modèles numériques de machines électriques[END_REF][46] [START_REF] Nguyen | Uncertainty Quantication Using Sparse Approximation for Models With a High Number of Parameters: Application to a Magnetoelectric Sensor[END_REF].

Perturbation methods (rst and second order) [48][49].

Most probable point (MPP) [START_REF] Du | A most probable point-based method for ecient uncertainty analysis[END_REF].

Worst case method [START_REF] Pugfelder | Worst case optimization: a method to account for uncertainties in the optimization of intensity modulated proton therapy[END_REF].

It is worth noting that each group is adapted to the performance metric needed.

Many performance metrics describing the behaviour of output under uncertainty have been suggested in engineering design. Perhaps the most common performance metric is the expected value of output under probabilistic uncertainties, also known as the mean. Another property of interest is a measure of the spread of the outputs, since if a design has a good expected value of performance but a high spread such that it can also give poor performance it could be undesirable. The standard measure of the spread is the variance. If a system is being designed for reliability, or the output is to be constrained, then the probability of failure will be of interest. A common measure of this sort is that the probability of failure is higher than some threshold value.

A Monte-Carlo simulation

One of the most widely used techniques for uncertainty propagation is the Monte Carlo (MC) sampling technique, which is based on a pseudorandom generator used to approximate the desired distribution. MC evaluates the deterministic model at the distribution realizations and, lastly, aggregates the model's outputs to get insights on the uncertainty of the performance. One could compute the mean, the variance and probability of failure with ease [START_REF] Hammersley | Monte carlo methods[END_REF]. MC sampling is a useful propagation technique because its convergence is independent of the dimensionality of the problem.

While this approach is attractive due to its simplicity, it has numerous deadfalls.

For example, a large number of samples must be taken in order to determine low probability events, e.g. thousands or millions of trial design variable vectors. The computational expense of this approach can quickly become intractable.

Other computational cost reduction techniques exist, such as Importance sampling, Latin Hypercube Sampling (LHS), and Hammersley Sequence Sampling. They enable one to sample the space in a smart way to reduce the number of sampled points. Nevertheless, sampling methods remain costly when considering computationally expensive models, such as a nite element analysis.

B Approximation methods

When the number of variables is low, more ecient methods are available than MC sampling. They primarily involve building a surrogate model as an approximation to the actual output, as presented for deterministic optimization, using relatively few model evaluations. Such a surrogate model can subsequently be sampled and analyzed at a negligible computational cost compared to evaluating the actual output or sometimes computing the performance metrics analytically. Polynomial Chaos expansion is widely used in engineering design [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF] [53] [START_REF] Nguyen | Uncertainty Quantication Using Sparse Approximation for Models With a High Number of Parameters: Application to a Magnetoelectric Sensor[END_REF], among other surrogate models.

C Perturbation methods

These methods have advantages when dealing with relatively small input variability and outputs that do not express high nonlinearity [54][48]. The model must usually be linearized by approximation to a rst-order or second-order Taylor series expansion.

The rst-order expansion of a function f with respect to a random variable X of mean µ X and variance σ 2

X is written. f (X) ≈ f (µ X ) + (X -µ X ) T ∇f (µ x ) (1.27)
The mean value µ f and the variance σ 2 f of the function f are computed as follows.

µ f = E(f (X)) ≈ E(f (µ X ) + (X -µ X ) T ∇f (µ x )) = f (µ x ) + E(X -µ x ) T ∇f (µ x ) µ f = f (µ x ) (1.28) σ 2 f = V(f (X)) ≈ V(f (µ x ) + (X -µ X ) T ∇f (µ X ))) = V(X -µ x ) T (∇f (µ X ))) 2 = V(X) T (∇f (µ X ))) 2 σ 2 f = σ 2 X T (∇f (µ X )) 2 (1.29)
The uncertainty propagation is done easily. The mean value of the function is obtained by evaluating the deterministic function at the variable mean value.

Similarly, the variance is the product of the gradient squared of the function and the variance of the variables.

As noted before, the gradient is not always a piece of accessible information;

thus, this method can be applied only if a gradient of good quality is available. It is also important to note that this formula is based on the linear characteristics of the gradient of f . Therefore, it is a reasonable estimation for the standard deviation of f as long as σ X are small enough.

D Most probable point MPP

The MPP concept is utilized to compute the whole cumulative distribution function (CDF) of model output by evaluating the probability estimates at a serial of limit states [START_REF] Du | A most probable point-based method for ecient uncertainty analysis[END_REF]. The MPP method was initially developed in the eld of reliability analysis. It requires that limit-state (constraints) functions be dened, and the probability of the limit-state functions less than zero can be evaluated approximately.

Generally, the MPP method is used to estimate the probability of an inequality function's failure and not the whole CDF. If a function of a random variable is needed to be less than zero, then the probability of failure is written as

P f = 1 -P (g(X) ≤ 0) = P (g(X) > 0) (1.30)
The rst-order reliability method (FORM) has been widely used in the probability of failure estimation. The method involves Taylor's expansion of the inequality function, i.e. the linearization of the limit-state function, not performed around the mean value of the function, but at the MPP. The selection of an appropriate linearization point is a crucial consideration and leads to an iterative solving procedure, i.e. an optimization algorithm.

E Worst-case method

With the worst-case analysis, all the uncertainties are assumed to lead to the worst possible combinations [START_REF] Pugfelder | Worst case optimization: a method to account for uncertainties in the optimization of intensity modulated proton therapy[END_REF]. Based on this assumption, the worst-case considers the most severe possible outcome that can occur for the uncertainty of the given variable. It is a non-probabilistic approach that is based on the estimation of the worst possible outcome. There is no need to determine the distributions of functions as opposed to the MPP method that identies a probability of failure. In the worst case, the probability of failure is zero.

Optimization under uncertainty

The formulation of an Optimization under Uncertainty (OuU) primarily concerns how to compare the model's behaviour under uncertainty between designs. The formulation denes the underlying landscape that an optimizer attempts to navigate and therefore inuences the design obtained from performing an optimization.

Therefore, the performances are related to two types of variables: rst, the optimization p, which are deterministic and secondly, the random uncertain variables X. Often, the deterministic variables are nothing but the mean value of the random variables. Thus for the sake of presentation simplicity, this assumption, i.e. p = µ X , is holding throughout this dissertation unless stated otherwise.

The original development of methods for taking account of uncertainty in the design process of engineering systems is often credited to Taguchi and his "robust design" approach. This framework initiated an interest in controlling robustness through the design of the system instead of trying to reduce the uncertainties on the inputs directly. Taguchi's method is not necessarily an optimization-based approach; instead, it makes use of the design of experiments. The design space is sampled, and the best design from these samples is selected. Since this initial progress, a mixture of methods for design under uncertainty has been proposed for which robustness is the priority.

The recent growth of interest in uncertainty quantication can perhaps be attached to the importance of rigorously quantifying the inuence of uncertainties for engineering. The eld of OuU grew from the desire to use this rigorous uncertainty quantication within design optimization. In the past few decades, OuU techniques have been applied to various engineering applications, and therefore a variety of approaches for formulating an OuU problem have been proposed.

A Robustness, reliability or worst-case

The appropriate formulation depends on the application. In engineering applications, the constrained quantities can represent an extreme failure of the system, and so the likelihood of constraint satisfaction is a key design driver; this is often known as "reliability-based design optimization" (RBDO). RBDO formulations often make use of the MPP method for uncertainty propagation [55] [START_REF] Liu | Optimization algorithms for structural reliability[END_REF].

min p µ f (p) s.t. P (g i (p, X) > 0) ≤ P f , i = 1 . . . n c (1.31)
In contrast, "robust design optimization" (RDO) formulations focus on the likelihood of obtaining excellent performance [START_REF] Wang | A unied framework for integrated optimization under uncertainty[END_REF]. Therefore RDO mainly deals with the inuence of uncertainties on the performances to be optimized instead of the constrained quantities. Often, "robust design" refers explicitly to simultaneously improving the expected performance represented by µ f (p) and reducing the variability represented by σ f (p) and σ g i (p).

min p [µ f (p), σ f (p), σ g i (p)] s.t.µ g i (p) ≤ 0, i = 1 . . . n c (1.32)
These "robust" and "reliability" paradigms are not mutually exclusive and can be combined.

On the other hand, worst-case optimization may be expressed as a minimax problem. This means the original values of the objective function and the constraint functions are substituted by their worst values, which is the maximum over the uncertainty set.

min p max x∈Ω(X) f (p, x) s.t. max x∈Ω(X) g i (p, x) ≤ 0, i = 1 . . . n c (1.33)

B Nested-loops to sequential

The general formulation of OuU is mathematically a nested two-level structure, consisting of a design optimization loop that repeatedly calls uncertainty quantication analyses in a series of inner loops. This approach is computationally costing and therefore has limited applicability for real engineering problems. In order to surpass the computational diculties, many formulations have been proposed in the literature [START_REF] Deng | Optimisation robuste pour des dispositifs électromagnétiques[END_REF].

The "sequential optimization and reliability assessment" (SORA) method decouples the RBDO process into a sequence of deterministic design optimization followed by a set of reliability evaluation loops [START_REF] Du | Sequential optimization and reliability assessment method for ecient probabilistic design[END_REF]. It uses the reliability information from the previous iteration to shift the deterministic constraints. The process continues until convergence of both sub-problems is reached. SORA has a structure of two iterative loops. However, the two loops are decoupled; hence it is signicantly more cost-eective than the nested two-level process.

In the same essence, a worst-vertex-based WCO was proposed. It aims to observe the values of the bounds of Ω(X) and predicts the worst vertex value [START_REF] Steiner | Managing uncertainties in electromagnetic design problems with robust optimization[END_REF]. The objective is to determine the directions of ascent in which the values of objective function and constraints increase and then evaluate at the bounds of Ω(X). Consequently, the inner maximization problems are replaced by a few evaluations on the uncertainty space's bounds.

Discussions

Well-known methods for uncertainty propagation have been reviewed in this section.

Methods based on the derivatives remain the best compromise between accuracy and eciency, especially when the coecients of variation of the random input variables are not too large, and the model is almost locally linear.

Derivatives can be used for dierent performance metrics, such as the perturbation method, MPP and, in some cases, the worst-case method. Nevertheless, they require models that implement the derivatives of the response quantities.

Approximation methods are preferred if the derivatives are not available, and the model is considered a black-box, as it is often the case of commercial electromagnetic software.

Electrical machines

An electrical machine is a broad term for machines using electromagnetic elds, such as motors, generators, and transformers. Some are electromechanical energy converters: an electric motor transforms electricity into mechanical power while a generator converts mechanical power into electricity. The moving components in a machine can be rotating (rotating machines) or linear (linear machines). The third category of electrical machines is transformers, which, although they do not have any moving components are also energy converters, changing the voltage level of an alternating current.

Electrical machines are based on the electromagnetic principles that generate force or induce voltages and currents. Five principles allow us to understand this link.

1. Current owing in a conductor will produce a magnetic eld around the conductor.

2. Magnetic materials, with high relative permeability, provide a mean to direct and focus magnetic ux.

3. Motor Action: if current ows in a conductor subject to an external magnetic eld, a force will be exerted on the conductor.

4. Generator Action: if a conductor is moved through a magnetic eld, a voltage will be induced in the conductor.

5. Transformer Action: if the magnetic ux passing through a coil changes, there will be a voltage induced in the coil.

The ability to produce a magnetic eld using coils and electric current can be coupled with ferromagnetic materials to direct and concentrate magnetic ux. Combined with transformer, generator and motor actions, these concepts form the foundation of electrical machine design and operation.

A rigorous physical model was given to model the electromagnetic phenomena, i.e. Maxwell's equations. This model is expressed as partial dierential equations (PDE). These PDE are generally dicult to solve analytically due to various reasons such as a complex studied domain, nonlinear material properties and dynamic time behaviour. Thus, numerical methods that discretize the space and time domain are used. The most used one is the nite element method. Such methods are costly in computational time. In this section, we will present the physical model, its discretization and the resolution using the nite element method.

Physical model

Maxwell's equations can describe electromagnetic phenomena. They are separated into two groups.

Firstly, fundamental equations contain only elds.

∇ × E = -∂B ∂t

(1.34) ∇ • B = 0 (1.35)
where E is the electric eld vector (V /m) and B is the magnetic ux density vector (T ).

CHAPTER 1. LITERATURE REVIEW

The second group links the elds to the sources.

∇ × H = ∂D ∂t

+ J + J s (1.36) ∇ • D = ρ (1.37)
where H is the magnetic eld vector (A/m), D is the electric ux density vector

(C/m
2 ), J and J s are the eddy and imposed current densities respectively (A/m 2 ) and ρ is the electric charge density (C/m 3 ).

For electrical machines, the quasi-static approximation is widely used. It states that displacement currents are negligible in comparison with the imposed currents.

Moreover, electric charges are also negligible.

∂D ∂t |J s |, |ρ| 1
Maxwell's equation in the quasi-static assumption are summarized as follows.

∇ × H = J + J s

(1.38) ∇ × E = - ∂B ∂t (1.39)
The system of equations (1.38 -1.39) is under-determined. To solve such a system, more equations should be given. Constitutive Laws link these quantities throughout electromagnetic material behaviour.

We particularly dene two equations

B = µ(H, t, T, ε, . . . )(H + H c ) (1.40) J = σ(t, T, . . . )E (1.41)
where H c is the coercive eld of magnets, µ is the magnetic permeability and σ the electric conductivity of the material. µ and σ depend on electromagnetic elds and other parameters like time t, temperature T , and mechanical strain ε.

In this work, all the parameters other than electromagnetic elds are considered constant. Thus, these laws are written as

B = µ(H)(H + H c ) (1.42) J = σE
∇ × H = J + J s (1.44) ∇ × E = - ∂B ∂t (1.45) B = µ(H)(H + H c ) (1.46) J = σE (1.47) 1.3. ELECTRICAL MACHINES 1.3.

Mathematical model

For the sake of presentation simplicity, we will only detail the two-dimensional magnetostatic case. The electromagnetic elds do not change in time and are invariant in the dimension perpendicular to the two-dimensional domain. For more information about the general formulation, we refer the reader to [60] [61]. The resolution of electromagnetic PDE was and still extensively studied.

Maxwell equation in the magnetostatic case decouple the magnetic and the electric ends and are reduced to

∇ × H = J s (1.48) B = µ(H)(H + H c ) (1.49) (1.50)
The quantities are then reduced to B

= [B 1 , B 2 , 0] , H = [H 1 , H 2 , 0] , H c = [H c1 , H c2 , 0] and J = [0, 0, J s3 ] .
We denote the studied domain as D ⊂ R 2 . Thus the quantities depend on the two coordinates of the two-dimensional domain. On the other hand, boundary conditions should be applied to solve Maxwell equations on the domain D. Furthermore, additional quantities are introduced to tackle issues. These quantities are called the potentials; they enable to reduce the number of unknowns of the original Maxwell equation and to have functionals that are suciently smooth since the elds (magnetic and electric) in the Maxwell equation present discontinuities on the material interfaces.

A Potential formulation

Since B is divergence-free i.e. ∇ • B = 0, we can nd a magnetic vector potential A such that

B = ∇ × A (1.51) Thus, equations (1.48 -1.49) can be combined in one equation ∇ × 1 µ ∇ × A = J s + ∇ × H c (1.52)
The solution A of this equation is not unique. Any quantity A is a solution for a scalar function u

A = A + ∇u (1.53)
On the other hand, we have

B =    ∂A 3 ∂x 2 -∂A 2 ∂x 3 ∂A 1 ∂x 3 -∂A 3 ∂x 1 ∂A 2 ∂x 1 -∂A 1 ∂x 2    (1.54) A is invariant in third dimension, then ∂A 2 ∂x 3 = 0 and ∂A 1 ∂x 3 = 0. Furthermore, by comparison with B = [B 1 , B 2 , 0] , it follows that ∂A 2 ∂x 1 - ∂A 1 ∂x 2 = 0 (1.55)
This equation has an innite number of solutions, the choice of A 1 and A 2 to be equal to zeros will enable reducing the number of unknowns to be solved.

The invariance in the third dimension and the last choice enable to ensure the uniqueness of the solutions. And, the quantities are written as

A =    0 0 A 3 (x 1 , x 2 )    (1.56) B = ∇ × A =    ∂A 3 ∂x 2 -∂A 3 ∂x 1 0    (1.57)
Therefore, only the third component of A (A 3 ) is considered as an unknown. For the sake of presentation homogeneity, we keep the notation A to denote the unknown function.

B Boundary and interface conditions

Let Γ be the boundary of D. Boundary condition model the behaviour of the system at the limit of the studied domains.

Boundary conditions on magnetic quantities Γ B and Γ H are complimentary on

Γ ( Γ B ∪ Γ H = Γ and Γ B ∩ Γ H = ∅) n × H| Γ H = 0 (1.58) n • B| Γ B = 0 (1.59)
n is the outer normal unit normal vector to Γ.

When passing from a domain to another one of dierent material properties, the electromagnetic elds present discontinuities. Thus, interface conditions ensure to have a certain level of continuity through material interfaces [START_REF] Griths | Introduction to Electrodynamics[END_REF].

Between two domains D i and D j , an interface denoted Γ i,j is dened and n an outer normal vector on Γ i,j . Interface conditions are stated as follows.

n × (H i -H j )| Γ ij = 0 (1.60) n • (B i -B j )| Γ ij = 0 (1.61)
With these conditions, the normal component of B and the tangential component of H are continuous.

C Finite element method

As mentioned before, the equation (1.52) is a partial dierential equation (PDE) (Strong formulation). This equation can be solved analytically for simple systems.

For complex ones, weak formulations are used instead [START_REF] Monk | Finite element methods for Maxwell's equations[END_REF], i.e. Finite Element Method.

The magnetostatic problem is reduced to the following equations.

∇ × 1 µ ∇ × A = J s + ∇ × H c in D (1.62) n × A = 0 on Γ B (1.63) 1 µ ∇ × A × n = 0 on Γ H (1.64)
These equations describe the electromagnetic elds and the corresponding boundary conditions.

The problem in 2D is further reduced since

∇ × 1 µ ∇ × A =    0 0 -∇ • 1 µ ∇A 3    and n × A = 0 is veried if A 3 = 0 and lastly 1 µ ∇ × A × n = 1 µ ∇A 3 • n.
Finally, the problem is wholly reduced to 2D, and by summarizing, we have

-∇ • 1 µ ∇A 3 = J s3 + ∇ • M in D (1.65) A 3 = 0 on Γ B (1.66) 1 µ ∇A 3 • n = 0 on Γ H (1.67) with M = -H c2 H c1 .
The interface conditions are implicitly veried in the two-dimensional magnetostatic problems.

Variational formulations

Multiplying the PDE from 1.65 by an arbitrary test function v such that v| Γ B = 0 and integrating over the domain leads to

- D ∇ • 1 µ ∇A 3 v = D J s3 v + D ∇ • M v (1.68)
Applying Green formulas on the l.h.s [START_REF] Girault | Finite element methods for navier-stokes equations: Theory and algorithms[END_REF].

-

D ∇ • 1 µ ∇A 3 v = D 1 µ ∇A 3 • ∇v - Γ 1 µ ∇A 3 • n v (1.69)
since v vanishes on Γ B (Dirichlet boundary) and 1

µ ∇A 3 • n vanishes on Γ H (Neumann boundary), the boundary integral vanishes too.

In the same manner, Green formula is applied to the term with M.

D ∇ • M v = D M • ∇v (1.70)
This leads to the following variational formulation.

a(A, v) = F, v (1.71) with a(A, v) = D 1 µ ∇A • ∇v (1.72) F, v = D J s3 v + D M • ∇v (1.73)
Using the Lax-Milgram theorem, we can guarantee the existence and the uniqueness of the solution of equation 1.71.

In the next section, we will show how we solve this equation using the nite element method.

Numerical model

Thanks to the Galerkin method, we assume having a nite-dimensional space V h , on which a discrete problem is associated.

Find A h ∈ V h such that a(A h , v h ) = F, v h ∀v h ∈ V h (1.74)
Since the space V h is of nite dimension denoted N h , we let α i denote the basis functions, i.e.

V h = span{α i /i = 1, . . . , N h }. Consequently, the solution A h is of the form A h (x) = N h i=1 u i α i (x), x ∈ D (1.

75)

Replacing A h by 1.75 and v by the basis functions α i , the problem becomes

Find u = (u i ) N h i=1 such that N h i=1 a(α i , α j )u i = F, α j j = 1, . . . , N h (1.76)
We can write Find u such that Ku = b

(1.77)

with

K ij = a(α i , α j ) = D 1 µ ∇α i • ∇α j (1.78) b j = F, α j = D J s3 • α j + M • ∇α j (1.79)
Now we face the problem of how to construct the space V h and corresponding basis functions. In the following paragraph, we are discussing this matter.

A Special choice of V h

The FEM is a particular case of the Galerkin method in the construction of V h in three main aspects:

1. A discretization T h of the domain D has to be established. 2. V h is constructed such that for each element T ∈ T h , α i ∈ V h is a polynomial.
3. There is a canonical basis of V h can be easily described and have small support.

The discretization T h of the domain D is called a triangulation or mesh. T h must satisfy the following properties Each T ∈ T h is a closed set with a nonempty, connected interior T .

Each T ∈ T h , its boundary ∂T is Lipschitz continuous.

D = ∪ T ∈T h T .
For each distinct T 1 , T 2 ∈ T h , the intersection is of measure zero. Dierent kinds of triangulation exist. Figure 1.2 shows possible triangulation of a squared domain using small triangles, rectangles could also be used.

The most simple but widely used elements are the linear and quadratic elements.

The domain is meshed, e.g. triangles for 2D domains, on which the functions α i should be linear or quadratic. The construction is as follows.

Let the vertices of the element (called nodes) be enumerated by an index set Id n and let x i , i ∈ Id n denote the corresponding nodes.

We set

α i (x j ) = δ ij , ∀i, j ∈ Id n ,
and with the restriction α i | T rst or second order polynomial, ∀T ∈ T h Such a basis is called the nodal basis since the unknowns are the function values on the nodes. Likewise, we can dene edge basis and facet basis in 2D and 3D.

Constructing V h in this manner enables one to make the matrix K sparse. Consequently, relatively easy to solve the system 1.77.

The stiness matrix K and the load vector b i are assembled by summing over the elements of the mesh.

K ij = D 1 µ ∇α i • ∇α j = T ∈T h T 1 µ ∇α i • ∇α j (1.80) b j = D J s3 α j + M • ∇α j = T ∈T h T J s3 α j + M • ∇α j (1.81)
The appearing integrals over the elements T can be evaluated eciently using the concept of the reference element.

B Reference element

In the 2D case, an element is a triangle. (x, y) = (x 0 , y 0 ) + (x r , y r )

x 1 -x 0 y 1 -y 0 x 2 -x 0 y 2 -y 0 (1.82)
Consequently, the integrals can be easily transformed into the reference element using the mapping. For eciently evaluating these integrals, numerical integration is used (Gauss quadrature).

Let us suppose that (x g , y g ) and w g are the points and weights of a Gauss quadrature of the reference element, then, the basis function and their derivatives are evaluated and written as

α g i = α i (x g , y g ) (1.83) ∇α g i = ∇α i (x g , y g ) (1.84)
Consequently, the integrals in (1.80) and (1.81) are computed as

K ij = T ∈T h g 1 µ ∇α g i J -1 T • ∇α g j J -1 T |J T |w g (1.85) b j = T ∈T h g J s3 α g j + M • ∇α g j J -1 T |J T |w g (1.86)
where J T is the Jacobian of the mapping from the reference element to the element T and |J T | is its determinant.

J T = x 1 -x 0 y 1 -y 0 x 2 -x 0 y 2 -y 0 (1.87)

C System resolution

Once the assembling of the stiness matrix and the load vector based on the contribution of each element of the triangulation is done, the system of equations (1.77) has to be solved. There exist multiple solvers. Direct solvers aim to decompose the matrix K to the product of two or more matrices that are relatively easy to solve;

for example, we can cite Cholesky factorization. However, iterative solvers are still the mostly used type of solvers since they usually lead to a considerable reduction of resolution time and memory.

It is worth noting that the system 1.77 may not be linear in the existence of ferromagnetic material, i.e. K is a function of u. Iterative solvers like Newton-Raphson or xed-point algorithms are then used for treating the nonlinearity of the system.

D Results interpretation

Usually, the solution of the system of equation is not of interest; nevertheless, other quantities computed from that solution are the ones of interest such as the magnetic ux density, the energy, the torque, the electromotive force, . . . These quantities are often used to make "colourful" visualization; these are easy to interpret and give insight about the electromagnetic phenomena in the studied domain.

MagFEM -2D magnetostatic FE code

We have developed a two-dimensional nite element code to investigate the challenges of using such a tool to simulate electrical machines. The developed code is far from being a complete interface for the modelling of electrical machines; however, it will enable us to test the developed approaches on some optimization benchmarks

from the literature [1][2].
The general structure of any nite element code is presented in Figure 1.4. At rst, the geometry of the device should be created, and material properties are assigned to their corresponding parts. Then a triangulation technique is used to subdivide the domains into small elements. Afterwards, the mathematical model is applied to each of these small elements. The latter are then assembled into a global system of equations. This system is solved to compute the quantities of interest.

We chose to explain these steps using a simple analysis of a device taken from FEMM website [START_REF] Meeker | Inductance calculation example[END_REF]. The analysis consist in computing the inductance of the gapped inductor shown in For this model, the ferromagnetic is considered to be linear with a relative permeability of 2500, and we consider that the current is distributed uniformally in the slots.

A Geometry and material properties

The geometry of the electromagnetic device needs to be created, and material properties are assigned to the corresponding domains. Generally, symmetries may exist in the studied device. In our example, we can model only half of the device since the other half is anti-symmetric. Figure 1.6 shows half of the modelled geometry.

Additionally, a boundary domain is dened to impose boundary conditions. Furthermore, we assign a magnetic permeability of 2500 to faces F 2 and F 3 while we assign a current density of 0.8184A/mm 2 to face F 4.

B Mesh

Now, we need to mesh the studied domain. In Figure 1.8, we show two meshes of the studied domain, one with elements that are bigger than the other one. Naturally, the one with smaller elements will lead to more accurate results than the other one.

However, it will be more expensive to evaluate. Once the material properties are dened and the domain is discretized, the assembling of the matrix and the second member is conducted. Afterwards, the equation is solved; since all the materials are linear, the solution consists in solving a linear system and u = K -1 b. We recall that the ux density B is expressed as follows.

B =    ∂A h ∂x 2 -∂A h ∂x 1 0    (1.88)
which can be computed from the solution of the linear system using equation (1.75).

D Post-processing

As noted before, the solution u is only an intermediate result; generally, other quantities are computed based on u are of interest. For example, we have made the plots in Figure 1.9 by computing the magnetic ux density from u. In this example, we are interested in computing the inductance of the magnetic device. To do so, we use the following formula.

W = 1 2 Li 2 (1.89)
where W is magnetic energy stored in the device, L is the inductance and i is the current owing in the coil.

We recall that the current was chosen to be equal to 1A, Thus, to compute the inductance, we need to compute the magnetic energy using the formula.

W = 1 2 D BHdV (1.90)
This quantity is nothing but W = 1 2 u Ku = 1 2 u b, where K and b are resulted from the assembling.

From the reference [START_REF] Meeker | Inductance calculation example[END_REF], a value of W = 0.865 mJ is found. In our analysis, we found W = 0.779 mJ for the coarse mesh and W = 0.868 mJ for the ne mesh.

The ne mesh gives the closest result to the one from the reference. Then, we can compute the inductance by equation (1.89), we nd a value of L = 1.736mH.

As we solve numerically Maxwell equation, the discretization of the domain, i.e. mesh, has an important impact on the solution quality. Ideally, a proper mesh is the one that gives suciently acceptable results without being very ne. Thus, some tuning of the mesh is necessary before conducting a nite element analysis.

Discussion

In this section, we reviewed the literature related to the electromagnetic modelling of electrical machines starting from the physical model and detailing the mathematical one to present the numerical model nally.

Finite element analysis (FEA) enables one to model the electromagnetic elds in the studied domain nely. FEA is now undeniably the most exible modelling methodology that is based on a sound mathematical foundation and able to handle dierent phenomena such as nonlinear material properties, time-varying behaviours and complex geometries. However, all these advantages come with the expense of the increasing computational cost needed to handle such phenomena. Thus, its usage for optimization, i.e. iterative processes, should be made with caution since only a limited number of evaluations of the FEA can be tolerated. Figure 1.10 shows an objective function f (function of the ux density) as a function of a geometry parameter p. An FEA was used to compute the response f by varying the variable p. Note that the response is both noisy and has several discontinuities. Numerically, nding the geometry parameter that minimizes the objective function is a design optimization problem that poses several challenges. These discontinuities are explained by the fact that when varying the geometry variable, the topology of the mesh is also changed.

Chapter Summary

This chapter is considered as the foundation of our research work; at the beginning, we presented the deterministic optimization. Afterwards, the integration of uncertainty related to variables in an optimization loop. Finally, electromagnetic modelling of the electrical machines using the nite element method.

We have highlighted the thesis's main challenges related to the computational cost of the FEA and the numerical error due to re-meshing. These two features limit the choices of well-suited optimization algorithms. Furthermore, we have seen that there are algorithms that are better than others in handling constraints and in treating multi-objective problems.

Consequently, two approaches can be considered, which depend on the capabilities of the model used. A non-intrusive approach where the model is considered as a black-box and surrogate model is t by evaluating the expensive model at few points, then, use the cheap to evaluate surrogate model for optimization and uncertainty quantication. This approach is detailed in Chapter 2. The second approach (detailed in Chapter 3), which is intrusive where the model evaluates not only the outputs but also their derivatives mainly the rst derivative, the derivatives are used either by a derivative-based algorithm or a perturbation method for uncertainty propagation.

The non-intrusive approach using surrogate models is widely used in the context of expensive model optimization. Many researchers had highlighted the challenges of using such a strategy. In this dissertation, we emphasize some of the very known issues and propose new methodologies of using the meta-models and how to speed up the optimization times.

On the other hand, for the intrusive approach, we propose how to compute the derivatives from a nite element model. In this approach, most of the work is done by the model rather than the algorithm. We consider that derivative-based algorithms are cost-ecient, and minor adjustments are needed for improvement, such as a multi-start strategy. The derivatives are eciently computed from the FEA using the adjoint variable method.

1.4. CHAPTER SUMMARY Chapter 2

Non-intrusive approach

The term "Black-box" has come into general use as a description or nomination of an undisclosed process or function about which we do not need to know how it operates -just that it does. A black-box is a device, system or object which can be viewed in only through its inputs and outputs (or transfer characteristics), without any knowledge of its internal workings. Its implementation is "opaque" (black). Through this dissertation, the black-box refers to the FEA that models the electromagnetic phenomena of a device, e.g. electrical machine.

The non-intrusive approach considers that the FEA as a black-box and all the eorts performed for device optimization is on the optimization algorithms. In the rst chapter, we detailed dierent algorithms for solving optimization problems; all these have pros and cons, but some are more well-suited than others for the problematic dealt with in this thesis. Approximation-based algorithms are of interest since they consider the model as a black-box and t a meta-model for the optimization process. This approach decouples the black-box model and the optimization process, which enable to reduce the number of calls to the black-box model (expensive) and renders the optimization tractable in a limited time. The meta-model also oers desired characteristics for the most ecient optimization algorithms, such as smoothness and the possibility to compute the gradients since the meta-model is an analytic function. Thus, derivative-based algorithms will excel in optimizing such models.

In this chapter, all the research has been done on how to exploit meta-models in the context of black-box optimization. A discussion of the limitations of the conventional way of usage is conducted, and we propose strategies to tackle these.

Meta-model design optimization MDO

Meta-modelling has led to new areas of research in simulation-based design optimization. Meta-modelling approaches have advantages over traditional techniques when dealing with the noisy responses and/or high computational cost characteristic of many computer simulations. Meta-models are used in many elds [START_REF] Jin | Comparative Studies of Metamodeling Techniques Under Multiple Modeling Criteria[END_REF], mainly to replace expensive black-box models [67] [68]. In an optimization problem, the objective function and/or constraints are not always cheaply available data. Thus these surrogate models aim to give a model able to approximate the expensive black-box models from a limited number of solutions.

As seen in the rst chapter, there is a multitude of approximation techniques in the literature; This dissertation focuses on one particular technique, Kriging. This decision is not baseless. Other researchers such as Jones et al. [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF] have made the argument that Kriging provides a suitable approximation for computer simulations and has led a thorough comparison between dierent approximation techniques.

Optimization using Kriging meta-models was rst introduced to tackle unconstrained optimization by Jones et al. [START_REF] Jones | Ecient Global Optimization of Expensive Black-Box Functions[END_REF]. This led to the famous EGO for the Ecient Global Optimization algorithm that forms the basis of this chapter. EGO is classied into a general class of optimization algorithms that is referred to as Bayesian analysis algorithms [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]. In this section, we will discuss dierent steps of the algorithm, starting from initial sampling to stopping criteria, with the corresponding literature and our contributions to the subject. We will use the one-dimensional function f (p) = cos(0.6p) + cos(3p) + sin(6p + 1.5)

Meta-model based optimization owchart is presented in

Design of experiment DoE

The choice of Design of Experiment (sampling points) plays a critical role in the accuracy of the meta-model and the subsequent use in prediction. The choice of DoE was addressed in the literature by performing empiric comparison on some test cases [START_REF] Simpson | A concept exploration method for product family design[END_REF]. Here, we detail two families of design of experiments: classical designs and space-lling designs; moreover, we discuss their adaptability to Kriging models.

A Classical DoE

The rst family of DoE consists of designs based on geometric considerations. Fullfactorial designs and central-composite designs belong to this category. They have been initially developed in the framework of linear regression, e.g. quadratic response surface. The idea is to choose the observation points that maximize the quality of statistical inference and minimize the uncertainty of the parameters of the meta-model.

For example, the samples of the full-factorial design are at the boundary of the domain. For the problem dened in Figure 2.2, the samples are p 1 = 0 and p 2 = 4, i.e.; borders of the design space.

In a two-dimensional space, the samples are on the corners of the design space as shown in Figure 2

B Space-lling DoE

A popular alternative to classical DoE is space-lling one. As the name suggests, they aim to spread the sample points evenly on design space. One of the most used designs for the construction of the Kriging response surface is Latin Hypercube Sampling (LHS) [START_REF] Simpson | A concept exploration method for product family design[END_REF]. Figure 2.4 shows a two-dimensional sampling using LHS using four samples. All dimensions are subdivided by four intervals, which produces 4 2 = 16 sub-spaces. The sub-spaces are sampled in a way that each column and each row contains only one point. Each point is randomly sampled in the sub-space. This enables one to have an acceptable spread of the points on the whole design space. The extension of LHS to n-dimensional space LHS is relatively easy.

There exist other sampling techniques such as Halton sequence, Hammersley set and Sobol sequence [START_REF] Santiago | Construction of space-lling designs using wsp algorithm for high dimensional spaces[END_REF]. These methods oer an aordable way of constructing good space-lling DoE, and they are suited for high dimensional problems.

Simpson has performed a comparison of dierent DoE methods on several test cases [START_REF] Simpson | A concept exploration method for product family design[END_REF]. It was found that space-lling designs oer better performances than classical design, and thus, for the rest of this dissertation, we only consider these designs unless stated otherwise.

Fitting Kriging

Kriging has got much interest as a method of curve tting and prediction. To show this, Figure 2.6 exposes the annual evolution of the usage of Kriging in research papers. The growing trend reects its usefulness in various engineering and science elds. Kriging oers some highly desired features such as data interpolation and the ability to compute a prediction error. This section is devoted to the development of the equations found in Kriging.

We will derive the Kriging formulas using a somewhat dierent approach. Readers interested in the standard derivation may consult [START_REF] Schonlau | Computer experiments and global optimization[END_REF] [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Much of the derivation presented below follows the presentation of Jones [START_REF] Jones | A taxonomy of global optimization methods based on response surfaces[END_REF].

A Understanding Kriging

Suppose we have two points p k and p l . Before we have sampled these points using the black-box model f (p), we don't know the function values at these points which can be considered as uncertain. Therefore, let us model this uncertainty by a random variable F (p) of a normal distribution of mean µ(p) and variance σ 2 (p). Intuitively, this means that we are saying that the function f (p) has a typical value of µ(p) and is expected to vary in some range like [µ(p) -3σ(p), µ(p) + 3σ(p)].

Furthermore, the function values f (p k ) and f (p l ) will tend to be close if the distance between the points p k and p l is small. Statistically, this means that F (p k ) and F (p l ) will be highly correlated if the points p k and p l are close to each other. Commonly, the exponential function is used to express the correlation between points, and it is given by

Corr[F (p k ), F (p l )] = exp - nv i=1 θ i |p k i -p l i | α (2.2)
where θ k and α are the correlation parameters. Figure 2.7 illustrates the impact of the correlation parameters on the distance between points.

On the left, α is held x while varying θ. The value of θ i determines how fast the correlation decay as one moves in i-th dimension. Large values serve to model rapid changes even over small distances. The higher θ, the smaller the correlation between two points. It means that the higher θ the faster the variation of the function can be between the two points.

-2 0 2 0 0.5 Going back to the random variable F (p), the uncertainty about the black-box model can be represented by sampling the random variable at N points.

1 (p k i -p l i ) θ = 2 θ = 1 θ = 0.5 (a) α = 2 -2 0 2 0 0.5 1 (p k i -p l i ) α = 1 α = 1.5 α = 2 (b) θ = 2
Y =    F (p 1 ) . . . F p N    (2.3)
The random vectorY has a mean value equal to M (vector of N × 1) and covariance matrix equal to

Cov(Y) = s 2 R (2.4)
where R is a N × N correlation matrix with (i, j) element given by equation (2.2) and s 2 is the variance.

The parameters M, s 2 , θ i and α characterize how we expect the function to vary when moving in dierent coordinate directions. These parameters need to be estimated in order to model the behaviour of the function.

The maximum likelihood estimation (MLE) oers one approach to estimate these parameters. The likelihood function is given by

1 (2πs 2 ) N 2 |R| 1 2 exp -(Y -M) T R -1 (Y -M) 2s 2 .
(2.5)

One should look for the parameters that maximize the likelihood function (LF).

Practically, this means that we look for the predictor that has the same behaviour as the sampled data.

Interestingly, one could easily deduce the values of M and s 2 by setting the derivative of the LF with respect to M and s 2 to zero. The optimal values can then be expressed as a functions of R

s 2 = (Y -M) T R -1 (Y -M) N (2.6) R -1 M = R -1 Y (2.7)
The equation (2.7 is written in that form because, generally, the mean cannot capture all the variations of the functions. Otherwise, s 2 is equal to zero, and consequently, there will be no covariance between the sampled random variables;

therefore, the Kriging model is meaningless. In the literature, the mean is often assumed to be some polynomial function that captures the global trend of the sampled data. By adopting this hypothesis, the mean M is written as

M(p) = f poly (p) β (2.8)
where f poly (p) is a polynomial vector (generally, polynomials of order two or less), and β are the coecients of the polynomial vector. Therefore these coecients β need to be estimated from equation 2.7,

R -1 F poly β = R -1 Y (2.9) β = F poly R -1 F poly -1 F poly R -1 Y (2.10)
where

F poly =    f poly (p 1 ) . . . f poly (p N )    .
(2.11)

Substituting equation (2.6) into equation (2.5), we get the expression of the LF that needs to be maximized.

1 (2πs 2 ) N 2 |R| 1 2 exp -N 2 .
(2.12)

Researchers usually refer to the concentrated log-likelihood function (CLLF), in which the log function is applied on LF, and constant terms are ignored. This function simplies the expression of the LF and changes its range of variation. This will improve upon the performances of the optimization algorithm that solves the problem. The CLLF is written as

-log(s 2 ) - 1 N log(|R|) (2.13)
The CLLF depends only on the correlation matrix R and, hence, on the correlation parameters θ i and α.

At this stage, we have modelled the sampled data by characterizing the parameters of the correlation function, i.e. θ and α. To understand how we can predict at some point p * , we follow the same principle as before. We assume that the parameters of the correlation are known, and we add an articial sample at p * to Y that is the unknown predictor f (p * ). This predictor is found by maximizing the updated LF and is given by

f (p * ) = f poly (p * ) β + r(p * ) R -1 (Y -F poly β) (2.14)
where r(p * ) is the correlation vector of length N between the untried p * and the sampled data points (p 1 , . . . , p N ).

The predictor f is made of two terms; the rst term is a regression model that enables to capture the global trend of the function, and the second term compensates the residuals produced by the regression; it is often referred as a stochastic term.

Therefore, the Kriging predictor interpolates the sampled points.

The mean squared error (MSE) of the prediction can be deduced from the following equation.

M SE(p

* ) = E f (p * ) -f (p * ) 2 (2.15)
The MSE is derived by developing this equation [START_REF] Schonlau | Computer experiments and global optimization[END_REF], and written as

M SE = s 2 1 + u F poly R -1 F poly -1 u -r (p * )R -1 r(p * ) (2.16)
where u(p * ) = F poly R -1 r(p * ) -f poly (p * ). First of all, we will look for the regression coecients of β. As we considered the rst-order polynomial then

f poly (p) = 1 p . (2.18)
Using the exponential correlation function for the stochastic part of the model, the correlation matrix is particularized for this example as:

R i,j = exp -|p i -p j | 2 (2.19)
The resulting correlation matrix is thus (2.20)

R =      1
The correlation matrix is symmetric and reects the basic idea discussed at the beginning of this section; close points have a signicant correlation, while distant points have a lower correlation.

The regression term of the global model can now be estimated using equation Now, new points are predicted using equation (2.14), which is repeated here.

f (p * ) = f poly (p * p * ) β + r(p * ) R -1 (Y -F poly β) (2.21)
where r(p) is a the 4 × 1 correlation vector between an untried point and the sampled data points

r i (p * ) = exp -|p * -p i | 2 (2.22)
Notice that p * needs to be scaled, as were the variables, before using the predictor, the response also needs to be transformed.

The M SE can be computed using (2.16), once the parameters of the Kriging model are identied.

A plot of the resulting Kriging model is shown in Figure 2.8 along with the prediction error σf , the original black-box function, and the four samples. Immediately evident from the gure is the fact that the Kriging model interpolates the sampled points. Furthermore, the prediction error σf (eq. 2.17) is equal to zero at those points and increases when going far from them.

B Parameter estimation

Previously, it was noted that Concentrated Log-Likelihood Function needs to be maximized (or minimize its opposite) to determine the parameters of the correlation function θ and α (2.2). Therefore, an optimization algorithm needs to be chosen for this purpose. The optimization problem is given by We can see clearly that the CLLF monotonically decreases with θ 1 , then the "optimal" solution is the maximum allowed value of θ 1 , i.e. θ 1 . With such a large value, the correlation between points will decrease, and the Kriging model will degenerate, producing the behaviour shown in Figure 2.10. We can distinguish that the Kriging model has a linear part corresponding to the rst-order polynomial regression term and "bell" shapes around the sampled points to interpolate these.

min θ l ,α log(s 2 ) + 1 N log(|R|) s.t. θ l ≤ θ l ≤ θ l 1 ≤ α ≤ 2
Obviously, this model is not well suited and should be avoided when optimizing.

Sasena explained that the points that are evenly sampled in the design space cause this behaviour and proposed to add new samples near the ones that are already sampled [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF] to prevent this erratic behaviour. Furthermore, he proposed a metric to assess if the CLLF is monotonic or not.

An asymptotic study was proposed when θ 1 approaches innity. Thus, the quantities mentioned earlier become:

R =I β = F poly F poly -1 F poly Y s 2 = (Y -F poly β) T (Y -F poly β) N (2.24)
Then, the asymptotic value limCLLF, when θ goes to innity, is equal to

limCLLF = log (Y -F poly β) T (Y -F poly β) N .
( The optimal value of the CLLF found by the optimization algorithm must be strictly less than limCLLF. Otherwise, we conclude that the meta-model is not suitable and special measures should be undertaken. In our example, we found that the optimal value of CLLF is equal to limCLLF = -0.55. This explains the poor tting of the Kriging meta-model.

As proposed in [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF], adding additional points near the already sampled ones could improve the meta-model quality. To verify this hypothesis, we considered three congurations, In cong. 1, we add a new point p = 3.85 near the point p = 3.9, in cong. 2, we add a point p = 2.6 near the point p = 2.65 and in cong.

3, we add a point p = 0.85. CLLF. We can see that adding points near to the sampled ones does not always solve the problem, as stated in [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]. In cong. 2, we can notice that indeed the CLLF is not monotonic because the optimal value of CLLF is less than limCLLF, and we have a minimal value at θ = 139.26. This will have the same behaviour as observed in Figure 2.10; A large value of θ is not very desired since it produces a model with as many local extrema as the number of samples.

Cong

Nevertheless, adding a sample point has solved the problem for the cong. 2, it has produced an acceptable value of θ while verifying the condition that the optimal value of CLLF is less than limCLLF.

The cong. 3 shows that points other than the ones near the already sampled point could prevent the monotonic behaviour of the CLLF. These observations beg the questions: What can one do to prevent the mono-tonic behaviour of the CLLF? Furthermore, where to add the sample points? Two measures should be considered. The rst one is to limit the upper value of the θ to guarantee a minimal correlation between points. The second one is a diagnosis to test the optimal value of CLLF compared to the limCLLF; Points should be added sequentially until the diagnosis is satised.

As we have shown before, adding a single point near sampled one is not sucient to improve upon the short-range variability, and therefore the CLLF may still be monotonic. Hence, additional points should be added until the condition mentioned above is veried. In the next section, we will discuss how to nd those points and add them to the sampled ones.

On the other hand, the maximal value of θ should be imposed, but how?

Considering the correlation function, we want at least two sampled points to be correlated with each other, and their correlation has to be higher than some threshold as shown in the equation below

exp -θ|p k -p l | ≥ c (2.26)
In most cases, the samples of the initial DoE are uniformly distributed on the design space, after scaling, all the variables will be bounded between -√ 3 and √ 3. Furthermore, they are almost equidistant in each dimension; thus, the distance between two samples is around

d = 2 √ 3 N + 1 (2.27)
where N is the number of samples. Thus, by xing the strength of the correlation c we can nd an upper value of θ. θ ≤ -log(c) d

(2.28)

A reasonable choice of c = 0.003 than will lead to an upper bound of θ which is equal to

1.68(N + 1).

It is worth noting that for n v -dimensional problem, the upper bound of θ in each dimension should be

θ = 1.68 N + 1 n v (2.29)
The upper bound depends only on the size N of the initial design and the dimensionality of the problem. Imposing an upper bound on θ will guarantee that, in the worst case, at least some small correlations between the samples are imposed.

We can follow the same reasoning for the estimation of the lower bound on θ. Imposing a maximal possible correlation between two samples to c = 0.997 will lead to the lower bound of θ in each dimension: since the CLLF can be multi-modal (many local minima) [START_REF] Mardia | On multimodality of the likelihood in the spatial linear model[END_REF].

θ = 10 -3 N + 1 n v (2.
The main idea for nding an adequate initial point is similar to the one used for identifying the upper and lower bounds of θ. We dene the initial point, such as the most distant samples are correlated. This will lead to an overestimation of the correlation strength with a small value of θ, where the CLLF exhibits large gradients.

The value of α = 2 is best suited for smooth functions but less numerically robust [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF]; that is why we choose the initial value of α = 1. Then, we assume all the samples are correlated by at least a small factor c. This small correlation corresponds to the most distant samples. As the variables are scaled, this distance between the farthest points is at most equal to 2 √ 3. Then, the initial point of θ 0

θ 0 = - log(c) 2 √ 3 (2.31)
A correlation factor c = 0.25 seems a good value since our experience have shown that underestimated values of θ 0 usually gives better results than overestimated ones.

θ 0 = 0.4 (2.32)
Finally the optimization problem that maximizes the CLLH to nd the correlation parameters is written as

min θ l ,α log(s 2 ) + 1 N log(|R|) s.t. 10 -3 N + 1 n v ≤ θ l ≤ 1.68 N + 1 n v , l = 1, . . . , n v 1 ≤ α ≤ 2 (2.33)
where N is the size of the initial design and n v is the number of variables.

After tting the Kriging model, additional sample points are added to improve its quality either by looking for the optimum of the Black-box model or globally improving its accuracy. Next, we discuss how we determine the sample points for both approaches.

Finding new samples

Optimization algorithms using the Kriging model do not follow a search path as conventional algorithms; it selects points from anywhere in the design space depending on the objective; it is referred to as inll sampling criterion (ISC). To illustrate this, we will start with the famous algorithm EGO developed by Jones et al. [START_REF] Jones | Ecient Global Optimization of Expensive Black-Box Functions[END_REF].

EGO exploits the Expected Improvement (EI) inll sampling criterion. Thus, EGO proposes to maximize the expected improvement function, shown in equation ( 2.34), rather than minimizing the meta-model of the objective function.

EI(p) = f min -f (p) Φ f min -f (p) σf (p) + σf (p)φ f min -f (p) σf (p) (2.34)
where f and σf are the Kriging meta-model and the prediction error of the objective function respectively. f min is the smallest sampled value of f . Maximizing EI(p) leads to the point x * with the most signicant expected improvement under the uncertainty of the Kriging model, either by sampling toward the optimum or improving the approximation of the meta-model. These characteristics can be justied by the fact that the derivative of EI with respect to f is negative, meaning that the smaller f , the higher EI (exploitation) and the derivative EI with respect to σf is positive, meaning that the bigger σf , the higher EI (exploration).

To demonstrate EGO's search strategy for unconstrained optimization problems, From this example, one can see that optimization using the Kriging model does not follow a search path. It chooses points from dierent regions in the design space, depending on where the sampling criterion is the highest.

In general, ISCs have strong inuence on how eciently and accurately the algorithm locates the optimum. They can be tuned to focus either on exploitation and exploration, or a compromise between the two. An exhaustive set of inll criteria was presented in [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF].

A Constraints handling

Constraints are omnipresent in most optimization problems, having ways to deal with these is mandatory. EGO was developed for the optimization of simply bounded problems, that is, with constraints only on the ranges of the design variable values to be considered. Some researchers extended the EGO algorithm to constrained problems. The probability of feasibility (P F ) criterion is one of the most signicant contribution. It quanties the probability that a constraint is satised (g(x) ≤ 0).

P F (x) = Φ -ĝ(x) ŝg (x) (2.35)
where ĝ(x) and ŝg (x) are the Kriging meta-model and the prediction error of the constraint g(x) respectively.

The ISC is then adapted to consider both objectives; to sample point that improves the actual solution and verify the constraints. Thus, both EI and P F should be considered. (2.36) Thus, another formulation was proposed as in equation (2.37) that considers the inll criterion as bi-objective to reduce the modality [START_REF] Durantin | Analysis of multi-objective Kriging-based methods for constrained global optimization[END_REF]. From the Pareto front solution, the point chosen is the one that maximizes the product of EI and P F . (2.37)

In the case of many constraints, both formulations consider the global P F as the product of the probability of feasibility of each constraint. However, this prevents the search close to the boundary of the constraint. So, if the optimum lies on the constraint boundaries, these inll criteria may fail to nd it. The maximization of P F leads to points inside the feasible region and less likely on the constraint boundaries.

Therefore, a new formulation was proposed in equation (2.38). It considers the ISC problem as a constrained one to reduce the modality of the inll criterion and to gain in precision of the solution in case of the optimum lies on the boundary of some constraints [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF].

max p EI(p) s.t. P F i (p) ≥ P tol , i = 1 . . . n c (2.38)
This formulation considers each constraint independently and calculates their respective probabilities of feasibility, ending up with the same number of constraints as the original problem.

A value P tol = 0.95 was recommended but it has an eect on the precision. It fails to locate the points on the constraint boundaries. In our opinion , a value P tol = 0.5 seems more reasonable because P F = 0.5 when ĝ(x) = 0.

To illustrate the sampling strategy for constrained problems, we will add a constraint to our test function. The optimization problem is now written as follows:

min p f (p) = cos(0.6p) + cos(3p) + sin(6p + 1.5) s.t. g(p) = cos(2p) -0.4 ≤ 0 0 ≤ p ≤ 4 (2.39)
In Figure 2.14, the problem is shown graphically. The objective function is the blue dashed line, the infeasible region, where g(p) > 0, is highlighted by the gray area, and an asterisk shows the global minimum at p = 3.721. On this gure, we can see that there are two separate feasible zones, one at the centre and the other at the right where the optimum lies. This constrained test problem is used to show the iterations of the constrained ISC presented in 2.36. In Figure 2.15, The blue dashed line is the true objective function we wish to minimize, while the solid red line is the Kriging approximation conditional to the sample points shown as green dots. The black plot at the bottom is the sampling criterion 2.36; its y-axis is on the right. The Kriging model of the constraints was not put in the gures for better legibility. Visually, we can see that the Kriging was able to approximate quite accurately the constraint after three iterations.

The ISC was able to locate the global solution quite accurately as for the unconstrained case. The new samples are found by maximizing the ISC. In the beginning, the algorithm explores the right region where the optimum lies. In contrast, at the last iterations, it explores the left region where there is a high probability that a better point can be found.

B Extensions to multi-objective optimization

Treating multi-objective problems with approximation methods is still an active research area, as noted by Rojas et al. [START_REF] Rojas-Gonzalez | A survey on kriging-based inll algorithms for multiobjective simulation optimization[END_REF]. There exist dierent ways to treat multi-objective optimization using an approximation method such as Kriging. We will distinguish here two signicant categories:

1. Inll criteria-based

Performance indicator-based

The inll criteria-based category aims to extend the inll sampling criteria from single-objective problems to multi-objective ones. The expected hypervolume improvement is commonly used; it extends the concept of the expected improvement to multi-objective optimization [START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by gaussian random eld metamodels[END_REF].

On the other hand, the performance indicator-based category is used as an assistant to some multi-objective optimization algorithm. For example, K-MOGA uses the prediction error to assist a genetic algorithm [START_REF] Li | A kriging metamodel assisted multi-objective genetic algorithm for design optimization[END_REF]. At each generation, a Kriging model is tted to each objective and used to evaluate each point in the population.

If the prediction error is higher than some threshold for any point in this population, the black-box model is used on that point to return the true response values.

It is worth recalling what we presented in Chapter 1, namely that a multiobjective problem could be transformed into a single-objective one (using weightedsum or ε-constraint) and treat it as a single-objective problem.

C Extensions to optimization under uncertainty

Optimization under uncertainty using approximation methods is scarcely treated in the literature because of the double complexity of optimization and uncertainty propagation, Deng has presented a thorough state-of-the-art on the subject [START_REF] Deng | Optimisation robuste pour des dispositifs électromagnétiques[END_REF].

Most of the formulations that deal with uncertainty are extended from the deterministic case. In fact, new ISCs are developed to explore only promising regions of the design space. This research topic has been gaining importance during the last two decades. Many researchers proposed their methods to reduce the computational burden and improve the design while taking into account the uncertainties [START_REF] Deng | Iterative Kriging-Based Methods for Expensive Black-Box Models[END_REF].

D Solving the inll criteria problem

It was noted that new samples are found by solving the ISC optimization problem.

Although dierent, ISCs often reveal striking similarities. One of the most common is multi-modal behaviour. Indeed the ISC has many local optima. Finding a local minimum is not necessarily a critical matter, but nding the global one is very desired. For example, in Figure 2.13, the expected improvement function has many optima and large regions where it is equal to zero. Optimization algorithms can nd the optimum in the rst iterations. However, in the last iterations, the non-zero region is tiny, which renders nding the global optimum a very delicate task.

Thus, the choice of the optimization algorithm is a crucial part of the success of the whole process. As noted, the ISC can have multiple local optima and can have zero-gradient on large parts of the design space. Thus, derivative-based algorithms are not well-suited. Nevertheless, the meta-model has an analytic formula and gradients can be computed with little overhead. We recall that derivative-based algorithms are the most ecient when the function and its derivatives are smooth enough, as detailed in Chapter 1. Consequently, we choose to bet on these algorithms for two main reasons:

1. The local search drawback, can be dealt with by adopting a multi-start strategy. Indeed, running multiple optimization with uniformly distributed initial points on the design space will necessary increase the chances of nding the global minimum.

Derivative-based algorithms deal eciently with constraints to enable us to treat constrained ISC problems.

We choose to use the SQP implementation of MATLAB for this purpose. The initial points are generated using an LHS to guarantee a uniform distribution. Additional measures were taken to avoid some specic problems, such as scaling the variables and the responses.

Stopping criteria

Proving the convergence of the approximation-based algorithm is very dicult. For most algorithms, the search is terminated by a limit on the number of black-box evaluations. Other stopping criteria have been investigated that involve some statistical property. For example, the expected improvement estimates how the sampled point will improve the actual solution. Thus, a threshold can be dened to stop the algorithm if the maximal value of the expected improvement is less than that threshold.

Nevertheless, these criteria might stop the algorithm prematurely since the optimization algorithm can sometimes fail to nd the global optimum of the ISC.

In the remaining of this chapter, the termination criterion is specied as the maximum allowable number of black-box model evaluations unless we state dierent conditions.

Analytical examples

The performances of the algorithms presented above are tested on two analytical examples.

A Example 1: Dimensionality

The rst example is the Rosenbrock function; H. Rosenbrock introduced it in 1960, which is used as a performance test problem for optimization algorithms. As shown in (2.41), the problem can be extended to many variables; thus, the objective, of using this function, is to study the impact of dimensionality on EGO.

min p f (p) = nv-1 i=1 (1 -p i ) 2 + 100 p i+1 -p 2 i 2 s.t. -2.4 ≤ p i ≤ 2.4, i = 1 . . . n v (2.41)
Rosenbrock problem has one global minimim at p o = (1, . . . , 1) and f (p o ) = 0. The second example aims to study constraint handling; therefore, we consider the optimization problem shown in (2.42). Sasena considered this problem as one of his test problems [START_REF] Sasena | Flexibility and Eciency Enhancements for Constrained Global Design Optimization with Kriging Approximations[END_REF].

min p f (p) = (p 1 -10) 2 + 100 (p 2 -15) 2 s.t. g(p) = p 2 - 5.1 4 * π 2 p 2 1 + 5 π p 1 -6 + 10 1 - 1 8π cos(p 1 ) -5 ≤ 0 -5 ≤ p 1 ≤ 10 0 ≤ p 2 ≤ 15 (2.42)
The problem has a global optimum at p o = (9.7775, 4.5238) and the constraint is active g(p o ) = 0.

Results assessment

We choose to consider these examples in order to test EGO algorithms and assess the quality of the solution and the computational burden expected when treating similar problems. Thus, some specic metrics are proposed to quantify the performances.

1. cost : The number of black-box function evaluations.

2. cost 1% : The number of black-box function evaluations before a point is sampled within a box with size 1% of the design space range around the true solution.

3. error : the normalized Euclidean distance from the best sampled point to the global solution. This quantity is a scaled to design size.

4.

êrror : The normalized Euclidean distance from the solution found by the Kriging model and the global solution.

CHAPTER 2. NON-INTRUSIVE APPROACH

To reduce the impact of the initial sampling on the result, we choose to use the Hammersley sequence rather than LHS as an initial DoE; it generates well-spread points. The size of the initial sample is chosen to be ve times the number of variables.

The Kriging meta-model is t at each iteration; the parameters of the correlation function are found using the Nelder-Mead algorithm. The upper and the lower bounds and the initial point of the correlation function are determined using the procedure described in the previous section mainly equation (2.33).

We consider the expected improvement criterion as an ISC for the rst example and for the second example, the expected improvement is subjected to the metamodel of the constraint. We use Matlab implementation of the SQP algorithm to solve the ISC problem. One thousand points were chosen using an LHS to run SQP using dierent initial points . This signicant number of initial points will enable us to determine a convergence rate. The latter is determined by computing the number of optimization that converged to the best solution (within a radius of 10 -8 ) divided by the total number of optimizations.

As for the stopping criterion, the iterations are kept going until the maximum value of the expected improvement is less than a threshold for two successive iterations. Imposing two iterations will prevent the algorithm from stopping prematurely in case the optimization algorithm fails to solve the ISC problem once. The threshold is taken to be equal to 10 -6 .

A Example 1

Four optimizations are considered with dierent number of variables n v . The results are summarized in Table 2 The total number of evaluations increases as the number of the variable does.

Interestingly, for this problem, the number of evaluations is almost proportional to the power of two (cost ≈ 6 2 nv ). As expected, this rate of increase is exponential.

The cost can be estimated for problems of higher dimensions.

The empty third column of the table shows that EGO could not sample points near the global solution because the algorithm expected no improvement in doing so. The values of error conrm this fact, the distance from the best-sampled point ( the sample having the minimum value of f ) to the global solution is higher than 1.5% for all the cases.

In the two-dimensional case, the algorithms were pretty accurate in the neighbourhood of the global optimum; this is justied by the low value of êrror. The minimum found by the Kriging model is very close to the global one.

For the three-dimensional case, the algorithm was misled since the best-sampled point was not the closest to the global solution. Consequently, the solution found by the meta-model is far from the global one -êrror = 11.5%.

For other cases, the solution found using the Kriging meta-model was not of good quality, as it was far from it by at least 3% of the design space range.

This example shows two drawbacks of EGO; the global solution is not determined reliably neither by :

1. sampling near enough to the global solution or,

having a good approximation near the global optimum

To investigate the causes of these drawbacks, we investigated the evolution of the algorithm iteration-by-iteration. In The decreasing trend of the convergence rate highlights the increased complexity of the ISC problem as the iterations progress. Indeed, we have seen for the onedimensional problem that the expected improvement criterion has zeros-gradient in the majority of the design space, and this explains why it is becoming cumbersome to solve the ISC as iterations continue.

Putting the chosen optimization algorithm to blame is legitimate, but we have considered one of the most ecient algorithms for local search and enhanced it with the multi-start strategy. In our opinion, no other algorithm could outperform the one we propose for solving these problems except an exact algorithm (which, by the way, explains why a branch and bound algorithm was used in the original paper for the optimization of the ISC [START_REF] Jones | Ecient Global Optimization of Expensive Black-Box Functions[END_REF]). Nevertheless, the estimation of the bounds to be used by an exact algorithm becomes computationally expensive as the number of samples increases; thus, the process becomes intractable and impractical.

All that being noted, we suggest that the EGO is not able to sample points near the global optimum because the algorithm was not able to nd the global solution of the ISC.

B Example 2

In this example, two Kriging models are t one for the objective function and the other for the constraint function. The result of the optimization is presented in Table 2 Admittedly, the constraint function is complicated; thus, the meta-model was not capable of approximating it accurately.

Some research has been done on how to locate all the feasible regions before starting optimization. We proposed to use the probability of feasibility of P F criterion [START_REF] Bechari | Enhanced Meta-Model-Based Optimization Under Constraints Using Parallel Computations[END_REF]. Maximizing P F will lead to solutions where the constraints are more likely to be satised.

We also proposed another criterion to look for multiple feasible solutions. This ISC allows sampling points that have the highest probability of feasibility while being distant from each other. This enables identifying multiple feasible zones (if any), and the optimization using another ISC for the optimum search can be conducted afterwards.

However, our experiences have shown that the optimization algorithm, maximizing P F , has the same problems as observed for the expected improvement ISC. Indeed, the P F function has a zero-gradient in large part of the design space; this renders the optimization problem harder to solve.

Discussion

We have noted some of the diculties that need to be dealt with when optimizing using Kriging models. Here we will try to summarize these and cite some others that were highlighted in the literature :

1. The ISC optimization problems are highly multi-modal, and their modality is highly correlated to the size of the sample of the black-box model; when the number of samples rises, the number of local optima increases too. Furthermore, some ISCs can express zero-gradient in large part of the design space.

Therefore, it becomes highly dicult to solve these optimization problems [START_REF] Picheny | A benchmark of kriging-based inll criteria for noisy optimization[END_REF].

2. The inll criteria that focus more on the local search tend to sample points close to each other, which decay the conditioning of the correlation matrix and, in consequence, the time needed for the construction of the meta-model [START_REF] Mohammadi | An analytic comparison of regularization methods for gaussian processes[END_REF][80].

3. The number of samples exponentially increases as the dimensionality of the optimization problem does (the curse of dimension). Thus, the size of the correlation matrix increases and consequently the time needed to t the metamodel. This time can even exceed the evaluation time of the black-box model.

4. When dealing with constrained optimization, the inll criteria tend to sample inside the feasible region but not on the constraint boundary which aects the solution accuracy [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF].

A Solving the ISC

We have seen in Figure 2.18 that solving the ISC becomes challenging as the iterations progress. Even at the beginning of the process, the convergence rate did not exceed 40% and decreased as the algorithm iterates. Solving the ISC is vital for the success of the process of nding a global solution. Thus, the choice of the ISC is crucial.

B Correlation matrix conditioning

As the correlation matrix is inverted many times in the process of tting the Kriging model, its conditioning is an indicator of how dicult model tting will be. The condition number is the quantity that measures the conditioning. A matrix with a low condition number is said to be well-conditioned, while a matrix with a high condition number is said to be ill-conditioned. In For both cases, the condition number increases, which means the the construction of the meta-model becomes since the correlation matrix is not guaranteed to be positive-denite due to rounding errors.

C Fitting cost

Aside from the conditioning, the number of samples also impact the time needed for tting the Kriging model. To illustrate this, we took the problem of Example 1 and tested it for dierent sample sizes and dierent number of variables (n v ). To eliminate the impact of the conditioning, the samples are drawn using the Hammersley sequence to get a space-lling design. We can see that the time needed for model tting increases as the sample size grows. The dimensionality seems not to have a signicant impact on tting time.

This can be explained by the fact that the inversion correlation matrix is the most Nevertheless, a problem of higher dimension needs more samples to approximate suciently the black-box model (the curse of dimension). Therefore, the increased cost related to dimensionality.

n v = 2 n v = 3 n v = 4 n v = 5

D Synthesis

The conventional way of using the Kriging model seems to be not reliable since the meta-model tting is a multifaceted process. Indeed, the conditioning of the correlation matrix and size of the samples profoundly impacts the time needed for model tting and the quality of the tted model. Additionally, traditional ISC problems increase in complexity as the sample size does.

To solve the problems mentioned above, "simple" meta-models need to be tted.

"Simple" stands for meta-models that are easier to t and lead to ISC problems that are tractable using classical optimization algorithms. Many optimization algorithms discussed in Chapter 1 solve optimization problems using many relaxed forms of the initial problem. This approach appears promising in the context of approximationbased algorithms and still an unexplored research path.

A divide and conquer strategy appears to be a good starting point; the initial design space is divided into sub-spaces with models easier to t and problems faster to solve. Then the sub-spaces are compared to further exploration for the most promising ones. The idea is based on the branch and bound (B&B) algorithm. This fundamental property of the Kriging model will enable estimating the bounds to be used by the B&B algorithm.

Proposed algorithm

We propose an algorithm that builds many meta-models each on a specied region of the design space; each one of them is relatively easy to t and easy to use to estimate the upper and lower bounds. Then, iteratively, prune the regions that are not promising. By adopting this algorithm, the drawbacks mentioned earlier are tackled, namely the numerical ones related to the correlation matrix. The owchart of the algorithm is shown in Figure 2.21.

Initialization

In the beginning, the initial space S is dened, and an empty list L is set, L will store the potentially optimal regions. Additionally, a maximum depth level l is dened to limit the exploration of the algorithm. The depth level can be seen as a metric of how many times the initial design space S can be subdivided.

Branching

Branching stands for subdividing the initial design space S. The most intuitive rule is splitting each dimension of S by two. This will lead to K = 2 nv sub-spaces s j .

Unfortunately, the number of sub-spaces increases exponentially as the number of variables grows. This is not desired, and a more practical rule is proposed.

For meta-models, the purpose behind branching is to improve upon it accuracy by reducing the size of the design space. Thus, the space S is split, such as the produced sub-spaces are of better accuracy. Therefore, the meta-model on S is used to identify the dimensions on which high inaccuracy is observed.

We recall that the stochastic term in the Kriging predictor enables us to model the error not captured by the regression term. This term is characterized by the correlation function and thus its parameters α and θ 1 , . . . , θ nv . The value of α characterizes the smoothness of the function while values of θ enable to model the changes over distances. Large values serve to model rapid changes even over small distances. Hence, the value of θ is an indicator of how the function varies.

If the value of θ i , corresponding to a variable p i , is the largest. Then, the change in the i-dimension is more important than other dimensions. Consequently, splitting the space S by dividing its dimension i will lead to meta-models on the sup-spaces that are better than the one on S.

A simple example is used to test this property of the Kriging model, the studied function is a three-dimensional function but that depends on only two of them, and one of the variables present more variation than the other 2.43. From this example, we expect that the value of θ 2 to be higher than θ 1 and θ 3 . So, a Kriging model is t to a DoE of 15 samples and using a linear regression term.

f (p 1 , p 2 , p 3 ) = p 1 p 2 2 , with 0 ≤ p i ≤ 10, i = 1, 2, 3 (2.43)
Indeed, the meta-model tting led to θ = (0.072, 0.171, 0.004) which conrms the idea we proposed. Ordering the value of θ will lead to the variables that have the highest variation. The function does not depend on p 3 then θ 3 = 0.004 which is the smallest allowed value for θ.

More than one variable can be chosen for the branching; these variables are picked in the decreasing order of their corresponding θ values. Then each variable is subdivided in two part which will lead to 2 d sub-space where d is the number of variables picked.

Meta-modelling

Once the sub-spaces are determined using the branching rule, Kriging models are t on each one of them. The initial DoE, such as Hammersley, is determined, and the black-box-model is evaluated.

Additional samples are added using inll criteria to improve the meta-model approximations and to deal with the monotonic CLLHF issue. The choice here is crucial for the estimation of the bounds; we rely on the prediction error σ, new samples are added where the prediction error is high; consequently, we improve the meta-model globally on the sub-space. Additionally, a limited number of points are sampled using expected improvement (EI) criterion to enhance local exploitation, in case the CLLF is monotonic. This will also improve the global performances of the algorithm.

The initial DoE size and number of the points added using σ are specied by the designer. In contrast, the number of points added using EI is conditioned by the monotonicity of CLLF. The algorithm slops if the CLLH function is not monotonic.

Bounding

Once the meta-model is t on each of the sub-spaces, we look for the lower and upper bounds on each of them. As Kriging is capable of giving the predictor and its error estimation, we exploit this information for the computation of bounds using the following formulation (2.44) and (2.45):

f = min p f (p) -kσ f (p) s.t. ĝi (p) -kσ g i (p) ≤ 0, i = 1, . . . , n c . (2.44) f = min p f (p) + kσ f (p) s.t. ĝi (p) + kσ g i (p) ≤ 0, i = 1, . . . , n c . (2.45)
where f and ĝi are the predictors of the objective function f and the constraints g i respectively and σf , σg i are their error estimations given by the Kriging meta- model. The parameter k can be dened as the condence level required, e.g k = 3, it corresponds to a condence level assuming that the prediction distribution is normal.

The optimization problems (2.44) and (2.45) has to be solved using appropriate optimization algorithms due to the multi-modal behaviour introduced by the kσ terms. We used sequential quadratic programming from Matlab for the task and improved by a multi-start strategy with points uniformly spread on the design space, the number of starting points is much higher than the number of samples of the black-box model. Additionally, the gradient information from the meta-model is exploited to increase the eciency of the algorithm.

The optimization problem (2.45) can present diculties in solving it since the objective function and constraints are mostly concave. Furthermore, our experience has shown that f is, most of the time, equal to the best-sampled point. Thus, this optimization can be discarded, and f is set to be equal to the sample having the smallest value of the objective function and satises the constraints on the sub-space.

For highly constrained problems, when no feasible sample is available, f is set to be equal a large value.

In most constrained optimization problems, the problem (2.44) may not have a solution on some sub-spaces. It means that these sub-spaces have no potentially feasible solution, and therefore the values of the bounds are set to some high value.

Elimination

The sub-spaces are added to L with their lower f and upper f bounds. A comparison between the newly added sub-spaces to L and the existing ones is conducted to eliminate the ones that cannot produce a better solution.

First, we set b as the smallest upper bound of all the spaces in L, and then, we prune the spaces from L that they have a lower bound that is bigger than b.

Selection

The remaining sub-spaces in L are all promising ones; thus, one of them has to be chosen for the next iteration to be subdivided. Many strategies could be considered; by preferring depth-rst for prioritizing local search or breadth-rst for prioritizing exploration of the search space. In our work, we use the best-rst; it is based on the lower bound of the objective function on the sub-spaces. Indeed, the sub-space chosen from L has the smallest lower bound, and its depth level is smaller than the maximum specied depth level. This strategy will prioritize sub-spaces that have the maximum potential for providing the best solution without any criteria on the size of the sub-space.

The sub-spaces that have attained the maximum allowed depth level are stored in a list L s .

Termination

The algorithm will naturally iterate until the list L is empty. It is challenging to come up with a robust termination criterion, but a heuristic can be used to stop the exploration of a sub-space if the dierence between its lower bound an upper bound is less than a small threshold specied by the designer. It remains possible to stop the algorithm by setting a maximum number of evaluations of the black-box model. The list L s contains the remaining sub-spaces that are at the maximum specied depth level, or their bounds are tight. One should investigate them, and search for the optimum in one them without considering the prediction error. Otherwise, if no solution is apparent, the algorithm could be restarted with the list L s and a higher maximum depth level.

Analytic Examples

The problems (2.41, 2.42) treated using the conventional approach are solved using the new algorithm. The setup of the algorithm parameters is as follows

The maximum depth level is two times the number of variables.

The parameter k in (2.44) and (2.45) is taken to be equal to three On each sub-space, the initial sample is chosen using Hammersley sequence with two times the number of variables, the same number of samples is added using the prediction error as an ISC and two points are added using EI.

Figure 2.22 shows the evolution of the sub-space division applied to the problem treated in Example 1 (n v = 2). As shown, on each sub-space, we compute the lower and upper bounds, i.e f and f , respectively, the new sub-spaces are compared with the ones in the list L to discard unpromising ones. The shaded area shows the sub-spaces discarded from the list L.

After only four iterations, the algorithm located the sub-space where the optimum lies. The green region is the only space that was left in the list L. When searching in this space without considering the prediction error, the solution found is close to the the global minimum of the problem.

The results for the remaining problems are summarized in table 2.5. Two of the metrics proposed before are not relevant here (cost 1% and error) since the algorithm itself does not aim to sample at the optimum but to have a meta-model of good quality around the optimum. The algorithm performed very well in locating the optimums for dierent problems; the relative error between the global solution and the solution found using the meta-model is less than 1% for all the cases.

n v cost
For Example 2, where the problem is constrained, The algorithm performed only one iteration because the dierence between the lower and upper bounds is

-2 -1 0 1 2 -2 -1 0 1 2 f = 1.003 f = 0.276 f = 0.107 f = -3.837 (a) Iteration 1 -2 -1 0 1 2 -2 -1 0 1 2 f = 0.00002 f = -0.004 f = 0.040 f = 0.035 (b) Iteration 2 -2 -1 0 1 2 -2 -1 0

Extension to multi-objective optimization

The algorithm can treat multi-objective optimization in many ways; an ε-constraint method transforms the multi-objective problem to multiple mono-objective ones.

Then, the algorithm could be used while storing the search tree in case one of the sub-spaces will be used by one of the remaining mono-objective problems.

Another strategy could be based on the concept of dominance. The bounds of each sub-space are now Pareto fronts, and a sub-space is eliminated if it cannot produce a better Pareto front than the actual best one. Nevertheless, some ecient algorithm has to be used for solving multi-objective problems to nd the bounds.

Extension to optimization under uncertainty

Uncertainty could be dealt with using the same algorithm by changing the bounding problems, to make them correspond to performance needed.

For example, A worst case optimization could be integrated by using a formulation from [START_REF] Deng | Optimisation robuste pour des dispositifs électromagnétiques[END_REF] p w,0 and p w,i for both problems are computed as stated in [START_REF] Deng | Optimisation robuste pour des dispositifs électromagnétiques[END_REF].

Interestingly, we have similar formulations as for the deterministic case, which enables us to treat dierent problems with the same algorithm an even to test multiple levels of uncertainty with little computational overhead since the explored sub-spaces could be reused.

Chapter Summary

In this chapter, we presented a thorough introduction to optimization using metamodels, mainly Kriging-related. We presented some of the challenges of tting "good" meta-models and how to assess their quality. The likelihood functions can sometimes lead to meta-models of lousy quality; thus, we proposed to limit the variation of the correlation parameters in a specied interval to have a minimum and maximum allowed correlation between the samples.

Furthermore, we have highlighted the drawbacks of using the conventional approach that consists of tting a single meta-model on the whole design space and enrich it sequentially using inll criteria. This approach becomes eventually unreliable since problems related to the conditioning of the correlation matrix, the complexity of tting the meta-models and the constraints handling surge.

Thus, we proposed a novel strategy that consists of building many meta-models each on a specied region of the design space, each one of them is relatively easy to t and easy to use for the purpose of optimization (less modality). Then, iteratively, prune the regions that are not promising. This process keeps a good trade-o CHAPTER 2. NON-INTRUSIVE APPROACH between exploration and exploitation of the design space; thus, it enables us to produce reliable solutions. Sometimes, this algorithm can be costly because of the B&B paradigm.

Nevertheless, the most expensive process is easily parallelizable, by either subdividing each sub-space to multiple smaller ones and explore them in a parallel manner or by taking multiple sub-spaces from the list and explore them likewise.

Chapter 3

Intrusive approach

As noted in the previous chapters, in this chapter, we will open the black-box to extract relevant information from the model for the optimization algorithms and uncertainty quantication methods. We have seen in the rst chapter that optimization algorithms based on the derivatives are the most ecient ones when the gradient information is available. Furthermore, the gradient could also be used for uncertainty quantication using perturbation methods. This double usefulness of the gradient is of high interest. Nevertheless, getting this information is not straightforward in the case of nite element analysis.

To illustrate the methodology and to point out the main diculties, we will treat a simple electromagnetic device shown in We chose to treat this problem thanks to its simple geometry, and also because of the explicit solutions of magnetic quantities can be easily calculated on the whole studied domain. This will enable us to validate the quantities computed using FEA and, most importantly, their derivatives.

The magnetic ux density inside the solenoid can be computed as follows.

B(x) = µ 0 Jd , x < R (3.1)
In the coil, we have a linear decrease (in the x-direction) of the magnetic ux density, B is written as

B(x) = µ 0 J(R + d -x) , R ≤ x < R + d (3.2)
Outside of the solenoid, the magnetic ux density is zero.

For computing B c , we take a small distance inside the coil at position x = R + e (e = 10 -3 ), then, B c (R, d, J) = µ 0 J(d + e).

The magnetic energy per unit length W is computed using the following formula

W = 1 2 1 µ 0 B 2 dx (3.3)
Therefore, the magnetic energy is reduced to the following form

W (R, d, J) = 1 2 µ 0 (Jd) 2 R + 1 3 d (3.4)
Now, we have a closed-form formula for the quantities of interest B c and W .

On the other hand, an FEA of the device is conducted, and the discrete solution is

shown in Figure 3.2. First order shape functions have been chosen to approximate the vector potential A. The exact solution of the problem is of the second order since the magnetic ux density B = curlA is of the rst order (see 3.2). Consequently, the Finite Element model can not t perfectly the exact solution and it will lead to numerical errors.

A comparison of the quantities computed with the analytic and FEA is conducted, and relative errors are then shown in Table 3 3.1 How to compute the gradient?

As noted in Chapter 1 of this dissertation, one could compute the gradient using nite-dierence by imposing a small perturbation on the variables.

nite-dierence

The nite-dierence method enables us to approximate the gradient by imposing a small perturbation on the considered variable. For example, the derivative of the energy with respect to the variable R can be approximated as follows

∂W ∂R ≈ W (R + ε, d, J) -W (R, d, J) ε (3.5)
where ε is a small value, this scheme implies running an additional FEA with the new value R + ε, generally speaking, we need as many simulations as there are variables. In our example, since we have three variables (R, d, J), we will need three additional FEA for computing the gradient.

The formula in (3.5) is called forward dierence; one could also use backward dierence or centred dierence. The latter requires twice as much number of FEA as the number of variables, but it can be more accurate than forward and backward dierence.

The choice of ε can be delicate, and some tuning is always necessary. In the literature, the estimation of the optimal values of ε is quite studied [START_REF] Iott | Selecting step sizes in sensitivity analysis by nite dierences[END_REF] [START_REF] Mathur | An analytical approach to computing step sizes for nite-dierence derivatives[END_REF]. Some prior knowledge of the model, such as an estimation of the second derivative, is often necessary for estimating a good value for ε. We can see that high values of ε (> 10 -5 ) have led to reasonably acceptable results for all the derivatives in our simple example. Nevertheless, it remains risky to use such a scheme since the optimal value of ε depends on the quantity computed (global or local) and the variables for which the derivatives are calculated.

Another behaviour that is seen in Figure 3.3 is the erratic variation when dealing with geometric variables R and d; the chaotic discontinuities are related to the change of the mesh when varying ε. In fact, when we apply the nite-dierence scheme, we have to calculate two values of the quantity of interest for two geometries (for example, geometry 1 with R and geometry 2 with R + ). The brute approach to obtaining the two FE solutions involves re-meshing the two geometries and then solving the FE problem. This approach often leads to a "discontinuity" between the two FE solutions, which are not then approximated in the same space. In Figure 3.4, we see that the same geometry could be meshed dierently depending meshing algorithm conguration. Even small variations in the geometry can lead to a change in the resulted mesh (re-meshing), which, as mentioned previously, impacts the FEA solution.

Figure 3.4: Same geometry, dierent meshes

To eliminate the impact of re-meshing, one would like to keep the same mesh topology while changing the geometry, as shown in Figure 3.5. This approach is commonly referred to as mesh morphing.

Mesh-morphing

Dierent approaches for mesh morphing have been developed in order to deform an initial mesh to take into account shape modication without changing the mesh topology. The concept consists in imposing a displacement for a set of mesh nodes and determining the new coordinates for all other ones by an interpolation approach.

In this context, the spring analogy [START_REF] Bottasso | The ball-vertex method: a new simple spring analogy method for unstructured dynamic meshes[END_REF] [84], Laplacian smoothing [START_REF] Hansbo | Generalized laplacian smoothing of unstructured grids[END_REF], linear elasticity [START_REF] Dwight | Robust mesh deformation using the linear elasticity equations[END_REF] and the Radial Basis Function [START_REF] Henneron | Mesh deformation based on radial basis function interpolation applied to low-frequency electromagnetic problem[END_REF] [88] interpolation method are largely treated in the literature. Explaining the ideas behind each method is out of the scope on this dissertation, we will limit ourselves to the application of the concept to our simple electromagnetic device.

Changing the geometry of an electromagnetic device by mesh morphing stands for moving the mesh nodes according to the change in the geometry. In our example when changing the parameter R, we proceed by moving all the mesh nodes on the boundary of the coil by ε in the x-direction while all the remaining nodes stand still.

In Figure 3.6, we show the resulted relative error in the derivative computation when using the mesh morphing. We can see that the erratic behaviour seen in Figure 3.3 has disappeared. Mesh morphing gives a piece of more reliable information about the derivatives when compared to the one given with re-meshing since its less sensitive to the value of ε. Still, it requires further attention to the choice of the best values of ε. Moreover, carrying out a mesh morphing is not necessarily evident for 3D geometries, especially in areas where the mesh is ne such as the airgap of an electrical machine.

Discussion

We have shown in this section how to compute the gradient of a quantity of interest using nite-dierence method. The nite-dierence has three shortcomings 1. The choice of the step size ε is not straightforward [START_REF] Berkani | Study on optimal design based on direct coupling between a fem simulation model and l-bfgs-b algorithm[END_REF].

2. Re-meshing error can highly deteriorate the information about the derivatives and mesh-morphing is not straightforward for most of geometries.

3. The computational cost is proportional to the number of variables which can be disadvantageous when dealing with computer expensive simulation such as FEA.

To tackle all these issues, some researchers have proposed to use the adjoint variable method for computing the gradient from the mathematical model of the electromagnetic phenomena [90] [91] [92]. This approach requires an intrusive manipulation of the FE code to calculate the gradient eciently and accurately.

The adjoint variable method is used in many elds that involve solving a large system of algebraic equations. The solution to such system is computationally expensive; thus, computing quantities of interest involving this system is also computationally expensive. The adjoint variable method can compute the gradient with a cost less or equal to the cost of solving the initial system. This last property makes the approach very attractive when dealing with many variables in the context of an FEA.

However, the implementation of such a method is not trivial and requires intrusive manipulation in the code of the simulation tool. For this reason, in this chapter, we will start by explaining the main ideas behind that adjoint variable method for a simple problem (an electrical circuit) before extending the concept to an FEA code.

Explanatory example

Let us consider a basic electrical circuit to explain the fundamental ideas of the adjoint variable method. Figure 3.7 shows an RL circuit composed of a resistor of resistance R and an inductor of reactance X driven by the sinewave voltage source V in .

For this circuit, we want to compute the Joule losses in the resistor. However, to be able to get to this end, the current owing through the circuit has to be computed since the Joule losses are equal to the resistance times the RMS current squared.

To show this, we will adopt the complex notation for the calculation.

The equation of the circuit is written as We suppose that the voltage source has only a real component (V in = v r ) and the current is written as

(R + jX)I = V in (3.6) V in R I jX
I = i r + ji i . Therefore, equation (3.6) becomes Ri r -Xi i = v r (3.7) Ri i + Xi r = 0 (3.8)
In a more compact matrix form, we write

Ku = b (3.9) where K = R -X X R , u = i r i i and b = v r 0
This system of equation has to be solved to identify u, i.e. the values of the current components, and then the Joule losses can computed by

P j = R(i 2 r + i 2 i ) = Ru 2 (3.10)
Thus, we could consider that R and X as the input variables of the model and an output P j, as shown in Figure 3.8. We denote this model as f (p) = f (p, u) = P j where p = (R, X). f and f are dierent since the rst one includes u implicitly in its denition, while f considers u as an independent variable from p. Practically, f is the quantity of interest when looking from outside of the box, while f is only seen inside the box.

u = K -1 b P j = Ru 2 R, X, v r P j Figure 3.8: Model of the circuit
For the sake of simplication, we will show the derivation of the gradient of the quantity of interest using this example, and interestingly, this derivation extends to the FE model.

Derivation of gradient

We denote the partial derivative of a function h with respect to a variable x as ∂ x h.

In the same manner, the total derivative is written as d x h.

We want to compute the gradient of the function f , it means, the partial derivatives of the function ∇f = (∂ R f, ∂ X f ).

The equation (3.11) expresses the relation between f and f

∂ p f = d p f = ∂ p f + ∂ u f d p u (3.11)
∂ p f shows the explicit dependence of the function f on the variable p while ∂ u f d p u

shows the implicit dependence of p through the variable u. We can note that if u is not dependent on p then d p u = 0. We retrieve then the fact that the total derivative is equal to the partial derivative.

To compute the gradient of the function f , we need to compute the total derivatives of f , i.e. d p f . In the following derivation, we will only consider f rather than f .

We rewrite the equation (3.9) as

g(u) = g(p, u) = Ku -b = 0 (3.12)
The dependence of u on the variable p is expressed in the equation g(p, u) = 0, Thus, the last term d p u of (3.11) is computed by deriving this equation

d p g = ∂ p g + ∂ u g d p u = 0 (3.13)
Then, d p u is the solution of the equation

d p u = -(∂ u g) -1 ∂ p g (3.14)
By replacing d p u in (3.11), we get

d p f = ∂ p f -∂ u f (∂ u g) -1 ∂ p g (3.15)
This formula enables to fully determine the gradient of the function f versus partial derivatives of the function f and the equation g.

For the example above f (p = (R, X), u) = Ru 2 and g(p = (R, X), u) = Ku -b, the quantities forming the gradient are written as

∂ u f = 2Ru (3.16) ∂ R f = u 2 (3.17) ∂ X f = 0 (3.18) ∂ u g = K = R -X X R (3.19) ∂ R g = ∂ R Ku= 1 0 0 1 u (3.20) ∂ X g = ∂ X Ku= 0 -1 1 0 u (3.21)
Thus, the derivative can be written as

d R f = u 2 -2Ru K -1 1 0 0 1 u (3.22) d X f = -2Ru K -1 0 -1 1 0 u (3.23)
We can notice that we have a general formula to compute the derivative for any value of R and X. The derivatives for each variable require solving a system of linear equations of the initial problem as we can see in (3.22) and (3.23) where we have to solve for X R and X

X KX R = 1 0 0 1 u (3.24) KX X = 0 -1 1 0 u (3.25)
Solving the equations is the most expensive operation; thus, the cost for computing the gradient is proportional to the time needed multiplied by the number of variables p. The adjoint variable method enables to tackle this drawback by solving the equations only once.

The adjoint variable method

As noted in the last section, the gradient of the function f is written as

∂ p f = ∂ p f -∂ u f (∂ u g) -1 ∂ p g (3.26)
By inverting the order of multiplication of the terms above, we get

∂ p f = ∂ p f -∂ u f (∂ u g) -1 ∂ p g (3.27)
The expression between the brackets does not depend on the parameters p; it depends only on the solution u of the equation g(p, u) = 0. Thus, this operation can be done once and be applied to all the variables p. We introduce

λ = -∂ u f (∂ u g) -1 (3.28)
where λ is called the adjoint variable. Then, the previous expression becomes

∂ p f = ∂ p f + λ∂ p g (3.29)
To summarize, for computing the gradient of a function f , the following steps has to be conducted :

1. Compute the partial derivatives for u :∂ u f and ∂ u g

2. Solve the linear system: λ ∂ u g = -∂ u f for λ 3. Compute the partial derivatives for p : ∂ p f and ∂ p g

4. Calculate the adjoint using 3.29

Step 1 is relatively easy to conduct. In fact, the function f is calculated from the solution u of the equation g(u) = 0; thus, its derivative should be no harder than calculating the function itself. On the other hand, the partial derivative of g with respect to u is sometimes already provided by the original model. If the equation g(u) = 0 is a linear system of equations, then its derivative is equal to the matrix of this system. Otherwise, if the equations g(u) = 0 are nonlinear, ∂ u g = ∂ u g is usually computed to solve the equation using an iterative process, e.g. Newton-Raphson, then the derivative is the Jacobian which are often available or at least easy to determine.

Step 2 is the most computationally expensive since it requires solving a system of linear equations of the same size of the system of equation dened by g(u) = 0. It is worth noting that even if g(u) = 0 is a nonlinear system of equations, the problem to solve in step 2 is linear. Therefore, solving this system will always be less expensive than solving g(u) = 0, That is to say, that the cost of computing the gradient will be equivalent or less than solving g(u) = 0 for any number of variables. This property makes the adjoint variable method very attractive when dealing with many variables.

Step 3 requires the calculation of the partial derivatives of the f and g with respect to p. ∂ p f calculation is done as for ∂ u f . Nevertheless, ∂ p g may not be straightforward since it involves the derivatives of the system of equations. Sometimes this link between the variables is not as obvious as it was in our example.

Step 4 is straightforward since the total derivatives is given by 3.29 and calculated from the product of terms determined in the previous steps.

Adjoint variable method for FEA

The development of the adjoint variable method for FEA has started in the eld of structural engineering, mainly in topology optimization for reducing the volume of devices [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF]. Afterwards, some researchers extend its usage for shape and topology optimization of electromagnetic devices [START_REF] Park | Design sensitivity analysis for nonlinear magnetostatic problems using nite element method[END_REF] [92] [START_REF] Mohamodhosen | Topology optimisation of electromagnetic devices[END_REF].

Despite the attractiveness of the method, to the author's knowledge, no commercial software dedicated to electromagnetic design has developed the method. At the same time, it has been present in almost all CFD and structural mechanics software for almost two decades.

In the following section, we will explain how the adjoint variable method can be implemented in a magnetostatic FEA code. The implementation can be extended to solve any other electromagnetic formulation derived from the Maxwell equations.

Finite element model

As detailed in Chapter 1, Maxwell equations are solved using FEM, and the derivation ends by solving a discrete system that reads as follows:

Find u : Ku = b (3.30)
where K is a sparse matrix, u is the state variable and b is the source term dened as follows (see 1.80 and 1.81)

K ij = T ∈T h g 1 µ ∇α g i J -1 T • ∇α g j J -1 T |J T |w g (3.31) b j = T ∈T h g J s α g j + M • ∇α g j J -1 T |J T |w g (3.32)
α g j and w g are, respectively, the basis functions evaluated in the Gauss quadrature points and their weights dened in the reference element, µ is the permeability of material, J s source current density, M is the magnetization of magnets and J T is the Jacobian of the mapping from the reference element to an element T (|J T | is determinant of J T ).

Additionally, we consider a quantity of interest f (p) for which we want to compute the gradient.

As for the electrical circuit, we dene two functions f and g for FEA

g(p, u) = K(p, u)u -b(p) = 0 (3.33) f (p) = f (p, u) (3.34)
This representation of the problem will simplify the derivation of the gradient using the adjoint variable method as for the electrical circuit; thus, we write

∂ p f = ∂ p f + λ ∂ p g (3.35)
where λ is the adjoint variable and is solution of the linear system of equations

λ∂ u g = -∂ u f (3.36)
The solution of the linear system 3.36 is the most computationally expensive operation, and it is independent of the variables p. Thus, the gradient computation cost is independent of the number of variables.

Computing the gradient using the adjoint method is reduced to compute the partial derivatives with respect to state variable u (∂ u f , ∂ u g) and those with respect to design variables p (∂ p f , ∂ p g).

Derivatives for state variables

Derivatives with respect to the state variables are relatively simple to express. ∂ u g is generally always computed in the initial FEA code.

If there is no nonlinear ferromagnetic material (∂ u µ = 0)

∂ u g = K (3.37) Otherwise (∂ u µ = 0) ∂ u g = ∂ u (Ku) (3.38)
For solving g(p, u) = 0, the Jacobian matrix ∂ u g is necessary for a nonlinear solver such as Newton-Raphson. Thus, this quantity could be retrieved without any additional eort.

On the other hand, ∂ u f can be computed eciently because the function f is given by an explicit expression of the state variable u. For example, if we consider that the quantity of interest is the magnetic energy that is written, as shown in

Chapter 1 f (p, u) = W = 1 2 u b, (3.39)
then, the derivative with respect to u is

∂ u f = 1 2 b , (3.40)
For other quantities of interest f , ∂ u f could be derived and implemented with ease. lation of derivatives of geometry with respect to design variables could be challenging within commercial CAD software.

The derivative of the quantities f and g with respect to a geometric design variable p can be written as follows

∂ p f = N i=1 ∂ x i f d p x i + ∂ y i f d p y i (3.41) ∂ p g = N i=1 [∂ x i g d p x i + ∂ y i g d p y i ] (3.42)
where (x i , y i ) for i = 1, . . . , N are the node coordinates used for the FEA, by writing the quantities in this form, we decouple the dependence the quantities related to the FEA from the geometry denition. Furthermore, we decouple the computation cost from the number of design variables since ∂ x i f , ∂ y i f , ∂ x i g and ∂ y i g are independent of p.

These quantities can be computed analytically. ∂ x i g and ∂ y i g can be calculated based on (3.31), (3.32) and (3.33). In fact, the only quantity that depends on the mesh coordinates is the mapping from the reference element to an arbitrary element. Thus, only the derivative of the mapping is needed for the computation of the quantities.

We have

∂ x i g = ∂ x i [Ku -b] (3.43) ∂ y i g = ∂ y i [Ku -b] (3.44)
Moreover, as detailed in Chapter 1, K and b are computed using the reference element concept, as shown in Figure 3.9. Equation (3.45) shows how the reference element coordinates (x r , y r ) are transformed into the coordinates (x, y) of an arbitrary element using the mesh nodes (x 0 , x 1 , x 2 , y 0 , y 1 , y 2 ) dening this element.

J T = x 1 -x 0 y 1 -y 0 x 2 -x 0 y 2 -y 0
The derivative of the mapping with respect to x 0 , x 1 , x 2 , y 0 , y 1 and y 2 is easy to compute; hence, ∂ x i g and ∂ y i g can be assembled as done for g. Now, let us go up to the other part of equations (3.41) and (3.42): d p x i and d p y i . These quantities somewhat link the variation in the node coordinates versus the change in geometry. Hence, they represent the real link to the CAD model after the latter has been meshed. In this dissertation, we will show how to compute these quantities without intrusive manipulation, neither in the CAD software nor in the mesh tool.

Geometric parameters are related to one or many edges that separate dierent regions (faces), as shown in Figure 3.10. This numbering is very useful for FEA since it enables to impose external conditions. For example, we used E2 to impose a Dirichlet boundary condition, and the face F 3 is used for imposing a current density J. We use this numbering for calculating the desired quantities. In this example, we have two geometric design variables R and d, as shown in Figure 3.10. These variables enable the positioning and sizing of the coil modelled by the face F 3; thus, they mainly act on edges surrounding the coil, namely E3, E4, E6 and E9. Hence, d p x i and d p y i are computed only for the mesh nodes that lie on these edges and not on all the other nodes. Therefore, we write

R d

d p x i =    d p x i , if node i is on E3, E4, E6 or E9.
d p y i =    d p y i , if node i is on E3, E4, E6 or E9. 0, otherwise. (3.47)
The edges E3, E4, E6 and E9 dene a rectangle. The objective is to dene the rectangle contour using a parametric equation as follows

P (t) = (x(t), y(t)) (3.48)
In general, one could propose many ways to parameterize a shape. For the rectangle, we could parameterize each edge independently and compute the derivative or use a single parametrization for the whole rectangle as below

x(t) = R + d 2 [1 + sign(cos(t))] y(t) = h 2 [1 + sign(sin(t))]
and t ∈ [0, 2π[.

The derivatives of x(t) and y(t) with respect to the design variables are computed with ease, d R x = 1

(3.49) d R y = 0 (3.50) d d x = 1 2 [1 + sign(cos(t))] = x -R d (3.51)
d d y = 0 (3.52)
Equations (3.49)-(3.52) are evaluated using the node coordinates previously identied, then, used in equations (3.41) and (3.42) to compute ∂ p f and ∂ p g. Now, we have everything for the computation of the gradient of a quantity of interest with respect to geometry variables using the adjoint variable method.

In this paragraph, we computed the shape sensitivities d p x i and d p y i for a rectangular shape using a parametric equation. This framework is extendable to dierent types of geometry and design variables. In appendix A, we derive the shape sensitivities for some recurrent types of geometries.

B Physical variable

The physical variables that can be considered in electromagnetic modelling are the permeability of ferromagnetic materials, the coercive eld of magnets, and the imposed current density. This kind of modelling is commonly used in the context of topology optimization, where the objective, in general, is to reduce the volume of material present in a device. Topology optimization was rst introduced in the context of structural optimization before being extended to other physical elds such as heat transfer, uid dynamics and electromagnetism.

In this dissertation, the objective is to derive the formula that enables to compute the gradient for parameters that act on the physical properties of electromagnetic devices.

Permeability

If a parameter p controls how µ varies, we can compute the derivatives ∂ p f , ∂ p g can be computed :

∂ p g = ∂ p [Ku -b] = [∂ p K] u (3.53) with ∂ p K ij = K ij -∂ p µ µ .
The quantity ∂ p µ highlights the sensitivity of µ with respect to the design variable p. In the literature, we can nd multiple ways of modelling. For example, the density method uses a polynomial mapping for the parameterization [START_REF] Mohamodhosen | Topology optimisation of electromagnetic devices[END_REF].

µ = µ air + (µ iron -µ air )p n (3.54)
with µ air and µ iron are respectively air and iron permeabilities, A3. This method enables us to chose where in the studied domain to put iron or air. In this case, ∂ p µ is written as

∂ p µ = n(µ iron -µ air )p n-1 (3.55)

Coercive eld

As for the permeability, if p controls the magnetization of a permanent magnet, one could, in a similar manner, deduce the derivatives

∂ p g = ∂ p [Ku -b] = -∂ p b (3.56) with ∂ p b j = T ∈T h g ∂ p M • ∇α g j J -1 T |J T |w g .
The quantity ∂ p M is computed based on the design variable; p can be the magnitude, the direction of the coercive eld, etc.

Imposed current density

Following the same principle, the derivative with respect to a variable that controls the current density is

∂ p g = ∂ p [Ku -b] = -∂ p b (3.57) with ∂ p b j = T ∈T h g ∂ p J s α g j . |J T |w g .

Discussion

In the precedent section, we showed the derivation of the gradient of a quantity of interest from a nite element code using the adjoint variable method. It is worth noting that two approaches exist to calculate the adjoint: the continuous approach and the discrete approach. The continuous approach aims to compute the gradient from the governing equations (Maxwell equations) before discretization, while the discrete approach calculates the gradient from the discrete equations. We focused and developed the discrete approach for many important reasons :

1. Relatively simple implementation (but tedious)

2. Implementation of the derivatives of each subroutine/process individually in analysis code 3. Build up larger components by chain rule dierentiation of analysis code 4. Ability to check the derivatives by nite-dierence while debugging

The discrete adjoint approach can be implemented in a modular fashion using the same data-structures/solution strategy as analysis as shown in Figure 3.11. The adjoint is composed of three main steps, at parameter setting, one should compute the sensitivities with respect to the design variables (geometric and physical), then, the assembling of the quantities

FE code

Adjoint code Start

∂ u f , ∂ u g, ∂ x i f , ∂ x i g, ∂ y i f , ∂ y i g and the 3.3. ADJOINT VARIABLE METHOD FOR FEA
ones related to the physical variables is conducted, afterwards the adjoint problem is solved, and nally, the gradient of the quantity of interest is computed as

∇f = ∂ x i f + λ ∂ x i g d p x i + ∂ y i f + λ ∂ y i g d p y i + ∂ µ f + λ ∂ µ g d p µ + ∂ M f + λ ∂ M g d p M + ∂ Js f + λ ∂ Js g d p J s
In this equation, one can see that all the dependencies on the design variables p are outside of the most computationally expensive operations, namely the solution of the adjoint problem and the assembling of the dierent quantities. The independence of the number of design variables makes the adjoint variable method very attractive for treating problems with large numbers of design variables and few numbers of objectives.

To show the utility and eectiveness of the adjoint variable method, in the next section, we treat some examples to highlight the advantages of the adjoint variable approach over the nite-dierence in computational cost and precision.

Examples

Solenoid model

We consider the device treated at the beginning of the chapter, where we have made a comparison between analytic and FEA computation of the quantities of interest.

We conduct the same comparison on the gradient of the quantities of interest as shown in Table 3 One can notice that the gradient of B c and W computed by the adjoint variable method are corresponding to the ones computed analytically with a relative error of the same magnitude as for the quantities B c and W .

TEAM workshop 22

This example is considered as one of the benchmarks for optimization problems of electromagnetic devices. The objective here is compare the gradient computed nite-dierence to the one computed using the adjoint variable method (further details in Chapter 4). To validate the eectiveness of the adjoint variable method, the approach is tested against the gradient computed using centred dierence (CD) 1 using mesh morphing technique and with remeshing technique. The mesh morphing (MM) technique enables "small" changes of the shape while maintaining the topology (node connectivity) of the mesh. The aim of using such a technique is to reduce the errors due to re-meshing in the comparison (see example treated at section 1).

The results of the comparison are summarized in Table 3 For nite-dierence schemes, we tried dierent values of ε, ranging from 10 -10 to 10 -2 . In the table, only the values that correspond to the best results in terms of relative error are shown.

When using mesh morphing (CD w. MM ), the gradients are coherent in less than 1% error relative to the gradient computed with the adjoint with a speedup of 3.2. This speedup is attained thanks to the main propery of the adjoint variable method that solves fewer equations. However, when not using mesh morphing (CD w/ MM), which means generating a new mesh for each geometry change, re-meshing errors appear to perturb the gradient information highly (relative error higher than 150%).

It is worth noting that the dierence in time for CD w. MM and CD w/ MM is related to the overhead of regenerating a new mesh for all geometry variations. In general, mesh morphing is not always a simple task; in this example, it was possible to morph to mesh nodes to correspond to the geometry change thanks to the simple shapes involved in the modelled device (rectangular regions). Thus, without using mesh-morphing, the nite-dierence method may lead to signicant errors when compared to the adjoint variable method. Furthermore, a speedup from 5.5 was attained for the computation of the gradient, which conrms the eciency of the approach. 1 dOF dpi (p) ≈ OF (p1,...,pi+ε,... )-OF (p1,...,pi-ε,... ) 2ε 3.4. EXAMPLES

Chapter Summary

In this chapter, we have presented the adjoint state variable method and how it can be applied to a FEM code. We have developed an ecient way to compute the derivatives of the shape sensitivities, which are vital for the computation of the gradient. Sometimes, nite-dierence could be used with a mesh-morphing strategy for computing the gradient but it can have redhibitory eect since the mesh displacement is not always straightforward and can lead to non-conforming mesh.

The adjoint variable method was compared to the nite-dierence one to validate and highlight its eectiveness in terms of precision and computational time.

Geometric parametrization of shape variables is still one of the shortcomings of the method. We presented an approach based on the parametric equations of the geometric shapes, however, for very complex shapes, this can be very cumbersome.

In the following Chapter, the gradient computed using the adjoint variable method is provided to a derivative-based optimization algorithm for the analysis of two benchmarks from the literature and comparison of the performances to the nonintrusive approach detailed in Chapter 2.

A FEM model

We have created a model for the benchmark with our FEA tool (MagFEM); as the model is symmetric with respect to the r-axis (the antisymmetric case is somewhat dierent than the 2D case. In appendix B, we detail the corresponding formulation).

In this part, we discuss some aspect of the model mainly related to the parameterization and the FEA.

Parameterization

The geometry of the device is simple; it is composed of three rectangles, as shown in

(R 1 , R 2 , h 1 /2, h 2 /2, d 1 , d 2 )
are used for the parameterization of the geometry, while the current densities (J 1 , J 2 ) are physical variables that are parametrized in the FE model.

Simulation

In the FEA, the geometry is meshed as shown in Figure 4.3 on the left. A ne mesh is imposed around the coils and coarse one far from them. On the right of The constraint (4.2) aims to limit the maximal ux density (|B|) in the coils to ensure the quench condition [1]. As shown in The second constraint aims to prevent both coils from overlapping. Unfortunately, optimization algorithms can sometimes during iterations violate some con-X X X As we are comparing dierent approaches, we chose to dene another optimization problem equivalent to the initial one that avoids the issue mentioned above.

We dene a variable A 2 as

A 2 = R 2 -R 1 -(d 1 + d 2 )/2
that replaces R 2 . Using interval arithmetic, we compute the bounds of this new variable [-3, 3.9]. By using the constraint (4.3) that is preventing the two coils from overlapping, we impose A 2 to be higher than 0.001.

R 1 + A 2 + (d 1 + d 2 )/2 ≤ 5 (4.8) -R 1 -A 2 -(d 1 + d 2 )/2 ≤ -1.8 (4.9)
where p are the design variables p = (R

1 , A 2 , h 1 /2, h 2 /2, d 1 , d 2 , J 1 , J 2 ).

TEAM Workshop Problem 25

In this test case, the device is used to orient the magnetic powder and produce anisotropic permanent magnets [2]. The magnetic powder is inserted in the cavity.

The orientation and strength of the magnetic eld should be controlled in order to obtain the required magnetization. A coil creates the magnetic eld. The current density is xed to 1.239219A/mm 2 . The die press and electromagnet are made of steel with a non-linear permeability. The geometry of the whole device is shown in

Simulation

In the FEA, the geometry is meshed as shown in Figure 4.6 on the left. A ne mesh is imposed around the moulds and coarse mesh far from them. At right, in the same gure, the distribution of the ux density is shown.

B Optimization problem

The objective of the shape optimization is to obtain a ux density that is radial in the cavity space and with a constant magnitude of 0.35T.

Algorithms settings

In this section, we present the conguration of the optimization algorithms. A set of parameters are dened based on our experience from using each algorithm.

B2M2 algorithm

The branch and bound based meta-model (B2M2) algorithm is detailed in Chapter 2, Section 2. Some options need to be congured, mainly the depth level and the number of FEA evaluations needed for tting the meta-models.

Initial design : The initial design of each sub-space is created using a spacelling LHS of size two times the number of variables. This may seem few, but as the algorithm progresses, the sample points of a parent space will be added to the initial design of its sub-spaces. Thus, the sample will grow as the sub-division of the space continues.

Inll points : Additional sample points are added in two dierent manners;

the rst one aims to further explore the design space by adding points where the 4.2. ALGORITHMS SETTINGS prediction error is high, while the second one enables the exploitation of the best solution by adding new sample points that maximize the expected improvement ISC.

The number of sample points to be added by each ISC (inll sampling criterion) is set to twice the number of variables. However, for expected improvement, it may be stopped before this limit if there is no improvement in sampling additional points.

Globally, in each sub-space, there will be sampled at most six times the number of variables.

Stopping criteria : The algorithm is set to stop if it attains a predened depth level. The maximum depth level is set to 10.

SQP algorithm assisted by adjoint variable method

We use SQP algorithm from the Matlab Optimization toolbox with the gradient of the quantities of interest computed using the adjoint variable method.

Multi-start : To cope with the local search of the SQP method, we use a multi-start strategy. We perform 100 runs with dierent initial points; These points are sampled using an LHS design to obtain a uniform distribution on the whole design space.

Stopping criteria : The stopping criterion is based on the StepTolerance option, meaning iterations end when the relative change in the solution is less than 10 -8 .

DIRECT

We use DIRECT implementation of [START_REF] Finkel | Direct optimization algorithm user guide[END_REF]. The default options are maintained except for the number of iterations.

Constraints : The constraints are handled by penalty method by minimizing the penalized objective function

OF (p) + i λ max(0, g i (p))
where λ is the penalty factor and has been chosen to be equal to 10.

Stopping criteria : As the algorithm does not have an appropriate stopping criterion except the number of iterations and the number of evaluations. We stop the algorithm manually when there is no improvement in the objective function in 50 successive iterations.

Genetic algorithm

We use GA from Matlab Global Optimization toolbox [START_REF]The MathWorks, Matlab, Global Optimization Toolbox[END_REF].

Multi-run : As the randomness inherent in genetic algorithm, meaning solving twice the same problem may not lead to the same solution, we adopted a multi-run strategy by performing 50 runs of the algorithm.

Population : The population options are set as follows Constraints : The constraints are handled by the penalty method as dened in [START_REF] Deb | An ecient constraint handling method for genetic algorithms[END_REF].

Stopping criteria : The algorithm is stopped if the value of the tness function did not change during 50 generations.

Results and comparison

Comparison protocol

The comparison of the performance of dierent algorithms will be based on two criteria: the cost and the quality of the solution.

The quality of the solution stands for the value of the objective function found.

The smaller, the better.

On the other hand, the cost of the optimization stands for the number of evaluations needed to attain the solution. As we compare dierent approaches, we propose a probabilistic metric to measure the cost.

The convergence rate (CR) is the percentage of the number of runs that converged to the best solution. The convergence rate concerns mainly SQP and GA since we perform multiple runs and the number of runs needed is unknown beforehand. To compute CR, we look, among the solutions of each run, for the best one (smallest 4.3. RESULTS AND COMPARISON 123 value), then we count the number of runs that converged to that solution within a tolerance of 1% of the design range.

Once the convergence rate is determined, we write the probability that at least one optimization lead to the best solution is

P = 1 -(1 -C.R.) n
where n is the number of optimization runs to be considered.

If we want to consider a high probability -(P = 99.73%) of nding the best solution, we can compute the expected number of runs needed n n = ln(1 -P ) ln(1 -CR)

(4.11)
Then, the expected total number of evaluations (Expected # FEM evals) is n times the total average number of evaluations performed.

Expected # FEM evals = n average(# FEM evals) (4.12)

A particular property related to SQP is the computation of the gradient using the adjoint variable method. The cost of computing the gradient is equivalent to evaluating the FEA once; thus, the # FEM evals are doubled for the SQP method.

TEAM Workshop Problem 22

The TABLE 4.3 summerizes the optimization results We notice that SQP is largely outperforming the other approaches in term solution quality, actually, 17% of the solutions found by this approach has an objective CHAPTER 4. APPLICATIONS AND BENCHMARKING function value less than 0.05 while being distinct in term of the variables, this explains the low CR of the method and, thus, the high expected number of evaluations.

The genetic algorithm performed better than DIRECT in terms of quality of the solution, but the former somewhat suers from its inherent randomness and thus the low convergence rate.

B2M2 approach was able to get a solution of better quality than GA and DI-RECT but less than SQP.

Discussion

TEAM Workshop Problem 22 is one of the most treated benchmarks for the optimization of electromagnetic devices. In the literature, researchers have reported their results in many papers. However, due to the dierences in the FE solver used by each researcher, the comparison to their results is not possible.

It is worth noting that some results from literature were taken into account and compared to ours [START_REF] Alotto | Stochastic algorithms in electromagnetic optimization[END_REF]. Nevertheless, those from literature, when evaluated in our FE model, have bigger values of the objective function .

Another characteristic that can be deduced from this problem is the multimodality. The low convergence rate of SQP highlights this fact, indeed, the solutions found tend to have similar values of the objective function while having design variables that are entirely dierent.

The SQP assisted the adjoint variable method was able to get the best solution.

However, due to multi-modality, the cost of obtaining this solution is high when compared to the one found by B2M2, which is of inferior quality and with a lower cost.

Two facts can explain GA and DIRECT failure; rst, the multi-modal behaviour of the problem treated, which rendered locating the best solution very delicate, and secondly, the constraint handling strategy. Indeed, the penalty approach used may not be best-suited, and the penalty factors may not be very adequate. Thus, the tuning of these factors will induce additional cost.

TEAM Workshop Problem 25

The results are summarized in In terms of solution quality, SQP and DIRECT methods lead to the best solutions (smaller value of W ), while SQP has a cost of almost seven times less than DIRECT.

B2M2 approach led to a competitive solution, even though the value of the objective function is a little bit higher. However, the design variables are very close to the ones found by SQP (dierence of less than 1% of the design range). Most importantly, B2M2 cost is the most eective.

Chapter Summary

This Chapter presents a comparison between intrusive and non-intrusive approaches for the optimization of electromagnetic devices using FEM. We treated two wellknown benchmarks from the literature [1] [2]. Then we used two metrics for the comparison, the rst one is the quality of the solution, and the second is the computational cost.

The choice of the algorithm from each category is based on what we can generally nd in the literature, the choice of a particular algorithm or implementation was

based on what we are working on, i.e. SQP and meta-model approaches, and the availability and simplicity of usage. These are usually the challenges that designers are facing when doing optimization.

The genetic algorithm (GA) is one of the most used algorithms to deal with noisy data. For both test cases, GA did not perform well. These performances might be slightly improved by doing some parameter tuning or using other implementations.

DIRECT performed very well for the second test case, but it was ineective for the rst one; this can be explained by the constraints handling strategy, the penalty method implemented in the algorithm may not be the best, and additional tuning is required.

In terms of performances, the SQP approach outperforms other strategies for both test cases; this was possible due to the computation of the gradient using the adjoint variable method. This improves the quality of the solutions drastically, but this comes with the expense of intrusive manipulation of the FEM code.

The B2M2 approach remains a good alternative in terms of implementation and solution quality. The developed approach was able to overcome some of the very well-known issues when using meta-models for optimization.

Chapter 5

Claw-pole machine

This last Chapter rst presents the main characteristics of the claw machine, the subject of this study. This machine diers from conventional electrical machines, in particular by the geometry of its rotor, which generates a three-dimensional magnetic ux and by its mechanical structure as the stator is enclosed in two brackets.

In this Chapter, we will focus on another aspect of the thesis related to the reliability of an electrical machine but closely related to the other aspects. The manufacturing processes, due to their imperfections, lead to dispersion in the material properties and the dimensions of the products. We focus on the impact of imperfections on the performance of the machine. The study presented here enables to characterize the imperfections and to conduct, for example, a robust optimization using the tools proposed in the previous Chapters.

The dispersion related to the manufacturing process induce "variability" on quantities of interest, which must be chosen among the preponderant performance factors for the product. For example, for a claw-pole machine (automotive alternator), the inuence of the dispersion on the output current or acoustic response must be studied. The aim is obviously to ensure that the product complies with the specications that are contractual between the supplier and the customer. Clawpole machines are commonly used as automotive generators due to their simplicity and low manufacturing cost. This machine is produced in mass. Deviations on the dimensions of the machines parts and also on the parts relative position versus the nominal values can occur due to the imperfections of the manufacturing and assembling processes. Moreover, these deviations are also subject to dispersions in mass production due to the variability of the processes with time.

The classical approach based on a deterministic model (the model inputs are entirely known) to address this problem is no longer sucient. It is necessary to consider the input quantities as uncertain quantities and to place oneself within the conceptual framework of uncertainties [START_REF] Liu | Study of the inuence of the fabrication process imperfections on the performance of a claw pole synchronous machine using a stochastic approach[END_REF]. The approach to uncertainties introduced by Taguchi and extended to various uncertainty propagation methods (as shown in Chapter 1) is today commonly used in dierent engineering domains.

In this Chapter, we present the claw-pole machine and its operating mode then introduce the metrology procedure adopted for the measurement of the electrical machine. These measurements are analyzed, and the variability of the manufacturing process is characterized. Afterwords, a parametrized nite element model of the machine is presented. Next, the uncertainty propagation using the variability characterized and the model of the machine is conducted employing Monte Carlo simulation assisted by a Kriging meta-model.

Electrical machine

The claw-pole alternator is a synchronous electrical machine with a wound rotor which aims to supply the on-board network and charge the vehicle battery. It is placed on the combustion engine's accessory side and is driven by a belt connected to the crankshaft pulley, as shown in Figure 5.1. The ratio between the rotation speed of the alternator and that of the combustion engine is generally between 2.5 and 3. The diode bridge and its heatsink, as well as the brushes and the regulator, are attached to the rear bearing. A plastic cover is tted to protect these components.

Therefore, the alternator is a complicated assembly comprising parts whose geometries are also complex, such as the rotor. The role and features of each part are The stator winding can be single, with one slot per pole and per phase, or double with two slots per pole and per phase. The stator phases can be star or delta coupled. Delta coupling is sometimes preferred since it does not require the connection of a neutral point. However, there is a third harmonic current owing in the delta in relation to the electrical frequency in contrast to star connection. This current generates additional Joule losses.

After the stator winding operation, the stator is impregnated with a varnish in order to improve the mechanical strength of the winding and to improve the thermal conduction between the stator and the winding (Figure 5.3).

B The rotor

The rotor consists in two main parts; the claws that form the machine's magnetic poles and an excitation coil. The machine contains 12 claws (6 in the front and 6 in the rear), each claw corresponds to a pole. These steel parts, of type SAE 1005, are obtained by forging. Surrounded by the claws, a coil of copper wire creates a magnetic eld in the axial direction (i.e. the direction of the axis of rotation). This magnetic eld is directed by ferromagnetic steel to the stator core.

C The brackets

The brackets refer to the parts on each side of the stator. The geometries of the two bearings ("front" side pulley and "back" on the heatsink side) are generally dierent (Figure 5.5). They ensure the containment of the stator by enclosing it in the axial direction (so-called "sandwich" conguration) with the help of four screws. The brackets also enable the alternator to be mounted on the combustion engine using mounting brackets. The brackets are die-cast aluminum alloy (EN AC-46200).

D The brushes and regulator

The excitation current of the rotor coil is supplied via two brushes (Figure 1.9), sliding on two rings placed on the shaft (see Figure 1.6). The brush holder is

E Bridge rectier

The bridge rectier enables to transform, using a diode bridge ( In their specications, car manufacturers generally require a minimum output current and eciency for dierent RPM (typically 1800, 3000, 6000 RPM). Various electrical tests for alternators are dened in ISO 8854 and internal Valeo specications.

Variability of output current in claw-pole machine

The curve shown in Figure 5.7 is obtained by measuring a machine issued from the production line. For our study, we collected three batches; each batch containing ten machines, and then we measured the outputs current. The variability in the measured current can attain 4 A; this leads to more than 5% of deviation at low speed. This variability is explained by the imperfection of the manufacturing process. In the next Section, we discuss the impact of the manufacturing process on the electrical machines' performances. 0 1,000 2,000 3,000 4,000 5,000 6,000 7,000 0

Impact of manufacturing process

There is a great deal of research into the design and optimization of electrical machines. New machine structures are emerging as well as new procedures to optimize their performance. Thus, the claw-pole machine's performance is analyzed and optimized by analytical and/or nite element methods. The inuence of magnetic materials on the behaviour of the claw machine has also been studied. However, it is noted that, in general, the manufacturability and imperfections of the electrical machine are not, or only to a limited extent, considered in most of these studies. Nevertheless, in the literature, there is some research on the impact of process imperfections and their eects, which we will summarize in the following Section.

A Impact on material properties

Dierent studies were conducted on the impact of the manufacturing process on the claw-pole machine's material properties. Some of the latest are studying the magnetic properties of the stator and the rotor [101] [102]. In [START_REF] Youssef | Slinky stator: The impact of manufacturing process on the magnetic properties[END_REF], El Youssef et al.

classify the impact of dierent operations in the manufacturing process (straightening, punching, rolling, ...) on the magnetic properties of the stator. They characterize qualitatively how each operation impacts the permeability and the iron losses of the stator. On the other hand, Borsenberger et al. highlighted how to measure the magnetic properties of a massive material (such as the rotor of the claw-pole machine) [START_REF] Borsenberger | Methodology for the analysis of the impact of the forging parameters on metallurgy and mechanical properties in case of solid electromagnetic manufactured parts[END_REF]. They revealed that the forging operation leads to disparities within the magnetic properties within the rotor.

Most manufacturing processes deteriorate the magnetic properties of materials. However, it may not always be the case; some can be benecial, El Youssef et al.

showed that compacting can improve the magnetic properties.

B Impact on geometry

The impact of the manufacturing process on the geometry of the machine is less addressed in the literature except on some limited studies dedicated to eccentricities. The tolerances induced by mass production can lead to variabilities on the performances of the electrical machine. In [START_REF] Liu | Prise en compte des incertitudes dimensionnelles introduites par les procédés de fabrication dans les modèles numériques de machines électriques[END_REF], S. Liu has applied a statistical approach for quantifying the impact of variability on the electromotive force and the torque of the claw-pole machine. He used a parametric nite element model with an approximation method (sparse approximation) to highlight the impact of geometry variability on the performances of the machine.

Summary

Our work can be seen as a succession of S. Liu Ph.D. thesis; We are interested in the geometry and its impact on the electrical machine's performances. We aim at studying how the manufacturing process induces variability on the electrical machine by performing metrological measurement series on machines withdrawn from the production line, then, using the metrology data, to model the variability using random variables. Afterwards, these variables are used in collaboration with a parametric nite element model of the machine to perform an uncertainty propagation using an approximation method to quantify the inuence of the variability of the geometry on the quantities of interest.

In the following Sections, we start with the metrology protocol to highlight how the electrical machine geometry is measured. Afterwards, we proceed by the variability modelling presentation. Then we present the parametric nite element model before explaining how we use the Kriging approximation for the uncertainty propagation to quantify the inuence of the variability of the geometry on the quantities of interest.

Metrology

In this Section, the aim is to measure the geometry of the claw-pole machine manufactured by VALEO. Metrology is performed using a Coordinate Measuring Machine (CMM). This device measurements the geometry of physical objects by detecting discrete points on the surface of the object using a probe. Dierent types of sensors are used in CMMs, including mechanical, optical, laser, ... In our case, we use a ZEISS UPMC CARAT machine, which has a mechanical probe (as shown in 5.9), The uncertainties on our machines due to this CMM are less than a few micrometers.

The main objective is to metrologically determine the geometry of air-gap zone between the stator and the rotor. The air-gap plays a key role in electromagnetic behaviour; the nominal size of the air-gap is equal to 325µm, which is inaccessible The position of the probe can be controlled manually or using a computer. CMMs typically specify a probe's position as a function of its movement relative to a reference position in a three-dimensional Cartesian coordinate system (XYZ axes). In addition to move the probe along the X, Y and Z axes, many machines also allow the angle control to measure surfaces that would otherwise be inaccessible. Figure 5.10 shows how the measurements are performed; indeed, the contact point is computed mathematically from the measured point in the CMM scale and the direction of the contact force. A surface can be reconstructed from limited measured points.

Methodology

As noted before, the machines are disassembled to measure the air-gap. However, some referencing elements are kept to virtually reassemble the measured data and then fully characterize the air-gap. This methodology was initially developed during S. Liu thesis in collaboration T. Coorevits of the MSMP lab [START_REF] Liu | Inuence of the stator deformation on the behaviour of a claw-pole generator[END_REF] [46] [START_REF] Liu | Prise en compte des incertitudes dimensionnelles introduites par les procédés de fabrication dans les modèles numériques de machines électriques[END_REF]. We adopted the a similar procedure for our thesis.

A Machine disassembling

The claw-pole machine is called a sandwich machine since the stator is maintained between the two brackets. To reduce the number of parts to be reassembled virtually, one of the brackets is glued to the stator (Figure 5.11 at left), then the pulley is removed using a tool as shown in Figure 5.11 at right. Furthermore, since the stator stiness is low, gluing it into a bracket makes it possible (as far as possible) to preserve its shape. Once the pulley and the clamping screws are removed, the machine is disassembled, as shown in Figure 5.12. This leads us to three parts to be measured The bracket on the right denoted front bracket in the rest of this dissertation.

The rotor and the shaft

The bracket and the stator on the left For measuring a workpiece using a CMM, a reference frame needs to be dened on that workpiece. In the following, we detail how these references are dened for each part and how the measurements are performed.

B Stator

In the stator, we are interested in measuring its internal surface to determine its deformation according to its nominal shape that is cylindrical. The stator teeth dene the internal surface. However, before proceeding with this measure, some preliminary ones are done to dene the workpiece's reference frame. The XY -plane is dened by using the plane on the top of the stator. Then the orthogonal vector Z is dened.

The center of the frame is dened using some measurement in the teeth of the stator by a least square method. We use the least square method because the aim is to simulate the assembly, not to verify a conformity in accordance with a standard.

The orientation of XY is dened by pointing the X -axis in the direction of the plastic positioning pin shown in Figure 5.13.

Measurement process

The measurement of the inner surface of the stator is carried out on the 36 stator teeth (Figure 5.14). On each tooth, 15 points are measured, which leads to 540 points that characterize the stator's inner surface.

C Rotor

For the rotor, we want to measure its external surface to dene the second surface dening the air-gap (the rst one is dened by the previous measurement of the stator). The surface of the claws denes the desired surface. The reference frame is dened in the following manner

Reference frame

The Z axis is dened as the axis of the shaft and the ball bearing

The XY -plane is dened by the surface (shoulder) retaining the axial movement of the rotor (as shown in Figure 5.15 at right).

The orientation of XY is dened by pointing the X -axis in a direction between two claws.

Measurement process

The measurement of the outer surface of the rotor is carried out on the 12 claws (Figure 5.16). As the claws from each side of the rotor are not identical, the number of measurement points is not the same; indeed, twelve points are performed on each claw of the upper ones and eight on the lower ones.

D Virtual assembly

We now have dealt with the measurements of the inner surface of the stator and the outer surface of the rotor. In order to be able to reconstruct the air-gap, we have to locate the stator versus the rotor. The air-gap is theoretically a zone of two coaxial cylinders, can be represented in a plane perpendicular to the axis as a curve in polar coordinates, which depends on the angular position of the rotor with respect to the stator.

First, we locate the front bracket versus the stator, then the frames of these two parts are confounded. Afterwards, the rotor is located in the frame dened on the stator.

Between stator and brackets

The bracket on which the stator is glued is positioned on the other bracket using the interface between the two parts. The center of the pin and the bore is then used for locating (Figure 5.17).

The assembling is done virtually by confounding the reference frame dened on the bracket with the one dened on the stator.

Summary

In this Section, we presented how the geometry of the claw-pole machine is measured. Because of a lack of accessibility to the air-gap, we disassembled the electrical machines and proceeded to measure the part individually. We showed how the reference frames are dened in each part. Reference frames are essential for the virtual assembly of the disassembled machine. It allows us to assemble the machine using measurement and calculation; this enables us to characterize the air-gap size and shape, although its small size.

Raw data and variability modelling

In this Section, we show the measurement data of each part of the claw-pole machine.

We propose a methodology to construct a stochastic model of the geometry from the measured points to exploit it in the context of uncertainty quantication.

Stochastic modelling aims to reduce the number of degrees of freedom; for example, on the internal surface of the stator, we measure 540 points. These points represent the degrees of freedom that describe the internal surface; however, all these cannot be taken into account for practical reasons (very chaotic surface) and computation (the curse of dimension). Thus, we deduce some variables representing the maximum information about the surface while being of small size.

Stator A measurementsmentss

As noted before, we measure the internal surface of the stator; In Figure 5.19, we show the measured data of one of the stators. The colours highlight the deviation from the nominal radius of the stator to the measured one.

R measured = R nominal + D (5.1)
where D is the deviation.

One can see the following 1. The deviation is not homogeneous on the surface.

2. The deviation varies between -50 µm and 20 µm.

3. The mean measured radius is smaller than the nominal one.

4. The surface tends to have an elliptic shape.

These observations are relatively the same for all the measured machines but with orders of magnitude that are slightly dierent. The purpose here is to propose a modelling technique for the shape of the stator's internal surface. As we have seen in (5.1), the measured radius of the stator is the dierence between the nominal radius and the deviation D, then, to model the internal surface of the stator we need to nd an appropriate model for the deviation

D.

From the observations mentioned above, we can deduce a simple model for the deviation. The rst part of this model is dedicated to the variation of the radius and the second part dedicated to the elliptic shape. This model is written as follows

D ≈ a 0 + a 2 cos(2θ) (5.2)
where a 0 model the variability in the radius, θ is the angular coordinate, a 2 is the amplitude of the elliptic deformation.

The parameters a 0 and a 2 can be calculated using the least square method.

The resulted shape is shown in Figure 5.20. One can visually see that this model reasonably approximates the deformation.

We adopt this approach since it reduces drastically the number of degrees of freedom to be taken into account in the numerical evaluation. These two variables (a 0 , a 2 ) model the shape of the stator's internal surface.

The same procedure is applied to all the machines to determine these variables for each one. Since we have made the measurement on 30 machines, we get at the end 30 values of a 0 and a 2 . Table 5 5.1: Variables modeling the stator shape for some machines Figure 5.21-a and 5.21-b show the histograms of the values of parameter a 0 and a 2 respectively while 5.21-c shows the distribution of variable a 2 with respect to variable a 0 . The variables a 0 and a 2 seem to be independent.

Rotor A measurementsments

Rotor measurements are somewhat tricky to represent graphically since they are performed on the twelve claws, six from the front of the rotor and six of its rear. To quantify the impact of each type on the global variability, we carry out an analysis of variance (ANOVA) study. ANOVA enables us to characterize the most signicant sources of variances and, most importantly, to eliminate the sources that do not contribute much in the total variability. The presentation of ANOVA is out of the scope of this dissertation; we refer the reader to the literature for further information [START_REF] Kaufmann | Analysis of variance anova[END_REF].

The factors of the ANOVA are : it means that the variation on a claw along the axis direction is signicant) 2. Claw (if this factor is important, it means that the variation of dimension along the angle position is signicant) 3. Rotor (if the factor rotor is signicant, it means that the mean of rotor radius changes from a machine to another)

The results of the analysis are summarized as follows:

Around 33% of variability is related to the factor "rotor" Around 27% of variability is related to the factor "claw" Less than 3% of variability is related to the factor "measurement point"

The remaining is related to the interactions between factors.

From these results, we notice that the variability related to the factor "measurement point" does not contribute much in the global variability; thus, it can be eliminated without a signicant impact on the analysis. It means that we do not need to develop a model accounting for position along the claw, i.e. along the z-direction.

The factor "measurement point" is responsible for the variability of measurement points on each claw independently. By eliminating this factor, we can consider that there is no variability on each claw. Using this argument, we can say that the radius of the rotor is constant for each claw but varies when going from a claw to another and from a rotor to another. From the ANOVA, we can see that the factor "rotor" and "claw" have an inuence, meaning that the mean radius will change signicantly from one machine to another and that the angular position also has an inuence.

In other words, based on the visual observation in Figure 5.22, One can notice that the measured cylindrical surface presents some eccentricity; in fact, the axis computed using a least square method is shifted from the axis that denes the nominal cylindrical surface (see the denition of the reference frame of the rotor).

Furthermore, the measured surface radius is smaller than the nominal one. Simply, we want to nd the deviation, which is the dierence between the measured radius and the nominal radius. The rotation axis is dened by the coordinate of the extremities X f and X r .

Each extremity is dened by three coordinate (x,y,z). Thus, we have six degrees of freedom for dening this axis. We can eliminate the z coordinate if we consider the x and y coordinates at the limit of the rotor at X n f and X n r .

B Variable modelling

We consider the x and y coordinates of the points addressed before as the variables to be controlled; they are dened as follows :

x and y coordinates at X n f are denoted c x and c y respectively

x and y coordinates at X n r are denoted d x and d y respectively

In Figure 5.25 and 5.26, we show the distribution of these variables.

Variability modelling

The measured data enables us to model the parameters that represent the variability of the electrical machine's shape. These parameters are calculated for each machine, and then we need to identify the probability distribution of each parameter. We consider the parameters as random variables with normal distribution. We have identied the best parameters (mean and standard deviation) that are the best tted to represent the sample of each parameter. We have also carried some Normality tests in order to validate our assumption. We used the toolbox developed by M.

Öner et al. to conduct ten of the most used Normality tests. As a result, at least ve out of the ten tests have shown that the normality hypothesis cannot be rejected (for a p-value of 0.05) for all the parameters. Figure 5.27 shows all the variable histograms with their relative normal distributions.

Summary

In this Section, we have shown the protocol developed for measuring the electrical machines; the machines were rst disassembled into three parts. The measurements were exploited to get intuitions of the parameters that model the variability. Each measurement point represents a degree of freedom of the geometry of the machine. Treating all the degrees of freedom is impractical and intractable due to their high number (540 in the stator, 120 in the rotor, ...). Therefore, we proposed some key parameters (a 0 , a 2 , b 0 ,...) to capture the maximum information about the shape of the measured surfaces. These parameters were mainly deduced using graphical representation and some statistical methods such as ANOVA.

Once the variables were dened for all the claw-pole machines and the variability were modelled as random variables by determining their marginal probability density function (these variables are considered to be independent), the probability distributions are considered as input for the FEA model that is used for the variability (uncertainty) propagation.

RAW DATA AND VARIABILITY MODELLING

FEA model

The modelling of the claw pole machine using nite element analysis is cumbersome, primarily because of the tridimensional geometry. The diculty related to the creation of the model is dealt with gradually to be able to run the uncertainty propagation. The robustness of the model is a key-ingredient to the success of this task. We start by creating a parameterized model of the electrical machine with its nominal geometry, material property, electrical circuit and functioning conditions (temperature, excitation current, ...). Then, the variables introduced in the previous Section are also added to extend the capabilities of the model.

For this analysis, we use JMAG as a nite element tool for the machine's electromagnetic modelling. Thanks to its numerous advantages in terms of geometry modelling, mesh morphing, eccentricity handling, etc

Geometry

The geometry of the electrical machine is built based on the 2D technical drawings (as shown in Figure 5.28) describing its dierent parts. The idea is to have a geometry the closest to the nominal one. The parameterization is not a straightforward task due to the complexity of the geometry and the invalid shapes generated when varying the parameters. Thus, the robustness of the model is critical when parameterizing a geometry. We can notice that the simulation results are in good agreement with the mean of the experimental measurements with a relative error of less than 1.3 %. It worth noting that the variability on the measured machines is around 4A for all the functioning points.

Variability parameters

In the following we aim at quantifying the eect of the random input parameter dened in the previous Section on the current. We proceed to take into account the variability that was identied using the metrology. We present how each of the parameters is taken into account using the developed model of the claw pole machine.

Uncertainty propagation

The uncertainty propagation aims to propagate the variability that we have characterized on the shape of the electrical machine and then compute the variability on the output current. For this purpose, we use a meta-model approach. The metamodel is constructed based on a small sample from the computer expensive FEA then the meta-model is used for statistical inferencing.

Sampling

The sampling aims at evaluating the model of the electrical machine at some sample points that represent the distribution of the variables that have been characterized.

As we considered the variables to be normally distributed, we use an LHS method adapted to normal distributions [START_REF] Stein | Large sample properties of simulations using latin hypercube sampling[END_REF], also implemented in MATLAB as lhsnormal.

We sample 30 points and ran the simulations for all the points. As each simulation take around 10 hours to complete, all the simulations were performed in around 13 days on a Windows core i7, 4 cores of 2.71Ghz.

Meta-modelling

After the evaluation of the FEA, we proceed to t a meta-model on the computed data; the Kriging model is used for this purpose. This variability needs to be reduced; therefore, identifying the variables that led to this variability is of high interest. In the following, we present a sensitivity analysis for this task.

Sensitivity analysis

In order to determine which input parameters inuence the most the variability of an output parameter, a sensitivity analysis is carried out. It involves estimating sensitivity indices that quantify the inuence of one input or group of inputs on the output. In our case, we want to study the impact of the variables that represent the deformation of the geometry on the output current of the machine. There exist many methods for sensitivity analysis [START_REF] Clénet | Approximation methods to solve stochastic problems in computational electromagnetics[END_REF]. One of the most known is the Sobol indices, which based on the decomposition of variance. If we want to reduce the variability on the output current, it will have almost no impact if we reduce the variability of other parameters than a 0 and b 0 . The variable a 0 models the deviation of the stator radius from the nominal one, and the variable b 0 models the deviation of the rotor radius from the nominal one. By reducing the variability on these two variables, the variability of the output current will be reduced. The other variables seem not to impact the output current, but they may have a signicant impact on other performances of the machine, such as acoustic noise.

Chapter Summary

In this Chapter, we presented a complete procedure of conducting a reliability analysis of the electrical machine; we started with the presentation of the electrical machine. Then, we proceeded by presenting dierent metrological measurements that were conducted to characterize the variability on the geometry. From these measurements, we have derived a simple stochastic model representing the deviation of the machine's geometry according to the nominal geometry. This model is dened by several random parameters, and their probability density functions have been identied. Afterwards, we presented a parametric nite element model of the machine and how the measurement can be taken into account in the numerical model.

The numerical evaluation of the impact of geometry variables on the performances of the machine was conducted using a meta-model approach to reduce the computational cost. The variability of the performance was then characterized, and the variables that inuence this variability were also identied.

Optimization under uncertainty of the electrical machine was not considered in this dissertation. Nevertheless, the nite element model is developed for this purpose. Other variables (shape, functioning temperatures, diode characteristics, ...) can be considered; thus, the optimization can be performed if further computational power is available.

5.6. CHAPTER SUMMARY

Concluding remarks

Conclusion

This dissertation dealt with the approaches for optimization and reliability analysis of electrical machines modelled by nite element method. The nite element method is the most sophisticated tool to model electromagnetic phenomena in electrical machines. It enables the ne modelling of the electromagnetic elds in the studied domain and handling complex geometries, such as the claw-pole machine. However, it is computationally expensive because of nonlinear materials, 3D geometries, time dependency. Thus, its use for optimization and reliability analysis (iterative processes) should be made with caution since only a limited number of evaluations of the simulation tool can be tolerated.

General optimization algorithms cannot be applied as they are; some necessary modications are needed to make them usable for optimizing electrical machines. In the rst chapter, we presented the literature about optimization algorithms, and we discussed their adaptability to this dissertation's subject. From the categories, two approaches were considered, which depend on the capabilities of the model used.

1. A non-intrusive approach based on the Kriging meta-model

2. An intrusive approach based on the adjoint variable method

The non-intrusive approach using surrogate models is widely used in the context of expensive model optimization. Many researchers had highlighted the challenges of using such a strategy. In this dissertation, we have emphasized some of the very well-known issues and proposed a new methodology of using the meta-models and how to speed up the optimization times. We have highlighted the drawbacks of using the conventional approach, which consists in tting a single meta-model on the whole design space and enriching it sequentially using inll criteria. We proposed a novel strategy that consists in building many meta-models on specied regions of the design space; each one of them is relatively easy to t and easy to use for optimization (less modality). Then, iteratively, prune the regions that are not promising. This process explores all the design space; thus, it enables us to produce reliable solutions. where

K ij = T ∈T h g 1 µ ∇α g i J -1 T • ∇α g j J -1 T |J T |w g b j =
T ∈T h g J s3 α g j |J T |w g A Derivatives w.r.t mesh nodes

The derivatives ∂ x k g and ∂ y k g are computed as follows

∂ x k g = ∂ x k Ku -∂ x k b (C.6) ∂ y k g = ∂ y k Ku -∂ y k b (C.7)
such as

∂ x k K ij = T ∈T h g 1 µ ∇α g i ∂ x k (J -1 T ) • ∇α g j J -1 T |J T |w g + 1 µ ∇α g i J -1 T • ∇α g j ∂ x k (J -1 T ) |J T |w g + 1 µ ∇α g i J -1 T • ∇α g j J -1 T ∂ x k (|J T |)w g ∂ x k b j = T ∈T h g J s3 α g j ∂ x k (|J T |)w g ∂ y k K ij = T ∈T h g 1 µ ∇α g i ∂ y k (J -1 T ) • ∇α g j J -1 T |J T |w g + 1 µ ∇α g i J -1 T • ∇α g j ∂ y k (J -1 T ) |J T |w g + 1 µ ∇α g i J -1 T • ∇α g j J -1 T ∂ y k (|J T |)w g ∂ y k b j = T ∈T h g J s3 α g j ∂ y k (|J T |)w g
These formula may seem intimidating, however, there are similar to the ones in the FE code. Only the derivatives of the jacobian J T are needed to be computed and then put into the formulas.

B Derivatives w.r.t state variables

The derivative ∂ u g is computed as follows

∂ u k B i = 0 1 -1 0 j ∇α i j J -1 T i δ jk (C.11)
such as δ jk is the Kronecker delta.

C Calculation of gradient

The objective function is written as follows

W (p) =

Figure 1 :

 1 Figure 1: Position of Claw-pole alternator and detailed view.

10 1 10 2 10 3 10 4 10 5 10 6 10 Figure 2 :

 102 Figure 2: Computational time versus the number of elements

Figure 1 . 1 :

 11 Figure 1.1: General process of optimization.

Figure 1 .

 1 Figure 1.1 shows the general framework of an optimization process; an optimization algorithm predicts where to evaluate the model to solve the optimization problem. It is an iterative process of model evaluation; The algorithm repeatedly evaluates the model at a variable and gets the values of the objective function and

  vergence. The main dierences among derivative-based algorithms are : How to look for the direction d k ? What is the step length α k to take in direction d k ? Algorithm 1: Derivative-based algorithms Choose p 0 and initialize k = 0; while Not converged do Find a descent direction d k ; Determine a step length α k that reduces the function f suciently; Compute the new iterate p k+1 = p k + α k d k ; Set k = k + 1;

  else replace the worst point p nv+1 with the reected point p r ; end else if f (p r ) ≥ f (p nv ) then Compute contracted point p c = p o -γ(p o -p nv+1 ); if the contracted point is better than the worst point, f (p c ) < f (p nv+1 ) then replace the worst point p nv+1 with the contracted point p e ; else replace all points except the best p 1 with p i = p 1 + δ(p i -p 1 ); end else replace the worst point p nv+1 with the reected point p r ; end end In n dimensions, the NM method maintains a set of n + 1 test points arranged as a simplex. It then extrapolates the behaviour of the objective function measured at each test point in order to nd a new point and to replace one of the old test points with the new one, and so the technique progresses. α, β, γ and δ are parameters of the algorithm, and they have a signicant eect on the performances of the algorithm. Standard values are α = 1, β = 2, γ = 1 2 and δ = 12.

Algorithm 3 :

 3 Genetic algorithm create a random initial population ; while Not terminated do evaluate the objective function at the current population ; rank the population based on the objective function values ; select the parents based on the ranking ; choose an elite from the current population to pass the next one ; produce children from the parents ; replace the parents with the children to form the next population; end Extension to constrained problems

(1. 43 )

 43 Note that equation 1.43 corresponds to the local form of Ohm law. Equation 1.42 transcribe the magnetic behaviour of materials.The physical model of the electromagnetic phenomenon is a combination of Maxwell's equation(1.38 -1.39) and constitutive laws(1.42 -1.43)

Figure 1 . 2 :

 12 Figure 1.2: Example of 2D square triangulation, (left) coarse triangulation, (right) ne triangulation.

Figure 1 .

 1 Figure 1.3: The reference element (left), an arbitrary element (right).

Figure 1 . 5 .Figure 1 . 4 :

 1514 Figure 1.4: Conventional structure of a nite element code

Figure 1 . 5 :

 15 Figure 1.5: Geometry of the inductor [65]

Figure 1 . 6 :

 16 Figure 1.6: Modelled geometry of the inductor (left), Studied domain with boundaries (right)

Figure 1 . 7 :

 17 Figure 1.7: Identiers of geometry

Figure 1 . 8 :

 18 Figure 1.8: Coarse mesh (left), Fine mesh (right)

Figure 1 .

 1 Figure 1.9 shows the distribution of the magnetic ux density B in the studied domain for both the meshes shown in Figure 1.8.

Figure 1 . 9 :

 19 Figure 1.9: Magnetic ux density in the studied domain. Coarse mesh (left), Fine mesh (right)

Figure 1 . 10 :

 110 Figure 1.10: Example of extreme re-meshing errors.

Figure 2 . 1 .Figure 2 . 1 :

 2121 Figure 2.1: Flowchart of Kriging meta-model based optimization.

(2. 1 )

 1 Figure 2.2: Test function

Figure 2 . 4 :Figure 2 . 5 :

 2425 Figure 2.4: Latin hypercube DoE

Figure 2 . 6 :

 26 Figure 2.6: The evolution of the number of research papers dealing with Kriging [source: dimensions.ai]

Figure 2 . 7 :

 27 Figure 2.7: Exponential correlation function

( 2 .

 2 10) : β = (0.2251 -0.7847)

Figure 2 . 8 :

 28 Figure 2.8: Kriging surrogate model

(2. 23 ) 2 . 1 . 1 CLLFFigure 2 . 9 :

 2321129 Figure 2.9: Evolution of CLLF with respect to θ 1

Figure 2 . 10 :

 210 Figure 2.10: Impact of the CLLF on the Kriging model

Figure 2 . 1 CLLFFigure 2 . 11 : 2 0Figure 2 . 12 :

 212112212 Figure 2.11 shows the evolution of CLLF with respect to θ in cong. 3, and the location of the optimal value. Furthermore, Figure2.12 shows the Kriging metamodel that corresponds to this conguration.

  EI was rst introduced to tackle unconstrained optimization. Its main advantage lies on sampling points that oer a good trade-o between local exploitation and global exploration of the design space. Indeed, It tends to choose the design points most likely to improve the accuracy of the model and/or have a better function value than the current best point.

  the one-dimensional test function shown in Figure2.2 is used. EI determines where the algorithm will evaluate the black-box function. It tends to choose the design points most likely to improve the accuracy of the Kriging model and/or have a better value than the current.In Figure2.13, The blue dashed line is the true objective function we wish to minimize, while the solid red line is the Kriging approximation conditional to the sample points shown as green dots. The black plot at the bottom is the sampling criterion (EI); its y-axis is on the right.After the initial sample of four points generated by LHS, the resulting Kriging model is of poor t compared to the true function. Nevertheless, the EI function leads the algorithm to sample points where it is expected to improve the actual smallest value of the objective function while taking into account the prediction error. After three iterations, the model has improved in the region of the optimum on the right. In iterations 5 and 6, EI nds other points in regions where there is a high probability of nding a better point. By iteration seven, all the region on the right has been explored. However, the left part is not accurate compared to the true function because the Kriging model expects no-improvement from sampling in that region. Nevertheless, the global minimum has been found quite accurately as shown in Figure2.13.

Figure 2 . 13 :

 213 Figure 2.13: Iterations of EGO

  max p [EI(p), P F (p)]

Figure 2 . 14 :

 214 Figure 2.14: Constrained test problem

( 2

 2 Sasena et al. have discussed that transforming the ISC to a constrained problem, as shown in (2.40), highly improves the accuracy of the solution found [74]. The ISC for constrained problems is then written as min p EI(p) s.t. ĝi (p) ≤ 0, i = 1 . . . n c

Figure 2 . 15 :

 215 Figure 2.15: Iterations of constrained EGO

Figure 2 .Figure 2 . 16 :

 2216 Figure 2.16 shows the contour plot of the function and the global solution is as asterisk.

Figure 2 .

 2 [START_REF] Lagarias | Convergence properties of the neldermead simplex method in low dimensions[END_REF] represents the cost function and the contraints. The objective function is shown as contour lines; the region inside the black curves is the feasible region while outside is not. This problem has three minima, one in each feasible region. A red asterisk shows the global one.

Figure 2 .

 2 Figure 2.17: Constrained function

Figure 2 .

 2 [START_REF] Jones | Lipschitzian optimization without the lipschitz constant[END_REF], at the left, we see the evolution of the maximum value of the expected improvement (max EI) found by the algorithm, and on the right, the convergence rate (CR) of the algorithm solving the ISC at each iteration. The erratic behaviour in both plots is related to the change in the correlation parameters of the meta-model. But, interestingly, after the tenth iteration, the maximum value of the expected improvement decreases and, at the same time, the rate of convergence.

Figure 2 . 18 :

 218 Figure 2.18: Evolution of the maximum expected improvement and the convergence rate of the algorithm

Figure 2 .Figure 2 . 19 :

 2219 Figure 2.19: Evolution of the condition number as function of iteration

Figure 2 .

 2 20 shows the evolution of tting time as a function of the sample size. The gure reports the mean values of multiple runs.

Figure 2 . 20 :

 220 Figure 2.20: Cost of tting the Kriging model for Example 1

Figure 2 . 21 :

 221 Figure 2.21: Flowchart of B&B Kriging

4 Figure 2 . 22 :

 4222 Figure 2.22: Space exploration by the algorithm

Figure 3 . 1 .

 31 It consists of the model of an innite solenoid powered by a coil of current density J = 0.01A/mm 2 . The position and the width of the coil are respectively R = 0.7m and d = 0.3m. A Dirichlet boundary condition is imposed on the edge shown in black.

Figure 3 . 1 :

 31 Figure 3.1: Innite solenoid model

Figure 3 . 2 :

 32 Figure 3.2: The mesh of the device (Top), and the distribution of magnetic ux density in T (Bottom)

Figure 3 .Figure 3 . 3 :

 333 Figure 3.3 shows the relative error of the derivative computed with the nitedierence compared to the exact derivative while varying ε from 10 -18 to 10 -2 . We can see clearly that the derivatives depend on the value of ε. Small values lead to more errors because of rounding errors. Moreover, high values can lead to some truncation errors, as seen for ∂W ∂J and ∂W ∂d .

Figure 3 . 5 :

 35 Figure 3.5: Same mesh topology, dierent geometries

10 - 18 Figure 3 . 6 :

 1836 Figure 3.6: Relative error of derivative computed by nite-dierence and meshmorphing

Figure 3 . 7 :

 37 Figure 3.7: Electrical circuit

Figure 3 . 9 :

 39 Figure 3.9: (left) The reference element, (right) an arbitrary element.

  (x, y) = (x 0 , y 0 ) + (x r , y r )J T (3.45) with

Figure 3 . 10 :

 310 Figure 3.10: Edges and faces numbering

Figure 3 . 11 :

 311 Figure 3.11: Flowchart of a modular adjoint implementation

Figure 4 . 1 :

 41 Figure 4.1: SMES Device [1]

Figure 4 . 2 ;

 42 two for the coils and one limiting the studied domain ([0, 15] × [0, 15]) which holds the two coils. The variables considered in the optimization problem are the ones related to the coils (position, size, current density imposed). The variables

Figure 4 . 3 ,

 43 Figure 4.3, the distribution of the ux density around the coils is shown.

Figure 4 . 2 :

 42 Figure 4.2: Modelled geometry of the SMES

Figure 4 . 3 ,

 43 the maximal values of |B| are located on the boundaries of the coils and specically on coordinatesP 1 = (R 1 -d 1 /2, 0), P 2 = (R 1 + d 1 /2, 0) and P 3 = (R 2 -d 2 /2,0) (shown by cross sign in Figure 4.3 on the right). Thus, we replace this constraint by three ones (4.5),(4.6) and (4.7).

Figure 4 . 3 :

 43 Figure 4.3: Mesh of studied domain (left), enlarged view of the ux density distribution around the coils (right)

Figure 4

 4 Figure 4.4.

Figure 4 . 4 :

 44 Figure 4.4: Model of die press with electromagnet, Whole view (left), Enlarged view (right) [2]

Figure 4 . 5 :

 45 Figure 4.5: Modelled geometry of the die press

W (p) = 10 i=1( 2 ( 4 . 10 )

 102410 The objective function W is the squared error between the Bx and By values sampled in 10 positions along the arc e-f shown in Figure 4.4. min p B xip -B xio ) 2 + (B yip -B yio ) where p = (R1, L2, L3, L4) are the design variables (their bounds are shown in Table 4.2) and B xio = 0.35 cos( π 40 i), B yio = 0.35 sin( π 40 i) and B xip and B yip are computed by the FEA of the device.

Figure 4 . 6 :

 46 Figure 4.6: The mesh of studied domain (left), The ux density distribution in the studied domain (right)

Figure 5 . 1 :

 51 Figure 5.1: Position of claw-pole alternator

Figure 5 .

 5 Figure 5.2 shows the structure and the main parts of a claw alternator. As with all rotating electrical machines, there are a rotor and a stator enclosed in two brackets.

Figure 5 . 2 :

 52 Figure 5.2: Detailed view of the claw-pole alternator

Figure 5 . 3 :

 53 Figure 5.3: Stator of claw-pole machine

Figure 5 . 4 :

 54 Figure 5.4: Rotor of claw-pole machine

Figure 5 . 5 :

 55 Figure 5.5: The brackets of claw-pole machine

Figure 5 . 6 )Figure 5 . 6 :Figure 5 . 7 :

 565657 Figure 5.6: Electrical diagram of claw-pole machine

Figure 5 .

 5 8 shows the bounds of variation (minimum and maximum) of the measured output current of 30 alternators. The gray area represents the possible values that the output current takes for dierent rotation speeds.

Figure 5 . 8 :

 58 Figure 5.8: Variation of measured output current on 30 machines

Figure 5 . 9 :

 59 Figure 5.9: The Coordinate Measuring Machine

Figure 5 . 10 :

 510 Figure 5.10: The measurement procedure using the CMM

Figure 5 .

 5 Figure 5.11: Claw-pole machine disassembling

Figure 5 .Figure 5 . 13 :

 5513 Figure 5.12: Claw-pole machine disassembled

Figure 5 . 14 :

 514 Figure 5.14: Measurement points of the inner surface of the stator

Figure 5 .

 5 Figure 5.15 shows the reference frame adopted for the rotor measurement.

Figure 5 . 15 :

 515 Figure 5.15: Reference frame of the rotor

Figure 5 . 16 :

 516 Figure 5.16: Measurement points of the outer surface of the rotor

5. 3 . 2 Figure 5 . 19 :

 32519 Figure 5.19: Deviation to nominal of the internal surface of the stator in mm

  .1 show the values of the variable for some machines and Figure 5.21 shows the statistical distribution of the values of the parameters.

2 Figure 5 . 20 :

 2520 Figure 5.20: Model of deviation to nominal of the internal surface of the stator

Figure 5 . 2 (c) Correlation plot between a 0 and a 2 Figure 5 . 21 :

 522521 Figure 5.22 shows the deviation of the radius of the external surface of the rotor to the nominal one. In Figure 5.22-a, we highlight the measure points (triangles) by colour with the scale shown on the colour bar. In contrast, in Figure 5.22-b, we show the deviations and compared to the nominal radius and the manufacturing tolerance allowed on the external surface of the rotor.

1 .

 1 Figure 5.22: Deviation from nominal of the external surface of the rotor

RFigure 5 . 23 :Figure 5 . 24 :

 523524 Figure 5.23: Synthesis of the measurements of the outer surface of the rotor

Figure 5 . 25 :

 525 Figure 5.25: Synthesis of the measurements of of virtual assembling (Part 1)

1 .- 5 - 5 • 10 Figure 5 . 26 :Figure 5 . 27 :

 15510526527 Figure 5.26: Synthesis of the measurements of of virtual assembling (Part 2)

Figure 5 . 28 :

 528 Figure 5.28: Technical drawing of the electrical machines

Figure 5 . 29 :Figure 5 . 30 :Figure 5 . 31 :

 529530531 Figure 5.29: Some parameters of the geometry

Figure 5 . 33 :

 533 Figure 5.33: Eccentricity module of JMAG software

5. 5 . 3 Figure 5 . 34 :

 53534 Figure 5.34: Variability of the output current resulted from variability on geometry

Figure 5 . 35 :

 535 Figure 5.35: Sobol indices of the variables

C. 2 . 1

 21 Figure C.1: Molds of the Die press (design region)

10 i=1( 2 (C. 12 )∂ x k W = 10 i=1 2 (∂ y k W = 10 i=1 2 (= 10 i=1 2 (

 10212102102102 B xip (p) -B xio ) 2 + (B yip (p) -B yio )Its gradient is written asd p W = ∂ p W + λ∂ p g = N k=1 [∂ x k W + λ∂ x k g] d p x k + N k=1 [∂ y k W + λ∂ y k g] d p y k (C.13)where∂ x k B xip (p))(B xip (p) -B xio ) + 2(∂ x k B yip (p))(B yip (p) -B yio ) ∂ y k B xip (p))(B xip (p) -B xio ) + 2(∂ y k B yip (p))(B yip (p) -B yio ) λ ∂ u B xip (p))(B xip (p) -B xio ) + 2(∂ u B yip (p))(B yip (p) -B yio ) (∂ u g) -1

  Linear optimization deals with optimization problems whose data, i.e. objective function and constraints, are linear functions. On the opposite side, nonlinear optimization copes mainly with optimization problems where the objective function or the constraints or both are nonlinear functions. Nonlinear problems are a complicated matter in the eld of optimization because it is challenging to converge to the best solution. Quadratic optimization is a particular case of nonlinear optimization where the objective function is quadratic, and constraints are linear functions.

	Quadratic optimization is widely used to solve nonlinear problems by solving a se-
	quence of quadratic problems.

On the contrary, in continuous optimization problems, variables can take any real value. Continuous optimization problems are usually easier to solve. An optimization problem mixing continuous variables and discrete variables is called a mixed-integer optimization problem.

Continuous optimization problems are more accessible to solve than discrete ones; the values of the objective function and constraint functions at a point x can be used to deduce information about points in a neighbourhood of x. However, advancements in computing technology have dramatically increased the size and complexity of discrete optimization problems that can be solved eciently. Continuous optimization algorithms are necessary for discrete optimization because many discrete optimization algorithms generate a sequence of continuous subproblems.

D Linear vs nonlinear problems E Explicit vs black-box models

Black-box optimization models are models for which the functions dening the objective and constraints are complex, dicult to evaluate, or their analytical expressions are unknown. For most industrial problems, this is a sophisticated computer program that is often inaccessible. They are exploited by analyzing the responses to given inputs; they are considered as black-boxes without any hypothesis on the properties of the functions. The classical mathematical methods based on rigorous hypotheses such as continuity, dierentiability, convexity are not straightforwardly applicable in black-box optimization.

  .3.

		1
	Variable 2	0.2 0.4 0.6 0.8					Samples
		0	0	0.2	0.4	0.6	0.8	1
					Variable 1
	Figure 2.3: Two variables Full-factorial DoE
	Although these designs remain reasonable in low dimensions, they require a large
	number of observations in high dimensions (2 nv , n v is the number of variables),
	making them impractical for computationally expensive problems, i.e. exponential
	number of evaluations is needed.

Table 2 .

 2 It will often be convenient to work with the square root of the M SE, it gives a standard error, or prediction error, for measuring uncertainty in the prediction and referred to as Before tting the Kriging model to the four sampled points, the variables and responses are scaled, as shown in table 2.1. They are scaled in a way to have a mean of zero and a standard deviation of one, other methods of scaling could be considered. When variables do not have the same order of magnitude, scaling becomes very important to treat all the variables in the same manner and to avoid

	√
	σf =	M SE	(2.17)
	To illustrate the Kriging predictor and its MSE, we use the test function men-
	tioned earlier (2.2). The process of tting a Kriging model is described step-by-step
	to foster a better understanding of what is involved in building a Kriging model.

Since we have no uncertainty about the point we have already sampled (because we have calculated it using the black-box), the formula for (2.16) has the intuitive property that it is zero at any sampled point. The MSE oers a metric for assessing the meta-model accuracy and can be used to improve it. Indeed, we can sample an additional point where the MSE is high and consequently expecting to improve the global approximation of the function.

After the function was sampled at the point dened in (2.5), we proceed to the construction of the Kriging meta-model by adopting the following simplifying hypotheses:

The regression term is rst order polynomial The correlation function is exponential with parameters θ 1 = 1 and α = 2 1: Variable and responses scaling

Table 2 .

 2 . 2: Impact of adding samples on the CLLF

		Optimal θ CLLF limCLLF
	1	139.26	-0.64	-0.54
	2	11.83	-0.82	-0.52
	3	7.57	-0.48	-0.44

Table 2 .

 2 2 shows the solution of the CLLF and the corresponding values of lim-

  [START_REF] Brisset | Pareto-based branch and bound algorithm for multiobjective optimization of a safety transformer[END_REF] Now, we fully identied the optimization problem formulated in(2.23), one should choose an optimization algorithm to solve it. As the evaluation of CLLF requires inverting the correlation matrix, which becomes computationally expensive for a large number of samples (size of R is equal to the number of samples). A derivativebased algorithm can be used, but the CLLF presents a plateau where the gradient is quasi-null, as shown in Figures 2.9 and 2.11. This will lead to slow convergence of derivative-based algorithms. As an alternative, we use the Nelder-Mead algorithm that requires only one or two evaluations per iteration and does not require derivatives.

The Nelder-Mead algorithm has two drawbacks; handling constraints and nding an adequate initial point. The problem formulated in (2.23) has only bound constraints; these can be easily integrated into the objective function by applying a variable transformation. Nevertheless, the initial point requires more attention

Table 2 . 3

 23 .3.

	n v cost cost 1% error %	êrror %
	2	27	-	1.5	0.06
	3	45	-	2.2	11.5
	4	87	-	8.2	3.0
	5	157	-	3.5	3.6

: Performances of EGO on the multidimensional Rosenbrock problem

Table 2 .

 2

	cost cost 1% error %	êrror %
	19	-	31.6	31.6

.

4

. The algorithm failed severely to solve the problem or even to locate where the global solution is. As shown in Table

2

.

4

, one can see that the best-sampled point is very far from the global solution (error = 31.6%), and even the solution found by the meta-model is also faraway ( êrror = 31.6%). 4: Performances of EGO on the constrained problem 2.17

In fact, the algorithm found a local solution and got stuck with it. It kept sampling points near to each other in one of the feasible regions found. It is worth noting that the inaccuracy of the constraint meta-model caused this behaviour.

  2.2 Branch and Bound assisted by Meta-models (B2M2)To recall, B&B algorithms consist of a systematic enumeration of candidate solutions by splitting the search space into smaller spaces. The algorithm explores each smaller space then estimates the upper and lower bounds of the objective function on each one of them. A space is pruned if it cannot produce a better solution than the best one found so far by the algorithm. B&B were mainly developed for discrete optimization by other extensions to continuous optimization were proposed. Messine et al. developed an interval B&B algorithm (IBBA)[START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF] [START_REF] Fontchastagner | Design of electrical rotating machines by associating deterministic global optimization algorithm with combinatorial analytical and numerical models[END_REF]. IBBA computes the bounds of the objective function using interval arithmetic. Its performances depend on the ecient estimation of the lower and upper bounds on the sub-regions. The lower and upper bounds are critical elements for IBBA, But, in the case of an expensive black-box model, the real bounds are, in most cases, dicult or even impossible to compute. Thus, it can be applied to only a limited number of situations with explicit objective function and constraints.

Nevertheless, Kriging enables us to model the uncertainty on a black-box function; This uncertainty is modelled by a random variable of a normal distribution of mean f (p) and variance σ 2 (p) often called Gaussian process. Intuitively, this means that we are saying that the function f (p) has a typical value of µ(p) and is expected to vary in some range like [ f (p) -3σ(p), f (p) + 3σ(p)].

  .t. ĝi (p w,i ) -kσ g i (p w,i ) ≤ 0, i = 1, . . . , m.

	(Section 1.1.1). The problems in (2.44) and (2.45) are transformed
	to
	f = min p	f (p w,0 ) -kσ f (p w,0 )
		(2.46)
	f = min

sp f (p w ) + kσ f (p w ) s.t. ĝi (p w ) + kσ p i (p w ) ≤ 0, i = 1, . . . , m.

(2.47)

Table 3 .

 3 .1. We notice that we have fewer errors for energy than for B c . Generally speaking, FEA leads to fewer errors for global quantities such as energy and higher errors for local quantities.

	Analytic	FEA	Rel. error %
	B c 3.757 mT	3.592 mT	4.39
	W 4.5239 J/m 4.5225 J/m	0.03

1: Comparison of analytic and FEA quantities

So now, we are interested in computing the gradient of the quantities of interest with respect to the variables R, d and J. The next section shows some approaches that are often considered for conducting this task.

  .2.

		Analytic	Adjoint	Rel. error %
	∂ R B c	0	0.006	0.60
	∂ d B c	0.0126	0.0125	0.38
	∂ J B c	0.376	0.359	4.39
	∂ R W	5.655	5.655	0.00
	∂ d W	32.044	32.032	0.04
	∂ J W	904.78	904.5	0.03
	Table 3.2: Comparison of analytic and and adjoint variable gradients

Table 3 .

 3 .3. It shows the gradient of the problem's objective function in the optimum. For more details, the reader may refer to reference [1]. 3: Gradient Comparison for TEAM workshop 22

	Method	ε	R 1	R 2	h 1 /2	h 2 /2	p	d 1	d 2	J 1	J 2	Rel. Error	Time (s)
	Adjoint	-	24.827	-19.222	13.725	-10.619	25.852	-81.824	8.418e-07	8.444e-07	-	1.27
	CD w. MM	10 -5 24.803 -19.212 13.726 -10.616 25.855 -81.817 8.418e-07 8.444e-07	0.1%	4.06
	CD w/ MM	10 -3 36.497 -23.592 15.756 -19.098 35.757 -93.534 8.418e-07 8.444e-07	155%	7.02

Table 4 .

 4 E ref = 180M J, B norm = 200µT and p are the design variables and their bounds are shown inTable 4.1. p R 1 R 2 h 1 /2 h 2 /2 d 1 d 1 J 1 J 2 1: Bounds of design variables for TEAM problem 22

						2 norm + |E(p) -E ref |/E ref	(4.1)
	s.t. |J| + 6.4|B| -54 ≤ 0					(4.2)
	R 1 -R 2 +	1 2	(d 1 + d 2 ) < 0				(4.3)
	min	1.0	1.8	0.1	0.1	0.1	0.1	10	-30
	max	4.0	5.0	1.8	1.8	0.8	0.8	30	-10

where

  Then, the interval of variation of A 2 is [0.001, 3.9]. this new variable is somewhat relaxing the initial optimization problem by allowing the variable R 2 = R 1 + A 2 + (d 1 + d 2 )/2 to vary in a more signicant interval[1.101, 8.7] than the one dened in the initial problem[1.8, 5].

	becomes
	min p	OF (p) = B 2 stray (p)/B 2 norm + |E(p) -E ref |/E ref	(4.4)
	s.t.	(J 1 -54)/6.4 + |B(p, P 1 )| ≤ 0	(4.5)
		(J 1 -54)/6.4 + |B(p, P 2 )| ≤ 0	(4.6)
		(-J 2 -54)/6.4 + |B(p, P 3 )| ≤ 0	(4.7)
	Introducing There-
	fore, two constraints (4.8) and (4.9) are added. Finally, the optimization problem

Table 4 .

 4 2: Bounds of design variables for TEAM problem 25

Table 4 .

 4

	3: TEAM Workshop problem 22 optimization results
	Approach	SQP	GA	DIRECT	B2M2
	R 1	1.336	1.457	1.543	1.369
	A 2	0.027	0.481	0.229	0.054
	h 1 /2	1.011	1.209	0.951	0.888
	h 2 /2	1.452	1.800	1.526	1.394
	d 1	0.677	0.347	0.374	0.791
	d 2	0.269	0.121	0.217	0.203
	J 1	15.579	19.834	22.346	14.099
	J 2	-15.069	-17.305	-13.441	-18.273
	OF	0.00197 0.03502 0.04881 0.00510
	# FEM evals	1902	160201	421995	106308
	CR	1 %	2 %	100 %	100 %
	Expected # FEM evals	529640	94276389	421995	106308

TABLE 4

 4

	.4.

Table 5 .

 5 2: Comparison of simulation result with respect to experimental one

1.1. DETERMINISTIC OPTIMIZATION

1.1. DETERMINISTIC OPTIMIZATION

1.2. UNCERTAINTY

1.2. UNCERTAINTY

1.3. ELECTRICAL MACHINES

1.3. ELECTRICAL MACHINES

1.3. ELECTRICAL MACHINES CHAPTER 1. LITERATURE REVIEW

2.1. META-MODEL DESIGN OPTIMIZATION MDO

2.1. META-MODEL DESIGN OPTIMIZATION MDO CHAPTER 2. NON-INTRUSIVE APPROACH

2.1. META-MODEL DESIGN OPTIMIZATION MDO

2.2. BRANCH AND BOUND ASSISTED BY META-MODELS (B2M2)

2.2. BRANCH AND BOUND ASSISTED BY META-MODELS (B2M2)

2.3. CHAPTER SUMMARY

2.3. CHAPTER SUMMARY

3.1. HOW TO COMPUTE THE GRADIENT?

3.1. HOW TO COMPUTE THE GRADIENT?

3.2. EXPLANATORY EXAMPLE

3.2. EXPLANATORY EXAMPLE

3.3. ADJOINT VARIABLE METHOD FOR FEA

3.3. ADJOINT VARIABLE METHOD FOR FEA

3.4. EXAMPLES

4.1. TEST CASES

4.1. TEST CASES

4.2. ALGORITHMS SETTINGS

4.4. CHAPTER SUMMARY

5.1. ELECTRICAL MACHINE

5.1. ELECTRICAL MACHINE

5.2. METROLOGY

5.2. METROLOGY

5.3. RAW DATA AND VARIABILITY MODELLING

5.4. FEA MODEL

B.2. VARIATIONAL FORMULATION

C.2. TEAM WORKSHOP PROBLEM 25

Acknowledgements

First of all, I would like to thank my thesis advisors, Pr. Stépahne Brisset, Pr.

Stépahne Clenet, Associate Pr. Frédéric Guyomarch, and Dr. Jean Claude Mipo, for their availability, great mastery of the dierent subjects treated in this thesis, and ability to formulate complex ideas in simple words.

I would like to thank the jury members for making the defense day possible and pleasant despite the Corona pandemic. Thank you Pr. Frédéric Messine and Pr.

Derivatives for design variables

The derivatives of the quantities f and g with respect to the design variables could be somewhat challenging due to the variety of parameters that could be considered.

Often, in the literature, researchers separate between geometric variables and physical variables where the rst are often referred to as shape sensitivity analysis [START_REF] Park | Design sensitivity analysis for nonlinear magnetostatic problems using nite element method[END_REF] and the latter are mainly used in the context of topology optimization [START_REF] Okamoto | 3-d topology optimization of single-pole-type head by using design sensitivity analysis[END_REF]. In this dissertation, we will treat both types of variables.

A Geometric variables

Geometric variables are the ones that involve changing the geometry of the device studied. A geometric denition of the problem must be made before starting optimization and uncertainty quantication processes. The choice of variables is of paramount importance since it is equivalent to dene the mathematical model of the problem studied. Clearly, it denes the nature, and the dimensions of the research space and possible solutions largely depend on it. The computation of the derivatives highly depends on the way the geometry is parameterized.

Parameterization

In the most straightforward conguration, the shape parameters can be considered as merely a collection of points in 3D space. The analysis, in this case, involves moving each of the points in the desired direction by some small amount and determining the eect on the function. This process is both complicated and computationally expensive. Shape parameterization exists as an eort to overcome these complexities and ineciencies. Shape parameterization is the method of determining a set of parameters that control the size and shape of the device to be designed. This is conducted by determining some driving parameters and reducing the number of design variables from the number of 3D space points to the number of parameters.

There exist many methods for geometry parameterization, such as free-form deformation [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF], polynomial and spline (NURBS) [START_REF] Piegl | The NURBS book[END_REF], iso-geometric analysis [START_REF] Bontinck | Recent advances of isogeometric analysis in computational electromagnetics[END_REF], and computer-aided design (CAD). Among these categories, CAD appears to have some capability that makes it very attractive for design engineers such as eciency, compactness and suitability for complex congurations. Besides, CAD-based parameterization methods can achieve signicant geometry changes. Another advantage of CAD-based is the availability of a comprehensive set of geometric functionalities provided by commercial CAD systems. Parameterizing a complex model is still a challenging task with today's CAD systems. Furthermore, they are not capable of calculating derivatives analytically. Computer codes for commercial CAD systems are huge; to dierentiate an entire system with automatic dierentiation tools may be a non-trivial task. Therefore, the calcu-Chapter 4

Applications and benchmarking

This chapter is dedicated to numerical tests of the developed approaches in the context of optimization. For this task, we adress two well-known benchmarks treated by electromagnetic community researchers. They are the ones of the TEAM (Testing Electromagnetic Analysis Methods) workshops [1] [2]. Since the comparison of results to other researchers is somewhat tricky because of the dierences related to the mesh and the numerical methods that can be used such as the formulation, the integration technique, ..., we solve these benchmarks with the approaches we have developed and some of the algorithms that are conventionally found in the literature.

The test cases are the TEAM workshop problems 22 and 25, there exist two publications describing both of them [1] [2]. We solve both problems using four methods; the non-intrusive approach developed in Chapter 2, an SQP algorithm assisted by the adjoint variable method for computing the gradient, the DIRECT algorithm [START_REF] Finkel | Direct optimization algorithm user guide[END_REF] and the genetic algorithm (GA of [START_REF]The MathWorks, Matlab, Global Optimization Toolbox[END_REF]).

In the following of this chapter, we detail the test cases and their related optimization problems, afterwards, we present the settings of each optimization algorithm for both test cases, and lastly, a comparison of the results is conducted.

Test cases 4.1.1 TEAM Workshop Problem 22

The Superconducting Magnetic Energy Storage (SMES) device in Figure 4.1 consists in two concentric superconducting coils fed by currents that ow in opposite directions [1]. The inner coil is used for storing magnetic energy E, while the outer one has the role of diminishing the magnetic stray eld B stray computed on line a and line b.

Pin Bore

Between rotor and brackets

To locate the position of the rotor with respect to the stator, we measure contact surfaces of the bearings and the shaft. By these measurements, we characterize the centre of the cylindrical surfaces. Then we confound the centers to obtain a virtual assembling of the rotor and the stator.

Assembling

When all the virtual assembling is done, we get the coordinate of the rotor axis in the stator reference frame. Figure 5.18 shows these axes, the z-direction of the rotor reference frame is confounded with the rotation axis and dened by the coordinates at the extremities of the red segment.

Stator axis Rotation axis The parameters a 0 and a 2 are used to model the interior surface of the stator;

thus, these variables used to deform the mesh of the teeth of the stator. On each tooth, we apply a radial deformation (as shown in red for one tooth in Figure 5.32) of amplitude dened by

(5.5)

with i is the number of the tooth and θ i = 10i JMAG is equipped with an eccentricity module (shown in Figure 5.33) that enables us to take into account these variables at once. The parameterization of the eccentricity on JMAG is somewhat dierent from the variables we use. Nevertheless, the transformation that links both of them is not very dicult (geometrical transformations such as translations and rotations).

Summary

In this Section, we presented the developed nite element model of the electrical machine with all the aspects of the geometry, the material properties and the electrical circuit. We veried the model validity with respect to experimental measurements.

We have shown how the parameters determined from the analysis of the metrology were taken into account in the numerical model. Due to the complexity of the deformations involved, such as eccentricity, the symmetries of the machine cannot 5.4. FEA MODEL On the other hand, for the intrusive approach, we propose how to compute the derivatives from a nite element model. Most of the work has been done on the model rather than on the algorithm (we used the SQP algorithm for this dissertation). We have presented the adjoint variable method and how the gradient can be derived from a FEM code. We have developed an ecient way to compute the derivatives of the shape sensitivities for geometric parameters, which are vital for the gradient's computation. The adjoint variable method was compared to nite-dierence to validate and highlight its eectiveness in terms of precision and computational time.

For the comparison between intrusive and non-intrusive approaches for optimizing electromagnetic devices using FEM, we treated two well-known benchmarks from the literature known as TEAM Workshop problems 22 and 25 [1] [2]. We used two metrics for the comparison, the rst one is the quality of the solution, and the second is the computational cost.

In terms of performances, the intrusive approach outperforms other strategies for both test cases; this was possible due to the computation of the gradient using the adjoint variable method. This improves the quality of the solutions drastically but comes with the expense of intrusive manipulation of the FEM code. The non-intrusive approach remains a good alternative in terms of implementation and solution quality. The developed approach was able to overcome some of the very known issues when using meta-models for optimization.

In the last Chapter, we presented a complete procedure for conducting a reliability analysis of the electrical machine; we started with the presentation of electrical machines. Then, we proceeded by presenting the dierent metrological measurements that were conducted to characterize the variability on its shape. Afterwards, we presented a parametric nite element model of the machine and how the measurement can be taken into account in the numerical model. The numerical evaluation of the impact of geometry variables on the performances of the machine was conducted using a meta-model approach to reduce the computational cost. The variability of the performance was then characterized, and the variables that inuence this variability were also identied.

In summary, the goals of this dissertation have been met. A comparison of the approaches appropriate for the optimization of electrical machines has been carried out. The advantages and shortcomings of each method have been highlighted.

The integration of reliability analysis in the design phase has been assessed on a real test case (a claw-pole machine). The impact of the manufacturing process on the geometry of the machine has been studied, and quantitative analysis has been accomplished.

Perspectives

Further work is currently underway to enhance the abilities of the non-intrusive approach. The multi-objective strategy is further enhanced for reducing the computational burden. The Pareto front can overlap many sub-spaces, which increases the computational cost of the method. A new technique for space division is being tested to tackle this issue.

The calculation of the gradient with respect to geometrical variables using the adjoint variable method relies mainly on the shape sensitivities. The geometric parameterization of shape variables is still one of the shortcomings of the method.

We have presented an approach based on the parametric equation of the geometric shapes; however, for very complex shapes, this can be very cumbersome. Some researchers use dedicated tools that couples the CAD software with the optimizer, such as The Computational Analysis PRogramming Interface (CAPRI) [START_REF] Brock | Adjoint-based design optimization using cad parameterization through capri[END_REF] [START_REF] Haimes | Capri cae gateway[END_REF].

CAPRI serves the purpose of providing custom communications from a computational software suite to the preferred CAD system. This allows designer access to the CAD system's geometry denitions and functionalities, providing the ability to query the CAD system whenever needed. Using this tool will enable us to compute the shape sensitivities related to variables automatically from the CAD software.

Another use of the adjoint variable method will be dedicated to topology optimization, where the number of variables is signicant. The gradient of the quantities of interest is computed at a minimal cost. Dierent methods might be considered, such as the SIMP method and the level-set method.

For the claw-pole machine, optimization under uncertainty will be considered to exploit the model full capabilities developed during the thesis. Taking into account the uncertainty related to the geometry of the machine while varying its shape can lead to highly ecient, while reliable, machines when considering the actual manufacturing process.

Appendix A Shape sensitivities

In this appendix, we derive the shape sensitivities for some recurrent types of geometries.

A.1 Rectangle

A parametric equation of the rectangle can be written as

Then, the shape sensitivities are calculated with respect to the variable dening the rectangle are

A.3 Other shapes

In a general manner, one can compute the shape sensitivities if the considered shape can be written as parametric equations of the edges dening it. We showed the shape sensitivities for the circular arc and the rectangle since there will be used in the application of the end of this chapter and the following chapter, where we will treat some optimization benchmarks.

Appendix B FEM formulation in cylindrical coordinates

This Appendix deals with the variational formulation for axisymmetric electromagnetic devices. There is few researchers who addressed this aspect. Here we present, the formulation that we used in MagFEM. The formulation is used for the modelling of the SMES device of the TEAM Workshop problem 22.

B.1 General equations

As detailed in Chapter 1, Maxwell equations are partial dierential equations (PDE) (Hard formulation). These equations can be solved analytically for simple systems.

For complex ones, numerical methods are best suited using weak formulations, i.e.

variational formulation.

The magnetostatic problem is written as follows.

∇ × 1 µ

These equations describe the electromagnetic elds and the corresponding boundary conditions.

B.2 Variational Formulation

Multiplying the PDE from B.1 by an arbitrary test function v such that n×v| Γ B = 0 and integrating over the domain leads to

Applying Green-Ostrogradski formulas on the l.h.s This leads to the variational formulation

In axisymmetric case, the potential vector A and the ux density B are reduced to the following thanks to the invariance around the symmetry axis.

Then, the problem is reduced to Except for the term 1 r in a(u, w), this formulation is similar to the 2D. Thus, the discretization is the same and the quantities of interest can be computed from the u, for example, the ux density is calculated as follows. This appendix treats the derivation of the gradient equations using the adjoint variable method for the TEAM Workshop problem 25 [2].

C.1 Gradient equations

In general, Finite element model solve a set of equations g(u, p) = K(u, p)u-b(p) = 0 for u parameterized by variables p (design parameters) and compute a quantity of interest f (u, p)

To recall from Section 3.2.2, for computing the gradient of a function f , the following steps has to be conducted :

1. Compute the partial derivative for u :∂ u f and ∂ u g The full description of the problem is detailed in Section 4.1.2.

B Variables L2 and L3

The variables L2 and L3 parameterize an ellipse arc (ij in Figure C.1). The equation of an ellipse is written as

Then, the arc ij can be parameterized by the following equations.

where y m = 10.5mm is y-coordinate of the vertex j in Figure C.1.

C Variable L4

The variable L4 parameterize the segment km in Then, the derivatives on the segment km are deduced:

The quantity ∂ u g is often computed in the nite element code. the subroutine could be reused in the adjoint variable code.

C.2.3 Quantity of interest

The aim is to minimize the function W (p) are ux densities computed by the FEA of the device, thus, we need to compute the derivatives of these quantities w.r.t. the mesh node coordinates (x i , y i ) and the state variables u.

The ux density at position i is computed as follows

where ∇α i j is the gradient of the base functions evaluated in the reference element and J T i is the jacobian of transformation from the reference element to the element T i where the the point i is located. Now, the derivatives of these quantities w.r. Optimisation et analyse de abilité des machines électrique modélisées par la méthode des éléments nis La méthode des éléments nis est l'outil le plus sophistiqué pour modéliser le phénomène électromagnétique. Cependant, elle est coûteuse en temps de calcul. Ainsi, son utilisation pour l'optimisation et l'analyse de abilité (processus itératifs) doit être faite avec prudence car seul un nombre limité d'évaluations du modèle peut être toléré. De plus, l'impact du processus de fabrication sur les machines électriques est peu étudié dans la littérature. L'intégration de cet aspect dans la phase de conception est l'un des apports de cette thèse aux côtés de la principale contribution, qui est le développement et la comparaison des approches d'optimisation pour les machines électriques. Nous exposons les approches adaptées au sujet de la thèse et en développons de nouvelles. D'une part, le modèle d'éléments nis peut être considéré comme une "boîte noire" pour laquelle nous développons une approche non intrusive basée sur des méta-modèles. D'autre part, nous considérons une approche intrusive, nous améliorons le modèle pour fournir les dérivées des quantités d'intérêt. Les dérivées sont calculées ecacement en utilisant la méthode de la variable adjointe. Finalement, les méthodes sont comparées pour donner un aperçu des avantages et des inconvénients de chacune d'entre elles. Enn, une étude de cas réel est abordée ; elle consiste à étudier l'impact du procédé de fabrication sur la machine à gries fabriquée par Valeo. Sur la chaîne de production, les machines sont prélevées pour mesurer leurs dimensions et caractériser leur écart par rapport aux dimensions nominales. Ensuite, une analyse statistique est menée pour évaluer la abilité et l'impact sur les performances.

Mots-clefs : Machines électriques, Meta-modèles, Méthode de la variable adjointe, Propagation d'incertitudes, Méthode des éléments nis.

Optimization and Reliability Analysis for Electrical Machines modeled by Finite Element Method

This dissertation deals with the approaches for optimization and reliability analysis of electrical machines modelled by the nite element method. The nite element method is the most sophisticated tool to model the electromagnetic phenomenon. However, it is computationally expensive. Thus, its usage for optimization and reliability analysis (iterative processes) should be made with caution since only a limited number of evaluations of the model can be tolerated. Furthermore, the impact of the manufacturing process on the electrical machines is scarcely studied in the literature. The integration of this aspect in the design phase is one of the contributions of this thesis alongside the main contribution, which is the development and comparison of optimization approaches for electrical machines.

We present the approaches adapted to the subject and develop new ones. On the one hand, the nite element model can be seen as a "black-box" for which we develop a non-intrusive approach based on Kriging meta-models. On the other hand, we consider an intrusive approach as we look inside the "black-box," we upgrade the model to provide the derivatives of the quantities of interest. The derivatives are essential to some optimization and reliability analysis tools. They are computed eciently using the adjoint variable method. Finally, the methods are compared to give insight into the advantages and the shortcomings of each of them.

Lastly, a real case study is considered; it consists of studying the impact of the manufacturing process on a claw-pole machine manufactured by Valeo. From the production line, machines are withdrawn to measure their dimensions and characterize their deviation from the nominal one. Then a statistical analysis is conducted to assess the reliability and impact on the performances.