Keywords: Optimisation mixte en nombres entiers, optimisation sans dérivée, simulation boîte noire, symétrie cyclique distance de collier, plan d Mixed integer non-linear programming, derivative free optimization, black-box simulation, cyclic-symmetry necklace distance, design of experiments, trust region. vi

Enfin, nous donnons des pistes d'amélioration de la méthode d'optimisation proposée pour une meilleure exploration de l'espace des variables de conception pour éviter d'être piégé dans des minima locaux.

First of all, I am grateful to my advisor Dr. Delphine Sinoquet for her constant guidance and support. Her passion for applied mathematics, especially optimization was an inspiration for me during my Ph.D. and fed my curiosity with state-of-the-art optimization techniques, supporting me technically and steering me in the right direction whenever she thought I needed it. I would also like to thank my director, Professor Marcel Mongeau, and my co-advisor, Dr. Sébastien Da Veiga, for their great mentorship in technique problems and research. Their enriching advice gave during the analysis of my numerical results part, or the state-of-the-art studying, or proofreading the manuscript is invaluable to my thesis.

I would like to thank the members of my thesis jury for having accepted to evaluate my work. In particular, Professor Charles Audet and Professor Laurent Dumas for their time of re-lecture my thesis in detail and their constructive comments. I also thank Professor Jean-Baptiste Hiriat-Urruty, Professor Ana Luísa Custódio, and Dr. Merlin Keller for the interest they showed by accepting to participate in this jury. I spent most of the time in Rueil Malmaison, surrounded by my colleagues in the department of applied mathematics where I quickly felt intimated. This integration was facilitated by Benjelloun-Touimi Zakia, the head of the department, who was always benevolent, who knew how to combine work and life. I especially mention Professor Quang Huy Tran, who had always answered all the questions regarding my thesis, conferences, visa, etc. My Vietnamese friends at IFP Energies Nouvelles: Duc Thach Son Vu, Quang Duy Le, Tri Dat Ngo, Manh Quan Nguyen, Ngoc Tu Trinh, Quang Long Nguyen, Van Quan Hoang, for their experiences sharing, life caring...they are all like my brothers. I also mention other colleagues with their sincerely constant help during my Ph.D.: Miguel Munoz Zuniga, Frederic Delbos for their valuable discussions about some parts of the work. Sylvie Hoguet, department assistant, for her help to fulfill the order missions, my salary time sheets, etc. Dominique Sfez, my proprietary, and her husband Pierre Doguin, for their support. Sylvie Pegaz-Fiornet, Henry Eon, Francoise Willien, Sylvie Wolf, Guillaume Enchery... for their great integration during our lunchtime, break times... and sometimes tried to teach French to me!! And of course, I can not forget to mention the former and current Ph.D. students in my department: Bastien Hamlat, Zakaria Jorti, Karine Laurent, Julien Coulet, Alexis Cousin, Sabrina Bassetto, Dongmo Nguepi Guissel Lagnol, Joelle Ferzly, Jingang Qu, Mohamed El Maarouf Aboul Karim, with them I had lots of great after-work times together. I would like to thank DSI, Security, cleaning, restaurant staff at IFP Energies Nouvelles for their informative support, safety insurance, making a clean working environment, and providing a variety of nutritious lunches.

iii I had the chance to visit ENAC, where I found a friendly environment while working with others. I would like to thank Ph.D. students at ENAC: Mahfoud Herraz, Sana Ikli, Hasnae Kasmi for their enthusiastic discussions and update ENAC's information for me. In particular, I mention Professor Nicolas Couellan, and former Ph.D. student Gabriel Jarry for their ideas that contributed to the reformulation part of the thesis.

During these three years, there were a huge number of administrative procedures that I had to do. As an international student, I would not be able to finish all those tasks without the help of Virginie Lemire, HR at IFP Energies Nouvelles, Helene Weiss, Optim team's assistant, Agnes Requis, EDMITT staff. I would also thank to Safran tech, the funding granted by them has allowed the development of this thesis.

I express my sincere gratitude to my boyfriend, Anshul Paigwar, a robotic researcher at INRIA Grenoble for his constant encouragement and his love. He is not only my boyfriend but also my teacher in life who has helped me to improve my communication skills, my fashion style, my general knowledge, and presenting skills. I would also thank my friends: Thi Thanh Yen Nguyen, Quoc Huy Vu,... for the time that we have shared.

Finally, I must express my very profound gratitude to my parents, my siblings, my nephews, my nieces for providing me with unfailing support. This accomplishment would not have been possible without them. Thank you.

Rueil Malmaison, 1 st , August, 2021,

Tran Thi Thoi

Résumé

Ces dernières années, un nombre considérable d'applications industrielles réelles impliquant des variables mixtes et des simulateurs coûteux en temps ont été réalisées, par exemple, chez Safran Tech et IFPEN, des conceptions optimales de turbine de moteur d'hélicoptère, de lignes d'amarrage d'éoliennes offshore, de stators et de rotors de moteurs électriques,. . . Dans ces problèmes d'optimisation non linéaire, les dérivées de la fonction objectif (et, éventuellement des fonctions contraintes) ne sont pas disponibles et ne peuvent pas être directement approximées. Une autre difficulté est que ces problèmes impliquent des variables de nature hétérogène: un nombre variable de composants (variables entières), différents matériaux (variables catégorielles, généralement non ordonnées), la présence ou non de certains composants (variables binaires) et des variables continues décrivant les dimensions/caractéristiques des pièces de la structure.

Cette thèse a pour but de développer et d'adapter des méthodes d'optimisation sans dérivées (ou DFO pour Derivative-Free Optimization) applicables à différents types d'applications, dont la conception optimale des moteurs d'hélicoptère. Dans la première partie, nous nous concentrons sur le développement et l'adaptation d'une méthode DFO aux problèmes avec variables mixtes continues et discrètes présentant une symétrie cyclique, caractéristiques présentes dans le problème d'optimisation des pales d'une turbomachine de moteur. À cette fin, nous introduisons une distance basée sur les colliers (necklace distance) et adaptons une distance basée sur les colliers pour une méthode d'optimisation du type région de confiance. Avant d'appliquer notre méthode à un cas applicatif simplifié fourni par Safran, nous construisons un ensemble de fonctions tests issues de la littérature que nous adaptons pour obtainir un emsemble de problèmes mixtes à symétrie cyclique. Notre méthode est évaluée sur ces cas tests et comparée à des méthodes d'optimisation sans dérivées de l'état de l'art. Nous donnons également une preuve de convergence locale de notre méthode adaptée.

Dans la deuxième partie, nous nous consacrons à la planification d'expériences dans un espace mixte (variables continues et discrètes) en étendant à cet espace mixte les approches basées sur des méthodes à noyaux pour estimer des distributions de probabilité. Cette partie de la thèse est motivée par le besoin d'améliorer la phase d'initialisation de l'algorithme d'optimisation. Le but étant de permettre une meilleure exploration de l'espace des variables mixtes, guidée par les informations a priori disponibles (types de variables, symétrie, corrélations,...). Nous illustrons également le potentiel de l'approche proposée dans le cadre plus classique de l'approximation d'une fonction par un méta-modèle pour des variables mixtes continues et discrètes mais aussi pour des séries temporelles. lytical problems for the five solvers (successful runs only) together with percentage of successes (in parentheses) . . .

The eight individual one-dimensional functions for DS-MM example. The global minimum is located at x " 0.808, y " p0, 0, 1, 0, 1q, the objective-function value is

33

