Hazaël Jones

Nicolas Sabouret

Atef Ben Youssef

Strategic Intentions based on an Affective Model and a simple Theory of Mind

Keywords: Theory of Mind, Strategic intentions, Affective model, Job interview

Cet article présente un modèle informatique pour raisonner sur les émotions de l'interlocuteur, en utilisant un paradigme de théorie de l'esprit (Theory of Mind, ou ToM, en anglais). Le système manipule des représentations sur des croyances à propos des émotions, des préférences et des buts de l'interlocuteur. Notre modèle affectif est conçu dans le contexte de la simulation d'entretien d'embauche mais il n'est pas lié à un ensemble d'affects spécifique. Il s'appuie sur des règles simples pour sélectionner les types de question lors de l'entretien en fonction de la personnalité de l'agent. Nous l'avons implémenté en utilisant une représentation de type OCC des émotions et un modèle dimensionnel PAD pour les humeurs.

Introduction

In order to build a credible interaction between a human and a virtual character, affective computing [Picard, 1995] proposes to simulate human affects in virtual agents, making them more realistic and engaging for interactions. In this context, one main challenge for Artificial Intelligence researchers is to make the virtual character adapt its behaviour to the perceived user's affective state, which will lead to a more natural and credible interaction for the user.

To this purpose, we claim that virtual characters must not only use reactive behaviour in answer to a wide range of affects (emotions, moods, social attitudes...) such as in [START_REF] Marsella | PsychSim: Agent-based modeling of social interactions and influence[END_REF], Kriegel and Aylett, 2008, Jones and Sabouret, 2013, Schröder, 2010]. It must also use strategic intentions about the human it interacts with. Strategic intentions can be seen as long term goals [Haddadi, 1996] for an agent. Indeed, in an interaction, people have intentions about the goal of a conversation, such as obtaining an information, finding an agreement, changing the interlocutor's point of view or having a fun and relaxing conversation. This paper proposes to analyse these strategic intentions and to use them in the reasoning model of an affective agent. To this purpose, we define a general model that can be adapted to different context. In our work, we apply this general model to a specific case of a formal interaction: job interviews in which the goal of the recruiter is to obtain concrete information about the interlocutor's social and technical skills, so as to select the best candidate. Our general model is based on logical rules and is inspired by the theory of mind [Baron-Cohen, 1995] paradigm. Based on affect perception from Social Signal Interpretation (SSI), our virtual agent's model derives beliefs about the interlocutor's self-estimation in the interview (job skills, importance of the salary, etc). These informations are confronted to the agent's goals so as to select the next course of actions in the interaction (in our case, to conduct the job interview). This paper is organized as follows. Section 2 makes a brief state of the art on theory of mind and shows how this has been used in the context of a virtual agent's reasonner. Section 3 briefly presents the job interview context and its specific features. Section 4 presents our architecture. In Section 5, we present in details our affective model that integrates the theory of mind for reasoning about affects in the context of interactions. The rules of this general model are illustrated on examples from the job interview context. The last section concludes on the model and its application to the job interview situation.

Related work on Theory of Mind

Theory of mind [Leslie, 1994, Baron-Cohen, 1995], or ToM, is the ability to attribute mental states (beliefs, intentions, desires, affects, . . . ) to others. The literature reports numerous ToM studies and implementations in agent-agent interaction [Bosse et al., 2007, Dastani andLorini, 2012]. In our work, we want to model the reasoning process of an agent that reason about the reasoning process of a human (the applicant). This particular configuration raises additional difficulties and leads to a original model for our representation of the ToM. For example, in [START_REF] Pynadath | PsychSim : Modeling Theory of Mind with Decision-Theoretic Agents[END_REF], an agent has beliefs about others in a subjective way. Agent A has belief about agent B following the real structure of agent B beliefs. However, in our work, since agent B is the human applicant, we do not have any information about its belief structure. We must guess them from the outputs of the affect recognition module.

Nevertheless, this model of influence and belief change [START_REF] Pynadath | PsychSim : Modeling Theory of Mind with Decision-Theoretic Agents[END_REF] is based on work in psychology: the authors use influence psychological factors in their simulation framework: consistency, selfinterest, speaker's self-interest and trust (or affinity). We believe that similar high-level reasoning structures must be proposed in reasoning models, to complement lowlevel reasoning on perceptions such as what is done in [Scassellati, 2002, Peters, 2006]. These papers focus on the perception aspects of ToM such as the desire of engagement, and are tailored for signal interpretation, not for the cognitive model of the virtual agent.

Several other applications have been studied with a Theory of Mind approach. For instance, [START_REF] Bosse | A two-level BDI-agent model for theory of mind and its use in social manipulation[END_REF] proposes a reasoner for task avoidance, the agent can change its behaviour in order to alter the other agent's desires, intentions and in fine, actions to occur. This work has been extended in a more generic version [START_REF] Bosse | A recursive BDI agent model for theory of mind and its applications[END_REF] that proposes a two-level BDI agent model: the first level is the agent's reasoner and the second one computes the ToM. Following a different approach, [START_REF] Dastani | A logic of emotions: from appraisal to coping[END_REF]] also propose a model based on modal logics that extend the BDI paradigm. Each agent has a set of actions and a set of formulas that represent the agent's mental state. A formula has a degree of desirability for the agent and a degree of plausibility. The use of modal logic allows researchers to model the recruiter beliefs, desires and intentions, but it seems difficult to represent a real humans' mental states based only on perceptions. This is the reason why we propose a model based on general rules that takes as inputs recognised affects from the interlocutor and strategic intentions for the virtual agents, and combines them in the ToM-based affective model. The goal of our model is to represent the reasoning process of an agent that reason about the reasoning process of a human. Our model will be applied and illustrated in the context of the TARDIS project 1 that considers a job interview simulation as an interaction.

1 TARDIS stands for Training young Adult's Regulation of emotions and Development of social Interaction Skills. url: http://www. tardis-project.eu/

Job interview context

In job interviews, the recruiter needs to reason about the potential behaviour of the applicant in front of him. This evaluation is done by selecting different questions in order to provoke particular reactions on the interviewee [START_REF] Rynes | Applicant reactions to alternative selection procedures[END_REF]. For example, to test the applicant's capability to manage his or her stress, the recruiter can be voluntarily aggressive during the interview. Our goal in the TARDIS project is to model that kind of recruiter strategic intentions. To this purpose, we propose in the next section a formalism to represent these strategies and to reason about the applicant's state of mind, based on perceived social and emotional signals.

It is important to note that job interview simulation is an interesting situation for studying multimodal affective interaction with a virtual agent. The literature shows that expressed emotions and non-verbal behaviour play a key role in job interview: 1) it is used by the recruiter for evaluating the candidate [Sieverding, 2009] and 2) it influences the applicant's behaviour [Sieverding, 2009]. In addition, [START_REF] Gatewood | Human resource selection[END_REF] shows that the recruiter uses different questions to assess the applicant's work performance, individual quality and specific regarding the job. In our work, we take into account these three aspects. First, the applicant's non-verbal performance is evaluated by comparing expected social cues to detected ones, using the SSI system. Second, we compute an affective reaction (see section 5.1) for the recruiter (that should influence the candidate's behaviour). Third, we select the next topic of interest for the recruiter that tries to evaluate the applicant's competencies (see section 5.2). We propose several strategies for a recruiter (in section 5.6), from provocative to helpful, which change this topic selection (and its affective reaction) so as to influence the applicant's behaviour.

TARDIS general architecture

Figure 1 shows our global architecture. The TARDIS architecture considers four main components:

• The SSI component provides the affective model with information about the applicant's affects and social attitudes that are detected by the system.

• The Interview Scenario component tells the virtual recruiter the expectation in terms of emotions and attitudes, depending on the interview progress. In TARDIS, the agent has no understanding of the applicant's actual answers to the questions. It follows a scenario, that can be influenced by the recruiter perception and internal states and focus on the affective recognition and adaptation.

• The Animation component is responsible for expressing the virtual recruiter's affective state through its behaviour and expressions.

• The virtual recruiter component which is composed of two modules: -The Decision module, that is the focus of this document. The goal of this module is to build a theory of mind for a cognitive agent in the context of a job interview. Our agent (the recruiter) will deduce intentions of the applicant considering its answers (based on SSI) in a particular context (the question that has just been asked by the recruiter). This model will also influence new questions.

A ToM-based model for a cognitive virtual agent

In this section, we will present our model, and then shows its application in the TARDIS project.

General model for theory of mind

Our main objective is to draw beliefs about the interlocutor's mental states, preferences and understanding of the situation in the course of human-agent interaction, based on the user's reaction in terms of non-verbal behaviour (and social signal interpretation). In human-agent interaction, the agent has to select a new question at each turn-taking.

To choose each new question, it is interesting that the agent has an idea of interlocutor's mental state regarding the past questions. This can be used in a wide spectrum of domains in which a human interacts with a virtual avatar in an interview simulation, such as teaching, training, . . . The common aspect of these simulations is the use of questions by the avatar. Our model considers the use of question in order to manage the context of the answers of the person in interaction with the simulation.

To summarize, our theory of mind model has three main properties:

• It is about a real person who interacts with the system,

• It is centred on the person's preferences, expectations and interest for the job,

• It uses the context of questions and the affective behaviour to analyse user responses.

Our ToM however does not include any memory, other that the immediate response of the person and the previous ToM. We do not represent the knowledge that the interlocutor might have acquired during the interaction: we focus on its (estimated) preferences and expectations.

Context Management