Ackley function with two continuous variables and four binary variables (corresponding 6 levels: y 1 " p0, 0, 0, 0q, y 2 " p0, 0, 0, 1q, y 3 " p0, 0, 1, 1q, y 4 " p0, 1, 0, 1q, y 5 " p0, 1, 1, 1q, y 6 " p1, 1, 1, 1q

Introduction

" It is through science that we prove, but through intuition that we discover.

-Henri Poincaré Contents This thesis stems from the collaboration between IFP Energies Nouvelles, Safran Tech and École Nationale de l'Aviation Civile (ENAC). The aim is to study and provide an algorithm to solve mixed continuous and discrete non-linear problems arising from industrial applications.

One of the motivating applications of the thesis is the optimal design of the turbine blades of a helicopter engine. The objective is to maximize the compressor efficiency under some stability constraints (minimizing vibrations). This involves very expensive mechanical simulations (solid and fluid mechanics), typically involving computer runs of the order of several hours for evaluating one design configuration set. Some continuous parameters describe the blade shapes, e.g., the thickness, length of the blades, and binary variables locate the different types of blade geometries that are considered. The lack of analytical form of the objective function, the appearance of discrete variables as well as the symmetry property provides challenging but also interesting problems.

This chapter is structured as follows: Section 1 introduces Safran's application context. In Section 2, we describe the idea of the reduced-order model which is currently used in practice. Section 3 presents the main contributions of the work. Finally, the outline of the thesis is presented followed by the list of the publications.

The motivating application

In the context of turbomachine design, manufacturing operations possibly affect the structural cyclic symmetry and can increase the forced response amplitudes. As a result, vibration appears in the system and destroys the engine. Therefore, minimizing vibration to obtain a stable engine is important.

Motivation

Air traffic is one of the most important means of transport nowadays, especially in Europe. There are around 8000 daily flights. Moreover, the amount of people travelling by airplane increases every year by around 5 %. Air traffic is associated with very high costs of fuel [3, 53], and also with a huge budget for the maintenance and manufacturing of the new engines.

Turbomachines

There are several ways to optimize cost in aviation: through optimizing trajectories, arrangements of passengers (through the design of seat arrangement), cargo storage, etc. Our study concentrates on optimizing the design of turbomachines by maximizing the efficiency (compressor) and by minimizing vibrations. Turbomachines Figure 1 or gas turbines are complex systems that are used in the aerospace, automotive, and power generation industries. Blades are important components that allow the exchange of energy with the flow. During the operation of the turbomachine, vibrations occur, mainly due to the excitation by modification of the aerodynamic flow and also resulting from a coupling between the flow and the movement of the blades [START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF].

Fig. 1.: Illustration of a turbomachine (from [START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF]).

In the particular application proposed by SAFRAN, the main idea is to optimize the so-called tuning and mistuning blade shapes in the engine turbine, and to find the optimal distribution of the resulting two shapes on the disk. The objective is to minimize the vibration of the compressor by changing the shapes of the compressor blades. Here, a single objective function evaluation may require several hours of computation time. This optimization problem involves a vector, x P R m , of m continuous variables, each of which describes one blade shape parameter, such as the thickness or the length of the blades. There are also integer variables that locate pre-defined possible blade geometries around the disk, as in [START_REF] Choi | Optimization of intentional mistuning patterns for the reduction of the forced response effects of unintentional mistuning: Formulation and assessment[END_REF]. In this study, we focus on the case involving only two different blade geometries; if we consider n blades, their relative positions are indicated with a binary vector y P t0, 1u n , where y i indicates whether the ith blade is of a given type (a) or of the other type, (b), for i " 1, 2, . . . , n.

Optimization challenges

Due to the complexity of fluid dynamic systems, we are not able to obtain an analytical form of the objective function, as well as physical-constraints functions. Therefore, the optimization problem that we are dealing with is a derivative-free optimization or blackbox simulator problem. Derivative-free optimization is a field of applied mathematics that concerns optimization problems without any information of derivatives of the objective function. We list in the following some mathematical challenges that must be addressed:

• No information of derivatives of the objective function and, possibly of some physicalconstraint functions.

• There is a combinatorial problem linked to the presence of discrete variables.

• The cyclic symmetry property of the problem yields a huge number of redundant arrangements.

Regarding this cyclic-symmetry challenge, remark that: two blade disks that differ only by a rotation of the pattern around the disk not only lead to a same value of the objective function, but also correspond to identical compressors. The number of equivalent solutions also rapidly increases with the number of blades, n, as illustrated in Table 3. Particularly, in this type of application, one simulation can cost several hours or more. Therefore, we aim to avoid recomputing during the optimization process, the costly objective-function value at equivalent configurations.

Tab. 3.: Number of distinct solutions for n blades (from [START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF]). The next section presents the Reduced Order Model (ROM) to reduce the computational cost of optimization for such highly combinatorial problems.

Reduced Order Model (ROM)

Two blade disks that differ only by a rotation of the pattern around the disk lead to the same value of the forced response on the disk. Such arrangements should obviously be considered as identical solutions. Recall that when the number of blades is increasing, the number of such equivalent solutions rapidly increases (Table 3).

For an instance of our problem involving n blades, we can roughly approximate the number of distinct arrangements by 2 n n . Reducing these redundancies will significantly reduce the combinatorial complexity of the optimization problem and, consequently, the computational cost of the envisaged optimization approach. Therefore, the methodology based on reduced-order modeling method is proposed. Suppose we have two types of blades, A and B, that need to be located at n locations on the disk. In [START_REF] Han | Optimization of intentional mistuning patterns for the mitigation of the effects of randon mistuning[END_REF], the authors present a physical discussion about the problem of redundancy: in the presence of strong coupling, i.e., if there are few switches (changes from A to B of from B to A) or if there exists a series of groupings of two consecutive blades of the form AA or BB (e.g., 7B5A: 7 consecutive blades B and 5 consecutive blades A), the vibration responses between blades tend to be globally uniform at the wheel. While, in the case of low coupling, i.e., if there are many switches or if the basic grouping forming the pattern are no longer AA and BB but rather AB and BA (e.g., 1A1B2A2B), the responses tend to be located on some blade neighbors.

If the optimal distribution of the blade geometry presents few alternations of the two different shapes, it can roughly be described by a distribution of patterns of the types AA or BB. On the other hand, if the optimal distribution presents many alternations of the two different shapes, one should rather use patterns AB or BA. In both cases, the idea is to consider the distribution by groups of two blades in view of obtaining a model of lower dimension that deals with n 2 groups instead of n blades. The idea of ROM is to consider n{2 groups of two blades with two possibilities for each group. We group either patterns AA and BB if there are few switches (high coupling), or AB and BA in the case where there are many switches (low coupling). For illustration, a disk with 10 blades consists of 108 distinct distributions of two types of blades, whereas a ROM approach leads to a problem involving 5 groups and only 8 distinct distributions.

A ROM optimization methodology consists of two main steps. It first optimizes in the two reduced spaces for patterns AA{BB, and then for patterns AB{BA considering on groups of 2 blades in order to provide good initial guesses for the orginal n blades problem.Then, a local search limited ROM reduces the dimension of the discrete space, which in turn tends to reduce the computation time required for each sub-problem. A disadvantage however is that it is likely to remove a large number of possibilities which will not be explored, leading thereby to a sub-optimal solution. For instance, for the problem involving n " 12 blades, ROM removes 324 disk configurations. Out of 352 distinct configurations as each of the two subproblem considers only 14 distinct configurations. This motivates our study to look for another method to address the abovementioned cyclic symmetry property to ensure that all the distinct arrangements are considered. In Chapter 5 of Part I, we will propose an adapted mathematical optimization method for this type of problem.

Thesis contributions

The three main contributions of the thesis are:

• Development of a new mixed continuous and discrete variables derivative-free optimization method based on trust-region method, [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] for cyclic-symmetry problems, called DFObd neck . The originality of our approach is the introduction of the necklace distance which avoids recomputing expensive simulations corresponding to the equivalent configurations.

• A theoretical convergence proof of the mixed-integer DFO method (evoked in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF] but never published) is provided. This proof is then further extended into a convergence proof of DFOb-d neck .

• A new design of experiments method based on a kernel-embedding probability distribution for mixed discrete variables problems.

A secondary contribution includes the proposition of benchmark cyclic-symmetry test functions and numeric comparison of our optimization approach with state-of-the-art methods such as NOMAD, RBFOpt, and DFLBOX.

Outline of the thesis

We begin by recalling in Chapter A the notation and some mathematical background necessary for the readers to follow the thesis. The rest of the thesis is composed of three parts that are detailed as follows:

Derivative-free mixed binary necklace optimization. In this first part, we focus on optimization methods for solving derivative-free problems:

• Chapter 2 gives an overview of blackbox optimization methods for mixed discrete variables problems.

• Chapter 3 carries out a review of the continuous DFO trust-region method.

• Chapter 4 focuses on the derivative-free trust-region methods for mixed binary problems.

Our contribution is the convergence proof of the original method proposed in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF].

• Chapter 5 presents our two main contributions. First, we introduce a new optimization method for cyclic-symmetry problems, named DFOb-d neck which is based on DFOb-d H .

Our second contribution, is the proposal of a set of benchmark functions adapted to cyclic-symmetry problems. These functions were used as tests for our method and for state-of-the-art methods. A simplified simulation provided by Safran is also tested.

• Chapter 6 shows an analysis result of methods about the impact of number of points in the interpolation set.

Design of experiments in mixed continuous and discrete variables space. Blackbox modelbased methods can be significantly impacted by the initial design (the points at which the objective function is evaluated before optimization). The choice of the initial design has a strong influence on how well the search domain will be explored.

• Chapter 7 reviews about design of experiments methods in continuous space and in mixed discrete space.

• Chapter 8 presents one of the main contributions of the thesis: the introduction of two extended kernel-embedding probability distribution DoEs to mixed discrete variables space. We also give examples with a large range of application types.

Escaping from local minima. This part is dedicated to escaping from local solutions by applying restart techniques.

Publications

The work of this thesis yielded the following publications:

• T. T.

Part I

Derivative-free mixed binary necklace optimization

Introduction to blackbox optimization methods for mixed integer variables

This part corresponds to the article entitled "Derivative-free mixed binary necklace optimization for cyclic-symmetry optimal design problems", that is accepted for publication in Optimization and Engineering (submitted on In this chapter, we give an overview of derivative-free optimization for mixed integer variables. We consider the following bound constrained optimization problem:

$ & % min x,y f px, yq x P rx, xs Ă R m , y P I n , (1)
where x P R m , y P I n are the continuous and discrete variables, respectively, and I denotes the discrete space (e.g., integer, binary or categorical variables). In this thesis, we focus on solving problem (1) with mixed continuous and binary variables, i.e., I " t0, 1u. For integer or categorical variables problems, encoding can be applied to reformulate the problem so as to deal with continuous and binary variables in many cases. There is a wide study in encoding categorical variables. To give the readers a summary, we make a list of encoding methods in Appendix B.

The objective function f is the output of a blackbox numerical simulator with two features:

1. The derivative, ▽f , is not available, we only have access to function values, and

It is computationally expensive to evaluate the function f .

There is a large body of works in the operation research community that regards Mixed-Integer NonLinear Programming problems (MINLP), see for instance [START_REF] Belotti | Mixed-integer nonlinear optimization[END_REF]. Most deterministic algorithms for solving MINLP are based on branch-and-bound methods. Briefly, the branch-and-bound algorithm is based on recursively sub-dividing the set of possible solutions during the branching step, and estimating bounds on the optimal objective-function value in each branch (the "cut"

or "bound" operation) to find a solution (see e.g., [START_REF] Lodi | Mixed integer nonlinear programming tools: An updated practical overview[END_REF]). In [START_REF] Shin | A penalty approach for nonlinear optimization with discrete design variables[END_REF], the authors do not create a search tree but relax the integrality constraints via a sine function that penalizes the variables for not being integers. Remark that relaxing the binary variables is not possible in our application context. While convex MINLPs can be tackled by several available software, for instance BONMIN [START_REF] Bonami | An algorithmic framework for convex mixed integer nonlinear programs[END_REF] or SCIP [START_REF] Achterberg | Constraint Integer Programming: A New Approach to Integrate CP and MIP[END_REF], nonconvex problems are more difficult and usually require convexification and reformulation strategies. Such strategies are either impossible (reformulation) or still need to be developed when dealing with black-box optimization. This is a real challenge, especially in our context of objective functions that are computationally expensive to evaluate (several hours or even days for a single evaluation). Therefore, we choose in this thesis not to focus on (meta) heuristic methods (e.g., evolutionary algorithms [START_REF] Gendreau | Handbook of Metaheuristics[END_REF] and simulated annealing) due to the large number of objective function evaluations that such approaches require.

Blackbox optimization algorithms only use the values of the objective function at evaluation points to perform the optimization. An overview introduction can be found in [START_REF] Audet | Derivative-Free and Blackbox Optimization[END_REF]. A review and comparison of algorithms and software for derivative-free optimization is presented in [START_REF] Ploskas | Review and comparison of algorithms and software for mixedinteger derivative-free optimization[END_REF]. Generally, derivative-free optimization methods can be classified in two main categories which are direct search methods and model-based methods. Direct Search Methods use systematic rules to explore the domain. Among direct search methods, we can count:

• Pattern search methods [START_REF] Torczon | On the Convergence of Pattern Search Algorithms[END_REF]. They consider the behavior of the objective function on a pattern of points over a rational lattice. A search exploratory strategy moves points in the vicinity of the current iterate followed by a so-called local poll step. A convergence proof for nonsmooth functions approaches can be found in [START_REF] Audet | Analysis of generalized pattern searches[END_REF]. Extension work to mixed variables problems is introduced in [START_REF] Audet | Pattern Search Algorithms for Mixed Variable Programming[END_REF].

• Heuristic methods. They include the popular genetic algorithms (GA) [START_REF] John | Adaptation in natural and artificial systems[END_REF] which are based on the evolution of a population of points. More precisely, by sorting the individuals for reproduction according to an adhoc so called fitness function, the best individuals (with highest fitness values) are selected to generate new candidate points. For example, crossover phase is operated to generate offsprings, i.e., new generations of solutions. Then possibly a mutation step is added to have complete new individuals. The algorithm terminates when the set of offsprings is not significantly different from the parents (no more improvement). This method generally requires a large number of evaluations to find a good solution. Some works extend genetic algorithm to mixed variables problem ( [START_REF] Wu | Genetic algorithms for solving mixed-discrete optimization problems[END_REF] is one example). This extension is usually straightforward by choosing adapted crossover and mutation operators.

• Simplex search methods. These heuristic methods (not to be confused with Dantzig's Simplex method for linear programming) include the famous Nelder-Mead algorithm [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] which is based on the idea of replacing the worst vertex in the set of evaluated points with a new one that gives a better value for the objective function by using geometrical transformations such as expansion, reflection, contraction.

A widely cited DFO algorithm in the family of pattern search methods, NOMAD [START_REF] Abramson | The NOMAD project[END_REF][START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear Optimization with the MADS algorithm[END_REF], implements the Mesh Adaptive Direct Search (MADS) algorithm [START_REF] Audet | Mesh Adaptive Direct Search Algorithms for Constrained Optimization[END_REF] for black-box optimization under general nonlinear constraints. MADS is an extension of Torczon's generalized pattern search algorithms [START_REF] Torczon | On the Convergence of Pattern Search Algorithms[END_REF][START_REF] Audet | Analysis of Generalized Pattern Searches[END_REF]. MADS principally relies on two main steps. The search step is flexible enough to allow local and global explorations with generic strategies such diverse as Latin-Hypercube Sampling (LHS), variable neighborhood search [START_REF] Audet | Nonsmooth optimization through Mesh Adaptive Direct Search and Variable Neighborhood Search[END_REF] or a Nelder Mead method [START_REF] Audet | Mesh-based Nelder-Mead algorithm for inequality constrained optimization[END_REF]. The poll step is critical to the local convergence proof. It involves evaluating the objective function on a discrete grid that is dynamically updated. More recently, [START_REF] Conn | Use of quadratic models with mesh-adaptive direct search for constrained black box optimization[END_REF] introduced a search strategy that automatically constructs quadratic models to try and find promising trial points. An extension of MADS to handle integer variables problems is presented in [START_REF] Audet | The mesh adaptive direct search algorithm for granular and discrete variables[END_REF]. The main idea is to harmonize the minimal granularity of variables with the finest mesh containing all trial points. A different approach for direct search methods for bound constrained problems, named DFLBOX is introduced in [START_REF] Luizzi | Derivative-free methods for bound constrained mixed-integer optimization[END_REF]: instead of searching in all the components of the variable, DFLBOX searches the solution component-wise. The algorithm consists of a continuous and a discrete search procedure that calculates the step sizes. It removes the sufficient decrease condition on the discrete variables and updates the iterate by choosing the coordinate that yields the largest objective function reduction. DFLBOX is extended to a more general constrained problems using the sequential quadratic penalty approach in [START_REF] Liuzzi | Derivative-free methods for mixed-integer constrained optimization problems[END_REF].

Among model-based methods, ones can identify two main categories in literature:

• Derivative-free trust-region methods. These include for instance the DFO algorithm of [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], and NEWUOA [START_REF] Powell | On trust region methods for unconstrained minimization without derivatives[END_REF]. These methods are based on local quadratic models that interpolate the objective function at some evaluated points. The models need to satisfy a so-called fully-linear or fully-quadratic property to ensure the local convergence [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. An exploitation phase minimizes the quadratic approximation within a trust region, while an exploration step explores the domain in view of diversification. An extension of this work was first presented in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF] at the SMAI conference.

• Response surface methods which include Radial Basis Function (RBF) based optimization methods and the Efficient Global Optimization (EGO) method or Bayesian Optimization based on Gaussian Processes (GP). The surrogate models used in these methods are global models, i.e., they use a single substitute of the objective function that aims to be sufficiently predictive in the whole search domain to detect areas of interest with good values of the objective function (exploration), and that can be refined in these areas (exploitation). Note that these exploitation and exploration objectives are similar to the goals of the poll and search steps of MADS. The RBF method for global optimization was introduced by Gutmann [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF], several variations followed [START_REF] Regis | Improved Strategies for Radial basis Function Methods for Global Optimization[END_REF][START_REF] Regis | A Stochastic Radial Basis Function Method for the Global Optimization of Expensive Functions[END_REF][START_REF] Costa | RBFOpt: An open-source library for black-box optimization with costly function evaluations[END_REF], and its extension to mixed discrete variables problems [START_REF] Nannicini | On the implementation of a global optimization method for mixed-variable problems[END_REF][START_REF] Holmström | An adaptive radial basis algorithm (ARBF) for expensive black-box mixed-integer constrained global optimization[END_REF]. EGO [START_REF] Jones | Efficient Global Optimization of Expensive Black-Box Functions[END_REF] is based on a Gaussian process surrogate and an adaptive strategy to propose new evaluation points based on the so-called expected improvement criterion, which balances between exploration and exploitation. These types of methods tackle two main disadvantages, first, the scale of dimension, and the difficulty of giving an explicit convergence proof. In [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF], Gaussian process kernels that are products of continuous and discrete kernels are integrated into an EGO method framework; the resulting mixed categorical (involving integer variables not related to effective quantities) optimization problem is then solved by NOMAD. The strengths and weaknesses of various types of kernels for Gaussian processes are discussed in [START_REF] Pelamatti | Efficient global optimization of constrained mixed variable problems[END_REF].

In recent years, there have been considerable studies dedicated to improving derivative-free optimization algorithms. Some of them focused on overcoming the limitation of Bayesian Optimization for high dimensional problems. The EGO-CMAS algorithm proposed in [START_REF] Mohammadi | Making EGO and CMA-ES Complementary for Global Optimization[END_REF] is one example. In this algorithm, an iteration starts with EGO search, then an interruption is set, to be followed by a CMA-ES search over whole starting points with the initial size (the Covariance Matrix is computed from simulated points). A Random Embedding BO (REMBO) method was proposed in [START_REF] Wang | Bayesian Optimization in a Billion Dimensions via Random Embeddings[END_REF] where instead of optimizing in the original high-dimension space, the function is optimized in a reduced-space variable defined by a random projection based on an assumption of a reduced effective dimensionality of f . In the same year, [START_REF] Kandasamy | High Dimensional Bayesian Optimisation and Bandits via Additive Models[END_REF] proposes the Add-GP-UCB approach which decomposes f as a sum of functions defined in lower-dimensional disjoint spaces and then use GP for each disjoint domain. SEGOKPLS, SEGOKPLS+K [START_REF] Bouhlel | Efficient global optimization for high-dimensional constrained problems by using the Kriging models combined with the partial least squares method[END_REF] use PCA (Principal Component Analysis) to find the most informative components of the variable set, and then apply EGO with an adapted user-defined kernel. A very recent work [START_REF] Tan | CobBO: Coordinate Backoff Bayesian Optimization[END_REF] A first stage involves the combination of RBF and trust-region method, for instance, ORBIT (Optimization by radial basis function interpolation in trust region), [START_REF] Wild | ORBIT: Optimization by Radial Basis Function Interpolation in Trust-Regions[END_REF], and BOOSTERS, [START_REF] Oeuvray | Boosters: A Derivative-Free Algorithm Based on Radial Basis Functions[END_REF].

The idea is to replace the quadratic model with more flexible radial basis functions. The second stage combines EGO and trust-region methods, for instance, TRIKE (Trust-Region Implementation in Kriging-based optimization with Expected improvement) [START_REF] Regis | Trust regions in Kriging-based optimization with expected improvement[END_REF], CYCLONE (CYClic Local search in OptimizatioN using Expected improvement) [START_REF] Regis | Trust regions in Kriging-based optimization with expected improvement[END_REF], TURBO (TrUst-Region BO solver) [START_REF] Eriksson | Scalable Global Optimization via Local Bayesian Optimization[END_REF], TREGO (Trust-Region Framework for Efficient Global Optimization) [START_REF] Diouane | TREGO: a Trust-Region Framework for Efficient Global Optimization[END_REF].

This part gave an overview of derivative-free optimization algorithms, a class of methods built to tackle expensive black-box optimization problems. In the next chapter, we present in more detail derivative-free optimization trust-region methods.

3

Derivative-free trust-region methods for continuous variables problems In this chapter, we address the bounded constrained problem of the form

# min x f pxq x P rx, xs, (2)
where x P R m is a vector of continuous variables. The objective function f is the output of a black-box numerical simulator. We further assume that f is smooth with respect to x and bounded from below to guarantee convergence of the algorithm (see [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF][START_REF] Audet | Derivative-Free and Blackbox Optimization[END_REF]). We also assume that derivatives of f exist but unavailable.

We now begin with the basic background required for the trust-region method framework. We first outline the basic framework of trust-region method, and then describe the derivative-free trust-region methods.

Trust-region framework basics

We begin by providing a review of the trust-region framework for continuous optimization when derivatives of f are available, see [START_REF] Nocedal | Numerical optimization[END_REF] for details. In trust-region methods, at each iteration k, we build a quadratic model r m k around the current iterate x k . Remark that this model is a tangent model with exact gradient at the center of the trust region. This model is assumed to approximate the objective function sufficiently well in a neighbourhood of center x k called the trust region, which is defined based on the center and radius pair px k , ∆ k q, where ∆ k ą 0:

Bpx k , ∆ k q " tx P R m : ∥x ´xk ∥ k ď ∆ k u.
The trust region norm ∥.∥ k can be taken from the standard 2´norm ∥.∥ 2 or the infinite norm. To determine the next iterate x k`1 , we solve a quadratic sub-problem of the form

min sPBp0,∆ k q r m k px k `sq, (3)
where r

m k px k `sq " f px k q `gT k s `1 2 s T H k s, g k is the gradient of f at x k , H k is
r m k px k q ´r m k px k `sk q ě κ d 2 ∥g k ∥ mint ∥g k ∥ ∥H k ∥ , ∆ k u,
where κ d P p0, 1s is a constant. Taylor's theorem guarantees the existence of a trust region that guarantees condition, but it does not give its precise radius. Therefore, updating and adjusting the trust-region radius after each iteration is necessary. Given an approximate solution s k of (3), the trust-region radius is updated depending on the ratio of actual improvement and predicted improvement

ρ k " f px k q ´f px k `sk q r m k px k q ´r m k px k `sk q .
If the model reduction matches well the actual reduction of the objective function (when ρ k ą 0 and near 1), the candidate s k is accepted and the trust region radius is possibly increased. Otherwise, the candidate is rejected and the trust region radius is decreased. We run the loop until convergence criteria on minimal gradient norm and minimal trust-region radius are reached.

Algorithm 1: Derivative based trust region algorithm, see [START_REF] Nocedal | Numerical optimization[END_REF] Input:

x 0 , 0 ă ∆ 0 ď ∆ max , 0 ă γ 0 ă 1 ă γ 1 , 0 ď η 0 ă η 1 ď 1 0. Initialization. Compute f px 0 q, ▽f px 0 q. Set k " 0 1. Model definition. Build the quadratic model r m k in Bpx k , ∆ k q 2. Solve sub-problem (3) s k P argmin sPBpx k ,∆ k q r m k px k `sq.
3. Center update. Compute

ρ k " f px k q ´f px k `sk q r m k px k q ´r m k px k `sk q .
Set:

x k`1 " $ & % x k `sk if ρ k ě η 0 , x k if ρ k ă η 0 .

Trust region update

∆ k`1 " $ ' ' ' & ' ' ' % mintγ 1 ∆ k , ∆ max u if ρ k ě η 1 , ∆ k if η 0 ď ρ k ă η 1 , γ 0 ∆ k if ρ k ă η 0 .
When the derivatives are not available, we build an interpolation model based on a set of given simulations (f values). In this case, the model is considered as a valid approximation of the objective function under some conditions which mainly depend on the geometry of the interpolation set. If the model does not satisfy these conditions, a new point is added to the interpolation set to improve the accuracy of the model. This is detailed in the next section.

2 Derivative-free trust-region methods

The basic idea of DFO trust region methods [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] is to replace the problem, which involves expensive simulations with no information about the derivatives, by a simpler problem (linear or quadratic form) for which we have derivatives. Then, we minimize the simpler problem using the idea of classical trust-region methods. At each iteration, we solve the sub-problem to find the possible candidate for the next iteration. Similarly to derivative-based trust region methods, if the model is qualified as valid (according to the ratio between the actual improvement and the predicted improvement), we decide to increase, decrease or keep the current trust-region radius and choose the solution of the sub-problem as the center of the new trust-region, or not.

We keep running the algorithm until stopping criteria are met, typically when the trust region and the norm of the gradient of the model are small enough.

The main differences between derivative free and derivative based trust region methods are:

• the quadratic models are based on a given interpolation set (available simulations),

• the Taylor expansion error bound is replaced by fully-linear or fully-quadratic model properties, or model qualification condition.

In the next section, we focus on the model construction and the necessary properties of the geometry of the interpolation set to control the error bounds of the models.

Geometry of the interpolation set

Consider a set of sample points given by X " tx 0 , x 1 , . . . , x p u, where

x i P R m , i " 1, 2, . . . , p. (4)
where p ă pm `1qpm `2q{2, we denote by p 1 the number of points of the sample set. Thus, p 1 " p `1.

We would like to build a quadratic model r mpxq which interpolates the points in X. If p 1 " pm `1qpm `2q{2, the problem (2) is a quadratic interpolation problem. If p 1 ą pm `1qpm `2q{2, the problem is over-determined, whereas p 1 ă pm `1qpm `2q{2, it is under-determined. In typical applications, the problem is under-determined since the simulations are computationally expensive.

Let r mpxq denote the polynomial of degree d (e.g., d " 1 with 1 ă p 1 " m `1 or with d " 2, m `1 ă p 1 ď pm `1qpm `2q{2) that interpolates f pxq at the points in X that: r mpx i q " f px i q, i " 0, . . . , p.

(

Let us consider P d m the space of polynomials of degree ď d in R m and a basis ϕ " tϕ 0 , ϕ 1 , ϕ 2 , . . . , ϕ p u. We can express r mpxq as

r mpxq " p ÿ i"0 α i ϕ i pxq, (6)
for some coefficients α i , i " 0, . . . , p,. It is clear that r mpxq is determined if the values of α 0 , . . . , α p are determined. From ( 5) and ( 6), the coefficients α are found by solving the following linear system ¨ϕ0 px 0 q ϕ 1 px 0 q . . . ϕ p px 0 q ϕ 0 px 1 q ϕ 1 px 1 q . . . ϕ p px 1 q . . . . . . . . .

ϕ 0 px p q ϕ 1 px p q . . . ϕ p px p q ‹ ‹ ‹ ‹ ' ¨α0 α 1 . . . α p ‹ ‹ ‹ ‹ ' " ¨f px 0 q f px 1 q . . . f px p q ‹ ‹ ‹ ‹ ' , ( 7
)
or, equivalently, in matrix form:

M pϕ, Xqα " f pXq,
where M pϕ, Xq and f pXq represent respectively the matrix and right-hand side vector of the system 7. If M pϕ, Xq is non-singular (i.e., has full rank) then the system has a unique solution.

In this case, the interpolation set X is said to be poised. Thus, if X is a poised set, then the interpolating polynomial r mpxq exists and is unique.

In [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] the authors prove that if the interpolation set is poised, then the condition number of M pϕ, Xq does not depend on the choice of the basis. The condition number of M pϕ, Xq can be considered as a measure of the poisedness of X, where X is a scaled interpolation set (detail in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF])

X " 1 ∆ rx 1 ´x0 , . . . , x p ´x0 s, with ∆ " ∆pXq " max 1ďiďp ∥x i ´x0 ∥.
In practice, we choose for ϕ the natural basis t1, x 1 , x 2 , . . . , x m , 1 2

x 2 1 , x 1 x 2 , . . . ,

1 pd ´1q! x d´1 m´1 x m , 1 d! x d m u, (8)
Definition 3.1. (Lagrange polynomials) Given a set of interpolation points X " tx 0 , x 1 , . . . , x p u, a basis of p 1 " p`1 polynomials l j pxq, j " 0, . . . , p, in P d m is called a basis of Lagrange polynomials if

l j px i q " $ & % 1 if i " j, 0 if i ‰ j.
If X is poised, then the basis of Lagrange polynomials exists and is unique. Definition 3.2. (Λ´poisedness) Let Λ ą 0 and a set B P R m be given. Let tϕ i pxqu p i"0 be a basis of P d m . A poised set X " tx 0 , x

λ i pxqϕpx i q " ϕpxq, with ∥λpxq∥ 8 ď Λ,
or, equivalently, 3. replacing any points in X by any x P B can increase the volume of the set tϕpx i q, i " 0, . . . , pu at most by a factor Λ, where the volume is defined as

volpϕpXqq " |detpM pϕ, Xqq| p 1 ! .
One note that if a set is Λ 1 -poisedness then it is also Λ 2 -poisedness with Λ 1 ă Λ 2 , but the reverse does not always hold.

Fully-linear and fully-quadratic models

As mentioned before, in trust-region methods, the quality of the model approximation should be controlled just as for Taylor expansion models. To formalize this idea, the class of fully-linear and fully-quadratic models is introduced and detailed in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF].

We suppose that x 0 is given as the initial iterate. We define the level set

Lpx 0 q " tx P R m : f pxq ď f px 0 qu.
We do not only consider the objective function f within the level set, but also extend to the enlarged region

L enl px 0 q " ď xPLpx 0 q
Bpx, ∆ max q.

Assumption 3.1. Suppose x 0 , ∆ max are given. Assume that f is continuously differentiable in an open domain containing the set L enl px 0 q and that ▽f is Lipschitz continuous on L enl px 0 q Definition 3.3. (Class of fully-linear models) Let a function f : R m Ñ R, that satisfies Assumption (3.1), be given. A set of model functions M " t r mpxq : R m Ñ R, r mpxq P C 1 u is called a fully-linear class of models if the following hold:

1. There exist positive constants κ ef , κ eg and ν m 1 such that for any x P Lpx 0 q and ∆ P p0, ∆ max s there exists a model function r mpx `sq P M, with Lipschitz continuous gradient and corresponding Lipschitz constant bounded by ν m 1 , and such that

• the error between the gradient of the model and the gradient of the function satisfies

∥▽f px `sq ´▽ r mpx `sq∥ ď κ eg ∆, @s P Bp0, ∆q, (9)
• the error between the model and the function satisfies

|f px `sq ´r mpx `sq| ď κ ef ∆ 2 , @s P Bp0, ∆q. (10)
Such a model r m is called fully linear on Bp0, ∆q.

2.

For this class M there exists an algorithm, which we will call a "model-improvement" algorithm, that in a finite, uniformly bounded (with respect to x and ∆) number of steps can

• either establish that a given model r m P M is fully linear on Bpx; ∆q

• or find a model m P M that is fully linear on Bpx; ∆q.

One example to illustrate the notion of fully linear model is the first-order Taylor approximation of a continuously differentiable locally Lipchitz function, i.e., f P C 1`.

Assumption 3.2. Suppose x 0 and ∆ max are given. Assume that f is twice continuously differentiable in an open domain containing the set L enl px 0 q and that ▽ 2 f is Lipchitz continuous on L enl px 0 q. Definition 3.4. (Class of fully-quadratic models) Let a function f : R m Ñ R, that satisfies Assumption (3.2) , be given. A set of model functions M " t r mpxq : R m Ñ R, r mpxq P C 2 u is called a fully-quadratic class of models if the following hold:

1. There exist positive constants κ ef , κ eg , κ eh and ν m 2 such that for any x P Lpx 0 q and ∆ P p0, ∆ max s there exists a model function r mpx `sq P M, with Lipschitz continuous Hessian and corresponding Lipschitz constant bounded by ν m 2 , and such that

• the error between the Hessian of the model and the Hessian of the function satisfies

∥▽ 2 f px `sq ´▽2 r mpx `sq∥ ď κ eh ∆, @s P Bp0, ∆q, (11)
• the error between the gradient of the model and the gradient of the function satisfies

∥▽f px `sq ´▽ r mpx `sq∥ ď κ eg ∆ 2 , @s P Bp0, ∆q, (12)
• the error between the model and the function satisfies |f px `sq ´r mpx `sq| ď κ ef ∆ 3 , @s P Bp0, ∆q.

Such a model r m is called fully quadratic on Bpx; ∆q.

2.

For this class M there exists an algorithm, which we will call a "model-improvement" algorithm, that in a finite uniformly bounded (with respect to x and ∆) number of steps can

• either establishes that a given model r m P M is fully quadratic on Bpx; ∆q

• or finds a model m P M that is fully quadratic on Bpx; ∆q.

It is clear that the class of fully quadratic models is better than the class of fully linear models in terms of Taylor-like bounded error, but we will see after that we need more functions evaluations to build a fully quadratic model.

Building the trust-region model

Next, we indicate how to construct a fully-linear (or fully-quadratic) model in the particular context of polynomial interpolation and regression.

Note that we only consider DFO problems whose initial sample set has at least m `1 points.

Case where the number of interpolation points is exactly equal to m `1.

We build a linear interpolation function L X pxq :" α 0 `αT 1 x where pα 0 , α 1 q is the solution of

r1 Xs ˜α0 α 1 ¸" f pXq. (14)
If the interpolation set X is poised, then there exists only one linear interpolation model whose coefficients satisfy [START_REF] Audet | Mesh-based Nelder-Mead algorithm for inequality constrained optimization[END_REF].

Case where the number of interpolation points is larger than m `1 and smaller than pm `1qpm `2q 2 .

To build the model in this case, we can use a least-square regression function L X pxq with the same form as in the case of linear interpolation, but pα 0 , α 1 q is now the solution of

min α 0 ,α 1 « p ÿ i"0 pα 0 `αT 1 x i ´f px i qq 2 ff . ( 15
)
We can also use the minimum Frobenius norm interpolation model M X pxq " α 0 `αx 1 2 x T Hx, H " ph i,j q where pα 0 , α, Hq is the solution of the following quadratic optimization problem

min pα 0 ,α,Hq 1 2 m ÿ i"1 m ÿ j"1 h 2 i,j α 0 `αT x i `1 2 px i q T Hx i " f px i q, i " 0, . . . , p, H " H T . . ( 16
)
Case where the number of interpolation points is exactly pm`1qpm`2q 2 .

We build the quadratic interpolation function M X pxq " α 0 `αT x `1 2 x T Hx, where pα 0 , α T , Hq is the unique solution of

α 0 `αT x i `1 2 px i q T Hx i " f px i q, i " 0, . . . , p. ( 17
)
One notes that the quadratic interpolation is fully linear if X is poised. Comparing with the linear interpolation, the error bounds of the quadratic interpolation are tighter, and the function M X provides both approximations of the gradient and the Hessian. Therefore, in terms of iterations, quadratic interpolation can converge much faster, but it requires more points in the interpolation set.

Case where the number of interpolation points is larger than pm`1qpm`2q In this case, we use regression to build the model as in the case where we have more than m `1 points for linear models.

To guarantee the fully-linear (or fully-quadratic) property of a model, we first have to check that the interpolation set is poised. If the sample set is not poised, we need to improve its poisedness, as explained in the next section.

Model improvement

It is essential to check and improve the model during the optimization iterations. Firstly, we need to "improve " the sample set to ensure its poisedness. In [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], the authors give two main ways to improve the sample set: one based on the Lagrange polynomials and one based on the LU factorization. In our study, we focus on LU factorization to improve the model. The following model improvement algorithm is applied before each model update.

Algorithm 2: Improving poisedness of X via LU factorization, see [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] 0. Initialization Initialize the pivot polynomial basis with some basis, e.g., the monomial basis (8), u i pxq " ϕ i pxq, i " 0, . . . , p. Select pivot threshold ξ ą 0. For i " 0, . . . , p 1. Point selection: If j i is found such that |u i px j i q| ě ξ, swap x j i and x i . Otherwise, recompute x i as

x i P argmax xPB |u i pxq|. Stop if |u i px i q| ă ξ.
2. Gaussian elimination: For j " i `1, . . . , p u j pxq Ð u j pxq ´uj px i q u i px i q u i pxq.

To illustrate the algorithm, we consider an example with a given initial set X 0 and three iterations of the improvement algorithm X 1 , X 2 and X 3 (see Figure 4):

X 0 " ¨0 0 0.1 0 0.5 0 0.6 0.1 0 0.3 0 0.7 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , X 1 " ¨0 0 0.6 0.1 0 0.7 0.5 0 0.1 0 ´0.7 ´0.7 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , X 2 " ¨0 0 ´0.7 ´0.7 0 0.7 0.6 0.1 0.1 0 0.9899 ´0.9899 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' , X 3 " ¨07 0 0.9899 ´0.9899 ´0.7 ´0.7 0 0.7 0.6 0.1 ´1.4 1.4 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . -2 -1 0 1 -1 0 1 2 -2 -1 0 1 -1 0 1 2 -2 -1 0 1 -1 0 1 2 -2 -1 0 1 -1 0 1 2
Fig. 4.: Improvement of the interpolation set via LU factorization: X 0 and X 1 (top); X 2 and X 3 (bottom).

Another way to improve the model is to use Lagrange polynomials. In [START_REF] Powell | On trust region methods for unconstrained minimization without derivatives[END_REF] the author uses the definition of Λ " max 0ďiďp max xPB |l i pxq| and maintains the sample set X by choosing a point entering or leaving the set X so that the value of Λ is reduced. To do that, the farthest point from the center of the trust region is removed, i.e., the point of X associated with the largest value of Lagrange polynomials in absolute value.

Convergence results:

The derivative-free trust-region method is a local method. The local convergence for continuous problems is proven in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. The validity of the convergence of our model is based on the fully-linear or fully-quadratic property in definition (3.3) and (3.4), which is guaranteed by the poisedness of the interpolation set. Besides, the model must satisfy inequality (10.10) in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] @k, r

m k px k q ´r m k px k `sk q ě κ 2 ∥▽ r m k px k q∥ minp ∥▽ r m k px k q∥ ∥H k ∥ , ∆ k q,
for some constant κ P p0, 1q. The first-order critical-point convergence lim kÑ`8

∥▽f px k q∥ " 0 is shown in Theorem (10.13) in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF].

Conclusion

In this chapter, we presented a basic framework of trust-region methods as well as DFO trustregion methods for continuous problems. We emphasized that the local convergence is guaranteed based on the so-called fully-linear and fully-quadratic properties. In the next chapter, we present an extension of DFO trust-region methods from continuous variables problems to mixed continuous and binary problems that was initially proposed in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF]. This chapter discusses the main ingredients of the DFOb-d H algorithm, an extension of DFO trust-region methods to mixed continuous and binary variables proposed in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF]. It will serve as the key building block when proposing our new algorithm in Chapter 5. After summarizing the algorithm, we focus on the theoretical proof of its local convergence, which was not explicitly provided in the original presentation [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF].

The algorithm DFOb-d H aims at improving iteratively a starting feasible solution by solving quadratic optimization "sub-problems" based on quadratic approximation models of the objective function. It starts with a set, Z, of interpolation points at each of which the objective function value is known. Each main iteration involves two main phases: exploitation and exploration.

In the exploitation phase, a quadratic model, m, of the objective function is built with fixed y, then some numerical condition (poisedness) for the interpolation set, Z, is verified, and otherwise the interpolation set is updated. A better current solution is sought by solving trust-region quadratic optimization subproblems yielding updates of Z, and of the radii of the trust regions.

The distance upon which is based the definition of the trust region for the discrete part, t0, 1u n , of the search space is the Hamming distance which can be formulated in a linear constraint:

d H pȳ, ỹq " ÿ j:ỹ j "0 ȳj `ÿ j:ỹ j "1 p1 ´ȳ j q, (18)
for ȳ, ỹ P t0, 1u n . Roughly speaking, this distance simply computes the minimal number of flips (from 0 to 1, or from 1 to 0) required to transform ȳ into ỹ. Then, an exploration phase is added to help the optimization explore wider the binary domain. The convergence result we are about to present in this section is in fact driven totally by the exploitation phase, which solves a continuous quadratic optimization subproblem by temporarily fixing the value of the discrete variables y, and by building fully-linear models of the continuous objective function f p¨, yq.

The outline of the chapter is as follows: Section 1 presents a mixed quadratic model construction based on the mixed interpolation set. Section 2 provides a review of the state-of-the-art initial design methods that are used in DFO. Sections 3 and 4 present the exploitation and exploration steps in DFOb-d H respectively. We end up with the complete DFOb-d H algorithm, some conclusions and discussion for this chapter in Section 4.

The quadratic model

This subsection details how the trust-region quadratic sub-problem model at iteration k, r m k , is built. We consider the vector of optimization variables z " px, yq, where x P R m , y P t0, 1u n . By adding binary variables, the notation of the interpolation set becomes Z " tz 0 , z 1 , . . . , z p u, where z i P R m ˆt0, 1u n .

The classification of problems based on the value of p is unchanged, but with dimension m

`n (for mixed continuous and binary) instead of m (for continuous). We keep all the notations and definitions from Chapter 2, but with mixed binary z " px, yq variables. Suppose that one is given a set of points z i " px i , y i q, x i P R m , y i P t0, 1u n , i " 0, 1, . . . , p, at which the objective function is evaluated with values f i :" f px i , y i q, i " 0, 1, . . . , p, where p ą m `n. The derivative-free trust-region algorithm for mixed binary variables is based on the local quadratic model

r m α,g,H pzq " α `gT z `1 2 z T Hz,
with z " px, yq, x P R m , y P t0, 1u n , and where the coefficients α P R, g P R m`n , and H, a pn `mq ˆpn `mq real symmetric matrix, are solutions of the regularized fitting problem

$ & % min α,g,H"H T 1 2 ∥H∥ 2 F r m α,g,H pz i q " f i , i " 0, 1, . . . , p, (19)
where ∥.∥ F is the Frobenius norm. From a computational perspective, subproblem ( 19) can be addressed by NLP solvers such as IPOPT (see details in [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF]). For the sake of notational simplicity, in the sequel, the model r m α,g,H will simply be denoted by r m.

The interpolation set needs to satisfy some conditions to ensure the uniqueness of the solution of the fitting problem [START_REF] Bremner | Necklaces, Convolutions, and X+Y[END_REF]. Let d " 1 or 2, and let tϕ i u h´1 i"0 be the natural basis of the space of polynomials of degree ď d in R m`n (h is then simply the dimension of this space). In our context where y P t0, 1u n , and when d " 2, the ϕ i elements of this basis are the components of the vector:

ϕpzq " p1, x 1 , . . . , x m , y 1 , . . . , y n , 1 2 x 2 1 , . . . , 1 2
x 2 m , . . . , x i x j , . . . , x i y j , . . . ,

x m y n , y 1 y 2 , . . . , y i y j , . . . , y n´1 y n q,

since the purely quadratic terms in the y i 's are discarded (since y 2 i " y i ). As a consequence, one has h " pm `n `1qpm `n `2q{2 ´n.

The proof of maintaining a Λ´poised sample set in a finite number of steps is evoked in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF]. We give the proof in the following lemma for the case of mixed continuous and binary variables, and d " 1 or 2. Remark that we only consider a proper mixed continuous and binary case, the pure continuous proof can be found in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], and the pure binary case can be asset easily. Lemma 4.1. Let v T ϕpzq be a quadratic polynomial of degree at most d (d " 1 or 2), where ϕpzq is defined above and ∥v∥ 8 " 1. Then, there exists a constant σ 8 ą 0, independent of v, such that

max xPBp0,1q,yPB |v T ϕpzq| ě σ 8 . ( 21
)
For d " 1, σ 8 ą 0 and for d " 2, σ 8 ě 1 4

, where Bpx k , ∆ k q is the ball centered at x k with radius ∆ k , and where B indicates the binary set t0, 1u n .

Proof. We first show that such a constant σ 8 exists. Let us consider

ψpvq " max xPBp0,1q,yPB |v T ϕpzq|. ( 22
)
We see that ψpvq is a norm in the space of vector v. Since the ratio of any two norms in finite dimensional spaces can be uniformly bounded by a constant, there exists a σ 8 ą 0 such that ψpvq ě σ 8 ∥v∥ 8 . The constant σ 8 is defined as:

σ 8 " min ∥v∥8"1 ψpvq. ( 23
)
Thus, if v has l 8 ´norm one then ψpvq ě σ 8 and there exits a continuous vector x P Bp0, 1q and a binary y P B such that |v T ϕpzq| ě σ 8 .

We demonstrate for the two cases: for d " 1, there exists σ 8 ą 0 and d " 2, σ 8 ě 1 4 . For d " 1, the natural polynomial basis of degree one is ϕpzq " r1, x 1 , x 2 , . . . , x m , y 1 , y 2 , . . . , y n s T .

Let u " rv 1 , v 2 , . . . , v m`1 s T and w " rv m`2 , v m`3 , . . . , v m`n`1 s T . Remark that we consider the mixed continuous and binary case, it means ∥u∥ 8 ‰ 0 and ∥v∥ 8 ‰ 0. Therefore, ψpvq becomes max xPBp0,1q,yPB |v 1

`uT x `wT y| which reaches an optimal value when x " u ∥u∥ and when y is chosen as

y i " $ & % 1, if w i ą 0, 0, otherwise,
or when x " ´u ∥u∥ and when

y i " $ & % 1, if w i ă 0, 0, otherwise.
Thus, the optimal value of ψpvq is

|v 1 | `∥u∥ `ř y i "1 |w i | ą 0.
This ends the proof of the first case.

For the case where d " 2: the natural polynomial basis of degree 2 is given as

ϕpzq " t1, x 1 , x 2 , . . . , x m , y 1 , y 2 , . . . , y n , 1 2 x 2 1 , . . . , 1 2 x 2 m , . . . , x i x j , . . . , x i y j , . . . , x m y n u.
Since ∥v∥ 8 " 1, by the definition of l 8 ´norm, there exits i : |v i | " 1. Thus, one of the coefficients of the polynomial qpzq " v T ϕpzq is equal to either 1 or ´1 (corresponding to linear part in the basis 1, x 1 , . . . , x m , y 1 , . . . , y n ) or 1 2 , or ´1 2 (corresponding to the quadratic part in the basis

1 2 x 2 1 , . . . , 1 2
x 2 m , . . . , x i x j , . . . , x i y j , . . . , x m y n ).