With a view to manage the context, labels are given to the questions/sentences of the virtual character in order to interpret the answer/reaction of the human in term of beliefs on some topics. A list of topics can be done for each specific application. The set of topics set topic contains N topics: {topic 1 , topic 2 , . . . , topic N }. Each subject is applica-tion dependent and based on the domain of the simulation. A question is concerned by 0 to n topics.

List of topics. In order to manage the context, some labels are given to the questions of the recruiter in order to interpret the answer of the applicant in term of beliefs on certain subject. Here is the set of topics set topic that can be tagged for a job interview:

• topic applicant : questions about the applicant (general questions),

• topic job : questions about the job,

• topic salary : questions about the salary,

• topic hours : questions about the working hours,

• topic skill : questions about the competencies, the skills of the applicant regarding the job,

• topic socialSkill : questions about the general social skills of the applicant, A question is concerned by 0 to n topics. For example, the question "In what position will you like to work in our enterprise?" can be tagged by two different topics: topic skill and topic job because it tells about the applicant skills (the position he thinks he can apply for in this job) and its knowledge of the job (organisation of the enterprise).

Beliefs build

In order to build beliefs about the human who interacts with the system, we consider the questions/sentences that were just expressed by the virtual agent (identified by labels about topics) and the quality of the answer of the human from an affective point of view (which is obtained by Social Signal Interpretation, or SSI). Based on that, the agent will update its beliefs about the human on a particular subject. We denote the beliefs of the agent about the human B Human (topic i ) for i in {1, . . . , N }.

According to the topic(s) raised by the question/remark of the agent, beliefs will be updated. In pursuance of building the beliefs of the human, we consider its answer (perceived via SSI) and decide if the answer is rather positive, negative or neutral. In order to determine if the global answer is positive or not, we use a performance index that compares the expected social cues (such as smile, large gestures, body movement, directed gaze...) with the detected ones. Expectations can be expressed as positive (signals that should be detected) or negative (behaviours that the user should avoid). This method in presented in [START_REF] Jones | Interpreting social cues to generate credible affective reactions of virtual job interviewers[END_REF].

The performance index pi is valued in the interval [0, 1].

Based on this value and the topic tags of the question/remarks just done by the agent, the beliefs can be computed. Updates of each belief are done with the following formula for each topic:

B Human (topic i ) ← B Human (topic i ) + α × pi
with α ∈ [0, 1] a value that can be altered if we want the recruiter beliefs about the human to evolve quickly (α = 1) or not (α near of 0). It can rely on the personality of the agent. An impulsive agent has an α near of 1 and a moderate one near of 0.

For instance, after a question about the job, with α = 1 (impulsive recruiter) and an actual belief value B Human (job) = 0.5, B Human (job) will become -0.3 for an AverageAns of -0.8 making the recruiter beliefs about the applicant change sign in one answer. A moderate recruiter (α = 0.2) will obtain a belief B Human (job) = 0.34 which will change the dynamic of our simulation. An impulsive recruiter will cause strong dynamics and a moderate one smooth ones, which is the expected behaviour.

List of beliefs. The list of beliefs we consider in the context of job interviews is the following:

• self-confidence B Human (young),

• knowledge about the job B Human (job),

• importance of the salary for the applicant B Human (salary),

• importance of scheduling and working hours for the applicant B Human (hours),

• qualities of job skills B Human (skills),

• qualities of social skills B Human (socialSkills),

According to the topic raised by the question of the recruiter, beliefs can mean different things. For instance B(young) is the belief of the applicant in himself. For a belief equal to 1, the applicant is very confident, and for -1, he has an important lack of self-esteem. Here the value quantifies the confidence of the applicant. If we look at B(salary), the signification is different, it is about the importance that the applicant put in the salary when applying for this job.

Desires and goals