Let us consider only the cases where one of the coefficients of qpzq is 1 or 1 2 (the cases ´1 or ´1 2 are analyzed similarly). The largest coefficient in absolute value in v corresponds to a term which is either a constant term, a linear term x i or y i , a quadratic term 1 2 x 2 i or x i x j or x i y j . Let us restrict all variables that do not appear in this term to zero. We will show that the maximum absolute value of qpzq is at least 1 4 by considering six cases of different terms that correspond to the largest coefficient. In each case we will evaluate qpzq at several points in the unit ball and show that at least at one of these points one has |qpzq| ě 1 4 :

• qpzq " 1: it implies directly |qpzq| ě 1 4

• qpzq " x i `1 2 αx 2 i `δ: in this case we have

qp1q " 1 `α 2 `δ, qp´1q " ´1 `α 2 `δ. If α 2 `δ ă 0 then |qp´1q| ą 1, otherwise, if α 2 `δ ě 0 then |qp1q| ě 1.
• , by adding these inequalities, we have

´1 2 ă 1 `2δ ă 1 2 implies that δ ă ´1 4
, and also we have qp0q " δ ă ´1 4

. Thus, we have |qp0q| ą 1 4 .

• qpzq " x i x j `1 2 αx 2 i `1 2 βx 2 j `γx i ``δx j `ϵ: in this case, we consider qpzq at four points

p 1 " p 1 ? 2 , 1 ? 2 q, p 2 " p 1 ? 2 , ´1 ? 2 q, p 3 " p ´1 ? 2 , 1 ? 2 q, p 4 " p ´1 ? 2
, ´1 ? 2 q. We have

qpp 1 q " α `β 4 `1 2 `γ `δ ? 2 `ϵ, qpp 3 q " α `β 4 ´1 2 ´γ ´δ ? 2 `ϵ, qpp 2 q " α `β 4 ´1 2 `γ ´δ ? 2 `ϵ, qpp 4 q " α `β 4 `1 2 ´γ `δ ? 2 `ϵ.
As a result, we obtain qpp 1 q´qpp 2 q " 1`?2δ and qpp 3 q´qpp 4 q " ´1`?2δ. If δ ď 0, we have

qpp 1 q ´qpp 2 q ď 1. If |qpp 1 q| ă 1 2 , then |qpp 2 q| ě 1 2
. Otherwise, δ ą 0, qpp 3 q ´qpp 4 q ě ´1.

Thus, if |qpp 3 q| ă 1 2
, then qpp 4 q ď ´1 2 .

• qpzq " x i y j `1 2 αx 2 i `βx i `γy j `δ: in this case, we evaluate qpzq at six points with the evaluation as follows

qp0, 0q " δ, qp´1, 0q " α 2 ´β `qp0, 0q, qp0, 1q " γ `δ, qp1, 1q " α 2 `p1 `βq `qp0, 1q, qp1, 0q " α 2 `β `qp0, 0q, qp´1, 1q " α 2 ´p1 `βq `qp0, 1q.
From the last two equations, we have qp1, 1q ´qp´1, 1q " 2 `2β.

We prove that with all the possible values of β, we can find at least one point which satisfies the lemma. If β ě 0, we have

qp1, 1q ´qp´1, 1q " 2 `2β ě 2. If |qp´1, 1q| ă 1 then qp1, 1q ą 1, otherwise, if |qp´1, 1q| ą 1 it is trivial. If ´1 2 ă β ă 0, we have qp1, 1q ´qp´1, 1q ą 1, then: if |qp1, 1q| ă 1 2 , we have ´qp´1, 1q ą 1 ´qp1, 1q ą 1 2 , so we obtain |qp´1, 1q| ą 1 2 . Otherwise, if |qp1, 1q| ą 1 2
, the result is trivial.

If β ď ´1 2
, from the third and the fourth equations, we have qp1, 0q ´qp´1, 0q " 2β ď ´1.

In this case, if |qp´1, 0q| ă 1 2 , then qp1, 0q ă ´1 2 , which implies that |qp1, 0q| ą 1 2 . Otherwise, if |qp´1, 0q| ě 1 2
, the result is trivial.

In order to define the poisedness of the interpolation set Z, we need first to define the corresponding pp `1q ˆh interpolation matrix:

M :" ¨ϕ0 pz 0 q ϕ 1 pz 0 q . . . ϕ h pz 0 q ϕ 0 pz 1 q ϕ 1 pz 1 q . . . ϕ h pz 1 q . . . . . . . . .

ϕ 0 pz p q ϕ 1 pz p q . . . ϕ h pz p q ‹ ‹ ‹ ‹ ' . ( 24
)
Let us consider the three possible cases for the dimensions of M (related to the number, p `1, of interpolation points and to the cardinality, h, of the basis -recall that m `n `1 ď p):

• h " p `1 (determined case): Following [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], the interpolation set Z is said to be poised if the determinant of M is non-zero.

• p `1 ă h (under-determined case): Again, as in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], Z is poised if M is full column rank (rankpM q " minpp, hq " p).

• h ă p `1 (over-determined case): In this case we propose to remove p ´h points from the interpolation set (we shall define precisely in the algorithm which points should be eliminated), so that one falls into one of the two previous cases.

A so-called ill-geometry situation leading to a non-poised interpolation set occurs when for instance at some iteration, two or more interpolation points collapse or are affinely dependent. This results in non-uniqueness of solutions of the fitting problem [START_REF] Bremner | Necklaces, Convolutions, and X+Y[END_REF]. There is also an illgeometry problem in the case of a near-singular interpolation matrix (when two interpolation points are too close to each other for example). To prevent this scenario, an improvement step based on LU factorization is set up in the mixed space R m ˆt0, 1u n , inspired from the continuous version of Section 2.2.4 from [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], and detailed in Algorithm 3. It involves solving a mixed-integer quadratic programming problem (MIQP) to be defined below.

Algorithm 3 provides a poised interpolation set such that when Gaussian elimination is applied to the interpolation matrix M , the absolute value of all pivots are not smaller than the chosen threshold ξ.

Algorithm 3: Improving poisedness of Z in the mixed trust region B, extended from [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] 0. Initialization Choose an initial pivot polynomial basis with some basis u i pzq, i " 0, 1, . . . , h, e.g., the monomial basis ϕpzq given by [START_REF] Burer | Non-convex mixed-integer nonlinear programming: A survey[END_REF]. Select a pivot threshold ξ ą 0. For i " 0, 1, . . . , h

Point selection:

• If there exists an index j P ti, i `1, . . . , |Z|u such that |u i pz j q| ě ξ, swap z j and z i in set Z,

• Otherwise, recompute z i as

z i P argmax zPB |u i pzq|, ( 25
)
where B is the trust region we are considering.

Stop if |u i pz i q| ă ξ.

2. Gaussian elimination: For j " i `1, i `2, . . . , p u j pzq Ð u j pzq ´uj pz i q u i pz i q u i pzq.

From a computational perspective, subproblem [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF] can be addressed by a MIQP solver such as CPLEX or BONMIN.

The introduction of binary variables requires an adapted trust-region definition. In the presentation [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF], the authors introduce an ∥.∥ 8 (l 8 -norm) trust region for the continuous variables, and a Hamming-distance trust region for the binary variables.

Assuming in the sequel that the current iterate under consideration is px 0 , y 0 q, the mixed trust region is defined as Bpx 0 , ∆ x q ˆBpy 0 , ∆ y q,

where

Bpx 0 , ∆ x q " tx P R m : ∥x ´x0 ∥ 8 ď ∆ x u, (27)
and Bpy 0 , ∆ y q " ty P t0, 1u n : d H py, y 0 q ď ∆ y u,

for some given trust-region radii ∆ x and ∆ y .

In what follows, we detail the major stages of the DFOb-d H algorithm.

Initial interpolation set

This subsection details the choice of the initial interpolation set (also often referred to as the initial design) in DFOb-d H .

In order to construct a first quadratic model, one requires an interpolation set that contains a sufficient number of points together with the corresponding objective function values. As indicated in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF], for DFOb-d H this number is often taken equal to m `n `1. Further, these points need to satisfy strict geometry conditions for the interpolation problem to be well posed. As remarked in [START_REF] Vu | Surrogate-based methods for black-box optimization[END_REF], a "good" design of experiments (DoE) not only needs to be affinely independent, but should additionally satisfy space-filling, non-collapsing properties [START_REF] Franco | Planification d'expériences numériques en phase exploratoire pour la simulation des phénomènes complexes[END_REF].

There are several methods to choose a given number of sample points in a continuous space, such as factorial designs, Latin Hypercube Sample (LHS), and Optimal LHS designs (see e.g. [START_REF] Vu | Surrogate-based methods for black-box optimization[END_REF]).

However, here we deal with mixed continuous and binary variables problems: we need to provide a DoE in the mixed space R m ˆt0, 1u n . When the dimension is small, one way to proceed is to sample among 2 n`m corner points of the boundary box: for example [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF] proposes a strategy that chooses m `n `1 corner points plus the central point of the box. For larger dimensions, a popular strategy is the Latin Hypercube Sample (LHS) [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code[END_REF][START_REF] Wild | Derivative-Free Optimization Algorithms for Computationally Expensive Functions[END_REF], originally used for generating samples for continuous variables in a bounded subset. However, points sampled by this strategy will surely not satisfy our binary constraints.

In our implementation, we therefore proceed as in [START_REF] Costa | RBFOpt: An open-source library for black-box optimization with costly function evaluations[END_REF] for the RBFOpt algorithm: we first construct a Latin Hypercube Design with maximin distance criterion of m `n `1 points in the considered bounded subset of R m`n , then we round the n components associated with the binary variables to zero or one. Remark that rounding recovers the binary domain but may destroy the desirable properties of LHS or, even worse, it may generate identical points. A future track of our research should therefore be dedicated to improving the method for generating the initial design of experiments.

Two main phases of DFOb-d H remain to be specified: the exploitation and the exploration phases. The exploitation phase attempts at finding locally-optimal solutions of the optimization problem with fixed binary variables. The exploration phase focuses on escaping from local minima when we cannot improve the current local solution, by exploring the binary domain.

Exploitation phase

The next step -that will be denoted Step 1a -involves solving a continuous quadratic-programming (QP) subproblem temporary fixing the binary variables y to the associated current values of the trust region center,

y k : # min x r m k px, y k q s.t. ∥x ´xk ∥ 8 ď ∆ x,k , (29)
where r m k is the current model at the k th iteration, px k , y k q is the current iterate, and ∆ x,k is the trust-region radius with respect to the continuous variables x at iteration k. Note that the infinity norm l 8 is used to define the trust region with respect to continuous variables for the sake of subproblem simplification (leading to bound constraints). In practice, IPOpt is applied to solve this QP subproblem.

The following step, Step 1b, tests whether the solution, x ˚, of (29) should be accepted, based on the ratio, ρ, of the true improvement in f brought by x ˚, over the improvement predicted by the model:

ρ " f px k , y k q ´f px ˚, y k q r m k px k , y k q ´r m k px ˚, y k q , ( 30
)
where one remarks that the denominator is always negative since x ˚is solution of [START_REF] Constantin | On the theory of binary asymmetric error correcting codes[END_REF]. We introduce η good , η ok , and η tol , some pre-defined acceptance threshold values such that η good ą η ok ą η tol ą 0.

If ρ ě η tol , the new iterate is accepted (successful iteration), If ρ ă η tol , the solution is rejected (unsuccessful iteration).

This exploitation phase (referred to as Step 1 in the sequel) is summarized in Algorithm 4.

Algorithm 4: Exploitation phase of DFOb-d H (Step 1) at iteration k

Step 1a (TR QP)

• Solve (29) for fixed y "

y k in B y k px k , ∆ x,k q to get x • Evaluate f px ˚, y k q; if n simu " n simu Ñ STOP • Add ppx ˚, y k q, f px ˚, y k qq to Z k ; p Ð p `1
Step 1b (Validation) Compute the acceptance ratio ρ via (30)

• If ρ ě η tol (successful Step 1):

x k Ð x else (unsuccessful Step 1): x k is rejected.

Exploration phase

After a successful Step 1 with fixed y k , the following step (which will be referred to as Step 1.5a) attempts to improve the current-iterate solution, px k , y k q, in the mixed-variable search space by solving the mixed binary quadratic subproblem:

$ ' ' & ' ' % min xPR m ,yPt0,1u n r m k px, yq s.t. ∥x ´xk ∥ 8 ď ∆ x,k , d H py, y k q ď ∆ y,k . (31)
In practice, this subproblem is addressed by MIQP solvers such as CPLEX or BONMIN.

Then, a validation step (referred to as Step 1.5b) checks if the solution of Step 1.5a provides a solution y ˚‰ y k , and whether the corresponding objective-function value, f px ˚, y ˚q, associated to this solution is smaller than the current best objective-function value.

In case of an unsuccessful Step 1.5a (i.e., no improvement in the minimization of f or failure in solving ( 31)), we continue with the same y k , with a trust-region management with respect to the x component, and a new Step 1 to improve the current solution with y fixed to the value y k :

• If ρ ą η good , then the solution px ˚, y k q is accepted and the model is considered as a "good" predictor of f , the trust-region size is then increased;

• If ρ P rη ok , η good s, then the solution is accepted and the model is considered as sufficiently predictive, the trust-region size remains unchanged;

• If ρ ă η ok , then px ˚, y k q is rejected and the model is not considered sufficiently predictive. The trust-region radius is then reduced.

This trust-region management (which will constitute Step 2) can be summarized as:

∆ x,k`1 " $ ' ' ' & ' ' ' % 2∆ x,k if ρ k ě η good , ∆ x,k if η ok ď ρ k ă η good , 1 2 ∆ x,k if ρ k ă η ok .
If the solution of (31) does not yield improvement with respect to the current center px ˚, y ˚q, and the pre-defined minimal value of the trust-region size, ∆ x , is reached, then px ˚, y ˚q is considered to be a locally-optimal solution. In this case, the algorithm explores in a Step 3 the binary search space using no-good cuts, analogous to those introduced in [START_REF] D'ambrosio | On interval-subgradient and no-good cuts[END_REF] for general mixed optimization problems. The aim of no-good cuts were to remove a point that was integer, but not feasible for all the constraints which we can informally call the "no-good" solutions. This leads in our case to relaxing the trust-region constraint:

d H py, y ˚q ď K,
for some K ą 0, and to force the algorithm to move away from the current locally-optimal solution by adding the extra (no-good cut) constraint:

ÿ j:y j "0 y j `ÿ j:y j "1 p1 ´yj q ě K ˚, (32)
where K ˚P N ˚is some user-defined discrepancy value strictly greater than 1. Note that for a given x k , several such no-good cut constraints are likely to cumulate, as there will be one constraint of the form (32) corresponding to each of the different y ˚values obtained. The set of no-good cut constraints at iteration k is denoted by Ω N GC k . If the exploration does not provide us a better solution, an add-points step is added. In this particular step, we allow to add an user-defined number of points generated by DoEs method. The exploration step (Step 1.5) is summarized in Algorithm 5.

Algorithm 5: Exploration phase of DFOb-d H (Step 1.5) at iteration k

Step 1.5a (MIQP subproblem)

• is_new_N GC " 0

• Solve MIQP (31) in Bpx k , ∆ x,k q ˆpBpy k , ∆ y,k q X Ω N GC k q to get px ˚, y ˚q • Evaluate f px ˚, y ˚q if n simu " n simu Ñ STOP • Add ppx ˚, y ˚q, f px ˚, y ˚qq to Z k ; p Ð p `1
Step 1.5b (Validation)

• If y ˚‰ y k and f px ˚, y ˚q ă min px,yqP|Z k |XpR m ˆΩNGC k q f px, yq (successful step 1.5): ∆ x,k Ð ∆ x,0 , px k , y k q else (unsuccessful step 1.5) If y ˚‰ y k : ∆ y,k Ð ∆ y,k ´1
The algorithm finally ends when the maximal budget of objective-function evaluations or the maximal number of no-good cuts is reached. The maximal number of possible no-good cuts is theoretically equal to 2 n ´1 (for n binary variables). We shall see later that in the context of cyclic symmetry, it is in fact approximately 2 n {n. But more importantly, with this type of property we will see that the definition of no-good cuts in [START_REF] D'ambrosio | On interval-subgradient and no-good cuts[END_REF] is not sufficient to discriminate equivalent configurations.

The complete DFOb-d H algorithm is outlined in Algorithm 6.

Algorithm 6: DFOb-d H algorithm Initialization

• Given initial TR radii 0 ă ∆ x ă ∆ x,0 ă ∆ x , 0 ă ∆ y ă ∆ y,0 ă ∆ y , and tolerances η good ą η ok ą η tol ą 0, a maximal budget of evaluations n simu ą p, maximal number of no-good cut constraints n N GC ą 0, and a corresponding discrepancy value K ˚ą 0.

• Initial interpolation set Z " tpz i , f i qu i"0,1,...,p , where z i " px i , y i q, f i " f pz i q

• Define initial iterate px 0 , y 0 q P argmin i"0,1,...,p f i

• Set k " 0, Ω N GC 0 " t0, 1u n , is_new_N GC " 0 Iteration k:
Step 0 (Model update and improvement)

• Build quadratic model r m k px, yq (cf. Section 1) • Improve poisedness of Z k by solving a TR MIQP in Bpx k , ∆ x,k q ˆBpy k , ∆ y,k q (Algorithm 3) a • If n simu " n simu Ñ STOP if is_new_N GC " 1 : go to Step 1.5
Step Step 2 (TR update and local convergence check)

• If ρ ď η ok : ∆ x,k Ð ∆ x,k {2 • If ρ ě η good : ∆ x,k Ð 2∆ x,k • If ∆ x,k ą ∆ x : k Ð k `1 and go to Step 0
Step 3 (Exploration after local convergence)

• Adding a new no-good cut :

Ω N GC k Ð Ω N GC k X ty P t0, 1u n : d H py, y k q ě K ˚u, n N GC Ð n N GC `1; is_new_N GC " 1 • Reinitialize TR radii: ∆ x,k Ð ∆ x,0 , ∆ y,k Ð ∆ y • If n N GC ď n N GC : k Ð k `1

Local convergence results

In order to avoid ill-conditioning and ensure the local convergence of the algorithm, we rely on a class of so-called fully-linear models [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] within the chosen trust region. As shown in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] for continuous problems, if an interpolation set is poised, then the model obtained by solving the minimal Frobenius fitting problem is fully linear in the trust region of size ∆ x " max i"1,2,...,p p∥x i ´x0 ∥ 8 q defined by the interpolation points, which ensures a control of the model error by controlling the size of the trust region and the interpolation set poisedness with the model improvement step.

In our case, we ensure the local convergence of our algorithm by considering the subproblem with fixed binary variables, and by checking that the model for fixed y " y 0 is fully linear in the trust region: B y 0 px 0 , ∆ x q " tpx, y 0 q : x P R m and ∥x ´x0 ∥ 8 ď ∆ x u.

Note that, in our implementation, the model improvement step (Algorithm 3) is performed in Bpx 0 , ∆ x q ˆt0, 1u n , where Bpx 0 , ∆ x q is defined by [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF]. This allows a larger exploration with respect to binary variables than an improvement step in the mixed trust region Bpx 0 , ∆ x q Bpy 0 , ∆ y q while still fulfilling the required assumptions for Lemma 4.2 below. In the following, we give the proof of fully-linear models for fixed y " y 0 . Assumption 4.1.

• f P C 1`i s a continuously differentiable function with respect to the x variables, and ▽ x f is Lipschitz-continuous in a closed subset, Ω, of the optimization domain, R m ˆt0, 1u n ;

• The interpolation set of p `1 points, Z, is poised in Bpx 0 , ∆ x q ˆt0, 1u n where p ą m `n and

∆ x " max i"1,2,...,p p∥x i ´x0 ∥ 8 q. Assumption 4.2.
At every iteration k of the Algorithm [START_REF] Audet | Derivative-Free and Blackbox Optimization[END_REF], the Frobenius norm of the model Hessian evaluated at iterate px k , y k q, H k , is bounded.

Lemma 4.2. Let px 0 , y 0 q be the initial iterate. Under Assumption 4.1 and Assumption 4.2, the model r mp¨, y 0 q which is constructed from r mpx, yq by fixing y " y 0 is fully linear in B y 0 px 0 , ∆ x q. In other words, for all x P B y 0 px 0 , ∆ x q, there exist κ f , κ g ą 0 such that: |f px, y 0 q ´r mpx, y 0 q| ď κ f ∆ 2

x , [START_REF] Lodi | Mixed integer nonlinear programming tools: An updated practical overview[END_REF] and ∥▽ x f px, y 0 q ´▽x r mpx, y 0 q∥ 2 ď κ g ∆ x .

(

) 34
Proof. The model constructed in mixed space is given as:

r mpzq " c `gT z `1 2 z T Hz,
where z " px, yq, g " pg Thus, the model with y fixed to y 0 is defined as follows:

r mpx, y 0 q " cx `ḡ x x `1 2 x T Hx x, ( 35
)
with cx " ´c `gT y y 0 `1 2 y T 0 H yy y 0 ¯, ḡx x " ´gT x `Hxy y 0 ¯and Hx " H xx .

The gradient of r mpx, y 0 q with respect to x is therefore:

▽ x r mpx, y 0 q " ḡx `H x x.
To be convenient, let us introduce the following notations: f 0 pxq " f px, y 0 q, r m 0 pxq " r mpx, y 0 q, ▽f 0 pxq " ▽ x f px, y 0 q , ▽ r m 0 pxq " ▽ x r mpx, y 0 q and B 0 p∆ x q " B y 0 px 0 , ∆ x q.

We define err f 0 pxq " r m 0 pxq ´f0 pxq, err g 0 pxq " ▽ r m 0 pxq ´▽f 0 pxq.

For all x i P B 0 p∆ x q, we develop px i ´xq T err g 0 pxq " px i ´xq T p Hx x `ḡ x ´▽f 0 pxqq " px i ´xq T Hx x `px i ´xq T ḡx ´f0 px i q `f0 pxq `rf 0 px i q ´f0 pxq ´px i ´xq T ▽f 0 pxqs " m 0 px i q ´m0 pxq ´1 2 px i ´xq T Hx px i ´xq ´f0 px i q `f0 pxq `rf 0 px i q ´f0 pxq ´px i ´xq T ▽f 0 pxqs " err f 0 px i q ´err f 0 pxq ´1 2 px i ´xq T Hx px i ´xq `rf 0 px i q ´f0 pxq ´px i ´xq T ▽f 0 pxqs.

Since f 0 is continuously differentiable, we have:

rf 0 px i q ´f0 pxq ´px i ´xq T ▽f 0 pxqs " ż 1 0 px i ´xq T p▽f 0 px `tpx i ´xqq ´▽f 0 pxqqdt,
which implies that px i ´xq T err g 0 pxq "

ż 1 0 px i ´xq T p▽f 0 px `tpx i ´xqq ´▽f 0 pxqqdt `err f 0 px i q ´err f 0 pxq ´1 2 px i ´xq T Hx px i ´xq. (37)
Using (37) with x " x 0 , we obtain:

px i ´x0 q T err g 0 pxq " px i ´xq T err g 0 pxq ´px 0 ´xq T err g 0 pxq "

ż 1 0 px i ´xq T p▽f 0 px `tpx i ´xqq ´▽f 0 pxqqdt `err f 0 px i q ´1 2 px i ´xq T Hx px i ´xq ´ż 1 0 px 0 ´xq T p▽f 0 px `tpx 0 ´xqq ´▽f pxqqdt ´err f 0 px 0 q `1 2 px 0 ´xq T Hx px 0 ´xq. (38)
First, note that err f 0 px 0 q " 0. Then, for each terms of [START_REF] El Amri | Set inversion under functional uncertainties with joint meta-models[END_REF], we obtain the following upper bounds:

• From the Lipschitz property of f 0 pxq (Assumption 4.1), one has:

ˇˇż 1 0 px i ´xq T p▽f px `tpx i ´xqq ´▽f pxqqdt ˇˇď 1 2 ν∥x i ´x∥ 2 ď 1 2 νp2∆ x q 2 ď 2ν∆ 2 x . (39)
• In the same way, we have:

ˇˇż 1 0 px 0 ´xq T p▽f px `tpx 0 ´xqq ´▽f pxqqdt ˇˇď 1 2 ν∥x 0 ´x∥ 2 ď 1 2 ν∆ 2 x . ( 40
)
• In the following two inequalities, note that ∥ Hx ∥ F is bounded from Assumption 4.2:

ˇˇ1 2 px i ´xq T Hx px i ´xq ˇˇď 1 2 ∥ Hx ∥ F ∥x i ´x∥ 2 ď 1 2 ∥ Hx ∥ F p2∆ x q 2 ď 2∥ Hx ∥ F ∆ 2 x . ( 41
) ˇˇ1 2 px 0 ´xq T Hx px 0 ´xq ˇˇď 1 2 ∥ Hx ∥ F ∥x 0 ´x∥ 2 ď 1 2 ∥ Hx ∥ F ∆ 2 x . (42)
• There exists ϵ 1 ą 0 such that

|err f 0 px i q| ď ϵ 1 ∆ 2 x , ( 43
)
which can be shown by contradiction. Indeed, suppose that we have

|err f 0 px i q| ą ϵ 1 ∆ 2 x @ϵ 1 ą 0. (44)
By definition of err f 0 and from the continuity assumption on f 0 and r m on B 0 p∆ x q, there exist ϵ 1 , ϵ 2 ą 0 such that:

|err f 0 px i q| " |f 0 px i q ´r m 0 py 0 q| " |f 0 px i q ´f0 px 0 q `r m 0 px 0 q ´r m 0 px i q| ď |f 0 px i q ´f0 px 0 q| `| r m 0 px 0 q ´r m 0 px i q| ď pϵ 1 `ϵ2 q∆ x . ( 45
)
Thus, setting 44) contradicts [START_REF] Gilles | Introduction to optimal vector quantization and its applications for numerics[END_REF].

ϵ 1 " ϵ 1 `ϵ2 ∆ x,min ě ϵ 1 `ϵ2 ∆ x in (
Thus, we find from [START_REF] El Amri | Set inversion under functional uncertainties with joint meta-models[END_REF] and the inequalities (39 -43):

|px i ´x0 q T err g 0 pxq| ď 5 2 ∆ 2 x ´ν `∥ Hx ∥ F `ϵ¯, (46)
with ϵ " 2 5 ϵ 1 . Using Cauchy-Schwarz inequality, we obtain:

∥err g 0 pxq∥ 2 ď 5 2 ∆ x ´ν `∥ Hx ∥ F `ϵ¯. ( 47
)
Consider the matrix X "

1 ∆ x " x 1 ´x0 , x 2 ´x0 , . . . , x p ´x0 ‰ .
We recall that the interpolation set Z is defined as

Z " ¨x0 y 0 x 1 y 1 . . . . . . x p y p ‹ ‹ ‹ ‹ ' . ( 48
)
Since Z is poised, Z is full rank, i.e., rankpSq " minpp, m `nq " m `n, based on the fact that p ą m `n (see Section 1), and the m column vectors x 0 , x 1 , . . . , x p are linearly independent. Therefore, X T is a non-singular matrix.

We have

X T err g 0 pxq " 1 ∆ x » - - - px 1 ´x0 q T . . . px p ´x0 q T fi ffi ffi fl err g 0 pxq. (49)
Then, we obtain from inequality (46):

∥X T err g 0 pxq∥ 8 " 1 ∆ x max i"1,...,p |px i ´x0 q T err g 0 pxq| ď 1 ∆ x « 5 2 ∆ 2 x ´ν `∥ Hx ∥ F `ϵ¯ff " ∆ x « 5 2 ´ν `∥ Hx ∥ F `ϵ¯ff . ( 50
)
Moreover, we have:

∥err g 0 pxq∥ 2 " ∥X ´T X T err g 0 pxq∥ 2 ď ∥X ´T ∥ 2 ∥X T err g 0 pxq∥ 2 . ( 51
)
Thus, we obtain:

∥err g 0 pxq∥ 2 ď ? m∥X ´T ∥ 2 ∥X T err g 0 pxq∥ 8 ď ? m∥X ´T ∥ 2 ∆ x « 5 2 ´ν `∥ Hx ∥ F `ϵ¯ff . ( 52
)
Recovering err f 0 pxq by equation ( 37), we have:

|err f 0 pxq| ď ∥err g 0 pxq∥∆ x `2ν∆ 2 x `2∥ Hx ∥ F ∆ 2 x `|err f 0 px i q| ď « ? m∥X ´T ∥ 2 5 2 ´ν `∥ Hx ∥ F `ϵ¯`2 pν `∥ Hx ∥ F q `ϵ1 ff ∆ 2 x . ( 53
)
We complete the proof by the definition of the two required constants:

κ g " 5 2
?

m∥X ´T ∥ 2 ´ν `∥ Hx ∥ F `ϵ¯, (54)
and

κ f " ´ν `∥ Hx ∥ F q ¯´5 2 ? m∥X ´T ∥ 2 `2¯`5 2 ϵ ´?m∥X ´T ∥ 2 `1¯. (55)
We can now state the local convergence of the algorithm.

lim kÑ8 ▽ x f px k , y k´1 q " 0, (56)
or all the limit points of the sequence of iterates in Algorithm 6 are first-order critical points.

The proof is obtained by following the same process as in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] (Theorem 10.13): from the results of Lemma 4.2 above, we can prove the local convergence (convergence for fixed y) of the algorithm with the additional assumption that f is bounded from below for all px, yq P Ω, a closed subset of R m ˆt0, 1u n . REMARK 4.1. As explained in [START_REF] Conn | Global Convergence of General Derivative-Free Trust-Region Algorithms to First-and Second-Order Critical Points[END_REF], an interpolation point outside B y 0 px 0 , ∆ x q has to be replaced in order to ensure a fully-linear model. However, in practice, in order to save expensive objectivefunction evaluations, we allow to go on with a model that is not certified to be fully linear when it yields effective progress in the minimization of the function.

Conclusion

This chapter provides the main ingredient of DFOb-d H , a derivative-free optimization method for mixed continuous and binary variables problems. Unfortunately, based on the fact that the Hamming distance cannot detect the equivalent configurations in the binary part of the variables, the method is not satisfying to address cyclic-symmetry problems. In the next chapter, we give our efforts to adapt DFOb-d H to this type of problem. We begin the chapter by recalling our target application: the optimal design of the blade shape of a turbine for aircraft's engine in order to minimize the vibrations. The aim is to optimize the blade arrangement on the disk of two different pre-defined shapes of blades: a reference shape called A and a mistuning shape B. A binary variable y i is associated with each blade location, taking value 0 for shape A and value 1 for shape B. The optimization solution provides the distribution of the two shapes around the turbine disk.

The reduced order model (ROM) is commonly used for addressing this problem. However, ROM removes a large number of possibilities for optimization, in the sense that ROM considers two reduced subproblems but the cyclic symmetry property remains and is not taken into account explicitly in the standard optimization workflow. For instance, in the case of n " 12 blades, ROM considers two subproblems with 6 binary variables, which ends up with only 28 distinct arrangements (14 for each sub-problem), while the total number of distinct arrangements in the original problem is 352, and the removed configurations could be "useful" candidates for the optimization.

Therefore, we define a new distance (to be used within the DFOb algorithm) to model the cyclic symmetry and avoid the redundancy. This new distance is expected to have a simple form. As noted in the previous chapter, the standard Hamming distance does not appear appropriate for this target.

In this chapter, we focus on defining this new distance, which is inspired by the concept of "necklace" in literature [START_REF] Fredricksen | An algorithm for generating necklaces of beads in two colors[END_REF][START_REF] Gabric | Constructing de Bruijn sequences by concatenating smaller universal cycles[END_REF]. In addition, the adaptation of the optimization method to this new distance is discussed.

Necklace context

To avoid useless costly evaluations of the numerous equivalent solutions for cyclic-symmetry problems, engineers typically resort to simplifications or adapted strategies (such as the reducedorder model methodology [START_REF] Choi | Optimization of intentional mistuning patterns for the reduction of the forced response effects of unintentional mistuning: Formulation and assessment[END_REF][START_REF] Choi | Pattern optimization of intentional blade mistuning for the reduction of the forced response using genetic algorithm[END_REF]) to reduce the optimization problem dimension. However, such simplifications are likely to discard interesting or even optimal configurations.

Concept of "necklace"

In order to avoid re-evaluating costly objective-function evaluations at equivalent blade arrangements, we propose to use the concept of necklace [START_REF] Fredricksen | An algorithm for generating necklaces of beads in two colors[END_REF][START_REF] Gabric | Constructing de Bruijn sequences by concatenating smaller universal cycles[END_REF].

The idea of arranging two different types of blade shapes on the disk is similar to distribute two different colors of beads on a necklace. In the application field of this concept, we identify several promising distances. Our application can be seen as a 2-colors necklace optimization with a fixed number of beads.