The desires are used to define the strategic intentions of the agent. We organize our desires in two categories: the highlevel ones and the more specific ones. The high level intentions are directly linked to social attitudes. Attitudes can be initialized with personality and can evolve during the simulation but with a dynamics slower than the emotional one which is quite reactive. For more detail about the computation of social attitudes, refer to [START_REF] Jones | An affective model for a virtual recruiter in a job interview context[END_REF]. The high-level intentions are about the general intentions of the agent for the interaction, the specific ones are about specific beliefs about the human that interest the agent during the interaction. The high-level desires are denoted: D(Attitude). The specific desires are denoted: D(B Human (topic i )) because specific desires in an interaction are about beliefs of the human on a particular topic. For instance D(B John (f ootball) is the desire of the agent to know if John has knowledge in the f ootball topic.

List of desires and goals. For a job-interview simulation, the recruiter will have a limited set of high-level intentions (provocative, pugnacious, friendly and helpful) and only one of them will be triggered in the same time.

The specific intentions are about subjects that the recruiter want to favour during the interview. The level of each subject will be adapted in function of the high-level intentions and will also consider the beliefs about the applicant. Actually, these specific goals for the recruiter are about the knowledge of applicant's beliefs on certain subjects. For instance, a question about the job will be associated to the goal G(B Human (job)) because this question will give more information to the recruiter about the belief B Human (job).

Dynamics of goals

The high-level desires evolve in function of the social attitude of the agent. Social attitudes used can be defined on Leary circumplex [Leary, 1996]. According to the application, some attitudes will be relevant and some not. As shown by Leary, attitudes can be separated in two categories, the positive ones (friendly, cooperative, extroverted, . . . ) and the negative ones (hostile, critical, . . . ).

Based on these two kind of attitudes, we define algorithm 1 in order to update the desires of the agent. This algorithm works as follows: if the agent has a negative attitude, he intends to select topics with a negative answer for the human. On the contrary, if the agent has a positive attitude, its desires are about topics with a positive answer from the human.

Goal selection

Several strategies can be defined for the selection of one desire in the list of possible desires. The most natural one is to select the desire with the maximum value in the available desires. At one moment of the dialogue, every possibilities (topics) cannot be approached in order to conserve the logical sequence of the conversation (a scenario for instance).

Goal selection based on recruiter's high level intentions.

The high level goals can be defined directly in the scenario or be computed on the personality of the recruiter. We define 4 main strategies (2 for the positive attitudes and 2 for the negative attitudes). A recruiter with positive attitudes will have positive desires on topics where the applicant has positive average answers. On the contrary, a recruiter with negative attitude will have positive desires on topics where the applicant has negative average answers. Here are some strategies that we use for the virtual recruiter:

• Provocative recruiter: the recruiter will have a negative attitude and will always select the worst topic for the user (the one with the maximum Desire for the negative agent). Intention = max(B Human (subject))

• Pugnacious recruiter: the recruiter will have a negative attitude but will randomly select one of the worst topic but not always the same. Intention = random(max n (B Human (subject))) with max n , the n worst subjects.

• Friendly recruiter: the recruiter will have a positive attitude and will randomly select one of the best topic for the user but not always the same. Intention = random(max n (B Human (subject))).

• Helpful recruiter: the recruiter will have a positive attitude and will always select the best topic (the one with the maximum Desire for the positive agent). Intention = max(B Human (subject))

These different strategies lead to different goals. At one moment of the interaction, only some subjects can be approached according to the possibilities of the scenario. The recruiter will select the maximum in the available choices.

Conclusion

In this article, we propose a theory of mind model for an affective virtual agent. The theory of mind is about a real person in interaction with the system. It is centred on the interpretation of the affective states perceived through Social Signal Interpretation. By building beliefs about the person in interaction with the simulation, we allow an interaction by understanding the subjects where the person is confident or not. Then, according to the virtual agent high level intentions, new questions will be selected in a coherent and credible way regarding the personality of the agent. This work is actually in the process of integration in the TARDIS platform. After integration, it will be evaluated in order to confirm that our model proposes a credible virtual recruiter for a job interview scenario. The theory of mind should provide coherent actions of the recruiter according to the reactions of the applicant and the personality of the recruiter. This aspect can be evaluated through literature and thanks to the applicants that will interact with the system. One current limit of our model is that it requires manual annotation of the scenario. The definition of an automated annotation process, based on the utterance's semantic and contextual information, would greatly increase the scalability of the model.

Figure 1 :

 1 Figure 1: Global architecture for a recruiter in a job interview -Affective and Decision modules

Algorithm 1

 1 Desires computation if (Attitude ∈ set(attitude -)) then for B Human (topic) ∈ set topic do if (AverageAns < 0) then D(topic) ← D(topic) + α × |AverageAns| else D(topic) ← D(topic) -α × |AverageAns| if (Attitude ∈ set(attitude + )) then for B Human (topic) ∈ set topic do if (AverageAns < 0) then D(topic) ← D(topic) -α × |AverageAns| else D(topic) ← D(topic) + α × |AverageAns|

Acknowledgment

This research is funded by the European Union Information Society and Media Seventh Framework Programme FP7-ICT-2011-7 under grant agreement 288578.


Reda El Bechari
email: reda.el-bechari@centralelille.fr

Stephane Brisset

Stephane Clenet

Jean-Claude Claude Mipo

Stéphane Brisset

Stéphane Clénet

Enhanced Meta-model Based Optimization under Constraints using Parallel Computations

Keywords: Constrained optimization, Expensive simulation, Kriging, Infill criteria, Parallelization strategy

published or not. The documents may come

L'archive ouverte

I. INTRODUCTION

M ETA-MODELS [1] are used in many fields, mainly to replace expensive black-box models [START_REF] Li | Kriging-assisted multi-objective design of permanent magnet motor for position sensorless control[END_REF] [START_REF] Kleijnen | Expected improvement in efficient global optimization through bootstrapped kriging[END_REF]. 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. Optimization using kriging meta-models were first introduced in [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] to tackle unconstrained optimization. Its main advantage is the reduction of the number of calls to the expensive model. However, for problem with high number of parameters the number of evaluations arises exponentially (curse of dimensionality). Thus, the purpose of this paper is to compare methods to handle constraints and propose a new strategy for parallelization, which enables to run several evaluations at each iteration. A brief review of meta-model based optimization and infill criteria for constrained optimization is presented. Then, the parallelization strategy is presented and tested on an analytical model and the TEAM workshop benchmark problem 22.

II. META-MODEL BASED OPTIMIZATION

An optimization problem can be formulated as follows min x y(x) s.t. g i (x) ≤ 0, i = 1 . . . n c [START_REF] Jin | Comparative Studies of Metamodeling Techniques Under Multiple Modeling Criteria[END_REF] where x are the design variables, y(x) is the objective function and g i (x) are n c constraints. In the case of electromagnetic modeling, finite element models are widely used, but the simulations can be expensive in terms of computational time.

An alternative consists in replacing them by cheaper models. These meta-models are constructed with a relatively small size sample obtained by the Finite Element model at first, but refined with other samples found using an infill criterion. Furthermore, the numerical noise due to discretization present a problem for gradient based algorithm. The optimization in this case may become difficult or even impossible.

Meta-model based optimization flowchart is presented in Fig. 1. The first step aims to determine initial set of parameter values (initial design) using a design of experiments, e.g. Latin Hypercube Sampling (LHS). The full (expensive) model is solved for each set of parameters. Afterwards a metamodel is built based on the initial design and the output data. Kriging is well suited for building the meta-model due to its statistical properties. The most important part in the process, on which the paper will focus, is to find the infill point which improves the actual best solution and increases the metamodel precision. This point will be evaluated using the Finite Element model in the next iteration. In Fig. 1, the arrow that goes out and comes back into step "Find Infill point" means that a sub-optimization problem is solved to find infill points. Finally some stopping criteria are evaluated to terminate the optimization.

A. Kriging

Kriging [START_REF] Martin | Use of Kriging Models to Approximate Deterministic Computer Models[END_REF] is an interpolation method based on a regression term and a stochastic term. The stochastic term aims to eliminate the error due to regression and is constructed based on the location of the sampled points. Furthermore Kriging characterizes the variance, or the precision, of the prediction. Thus, Kriging considers the response a normally distributed random variable with given expected value and standard deviation.

The Kriging method, also called Gaussian process regression, was originally developed by D. Krige a mining engineer in South Africa. Afterwards the method was introduced into the field of numerical design.

In [START_REF] Lophaven | DACE A Matlab Kriging Toolbox[END_REF] (Dace a Matlab Kriging Toolbox) an exhaustive presentation and the implementation of kriging predictor are detailed. The developed toolbox is used in the numerical evaluation to construct a response surface of the objective function and of the constraints from a set of sampled points.

III. INFILL CRITERIA A. Expected Improvement

The most used infill criterion to deal with unconstrained optimization is Expected Improvement (EI) criterion, presented in [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF]. This criterion enables a trade-off between exploitation and exploration of the design space.

EI(x) = (y min -ŷ(x))Φ(u(x)) + ŝy (x)φ(u(x)) (2)
where u(x) = ymin-ŷ(x) ŝy(x)

, ŷ(x) and ŝy (x) are the expected value and the standard deviation of the kriging predictor for the objective function y and y min is the smallest sampled value of y.

Maximizing EI(x) leads to the point x * with the highest probability of improvement, either by sampling toward the optimum or improving the approximation of the meta-model. These characteristics can be justified by the fact that the derivative of EI(x) with respect to ŷ(x) is negative, meaning that the smaller ŷ(x) the higher EI(x) (exploitation) and the derivative EI(x) with respect to ŝy (x) is positive, meaning that the bigger ŝy (x) the higher EI(x) (exploration).

There exist other infill criteria, that often reveal striking similarities. An exhaustive set of infill criteria was presented in [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF] for unconstrained and constrained optimization.

B. Constraint handling

The probability of feasibility (P F ) criterion is widely used for constrained optimization. It quantifies the probability that a constraint is satisfied

P F (x) = Φ -ĝ(x) ŝg (x) (3)
where ĝ(x) and ŝg (x) are the expected value and the standard deviation of the kriging predictor for a constraint.

To sample point that improves the actual solution and respects constraints both EI and P F should be considered. Table I summarizes the main formulations that aim to satisfy these requirements. The formulation (A) has a statistical derivation (assuming that the objective function and the constraint are statistically independent). It aims to sample points that maximize the expectation of improvement and fulfill constraints. However the product of EI and P F reveals high modality and algorithms maximizing this criterion often fail to find the global optimum.

Thus, the formulation (B) was proposed in [START_REF] Durantin | Analysis of multi-objective kriging-based methods for constrained global optimization[END_REF] to consider the infill criterion as bi-objective to reduce the modality. From Pareto front, the point chosen is the one that maximizes the product of EI and P F . 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 impacts the search close to the constraints boundary, so if the optimum lies on the constraint boundaries these infill criteria may fail to find it. This is due to the maximization of P F that leads to point inside the feasible region and less likely on the constraint boundaries.

The formulation (C) was proposed in [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF] and considers the problem as a constrained one to reduce the modality of the infill criterion and to gain in precision of the solution in case of the optimum activates some constraints. P tol = 0.95 was recommended but it has an effect on the precision, it fails to locate the points on the constraint boundaries. The authors opinion is that P tol = 0.5 seems more reasonable because P F = 0.5 when ĝ(x) = 0. 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.

IV. PARALLELIZATION STRATEGY

In the sequential process, only one point is found by maximizing EI. The aim of the parallelization strategy is to find multiple promising points at each iteration to take advantage of clusters by running distributed expensive evaluations.

In [START_REF] Ginsbourger | Kriging Is Well-Suited to Parallelize Optimization[END_REF] the straightforward extension of EI from sequential to parallel was presented. It computes an analytic expression for two points. For more than two points, the authors propose its estimation through Monte Carlo simulation.

In [START_REF] Schonlau | Computer experiments and global optimization[END_REF] a hybrid method was proposed. It adds artificially the point found at each iteration to a subset, evaluates it with the kriging predictor, and reconstruct the meta-model. Afterwards, the next point is found until reaching the number of required points. When the number of evaluated points increases, the construction of the meta-model become a time consuming step which penalize the whole process.

The proposed method is based mainly on the multimodal behavior of EI and searches the points that maximize EI and exclude its vicinity to find another point. The process is repeated until the number of excluded points equals the number of required parallel evaluations. The exclusion area is defined by the distance from that point.

In the infill criterion, the exclusion area is defined by the constraint x-xi range(x) 2 ≥ 1 4n 2 , where range(x) is the range of variable x and n is the number of parallel evaluations wanted, this constraint defines an area around the points already found in the current iteration to look for other promising regions in the design space. The advantages of this strategy over the two abovementioned is that we have an analytical form for the infill criterion and the reconstruction of the meta-model is not needed multiple times at the same iteration. However, it presents a small drawback, at the last iterations of the optimization EI becomes 0 in a large part of the design space. The points found, at an iteration, are close to each other. This drawback will be further investigated, an initial idea is to vary the size of the exclusion area, or using another infill criterion (e.g. ŷ(x) -kŝ y (x) [START_REF] Sasena | Exploration of Metamodeling Sampling Criteria for Constrained Global Optimization[END_REF] [11] k = 2.5, 3) which keeps the multimodal behavior.

V. ENHANCEMENTS A. Pre-Optimization

Frequently, in the case of constrained optimization the initial sample (generated by LHS) has no feasible point. Thus, before starting the optimization. An initial infill criterion is used for this purpose, it aims to find feasible points.

max x P F (x) (4)
Maximizing P F means that we look for the point that is most likely feasible.

In the case were multiple feasible points are wanted, another infill criterion is proposed

max x P F (x) × D(x) (5)
where the distance D(x) = min d∈X f eas

x -d and X f eas is the set of feasible sampled points found until the current iteration. This infill criterion enables to sample points that have the highest probability of feasibility while being distant from previously sampled points.

B. Cheap Constraints handling

Cheap constraints are the constraints that can be evaluated rapidly without calling the FE model. Generally these of constraints are explicit functions of design variables. Sometimes these constraints express the feasibility of design, it means that if this constraint is not fulfilled the modeled phenomena have no physical meaning (e.g. overlapping constraint of SMES Device). These constraints are handled differently, they are embedded in the infill criterion directly.

Given g(x) ≤ 0 a cheap constraint, the infill criterion problem is written Adopting this procedure ensures that the point found by the infill criterion satisfy the constraint.

max x EI(x) s.t. P F (x) ≥ P tol (6) g(x) ≤ 0

VI. ANALYTIC TEST

The analytic problem from [START_REF] Kim | A single-loop strategy for efficient reliability-based electromagnetic design optimization[END_REF] is taken for comparison purpose.

min x -(x1+x2-10) 2 30 -(x1-x2+10) 2 120 s.t. 1 - x 2 1 x2 5 ≤ 0 1 -(x1+x2-5) 2 30 -(x1-x2-12) 2 120 ≤ 0 (7) 1 - 80 x 2 1 +8x2+5 ≤ 0 0 ≤ x 1 , x 2 ≤ 10
To compare formulations, 10 different initial designs of experiments are generated by the LHS are used. Furthermore, SQP algorithm with 10 start points uniformly distributed in the design space, is presented.

The averages of the results are shown in Table II. The metrics used for comparison are the convergence rate (C.R.) that is the percentage of the 10 tests that converge to the known solution for less than 1% of the range of variable. dist is the normalized Euclidean distance to the known solution for the tests that converged, evals is the average number of evaluations of the exact (supposed expensive) model, and iters is average number of iterations, this metric is only considered for parallel.

The table shows that SQP and the formulation (C) with P tol = 0.95 has the lowest convergence rate. The formulations (A) and (B) have good convergence rates but the number of evaluations is higher. Due to the multimodality of the problem, the algorithm often fails to find the global optimum. The formulation (C) with P tol = 0.5 shows the best results among the sequential formulations.

For the parallelization strategy, three evaluations of the exact model at each iteration were done. The results show that the number of exact model evaluations has increased by 40%, however the number of iterations has decreased by 52%.

VII. TEAM WORKSHOP PROBLEM 22

The SMES device [START_REF] Alotto | Smes optimization benchmark: Team workshop problem 22[END_REF] consists of two concentric superconducting coils fed with currents that flow in opposite directions. The inner coil is used for storing magnetic energy E, while the outer one has the role of diminishing the magnetic stray field B stray .

The goal of the optimization problem is to find the design configurations (8 parameters) that give a specified value of stored magnetic energy and a minimal magnetic stray field. Mathematically, this is formulated as

min x OF (x) = B 2 stray (x) B 2 norm + |E(x)-E ref | E ref s.t. |J| + 6.4|B| -54 ≤ 0 (8) R 1 -R 2 + 1 2 (d 1 + d 2 ) < 0 where E ref = 180M J, B norm = 200µT and x are the design variables x = (R 1 , R 2 , h 1 /2, h 2 /2, d 1 , d 2 , J 1 , J 2 ).
The FE simulation is considered as a Black-Box model with inputs and outputs and response surfaces using Kriging for the objective function and the constraints are constructed.

Thus, three response surfaces are constructed, one for the objective function, and two for the quench constraint. These last two response surfaces represent each one the quench condition in a coil (coil 1 and coil 2). The third constraint is handled as a cheap constraint because it depends directly on the design variables.

The results of optimization are summarized in Table III. The comparison is done between the reference [START_REF] Alotto | Smes optimization benchmark: Team workshop problem 22[END_REF] and two formulation of infill criterion, sequential formulation C with P tol = 0.5 and the parallel formulation (also with P tol = 0.5) to find n = 8 points at each iteration.

The number of evaluations (evals) has increased by 21% on the other hand the overall time for the optimization has decreased 36%. These performances are due to two reasons. Firstly, the parallelization of evaluations. Secondly, the number of reconstructions of the kriging meta-model (3403 times for sequential, and 516 for parallel).

VIII. CONCLUSION

In this communication, we have developed a strategy of optimization based on kriging meta-model and the parallelization of the computations of the full model. The results obtained on an analytical example were promising.

This strategy was assessed on the TEAM Workshop problem 22 with 8 variables, the results show that the proposed parallelization strategy enables to reduce significantly the computational time (a speedup of 1.56). However, the expected speedup was to be 8 (8 points in parallel), this will be further investigated, to enhance the performances.

Fig. 1 .

 1 Fig. 1. Flowchart of Kriging meta-model based optimization.

Fig. 2 .

 2 Fig. 2. Infill point determination

Fig. 3 .

 3 Fig.3. SMES Device[START_REF] Alotto | Smes optimization benchmark: Team workshop problem 22[END_REF]

TABLE III OPTIMIZATION

 III

	Enegy(M J)	178.73	180.14	180.0
	Bstray(µT )	14.56	8.81	8.93
	OF	0.0123	0.0029	0.0020
	dist	0	2.85e-3	3.41e-3
	evals	-	3403	4128
	iters	-	3403	516
	time(s)	-	10254	6541

RESULTS

Reference C (P tol = 0.5) Parallel


Sajid Hussain

Abdelkader Benabou

Stéphane Clénet

David A Lowther

Temperature Dependence in the Jiles-Atherton Model for Non-Oriented Electrical Steels: An Engineering Approach

Keywords: Electrical machines, electrical steels, ferromagnetic materials, Jiles-Atherton (JA) model, temperature

High operating temperatures modify the magnetic behavior of ferromagnetic cores which may affect the performance of electrical machines. Therefore, a temperature-dependent material model is necessary to model the electrical machine behavior more accurately during the design process. Physics-inspired hysteresis models, such as the Jiles-Atherton (JA) model, seem to be promising candidates to incorporate temperature effects and can be embedded in finite element simulations. In this paper, we have identified the JA model parameters from measurements for a temperature range experienced by non-oriented electrical steels in electrical machines during their operation. Based on the analysis, a parameter reduction has been performed. The proposed approach simplifies the identification procedures by reducing the number of model parameters and does not require any additional material information, such as the Curie temperature. The resulting temperature-dependent JA model is validated against measurements, and the results are in good agreement.

I. INTRODUCTION

M ODERN electrical machines may run for a prolonged time at extremely high temperatures which affect the magnetic behavior of the electrical steels used in their cores. Iron losses are usually computed using the Steinmetz equation [START_REF] Steinmetz | On the law of hysteresis[END_REF]-based formulae in computer-aided design (CAD) simulations, and the effects of the temperature on the magnetic properties of the ferromagnetic material and iron loss are often neglected. These effects should be taken into account during the CAD of the modern electrical machines to optimize their performance.

The magnetic properties (B-H relationship and iron losses) of the ferromagnetic materials are highly affected by the variation of temperature [START_REF] Bertotti | Hysteresis in Magnetism: For Physicists, Materials Scientists, and Engineers[END_REF]. Physics-inspired hysteresis models, such as the Jiles-Atherton (JA) model [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF], are promising candidates to incorporate temperature effects and can be embedded in finite element (FE) simulations. Although the physical aspects of the JA model and its parameters have been widely discussed [START_REF] Zirka | On physical aspects of the Jiles-Atherton hysteresis models[END_REF], this is the most computationally efficient and easy-to-implement hysteresis model available in the literature and can be coupled with the FE simulations to incorporate the effects of hysteresis in field solutions. The JA model has five input parameters (M S , a, α, c, and k) and the analytical expressions for the identification of these parameters are given in [START_REF] Jiles | Numerical determination of hysteresis parameters for the modeling of magnetic properties using the theory of ferromagnetic hysteresis[END_REF]. The details of the JA model and its parameters are given in Section II.

The temperature-dependent JA model was recently presented in the literature [START_REF] Raghunathan | Modeling the temperature dependence of hysteresis based on Jiles-Atherton theory[END_REF], [START_REF] Raghunathan | Theoretical model of temperature dependence of hysteresis based on mean field theory[END_REF]. Two additional parameters were introduced, and analytical expressions were proposed for the evolution of the JA model parameters with temperature. In addition to this, some numerical adjustments were applied, such as the value of parameter "c" was constrained in order to keep it within the physical limit, i.e., 0-1. The model was applied to cobalt ferrite and requires parameters to be determined at 0 °K, which is not possible; therefore, mathematical extrapolations were performed. The accuracy of the proposed model is limited, especially around the Curie temperature. A similar model was proposed in [START_REF] Ladjimi | Modeling of thermal effects on magnetic hysteresis using the Jiles-Atherton model[END_REF] which required parameters to be measured at room temperature instead of 0 °K. The model was applied to the 3F3 material.

A new method was presented in [START_REF] Messal | Temperature dependent extension of the Jiles-Atherton model: Study of the variation of microstructural hysteresis parameters[END_REF] based on the direct identification of parameters from the measured data at different temperatures and was compared with the model presented in [START_REF] Raghunathan | Modeling the temperature dependence of hysteresis based on Jiles-Atherton theory[END_REF] and [START_REF] Raghunathan | Theoretical model of temperature dependence of hysteresis based on mean field theory[END_REF]. The proposed method was applied to NiFe 80/20 material. Although the model predicted good agreement with the measurements, the temperature variations of the parameters could not be explained. Also, the direct relation of the model parameters with the physics of the material could not be established. None of the materials reported in [START_REF] Raghunathan | Modeling the temperature dependence of hysteresis based on Jiles-Atherton theory[END_REF]- [START_REF] Ladjimi | Modeling of thermal effects on magnetic hysteresis using the Jiles-Atherton model[END_REF] is electrical steel, and all have Curie temperatures less than 369 °C. The goal of this paper is to create a temperature-dependent JA approach which can represent the magnetic behavior of non-oriented (NO) electrical steels with suitable accuracy for the temperature range experienced in electrical machines during their operation. The JA model has five parameters (M S , a, α, c, and k) which a priori depend on the temperature [START_REF] Raghunathan | Modeling the temperature dependence of hysteresis based on Jiles-Atherton theory[END_REF]- [START_REF] Ladjimi | Modeling of thermal effects on magnetic hysteresis using the Jiles-Atherton model[END_REF]. The idea is to identify a way to reduce the number of parameters needed to model the dependence on temperature. The approach is applied to 35WW300 NO electrical steel (B H =2500 A/m = 1.68 T and P loss,B=1.5 T = 1.96 W/kg) and the results are compared with measurements in Section IV.

II. JILES-ATHERTON MODEL

The JA model is one of the most popular hysteresis models inspired by the physics of ferromagnetism. It explains the hysteresis mechanism with the help of domain wall motion. The two modes of domain wall transitions (bending and translational motion) contribute to reversible and irreversible components of magnetization, respectively. The total magnetization is computed using a differential equation

d M d H = c (1 + c) d M an (M s , a, α) d H + 1 (1 + c) (M an (M s , a, α) -M) βk μ o -α(M an (M s , a, α) -M) (1) M an (M S , a, α) = M s coth H +αM a - a H +αM (2)
where M an (M S , a, α) is the anhysteretic magnetization which is computed using Langevin's function ( 2), M S is the saturation magnetization, α is the interdomain coupling coefficient, a is a parameter that determines the shape of the anhysteretic curve and has the units of magnetic field, k is the pinning coefficient, and c is the domain wall flexibility coefficient. These are known as the JA model parameters and are identified from the measured major B-H loop. The details of these parameters and the differential equation are given in [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF]. One limitation of the JA model is that the five parameters are identified from the major loop, and thus cannot model the inner loops accurately. However, this issue can be resolved using a combination of physical and empirical approaches to represent the experimental data, as presented in [START_REF] Hussain | The modified Jiles-Atherton model for the accurate prediction of iron losses[END_REF].

III. PROPOSED APPROACH

An engineering approach is presented here to predict the magnetic behavior of the NO electrical steels at different values of operating temperature. The B-H loops (B max = 1.5 T) were measured for 35WW300 NO electrical steel (sample size: 150 mm × 150 mm) at different values of temperature using a high-temperature Brockhaus Single Sheet Tester placed in an electric oven. The electric oven has a maximum operating temperature of 330 °C. Some of these measured B-H loops are shown in Fig. 1 demonstrating the effect of the temperature on the B-H loop's characteristics. The final values of the JA model input parameters are obtained for each of the measured B-H loops using a curve fitting technique, i.e., the nonlinear least squares (NLS) method [START_REF] Kis | Parameter identification of Jiles-Atherton model with nonlinear least-square method[END_REF]. The optimization function, used in this case, is the square of the difference between the measured and computed data points. The NLS method is a local search optimization method, which can lead to the solution faster if the initial guess of the JA parameters is carefully estimated. Jiles et al. [START_REF] Jiles | Numerical determination of hysteresis parameters for the modeling of magnetic properties using the theory of ferromagnetic hysteresis[END_REF] have proposed a method for determining these parameters using the characteristics of the B-H loop, and we have used the same technique in this paper to obtain the initial guess for the NLS method.

This paper targets the application of electrical steels in electrical machines at operating temperatures. Therefore, we have measured the B-H loop at seven temperatures within a limited temperature range (i.e., 50 °C-330 °C) which is well below the Curie temperature of electrical steels (T C > 700 °C). As a result, we obtain a total of 35 values for the five JA model parameters, i.e., five at each temperature. Instead of developing the JA model parameters with analytical formulae, we consider the variation of the JA model parameters from the measured B-H loops and perform a primary analysis to reduce the number of parameters based on their thermal evolution.

A. Primary Analysis-Saturation Magnetization (M S ) and Reversibility Factor (c)

The first parameter in our analysis is M S , which represents the saturation magnetization. It can be seen in Fig. 1 that the saturation magnetization (for constant applied magnetic field intensity, H ) decreases with the increase in temperature. The evolution of the parameter M S identified from measurements at different temperatures for 35WW300 NO electrical steel is shown in Fig. 2(a) and the variation is less than 4%. The temperature dependence of M S in the studied steel sample correlates with the Weiss molecular theory [START_REF] Bertotti | Hysteresis in Magnetism: For Physicists, Materials Scientists, and Engineers[END_REF] which states that if the operating temperature is well below the Curie temperature, T C (i.e., T /T C < 0.5), as in our case, the variation of M S with temperature could be neglected [START_REF] Bertotti | Hysteresis in Magnetism: For Physicists, Materials Scientists, and Engineers[END_REF]. Therefore, in this paper, we have assumed that the value of M S to be constant, i.e., 1.2678 × 10 6 A/m, which is the value, obtained from the B-H loop measured at 50 °C.

The reversibility factor, c, in the JA model, represents the flexibility of the domain walls [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF] and depends on the surface energy of the domain walls, pinning field, and the applied magnetic field intensity, H . The values of c, obtained using the NLS method, are shown in Fig. 2(b) and seem to slightly increase at relatively high temperatures. However, it has been observed that a small variation in c does not have a significant effect on the area of the B-H loop. Hence, the parameter c was also assumed to be temperature independent, i.e., 0.02.

It is important to mention here that M S and c are identified using the B-H loop measured at the lowest temperature of interest, i.e., 50 °C. The rest of the JA model parameters (i.e., a, α, and k) are then identified from the measured B-H loops at other temperatures. The reduction in the number of unknown parameters makes the identification procedure simpler and faster.

B. Domain Density (a) and Domain Coupling Coefficient (α)

The domain density a is related to Boltzmann's energy and has a direct dependence on temperature, i.e., a = k B T /μ o m, where k B is Boltzmann's coefficient, T is the temperature, μ o is the permeability of free space, and m is the magnetic moment of domains per unit volume. The domain coupling factor α in the JA model is responsible for the interdomain coupling which is due to the exchange interaction at the domain level [START_REF] Jiles | Theory of ferromagnetic hysteresis[END_REF].

Two of the three remaining JA model parameters (i.e., a and α), identified from the measurements, are plotted against temperature in Fig. 2(c) and(d), respectively. It can be seen that both parameters show a linear dependence on the temperature for the given temperature range and can be modeled using the first-order polynomials (i.e., with two unknowns only).

C. Pinning Factor (k)

The coercivity of a ferromagnetic material decreases with temperature [START_REF] Hauschild | Determination of the temperature dependence of the coercivity in Fe/Cr (110) multilayers[END_REF], and it loses this property completely above the Curie temperature. The parameter k in the JA model represents the hysteresis effect, i.e., the energy required by domain walls to overcome the pinning sites. In the case of soft magnetic materials, the pinning factor k is similar to the coercivity [START_REF] Jiles | Introduction to Magnetism and Magnetic Materials[END_REF] which hints about the decrease in the value of the parameter k with the increase in temperature. The parameter k was identified from the measured B-H loops at different temperatures. For the given temperature range, a quadratic dependence (i.e., three unknowns only) of the parameter k on temperature was observed for this material and is shown in Fig. 2(e). The measured values of coercivity at different temperatures are plotted in Fig. 2(f). Both the coercivity and the parameter k show similar decreasing trends against the temperature in the studied steel.

IV. RESULTS AND DISCUSSION

The temperature evolution of the JA model parameters, described in Section III, is fitted to their respective polynomial equations, and the details of the number of unknowns for each parameter are given in Table I. The unknowns' values and the quality of the fit for all the parameters are also given. It can be seen that the proposed approach not only eases the process of identification of the JA model parameters but also helps in reducing the number of unknowns. No additional material information is needed.

The next step is to use these unknowns in their respective polynomial equations to compute the values of the JA model parameters which can then be used to predict the B-H loops of ferromagnetic materials at different temperatures. The measured and computed hysteresis loops of 35WW300 NO electrical steel are compared in Fig. 3 for different temperatures (50 °C, 150 °C, and 300 °C). We can see that as the temperature increases, the area and the width of the B -H loop decrease.

An error metric is defined to examine the accuracy of the predictions and is given in the following equation:

Error (%) = ELD computed -ELD measured ELD measured × 100% ( 3
)
where ELD is the energy loss density that is computed by calculating the area of the loop. The proposed approach gives reasonable estimates for the given temperature range, and the errors in ELD are well within 6%. The errors obtained using the proposed model are slightly higher than that of the original model, because the proposed approach approximates the JA parameters requiring fewer measurements which affects the accuracy. The same approach was applied to a similar material, i.e., B35AV1900 NO electrical steel (B H =2500 A/m = 1.68 T and P loss,B=1.5 T = 1.91 W/kg) and similar results were obtained (Fig. 4 shows a sample result at 150 O C).

A. Application Example-TEAM 32 Problem

Case 3 of TEAM 32 problem (details can be found in [START_REF] Bottauscio | A test-case for validation of magnetic field analysis with vector hysteresis[END_REF]) has been solved in 2-D using the vector JA model [START_REF] Bergqvist | A simple vector generalization of the Jiles-Atherton model of hysteresis[END_REF] in Infolytica's MagNet [16]. We have computed and used two sets of the JA parameters for 35WW300 NO electrical steels at two different temperatures (i.e., 50 °C and 300 °C) using the data given in Table I. The B-field solutions at the point x = 0 mm, y = 61.5 mm (the model is centered at the origin) in the "T" section of the ferromagnetic core are plotted for two temperatures in Fig. 5. We have applied the proposed approach in the TEAM problem to demonstrate that the proposed empirical approach can be used in thermally coupled electromagnetic simulations for the accurate design of modern electrical machines.

V. CONCLUSION

The temperature dependence of the JA model parameters was studied for NO electrical steels. A simple engineering approach has been presented to extend the JA model to incorporate the effects of temperature. The proposed approach makes use of a minimal amount of the experimental data, i.e., B-H loops measured at only three different temperatures are sufficient for model predictions within a given temperature range with reasonable accuracy. The JA parameters are identified from these measured B-H loops. The values of the parameters M S and c are obtained at the lowest temperature of interest from the measured data and are kept constant. The remaining JA model parameters (a, α, and k) are then identified at different temperatures to fit the respective polynomials.

The proposed approach offers two advantages. The identification of the JA model parameters has been simplified by reducing the number of parameters. The JA model parameter information at extremely low temperatures (i.e., 0 °K) is not required. Since this is an experimental approach, the identification of the JA model parameters using analytical formulae, as previously described in the literature, may be ignored. The proposed approach has only been tested for NO electrical steels for a limited temperature range which is a typical range of thermal operation of an electrical machine.

Fig. 1 .

 1 Fig. 1. Effect of temperature on the B-H loops measured at B max = 1.5 T. Inset: variation of coercivity with temperature. The operating frequency is 50 Hz.

Fig. 2 .

 2 Fig. 2. Variation of the (a) saturation magnetization, M S , (b) reversibility factor, c, (c) domain density, a, (d) domain coupling coefficient, α (e) pinning factor, k, and (f) measured coercivity of 35WW300 NO electrical steel at seven discrete temperatures. Red dots in (a)-(e): identified values of the JA parameters from the measured data.

Fig. 3 .

 3 Fig. 3. Measured and computed B-H loops at (a) 50 °C, (b) 150 °C, and (c) 300 °C, using the identified JA model parameters (the original JA model) and the computed JA model parameters (the proposed JA model).

Fig. 4 .

 4 Fig. 4. Measured and computed B-H loops for B35AV1900 NO electrical steel at 150 °C, using the identified JA model parameters (the original JA model) and the computed JA model parameters (the proposed JA model). The parameters of the original JA model are M S = 1.2594 × 106 A/m, c = 0.02, a = 27.1854 A/m, α = 7.7459 × 10 -5 , and k = 50.443 A/m, and for the proposed approach, M S = 1.2594 × 106 A/m, c = 0.02, a(T ) = 13.92 + 0.1189 T A/m, α(T ) = 4.8452 × 10 -5 + 2.549 × 10 -7 T, and k(T ) = 50.2 + 0.00842T --6.507 × 10 -7 T 2 A/m.

Fig. 5 .

 5 Fig. 5. Computed rotational fields at the point (marked blue) in the T section of the ferromagnetic core at 50 °C and 300 °C.

TABLE I CURVE

 I FITTING OF THE JA PARAMETERS

ACKNOWLEDGMENT

This work was supported in part by the Automotive Partnership Canada under Grant APCPJ 418901-11 and in part by the National Science and Engineering Research Council of Canada.


Siyang Deng
email: siyang.deng@centralelille.fr

Stéphane Brisset

Stéphane Clenet

COMPARATIVE STUDY OF METHODS FOR OPTIMIZATION OF ELECTROMAGNETIC DEVICES WITH UNCERTAINTY

Keywords: Reliability-based design optimization, uncertainty, reliability, safety transformer. 1

This paper compares different probabilistic optimization methods dealing with uncertainties. Reliability-Based Design Optimization is presented as well as various approaches to calculate the probability of failure. They are compared in terms of precision and number of evaluations on mathematical and electromagnetic design problems to highlight the most effective methods.

INTRODUCTION

In most optimization problems, the variables are usually considered as deterministic, i.e. without any variability. This traditional Deterministic Design Optimization (DDO) addresses only the performances but not the reliability and robustness. Since, the manufacturing process, the characteristics of materials, and the dimensions undergo variability, the device performances are in practice not deterministic but uncertain. Therefore, uncertainties related to design parameters are more and more taken into account in all engineering fields.

Various formulations are available in the literature to express optimization problems with uncertainty, which can be mainly divided into Worst-Case Optimization (WCO), Robust Design Optimization (RDO), Reliability-Based Design Optimization (RBDO) and Reliability-Based Robust Design Optimization (RBRDO). WCO is a non-probabilistic approach that is based on minimax problem formulation. For instance, [START_REF] Xiao | Considerations of uncertainty in robust optimisation of electromagnetic devices[END_REF] solves a multi-objective problem that aims to minimize the objective function, its maximum feasible value in a surrounding box, and the greatest component of the objective function derivative in the surrounding box.

The three others are probabilistic approaches that quantify the uncertainty of quantities of interest by probability distribution functions. RDO minimize a weighted sum of the mean value and variance of the objective function subject to deterministic constraints. RBDO minimize the mean value of the objective function subject to constraints on the probability of failure, i.e. constraint violation. Finally, RBRDO [START_REF] Kim | A Comparative Study on Probabilistic Optimization Methods for Electromagnetic Design[END_REF] integrates both last formulations by changing the objective function and constraints at the same time.

A comparative study [START_REF] Ren | Investigation of reliability analysis algorithms for effective reliability-based optimal design of electromagnetic devices[END_REF] of two RBDO Double-Loop Methods (DLM) with Monte Carlo Simulation (MCS) shows the interest of RBDO-DLM compared to MCS that requires a very large sampling to be accurate. However, other RBDO approaches such as Single-Loop (SLM) and Sequential Decoupled Methods (SDM) were not simultaneously investigated and this paper proposes to compare 6 algorithms belonging to the three aforementioned RBDO approaches with MCS in order to highlight the most accurate and the less time consuming. This is performed with a simple mathematical model and the multidisciplinary optimization problem of a safety transformer with uncertainty.

The paper is organized into three parts. Chapter 2 introduces the different categories of RBDO approaches. Two examples are detailed in the chapter 3 and used to compare the different methods. Last chapter is the conclusion.

RELIABILITY-BASED DESIGN OPTIMIZATION

The original formulation without any uncertainty or DDO is expressed as:

(1)

min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑔(𝑑) ≤ 0 𝑑 𝐿 ≤ 𝑑 ≤ 𝑑 𝑈
where 𝑑 is the input design variable, 𝑓(•) and 𝑔(•) are the objective function and the inequality constraint, and 𝑑 𝐿 , 𝑑 𝑈 represent the lower and upper bounds of 𝑑, respectively.

As the variability of the design variables is taken into account, the original deterministic input parameter 𝑑 should be replaced by a random input parameter 𝑋, which follows the normal law in this paper for simplicity. The mean value of 𝑋 is denoted 𝑑 and is the unknown of the new design variables, the standard deviation is denoted 𝜎 and is considered constant which means that the variability of the input parameter X doesn't depend on its mean.

The method RBDO aims to find the optimal design with the allocation of a target reliability level. The probability of failure is close to the number of points in a sampling around the mean value that fall within the reliable domain. This means that RBDO attempt to find the optimal design with the reliability that the probability of failure must be smaller than a given target value. The formulation of RBDO is as follows:

(2)

min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑃 𝑓 (𝑔(𝑥) > 0) ≤ 𝑃 𝑡 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
where 𝑃 𝑓 is the probability of failure, 𝛽 𝑡 is the reliability index, and 𝑃 𝑡 = 𝛷(-𝛽 𝑡 ) the target value for the probability of failure. The RBDO uses probabilistic constraints to make sure that the design variables satisfy a desired reliability level while minimizing the function objective.

Probability of failure

The probability is the likelihood of an event, estimated by a real number between 0 and 1. The Probability Density Function (PDF) and the Cumulative Distribution Function (CDF) define the occurrence of stochastic quantities inherently uncertain. The statistical description of a random variable 𝑋 given by the CDF Fx or PDF fX is expressed as follows:

(3)

𝑃[𝑋 ≤ 𝑥] = 𝐹 𝑋 (𝑥) = ∫ 𝑓 𝑋 (𝜏)𝑑𝜏 𝑥 -∞
where 𝑃 is the probability of occurrence of an event.

In the field of electromagnetic device manufacturing, the system ability to satisfy consumer's demand or operating constraints is important. The reliability means that designers should reduce the probability of failure as much as possible. The determination of the reliability of the system is based on the limit state function. Each constraint 𝑔(𝑋) ≤ 0 can separate the domain of 𝑋 into three parts: the limit state curve is 𝑔(𝑋) = 0, the domain where 𝑔(𝑋) > 0 is the failure domain, on the contrary the security domain represents the area 𝑔(𝑋) < 0. where 𝐷 𝑓 is the failure domain. Because the computational burden is heavy with numerous random parameters, direct integration is almost impossible and thus Monte-Carlo Simulation (MCS) or other techniques such as First-Order Reliability Method (FORM) [START_REF] Liu | Optimization algorithms for structural reliability[END_REF], [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF] is often used to calculate an approximation of 𝑃 𝑓 .

FORM and inverse FORM are based on an isoprobalistic transformation to have a normalized vector of statistically independent random variables 𝑈 instead of the initial input parameter 𝑋.

For the Gaussian vector 𝑋 in this paper, the transformation 𝑇 is as follows: [START_REF] Hasofer | Exact and invariant second-moment code format[END_REF] 𝑇: 𝑈 = (𝑋 -𝑑) 𝜎 ⁄ Figure 2. The first order reliability method [START_REF] Lopez | Reliability-based design optimization strategies based on FORM: a review[END_REF].

Then the limit state function changes from 𝑔(𝑑) = 0 to 𝐺 𝑈 (𝑢) = 0 and the CDF from 𝐹 𝑋 to 𝐹 𝑈 . 𝐺 𝑈 (𝑢) is defined as the performance function, and the FORM or inverse FORM method uses a linear approximation to replace the real performance function at the Most Probable Point of failure (MPPF) in U-space. The MPPF 𝑢 * is the one that minimize the distance between the origin 𝑂 and 𝐺 𝑈 (𝑢) = 0. After 𝑢 * is found, the limit-state function is replaced by a tangent hyperplane crossing 𝑢 * . So the probability of failure is calculated by: ( 6)

𝑃 𝑓 ≈ ∫ 𝑓 𝑈 (𝑢)𝑑𝑢 = 𝛷(-𝛽) 𝐺 ̃𝑈(𝑢)>0
where 𝐺 ̃𝑈(𝑢) = 0 is the hyperplane that approximates the limit-state function 𝐺 𝑈 (𝑢) = 0, 𝛷 is the standard Gaussian cumulative distribution function, and 𝛽 = ‖𝑢‖ is the reliability index which is equal the norm of 𝑢.

The formulation ( 2) is the basis of all RBDO methods and can be solved by different approaches which are usually separated in three main categories: double-loop, single-loop, and sequential decoupled methods. The following sections introduce the principle of these different categories and present some approaches for each of them.

Double-loop methods

Double-loop methods use two loops to solve RBDO problems: the inner loop aims to analyze reliability of the chosen configuration and to calculate the probability of failure using FORM or inverse FORM; the outer loop seeks the mean values of input designs that minimize the objective function and constrain the probability of failure computed by inner loop.

There are several approaches for double-loop methods, the most popular are Reliability Index Approach (RIA) and Performance Measure Approach (PMA) [START_REF] Ren | Investigation of reliability analysis algorithms for effective reliability-based optimal design of electromagnetic devices[END_REF], [START_REF] Tu | A new study on reliability-based design optimization[END_REF].

Reliability Index Approach (RIA)

RIA uses the FORM to calculate the reliability index in the inner loop:

(7) 𝛽 = min 𝑢 ‖𝑢‖ 𝑠. 𝑡. 𝐺 𝑈 (𝑢) = 0
The outer loop of RIA, any constrained non-linear algorithm like SQP or others can be chosen to minimize the objective function 𝑓 with the index 𝛽 given by the FORM:

(8) min 𝑑 𝑓(𝑑) 𝑠. 𝑡. 𝑔 ≤ 0 𝛽 ≥ 𝛽 𝑡 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
where 𝑔 < 0 is used to restrict the because the definition in equation ( 6) only validate if the origin is located in the security domain.

Performance Measure Approach (PMA)

In former formulation of RBDO, optimization is carried out with the limitation of the reliability index that must be greater than or be equal to a target value. The calculation of this index leads to the search for the MPPF. In contrast, for PMA formulation, optimization is formulated with the limitation of maximum performance with a given value of reliability index. The search of this maximum performance is to maximize the function 𝐺 𝑈 with the limitation that reliability index must achieve the target value. This approach is considered as the inverse of the FORM approximation [START_REF] Wu | Advanced probabilistic structural analysis method for implicit performance functions[END_REF], [START_REF] Youn | Hybrid analysis method for reliability-based design optimization[END_REF].

So the outer loop becomes: where 𝑢 * is the Maximum Performance Target Point (MPTP) that corresponds to the target index 𝛽 𝑡 .

Single-loop methods

For the so-called single-loop methods, the main point is that the inner loop is replaced by an approximation to avoid the iterative evaluations for reliability analysis in order to accelerate the convergence to the optimum.

Approximate Moment Approach (AMA)

AMA is based on statistical moments. The first order Taylor expansion is used to calculate the mean 𝜇 𝑔 𝑖 and the standard deviation 𝜎 𝑔 𝑖 for all constraint functions 𝑔 𝑖 (𝑑) by using the following expressions [START_REF] Putko | Approach for input uncertainty propagation and robust design in CFD using sensitivity derivatives[END_REF]:

(11) { 𝜇 𝑔 𝑖 = 𝑔 𝑖 (𝑑) 𝜎 𝑔 𝑖 2 ≈ (𝛻 𝑔 𝑖 𝑇 ) 2 • 𝜎 2
where 𝜎 is the standard deviation of the input variables 𝑋 and 𝛻 is the gradient operator. With these expressions, an approximation 𝑔 𝑖 (𝑑) + 𝛽 𝑡 𝜎 𝑔 𝑖 is used to replace the original constraints 𝑔 𝑖 (𝑑′) as shown in Figure 3. In this way, the margin of security is shifted in order to keep the reliability as desired. Figure 3 shows the basic idea of AMA approach, 𝑔(𝑑) = 0 is the deterministic limit-state function. As the optimum point 𝑑 * is on the limit-state curve, the probability of failure is about 50%. The situation which could cause the failure around this point with a given probability is shown as shaded area on the left side. This approach is based on a local linear approximation of the constraint functions at the mean value of the design parameters and as the probabilistic distribution of the performances depends on two moments that are also approximated, the reliability assessment could produce significant numerical error [START_REF] Youn | Selecting probabilistic approaches for reliability-based design optimization[END_REF].

Single Loop Approach (SLA)

SLA uses the same approximation as AMA for calculating the moments and the position of MPPF (marked as 𝑑 ′ * in Figure 3 and 𝑥 𝑝 in Figure 4). The difference between these two methods is that SLA evaluates the constraints 𝑔 at the approximate MPPF unlike AMA that uses a first-order approximation of the constraint function around the mean value 𝑑 directly. The formulation is given as [START_REF] Liang | A single-loop method for reliability-based design optimization[END_REF]: where 𝛼 𝑖 is the normalized gradient of the 𝑖th constraint, 𝑥 𝑝 𝑖 is the approximation of the MPPF for the constraint 𝑔 𝑖 , and ∘ is the Hadamard operator (element-wise) product.

Figure 4. Principle of the SLA approach

The principle of this approach is similar to AMA. As shown in Figure 4, it does not search for the MPPF of each constraint by using an inner loop but approximate its position.

Sequential decoupled methods

Sequential decoupled methods aim to change the initial problem into a series of optimization cycles. The cycles are sequential, each individual optimization is deterministic and uses the optimum given by the former optimization as an initial point. At the first iteration, the algorithm searches a deterministic solution without considering uncertainty and then compute the reliability index of this solution to deduce a shift in order to achieve a given probability of failure. The next iterations refine the shift.

Sequential Optimization and Reliability Assessment (SORA)

SORA employs a series of cycles of deterministic optimizations and reliability assessments. In each cycle, optimization and reliability assessment are decoupled from each other, the reliability assessment is only conducted after the deterministic optimization to verify constraint feasibility under uncertainty. The main point of this method is to shift the boundaries of violated constraints to the feasible direction based on the reliability information obtained in the former cycle [START_REF] Du | Sequential optimization and reliability assessment method for efficient probabilistic design)[END_REF]. The updated point is used in the next cycle of the deterministic optimization. This cycle is repeated until the fulfillment of the convergence criteria.

The process of SORA method is presented as follows. First of all, an initial shift 𝑠 0 = 0 allows finding the solution 𝑑 * 0 of the deterministic problem. At each iteration, the optimization problem is defined by: where 𝑘 indicates the iteration number. The optimal value is set as 𝑑 * 𝑘 . After each optimization, some of the constraints may become active. For an active constraint, the optimal point 𝑑 * 𝑘 is on the limit-state curve. When considering the randomness of 𝑋, the actual probability of failure is about 0.5, so a reliability assessment is implemented at the deterministic optimum solution to locate the MPTP 𝑥 𝑘 that corresponds to the desired probability of failure. To ensure the MPTP onto the deterministic boundary, a shifting vector 𝑠 𝑘+1 is derived: [START_REF] Tran | A Benchmark for Multi-Objective, Multi-Level and Combinatorial Optimizations of a Safety Isolating Transformer[END_REF] 𝑠 𝑘+1 = 𝑑 * 𝑘 -𝑥 𝑘 Therefore, when establishing the equivalent deterministic optimization model in the next cycle, the constraints is modified to shift the MPTP at least onto the deterministic boundary.

Sequential Approximate Programming (SAP)

SAP is another approach based on Taylor expansion at the first order. The original optimization problem is decomposed into a sequence of sub-optimization problems. Each sub-optimization consists of an approximate objective function subjected to a set of approximate constraint functions [START_REF] Yi | A sequential approximate programming strategy for performance-measure-based probabilistic structural design optimization[END_REF]. The details are presented as follows:

For considering the PMA formed optimization problems, the expressions are as equation ( 8), following the idea of SAP, a sequential approximate formulation is constructed as: where 𝐺 ̂𝑝(𝑑 𝑘-1 ) = 𝐺 𝑢 (𝑢 𝑘-1 ) is the approximate probabilistic performance measure and 𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 ) is its gradient. To avoid another optimization to calculate 𝑢 𝑘 , it is updated by the function below:

(19) 𝑢 𝑘 = -𝛽 𝑡 𝛻 𝑈 𝐺 𝑈 (𝑢 𝑘-1 )
‖𝛻 𝑈 𝐺 𝑈 (𝑢 𝑘-1 )‖

where 𝑢 0 = 0 is usually chosen as the initial estimation.

With the first order approximation of 𝐺 𝑝 (𝑑 𝑘 ), new formulation is obtained:

(20)

min 𝑑 𝑓 𝑘 (𝑑) 𝑠. 𝑡. 𝐺 ̂𝑝(𝑑 𝑘-1 ) + (𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 )) 𝑇 (𝑑 -𝑑 𝑘-1 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎
This formulation of SAP is similar to the outer loop of PMA. The same principle is used to convert a RIA formulation into a sequential approximate programming [START_REF] Cheng | A sequential approximate programming strategy for reliability-based structural optimization[END_REF].

CASE STUDIES

In this chapter, two examples are tested to compare the performances of all approaches.

Numerical example

To assess the efficiency of these methods, the simple two-variable problem in [START_REF] Kim | A Single-Loop Strategy for Efficient Reliability-Based Electromagnetic Design Optimization[END_REF] is analyzed. The optimization problem of this numerical example is: In order to understand the implications of the choice of different formulations and also to quantify the impact of the choice of approaches and algorithms on the accuracy and the number of evaluations, the results of the 7 mentioned methods are given in Table 1. The different optima are compared with each other in Figure 5, where the dotted curves present contours of 𝑓 and constraint boundaries 𝑔 𝑖 = 0 are depicted as solid lines. The target probability of failure 𝑃 𝑡 for RBDO is 2.28%, 𝑃 𝑓 is computed by Monte-Carlo simulation with a sampling of 10 6 realizations. For the simplicity and without loss of generality, all the uncertainties are modeled with the normal law. The column named convergence in Table 1 means that we run the algorithms 100 times with different initial points to see how many times it converges to the same optimum. From Table 1 and Figure 5, we can see that most of the methods can find a result that satisfies the constraints. DDO leads to a solution with the lowest objective function mean value, the highest standard-deviation of objective function and a probability of failure around 50%. On the contrary, the probabilities of failure of RBDO methods are close to the target probability except for the single-loop method AMA because of the approximation used to increase the speed of convergence by sacrificing the accuracy. Among RBDO methods, single-loop strategies are fast but inaccurate, double-loop and sequential decoupled methods lead to the same results but sequential decoupled methods are greatly faster.

Example of safety transformer optimization

The example of a safety isolating transformer [START_REF] Tran | A Benchmark for Multi-Objective, Multi-Level and Combinatorial Optimizations of a Safety Isolating Transformer[END_REF] is also investigated. This is a single-phase transformer with grain-oriented E-I laminations. The primary and secondary windings are both wound around the frame surrounding the central core. The model inputs consist in 7 random design parameters: three parameters 𝑎, 𝑏, 𝑐 for the shape of the lamination, one for the frame 𝑑, two for the section of conductors 𝑆 1 , 𝑆 2 and the last one for the number of primary turns 𝑛 1 . These design variables are shown in Figure 6. The range of input variables are:

(22) 3𝑚𝑚 ≤ 𝑎 ≤ 30𝑚𝑚 14𝑚𝑚 ≤ 𝑏 ≤ 95𝑚𝑚 6𝑚𝑚 ≤ 𝑐 ≤ 40𝑚𝑚 10𝑚𝑚 ≤ 𝑑 ≤ 80𝑚𝑚 200 ≤ 𝑛 1 ≤ 1200 0.15𝑚𝑚 2 ≤ 𝑆 1 ≤ 19𝑚𝑚 2 0.15𝑚𝑚 2 ≤ 𝑆 2 ≤ 19𝑚𝑚 2
There are 8 inequality constraints in this problem. The copper and iron temperatures 𝑇 𝑐𝑜𝑛𝑑 and 𝑇 𝑖𝑟𝑜𝑛 should be less than given temperatures 120℃ and 100℃ respectively. Both the magnetizing current should be less than 10%. The filling factors of the primary coil 𝑓 1 and secondary coil 𝑓 2 should be lower than 1, and the efficiency 𝜂 should be larger than 0.8. At last, the residue must be less than 10 -6 .

The constraint functions are: The aim is to minimize the total mass 𝑀 𝑡𝑜𝑡 with a target probability 𝑃 𝑡 = 0.13%, so that the target reliability index is 𝛽 𝑡 = 3.

Table 2 shows the results of some methods. For each constraint, there are two probabilities of failures, upper ones are calculated by the methods themselves and the below ones are calculated by the Monte-Carlo simulation for comparison. For the reason that the 8 th constraint is a residue, we consider it as a deterministic one, so there are only 7 probabilities of constraints presented in the table 2. Note that PMA can find an optimum that satisfies all the constraints even if the number of evaluations of the model is high. For SLA, it has a smaller number of evaluations but there is one constraint violated. The reason is that SLA sacrifices the accuracy in order to reduce the number of evaluations, leading to a coarse computation of the probability of failure. The convergence rate of SORA is not as good as the other two but it has the smallest number of evaluations among them. It can be seen that this number is nearly 15 times less than PMA and 3 times less than SLA. Unfortunately, the rate of convergence for SORA is 4 times lower. This means that a multi-start process is required and the number of evaluations will increase consequently. We obtain almost the same conclusions as with the numerical example: the single-loop method SLA is the most inaccurate method and doubleloop method PMA has the highest convergence rate. For this more complicated example, the number of evaluations of sequential decoupled method SORA is even smaller than of SLA, so that the more complex the model is, the more effective SORA may be. Other methods fail to find a solution, probably because this problem is hardconstrained and the solution of DDO is on the limit-state of four constraints. So, it also indicates that not all the aforementioned approaches can handle complicated models.

CONCLUSION

RBDO methods change the initial constraints into probabilistic ones and use different approaches to approximate the limit state function or the most probable point of failure to calculate the probability of failure or the reliability index.

The mathematical example shows that the 6 RBDO approaches have almost the same results except AMA that is less accurate. The optimization of the safety transformer highlights that not all the methods can converge to the global solution. PMA, SLA, and SORA appear to be more stable. Considering both numerical examples, SORA is the most effective method among all RBDO approaches.

Figure 1 .

 1 Figure 1. The failure domain, security domain, and limit-state curve. The probability of failure is the probability of the event 𝑔(𝑋) > 0. It is calculated with the integral: (4) 𝑃 𝑓 = 𝑃[𝑔(𝑋) > 0] = 𝑃[𝑋 ∈ 𝐷 𝑓 ] = ∫ 𝑓 𝑋 (𝑥)𝑑𝑥 𝐷 𝑓

  . 𝐺 𝑝 (𝑢 * ) ≤ 0𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎where 𝐺 𝑝 (𝑢 * ) is the maximal performance measurement. The purpose of the inner loop is to find 𝑢 * in U-space such as:[START_REF] Youn | Selecting probabilistic approaches for reliability-based design optimization[END_REF] 𝐺 𝑝 (𝑢 * ) = max 𝑢 𝐺 𝑈 (𝑢) 𝑠. 𝑡. ‖𝑢‖ = 𝛽 𝑡

Figure 3 .

 3 Figure 3. Principle of the AMA approach

  . 𝑔 𝑖 (𝑥 𝑝 𝑖 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎 with (14) 𝑥 𝑝 𝑖 = 𝑑 -𝛽 𝑡 𝛼 𝑖 ∘ 𝜎 𝛼 𝑖 = 𝛻𝐺 𝑈 𝑖 (𝑑 ) ‖𝛻𝐺 𝑈 𝑖 (𝑑 )‖ = 𝜎∘𝛻𝑔 𝑖 (𝑑 ) ‖𝜎∘𝛻𝑔 𝑖 (𝑑 )‖

  . 𝑔(𝑑 𝑘 -𝑠 𝑘 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎

  . 𝐺 𝑝 (𝑑 𝑘 ) ≤ 0 𝑑 𝐿 + 𝛽 𝑡 • 𝜎 ≤ 𝑑 ≤ 𝑑 𝑈 -𝛽 𝑡 • 𝜎The approximate function 𝐺 𝑝 (𝑑 𝑘 ) is built with a first order Taylor expansion with respect to the design variables at the current point:(18)𝐺 𝑝 (𝑑 𝑘 ) ≈ 𝐺 ̂𝑝(𝑑 𝑘-1 ) + (𝛻 𝑑 𝐺 ̂𝑝(𝑑 𝑘-1 )) 𝑇 (𝑑 𝑘 -𝑑 𝑘-1 )

Figure 6 .

 6 Figure 6. 3D geometric simulation of the transformer

Table 1 . Results for different optimizations

 1

	Methods	𝑑 1	𝑑 2	𝜇 𝑓	𝜎 𝑓	Max of 𝑃 𝑓 (%)	Convergence (100 times)	Function evaluations of 𝑓 and 𝑔
	DDO	2.4398	0.8400	-2.6266	0.2072	50	20%	84
	RIA	2.2513	1.9691	-1.9945	0.1792	2.28	29%	4539
	PMA	2.2513	1.9691	-1.9945	0.1792	2.28	82%	3183
	AMA	2.3524	2.1568	-1.8712	0.1715	1.87	46%	132
	SLA	2.2512	1.9677	-1.9953	0.1793	2.29	32%	165
	SORA	2.2513	1.9691	-1.9945	0.1792	2.28	45%	531
	SAP	2.2513	1.9691	-1.9945	0.1792	2.28	37%	181

Figure 5. The results of different methods plotted in search space.

Table 2 . Results on the safety transformer optimization problem

 2

		Transformer		DDO	PMA	SLA	SORA
		𝑎 (𝑚𝑚)		13.004	13.078	13.045	13.078
		𝑏 (𝑚𝑚)		50.1	52.274	51.613	52.278
		𝑐 (𝑚𝑚)		16.537	16.758	16.837	16.757
		𝑑 (𝑚𝑚)		43.05	42.182	42.603	42.180
		𝑛 1 𝑆 1 (𝑚𝑚 2 ) 𝑆 2 (𝑚𝑚 2 )		639.76 0.3238 2.9026	658.36 0.3256 2.9296	653.56 0.3254 2.9273	658.39 0.3256 2.9296
		𝜇 𝑀 𝑡𝑜𝑡		2.3112	2.3546	2.3534	2.3546
		𝜎 𝑀 𝑡𝑜𝑡		0.0079	0.0081	0.0081	0.0081
		Convergence		0.12	0.84	0.72	0.23
	𝑃 𝑓 1 (𝑇 𝑐𝑜𝑛𝑑 > 120℃) (%)	by MC by itself	0 0	0 0	0 0	0 0
	𝑃 𝑓 2 (𝑇 𝑖𝑟𝑜𝑛 > 100℃) (%)	by MC by itself	50 50.658	0.13 0.1337	0.13 0.3754	0.13 0.1389
	𝑃 𝑓 3 (∆𝑉 2 𝑉 20 ⁄	> 0.1) (%)	by MC by itself	0 0	0 0	0 0	0 0
	𝑃 𝑓 4 (𝐼 10 𝐼 1 ⁄ > 0.1) (%)	by MC by itself	50 50.203	0.13 0.1356	0.13 0.1357	0.13 0.1370
	𝑃 𝑓 5 (𝑓 1 > 1) (%)	by MC by itself	50 49.987	0.13 0.1323	0.13 0.1362	0.13 0.1388
	𝑃 𝑓 6 (𝑓 2 > 1)(%)	by MC by itself	50 50.021	0.13 0.1415	0.13 0.1368	0.13 0.1382
			by MC	50	0.13	0.13	0.13
			by itself	49.757	0.1324	0.1376	0.1390
	𝑃 𝑓 7 (0.8 -𝜂)(%)
	Model evaluation		242	30607	6048	2171