The number of distinct arrangements is approximately 2 n n as mentioned in the turbine application [START_REF] Han | Optimization of intentional mistuning patterns for the mitigation of the effects of randon mistuning[END_REF][START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF]. By using the result from [START_REF] Fredricksen | An algorithm for generating necklaces of beads in two colors[END_REF], we obtain that the exact number of necklaces for a given number of beads n is equal to 1 n ř d|n ϕpdq2 n{d , where ϕ is Euler's totient function, i.e., the function that counts the positive integers up to n that are relatively prime to n, and where the summation is taken over all divisors, d, of n.

There are numerous applications of the concept of necklace, many of them are based on necklace distances such as similarity of different types of music [START_REF] Toussaint | A comparison of rhythmic similarity measures[END_REF][START_REF] Toussaint | A mathematical analysis of African, Brazilian, and Cuban clave rhythms[END_REF][START_REF] Toussaint | Computational geometric aspects of rhythm, melody, and voice-leading[END_REF], the dissimilarity of DNA in biology [START_REF] Alsmadi | String matching evaluation methods for DNA comparation[END_REF], calculating the leap years and design calendars [START_REF] Demaine | The distance geometry of music[END_REF], painting car in manufacturing [START_REF] Meunier | Computing solutions of the paintshop-necklace problem[END_REF]. In the the next section, we review some usual necklace distances.

Survey of necklace distances

The study of distances in continuous space is well known and widely studied in the literature. Nevertheless, in the discrete space, there is a limitation of defining a distance: the choice of a distance depends on the problem formulation. Taking into account cyclic symmetry property is a challenge of our problem as it can save a large number of simulations in the optimization procedure. There exists several distances for measuring the difference between two discrete strings. We list the main ones in the following. Definition 5.2. (The Hamming distance) Given two binary strings y " py 1 , . . . , y n q and y1 " py 1 1 , . . . , y 1 n q, the Hamming distance between y and y 1 is given by

d H py, y 1 q " n ÿ i"1 |y i ´y1 i |. ( 57
)
The Hamming distance is easily computed in Opnq operations. The constraint d H py, y ref q ď c, for a given binary string y ref and constant c, is linear.

• A deletion, delpiq, changes y i from 1 to 0 with cost c del ;

• A shift, shpi, jq, changes y i from 1 to 0 and y j from 0 to 1 with cost ∥i ´j∥c sh .

The Hamming distance with shifts between y and y 1 is the minimum cost of a sequence of operations that transform y to y 1 .

The Hamming distance with shifts was mentioned in [START_REF] Minghui | On the sum of distances along a circle[END_REF]. It measures not only the number of mismatches but also how far apart the mismatches occur.

We note that with the values of c ins " c del " 1, c sh " 8, the Hamming distance with shifts reduces to the Hamming distance, and when c ins " c del " 8, c sh " 1, it reduces to the swap distance.

In [START_REF] Minghui | A Linear-Time Algorithm for Hamming Distance with Shifts[END_REF] the author states that the Hamming distance with shifts can be computed in Opnq operations.

Definition 5.6. (The Euclidean interval vector distance) Given y, y 1 two binary sequences of the same length, we represent them as their interval vectors

U 1 " pu 1 1 , . . . , u 1 k q, V 1 " pv 1 1 , . . . , v 1 k q where u 1 i , v 1
i are the number of 0 elements between two 1s values (taking into account the rotations). For instance, X " p1100100000q can be represented as X " p0, 2, 5q. Then, the Euclidean interval vector distance for two binary strings with the same length and the same number of 1s is given by

d E py, y 1 q " g f f e k ÿ i"1 pu 1 i ´v1 i q 2 . ( 60
)
The interval vectors can be computed from the given binary sequences in Opnq operations, and d E py, y 1 q can be computed in Opkq operations, where k is the number of 1s in the sequence. We note that this distance does not take into account the cyclic symmetry property.

Definition 5.7. (The interval-difference vector distance) Given y, y 1 two binary sequences and

U 1 " pu 1 1 , . . . , u 1 k q, V 1 " pv 1 1 , . . . , v 1 k q their interval vectors.
For each sequence, we define new sequences, Z " pz 1 , . . . , z k q, with z i "

u 1 i`1 {u 1 i , i " 1, . . . , k ´1, z k " y k {u 1 1 and Z 1 " pz 1 1 , . . . , z 1 k q, with z 1 i " v 1 i`1 {v 1 i , i " 1, . . . , k ´1, z 1 k " y 1 k {v 1 1 . The interval difference vector distance is d ID py, y 1 q " ´k ÿ i"1 maxpz i , z 1 i q minpz i , z 1 i q ¯´k. (61)
In [START_REF] Toussaint | A comparison of rhythmic similarity measures[END_REF] d ID py, y 1 q is computed in Opnq operations.

Definition 5.8. (Minimum distance of pair assignment (MDPA)) Given y, y 1 two integer sequences, MDPA is defined as

d M DP A py, y 1 q " min Y,Y 1 ´n ÿ i,j"1 |y i ´y1 j | ¯. (62)
The definition can be applied to two binary strings.

There are other distances related to the necklace application field, such as the chronotonic distance or the geometry distance [START_REF] Toussaint | The Geometry of Musical Rhythm[END_REF][START_REF] Toussaint | A comparison of rhythmic similarity measures[END_REF]. We listed above the main formulations to see how the cyclic symmetry could possibly be taken into account. In the next section, we present a new necklace distance tailored to our application.

An adapted distance for cyclic-symmetry problems

To avoid useless costly evaluations of the numerous equivalent solutions for cyclic-symmetry problems (as illustrated in Table 3), scientists and engineers typically resort to simplifications or adapted strategies (such as the reduced-order model methodology [START_REF] Choi | Optimization of intentional mistuning patterns for the reduction of the forced response effects of unintentional mistuning: Formulation and assessment[END_REF][START_REF] Choi | Pattern optimization of intentional blade mistuning for the reduction of the forced response using genetic algorithm[END_REF]) to reduce the optimization problem dimension. However, such simplifications are likely to discard interesting or optimal configurations. This section first defines a new distance to be used in algorithm DFOb so as to avoid to evaluate the costly blackbox objective function at configurations that were previously evaluated (equivalent solutions), without arbitrary removal of (potentially good) candidate configurations. The new distance should lead to constraints easily manageable in efficient optimization methods, just like the Hamming distance which leads to linear constraints (cf. ( 18) and ( 32)). Then, in Subsection 2.2 we propose a reformulation of the algorithm optimization subproblems with this new distance for both the exploitation and the exploration phases. Subsection 2.3 summarizes the adapted algorithm, named DFOb-d neck , and provides the local convergence statement.

The necklace distance

Our blade design application can be seen as a 2-color (or binary) necklace optimization problem involving a fixed number, n, of beads (the number of reference blade shapes). The number of distinct arrangements in our applicative context is therefore given by the number of n-bead necklaces:

1 n ř d|n ϕpdq2 n{d , where ϕ is Euler's totient function (the number of positive integers between 1 and n that are relatively prime to n) and the summation is taken over all divisors d of n.

The new distance we shall use is inspired from the particular l p necklace alignment distance (denoted l p NAD) where p " 1 [START_REF] Bremner | Necklaces, Convolutions, and X+Y[END_REF]. Given two vectors of n real numbers, v " pv 1 , v 2 , . . . , v n q, v 1 " pv 1 1 , v 1 2 , . . . , v 1 n q, v i , v 1 i P r0, 1q, the l p NAD is defined as:

d lpN AD pv, v 1 q " min c,s n ÿ i"1 pd 0 ppv i `cq mod 1, v 1 pi`sq mod n qq p , ( 63
)
where c P r0, 1q is a clockwise rotation angle of the first necklace relative to the second necklace, s P t0, 1, . . . , n ´1u is the perfect matching (best possible shift) between beads, and d 0 is the distance:

d 0 pv i , v 1 j q " mint|v i ´v1 j |, p1 ´|v i ´v1 j |qu
(see [START_REF] Bremner | Necklaces, Convolutions, and X+Y[END_REF] for more detail).

Taking into account the fact that our applications involve uniformly-distributed discrete locations, we set the rotation angle to the constant value c " 0, and we replace d 0 with the simple univariate Euler distance (|v i ´v1 j |). This yields the discrete necklace distance: Definition 5.9. Given two k-ary necklaces of length n: u " pu 1 , u 2 , . . . , u n q and u 1 " pu 1 1 , u 1 2 , . . . , u 1 n q, where u i , u 1 i P ta 1 , a 2 , . . . , a k u, i " 1, 2, . . . , n, the discrete necklace distance between u and u 1 is:

d neck pu, u 1 q " min s"1,2,...,n n ÿ i"1 |u i ´u1 i`s |. (64)
For the purpose of the present study which considers only two possible types of blade design, we focus on the case where k " 2 and the alphabet ta 1 , . . . , a k u reduces to t0, 1u. This leads to the binary necklace distance, denoted d neck , on which our algorithm DFOb-d neck will be based: Definition 5.10. Given y, y 1 P t0, 1u n , the binary necklace distance between y and y 1 is:

d neck py, y 1 q " min i"1,2,...,n d H py, Rot i py 1 qq, (65)
where d H denotes the Hamming distance, and Rot i pyq is the rotation of y by i positions.

For example, considering two binary strings: y " p0, 1, 0q and y 1 " p1, 1, 0q. One has: d H py, y 1 q " 1, d H py, Rot 1 py 1 q " d H py, p0, 1, 1qq " 1, d H py, Rot 2 py 1 qq " d H py, p1, 0, 1qq " 2. Using Definition (65): d neck py, y 1 q " 1.

Metric properties of the necklace distance d neck .

It is clear that d neck is a distance since, for any y, y d neck py, y 1 q ď d neck py, y"q `dneck py", y 1 q.

Proof. Direct consequence of the triangle inequality property of the Hamming distance.

To state the next metric property, we need to introduce the following definition: Definition 5.11. (Cyclic symmetry property in distance) We say that a distance d has the cyclic symmetry property if dpy, y 1 q " 0 ðñ y P Rotpy 1 q, for any given two distinct vectors y, y 1 .

Proposition 5.5. (Cyclic symmetry property) d neck has the cyclic symmetry property. In other words d neck py, y 1 q " 0 ðñ y P Rotpy 1 q, @y, y 1 P t0, 1u n .

Proof. To prove this statement, we prove both implications:

1. Proof of the first implication " ùñ ": "d neck py, y 1 q " 0 ùñ y P Rotpy 1 q Assume d neck py, y 1 q " 0 then min r pd H py, Rot r py 1 qqq " 0, which implies that there exists i P t1, 2, . . . , nu such that: d H py, Rot i pxqq " 0, by the metric property of Hamming distance d H py, Rot i py 1 qq " 0 ðñ y " Rot i py 1 q.

2. Proof of the second implication " ðù ": y P Rotpy 1 q ùñ d neck py, y 1 q " 0.

If y P Rotpy 1 q then there exists i P t1, 2, . . . , nu such that y " Rot i py 1 q, which implies d H py, Rot i py 1 qq " 0. Then, min r d h py, Rot r py 1 qq " 0, hence d neck py, y 1 q " 0.

This invariance property will ensure that equivalent solutions are considered as identical solutions.

Unfortunately, contrary to the Hamming distance (see Equations ( 18) or (32) for instance), a constraint involving the binary necklace distance cannot be straightforwardly expressed as linear constraints (due to the "min" operator involved in the definition of d neck ). The next section proposes a way to address this critical issue for adapting Algorithm 6 to the new distance d neck (which will replace the Hamming distance) so as to obtain an algorithm, that we shall name DFOb-d neck , that deals only with linear constraints.

Reformulation of the QP subproblems involving the necklace distance

The incorporation of the new distance d neck in the QP subproblems involves specific modifications in the formulation of the no-good cuts and of the trust-region constraints.

Necklace-distance based no-good cuts. In order to replace the Hamming distance by the necklace distance in the formulation of no-good cuts, first note that for any real numbers a 1 , a 2 , . . . , a n and any positive integer K ˚, one has:

min i"1,2,...,n ta i u ě K ˚ðñ a i ě K ˚, i " 1, 2, . . . , n. (66)
Now, letting y, y 0 P t0, 1u n and using the above equivalence with a i " d H py, Rot i py 0 qq, i " 1, 2, . . . , n, one straightforwardly obtains:

min i"1,2,...,n td H py, Rot i py 0 qqu ě K ˚ðñ d H py, Rot i py 0 qq ě K ˚, i " 1, 2, . . . , n, (67)
or

d neck py, y 0 q ě K ˚ðñ d H py, Rot i py 0 qq ě K ˚, i " 1, 2, . . . , n. (68)
To summarize, one can formulate a no-good cut that avoids useless costly evaluations by using n linear constraints since [START_REF] Meunier | Computing solutions of the paintshop-necklace problem[END_REF] involves n Hamming-distance inequalities, each of which can be written under the form of a linear inequality following [START_REF] Bouhlel | Efficient global optimization for high-dimensional constrained problems by using the Kriging models combined with the partial least squares method[END_REF].

Necklace-distance based trust regions. The way we replace Hamming distances by binary necklace distances in the exploration phase (more precisely, in the trust-region mixed binary quadratic sub-problem (31) of Step 1.5 of Algorithm 6) is less straightforward.

We consider the mixed binary optimization problem:

$ ' ' ' ' & ' ' ' ' % min x,y r m k px, yq s.t. ∥x ´xk ∥ 8 ď ∆ x,k , d neck py, y k q ď ∆ y,k , y P t0, 1u n , ( 69
)
where r m k : R m ˆt0, 1u n Ñ R is a quadratic function, and x k P R m , y k P t0, 1u n , ∆ x,k , ∆ y,k P R are given. functions satisfying 0 ď g i pzq ď M , for all z P Ω, for some M ą 0. Then, the two following optimization problems are equivalent: Proof. We prove the proposition in two steps:

$ & % min z,
• firstly, we show that, (P 2 ) is a relaxation of (P 1 ) in the sense that if pz, tq is a feasible solution of (P 1 ), then pz, ȳ, tq is a feasible solution of P 2 ;

• secondly we prove that any optimal solution pz ˚, y ˚, t ˚q of (P 2 ) is feasible for (P 1 ).

Let us consider the first assertion: (P 2 ) is a relaxation of (P 1 ).

Let pz, tq be a feasible solution of (P 1 ). Consider now the point pz, ȳ, tq where, for i " 1, 2, . . . , n:

ȳi :"

$ & % 0, if i is the smallest index such that t " min i"1,2,...,n tg i pzqu, 1, otherwise. (72)
Let I be the unique index i such that ȳi " 0.

From the definition of ȳ , we note that:

• t " g I pzq,

• ȳI " 0,

• ȳi " 1 for all i ‰ I.

Then, for i ‰ I, the constraint t ě g i pzq ´M holds since t " g I pzq " min i"1,2,...,n tg i pzqu ě 0 ě g i pzq ´M.

And for i " I, t ě g I pzq ´M holds also since t " g I pzq ě g I pzq ´M ȳI " g I pzq.

Then, pz, ȳ, tq is feasible for (P 2 ).

For the second step, let us now show that: if pz ˚, y ˚, t ˚q is an optimal solution of (P 2 ), then pz ˚, t ˚q is feasible for (P 1 ), i.e., we want to prove that t ˚" min i"1,2,...,n tg i pz ˚qu.

By contradiction, we shall suppose that this optimal solution of (P 2 ) is such that Let I y ˚denote the unique index i such that y i " 0.

t ˚‰ min i"1,
Then, y I y ˚" 0 and y i " 1, for all i ‰ I y ˚.

Using the fact that pz ˚, y ˚, t ˚q is a feasible solution for (P 2 ), we have

t ˚ě g I y ˚pz ˚q. (73)
In the first case, with t ˚ă min

where I ˚" ti : g i pz ˚q " min i"1,2,...,n tg i pz ˚quu.

We have:

• This new solution pz, ȳ, tq is feasible for (P 2 ). Indeed, for i ‰ I ˚, the i th constraint, t ě g i pzq ´M y i , is satisfied for pz, ȳ, tq, since M is an upper bound for the function g i pzq and ȳi " 1.

If i " I ˚, then on the one hand min

r m ϵ px, y 0 q " r mpx, y 0 q `ϵ, (76)
with ϵ "

ϵ 1 ∆ 2 x n ď ϵ 1 ∆ 2
x , where ϵ 1 ą 0 is some small pre-defined value (in our computational results, we choose ϵ 1 " 10 ´8). Lemma 5.1. Under Assumption 4.1 and Assumption 4.2, the perturbed model r m ϵ p¨, y 0 q (defined with fixed y " y 0 ) is fully linear in B y 0 px 0 , ∆ x q. In other words, for all x P B y 0 px 0 , ∆ x q, there exist positive constants κ ϵ f , κ ϵ g such that:

|f px, y 0 q ´r m ϵ px,

y 0 q| ď κ ϵ f ∆ 2 x , ( 77
)
and

∥▽ x f px, y 0 q ´▽x r m ϵ px, y 0 q∥ 2 ď κ ϵ g ∆ x , ( 78
)
with κ ϵ f " κ f `ϵ1 , and κ ϵ g " κ g , and where κ f and κ g are the constants of Lemma 4.2.

The proof is a straightforward adaptation of the proof of Lemma 4.2.

As for the convergence proof (detailed in Lemma 4.2) for the DFOb-d H algorithm, from Lemma 5.1 and with the additional assumption that f is bounded from below, we can prove that the algorithm DFOb-d neck is locally convergent, following the lines of the proof of convergence of the (continuous) DFO algorithm in [START_REF] Conn | Introduction to Derivative-Free Optimization[END_REF] (Theorem 10.13).

Numerical results

This section presents comparative numerical results. After briefly presenting the comparison methodology, we propose in Subsection 3.1 a set of benchmark mixed binary optimization problems that features cyclic symmetry. It consists of a set of 25 instances constructed by transforming existing analytical test problems from the literature, plus one completely original problem related to the design of compressor blades in a helicopter turbomachine. Subsection 3.2 reports numerical results on the 25-instance set, while Subsection 3.3 presents comparative results on the helicopter application problem.

We compare the two versions of our DFOb method (denoted DFOb-d H for the version involving the Hamming distance, and DFOb-d neck for the one with the necklace distance) with three state-of-the-art mixed integer derivative-free optimization methods:

• The mesh adaptive direct search algorithm implemented in NOMAD software [START_REF] Abramson | The NOMAD project[END_REF][START_REF] Digabel | Algorithm 909: NOMAD: Nonlinear Optimization with the MADS algorithm[END_REF],

• The surrogate-based optimization method implemented in RBFOpt [START_REF] Costa | RBFOpt: An open-source library for black-box optimization with costly function evaluations[END_REF] (based on radial basis functions),

• The derivative-free line-search bound constrained method, DFLBOX [START_REF] Luizzi | Derivative-free methods for bound constrained mixed-integer optimization[END_REF].

Following the methodology proposed in [START_REF] Moré | Benchmarking Derivative-Free Optimization Algorithms[END_REF], we compare the solvers' performances in terms of number of evaluations of the objective function. This is a classical indicator for applications involving expensive objective-function evaluations where a solver is evaluated through its capacity to achieve a given function reduction within a limited budget of simulations (evaluations of the objective function). Remark that we did not integrate techniques to overcomes cyclic symmetry property of other DFO methods. One technique can be use to avoid recomputing equivalent solutions is to check if one equivalent necklace is evaluated, if so we neglect this solution and go to other candidates.

In our comparisons we consider that a method solves a problem if it provides a solution x satisfying the following criterion on the objective-function value: f px 0 q ´f pxq ě p1 ´τ qpf px 0 q ´f ˚q, [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] where, in the sequel, f ˚denotes the best function value found by any solver (or the globalminimum value, if known), x 0 is the starting point for each solver (or the best point of the initial interpolation points), and τ is the desired accuracy, a user-defined tolerance value (in our tests, τ " 10 ´3 or 10 ´5). If a solver does not provide a solution that satisfies [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF], we consider that it fails.

Performance and data profiles (see [START_REF] Moré | Benchmarking Derivative-Free Optimization Algorithms[END_REF]) are complementary tools to compare solvers on a collection of problems.

For a set of n p problems P " tp For η ě 1, we define a distribution function ρ s for the performance ratio for a solver s as ρ s pηq " 1 n p cardtp P P, r p,s ď ηu ď 1, for a given solver s, with card denotes the cardinal of a set. This distribution computes the number of problems p that are solved with a performance ratio below a given threshold η. If the value of ρ s pηq is near 1, it means that the solver is good and can solve almost all the problems. The performance profiles is a graph of the function ρ s pηq, η ą 1.

The performance profile of a solver depends therefore on the other solvers tested. For instance, the value of the performance profile of a given solver for a performance ratio of 2 is the number of problems solved by this solver within less than twice the number of evaluations required by the most efficient solver for each problem. However, this does not give an accurate information on the number of evaluations required by a solver to solve a whole collection of problems. This is the reason why data profiles are widely used to compare DFO methods. Data profiles give the fraction of problems that can be solved within a given number of objective-function evaluations. Data profiles therefore provide the performance of the solvers for any given simulation budget.

With expensive optimization problems, we are interested in the performance of solvers as a function of the number of function evaluations which leads to the definition of data profiles. The data profile for a solver s is the fraction of problems that are solved within a fixed simulation budget κ d s pαq "

1

20, otherwise

No-good cut parameter K ˚" 1

New mixed binary optimization test problems featuring cyclic symmetry

This subsection details how we build 25 mixed binary cyclic-symmetry analytical benchmark problems from instances of the literature, and it also presents a simplified real-life application from Safran for designing the compressor blades of a helicopter turbomachine.

The 25 analytical benchmark problems. To our knowledge, no instance of cyclic-symmetry mixed binary optimization problems are proposed in the literature. The last two columns of Table 5

$ ' ' ' ' ' ' & ' ' ' ' ' ' % f 1 pxq, if w " 1, f 2 pxq, if w " 2, . . . f l pxq, if w " l, ( 80
)
for which the integer w is the category variable.

For each of the 15 remaining benchmark problems (the 12 problems from [START_REF] Hock | Test examples for nonlinear programming codes[END_REF][START_REF] Dixon | The global optimization problem: An introduction[END_REF][START_REF] Neumaier | Neumaier's collection of test problems for global optimization[END_REF] and the three problems from [71]), the transformation into a mixed categorical problem goes as follows.

where x end and xend are respectively the lower and upper bounds of the variable x end in the original problem, and where the number, l, of categories is to be set by the user. One thereby obtains a mixed categorical optimization problem involving an objective function of the form F px, wq, where w is a category variable (w P t1, 2, . . . , lu). In our numerical tests, we set the number of categories to l " 4 for these 15 benchmark problems (i.e., other than those that originate from minimax problems, for which l is given by the original instance).

After transforming the 25 problems into mixed categorical problems as above, we then artificially introduce a cyclic symmetry by associating to each value of the categorical variable, w, a necklace (with all the solutions corresponding to its rotations). For instance, if a categorical variable w takes three values in t1, 2, 3u (case where l " 3), the mixed categorical problem of minimizing F px, wq is transformed into a mixed binary problem with cyclic symmetry by considering the following new objective function:

min xPrx,xs,yPt0,1u 2 f px, yq :" $ ' ' ' & ' ' ' % F px, 1q, if y " p0, 0q, F px, 2q, if y " p0, 1q, p1, 0q, F px, 3q, if y " p1, 1q. (81)
The resulting new mixed binary necklace-optimization instances are listed in Table 5 together with the number of continuous variables (m), the number of binary variables (n), the name of the original instance from which it was constructed together with the literature reference. The particular choice of problems is motivated by the sake of comparing methods on a diverse range of problem dimensions and difficulties, including the presence of multiple local minima.

Safran's helicopter application. As mentioned in Chapter 1, the present study is motivated by an application provided by Safran: proposing a design of the turbine blades of a helicopter engine that minimizes the vibrations of the compressor. This application involves m " 1 continuous optimization variable controlling the frequency amplitude, and a vector of n " 12 binary decision variables describing the layout of two reference types of blades on the turbine disk. Safran's engineers provide us a surrogate model built from costly real simulations (several hours of computer time are required for one real single simulation) to allow the computational comparison of the five optimization solvers under study.

Results obtained on the 25 analytical problems

The results obtained with DFOb-d H , DFOb-d neck , NOMAD, RBFOptand DFLBOX over the 25 analytical problems are presented under the form of performance profiles (Figure 5) and data profiles (Figure 6). The required accuracy is set to τ " 10 ´3. The number of objective-function evaluations necessary to reach the best known solution for each problem is also displayed in Figure 7. If the solution is not reached by a solver, no point is displayed but the percentages of successes are indicated for each solver. For all the numerical results, DFOb-d H results are displayed in blue, DFOb-d neck results in red, NOMAD results in black, RBFOpt in magenta, and DFLBOX in cyan.

These three figures show that for this benchmark of 25 analytical problems, DFOb-d neck performs overall the best among the four other optimization methods.

The new method DFOb-d neck succeeds to solve 72% of the 25 problems, whereas RBFOpt succeeds to solve 68%, DFOb-d H 64%, NOMAD and DFLBOX 60%.

The performance profiles of Figure 5 reveal that for a (number-of-simulation) performance ratio equal to 1 (to indicate the percentage of problems for which the solver does as well as the best solver), DFLBOX, DFOb-d neck , and RBFOpt reaches the best value: 36%. In addition, the data profiles of Figure 6 illustrate that DFOb-d neck solves most problems with a relatively small number of objective-function evaluations. For instance, 62% of the problems are solved within a number of simulations less than 15 times the number of variables.

Finally, Figure 7 shows that: DFLBOX requires a small number of simulations (less than 200) for the 60% of the problems it succeeds to solve, the two DFOb methods need less than 200 function evaluations for all successful problems but one (only one successful problem requires 270 function evaluations for each solver), whereas RBFOpt and NOMAD require more than 200 objective-function evaluations for 29.4% and 66.7% of their successful problems, respectively.

To summarize, the results obtained on this benchmark of 25 analytical functions demonstrate the efficiency of the DFOb methods within a limited budget of function evaluations. Moreover, 18 out of the 25 problems are solved by DFOb-d neck with, generally, a smaller number of function evaluations compared with the two state-of-the-art solvers NOMAD and RBFOpt. DFLBOX method shows good performances in terms of number of simulations but is able to solve fewer problems than the DFO methods and RBFOpt.

Design of compressor blades in a helicopter turbomachine

In this subsection, we present results for a simplified optimal design application provided by Safran. Contrary to the benchmark of the previous subsection, this is a real-life cyclic-symmetry application. For the sake of fair comparison, we repeat 50 runs of each of the four solvers. Each run starts with a different design of experiments (of cardinality n `m `1 " 14), following the construction process described in Section 2. Again, we report performance (Figure 11) and data profiles (Figure 12), as well as the number of iterations (Figure 13) to reach, this time, a reduced accuracy of τ " 10 ´5 is mainly considered (due to the scale of the objective function) on the function reduction (79) associated with these 50 runs, we also show the results with the accuracy of τ " 10 ´3 after.

Considering first the results with the accuracy τ " 10 ´5. DFLBOX shows very good performances on this test case: it can solve all of the 50 repetitions within a very small number of simulations. As already noticed on the benchmark function results, DFLBOX is more efficient on small-dimensional problems than on large ones. For this application, which involves only one continuous variable, it clearly outperforms the four other methods. DFLBOX algorithm is based on line searches in each coordinate direction for both continuous and discrete variables.The analysis of the results of this blackbox optimization problem revealed its particular structure: for all feasible combinations of values of the binary variables tested (a few hundreds), the minimum of the function corresponds to the same value of the continuous variable. In other words, this simplified optimal design application problem appears to be separable. Indeed, we observed that for all tested combinations of the n " 12 binary variables, the resulting univariate continuous optimization problem yields the same optimal value for the continuous variables. As a consequence, one could therefore separately minimize the univariate continuous problem for one combination of values of the binary variables, then fix the continuous variables to the corresponding value, and then finally address the resulting combinational combination problem (in 12 binary variables) The relaxation used in DFLBOX works very well on such a problem: one continuous variable and 12 binary variables; this methodology is indeed quite efficient as shown in the presented results.

Concerning the four other methods, the results show the very good performance of DFOb-d neck , as it reaches the value f ˚(up to τ accuracy) for 88% of the runs, compared with 84%, 56% and 78% for DFOb-d H , NOMAD and RBFOpt, respectively. One therefore observes that, depending on the initial design of experiments, some runs do not achieve to reach the required reduction of the objective function. However, the two versions of DFOb are robust methods with regard to the initial design, with more than 80% of success within the budget of simulations. Data profiles in Figure 12 show the robustness of the DFOb methods and their good performances in terms of number of simulations compared to RBFOpt and NOMAD methods. Figure 13 illustrates the efficiency of DFOb-d neck , with only one successful run that requires more than 150 simulations. Finally, Figure 16 displays the distribution of the 50 solutions found by DFOb-d neck for each of the 50 runs. For 88% of the runs, DFOb-d neck converged to a point, denoted px ˚, y 1 q, corresponding to the minimal objective-function value f ˚(probably a global minimum); the remaining runs terminated with three other (locally-optimal) solutions, denoted px ˚, y 2 q, px ˚, y 3 q and px ˚, y 4 q.

All the runs converged to the same value of the continuous-variable component: x ˚" 0.03 (corresponding to some geometry feature of the blades). The discrete-variable components of these four different solutions found by DFOb-d neck are displayed in Figure 17, where the 0 and 1 values correspond to different pre-defined types of blades. The results with the accuracy τ " 10 ´3 are shown in the Figures [START_REF] Bartoli | Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design[END_REF][START_REF] Audet | Mesh-based Nelder-Mead algorithm for inequality constrained optimization[END_REF]. The results are slightly different compare to the above results with the accuracy of 10 ´5: RBFOpt shows a better performance with 90% of successes runs. The results stay the same for four other methods: DFOb-d H , DFOb-d neck , NOMAD, and DFLBOX.

Conclusion

In this chapter we addressed derivative-free mixed binary optimization problems involving a cyclic symmetric property, by proposing an adapted distance, d neck (necklace distance), for the binary search space. We presented theoretical results related to the linear formulation of constraints involving this necklace distance that allowed us to integrate d neck in the trust-region derivative-free method DFOb-d H proposed by [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF] for mixed binary problems, in place of the Hamming distance. The convergence of both DFOb-d H and that of the adapted algorithm, named DFOb-d neck , to a locally-optimal solution was proved.

We proposed 25 analytical mixed binary cyclic-symmetry test problems built from a collection of continuous-optimization instances from the literature. The DFOb-d neck method was evaluated on these analytical instances as well as on a surrogate approximation of a real optimal design application. Three state-of-the-art derivative-free mixed binary optimization solvers, NOMAD, DFLBOX and RBFOpt, were also applied for comparison.

Further analysis on the optimization parameters

This chapter provides an analysis of the impact of the number of points, nDoE, in the interpolation set on the performance of the DFOb-d neck algorithm in Chapter 5. Two options are studied: nDoE " n `m `1 and nDoE " 2pn `mq `1 over the 25 analytical problems detailed in Table 6. The initial designs are generated by the adapted greedy algorithm which will be described in detail in Part II.

We repeat 20 runs for each problem, so in total we have 500 runs. We report the results in performance profiles (Figure 18), data profiles (Figure 19), and the number of iterations (Figure 20). Results for DFO-d neck with nDoE " n `m `1 are displayed in blue, whereas the results for DFO-d neck with nDoE " 2pm `nq `1 are in red.

These three figures show that, for this 25 analytical problems, the results are similar for the two sizes of DoE: nDoE " 2pm`nq`1 and nDoE " m`n`1: DFOb-d neck with nDoE " 2pm`nq`1 succeeds to solve 86.5% of the problems, whereas DFOb-d neck with nDoE " m `n `1 succeeds to solve 86.3%. The performance profiles in Figure 18 show that for a ratio equal to 1, DFOb-d neck with nDoE " m`n`1 reaches the best value: 43.9%, while DFOb-d neck with nDoE " 2pm`nq`1 yields a value of 58.3%. The data profiles of Figure 19 and the number of simulations to reach the best solutions in Figure 20 show very similar results for nDoE " m `n `1 and for nDoE " 2pm `nq `1. The results show that the number of points in the interpolation set has a small impact in the optimization procedure (when this number is relatively small, since we are considering an expensive-simulation context). It is likely that, the distribution of the points itself is more impactful. We suspect that an initial set of well-distributed points, especially for the set of discrete variables, provides useful information to build good approximation models for the optimization and should help to avoid local solutions. Therefore, the second part of the thesis is dedicated to the design of experiments in which we can generate a "good" initial design for the optimization.

Conclusion

Part I presents the main materials for the first subject of the thesis which is dedicated to derivativefree optimization trust-region methods. The work is an extension of [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF], adapted to cyclic symmetry problems. The contribution of the author in this part is, first, the convergence proof of the original work in [START_REF] Conn | A trust region method for solving grey-box mixed integer nonlinear problems with industrial applications[END_REF], a new method, DFOb-d neck , and its convergence proof. This method is an extension of the original method DFOb-d H for cyclic symmetric problems with mixed continuous and binary variables. The introduction of the proposed necklace distance allows to deal with the symmetry and to define new sub-problems for the trust region method so that the problem can be addressed by standard MINLP solvers. The preliminary comparative results obtained on a benchmark of test functions and on the application provided by Safran are very encouraging in terms of both robustness and of simulation costs when compared to state-of-the-art methods. In recent years, the need for efficient design of experiments has emerged as a key research area for analyzing complex physical numerical models. Indeed, such models are usually expensive to evaluate, sometimes several hours or even days are required to run one single simulation, which means that using them to conduct optimization studies or uncertainty quantification and sensitivity-analysis investigations is in general too computationally demanding. In the computer experiments community, the standard practice thus consists in building a surrogate model of the numerical code and use it as a proxy for all further intensive computations required by optimization or uncertainty propagation. However, the final surrogate model accuracy heavily depends on the available samples of the computer code inputs/outputs relationship. For this task, it is critical to build a Design of Experiments (DoE) which provides information about all portions of the experimental region. For example, two appealing concepts are the so-called space filling, [START_REF] Vu | Surrogate-based methods for black-box optimization[END_REF][START_REF] Santner | The design and analysis of computer experiments[END_REF] and the Latin Hypercube Sampling (LHS) [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. If designs of experiments for continuous variables have been extensively studied with several space-filing criteria minimax and maximin [START_REF] Pronzato | Minimax and maximin space-filling designs: some properties and methods for construction[END_REF], discrepancy [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF], maximum projection [START_REF] Joseph | Maximum projection designs for computer experiments[END_REF], the case where models involve both continuous and discrete input variables has been much less investigated and tested.

Design of experiments for continuous variables problems

This chapter focuses on design of experiments in the continuous space D " R m . Remark that as mentioned in [START_REF] Vu | Surrogate-based methods for black-box optimization[END_REF], the two most important requirements for a good DoE are:

• Space-filling property: The design points should be uniformly spread over the entire design space

• Noncollapsing property: Two design points should not share any coordinate values if we not know a priori which dimensions are important This chapter is organized as follows. In Sections 1, 2, 3, and 4 we present the geometrical, marginal, low discrepancy and super-sample from kernel herding designs for continuous variables problems respectively.

Geometrical designs

In this class of DOEs, one uses a deterministic procedure to create initial designs. We highlight some space-filling designs:

Full factorial designs [START_REF] Forrester | Engineering Design via Surrogate Modelling, a Practical Guide[END_REF]. It is the most straightforward method to create initial sets uniformly. The main idea is to divide each dimension i into n i intervals, called cells, of equal length, then select the centers of these cells to constitute the DoE n 1 . . . n m cells. It is clear that the sample set obtained by full factorial designs are highly uniform, but the noncollapsing property is not satisfied. Moreover, a full-factorial design can only be used if the number of sampling points is of the form n 1 . . . n m .

Fractional factorial designs [START_REF] Forrester | Engineering Design via Surrogate Modelling, a Practical Guide[END_REF]. It is a variation of the full factorial designs in which one only uses a fraction of the sample set created by the full factorial designs.

Marginal designs

Latin Hypercube Sample (LHS) designs [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]. LHS divides each axis into n intervals of equal length so that one considers n m cells, then selects n points in each of these n m cells so that there are no two points with the same coordinate in any dimension. To do that, one can for instance simply perturb the coordinates of the points as needed, starting from the beginning.

Optimal LHS designs [START_REF] Johnson | Minimax and maximin distance designs[END_REF]. This type of design aims at improving the space-filling property of LHS by using either the minimax distance:

min |S|"p max xPD dpx, Sq,
where dpx, Sq " min

x 1 PS
dpx, x 1 q, p is the number of required points; or by using the maximin distance: max

|S|"p min x,x 1 PS dpx, x 1 q.
3 Low discrepancy designs

Discrepancy is a concept for measuring the "difference" between uniform point distribution and a given point distribution. For instance, in the case of dimension 1, given the empirical distribution F N of a set X of N points x 1 , x 2 , . . . , x N , and F U a uniform distribution on r0, 1s, we define the discrepancy, D N pXq, by:

D N pXq " sup xPr0,1s
|F N pXq ´FU pXq|.

From the above description, if a sample set has lower discrepancy, it approximates the empirical distribution better. The low discrepancy design can be based on Sobol sequences, Halton sequences, Hammersly sequences or Faure sequences (details in [START_REF] Franco | Planification d'expériences numériques en phase exploratoire pour la simulation des phénomènes complexes[END_REF]). However, in the implementing scope of this study, we focus on adapting the low discrepancy design based on Sobol sequences to a mixed binary and continuous space, which will be presented in the next subsection.

Low discrepancy designs based on Sobol sequences

Low discrepancy designs based on Sobol sequences are defined from primitive polynomials over the finite field Z 2 " t0, 1u (i.e., with coefficients either 0 or 1):

Definition 7.1. A polynomial, pptq, of degree s of the form t s `u1 t s´1 `. . . `us´1 t `us is primitive over Z 2 if it is irreducible and it has order 2 s ´1.

The term of irreducible signifies polynomials which cannot be factored into nontrivial polynomials over the field. The order of polynomial is the smallest positive number i such that pptq divides t i ´1 (or t i `1).

Definition 7.2.

A Sobol sequence S " tx 0 , x 1 , . . . , x N ´1u in dimension 1 is defined as:

x i " 1 2 h p' h i k"1 p k l k q,
where pp 1 , p 2 , . . . , p h q is the binary representative of i, h i "

$ & % 1 `tlog 2 u if i ‰ 0, 1 if i " 0,
and ' refers to the addition in Z 2 , and each l k , k ą s is obtained using the following recurrence relationship:

l k " 2u 1 l k´1 ' 2 2 u 2 l k´2 ' . . . ' 2 s´1 u s´1 l i´s`1 ' p2 s l i´s ' l i´s q,
where u k are the coefficients of a primitive polynomial degree s over Z 2 , and l 1 , . . . , l s must be odd integers such that 1 ď l k ď 2 k , k " 1, 2, . . . , s. In order to construct a Sobol sequence in dimension d, we choose d distinct primitive polynomials.

A low discrepancy design based on Sobol sequences is rapid to build given the binary's representative. However, due to the recurrence relation, the new points constructed tend to be in the vicinity of the previous iteration.

Super-sample from kernel herding

Another type of method in this class DoEs is based on Maximum Mean Discrepancy (MMD) designs, called kernel herding [START_REF] Chen | Super-Samples from Kernel Herding[END_REF]. MMD (or kernel distance) between two probability distributions P X 1 and P X 2 is defined as follows:

M M D 2 pP X 1 , P X 2 q " ∥µ P X 1 ´µP X 2 ∥ 2 H k , ( 82
)
where H k is a Reproducing Kernel Hilbert Space (RKHS) with kernel k and µ P " ş kpx, .qdPpxq is the representative of the probability distribution P. Recall that, a RKHS is a Hilbert space in real-value functions in a closed bounded subset X of R m with the property that for each t P X the evaluation functional is a bounded linear functional.

We start with two important definitions to follow up this part: Definition 7.3. (Positive definite kernel, [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]) Consider the domain I. A function k : I ˆI Ñ R, is a positive definite (p.d) if it satisfies the two following criterion 1. symmetric: for any pair of inputs x i , x j P I, k satisfies kpx i , x j q " kpx j , x i q, 2. positive definite: for any p inputs x 1 , x 2 , . . . , x p P I and any p constants c 1 , c 2 , . . . , c p P R, k satisfies

p ÿ i"1 p ÿ j"1 c i c j kpx i , x j q ě 0.
Definition 7.4. (Stationary kernel (or translation invariant kernel), [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]) A stationary kernel k can be expressed in terms of a multivariate function k : R m Ñ R of a single m´dimensional variable, i.e., kpx 1 , x 2 q " kpx 1 ´x2 q, x 1 , x 2 P R m .

A standard example is given by the anisotropic Gaussian kernel:

kpx 1 , x 2 q " m ź l"1 e ´λ2 l px 1 l ´x2 l q 2 .
If X 1 " P X 1 and X 2 " P X 2 we can expand M M D 2 to get an expression involving kp¨, ¨q explicitly:

M M D 2 pP X 1 , P X 2 q " E X 1 ,X 1 pkpX 1 , X 1 qq `EX 2 ,X 2 pkpX 2 , X 2 qq ´2E X 1 ,X 2 pkpX 1 , X 2 qq, (83)
where X 1 " P X 1 and X 2 " P X 2 are independent copies of X 1 and X 2 respectively.

To build DoEs based on MMD, suppose that we have a target probability distribution P X given as a sample, i.e.,

P X " 1 N N ÿ i"1 δ x i , (84)
where px i q i"1,2,...,N is the given sample of size N , δ x is the mass function at x. Typically N will be large.

We want to find a probability distribution P X 1 which approximates P X , in the sense of the MMD, with the following form:

P X 1 " 1 n DoE n DoE ÿ j"1 δ x 1 j , x 1 j P X 1 , (85)
where n DoE is the number of required sample points and is strictly smaller than N . One can choose the x 1 j , j " 1, 2, . . . , nDoE so that P X 1 represents the target probability distribution P X .

To find px 1 j q j"1,2,...,n DOE we have to solve the following optimization problem:

min x 1 1 ,x 1 2 ,...,x n 1 DOE M M D 2 p 1 n DOE n DOE ÿ j"1 δ x 1 j , 1 N N ÿ i"1 δ x i q. ( 86
)
If nDoE is large, the optimization problem might be difficult to solve, since it is a nonlinear and non-convex problem. A possible work-around is to use a greedy approach where we optimize x 1 1 first, then x 1 2 , and so on. In kernel-herding [START_REF] Chen | Super-Samples from Kernel Herding[END_REF], a greedy sequential algorithm is proposed, which consists of the following steps:

x 1 " argmax xPU 1 N N ÿ j"1 kpx, u j q, x t`1 " argmin xPU 1 t `1 t ÿ i"1 kpx, x t q ´1 N N ÿ j"1 kpx, u j q.
From a computational perspective, in this case we only have to compute kp¨, ¨q for the samples px j q, so we can precompute

1 N N ř j"1 kpx 1 , x j q for all x 1 P tx 1 , x 2 , . . . , x N u.
Regarding the choice of the kernel, we assume that px j q are continuous variables that lie in R m ; we have to specify a kernel function kp., .q; k :R m ˆRm ÝÑ R px, zq ÝÑ kpx, zq

Starting from a kernel k : R ˆR ÝÑ R defined in dimension 1, we can typically define the kernel in dimension m by two ways; either "isotropic": kpx, x 1 q " kp∥x ´x1 ∥q, or by "tensorization":

kpx, x 1 q " m ź l"1 kpx l , x 1l q.
We use tensorization to construct the multidimentional kernel in our case.

Besides the kernel-herding approach that tries to minimize MMD, the Kernel Stein Discrepancy (KSD) [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF] is also a tool to measure MMD. The main difference between these two approaches is that the greedy algorithm (see Algorithm 7 above) requires that the target distribution µ is Algorithm 7: Greedy algorithm for continuous-variable problems Input: px j q j"1,2,...,N , n DoE , and a p.d kernel k

• Precompute 1 N N ÿ j"1
kpx 1 , x j q, @x 1 P px j q.

• Compute

x 11 " argmax x 1 PX 1 1 N N ÿ j"1 kpx 1 , x j q. (88)
• For t " 1, 2, . . . , n DoE :

x 1t `1 " argmin x 1 Ptx 1 ,x 2 ,...,x N u 1 t `1 t ÿ i"1 kpx 1 , x 1t q ´1 N N ÿ j"1 kpx 1 , x j q. (89)
Return x 11 , x 12 , . . . , x 1n DoE either explicit or can be easily approximated by sampling, whereas KSD can be used in the case where µ is not explicit, such as when it arises as an intractable posterior in a Bayesian context. Thus, the two methods have the same goal, but a different way of approach. As mentioned in the introduction section, there exists an extended KSD measure approach for discrete space with a new choice of "discrete" Stein operation [START_REF] Yang | Goodness-of-Fit Testing for Discrete Distributions via Stein Discrepancy[END_REF], which is considered as a state-of-the-art method. However, in our case, when the experimental region is given, it is possible to generate a target distribution by uniformly sampling points over the domains. That is why we decided to use kernel herding in this context.

Remark that to apply the kernel-herding algorithm, the kernel to be used must necessarily be positive definite.

Besides, we can count random design which is one of the simplest way to sample points. It generates points randomly in the design space. This type of design has low uniformity, i.e., design points are not uniformly spread over the domain, and it is possible to have duplicated points.

There are model-based DOEs which assume the form of the function to evaluate based on the correlation definition, also known as statistical designs. Some popular methods in this class can be counted such as D-optimality, I-optimality and G-optimality designs, [START_REF] Vu | Surrogate-based methods for black-box optimization[END_REF]. Nevertheless, in this framework, we do not concentrate on this class of method so we only list them here for the readers to have a larger view of DOEs.

Introduction to mixed continuous and discrete variables DoEs

For handling mixed DoEs, the first proposed approaches mainly consisted on simple extensions of continuous LHS, either by randomly discretizing continuous values [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] or by rounding continuous DoEs to obtain feasible integer candidates [START_REF] Costa | RBFOpt: An open-source library for black-box optimization with costly function evaluations[END_REF][START_REF] Müller | SO-MI: A surrogate model algorithm for computationally expensive nonlinear mixed-integer black-box global optimization problems[END_REF]. These techniques can recover the integer domain but may destroy the desirable properties of the original DoEs, or even worse, they may generate identical points which must be avoided. On the other hand, it is possible to generate independently a LHS for every given discrete possible values [START_REF] Santner | The design and analysis of computer experiments[END_REF] or use a so-called Sliced LHS, with LHS for the continuous factor which is sliced into groups of smaller LHS designs associated with different discrete levels [START_REF] Qian | Sliced Latin Hypercube Designs[END_REF]. Although popular, these last two approaches usually require a very large number of samples when the number of discrete variables increases. A recent work [START_REF] Pelamatti | Bayesian optimization of variable-size design space problems[END_REF] proposes to sample continuous variables through a single continuous LHS, while the discrete variables are obtained by randomly assigning an equal number of data samples to each discrete value. A different line of work mainly studied in the machine learning research area relies on the approximation of probability measures with an empirical probability distribution supported by a small number of points, also called quantization [START_REF] Graf | Foundations of quantization for probability distributions[END_REF]. Interestingly, this point of view will be the main ingredient of our proposal.

More precisely, we will focus on the framework of kernel-embedding of probability distributions [START_REF] Smola | A Hilbert space embedding for distributions[END_REF], which offers mathematical tools to define distances between probability distributions which can be very efficiently computed. It is then straightforward to recast the problem of finding a DoE with good space-filling properties as an optimization problem, where one wants to find the empirical probability distribution (the DoE) which is as close as possible to a target probability measure (for example the uniform measure on the feasible experimental region). For continuous probability measures, this point of view has already been proposed with the name of kernel herding [START_REF] Chen | Super-Samples from Kernel Herding[END_REF] or support points [START_REF] Mak | Support points[END_REF]. Note that the target distribution here can either be given explicitly or as an empirical distribution itself but with a large number of points. Here, our goal is to extend this framework to the mixed case.

This chapter is organized as follows. First, we present in Section 2 a naive extension of continuous DoEs which will serve as a baseline for numerical comparisons. In Sections 3 and 4 two novel methods based on kernel-embedding are introduced and discussed with their respective advantages and imitations. Extensive numerical experiments are finally conducted in Section 6, where we show that on a variety of problems the kernel point of view clearly surpasses standard approaches.

Straightforward transformation of a continuous DoE into a mixed DoE

Recall that our experimental mixed region is

D " tz " px, yq P R m ˆIn u. (90)
where x P R m , y P I n are the continuous and discrete variables, respectively, and I denotes the discrete space. Our goal is to select n DoE representative points from D. Inspired from previous works [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF], [START_REF] Tran | Derivative-free mixed binary necklace optimization for cyclic-symmetry optimal design problems[END_REF]), we first propose a straightforward algorithm which is based on the projection from continuous space to integers representing the distinct indices of discrete variables. We first generate a DoE for the discrete variables as detailed in the pseudo-code given in Algorithm 8. if u i P ppj ´1qw, jws, i " 1, . . . , n DoE , j " 1, . . . , n levels then assign the j´th row of L to y i Return y 1 , . .

. , y n DoE

The choice of the method to generate the continuous DoE is left to the user, but for the sake of having a uniform distribution in the discrete space, it is recommended to either use a low discrepancy method such as Sobol sequences, Halton sequences or others (see details in [START_REF] Franco | Planification d'expériences numériques en phase exploratoire pour la simulation des phénomènes complexes[END_REF]), or standard space-filling DoEs such as the minimax or maximin ones for example [START_REF] Mckay | A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF][START_REF] Pronzato | Minimax and maximin space-filling designs: some properties and methods for construction[END_REF].

Remark that this adaptation focuses on picking distinct discrete candidates. If it yields badlydistributed points (this will happen if the continuous DoE is of poor quality), then the algorithm can post process duplicated discrete points (which need to be avoided). Algorithm 8 only provides the DoE for discrete variables: in order to obtain the full DoE in mixed space, the continuous parts are generated separately by another continuous DoE and added to the output of Algorithm 8. In the numerical comparisons of Section 5, we shall use a Sobol sequence and a classical LHS for both continuous DoEs.

3 Greedy-MDS: kernel herding and multi-dimensional scaling

In this section, we propose a first extension of the previous greedy algorithm to mixed continuous and discrete space by using Multi-Dimensional Scaling (MDS). The main idea is to apply MDS to build a continuous encoding of the discrete variables.

Continuous encoding of discrete variables is a standard way to handle mixed variables in regression or classification problems. Successful strategies rely on a data-driven approach, where the encoding is optimized during the supervised learning task [START_REF] Guo | Entity embeddings of categorical variables[END_REF][START_REF] Zhang | A latent variable approach to Gaussian process modeling with qualitative and quantitative factors[END_REF]. However in our unsupervised setting, such approaches cannot be envisioned. In order to circumvent this limitation, we take a different road by assuming that an user-defined distance in the discrete space characterizing any prior information is available. Such a distance can include symmetry properties such as in the numerical experiments presented below, or hand-crafted correlation given by experts of the phenomenon under study. Once this distance is available, it can be used to compute all the pairwise distances between the points from the target distribution u i .

MDS can then be applied to this distance matrix. Indeed, MDS is a technique that creates a map displaying the relative positions of a number of objects, given only a table of distances or dissimilarities between them. In another words, MDS translates the information about the pairwise distances among a set of N objects (normally in a high dimensional space) into a configuration of N points mapped into a smaller dimensional space preserving the distances as much as possible. Suppose that we have a set of N points tu 1 , u 2 , . . . , u N u for which we know all the pairwise distances d ij " dpu i , u j q. The output of MDS is a new set of points v 1 , v 2 , . . . , v N in a Euclidean space such that their Euclidean distance approximates d ij , i.e., d ij « dij " }v i ´vj }.

The dimension of this Euclidean space is chosen by the user, with a maximum equal to N ´1.

To measure how well the original set of data are represented by the MDS model, it is typical to use the goodness-of-fit values, or stress values [START_REF]Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF], based on the differences between the actual distances and their predicted values:

stress " d ř pd ij ´d ij q 2 ř d 2 ij . ( 91
)
Based on the stress values, [START_REF]Multidimensional scaling by optimizing goodness of fit to a nonmetric hypothesis[END_REF] asserts how the MDS model performs: if stress " 0 the MDS model is said to be perfect, if stress " 0.025 the MDS model is said to be excellent, if stress " 0.05 the MDS model is said to be good, if stress " 0.1 the MDS model is said to be fair, and finally if stress " 0.2 the MDS model is said to be poor.

Our proposal, called greedy-MDS, is then to apply a preliminary MDS step for all the discrete variables Y " ty 1 , y 2 , . . . , y N u in the large target sample for which we have all the pairwise distances defined by the user. After this step, we have a continuous encoding V " tv 1 , v 2 , . . . , v N u of the discrete variables, which can be stacked with the continuous variables X " tx 1 , x 2 , . . . , x N u. This gives a new representation of the large target sample with only continuous values. Then, we use kernel-herding from Algorithm 7 with a standard continuous kernel to obtain a small subset of points which represents as well as possible the target distribution. The last step is finally to invert the encoding: this is an easy task, since by bijection we exactly know the correspondence between the original discrete values and their encoding. A pseudo-code is given in Algorithm 9 below.

Algorithm 9: Greedy-MDS algorithm. Input: nDoE, d, Y

1. Pre-processing step

• Compute d, the matrix of distances d ij " dpy i , y j q for all y i , y j P Y .

• MDS step: apply MDS to d, obtain accordingly continuous points, u 1 , . . . , u N

• Stacking with continuous variables: for each u i , i " 1, . . . , N , add a large number of continuous values, save in matrix R 2. Greedy algorithm Algorithm 7 with inputs R, nDoE.

Retrieve the original values by inverting MDS

A strong potential limitation of Greedy-MDS lies in the choice of the user-defined distance. Indeed, very poor distance reconstruction performance can occur in practice, meaning that MDS cannot find an encoding which is able to well preserve the pairwise distances even with a maximum dimension of N ´1. To circumvent this issue, we observed numerically that in practice MDS has a much better reconstruction performance if instead of relying on a user-defined distance, we employ a user-defined kernel k for discrete variables and build the corresponding kernel-induced pseudo-distance d k pu, vq " a kpu, uq `kpv, vq ´2kpu, vq.

Note however that requiring a user-defined kernel is much more demanding than just a distance, and that if a kernel is available the following approach should be preferred.

4 Adapted-Greedy: directly using a kernel defined in the mixed space

In this section, we propose an adapted greedy algorithm for mixed continuous and discrete variables problems that takes into account directly an appropriate mixed kernel inside the greedy algorithm 7, appropriate being related here to the distance and type of information that are given by the user. The first step is to build a positive definite kernel on the discrete variables from a given distance dp¨, ¨q. The naive approach consists in generalizing the Laplace kernel as follows:

kpy, y 1 q " e ´γdpy,y 1 q for any values y and y 1 of the discrete variables. Unfortunately this does not provide a positive definite kernel in general [START_REF] Van Den | Harmonic analysis on semigroups: theory of positive definite and related functions[END_REF]. A particular case of interest where this result holds concerns binary variables with the Hamming distance, d H py, y 1 q, where the kernel k H py, y 1 q " e ´γd H py,y 1 q (92)

is positive definite for any γ ą 0, see [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]. For the general case, we have to follow a different road and build upon the work of [START_REF] Wu | Efficient Global String Kernel with Random Features: Beyond Counting Substructures[END_REF] where they introduce the so-called soft string kernel. For any distance d between any variables u and v in a space I, it is given by k sof t pu, vq :" ż ωPI ppωqe ´γpdpu,ωq`dpv,ωqq dω [START_REF] Shin | A penalty approach for nonlinear optimization with discrete design variables[END_REF] where ppωq : I Ñ R is a probability distribution over I, and it can be proved that it is positive definite [START_REF] Wu | Efficient Global String Kernel with Random Features: Beyond Counting Substructures[END_REF]. When I is a discrete space it boils down to the discrete sum k sof t pu, vq :" ÿ ωPI ppωqe ´γpdpu,ωq`dpv,ωqq . (

Interestingly, for a uniform distribution, this kernel is an approximation of the naive proposal above, by using well-known approximation results for the LogSumExp function (LSE). More precisely, the LSE function reads LSEpx 1 , x 2 , . . . , x n q " logpe x 1 `ex 2 `. . . `exn q and is an approximation to the maximum maxpx 1 , . . . , x n q [START_REF] Lange | Applications of lp-Norms and their Smooth Approximations for Gradient Based Learning Vector Quantization[END_REF]. Consequently, we also have minpx 1 , . . . , x n q « ´logpe ´x1 `e´x 2 `. . . `e´xn q.

Now if p is a uniform distribution on I, meaning that @ω P I, ppwq " p 0 the soft string kernel simplifies as follows k sof t pu, vq "

ÿ ω i PI
p 0 e ´γpdpu,ω i q`dpv,ω i qq , "

ÿ ω i PI e p´γη i `logpp 0 qq , (95)
with η i " dpu, ω i q `dpv, ω i q and applying ´log to both sides leads to ´logpk sof t pu, vqq " ´log ˜ÿ ω i PI e p´γη i `logpp 0 qq « min i pγη i ´logpp 0 qq using the approximation of the minimum given above. From the distance triangle inequality we also have η i ě dpu, vq, which means that min i pγη i ´logpp 0 qq " γdpu, vq ´logpp 0 q. This finally gives k sof t pu, vq " e logpk sof t pu,vqq « e ´pγdpu,vq´logpp 0 qq " p 0 e ´γdpu,vq .

This result can help interpret k sof t as a positive-definite approximation of the Laplace kernel where we would inset the distance d.

Combining one of the above kernels [START_REF] Santner | The design and analysis of computer experiments[END_REF] or [START_REF] Smola | A Hilbert space embedding for distributions[END_REF] for the discrete variables and any kernel k cont for the continuous variables is finally straightforward with tensorization, and leads to a mixed kernel that can directly be used in kernel herding.

Extension beyond mixed variables problems.

Interestingly, the Adapted-Greedy approach presented before is generic, in the sense that we discussed building a kernel adapted for discrete variables, but the same principle can be applied to any other type of variables provided the user can provide a kernel or a distance. A prominent example in computer experiments involves time series variables, for which it is usually a challenge to build a space-filling DoEs. Some previous works are for example based on decomposition of functional data on orthogonal bases such as in [START_REF] Nanty | Sampling, metamodeling, and sensitivity analysis of numerical simulators with functional stochastic inputs[END_REF], but here we can directly make use of the literature on kernels for time-series data. We can for example propose a normalized kernel between two time series uptq and vptq involving the recently proposed global alignment kernel [START_REF] Cuturi | A Kernel for Time Series Based on Global Alignments[END_REF]:

K GAK puptq, vptqq " Kpuptq, vptqq ´1 2 pKpuptq, uptqq `Kpvptq, vptqqq, ( 96
)
where K is the global alignment kernel for same length time series Kpuptq, vptqq "

| ź i"1 e ´ϕσpuptiq,vptiqq , (97)
with ϕ σ pupt i q, vpt i qq " 1 2σ 2 ∥upt i q ´vpt i q∥ 2 `logp2 ´e´∥ upt i q´vpt i q∥ 2 2σ 2

q. We will see in Section 5 how it performs numerically on a time series example.

Numerical experiments

The aim of this section is to demonstrate the large range of applicability of the proposed design of experiment methods. The first application is the surrogate modelling of the simulated performances of an electric engine with respect to design parameters of the rotor component. This application deals with continuous, integer and categorical variables. Then, we consider the case of cyclic-symmetry problems, with continuous and binary variables, motivated by an application of optimal design of the turbomachine of an helicopter engine [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF]. For these specific problems, we address two different operational objectives: (a) surrogate modelling illustrated on various benchmark functions and (b) the choice of initial points for the shape optimization problem of the turbomachine. Finally, our DoE methods are applied to two examples involving time series.

Design of experiments for a surrogate model with continuous, integer and categorical variables

The first application relates to the optimal design of an electric engine with respect to the rotor component. The purpose of this study is to choose a set of values of the design variables that will be used to perform simulations of the engine operation in order to build surrogate models of the maximal power and of the maximal torque of the engine. The simulations are performed with the FEMM simulator1 designed for electromagnetic problems. The design variables are composed of one continuous variables (the rotor length), 2 integer variables (number of wires and number of coils) with respectively 4 and 11 feasible values and one categorical variable that characterizes the type of rotor geometry (with 16 possible types of rotor). For this problem, we used a tensorized kernel composed of a Gaussian kernel for continuous variables and the kernel (92) based on the Hamming distance applied to mixed binary variables resulting from the encoding of integer and categorical variables. Here we compare designs of experiments built with adapted-greedy approach to the projected LHS (algorithm 8) with a number of points ranging from 20 to 100 points. Surrogate models based on Gaussian processes and adapted kernel for mixed continuous and categorical variables ( [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF], [START_REF] Qian | Sliced Latin Hypercube Designs[END_REF] ) are built from the simulations performed on these designs. The accuracy of the obtained surrogate models is assessed by

1 ´Q2 pF, F q " ř Nv i"1 pF i ´F i q 2 ř N i"1 pF i ´F q 2 , ( 98
)
where Q 2 , the predictive accuracy coefficient, is computed on a new set of validation points of size N v , N " N v `nDOE is the total number of simulations including the points of the DoE, F is the vector of simulated values, F is the vector of predicted values (via the surrogate model) and F is the mean of F . In our experiments, we choose N v " 300 validation points obtained by a projected LHS method. In Figure 21, we study the predictivity of each model built from the designs of experiments for the two simulator outputs of interest, the maximal torque and the maximal power of the electric machine. We achieve a better predictivity (smaller values of 1 ´Q2 ) of the surrogate models for the adapted-greedy method compared to projected LHS for both responses of interest and for all the design sizes.

Design of experiments for cyclic-symmetry problems involving continuous and binary variables

The study of such cyclic-symmetry problems is driven by an application of optimal design for turbine blades in an helicopter engine [START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF], which involves mixed continuous and binary input variables.

We first apply our novel DoE methods on a set of benchmark functions (all with cyclic symmetry) with the aim of estimating the expectation of the actual functions. The second application presented in this section focuses on the choice of the initial points of an optimization method.

For both applications, the functions of interest are not affected by cyclic permutations of the vector of binary variables:

f px, y 1 q " f px, yq, @px, yq P R m ˆt0, 1u n (99) @y 1 P t0, 1u n such that Di P t1, 2, . . . , nu Rot i py 1 q " y, with Rot i pyq is the rotation of y by i positions.

For the MDS-greedy and the adapted-greedy methods, we use a tensorized kernel with, for the binary variables, the soft kernel [START_REF] Smola | A Hilbert space embedding for distributions[END_REF] with the necklace distance ( [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF])

d neck py, y 1 q " min i"1,2,...,n d H py, Rot i py 1 qq, (100)
where d H denotes the Hamming distance. A necklace of binary variables y P t0, 1u n is defined as the equivalence class of binary vectors considering all rotations as equivalent vectors (analogy with a necklace of n beads of 2 colors represented by 0 and 1 values): N y " tRot i pyq @i P t1, 2, . . . , nuu.

Results for expectation estimation

The benchmark consists of 6 different functions chosen from [START_REF] Dixon | The global optimization problem: An introduction[END_REF][START_REF] Luksan | Test Problems for Nonsmooth Unconstrained and Linearly Constrained Optimization[END_REF][START_REF] Neumaier | Neumaier's collection of test problems for global optimization[END_REF][START_REF] Chen | Super-Samples from Kernel Herding[END_REF], listed in table 6. These functions were set up originally for continuous optimization purpose. Therefore, we transform them into mixed binary functions with cyclic symmetry following the methodology of [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF]. The problem dimensions are listed in For each DoE δ n DoE , the accuracy of the estimations is rated by the root mean squared errors of the function expectations at the different necklaces

Errorpδ n DoE q " g f f e 1 n neck n neck ÿ h"1 pF ref h ´F DoE h q 2 , ( 101
)
where n neck is the number of necklaces, F ref The results of 20 repeated runs of the 6 methods with a design size n DOE " k n neck , k ranging from 1 to 27, are shown in Figure 22. Figure 23 provides the distributions of estimation errors for each necklace obtained by one run of the 6 methods on Branin test case with a design size of 100 points (k " 5). We observe on both Figures that the approximation errors are much smaller for greedy-MDS and adapted-greedy methods compared to the four others whereas the two projected methods based on LHS and Sobol sequences provide slightly better global results on the benchmark compared to standard methods with rounded values (Figure 22).

Results for initial design of optimization

The aim of the turbine application from [START_REF] Moustapha | Conception robuste en vibration et aéroélasticité des roues aubagées de turbomachines[END_REF][START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF] is the design of the turbine blades of an helicopter engine in order to minimize the vibrations of the compressor. This optimization problem involves one continuous parameter controlling the frequency amplitude and a vector of 12 binary variables describing the repartition of two reference blades on the turbine disk. In [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF], the authors propose an adapted optimization method based on a derivative-free trust-region method that uses the necklace distance (100) to take into account the cyclic symmetry of the problem. In this work, the initial set of points are chosen by a LHS procedure with rounded values for binary variables. We compare here these optimization results with the results obtained with the same optimization method coupled with two types of initial DoE: the projected LHS and adapted-greedy methods. As for previous benchmark functions, the chosen kernel for the latter is a tensorized kernel of a Gaussian kernel for continuous variables and the softmax kernel with necklace distance. The size of the initial DoE is n `m `1 " 14 and 100 repetitions of each DoE methods and associated optimizations are run. The results are reported in data and performance profiles in Figures 24 and25, respectively. As explained in [START_REF] Moré | Benchmarking Derivative-Free Optimization Algorithms[END_REF], the data and performance profiles with respect to the number of function evaluations are convenient tools to compare the performances of derivative-free optimization methods, counting the ratio of successful runs for each solver with respect to a chosen criterion, the number of simulations in our study. We consider that an optimization run is successful if the best solution z obtained within a given number of function evaluations satisfies f pz 0 q ´f pzq ě p1 ´τ qpf pz 0 q ´f ˚q, [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF] where, in the sequel, f ˚denote the best function value found by any solver (or the globalminimum value, if known), z 0 is the best point of the initial point for each run and τ is the desired accuracy, a user-defined tolerance value (in our tests, τ " 10 ´5). If a run does not provide a solution that satisfies [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF], we consider that it fails. The data profile displays the ratio of successful runs of each solver over the total number of runs with respect to the number of simulations (normalized by the number of variables: n `m `1 " 14), whereas the performance profile displays this ratio with respect to a performance ratio that is, in our case, the ratio between the number of evaluations needed by a given method to satisfy the condition [START_REF] Tran | An Adapted Derivative-Free Optimization Method for an Optimal Design Application with Mixed Binary and Continuous Variables[END_REF] and the smallest number of evaluations for all the compared methods.

Figures (24 and 25) show that the optimization method coupled with DoEs generated by projected LHS is the most robust method with 90% of successful runs, while the optimization coupled with adapted-greedy method succeeds to solve 88% of problems. The optimization method coupled with LHS with rounded values for binary variables is the least robust method with only 82% of successful runs. The optimization method coupled with adapted-greedy DoE appears to be less efficient for small budgets of simulations. For this kind of derivative-free optimization methods coupled with expensive simulators, the sizes of initial set are generally small in order to limit the number of expensive simulations. In the presented application, the size of the initial set ( 14) is small compared to the number of necklaces (equal to 352). This can explain the small differences between the results of the optimization method coupled with the three DoE methods: the exploration of the input space is done essentially along the optimization iterations and not during the initial phase.

Design of experiments for times series

To illustrate the ability of the proposed DoEs to be applied to kinds of variables other than the mixed discrete case, we propose in this section to apply the greedy-MDS and the adapted-greedy approaches coupled with the global alignment kernel (see [START_REF] Teymur | Optimal quantisation of probability measures using maximum mean discrepancy[END_REF] in Section 4) to two time series examples and compare the results with those obtained with some state-of-the-art approaches.

We consider two analytical functions with a functional input random variable V ptq that is known through a sample of 200 realizations. Two cases are studied: V is either a standard Brownian motion (Figure 26) or a Max-stable process (Figure 27). The aim of these experiments is to estimate the expectation of the function with respect to the functional variable V ptq with a limited number of samples of V . To achieve this goal, [START_REF] Gilles | Introduction to optimal vector quantization and its applications for numerics[END_REF][START_REF] El Amri | Data-driven stochastic inversion via functional quantization[END_REF] propose a greedy functional data-driven quantization approach coupled with a functional principal component analysis (FPCA) to reduce the dimension of V to a finite (small) dimension. The impact of the number of components in FPCA on the explained variance of the functional data for two types of uncertainties (Brownian and Max-stable processes) is discussed in [START_REF] El Amri | Set inversion under functional uncertainties with joint meta-models[END_REF]. For the presented numerical experiments, we choose 8 components that explain respectively 97.6% and 70% of the data variance. We propose here to compare this quantization method to our two DoE methods applied directly in the functional domain with the global alignment kernel K GAK introduced in Subsection 3 for time series. The greedy-MDS method relies on the associated kernel distance defined as

d 2 GAK puptq, vptqq " K GAK puptq, uptqq `KGAK pvptq, vptqq ´2K GAK puptq, vptqq. (103)
A standard LHS method applied in the reduced space obtained by FPCA is also evaluated. The two studied functions are defined as where we fix x 1 " 50 and x 2 " ´80, and

f 1 : px, V q Ñ x 2 1 `2x 2 2 ´0.3 cosp3πx 1 q ´0.4 cosp4πx 2 q `0.7 `T ż 0 e Vt dt, (104)
f 2 : px, V q Ñ 2x 2 1 `x2 2 ´0.3 cosp3πx 1 q ´0.4 cosp4πx 2 q `0.7 `T ż 0 2 sinpV t qdt, (105)
where x 1 " 50 and x 2 " ´80.

We apply the 4 DoE methods for the functions f 1 and f 2 and the two functional data, Brownian motion and max-stable process, with various sizes of DoE equal to 20, 40 and 60 chosen among the 200 available realizations of V . The performances of the methods are measured in terms of accuracy and computational time in Tables 7,8 The quantization method which is dedicated to the expectation estimation gives often the best results but our two methods provide close results in terms of accuracy whereas the results of the LHS method applied in the reduced space are not as accurate in general. In terms of computational time, our two methods are more efficient than the quantization method implementation of [START_REF] El Amri | Data-driven stochastic inversion via functional quantization[END_REF]. Moreover, these methods directly consider the functional variable and do not require to apply a dimension reduction method. DoE techniques as low discrepancy methods or standard space filling approaches coupled with the projection from continuous space to integers representing the distinct indices of discrete variables. A second proposed approach relies on kernel-embedding of probability distributions extended to mixed continuous and discrete space, which leads to two new methods: greedy-MDS and adapted-greedy. These methods are generic and can address other types of objects than mixed discrete and continuous variables, provided that a suitable kernel is available. A variety of useful kernels are presented: in particular, the softmax kernel that can be built from any distance and ensures the positive-definite property and the global alignment kernel that is suited to time series.

We illustrate the performances of the proposed DoE methods in various context: building a training data set for a surrogate model handling continuous, integer and categorical variables, and for the estimation of the expectation of a function with cyclic symmetry for a vector of binary variables, choosing initial points of a mixed binary and continuous shape optimization problem for an helicopter engine, and finally DoE for time series for expectation estimation of a function that depends on a functional variable.

Future studies can be carried on with other types of objects and associated kernels (images, graphs ...) and also with different operational objectives. For instance, the exploration step Global optimization, is an attractive objective for many applications dealing with complex multimodal functions. We mentioned in Part I that DFO methods, for instance, DFOb-d H or adapted method DFOb-d neck , guarantee convergence but only toward a local minimum. As we test and compare DFOb-d neck with RBFOpt (which is declared as a global method) in the previous part, we observe that for some instances, our DFO method is comparable, or even more efficient than RBFOpt with the same accuracy. Nevertheless, DFO trust-region methods may be trapped in a local minimum without exploring regions that includes a globally optimal solution. There are techniques have been used in DFO literature for escaping from local minima, for example multi-start mechanism using in RBFOpt [START_REF] Costa | RBFOpt: An open-source library for black-box optimization with costly function evaluations[END_REF] or in the work of [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF], restart strategy in [START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF], or multistart using multi-level single linkage or maximum information from previous evaluations in GORBIT [START_REF] Wild | Derivative-Free Optimization Algorithms for Computationally Expensive Functions[END_REF].

Our idea to circumvent this drawback is a restart strategy coupled with an augmented DoE strategy. Section 1 presents the principle of an augmented DoE based on the sampling of the discrete variable space. Section 2 introduces three restart strategies including two original proposals. Section 3 gives preliminary results. Section 4 concludes this third part of the thesis.

1 Augmented DoE strategies for sampling a discrete variable space

The purpose of this section is the improvement of the standard "add-points" step (see Chapter 3, Part I) as a diversifying strategy to help the exploration step for new values of the binary set y.

As we mentioned, the current "add-points step" chooses a user-defined number of new points which only need to satisfy the trust-region and no-good cuts constraints without any requirement regarding the available simulated points and the spreading of the points. We propose here a strategy that uses the information of the available simulated points to anticipate whether we should explore more or less the continuous search space for a given binary configuration (or necklace). The proposed method is an extension of the DOE methods proposed in Part II integrating a given initial design of experiments (the available simulated points) to be completed.

Our strategy is based on the sampling method proposed in [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF] for general mixed continuous and categorical variables in the Bayesian optimization context. In the context of cyclic-symmetry problems, we adapt the sampling method by considering each necklace as one category of a categorical variable, called level in [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF].

The authors of [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF] define a discrete probability distribution in the categorical space that takes into account the density of available evaluated points at each level. The discrete probability of level c i is defined as

p k,i " α k p g k,i `p1 ´αk qp m k,i , @i " 1, 2, . . . , n levels , ( 106
)
where k is the optimization iteration counter, α k P r0, 1s is a weighting parameter, the quantity p g k,i stands for the probability that the level c i has not been fully explored (it is large when the number of simulations associated with this level is small) and is computed as:

p g k,i " 1 ´ˆn k,i n k ˙l , ( 107
)
where n k,i and n k are respectively the number of evaluated points in the level c i and the total number of evaluations of the objective function over all the levels at iteration k; l ą 0 is a user-defined parameter (that we shall set to 1 2 in the implementation); the term p m k,i is the probability that the level c i has a high potential of containing the minimum of the objective function, which is approximated by: We realize two experiments: one with five adding points (see Figure 29), the other one with ten points (see Figure 30). On the left-hand side, the points provided by the standard procedure are displayed, whereas on the right-hand side points added by the augmented DoE algorithm are plotted. We observe, for two experiments, that the new procedure proposes points on the levels that were not already explored and with a minimization potential learnt from the available simulated points.

p m k,i " S R k,i ř n levels j"1 S R k,j , (108)

Restart strategies

This section presents different restart strategies that allow the DFO trust-region methods DFOb to escape from local minima based on ideas proposed in [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF][START_REF] Nannicini | On the implementation of a global optimization method for mixed-variable problems[END_REF][START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF] for continuous derivativefree trust-region method. In [START_REF] Gutmann | A Radial Basis Function Method for Global Optimization[END_REF][START_REF] Nannicini | On the implementation of a global optimization method for mixed-variable problems[END_REF], the authors use a mechanism to restart the algorithm from completely new initial points whenever certain conditions are satisfied. Beside this hard- restart strategy, the authors in [START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF] propose an adapted soft restart which allows one to use the information from the evaluated points and restart with the algorithm scale of 1.1∆ 0 : the algorithm starts with a trust-region radius set to 1.1∆ 0 , where ∆ 0 is the initial radius.

Another attractive track is the multistart strategy, proposed for example in the method GORBIT [START_REF] Wild | ORBIT: Optimization by Radial Basis Function Interpolation in Trust-Regions[END_REF]. The principle of this tactic is to start the algorithm from different starting points in parallel.

In our work, we integrate the above hard restart in the DFOb method, and we propose two other restart methods: complete soft restart and partial soft restart. The next three subsections will give more details on the three approaches.

We choose two conditions proposed in [START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF] to trigger the restart strategies, we involve restart when either:

• The trust-region radius at the current iteration k, ∆ k , is small,

• The decrease of the objective function over consecutive iterations is too small.

Remark that we also use the re-initialized value of the trust-region radius followed by the adapted soft restart in [START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF]: we increase the trust-region radius to 1.1∆ 0 .

Hard restart

The idea of hard restart is to restart the algorithm from totally new starting points which are generated from the same mechanism as generating the initial interpolation set. The advantage of this restart type is its simplicity of implementation. When the algorithm detects that it is trapped in the valley of a local solution, the algorithm stops and restarts from a new initial interpolation set.

However, this strategy has some drawbacks also. First, it is not guaranteed that the algorithm will not converge towards in the same local minimum again. Second, it is clear that with this strategy, the users do not use any information, i.e., objective function evaluations, from the previous runs, which is not efficient, especially in the context of blackbox optimization where we want to save expensive simulations.

Complete soft restart

Our first proposal restart strategy restarts the algorithm with a subset of starting points that are chosen among already-evaluated objective-function points. The intuition comes from the fact that, ordinarily, the condition of restarting is met after a certain number of iterations that may provide a good exploration of the domain to set up a new initial interpolation set.

The advantage of this method is that we have no waste of new simulations to set up the new initial interpolation set. However, it can happen that the set of evaluated points are not very well distributed over the domain, yielding a set of starting points that can be biased. Also, like the hard restart, it is possible that the algorithm converges towards a local solution that was already explored.

Partial soft restart

Our second proposal improves the previous method by starting from a set of points that includes both already-evaluated points and new points using the augmented DOE method proposed in section 1.

The advantage of this proposal is that we use the available information from previous simulations, while ensuring the good properties of the new starting point set.

Numerical results

This section presents preliminary comparative numerical results. Subsection 3.1 and Subsection 3.2 report numerical results on the DS-MM and Ackley functions respectively.

We compare five options of DFOb-d neck :

• The standard DFOb-d neck with the standard option of add-point step (named Standard),

• DFOb-d neck with the augmented DoE add-points step (Augmented),

• DFOb-d neck coupling the augmented DoE add-point step and hard restart (A-hard restart),

• DFOb-d neck coupling the augmented DoE add-point step and complete soft restart (Acomplete-soft restart),

• DFOb-d neck coupling the augmented DoE add-point step and partial soft restart (A-partialsoft restart).

We set the maximum number of evaluation to 300. The five experiments both share two reference options and three restart options with the same initial design (generated by adapted-greedy algorithm, nDoE " m `n `1). The value of the parameter α k in the augmented DOE method is set to 0.01. The maximum number of consecutive discarded iterations is 50.

Results for the DS-MM function

The first test function is the piecewise-defined DS-MM function (see equation ( 111)) involving one continuous and five binary variables (corresponding to 8 individual one-dimensional functions) as described in section 1.

Figure 31 shows the best current objective-function value achieved by the 5 options, with respect to the number of objective-function evaluations. The standard, augmented, and coupling augmented and hard restart DFOb-d neck options yield the same solution with the associated objective-function value: ´1.948 within 43, 38 and, 38 evaluations respectively. The A-completesoft restart method yield the globally optimal solution with associated objective-function value ´2.328 within 66 evaluations. The A-partial-soft restart finds a solution with objective-function value ´2.319 within 67 evaluations, a minimum local solution at the level (y value) corresponding to the global solution. Figure 32 shows the solutions obtained by the five options.

Results for the Ackley function

We now perform a similar study for the well-known Ackley function. The Ackley function is originally for continuous variables and given by:

f pxq " ´ae ´bd 1 m m ř i"1 x 2 i ´e 1 m m ř i"1 cospcx i q `a `e, (112)
where the parameter values are set to a " 20, b " 0.2, c " 2π and m " 2. The function is evaluated on the square x i P r´5, 5s, i " 1, 2. We use the cyclic-symmetry form construction from Chapter 5 of Part I with n " 4 binary variables (corresponding to 6 individual twodimensional functions to define the resulting piecewise-defined function). The global minimizer is x " p0, 0q, y " p0, 0, 0, 0q. The associated objective-function value is ´1.7634. The plots of the 6 individual two-dimensional functions associated with the 6 necklaces are shown in Figure 33.

Conclusion

In this part of the thesis, we presented several original contributions to escape from local minima. A first contribution improves the "add-point" step of DFOb methods based on an original augmented design of experiment based on a discrete distribution taking into account the values of the available objective-function evaluations and the associated point distribution.

Another contribution is the introduction of two restart strategies based on new initial interpolation sets using previously simulated points. These methods were compared with the standard DFOb method, and with some reference restart strategies implemented in DFOb. The preliminary results obtained on two analytical functions with cyclic symmetry involving several local minima are encouraging. The proposed methods can also be used with DFOb-d H for general mixed continuous and binary problems. y 1 " p0, 0, 0, 0q, y 2 " p0, 0, 0, 1q, y 3 " p0, 0, 1, 1q, y 4 " p0, 1, 0, 1q, y 5 " p0, 1, 1, 1q, y 6 " p1, 1, 1, 1q).

General conclusion and perspectives

Conclusion

In this work, we presented several contributions to solve blackbox optimization problems involving computationally expensive simulators with mixed continuous and discrete variables. Numerous applications of this type of problems arise in the domain of optimal design: in particular, our research work was motivated by an application provided by Safran Tech for the design of turbomachine of aircraft.

• The first part of the thesis considers optimization methods for solving derivative-free problems. We presented an in-depth review of DFOb-d H , a derivative-free trust-region method for mixed continuous and binary problems based on the Hamming distance. A first contribution was the convergence proof of this method which was not presented in the original work that proposed this method. The second contribution is a derivative-free trust-region method adapted to cyclic-symmetry mixed continuous and discrete problems: DFOb-d neck method. This method relies on an adapted distance for cyclic-symmetry problems based on the "necklace" concept from combinatorics. Reformulations of the subproblems involved in the trust-region method DFOb were built to preserve the MIQP structure of the subproblems. The convergence proof of this new method is given. The proposed method was compared to three state-of-the art optimization methods: NOMAD, RBFOpt and DFLBOX and to the original DFOb-d H method based on usual Hamming distance. This comparison was performed on a benchmark of 25 functions with cyclic symetry that was set up from various published benchmarks. The comparison is also done on an application problem provided by Safran Tech. DFOb-d neck shows robust results within a limited number of function evaluations.

• The second part of this thesis is dedicated to the design of experiments for mixed continuous and discrete variables. We propose two extensions of a method based on kernel embedding of probability distribution for mixed variables space. The proposed methods can be applied to various types of variables, provided a positive defined kernel. The performances of the methods are illustrated on examples with mixed discrete variables, cyclic-symmetry problems and also time series with dedicated kernels.

• Finally, in the third part of the thesis, new ideas are proposed and evaluated to improve the robustness of the optimization method DFOb-d neck with respect to not getting trapped in local minima. Some encouraging preliminary results are presented.

Perspectives This work opens up several research avenues for future studies. From optimization perspective:

• First, apply the proposed method to an application of the optimal design problem in Safran's industrial context with the real complex black-box simulators (rather than on the simplified version of it).

• Regarding the problem of escaping from local minima:

1. Extend the preliminary study described in Part 3: coupling the DOE methods proposed in Part 2 with restart strategies proposed in [START_REF] Cartis | Escaping local minima with derivative-free methods: A numerical investigation[END_REF].

2. Couple an exploration step based on global surrogate methods (EGO, RBFOpt) with the proposed trust-region method.

• Extend the proposed optimization methods DFOb for noisy problems: extension of works of [START_REF] Wild | Derivative-Free Optimization Algorithms for Computationally Expensive Functions[END_REF] for mixed continuous and discrete variables.

• Extend the proposed optimization method DFOb-d neck to general constrained problems. One technique can be used is a progressive barrier derivative-free trust-region in [START_REF] Audet | A progressive barrier derivative-free trust-region algorithm for constrained optimization[END_REF]: by using a constraint violation function to ensure the algorithm to find the feasible solutions.

From DoEs points of views:

• Evaluate the proposed DOE methods on various types of problems: for instance, images or applications with prior information.

A Mathematical Background

We revise some concepts necessary to follow the work in this thesis. Some results in Section 1 and 2 are used in the proof of convergence of our method (see Chapters 4 and 5, Part I). Section 3 presents quadratic programming which is the core of our algorithm.

Vector and matrix norms

We begin by introducing the notations used in this thesis: vector and matrix norms, condition number, quadratic programming and mixed integer quadratic programming. Besides satisfying the three norm properties (non-negativity, scalar multiply, triangle inequality), they also satisfy the "submultiplicative" property ∥AB∥ ď ∥A∥∥B∥.

B Encoding techniques and list of benchmark functions 1 Categorical variables encoding techniques

There is a large number of engineering optimization problems with both continuous and discrete variables. In applications, the discrete variables are sometimes categorical ones; i.e., they refer to a list or set of categories. Many optimization algorithms cannot deal directly with this type of variables: the categories must often be converted into numbers.

There are two types of categorical variables: nominal and ordinal variables. The term nominal refers to no intrinsic ordering, for instance, gender (male and female). Unlike nominal variables, ordinal variables have a clear ordering.

This section aims to give the readers a short review about different ways of encoding categorical variables. Since the work is more focused on mixed continuous and binary variables, we will give more detail on the techniques that convert categorical variables to binaries.

Binary encoding. This is one of the simplest way to encode categorical variables into quantitative variables (numbers). First, the categorical variables are converted to numeric (integer) numbers, denoted by l 1 , l 2 , . . . , l N , then one applies the following binary encoding to l k , k " 1, 2, . . . , N :

l k " tlog 2 pnqu`1 ÿ i b k i 2 i´1 ,
where b k i P t0, 1u. Thus, the encoding of l k is pb k 1 , b k 2 , . . . , b k tlog 2 pnqu`1 q, k " 1, 2, . . . , N .

One-hot encoding. This includes encoding each categorical variable with a different binary or boolean variable (also called dummy variable) which takes values 0 or 1, indicating whether a category is present in an observation. `x4 2 , f 2 pxq "p2 ´x1 q 2 `p2 ´x2 q 2 , f 3 pxq "2e `x2 2 , f 2 pxq "p2 ´x1 q 2 `p2 ´x2 q 2 , f 3 pxq "2e f i pxq, f 1 pxq " ´apx 1 ´x3 q 2 `px 2 ´x4 q 2 , f 2 pxq " ´apx 3 ´x5 q 2 `px 2 ´x6 q 2 , f 3 pxq " ´apx 5 ´x1 q 2 `px 6 ´x2 q 2 . Problem B.6. ROSEN-SUZUKI.

F pxq " maxtf 1 pxq, f 1 pxq `10f 2 pxq, f 1 pxq `10f 3 pxq, f 1 pxq `10f 4 pxqu, 3 ´7x 4 ´120q, f 3 pxq "f 1 pxq `10p5x 2 1 `8x 2 `px 3 ´6q 2 ´2x 4 ´40q, f 4 pxq "f 1 pxq `10p0.5px 1 ´8q 2 `2px 2 ´4q 2 `3x 2 5 ´x6 ´30q, f 5 pxq "f 1 pxq `10px 2 1 `2px 2 ´2q 2 ´2x 1 x 2 `14x 5 ´6x 6 q, f 6 pxq "f 1 pxq `10p´3x 1 `6x 2 `12px 9 ´8q 2 ´7x 10 q, f 7 pxq "f 1 pxq `10px 2 1 `5x 11 ´8x 12 ´28q, f 8 pxq "f 1 pxq `10p4x 1 `9x 2 `5x 2 13 ´9x 14 ´87q, f 9 pxq "f 1 pxq `10p3x 1 `4x 2 `3px 13 ´6q 2 ´14x 14 ´10q, f 10 pxq "f 1 pxq `10p14x F pxq " 100px 2 ´x2 1 q2 `p1 ´x1 q 2 . Problem B.12. HS3. F pxq " x 2 `10 ´5px 2 ´x1 q 2 . Problem B.13. HS29log.

f 1 pxq "x 2 1 `x2 2 `2x
F pxq " log 10 p100px 2 ´x2 1 q 2 `p1 ´x1 q 2 q.

Dixon-Szegö benchmark Problem B.14. Branin.

F pxq " px 2 ´p 5.1 4π 2 qx 2 1 `5 π x 1 ´6q 2 `10p1 ´1 8π q cospx 1 q `10.

Problem B.15. Camel.

F pxq " p4 ´2.1x

´2. 329 .

 329 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Five points (blue stars) to complete the set of ten evaluated points (green circles) by the standard procedure (left) and by the augmented DoE method (right). . . . . . . Ten points (blue stars) to complete the set of ten evaluated points (green circles) by the standard procedure (left) and by the augmented DoE method (right). . . . . . . 31 The best current objective-function values versus number of evaluations of DFObd neck with standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the DS-MM function. . . . . . . . . . . . . . . . . . . . . 32 Solutions obtained by DFOb-d neck with standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the DS-MM function. . . . .

Figure 2 (

 2 left) illustrates the case with n " 23 blades with the two possible types of blade geometry. Figure2(right) displays, for the case of n " 6 blades, all the distinct arrangements and their equivalent configurations obtained by rotation.

Fig. 2 .

 2 Fig. 2.: Examples of blade configurations for two different geometries (from [74]).

Fig. 3 .

 3 Fig. 3.: Illustration of a case with 23 blades with groups AA ´BB (from [74]).

1 ( 2 Step 1 . 5 (

 1215 Exploitation phase): Algorithm 4 (cf. Section 3) • If unsuccessful Step 1: go to Step Exploration phase): Algorithm 5 (cf. Section 4) • If successful Step 1.5: k Ð k `1 and go to Step 0

  and go to Step 0 else: Ñ STOP a Algorithm 3 adds possibly in Z k new points to be simulated.

Theorem 4 . 1 .

 41 Let Assumptions 4.1 and 4.2 hold. Then,

Definition 5 . 1 .

 51 (Necklace) Let k and n be strictly positive integers. In combinatorics, a kary necklace of length n is an equivalent class of n-character strings over an alphabet ř k " ta 1 , . . . , a k u of size k, taking all rotations as equivalent. It represents a structure with n circularly connected beads with k available colors.Let N eckpnq be the set containing all the necklaces of length n. For instance, considering the necklace with 2 colors and 4 beads (n " 4, k " 2, ř k " t0, 1u), we have the following list of representative necklaces (in red color) and all rotations (in black):

Figure 5

 5 Figure 5 also shows that the new necklace distance improves the efficiency of the DFOb method when addressing cyclic-symmetry optimization problems, as DFOb-d H solves only 28% of the problems with the minimal number of objective-function evaluations.

Fig. 5 .

 5 Fig. 5.: Performance profiles of the five solvers for the 25 analytical problems

Fig. 6 .Fig. 7 .

 67 Fig. 6.: Data profiles of the five solvers for the 25 analytical problems

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8.: Performance profiles of the five solvers for the 10 highest dimension analytical problems

Fig. 11 ." 10 ´5Fig. 12 . 10 ´5Fig. 13 .

 1110121013 Fig. 11.: Performance profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10´5

Fig. 14 ." 10 ´3Fig. 15 . 10 ´3Fig. 16 .

 1410151016 Fig. 14.: Performance profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10´3

Fig. 17 .

 17 Fig. 17.: The four solutions found by DFOb-d neck for the blade design application (the continuous component has the same common value: x ˚" 0.03)

Fig. 18 .

 18 Fig. 18.: Performance profiles of DFOb-d neck with nDoE " m `n `1 and nDoE " 2pm `nq `1 for the 25 analytical problems

Fig. 19 .Fig. 20 .

 1920 Fig. 19.: Data profiles of DFOb-d neck with nDoE " m `n `1 and nDoE " 2pm `nq `1 for the 25 analytical problems

Algorithm 8 :

 8 Projected DoE for discrete variables Input: n, nDoE 0. Pre-processing • Compute all n levels distinct discrete arrangements of size n and arrange them in the matrix L • Compute the weight, w " 1 n levels 1. Main computation • Build a continuous DoE of size n DoE in dimension 1: U " tu 1 , . . . , u n DoE u • Assign discrete values according to the following rule:

Fig. 21 .

 21 Fig. 21.: Accuracy of the surrogate models of the maximal torque (left) and the maximal power (right) built with designs of experiments obtained by the projected LHS method (in blue) and by the adapted-greedy algorithm (in red). The size of the designs, n DOE , ranges from 20 to 100.

h and F

 F DoE h are respectively the empirical expectations of the function for the Monte-Carlo sampling points and the points of the DoE for the necklace h.

Fig. 22 .

 22 Fig. 22.: Median RMSE (in log-scale) of the estimation of benchmark function expectation obtained with 5 repetitions of the 6 methods: a random sampling and a standard LHS with rounded values for binary variables, 2 projected methods based on Sobol sequence and projected LHS, the greedy-MDS and adapted-greedy methods.

Fig. 23 .

 23 Fig. 23.: Distributions of estimation errors for each necklace of the Branin function obtained with 6 methods: a random sampling and a standard LHS with rounded values for binary variables, 2 projected methods based on Sobol sequence and projected LHS, the greedy-MDS and adaptedgreedy methods. The size of the DoE is 100 points.

Fig. 24 .

 24 Fig. 24.: Safran's application: data profiles of the optimization runs with 100 initial DoE obtained with 3 methods: LHS with rounded values for binary variables, projected LHS, and adapted-greedy methods.

Fig. 25 .

 25 Fig. 25.: Safran's application: of the optimization runs with 100 initial DoE obtained with 3 methods:LHS with rounded values for binary variables, projected LHS, and adapted-greedy methods.

, 9 ,

 9 and 10. The errors are computed as the relative errors of the expectation estimations based on the obtained DoE and the estimation based on the 200 realizations. The results of LHS and quantization methods are the means of the results obtained for 50 repeated runs.

Fig. 26 .

 26 Fig. 26.: 200 realizations of a Brownian motion. The thick blue curves are the 10 curves selected by the first iterations of the adapted-greedy algorithm.

Tab. 7 .:

 7 Computational times and relative errors of the expectation estimation obtained for the function f 1 with the functional data of Brownian motion. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold.

Fig. 27 .

 27 Fig. 27.: 200 realizations of a Max-stable process. The thick blue curves are the 10 curves selected by the first iterations of the adapted-greedy algorithm.

Fig. 28 .

 28 Fig. 28.: The eight individual one-dimensional functions for DS-MM example. The global minimum is located at x " 0.808, y " p0, 0, 1, 0, 1q, the objective-function value is ´2.329.

Fig. 29 .

 29 Fig. 29.: Five points (blue stars) to complete the set of ten evaluated points (green circles) by the standard procedure (left) and by the augmented DoE method (right).

Fig. 30 .

 30 Fig. 30.: Ten points (blue stars) to complete the set of ten evaluated points (green circles) by the standard procedure (left) and by the augmented DoE method (right).

Fig. 31 .

 31 Fig. 31.: The best current objective-function values versus number of evaluations of DFOb-d neck with standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the DS-MM function.

Figure 34

 34 Figure 34 shows the best objective-function value achieved by each of the five options on the Ackley function, with respect to the number of objective-function evaluations. The solution found by the standard DFOb-d neck has objective-function value: ´1.336 and required 25 evaluations. The solution found by the augmented DFOb-d neck has objective-function value: ´1.7604 and required 92 evaluations. For A-hard restart the value is ´1.7521, within 80 evaluations, Acomplete-soft restart yields the objective-function value ´1.7541 within 90 evaluations, and the solution found by A-partial-soft restart gives objective-function value ´1.7607 and 110 evaluations.Figure 35 zooms on the main best objective-function values of the Figure 34.Figure 36 plots the solutions found by the five options; they all correspond to a same value of the binary-variable vector y corresponding to the first individual two-dimensional function.

Figure 35 zooms

 35 Figure 34 shows the best objective-function value achieved by each of the five options on the Ackley function, with respect to the number of objective-function evaluations. The solution found by the standard DFOb-d neck has objective-function value: ´1.336 and required 25 evaluations. The solution found by the augmented DFOb-d neck has objective-function value: ´1.7604 and required 92 evaluations. For A-hard restart the value is ´1.7521, within 80 evaluations, Acomplete-soft restart yields the objective-function value ´1.7541 within 90 evaluations, and the solution found by A-partial-soft restart gives objective-function value ´1.7607 and 110 evaluations.Figure 35 zooms on the main best objective-function values of the Figure 34.Figure 36 plots the solutions found by the five options; they all correspond to a same value of the binary-variable vector y corresponding to the first individual two-dimensional function.

Figure 36

 36 Figure 34 shows the best objective-function value achieved by each of the five options on the Ackley function, with respect to the number of objective-function evaluations. The solution found by the standard DFOb-d neck has objective-function value: ´1.336 and required 25 evaluations. The solution found by the augmented DFOb-d neck has objective-function value: ´1.7604 and required 92 evaluations. For A-hard restart the value is ´1.7521, within 80 evaluations, Acomplete-soft restart yields the objective-function value ´1.7541 within 90 evaluations, and the solution found by A-partial-soft restart gives objective-function value ´1.7607 and 110 evaluations.Figure 35 zooms on the main best objective-function values of the Figure 34.Figure 36 plots the solutions found by the five options; they all correspond to a same value of the binary-variable vector y corresponding to the first individual two-dimensional function.

Fig. 32 .

 32 Fig. 32.: Solutions obtained by DFOb-d neck with standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the DS-MM function.

Fig. 33 .

 33 Fig. 33.: Ackley function with two continuous variables and four binary variables (corresponding 6 levels:y 1 " p0, 0, 0, 0q, y 2 " p0, 0, 0, 1q, y 3 " p0, 0, 1, 1q, y 4 " p0, 1, 0, 1q, y 5 " p0, 1, 1, 1q, y 6 " p1, 1, 1, 1q).

Fig. 34 .

 34 Fig. 34.: The best current objective-function values versus number of evaluations of DFOb-d neck with the five options standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the Ackley function.

Fig. 35 .

 35 Fig. 35.: Zoom on the best current objective-function values versus number of evaluations of DFOb-d neck with the five options.

Fig. 36 .

 36 Fig. 36.: The two solutions obtained with DFOb-d neck and the five options on the Ackley function.

For a vector

  x P R n , p´norms (p ě 1) of x are given by ∥x∥ p "´n ÿ i"1 |x i | p ¯1 p ,and the 8´norm by ∥x∥ 8 " max1ďiďn |x i |.Similarly, we define the p´norm of matrices. Let A " pa ij q mˆn be a m ˆn matrix. The p´norms of A, p ě 1, are defined by ∥A∥ p " sup x‰0

1 a ji px j ´Pji q 2 . 1 a ji px j ´Pji q 2 .

 1212 [0.36890 , 0.46990 , 0.10910 , 0.03815] , [0.11700 , 0.43870 , 0.87320 , 0.57430] , [0.26730 , 0.74700 , 0.55470 , 0.88280]] c = [ 1 . 0 , 1 . 2 , 3 . 0 , 3 . 2 ] 0.8307 , 0.3522 , 0 . 8 7 3 2 ] , [ 0 . 0 1 2 4 , 0.3736 , 0.2883 , 0 . 5 7 4 3 ] , [ 0 . 8 2 8 3 , 0.1004 , 0.3047 , 0 . 1 0 9 1 ] , [ 0 . 5 8 8 6 , 0.9991 , 0.6650 , 0.0381] ] c = [ 1 . 0 , 1 . 2 , 3 . 0 , 3 . 2 ] Problem B.19. Shekel7. a = [ [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 3 . 0 , 2 . 0 , 5 . 0 ] , [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 7 . 0 , 9 . 0 , 5 . 0 ] , [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 3 . 0 , 2 . 0 , 3 . 0 ] , [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 7 . 0 , 9 . 0 , 3 . 0 ] ] c = [ 0 . 1 , 0 . 2 , 0 . 2 , 0 . 4 , 0 . 4 , 0 . 6 , 0 . 3 ] Shekel10. a = [ [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 3 . 0 , 2 . 0 , 5 . 0 , 8 . 0 , 6 . 0 , 7 . 0 ] ,

  ´5 for each of the 50 runs of the blade design application for the five solvers (successful runs only) together with percentage of success (in parentheses) . . . . . . . . . . . . . . . . . . . . . . 14 Performance profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10 ´3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Data profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10 ´3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . found by DFOb-d neck for the blade design application (the continuous component has the same common value: x ˚" 0.03) . . . . . . . . . . . . . . 18 Performance profiles of DFOb-d neck with nDoE " m`n`1 and nDoE " 2pm`nq`1 for the 25 analytical problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 Data profiles of DFOb-d neck with nDoE " m `n `1 and nDoE " 2pm `nq `1 for the 25 analytical problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Number of function evaluations to reach f ˚up to τ " 10 ´3 for each of the 25 analytical problems for DFOb-d neck with nDoE " m`n`1 and nDoE " 2pm`nq`1 (successful runs only) together with percentage of success (in parentheses) . . . . . 21 Accuracy of the surrogate models of the maximal torque (left) and the maximal power (right) built with designs of experiments obtained by the projected LHS method (in blue) and by the adapted-greedy algorithm (in red). The size of the designs, n DOE , ranges from 20 to 100. . . . . . . . . . . . . . . . . . . . . . . . . .

	List of Figures	xi
	17	The four solutions
	x

. . . . . . . . . . . . . . . . . . . . . . . . . . 8 Performance profiles of the five solvers for the 10 highest dimension analytical problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Data profiles of the five solvers for the 10 highest dimension analytical problems . . 10 Number of function evaluations to reach f ˚up to τ " 10 ´3 for each of the 10 highest dimension analytical problems for the five solvers (successful runs only) together with percentage of success (in parentheses) . . . . . . . . . . . . . . . . . . . . . . 11 Performance profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10 ´5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Data profiles of the five solvers for the blade design application with 50 repetitions, accuracy τ " 10 ´5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Number of function evaluations to reach f ˚up to τ " 10 16 Distribution of the binary component of the 50 solutions found by DFOb-d neck for the blade design application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 Median RMSE (in log-scale) of the estimation of benchmark function expectation obtained with 5 repetitions of the 6 methods: a random sampling and a standard LHS with rounded values for binary variables, 2 projected methods based on Sobol sequence and projected LHS, the greedy-MDS and adapted-greedy methods. . . . . 23 Distributions of estimation errors for each necklace of the Branin function obtained with 6 methods: a random sampling and a standard LHS with rounded values for binary variables, 2 projected methods based on Sobol sequence and projected LHS, the greedy-MDS and adapted-greedy methods. The size of the DoE is 100 points. . . 24 Safran's application: data profiles of the optimization runs with 100 initial DoE obtained with 3 methods: LHS with rounded values for binary variables, projected LHS, and adapted-greedy methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Safran's application: of the optimization runs with 100 initial DoE obtained with 3 methods: LHS with rounded values for binary variables, projected LHS, and adapted-greedy methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 200 realizations of a Brownian motion. The thick blue curves are the 10 curves selected by the first iterations of the adapted-greedy algorithm. . . . . . . . . . . . 27 200 realizations of a Max-stable process. The thick blue curves are the 10 curves selected by the first iterations of the adapted-greedy algorithm. . . . . . . . . . . .

  ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 The best current objective-function values versus number of evaluations of DFObd neck with the five options standard, augmented add-point step, hard restart, complete soft restart, and partial soft restart on the Ackley function. . . . . . . . . . . . Computational times and relative errors of the expectation estimation obtained for the function f 1 with the functional data of Brownian motion. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are Computational times and relative errors of the expectation estimation obtained for the function f 1 with the functional data of Max-stable process. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Computational times and relative errors of the expectation estimation obtained for the function f 2 with the functional data of Brownian motion. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Computational times and relative errors of the expectation estimation obtained for the function f 2 with the functional data of Max-stable process. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are List of Algorithms Derivative based trust region algorithm, see [81] . . . . . . . . . . . . . . . . . . . 19 Improving poisedness of X via LU factorization, see [27] . . . . . . . . . . . . . . 26 Improving poisedness of Z in the mixed trust region B, extended from [27] . . . . 36 Exploitation phase of DFOb-d Bpx 0 , ∆ x q: continuous l 2 ball centered at x 0 of radius ∆ x Bpy 0 , ∆ y q : binary l 2 ball centered at y 0 of radius ∆ y B y 0 px 0 , ∆ x q : continuous l 2 ball with fixed binary part y " y 0 of radius ∆ x

		List of Tables	xiii
	M	8	: class of fully-linear (or fully-quadratic) model
	d H		: Hamming distance
	ACRONYMS d neck : necklace distance NOTATIONS Rot r pyq : rotation of y by r positions
	nDoE	: number of points in the interpolation set
		f m	BO BONMIN : objective function : Bayesian Optimization : Basic Open-source Nonlinear Mixed INteger programming DFO : Derivative Free Optimization : number of continuous variables
		x	DFOb	: DFO trust-region method with mixed binary variables : continuous-variable vector
		x	DFOb-d H : upper bound for x : DFOb with Hamming distance
		x	DFOb-d neck : DFOb with necklace distance : lower bound for x
		X X	DoEs EGO GA	: Design of Experiments : continuous interpolation set : Efficient Global Optimization : Genetic Algorithm : continuous normalized interpolation set
		m k	GP	: Gaussian Process : continuous quadratic model at iteration k
		y	IPOPT	: Interior Point OPTimizer : binary-variable vector
		z	LHS	: Latin Hypercube Sampling : mixed-variable vector z " px, yq
	List of Tables MADS MINLP n : number of binary variables : Mesh Adaptive Direct Search : Mixed-Integer NonLinear Programming MIQP : Mixed Integer Quadratic Programming f pxq : objective function with continuous variables
		f px, yq	NLP	: NonLinear Programming : objective function with mixed continuous and binary variables
		▽ x f	NOMAD : derivatives of f with respect to x : Nonlinear Optimization by Mesh-Adaptive Direct Search
		Z	QP	: Quadratic Programming : mixed variables interpolation set
		Ẑ mk	RBF ROM SA	: Radial Basis Function : normalized mixed variables interpolation set : Reduce Order Model : Simulated Annealing : mixed variables quadratic model at iteration k
		p	SCIP	: Solver for Constraint Integer Programming : cardinality of Z
		∥.∥ F		: Fröbenius norm
		∥.∥ 8		: l 8 norm
		∥.∥ p		: l p norm
		∥.∥ 2		: l 2 norm
		∆ x		: trust-region radius relative to the x search subspace
		∆ y		: trust-region radius relative to the y search subspace
		∆ x,0		: initial value of ∆ x
		∆ y,0		: initial value of ∆ y
		H k		: Hessian matrix at iteration k
		ϕ		: natural basis of polynomials
	xii			xv xvii xix

35 Zoom on the best current objective-function values versus number of evaluations of DFOb-d neck with the five options. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 The two solutions obtained with DFOb-d neck and the five options on the Ackley function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Fuel consumption for aircraft and impact of the small reduction from 2016. . . . . 3 Number of distinct solutions for n blades (from [74]). . . . . . . . . . . . . . . . . 4 Solver parameters and options used for the benchmark. . . . . . . . . . . . . . . . . 5 Analytical necklace-optimization benchmark problems . . . . . . . . . . . . . . . . 6 List of benchmark functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 indicated in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . indicated in bold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . H (Step 1) at iteration k . . . . . . . . . . . . . . . . 38 Exploration phase of DFOb-d H (Step 1.5) at iteration k . . . . . . . . . . . . . . . 40 DFOb-d H algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Greedy algorithm for continuous-variable problems . . . . . . . . . . . . . . . . . . 97 Projected DoE for discrete variables . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Greedy-MDS algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Table 2

 2

	∽ 240.10 12	∽ 83.10 9
	0.5% Ó	410.10 6 Ó

shows the fuel consumption for aircraft: if we can reduce by 0.5% the amount of fuel bust, we can gain around 410 million US $ per year. Therefore, reducing fuel consumption (by increasing engine efficiency) and maintenance costs (by decreasing vibrations) are two major concerns of the aeronautics industry.

Tab. 2.: Fuel consumption for aircraft and impact of the small reduction from 2016.

World fuel consumption (liters) Corresponding price (US $)

  15.03.2021, revision submitted on 18.06.2021, accepted on 27.08.2021), and a technical report entitled "Nonlinear optimization of mixed continuous

and discrete variables for black-box simulators", https://hal-ifp.archives-ouvertes.fr/ hal-02511841/document.

  a symmetric approximation of the Hessian of f at x k , and s P R m . The approximate solution obtained, s k , should satisfy the condition

4

  DFOb-d H : A derivative-free trust-region method for mixed binary problems Contents 1 The quadratic model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2 Initial interpolation set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3 Exploitation phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4 Exploration phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5 Local convergence results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

  Necklace context . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 1.1 Concept of "necklace" . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 1.2 Survey of necklace distances . . . . . . . . . . . . . . . . . . . . . . . . 53 2 An adapted distance for cyclic-symmetry problems . . . . . . . . . . . . . . . . 56 2.1 The necklace distance . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

	DFOb-d neck : A DFOb method for	5
	cyclic-symmetry problems
	Contents
	1
	2.2	Reformulation of the QP subproblems involving the necklace distance . 60
	2.3	Algorithm DFOb-d

neck . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.1 New mixed binary optimization test problems featuring cyclic symmetry 69 3.2 Results obtained on the 25 analytical problems . . . . . . . . . . . . . . 71 3.3 Design of compressor blades in a helicopter turbomachine . . . . . . . 74 4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

  1 , y 2 P t0, 1u n , it satisfies the following properties:

	Proposition 5.1. (Non-negativity property) d neck has the non-negativity property
	d neck py, y 1 q ě 0, @y, y 1 P t0, 1u n .

Proof. The statement is a corollary of the non-negativity property of the Hamming distance. Proposition 5.2. (Reflexivity property) d neck has the reflexivity property d neck py, yq " 0. Proof. Since d H py, yq " 0 holds true, 0 is the minimum bound of d neck py, y 1 q. Therefore d neck py, yq " 0. Proposition 5.3. (Commutativity property) d neck has the commutativity property d neck py, y 1 q " d neck py 1 , yq. Proof. Based on the fact that d H py, Rot i py 1 qq " d H py 1 , Rot n´i pyqq, thus min i d H py, Rot i py 1 qq " min i d H py 1 , Rot n´i pyqq " min i d H py 1 , Rot i pyqq. Proposition 5.4. (Triangle inequality property) d neck satisfies the triangle inequality property:

  Consider the special case ofProposition 5.6 where Ω " R m ˆt0, 1u n , z " px, yq, and N " m `n, and restrict both feasible sets of (P 1 ) and (P 2 ) by adding the two constraints: ∥x ´xk ∥ 8 ď ∆ x,k and t ď ∆ y,k .2.3 Algorithm DFOb-d neckThe algorithm we are introducing in this work to deal with derivative-free mixed binary optimization problems is a modified version of DFOb-d H in which we replace the Hamming distance, d H , with the necklace distance, d neck . To do so we use the above formulation of the no-good cut constraints as linear constraints, and the reformulated MIQP subproblem. This subsection presents the new algorithm DFOb-d neck and establishes its local convergence.

	The new algorithm is named Derivative-Free trust-region method for mixed binary necklace
	optimization, and is denoted DFOb-d neck . It follows exactly the steps of DFOb-d H (Algorithm 6,
	given at the end of Chapter 3), except for the following specific changes:
	• In Step 1.5a:
	Solve MIQP subproblem (71), instead of MIQP subproblem (31).

i"1,2,...,n tg i pz ˚qu " t, and on the other hand g I ˚p zq " g I ˚pz ˚q " min i"1,2,...,n tg i pz ˚qu by definition of I ˚. Therefore, the I ˚th constraint of (P 2 ) is satisfied for pz, ȳ, tq.

• In terms of objective-function values, it is clear that f pz ˚q `µt ˚ą f pzq `µt , since by hypothesis t ˚ą min i"1,2,...,n tg i pz ˚qu, while t " min i"1,2,...,n tg i pz ˚qu. This contradicts the optimality of pz ˚, y ˚, t ˚q.

Corollary 5.1. The two problems

[START_REF] Minghui | On the sum of distances along a circle[END_REF]

and (

71

) are equivalent.

Proof.

• In Step 3:

Replace, in the new no-good cut, the Hamming-distance inequality:

d H py, y k q ě K ˚,

by the n inequalities: d H py, Rot i py k qq ě K ˚, i " 1, 2, . . . , n, which are linear constraints equivalent to d neck py, y k q ě K ˚by

[START_REF] Meunier | Computing solutions of the paintshop-necklace problem[END_REF]

.

Let us now derive a result of local convergence for the new algorithm DFOb-d neck , analogous to that established in Lemma 4.2 and Theorem 4.1 for DFOb-d H . First, let us consider the perturbed model:

  with n v the number of variables of problem p. It is normalized by n v `1 since the number of simulations grows when the number of variables increase. Table 4 summarizes the options and the main parameter values used in our comparison for the five solvers under study: DFOb-d H , DFOb-d neck , NOMAD, RBFOpt and DFLBOX.

		n p	cardtp P P :	t p,s n v `1 ď κu,
		Maximal number of	300
		objective-function
		evaluations
	All	Initial design description	LHS design with rounding of
	solvers			discrete variables (default op-
				tion of RBFOpt)
				For NOMAD and DFLBOX: the
				initial point is the best point of
				the initial design
		Initial design size		m `n `1 points
		Trust-region radii for	∆ x,0 " 1 ∆ x " 10 ´3, ∆ x "
		continuous variables	10 2
	DFOb Trust-region radius for bi-	∆ y,0 " 2
		nary variables
		Maximal number of	minp14, 2
		no-good cuts

Tab. 4.: Solver parameters and options used for the benchmark. n ´1q, if n ď 6

  list the name and the source of 25 optimization problems that we selected to build our 25-instance benchmark problem set, more details in Appendix 2. Originally these benchmark problems involve continuous optimization variables only. Some of these optimization problems also include constraints. The present study considers mixed binary optimization problems featuring cyclic symmetry and involving only bound constraints. Thus, we describe now how we transform these instances from the literature into new mixed binary bound-constrained optimization problems featuring cyclic symmetry. A first step in this transformation process is to propose intermediate mixed categorical (involving integer variables not related to effective quantities) optimization problems.

	The 10 instances from [64] are associated with continuous-optimization minimax problems of the
	form:
	min xPrx,xs	F pxq :" max

wPt1,2,...,lu f w pxq, where l ě 2, and f 1 , f 2 , . . ., f l are given functions. We transform these instances into mixed categorical problems of the form: min xPrx,xs,wPt1,2,...,lu F px, wq :"

Contents 1

 1 Geometrical designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 2 Marginal designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 3 Low discrepancy designs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

4 Super-sample from kernel herding . . . . . . . . . . . . . . . . . . . . . . . . . 94 Design of experiments (DoE) are used in various contexts such as optimization or uncertainty quantification based on a time-consuming numerical simulator. It aims to select a limited number of values to assign to the simulator input variables that give a maximal knowledge on the simulator outputs of interest.

Table 6

 6 These results are compared with random sampling and standard LHS with rounded values for binary variables. The empirical means of the benchmark functions for each necklace are computed from the points of the DoE and compared with the empirical means computed with a very large set of points obtained by Monte-Carlo sampling (10 4 evaluations).

	with a number of continuous variables

  The aim of this chapter is to propose design of experiments methods adapted to mixed discrete and continuous variable problems. A first proposal is the straightforward extension of continuous

		Adapted-greedy	Greedy-MDS	LHS + FPCA	Quant.+FPCA
	nDoE time(s)	errors	time(s)	errors	time(s)	errors	time(s)	errors
	20	0.04	8.3.10 ´8	0.02	1.9.10 ´7	0.08	1.6.10 ´7	0.36	1.3.10 ´7
	40	0.072	9.3.10 ´8	0.027 2.5.10 ´7	0.096 2.9.10 ´8	0.94	6.8 ´9
	60	0.109 1.6.10 ´8 0.030 4.1.10 ´8	0.11	5.7.10 ´8	2.04	1.6 ´8
	6 Conclusions

  Computational times and relative errors of the expectation estimation obtained for the function f 1 with the functional data of Max-stable process. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold. Computational times and relative errors of the expectation estimation obtained for the function f 2 with the functional data of Brownian motion. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold. Computational times and relative errors of the expectation estimation obtained for the function f 2 with the functional data of Max-stable process. 4 methods are evaluated: the adapted-greedy and MDS-greedy methods coupled with the global alignment kernel, the LHS method and the quantization method coupled with the FPCA dimension reduction method. The best values among the 4 methods are indicated in bold.

	Adapted-greedy nDoE time(s) errors 20 0.036 2.3.10 ´7 0.028 Greedy-MDS time(s) errors 2.10 ´6 40 0.073 1.8.10 ´6 0.034 1.6.10 ´8 " Life is like riding a bicycle. To keep your balance you LHS + FPCA Quant.+FPCA time(s) errors time(s) errors 0.071 1.7 ´6 0.353 5.4.10 ´7 0.09 8.10 ´7 0.93 3.2 ´7 must keep moving.
	60	0.11		5.10 ´7	0.039	4.1.10 ´7	0.11	3.8.10 ´7 -Albert Eistein 1.93	2.10
	Contents								Part III
		Adapted-greedy		Greedy-MDS	LHS + FPCA	Quant.+FPCA
	nDoE time(s)	errors		time(s)	errors	time(s)	errors	time(s)	errors
	20	0.036		3.10 ´7		0.022	6.3.10 ´7	0.05	5.5.10 ´7	0.33	2.3.10 ´7
	40	0.072 2.7.10 ´7 0.027 Adapted-greedy Greedy-MDS 7.6.10 ´7 Escaping from local minima 0.07 1.8.10 ´7 0.89 6.1.10 ´8 LHS + FPCA Quant.+FPCA
	60 nDoE time(s) 0.11	3.10 ´8 errors time(s) 0.03	1.1.10 ´8 errors time(s) 0.09	1.1.10 ´7 errors time(s) 1.9	2.5.10 errors
	20	0.036 5.2.10 ´7 0.028 1.1.10 ´6	0.071 2.8.10 ´6	0.35	1.6.10 ´6
	40	0.073	1.9.10 ´6	0.034 1.9.10 ´6	0.091	1.10 ´6	0.91	8.9.10 ´7
	60	0.11	4.10 ´7		0.04	2.1.10 ´6	0.11	8.1.10 ´7	1.91	6.3.10 ´7

Tab. 8.:

´7 Tab. 9.: ´8 in global optimization methods such as surrogate optimization or direct search methods could exploit the proposed DoEs for mixed continuous and discrete problems or other type of variables provided that some prior information is available and can be associated to a kernel or a distance (as we did for cyclic symmetry). Tab. 10.: 1 Augmented DoE strategies for sampling a discrete variable space . . . . . . . . 122 2 Restart strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 2.1 Hard restart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 2.2 Complete soft restart . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 2.3 Partial soft restart . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 3.1 Results for the DS-MM function . . . . . . . . . . . . . . . . . . . . . . 128 3.2 Results for the Ackley function . . . . . . . . . . . . . . . . . . . . . . . 128 4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

  x 2 ´x1 .

	Problem B.2. CB3.
	F pxq " max

1ďiď3 f i pxq, f 1 pxq "x 4 1

  x 2 ´x1 .

	Problem B.5. PENTAGON.
	F pxq " max 1ďiď3
	Problem B.3. QL.
	F pxq " max 1ďiď3	f i pxq,
	f 1 pxq "x 2 1	`x2 2 ,
	f 2 pxq "x 2 1	`x2 2 `10p´4x 1 ´x2 `4q,
	f 3 pxq "x 2 1	`x2 2 `10p´x 1 ´2x 2 `6q.
	Problem B.4. WF.
	F pxq " max 1ďiď3	f i pxq,
	f 1 pxq "	1 2	px 1	`10x 1 x 1 `0.1	`2x 2 2 q,
	f 2 pxq "	1 2	p´x 1	`10x 1 x 1 `0.1	`2x 2 2 q,
	f 3 pxq "	1 2	px 1	´10x 1 x 1 `0.1	`2x 2 2 q.

  `x1 x 2 ´14x 1 ´16x 2 `px 3 ´10q 2 `4px 4 ´5q 2 `px 5 ´3q 2 2px `x1 x 2 ´14x 1 ´16x 2 `px 3 ´10q 2 `4px 4 ´5q 2 `px 5 ´3q 2 2px 6 ´1q 2 `5x2 7 `7px 8 ´11q 2 `2px 9 ´10q 2 `px 10 ´7q 2 `px 11 ´9q 2 10px 12 ´1q 2 `5px 13 ´7q 2 `4px 14 ´14q 2 `27px 15 ´1q 2 `x4 16 `px 17 ´2q 2 13px 18 ´2q 2 `px 19 ´3q 2 `x2 20 `95, f 2 pxq "f 1 pxq `10p3px 1 ´2q 2 `4px 2 ´3q 2 `2x 2

	2 3 `x2 `x2 4 ´5x 1 ´5x 2 ´21x 3 `7x 4 , 4 `x1 ´x2 `x3 ´x4 ´8, `x2 `x2 3 1 `2x 2 f 2 pxq "x 2 1 `x2 2 f 3 pxq "x 2 2 3 `2x 2 4 ´x1 ´x4 ´10, f 4 pxq "x 2 1 `x2 2 `x2 3 `2x 1 ´x4 ´5. Problem B.7. WONG 2. F pxq " max 1ďiď6 f i pxq, f 1 pxq "x 2 1 `x2 F pxq " max 1ďiď14 f i pxq, 2 Problem B.8. WONG 3. f 1 pxq "x 2 1 `x2 2

  2 1 `35x 15 ´79x 16 ´92q,f 11 pxq "f 1 pxq `10p15x 2 2 `11x 15 ´61x 16 ´54q, f 12 pxq "f 1 pxq `10p5x 2 1 `2x 2 `9x 4 17 ´x18 ´68q, f 13 pxq "f 1 pxq `10px 2 1 ´x9 `19x 19 ´20x 20 `19q, f 14 pxq "f 1 pxq `10p7x21 `5x 2 ´1, f 2 pxq " sinpx 1 q, f 3 pxq " ´cospx 2 q.

	Problem B.11. HS2.
	2	`x2 19 ´30x 20 q.
	Problem B.9. MAD1.
		F pxq " max 1ďiď3	f i pxq,
	f 1 pxq "x 2 1 2 `x1 x 2 Problem B.10. MAD4. `x2
		F pxq " max 1ďiď3	f i pxq,
		f 1 pxq " ´ex 1 ´x2 ,
		f 2 pxq " sinhpx 1 ´1q ´1,
		f 3 pxq " ´logpx 2 q ´1.
	Hock and Schittkowski benchmark

  F pxq " p1 `px 1 `x2 `1q 2 p19 ´14x 1 `3x 2 1 ´14x 2 `6x 1 x 2 `3x 2 2 qqp 30 `p2x 1 ´3x 2 q 2 p18 ´32x 1 `12x 2 1 `48x 2 ´36x 1 x 2 `27x 2 2 qq.

	2 1	`x4{3 1 qx 2 1 `x1 x 2 `p´4 `4x 2 1 qx 2 1 .

Problem B.16. Goldstein-Price.

The motivating application

Reduced Order Model (ROM)

Outline of the thesis

Chapter 2 Introduction to blackbox optimization methods for mixed integer variables

Chapter 3 Derivative-free trust-region methods for continuous variables problems

Derivative-free trust-region methods

2 Derivative-free trust-region methods

The quadratic model

Initial interpolation set

Exploration phase

Local convergence results

Conclusion

Chapter 5 DFOb-d neck : A DFOb method for cyclic-symmetry problems

Necklace context

An adapted distance for cyclic-symmetry problems

Numerical results

Conclusion

Low discrepancy designs

Super-sample from kernel herding

Greedy-MDS: kernel herding and multi-dimensional scaling

Adapted-Greedy: directly using a kernel defined in the mixed space

Numerical experiments

https://www.femm.info/wiki/HomePage

Conclusions

Augmented DoE strategies for sampling a discrete variable space

Restart strategies

Quadratic programming and mixed integer quadratic programming

´1q 2 `5x 2

`7px

´11q 2 `2px

´10q 2 `px

´7q 2 `45, f 2 pxq "f 1 pxq `10p3px 1 ´2q 2 `4px 2 ´3q 2 `2x 2 3 ´7x 4 ´120q, f 3 pxq "f 1 pxq `10p5x 2 1 `8x 2 `px 3 ´6q 2 ´2x 4 ´40q, f 4 pxq "f 1 pxq `10p0.5px 1 ´8q 2 `2px 2 ´4q 2 `3x 2 5 ´x6 ´30q, f 5 pxq "f 1 pxq `10px21 `2px 2 ´2q 2 ´2x 1 x 2 `14x 5 ´6x 6 q, f 6 pxq "f 1 pxq `10p´3x 1 `6x 2 `12px 9 ´8q 2 ´7x 10 q.2 Benchmark functions

Benchmark functions

Acknowledgments

I would like to take this opportunity to express my sincere gratitude and appreciation to everyone who supported me during this thesis, without them, this thesis would not be completed.

Definition 5.3. (Varshamov distance) Given two binary strings y " py 1 , . . . , y n q and y 1 " py 1 1 , . . . , y 1 n q, the Varshamov distance between y and y 1 is given by d V py, y 1 q " maxpN 01 py, y 1 q, N 10 py, y 1 qq, [START_REF] Kandasamy | High Dimensional Bayesian Optimisation and Bandits via Additive Models[END_REF] where N 01 py, y 1 q " #tpy i , y 1 i q : y i " 0, y 1 i " 1u and N 10 py, y 1 q " #tpy i , y 1 i q : y i " 1, y 1 i " 0u, where # denotes the cardinality of the set.

One has:

N 01 py, y 1 q " n ÿ i"1 p1 ´yi qy 1 i "

and

Therefore, we have: d V py, y 1 q " max `∥y∥ 1 , ∥y 1 ∥ 1 ˘´y T y 1 .

For more details, see [START_REF] Constantin | On the theory of binary asymmetric error correcting codes[END_REF]. Varshamov distance can detect cases where the number of 1s in string y is larger than the number of 1s in string y 1 .

Definition 5.4. (The swap distance) y and y 1 are as above. Assume that y and y 1 have the same number of 1, U " pu 1 , . . . , u k q, V " pv 1 , . . . , v k q where u i , v i are the indices such that y u i " 1 and

The swap distance is given by

In [START_REF] Ardila | Necklace Swap problem for rhythmic similarity measure[END_REF] the authors use the swap distance to measure the similarity of two rhythms (application in classification of types of music) or in bioinformatics where the two strings to be compared are chain polymers. Computing U, V from y, y 1 costs Opnq operations and Toussaint [START_REF] Toussaint | The Geometry of Musical Rhythm[END_REF] states that the swap distance can be computed in Opn 2 q. However, in [START_REF] Ardila | Necklace Swap problem for rhythmic similarity measure[END_REF], the authors introduce algorithms computing the swap distance in Opk 2 q and Opk 3 q operations. Definition 5.5. (The Hamming distance with shifts) Given y, y 1 two binary strings, we define three types of operations on y:

• An insertion, inspiq, changes y i from 0 to 1 with cost c ins ;

We propose to replace [START_REF] Minghui | A Linear-Time Algorithm for Hamming Distance with Shifts[END_REF]

where M is some large-enough positive constant (one can easily verify that in fact it suffices to set M to the value n `1) and ỹ is a vector of n auxiliary binary variables.

In the sequel we shall write that two optimization problems are equivalent if an optimal solution of one problem straightforwardly provides an optimal solution of the other problem, and vice versa. Proposition 5.6 below is introduced in order to show that the new problem (71) (involving only linear constraints) is equivalent to problem [START_REF] Minghui | On the sum of distances along a circle[END_REF]. Since the essential difficulty resides in the "min" constraint of ( 70) and in the Hamming-distance constraints of (71), the proposition statement disregards the trust-region constraint on x and the constraint t ď ∆ y,k (both of which are straightforwardly modeled in an MIQP). Corollary 5.1 below will establish the equivalence of problems ( 70) and ( 71), as a special case of Proposition 5.6.

Proposition 5.6. (Mini-min reformulation) Let µ ą 0 be a given constant, N, n be positive integers, and let f :

, n, be real-valued

Tab. 5.: Analytical necklace-optimization benchmark problems

Test problem m n source instance reference

Branin-nl 1 3 Branin [START_REF] Dixon | The global optimization problem: An introduction[END_REF] Camel-nl 1 3 Camel [START_REF] Dixon | The global optimization problem: An introduction[END_REF] Goldstein-Price-nl 1 3 Goldstein-Price [START_REF] Dixon | The global optimization problem: An introduction[END_REF] Hartman3

Sporttournament06-nl 14 3 Sporttournament06 [71] We first restrict the last continuous variable, say x end , to take only a finite number of values in the discretized-interval set: In this section, we focus on the specific case of models with mixed experimental region defined as D " tz " px, yq P R m ˆIn u, where x P R m , and y P I n are the continuous and discrete variables, respectively, and where I denotes the discrete space.

where S R k,i is the sigmoid function calculated by:

and

with f min , fk,i and σ k,i are respectively the current minimum objective-function value, the mean of the evaluated objective function values and the corresponding standard deviation in the level c i ; p k,i are for all levels i " 1, 2, . . . , n levels which give, for each level, the probabilities that it has not been sufficiently explored. This probability distribution p " pp k,1 , p k,2 , ..., p k,n levels q is sampled to obtain a target distribution for our DOE method: MDS greedy or adapted-greedy DOE method (see Part II), completed with the available simulated points. The greedy algorithm is applied from the given set of simulated points instead of starting with an empty set as the original proposal.

We illustrate the proposed method on the toy problem introduced in [START_REF] Zuniga | Global optimization for mixed categorical-continuous variables based on Gaussian process models with a randomized categorical space exploration step[END_REF], called DS-MM function in the sequel, with m " 1 continuous and n " 5 binary variables. Remark that the original problem is adapted to a mixed cyclic binary problem, with the following piecewise-defined function involving 8 univariate functions

We shall use this function as a test function for optimization in Section 3. As we can see in Figure 28, this function has several local minima and some of the individual one-dimension functions (corresponding to each level) are strongly correlated or anti-correlated. The global minimum is marked as a circle point in the individual function associated to level 4 (noted OF 4). Finally, we observe that four individual functions have the local minima OF1, OF4, OF7, OF8, whereas the four others (OF2, OF3, OF5, OF6) are almost flat. As a result, we would expect that the augmented can put more points in these levels.

There are several useful inequalities for matrix norms, namely

and for vector norms ∥x∥ 8 ď ∥x∥ 2 ď ? n∥x∥ 8 .

Conditioning

One important property that we use to study a matrix is its condition number which gives a quantitative information as to whether a matrix is ill-conditioned or well-conditioned. The condition number of a n ˆn non-singular matrix A is noted χpAq and is given by

which depends on the chosen norm. The value of the condition number gives an indication on the difficulty of solving numerically a linear system: the system is easier to solve if its condition number is small (ideally χpAq equal to 1) and it becomes difficult or unsolvable if the condition number is going to infinity. More precisely, a small perturbation (due for instance to the finite precision of the machine or to the imprecision/uncertainty in the data) in the right-hand side of a system of linear equations can yield an enormous change in its solution when the matrix is ill-conditioned.

where f, and g are quadratic and linear functions, respectively. To solve MIQP, the most common class of methods is Branch-and-bound algorithms and their extended versions, such as spatial branch-and-bound, branch-and-reduce, α branch-and-bound (for details see [START_REF] Burer | Non-convex mixed-integer nonlinear programming: A survey[END_REF][START_REF] Lodi | Mixed integer nonlinear programming tools: An updated practical overview[END_REF]).

[ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 7 . 0 , 9 . 0 , 5 . 0 , 1 . 0 , 2 . 0 , 3 . ] , [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 3 . 0 , 2 . 0 , 3 . 0 , 8 . 0 , 6 . 0 , 7 . ] , [ 4 . 0 , 1 . 0 , 8 . 0 , 6 . 0 , 7 . 0 , 9 . 0 , 3 . 0 , 1 . 0 , 2 . 0 , 3 . ] ] c = [ 0 . 1 , 0 . 2 , 0 . 2 , 0 . 4 , 0 . 4 , 0 . 6 , 0 . 3 , 0 . 7 , 0 . 5 , 0 . 5 ] ppi `1q `βqpx k i ´p 1 pi `1q q k q 2 `1000.

MINLPLib2 benchmark

2 Benchmark functions