text:
Elizabeth Montoya Delgadillo

Fernanda Viola
email: fviola@famaf.unc.edu.ar

Laurent Vivier
email: laurent.vivier@univ-paris-diderot.fr

Choosing a Mathematical Working Space in a modelling task: The influence of teaching

Keywords: Mathematical Working Space, paradigms, modelling task

come   L'archive ouverte pluridisciplinaire

Introduction

Nowadays, modelling is fostered in secondary school and at the university level, and one expects teachers and engineers to be trained to achieve the modelling competence. For instance, modelling constitutes a specific competence that has to be developed by civil engineers in Latin America, according to Tuning project [START_REF] Guerrero | Educación superior en América Latina: Reflexiones y perspectivas en Ingeniería Civil: Proyecto Tuning Latinoamérica[END_REF]. But, what modelling problems have to be considered in upper secondary school and in university? Are the modelling problems used in high school, and in the university, standardized by the institutions? Do they refer to real situations? We think that mathematical objects, different representations used, and properties involved in the modelling process give a richer sense for students in term of the knowledge they can construct. Teachers have of course an important role to play that highlights the importance of the teachers training. In that way, the population of this study was mainly choosen in secondary teacher training.

In this study we present a relation between the Mathematical Working Space research model [START_REF] Kuzniak | Spaces for mathematical work. Viewpoints and perspectives[END_REF][START_REF] Kuzniak | Mathematical Working Spaces in schooling: An introduction[END_REF] and the modelling cycle [START_REF] Blum | Mathematical modelling: Can it be taught and learnt?[END_REF] that could be understood as a first cycle for the resolution of the modelling task. Mathematical Working Space (MWS) is a model that is used in research in mathematics education, first developed in the field of geometry. When a student starts with the given situation, one assumes that he/she begins a horizontal mathematization process which is a foundation for bringing the situation problem into a mathematical domain. Then a vertical mathematization process takes place where the MWS framework and the modelling cycle can interact with each other. That is our approach in this study.

We first present how to use the MWS model to analyse a modelling activity. Then, we study a modelling task given to pre-service teachers in Argentina and in France. Our aim is not to make a comparative study, since both populations are different, but rather to identify the personal MWS of the students when solving this modelling task. The focus is mainly on three levels of teachers training in Argentina. We also looked at primary teachers training in France in order to have an idea of no scientific students' answers.

MWS and Modelling

We consider a MWS that depends on a specific mathematical field [START_REF] Kuzniak | Mathematical Working Spaces in schooling: An introduction[END_REF] such as, among others, analysis, geometry, algebra or statistics. Paradigms in an MWS, depending of the domain, serve to characterize the work according to a community or an institution. We develop the paradigms for analysis [START_REF] Montoya-Delgadillo | Mathematical Working Space and paradigms as an analysis tool for the teaching and learning of analysis[END_REF].

The MWS model

Three types of MWS may be distinguished: (i) MWS of reference, which is defined according to the relation to knowledge, ideally under mathematical criteria; (ii) suitable MWS, which depends on the institution involved, and is defined according to the way that this knowledge is supposed to be taught in the institution; (iii) personal MWS, which depends on the individual and is defined by the way in which the individual handles a mathematical problem with his or her own knowledge and cognitive capacities.

MWS is an environment in which reflection results from the interaction between an individual and a problem in a mathematical domain. It is an environment organized for an expert of this domain, by means of two interconnecting planes: the epistemological and the cognitive planes [START_REF] Kuzniak | Spaces for mathematical work. Viewpoints and perspectives[END_REF][START_REF] Kuzniak | Mathematical Working Spaces in schooling: An introduction[END_REF].

The epistemological plane is composed of three poles (Figure 1a), namely referential (properties, theorems, definitions…), representamen (semiotic signs), and artefacts (material or symbolic). The cognitive plane consists of the processes of visualization, construction and proof. The functioning of a MWS must not be understood as a union of single components lying on the epistemological and cognitive planes, but rather as links activated by two or three geneses, semiotic, instrumental and discursive genesis, that articulate the two planes.

Paradigms of analysis

The situation we propose (see following section) is an optimization task. Hence, using (mathematical) analysis is quite natural, even if the problem can be solved in various ways. In order to identify the paradigms of analysis in the answers of the students, we present the three working paradigms of analysis identified by Montoya Delgadillo and Vivier (2016):

-Arithmetic/Geometric Analysis (GA): it enables interpretations with implicit assumptions based on geometry, arithmetic calculations or the real world. -Calculation Analysis (CA): the rules of calculation are defined more or less explicitly and are applied independently of reflection on the existence and nature of the objects. -Real Analysis (RA): it is characterized by work involving approximation and neighbourhoods; definition and properties are set theoretically; an "ε work".

Modelling

The development of the modelling skills [START_REF] Blum | Mathematical modelling: Can it be taught and learnt?[END_REF] mobilizes notions and mathematical objects of different mathematical domains, where the knowledge that the students can learn is grounded on arguments that belong to different domains (analysis, probabilities, geometry, etc.). This gives rise to different MWS and paradigms.

Recently, two doctoral theses in probability in Paris Diderot university proposed a use of MWS in modelling processes1 . The whole modeling cycle (Figure 1b) is not taken into its totality: the focus is on phases 3 to 5 of the cycle and our aim is to analyse, with the MWS model, the mathematization process when students solve a modelling task. We expect to identify, in students' work, the solving mathematical domain and mathematical objects, representations or signs, artefacts, mathematical knowledge, and the working paradigm. We show below how the two frameworks may be used together and their complementarity.

In the following, we do not consider all the modelling process with the lens of the MWS model, but rather how is it possible to analyse mathematical activity. In particular, we set the following questions: Given a specific task, what domain and (personal) MWS will a student choose during the modelling process? Is it possible to see an influence of a suitable MWS?

Experimental study

The "gutter" situation (see below) was given to three groups of future secondary mathematics teachers in Argentina: two groups of a private training institute, namely 1PC of 1 st year and 2PC of 2 nd year, and one group of 4 th university year, 4PF, of the Universidad Nacional de Córdobaunfortunately, a 3 rd level group was not available. The study focuses on these pre-service teachers, but the task was also posed in France, in an examination of 3 rd year future primary teachers, at the Université Paris Diderot. This extra population helps to understand what kind of solving processes one can expect from non-scientific students. Students of the study are named L1 to L28, 1PC1 to 1PC24, 2PC1 to 2PC15, and 4PF1 to 4PF12 accordingly.

As announced above, it is an optimization task, quite classic except for modelling aspects. Here, we make the hypothesis that the majority of students who followed a calculus course that provided a method for solving a class of optimization problems (derivation or optimum of quadratic functions) will work in paradigm CA, and that they will resort to functions.

The aforementioned was the case for 4PF students trained at the university, since these students had courses on mathematical analysis. But CA paradigm was not expected for 1PC students and 2PC students, since in Argentina the first calculus course takes place during the second year and the experimentation was at the very beginning of the academic year. We then expected more GA methods (see below).

On the other hand, French students were in a multidisciplinary third university year after the validation of two disciplinary university years (mathematics, or biology, or history, or English…). Students' profiles were very diverse and a few of them studied mathematics, and specifically calculus.

We make the hypothesis that, spontaneously, most of L-students work in GA paradigm, whether or not using functions.

The gutter situation: a priori analysis

First, we give the statement of "the gutter", then we make an a priori analysis.

We have a rectangular metal sheet of 30 cm width and of big length. We fold up perpendicularly edges on each side to make a gutter (see dotted lines on the figure below). For obvious reasons, both side edges of the gutter will have the same size.

How should we fold up the metal sheet in order to obtain a gutter with a maximum flow?

This is a partially modelled task, since many parameters are fixed and, moreover, the geometrical mathematization is given. The choice to fold perpendicularly allows simplifying the task by eliminating other cases. There is a lack of parameters to study the flow, among them, the slope of the gutter (that we shall not discuss) and the length of the sheet.

The latter is not necessary because we can replace the study of the flow by the study of the area of a cross-section of the gutter supposing that the flow through this section is constant. Nevertheless, it is simpler to work with the volume than with the flow which is a quotient magnitude. Thus, some students may choose a length to make a calculation of volume. Is this length considered as a parameter or either as a numerical value added to the statement? This length is used in calculations or is it only a useful intermediary to think of the situation?

Let us note that there may be problems of understanding of the proposed geometrical modelling. An inadequate understanding of the situation corresponds to the addition of an idea or a belief, which is not correct regarding the situation. We expect in particular to have equitable answers with 10 cm for each of the three edges, as well as the assertion that the area, and thus the flow, does not vary because of the constant 30 cm.

On the other hand, there is no indication on the mathematization allowing to make calculations for solving. We focus our analyses on this phase of mathematization, which we interpret as the choice of a MWS by students, and the phase of solving. We then look at students' answers with an adequate understanding of the modelling situation: good geometrical shape (specified in the statement), and considering the area of a cross-section, or a volume, to study the various values (the variations) of the flow. We also look at a schema of the gutter: cross-section, in 2D or in 3D.

In the following, we present the types of answers expected for solving the modelling task. Before that, let us notice that the task can be solved in the geometry field of the statement, with magnitude and a knowledge on the areas of rectangles of given perimeter: By taking two "gutters", one forms a rectangular pipe the section of which has a perimeter of 60 cm. Since the area is maximal for a square, the solution is obtained for a 60 cm/4 side, that is, 15 cm. So, the basis has 15 cm of length and one has to fold in 15 cm/2 = 7.5 cm.

The resolution can be made in a numerical MWS, that is, by doing several calculations of the area or the volume. Several levels are possible: choice of some values of the length of the edge, calculations, and then decision-making. Numbers to be considered are essentially whole or decimal numbers. The solution being 7.5 cm, considering only integers multiples of cm does not allow to find the optimal value (unless changing the unit of length). These values can be grouped in a table or not, obtained in an organized way or not (for example, with all integers from 0 to 15, or by an oriented search), with a sign (as the letter x) to denote the length of the edge or not. A formula allows to automate calculations and, possibly, the implementation in an instrument of calculation. Obtained values can be also put in a graph.

The production of a formula can allow to change of MWS to a functional MWS, in calculus. One expects the introduction of a quadratic function2 :

 Using only algebra, with the expression of a global variational principle:

f(x)  f(a)
where a is the candidate value, that has to be find, to be the optimal value;

 With recognition of a quadratic function, its properties, vertices or symmetry axis of a parabola, allow to solve the problem. This is a work in CA paradigm.

 Using the derivative of the function. Is the theorem quoted? Is the change of sign of the derivative evoked or forgotten? This is mainly a work in CA paradigm.

 Setting a values table or a graphic, recognizing or not a parabola, in GA paradigm. The difference with numerical MWS mentioned above may be difficult to identify.

Results

Let us note that 12 L-students do not answer, and utterances from 5 students of 1PC are not classified. Three of these 1PC-students gave the solution 7.5 cm without any justification: possibly an exchange of information between students or, for at most two students, the divisions of 30 cm by 4 like for a square3 . There is also the atypical 1PC24's production, which will be mentioned later on.

In Table 1 there is a summary of the results with four groups of indicators:

1. Representamen: cross-section of the gutter (CS), letter for the edges (Ledg), letter (L) or value (V) for the fixed length of the gutter (Par), table (T) and graph (G);

2. Objects used: function (Fu) or formula (Fo);

3. Knowledge: derivation (Der), vertex of a parabola or quadratic function (Ver) in a MWS of function, CA paradigm; numerical calculations (Num) in a numerical MWS;

4. Non adequate modelling: 10 cm for each edge (10), no variation (noV) and also extreme folding (ExtFold) for "the edge must be the lower" or "the greater". There is a lot of 2D schema, close to the statement schema, but only four 3D schema. On the other hand, a cross-cutting (see Table 1) seems characteristic of 4PF-students, more advanced mathematically. It is the same with formalization, introduction of letters (variables and parameters) or using functions: the 12 4PF-students used a modelling with functions. It is also the same with the use of knowledge: 4PF-students work in CA paradigm, either by using derivation or the formula giving the vertex of a parabola (rarely justifying that it is a maximum; only two students justified it, by calculating the second derivative). There were a few students that used this knowledge in the three other groups, or in an empirical way by working in GA paradigm (by numerical calculations, more or less organized, or graphs).

Although the work seems largely guided by the knowledge the 4PF-students have, with few variations, modelling was adequate, what is not the case in other groups. For students who do not well understand the situation, either no MWS rises from modelling, or a very poor MWS appears, not allowing a substantial mathematical work. Answers of inadequate modelling are:

 Equitable length of the three sides: 10 cm each (5 students);

 There is no variation, all is constant (4 students);

 An extreme folding, the littlest lateral side possible (5 students), but sometimes with the (real) constraint that 0 is forbidden, so that the water remains in the gutter;

 L27: a pyramid (to have a pipe?) without understanding the geometry of the statement;

 1PC24 drew a roof explaining: "if there is a few water that does not go strongly", and "if there is a lot of water and with great strength".

Although the productions of the group 1PC are very poor, with not much formalization, in the groups 2PC and L we find a various types of rich productions: trying values with calculations, sometimes leading to the accurate solution when the student thinks out of the integers domain (some students remain in the field of the arithmetic of integers, sometimes with answers with two solutions, 7 and 8, and even the constraints to have a basis of even length in order to be able to divide by two), graphs and tabulations which can lead to the awareness of a symmetry (L6, L28 below).

L14 proposes a formula with a letter x for the length of the edge and l=10 for the length of the gutter, to calculate what he calls Dmax (figure 2). He uses this formula to find a values table for a whole number between 1 and 15. He writes: "we notice that the maximal flow is achieved for an edge between 7 and 8 cm". Other students stop at this stage, while L14 continues with another table for the values of x between 7 and 8 with a step of 0,1. This allows him to conclude correctly, in GA paradigm, that "to have a maximal flow, edges have to be 7.5 cm each".

Student L6, after the introduction of the function f(x) = (30-2x)x, establishes a valuable table for integers and draws the graph of f (figure 3). She concludes with visualization on the graph to determine the solution by drawing what seems to be the axis of symmetry of the curve. Finally, 2PC4 has an atypical production with a modelling by means of a paper model: he takes two paper strips, which he folds to form a pipe and gives, by a reasoning close to the solution in the MWSGeometry described above, the right solution.

Conclusion

The recognition of a function in the situation, like the introduction of a letter, is an important modelling activity, because a concept is introduced, which a priori has no relation to the situation. But more than a letterby using it the work can remain in a numerical MWS -, the introduction of a

Proceedings of CERME10

function switches the work to a MWS of functions with specific techniques (e.g. derivation) and representations (graphic, table value). However, the work done can be very different and the personal MWS involved can also be different from a student to another.

The understanding of the situation and the asked question made each student choose a MWS, that is to say, mathematical objects, knowledge, theorems, signs (letters, graphics,…), artefacts (calculator, spreadsheet,…) and the processes of visualization, construction and proving in various domainsin other words, the MWS... -. Those are complex issues for a teacher to control the various knowledge involved.

Then we wonder whether the modelling process for some students is a routine that has been "normalized" by the institution, mostly for a specific class of problems, or a mathematical activity which has been taught in their institution. It seems that for 4PF group, students' answers are standardized with only a choice between two knowledge on functions. This is a first influence of the suitable MWS that we point out. The diversity of the L-students productions reinforce this standardization and we set the problem of the possibility for teachers to be aware and open to alternative answers that their students may have.

MWS model allows us to analyse in depth the concepts and mathematical objects coming from mathematization identified in the modelling process. However, the possibility that different MWS from this process are generated makes an attractive but complex task to control for the teacher, overall for a second modelling cycle. Indeed, this study focus on a first cycle where one asks to students to solve a modelling task. Then, teachers may use a second cycle in order to reorganize and to point out knowledge that arose during the first cycle.

Figure 1a :Figure 1b :

 1a1b Figure 1a:The Mathematics Working Space, geneses and vertical planes[START_REF] Kuzniak | Spaces for mathematical work. Viewpoints and perspectives[END_REF]

Figure 2 :

 2 Figure 2 : Answer by student L14 Figure 3 : L6's graphL28 declares the same function but in a more formal way: "x(0,15) f(x)=(30-2x)x". He then calculates the values for the whole values of x that he associates two by two, which is an organization adapted to the function at stake. "We notice a symmetry in the values of f(x) when x varies between 1 and 15. We deduct that the maximum is in the middle of the values, that is: fmax(x)=f(7,5)=112,5 cm 2 ."

Table 1 :

 1 Students' answers to the "gutter task"

		CS	Ledg Par	T/G	Fu/Fo	Der Ver Num 10	noV	ExtFold
	4PF (12) 8	12	6L,0V 0T, 3G 12Fu, 0Fo 7	5	0	0	0	0
	2PC (15) 3	2	0L,1V 2T, 0G 0Fu, 1Fo	0	1	5	4	5	2
	1PC (24) 7	2	1L,2V 0T, 0G 0Fu, 0Fo	0	1	11	5	2	0
	L (16)	0	7	2L,3V 2T, 2G 4Fu,2Fo	2	0	5	2	0	4

See also the poster "Modelling tasks and mathematical work" in TSG6 of CERME 10.Thematic Working Group 06Proceedings of CERME10

Proceedings of CERME10

A quadratic function limits the technical difficulties. The task was also chosen for this reason.

That method is not expected here: since (a,b)ab is bilinear, optimizing the area of rectangles with fixed perimeter p is the same than optimizing the area of rectangles of sizes a and b with a + b + a = p.

Acknowledgments

This study was supported by mobility project of USPC and ECOS-Sud project C13H03.


Martine Labbé

Marina Leal

Justo Puerto

New models for the location of controversial facilities: A bilevel programming approach

Keywords: Bilevel optimization, locational analysis, combinatorial optimization

Motivated by recent real-life applications in Location Theory in which the location decisions generate controversy, we propose a novel bilevel location model in which, on the one hand, there is a leader that chooses among a number of fixed potential locations which ones to establish. Next, on the second hand, there is one or several followers that, once the leader location facilities have been set, chooses his location points in a continuous framework. The leader's goal is to maximize some proxy to the weighted distance to the follower's location points, while the follower(s) aim is to locate his location points as close as possible to the leader ones. We develop the bilevel location model for one follower and for any polyhedral distance, and we extend it for several followers and any p -norm, p ∈ Q, p ≥ 1. We prove the NP-hardness of the problem and propose different mixed integer linear programming formulations. Moreover, we develop alternative Benders decomposition algorithms for the problem. Finally, we report some computational results comparing the formulations and the Benders decompositions on a set of instances.

Introduction

Location is a research area devoted to the optimal placement of facilities [START_REF] Albareda-Sambola | The reliable p-median problem with at-facility service[END_REF][START_REF] Boonmee | Facility location optimization model for emergency humanitarian logistics[END_REF][START_REF] Hammad | A Bi-level Mixed Integer Programming Model to Solve the Multi-Servicing Facility Location Problem, Minimising Negative Impacts Due to an Existing Semi-Obnoxious Facility[END_REF][START_REF] Kalcsics | The multicriteria p-facility median location problem on networks[END_REF][START_REF] Labbé | Location on networks[END_REF][START_REF] Nickel | Location theory: a unified approach[END_REF][START_REF] Owen | Strategic facility location: A review[END_REF], including among many others emergency systems [START_REF] Bélanger | Recent optimization models and trends in location, relocation, and dispatching of emergency medical vehicles[END_REF][START_REF] Boonmee | Facility location optimization model for emergency humanitarian logistics[END_REF][START_REF] Schmid | Ambulance location and relocation problems with time-dependent travel times[END_REF], service providers [START_REF] Albareda-Sambola | The multi-period incremental service facility location problem[END_REF][START_REF] Berman | Generalized coverage: New developments in covering location models[END_REF], infrastructures, etc., and it is a basic building block of most transportation, communication or logistic problems. An optimal location can be chosen according to different criteria depending on the rationale behind the considered model. The most popular ones are the minimization of the total or maximum transportation cost [START_REF] Albareda-Sambola | The reliable p-median problem with at-facility service[END_REF][START_REF] Boonmee | Facility location optimization model for emergency humanitarian logistics[END_REF], the maximization of some coverage goal [START_REF] Albareda-Sambola | The multi-period incremental service facility location problem[END_REF][START_REF] Bélanger | Recent optimization models and trends in location, relocation, and dispatching of emergency medical vehicles[END_REF][START_REF] Berman | Generalized coverage: New developments in covering location models[END_REF][START_REF] Boonmee | Facility location optimization model for emergency humanitarian logistics[END_REF], or the minimization of the undesirable effects induced by the facilities [START_REF] Brimberg | A bi-criteria model for locating a semi-desirable facility in the plane[END_REF][START_REF] Church | Protecting critical assets: the r-interdiction median problem with fortification[END_REF][START_REF] Hammad | A Bi-level Mixed Integer Programming Model to Solve the Multi-Servicing Facility Location Problem, Minimising Negative Impacts Due to an Existing Semi-Obnoxious Facility[END_REF].

Location Theory includes a number of real-life applications in which the location decisions generate controversy. This controversy must be understood as a disagreement among users with different, non-aligned or opposite interests. Examples of this controversial location can be found in the literature, for example, in the areas of semiobnoxious facility location or in problems that involve the location and protection of critical infrastructures or facilities sensitive to intentional attacks.

In the last decades, the consciousness-raising in environmental issues has grown substantially, specially in those aspects that affect human health or have adverse effects on people. As a consequence of this awareness-raising, the location of (semi-) obnoxious facilities has been extensively studied. Obnoxious facilities are those that generate a disservice to the people nearby while producing an intended product or service [START_REF] Erkut | Analytical models for locating undesirable facilities[END_REF]. However, if only these undesirable effects are taken into account when locating them, these facilities would never be opened or would be located too far from the population centers making use of the produced services, thus generating huge costs. For that reason, in the last years, there has been an increasing focus in analyzing the problem of locating semi-obnoxious facilities [START_REF] Brimberg | A bi-criteria model for locating a semi-desirable facility in the plane[END_REF][START_REF] Hammad | A Bi-level Mixed Integer Programming Model to Solve the Multi-Servicing Facility Location Problem, Minimising Negative Impacts Due to an Existing Semi-Obnoxious Facility[END_REF][START_REF] Melachrinoudis | Semi-obnoxious single facility location in Euclidean space[END_REF]. Semiobnoxious facilities has been defined as useful but unwelcome facilities that produce environmental concerns. That is, facilities that population centers (demand points) want them away, but there are some interests (political, economical ...) in locating them close the demand points, generating in this way, location controversy. Classical examples of this kind of facilities are chemical and power plants, airports, waste dumps, detoxification centers, etc., as listed in [START_REF] Melachrinoudis | Semi-obnoxious single facility location in Euclidean space[END_REF].

Another area that has also attracted increasing attention of researchers in the last years is the location and protection of vulnerable facilities (with high risk of disruption) and the protection of critical facilities, including not only those related to disruptions produced by natural disasters or natural failures, but also those referred to disruptions produced by man-made attacks [START_REF] Church | Protecting critical assets: the r-interdiction median problem with fortification[END_REF][START_REF] Aksen | A bilevel fixed charge location model for facilities under imminent attack[END_REF][START_REF] Scaparra | A bilevel mixed-integer program for critical infrastructure protection planning[END_REF]. Critical infrastructure is a term to describe assets that are essential for the functioning of a society and economy.

Most commonly associated with the term are facilities for heating, water supply, public health, security services, telecommunication, economic sector, etcetera. Clearly, the location and protection of these types of facilities generates also controversy, since the population is aware of the risks that may appear close to those facilities due to the confrontation between two antagonist parties: attackers and defenders with visibly opposite goals.

The above-mentioned problems have been usually addressed via biobjective (multiobjective) approaches, difference of objective functions, maximin optimization and, if there exists a hierarchical structure in the decision-making process, by means of bilevel optimization.

Motivated by the increasing interest in real-world applications generating location controversy, we introduce a new model for its study and analysis.

The situation that we want to address models the existence of two parties acting sequentially in a decision-making process. On the one hand, there is a leader who wants to locate some primary facilities and must choose among a number of fixed potential locations where to establish them. On the other hand, there is one (or several) follower(s) that, once the primary facilities have been set, chooses the placement of some secondary facilities, in a continuous environment. The leader's goal is to maximize some proxy of the overall weighted distance to the follower's secondary facility locations. Meanwhile, the follower(s) aim is to locate their secondary facilities as close as possible to the primary ones, minimizing a cost proportional to the distance from the secondary facilities to the primary ones set by the leader.

The reader may observe that this model fits perfectly to the cases mentioned above. The paper is structured as follows. The following section introduces the model, sets the notation and proves the NP-hardness of the considered problem. In Section 3, we develop the mathematical programming formulations and resolution algorithms for the problem with one follower and any block norm. Two different approaches, based on the representation of the norms, have been considered. Furthermore, due to their importance, they have been applied to the case of the 1 and ∞ norm. In the next Section 4, we extend the model to several followers and non-polyhedral norms. Section 5 is devoted to the computational study of the different methods discussed in the previous sections. Finally, Section 6 concludes the paper.

The model

We consider a situation with two different types of location entities: the primary facilities (critical infrastructures, goods to protect, demand-points, etc.), and the secondary facilities (terrorists nets, thefts, semi-obnoxious facilities, detoxification centers, recycling or power plants, etc.). The primary facilities wish to be located as far as possible from the secondary facilities, meanwhile the secondary facilities aim to be located as close as possible to the primary ones. The model we present consists in choosing the location of the primary facilities (these are set first), taking into account that, the secondary facilities will be located afterwards knowing their location. For the ease of presentation, we restrict ourselves to the case where a unique secondary facility will be located. The reader is referred to Section 4 for the extension to several facilities.

We will model this hierarchical structure using Bilevel Optimization. Bilevel programming targets hierarchical optimization problems in which part of the constraints translate the fact that some of the variables constitute an optimal solution of another optimization problem. There exist a leader that acts first, and then a follower that reacts to the leader's decision.

We assume that there is a leader (setting the primary facilities) that chooses among a set of potential locations B the placement of some new primary facilities. We also consider that there is a set N B of primary facilities already established, and there exists a budget constraint on the overall investment for the location of the new primary facilities. On the other hand, once the primary facilities are set, the follower chooses the location of the secondary facility in a continuous framework. The proximity between the primary and secondary facilities is measured as a weighted sum of a distance to all primary facilities.

We denote by c j the cost of opening the primary facility j, for all j ∈ B, by C the maximum budget, by f j ∈ R n the given location point j ∈ B ∪ N B, and by w j the weight factor that scales the distance from the secondary facility to f j according to its importance. We define the binary decision variables y j = 1 if f j , j ∈ B, is open, and y j = 0 otherwise.

For the follower problem we define the decision variable x ∈ R n that specifies the location of the secondary facility.

Therefore, the bilevel problem can be modeled as:

max j∈B w j d(x, f j )y j + j∈N B w j d(x, f j ) (BLP) st. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1} j ∈ B, (2)
x ∈ arg min

x j∈B

w j d(x, f j )y j + j∈N B w j d(x, f j ),
where d(x, f j ) denotes any distance induced by some norm:

d(x, f j ) = x -f j .
Observe that BLP is a bilevel max-min problem, in which the objective function of both levels is a proxy of the distance between the primary and the secondary facilities. The resulting bilevel problem contains a knapsack problem at the upper level, and a continuous single-facility location problem at the lower level.

To state the complexity of the problem, we provide the following result.

Theorem 1. The bilevel location model BLP is NP-hard.

Proof. Let us consider the distance induced by a norm • , and an instance such that |N B| = 1 and w j 0 > j∈B w j for j 0 ∈ N B.

For this instance, we know, using the majority theorem, see for example [START_REF] Love | Facilities location: models and methods[END_REF], that the optimal solution of the continuous location problem is x * = f 0 j . Then, if we denote r * j = f 0 j -f j , Problem BLP can be written as:

max j∈B w j y j r * j s.t. j∈B c j y j ≤ C, y j ∈ {0, 1} ∀j ∈ B,
which is a knapsack problem, known to be NP-hard [START_REF] Garey | Computers and intractability. A guide to the theory of NP-Completeness[END_REF].

Mathematical programming formulations and resolution algorithms

This section is devoted to present useful mathematical programming formulations for BLP in order to solve it with off-the-self solvers. In addition, we will present alternative add-hoc algorithms, based on decompositions, that prove to be more efficient than the solvers acting on the Mixed Integer Linear Programming (MILP) above-mentioned formulations.

For the sake of presentation, we assume in this section that we measure distances via block norms. The family of block norms, also called polyhedral norms, includes all norms ||.|| P whose unit ball P is a symmetric with respect to the origin, convex bounded polyhedral set containing the origin in its interior. We will denote then by ext(P ) the set of extreme points of P and by P o : the polar set of P , that is, P o := {x ∈ R n : x, p ≤ 1, ∀p ∈ P }. The reader may note that the commonly used 1 and ∞ norms belong to this family.

In order to deal with the problem we develop two different procedures: the first one is based on the evaluation of the norm through its primal expression, using its unit ball defined by P , and the second one evaluates the norm through its dual expression, using its dual unit ball P o .

These two different forms used to handle the problem are justified, as we will see, by the fact that depending on the cases one can be more efficient than the other due to the structure of the set of extreme points of P and P o . We will illustrate this behavior in the following sections with the 1 and ∞ norms.

3.1. First approach: Evaluating norms with its primal expression. It is well-known (see for example [START_REF] Nickel | A unified approach to location analysis: Ordered Median Problems[END_REF] and [START_REF] Rockafellar | Convex analysis[END_REF]), that the value of ||x|| P is given as:

||x|| P = min b∈ext(P ) µ b , (PrimalNormP) s.t. x = b∈ext(P ) µ b b, µ b ≥ 0, b ∈ ext(P ).
This representation of the norm gives rise to a trilevel problem. Using such representation we develop a MILP Formulation and a Benders like algorithm in order to solve BLP, using off-the-shelf solvers.

3.1.1. A MILP formulation. Let us assume that x = (x 1 , . . . , x n ), f j = (f j1 , . . . , f jn ), for all j ∈ B ∪ N B and b = (b 1 , . . . , b n ) for all b ∈ ext(P ). By representing ||x|| P as in PrimalNormP, BLP can be written as the following bilevel problem:

max j∈B w j y j r j + j∈N B w j r j (BLP-P) s.t. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1} j ∈ B, (2) r ∈ arg min x,r j∈B w j y j r j + j∈N B w j r j , (3)
r j = b∈ext(P ) µ j b , j ∈ B ∪ N B, (4)
x i = b∈ext(P ) µ j b b i + f ji , j ∈ B ∪ N B, i = 1, ..., n, (5)
µ j b ≥ 0, b ∈ ext(P ), j ∈ B ∪ N B, (6)
r j ≥ 0, j ∈ B ∪ N B, (7)
x i ∈ R, i = 1, ..., n. (8)
In the above formulation, variables r j , defined in constraint (4), represent the distance between x and f j and constraints (5) set the correct representation of coordinates of the secondary facility in terms of the reference system induced by f j and the extreme points of P . Constraints ( 1) and ( 2), as in BLP, are knapsack constraints, corresponding to the choice of the location of the primary facilities, according to a budget constraint. Constraints (3)-( 8) define the lower level problem, the continuous location problem, in which the representation of the norm has been included.

Proposition 1. Problem BLP-P can be reformulated as the following single level problem BLP-P'. max j∈B w j y j r j + j∈N B w j r j (BLP-P')

s.t. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1} j ∈ B, (2)
j∈B w j y j r j + j∈N B w j r j = n i=1 j∈B∪N B β ji f ji (9) r j = b∈ext(P ) µ j b , j ∈ B ∪ N B, (4)
x i = b∈ext(P ) µ j b b i + f ji , j ∈ B ∪ N B, i = 1, ..., n, (5)
µ j b ≥ 0, b ∈ ext(P ), j ∈ B ∪ N B, (6)
r j ≥ 0, j ∈ B ∪ N B, (7)
x i ∈ R, i = 1, ..., n, (8)
α j ≤ w j y j , j ∈ B, (10)
α j ≤ w j , j ∈ N B, ( 11
) j∈B∪N B β ji = 0, i = 1, ..., n, (12)
-α j - n i=1 b i β ji ≤ 0, j ∈ B ∪ N B, b ∈ ext(P ). ( 13
)
Proof. Given a solution y, representing a feasible set of locations for the primary facilities in BLP-P, the inner problem in BLP-P is a feasible Linear Program (LP) with finite solution, and its dual is:

max n i=1 j∈B∪N B β ji f ji (Dual-P) s.t. α j ≤ w j y j , j ∈ B, (10)
α j ≤ w j , j ∈ N B, (11) j∈B∪N B β ji = 0, i = 1, ..., n, (12)
-α j - n i=1 b i β ji ≤ 0, j ∈ B ∪ N B, b ∈ ext(P ). ( 13
)
Then, Problem BLP-P is equivalent to the single level formulation BLP-P' since constraint ( 9) is the strong duality condition stating that the primal and dual objectives of the lower level problem must be equal, and the blocks of constraints ( 1)-( 2), ( 4)-( 8) and ( 10)-( 13) represent, respectively, the upper level problem constraints, the lower level primal problem constraints and the lower level dual problem constraints.

We can observe that the above formulation contains some bilinear terms: r j y j . In order to transform that formulation into a mixed integer linear problem, the bilinear terms can be linearized (see [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I Convex underestimating problems[END_REF]) giving rise to an exact MILP formulation for the bilevel problem. To this end, we substitute the terms r j y j by the variables rj ; ∀j ∈ B ∪ N B and add the following set of constraints:

(14) rj ≤ r j , j ∈ B ∪ N B, rj ≤ M j y j , j ∈ B ∪ N B rj ≥ r j -M j (1 -y j ), j ∈ B ∪ N B rj ≥ 0, j ∈ B ∪ N B.
The previous block of constraints requires to set a valid value for the "big-M "constants. It is easy to observe that M j can be chosen equal to the maximum distance between f j and any other point in B ∪ N B.

3.1.2.

Benders like algorithm for solving BLP. Now, we propose an alternative method to solve the bilevel location problem under a block norm which is based on a decomposition of the problem.

For a given solution y, the inner problem in BLP-P is an LP whose set of constraints does not depend on the variables associated to the master (leader) problem (does not depend on y). Then, if we denote by P the set of extreme points of the inner problem, solving such problem is equivalent to evaluate the objective function at the points in P and to take the minimum objective function value. Then, the continuous location inner problem can be rewritten as the following optimization problem:

max q s.t. q ≤ j∈B n i=1 w j y j r τ j + j∈N B n i=1 w j r τ j , ∀ r τ ∈ P.
In order to apply Benders decomposition, and using the above formulation, Problem BLP-P can be reformulated as:

max q s.t. j∈B c j y j ≤ C, y j ∈ {0, 1}, j ∈ B, q ≤ j∈B n i=1 w j y j r τ j + j∈N B n i=1
w j r τ j , ∀ r τ ∈ P.

Our approach to solve the above problem is to sequentially identify and add extreme points in P to the problem until a certificate of optimality is fulfilled (eventually in the worse case after adding all extreme points).

To describe the algorithm, we denote by P a subset of points in P. With the purpose of obtaining upper bounds for BLP-P, in the algorithm, we define the following Master Problem:

max q (MP) s.t. q ≤ j∈B n i=1 w j y j r τ j + j∈N B n i=1 w j r τ j , r τ ∈ P, j∈B c j y j ≤ C, y j ∈ {0, 1} ∀j ∈ B.
Lower bounds for the Problem BLP-P are obtained in the algorithm by solving the following subproblem:

q(ȳ) = min j∈B w j ȳj r j + j∈N B
w j r j , (PP-P (ȳ))

r j = b∈ext(P ) µ j b , j ∈ B ∪ N B, x i = b∈ext(P ) µ j b b i + f ji , j ∈ B ∪ N B, i = 1, ..., n, µ j b ≥ 0, b ∈ ext(P ), j ∈ B ∪ N B, r j ≥ 0, j ∈ B ∪ N B, x i ∈ R, i = 1, ..., n.
If r is an optimal solution of the above problem for a given solution ȳ feasible to the master problem MP, the inequality q ≤ j∈B n i=1 w j y j rj + j∈N B n i=1 w j rj either generates a new lower bound for MP or, if the optimal solution coincides with the previous one, it is a certificate of optimality. Based on this recursion, we propose the following algorithm:

Algorithm 1: Benders decomposition Algorithm

Initialization: Choose a solution y 0 satisfying the knapsack constraint, and solve the problem PP-P (ȳ) for ȳ = y 0 . Let r 0 be an optimal solution for PP-P (ȳ). Take P = {0} and go to iteration ν = 1.

Iteration ν = 1, 2, . . .: Solve the Master Problem MP. Let y * be an optimal solution of such problem and q * the corresponding optimal value.

• Solve PP-P (ȳ) for ȳ = y * .If q * = q(y * ). END.

• Otherwise, let r * be an optimal solution of PP-P (ȳ) . Take r ν = r * , P := P ∪ {ν}, and go to iteration ν := ν + 1.

3.1.3. The case of the 1 -norm. In this section, we apply the above reasoning to the particular important case of problem BLP-P under the rectangular distance, that is, the distance induced by the 1 -norm. We take advantage of some specific properties of this norm to exploit further its algorithmic implications. As before, n denotes the dimension of the space. The set of extreme points of the unit ball of the 1 norm is ext(P ) = {e 1 , ..., e n , -e 1 , ..., -e n }, where e i the i-th vector of the canonical basis. Further, the 1 -norm of a vector x is given by x 1 = n i=1 |x i |. By introducing variables r ji representing the non linear terms |x i -f ji |, we adapt BLP-P to the 1 -norm case.

r ji ≥ x i -f ji , j ∈ B ∪ N B, i = 1, ..., n, (16)
r ji ≥ f ji -x i , j ∈ B ∪ N B, i = 1, ..., n, (17)
As in Subsection 3.1.1, we can derive a MILP by using the primal dual optimality conditions and then linearizing the bilinear terms y j r ji by introducing new variables rji . In this formulation, dual variables α ji correspond to contraints [START_REF] Luenberger | Linear and nonlinear programming[END_REF]. The dual variables associated to constraints [START_REF] Mccormick | Computability of global solutions to factorable nonconvex programs: Part I Convex underestimating problems[END_REF]

s.t. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1}, j ∈ B, (2)
j∈B n i=1 w j rji + j∈N B n i=1 w j r ji = j∈B∪N B n i=1 -f ji α ji + + j∈B n i=1 f ji (w j y j -α ji ) + j∈N B n i=1 f ji (w j -α ji ), ( 18
)
r ji ≥ x i -f ji , j ∈ B ∪ N B, i = 1, ..., n, (19)
r ji ≥ f ji -x i , j ∈ B ∪ N B, i = 1, ..., n, (20)
x i ∈ R, i = 1, ..., n, (8) rji
≤ M ji y j , j ∈ B, i = 1, ..., n, (21) rji
≤ r ji , j ∈ B, i = 1, ..., n, (22)
rji ≥ r ji -(1 -y j )M ji , j ∈ B, i = 1, ..., n, [START_REF] Scaparra | A bilevel mixed-integer program for critical infrastructure protection planning[END_REF] rji ≥ 0, j ∈ B, i = 1, ..., n, [START_REF] Schmid | Ambulance location and relocation problems with time-dependent travel times[END_REF] α ji ≤ w j y j , j ∈ B, i = 1, ..., n, (25)

α ji ≤ w j , j ∈ B, i = 1, ..., n, (26)
α ji ≥ 0, j ∈ B ∪ N B, i = 1, ..., n, (27) j∈B
(-2α ji + w j y j ) + j∈N B (-2α ji + w j ) = 0, i = 1, ..., n. (28)
The reader can observe that valid big-M constant in this formulation are M ji = max k∈B∪N B |f ki -f ji |, for all i = 1, . . . , n and j ∈ B.

An alternative formulation can be derived for Problem BLP-1 by using the fact that the inner location problem can be decomposed into n independent linear programs, one for each coordinate. Using the optimality conditions for each such problem and the linearization technique described above, we obtain a formulation BLP-1 -2 identical to BLP-1 -1 except that constraint [START_REF] Melachrinoudis | Semi-obnoxious single facility location in Euclidean space[END_REF] is replaced by the group of contraints:

j∈B w j rji + j∈N B w j r ji = j∈B∪N B -f ji α ji + j∈B f ji (w j y j -α ji )+ j∈N B f ji (w j - α ji ), i = 1, ..., n.
Algorithm 1 can also be adapted to the case of the 1 -norm. Then, q(ȳ) is obtained by solving the lower level problem defined by ( 15) -( 17) but it can be solved in O(n|B ∪ N B|) time since, for each coordinate, it amounts to find the median of a discrete distribution.

We can also use the separability property in the proposed Benders Algorithm, by solving in each iteration n subproblems q i (y) (one for each coordinate), and considering the following Master Problem:

(MP-1 -i) max n i=1 q i s.t. q i ≤ j∈B w j y j z τ ij + j∈N B w j z τ ij ∀τ ∈ P, ∀i = 1, ..., n, j∈B c j y j ≤ C, y j ∈ {0, 1} ∀j ∈ B,
We will compare the performance of the four approaches in the computational study presented in Section 5.

3.2.

Second approach: Evaluating the norm by its dual expression. Since the polar set of a polyhedron is a polyhedron, P o induces the so-called dual norm of || • || P that can also be used to evaluate || • || P . In this case, || • || P is the optimal solution of the following linear program (see for example [START_REF] Nickel | A unified approach to location analysis: Ordered Median Problems[END_REF] or [START_REF] Rockafellar | Convex analysis[END_REF]):

||x|| P = min r (NormP 0 ) s.t. n i=1 u i x i ≤ r, u ∈ ext(P 0 )
Depending on the number and structure of the set of extreme points of P and P o , it may be more convenient to compute || • || P , by using its primal or dual expression. Fruther, this dual representation leads to different MILP formulations and the Benders approach can also be adapted.

MILP formulation:

Using the dual representation of the norm, BLP can be formulated as:

max j∈B w j y j r j + j∈N B w j r j (BLP-P 0 ) s.t. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1} j ∈ B, (2) r ∈ arg min x,r j∈B w j y j r j + j∈N B w j r j , (9)
r j ≥ n i=1 u i (x i -f ji ), u ∈ ext(P 0 ), j ∈ B ∪ N B, (29)
r j ≥ 0, j ∈ B ∪ N B, (7)
x i ∈ R, i = 1, ..., n, (8)
where variables r j , defined in constraint (29), represent the distance between x and f j . Constraints (1) and (2) relate the choice of the location of the primary facilities, according to a budget constraint, and constraints ( 9), (29), ( 7) and ( 8) define the inner subproblem, in which the representation of the norm via its dual expression has been included. Proposition 2. Problem BLP-P 0 can be formulated as the following single level problem BLP-P 0 '. max j∈B w j y j r j + j∈N B w j r j (BLP-P 0 ')

s.t. j∈B c j y j ≤ C, y j ∈ {0, 1} j ∈ B, j∈B w j y j r j + j∈N B w j r j = u∈ext(P 0 ) j∈B∪N B n i=1 -f ji u i γ uj , r j ≥ n i=1 u i (x i -f ji ), u ∈ ext(P 0 ), j ∈ B ∪ N B, r j ≥ 0, j ∈ B ∪ N B, x i ∈ R, i = 1, ..., n, u∈ext(P 0 ) γ uj ≤ w j y j , j ∈ B, u∈ext(P 0 ) γ uj ≤ w j , j ∈ N B, u∈ext(P 0 ) j∈B∪N B (-u i )γ uj = 0, i = 1, ...n, γ uj ≥ 0 u ∈ ext(P 0 ), j ∈ B ∪ N B
The proof of this proposition follows the same lines as that of Proposition 1 and is thus omitted.

We can observe that in the above formulation there appear the same bilinear terms that we have already obtained in Section 3.1.1. Therefore, the same linearization [START_REF] Labbé | Location on networks[END_REF] can be applied to obtain the corresponding MILP formulation.

3.2.2.

Benders like algorithm for solving BLP-P 0 . The Benders Algorithm proposed in Section 3.1.2 can also be applied when the norm is induced by the polar polyhedron, with the same Master Problem MP and the following primal problem:

min q(y) = j∈B w j y j r j + j∈N B w j r j , (PP-P 0 ) s.t. r j ≥ n i=1 u ki (x i -f ji ), u ∈ ext(P 0 ), j ∈ B ∪ N B, r ≥ 0, x ∈ R n , x i ∈ R, i = 1, ..., n.

3.2.3.

The case of the ∞ -norm. This is Section, we apply the above results to the important case of the infinity norm. The set of extreme points of the infinity norm is ext(P ) = (a 1 , ..., a n ) ∈ R n : a i ∈ {1, -1}, i = 1, ..n, , so that |ext(P )| = 2 n . Then formulation BLP-P' would include 2 n (|B| + |N B|) µ j e variables, and more than 2 n (|B| + |N B|) constraints. However, the number of extreme points of the polar polyhedron is much smaller: ext(P 0 ) = {e 1 , ..., e n , -e 1 , ..., -e n } and |ext(P 0 )| = 2n. Further, the ∞ -norm of a vector x is given by x ∞ = max i=1,...,n |x i |. This allows to adapt BLP-P to the ∞ -norm case as follows:

max j∈B w j y j r j + j∈N B w j r j (BLP-∞ ) s.t. j∈B c j y j ≤ C, y j ∈ {0, 1}, j ∈ B, x ∈ arg min x j∈B w j y j r j + j∈N B w j r j , r j ≥ x i -f ji , j ∈ B ∪ N B, i = 1, ..., n, r j ≥ f ji -x i , j ∈ B ∪ N B, i = 1, ..., n.
Further, a MILP formulation can also be derived.

max j∈B w j rj + j∈N B w j r j (BLP-∞ -1) s.t. j∈B c j y j ≤ C, y j ∈ {0, 1}, j ∈ B, j∈B w j rj + j∈N B w j r j = n i=1 j∈B∪N B f ji (-γ e i j + γ -e i j ), r j ≥ x i -f ji , j ∈ B ∪ N B, i = 1, ..., n, r j ≥ f ji -x i , j ∈ B ∪ N B, i = 1, ..., n, r j ≥ 0, j ∈ B ∪ N B, n i=1
(γ e i j + γ -e i j )≤w j y j , j ∈ B, n i=1

(γ e i j + γ -e i j )≤w j , j ∈ N B, j∈B∪N B (-γ e i j + γ -e i j ) = 0, i = 1, ..., n,

γ e i j ≥ 0, j ∈ B ∪ N B, i = 1, ..., n, γ -e i j ≥ 0, j ∈ B ∪ N B, i = 1, ..., n, rj ≤ M y j , j ∈ B rj ≤ r j , j ∈ B, rj ≥ r j -(1 -y j )M, j ∈ B, rj ≥ 0, j ∈ B.
Finally, The proposed Benders algorithm can also be applied to the problem under the ∞ norm. The resulting inner subproblem is given by the lower level problem of BLP-∞ .

Extensions

This section is devoted to present extensions of the problem BLP to several secondary facilities and non-polyhedral norms. We analyze the problem with K > 1 secondary facilities which means to locate K new facilities also in the lower level problem. Moreover, we extend the problem BLP to deal with norms τ for τ ∈ Q, τ ≥ 1 which requires to apply conic programming and conic duality to obtain results similar to the ones presented in previous sections. [START_REF] Bélanger | Recent optimization models and trends in location, relocation, and dispatching of emergency medical vehicles[END_REF].1. The model with K secondary facilities (independent followers) . We are interested to incorporate to the problem BLP the possibility to locate several secondary facilities rather than only one, and the goal of each secondary facility is to minimize the overall distance to the primary facilities. In the following, we analyze problem BLP with K secondary points to be located in the lower level problem, that is, we consider that instead of locating one secondary facility, K of these points must be located. For this extension we assume that we are given vectors of weights w k ∈ R n + , for k = 1, . . . , K, and we define K vectors of decision variables x k ∈ R n , for k = 1, . . . , K; where x k are the coordinates of the location of the k-th secondary point. With this notation, the new problem can be written as:

max k∈K j∈B w k j d(x k , f j )y j + j∈N B w k j d(x k , f j ) (BLP-K) s.t. j∈B c j y j ≤ C, y j ∈ {0, 1}, j ∈ B, x k ∈ arg min x k j∈B w k j d(x k , f j )y j + j∈N B w k j d(x k , f j ) ∀k = 1, ..., K.
In the particular case in which w 1 = w 2 = ... = w k , we observe that by symmetry, there is an optimal solution where the secondary facilities co-locate.

Coming back to the general problem BLP-K, the evaluation of the norm can be done via the primal or dual expression. In both cases, in order to develop a MILP formulation for the model with K secondary facilities, we can apply the same technique that in the previous section. Given a solution y of the upper level problem, the continuous location problem of each follower is linear and thus, the strong duality theorem can be applied as before. This implies that K different one-secondary facility problems are added to the leader problem. In conclusion, the same approach used with the one-secondary facility location problem is replicated K times and the same results follow.

Furthermore, the Benders algorithm can also be extended to the case with K followers. The Master Problem for this extension must be slightly modified:

max k∈K q k (MP-K) s.t. q k ≤ j∈B n i=1 w k j y j r kτ j + j∈N B n i=1 w k j r kτ j ∀τ ∈ P, ∀k ∈ K, j∈B c j y j ≤ C, y j ∈ {0, 1} ∀j ∈ B.
In addition, in this formulation, there are K primal subproblems with the same structure but with different set of w weights. Therefore, in each iteration of this Benders approach, K primal subproblems must be solved.

4.2.

The problem under the τ -norm . This section extends the analysis of the problem to the case where the inner subproblem measures distances with τ -norms

with τ ∈ Q, τ ≥ 1. Recall that ||x|| τ = ( n i=1 |x i | τ ) 1/τ . The problem to be considered in this case is max j∈B w j ||x -f j || τ y j + j∈N B w j ||x -f j || τ (BLP-τ ) st. j∈B c j y j ≤ C, (1)
y j ∈ {0, 1} j ∈ B, (2)
x ∈ arg min

x j∈B

w j ||x -f j || τ y j + j∈N B w j ||x -f j || τ , Let ρ ∈ Q be such that 1/τ + 1/ρ = 1.
In order to reformulate Problem BLPτ we give the following proposition, using the representation given in Blanco et al [START_REF] Blanco | Revisiting several problems and algorithms in continuous location with l τ norms[END_REF]. Proposition 3. The problem BLPτ can be reformulated as a the single level mixed integer conic program.

max j∈B w j ||x -f j || τ y j + j∈N B w j ||x -f j || τ (ConicP) s.t. j∈B c j y j ≤ C, y j ∈ {0, 1} j ∈ B, j∈B w j ||x -f j || τ y j + j∈N B w j ||x -f j || τ = j∈B∪N B n i=1 V ji f ji , x + -x --Z j = f j , ∀ j ∈ B ∪ N B, ||Z j || τ ≤ r j , ∀ j ∈ B ∪ N B, x + , x -∈ R n + , Z j ∈ R n , r ∈ R n , j∈B∪N B V ji + λ 1 i = 0, ∀ i = 1, . . . , n, - j∈B∪N B V ji + λ 2 i = 0, ∀ i = 1, . . . , n, -V ji + µ ji = 0, ∀ j ∈ B ∪ N B, i = 1, . . . , n, ||µ j || ρ ≤ γ j , ∀ j ∈ B ∪ N B, γ j = w j y j , if i ∈ B, γ j = w j , if i ∈ N B, λ 1 , λ 2 ∈ R n + , µ j ∈ R n , ∀j ∈ B ∪ N B, γ ∈ R n .
Proof. We observe that the inner location problem can be formulated as the following conic linear program in standard form [START_REF] Luenberger | Linear and nonlinear programming[END_REF]:

min j∈B w j y j r j + j∈N B w j r j s.t. x + -x --Z j = f j , ∀ j ∈ B ∪ N B, ||Z j || τ ≤ r j , ∀ j ∈ B ∪ N B, x + , x -∈ R n + , Z j ∈ R n , r ∈ R n .
Therefore, its conic dual can be written in the following form:

max j∈B∪N B n i=1 V ji f ji s.t. j∈B∪N B V ji + λ 1 i = 0, ∀ i = 1, . . . , n, , - j∈B∪N B V ji + λ 2 i = 0, ∀ i = 1, . . . , n, , -V ji + µ ji = 0, ∀ j ∈ B ∪ N B, i = 1, . . . , n, ||µ j || ρ ≤ γ j , ∀ j ∈ B ∪ N B, γ j = w j y j , if j ∈ B, γ j = w j , if j ∈ N B, λ 1 , λ 2 ∈ R n + , µ j ∈ R n , ∀j ∈ B ∪ N B, γ ∈ R n .
Clearly, the inner primal and dual problems satisfy Slater condition so that strong duality applies. This allows us to insert the optimality conditions in BLPτ to obtain the final single level program ConicP.

Computational Results

In the following we report some numerical results conducted to compare the efficiency of the different methods proposed to solve BLP, and to check experimentally their scope.

The computational experiments were carried out on a personal computer with Intel B. Core (TM) i7-4720HQ, 2.60 gigahertz with 16384 megabytes RAM. The MILP formulations and algorithms were implemented and solved by using Xpress Version 8.0.

The distances considered for the numerical experiments were computed using the 1 and ∞ norms. Therefore, we implemented the MILP formulations and algorithms proposed in Sections 3.1.3 and 3.2.3, in which we adapted the general methods in Section 3 to the models with these two particular distances.

For the computational study we generated different random instances taking into account the following factors: the dimension of the space, n, the cardinality of B and N B, which are the set of potential locations for the new primary facilities and the set of existing primary facilities, respectively, and also the maximum budget C. We considered the following levels for each factor:

• n = 2, 3, 10, 20, • |B| = 1000, 2000, 5000, 10000,

• |N B| = 1 4 |B|, 1 3 |B|, 1 2 |B|, • C = 1 C |B| = 1 4 |B|, 1 3 |B|.
The weights w j were generated randomly in the interval [0, 1], and each coordinate, f ji , of the location of the primary facilities f j was generated randomly in the interval [-1000, 1000], for all the instances.

For each combination of levels, 5 different instances were generated and solved. The CPU time limit to solve the problems was set to 1800 seconds.

In Figures 1 and2 we show the performance profile graphs of the number of solved instances for the different proposed models for the 1 -norm (Figure 1) and ∞ -norm (Figure 2). We represent in the abscissa axis the time (in seconds) and in the ordinate axis the number of solved instances. Figure 1 reports the results for the 1 -norm and it compares the two MILP formulations BLP-1 -1 and BLP-1 -2, the basic Benders algorithm, that we denote by Bend-1 , and the Benders algorithm using the separability property, denoted as Bend-1 -sep. Figure 2 shows the results for the ∞ -norm and it compares the MILP formulation, BLP-∞ -1, and Benders algorithm, denoted as Bend-∞ .

We can observe in Figures 1 and2 that the Benders algorithms are more efficient than the MILP formulations, in both cases with the 1 -and-∞ -norm cases. The Benders algorithms solve all the instances in very short time, whereas none of the MILP formulations could solve to optimality all the instances. We can see in the figures that the formulations BLP-1 -1 and BLP-1 -2 solve around 300 out of the 480 instances in 1800 seconds, and formulation BLP-∞ -1 solves around 400 instances in the same time.

The average number of cuts added in the Benders algorithm is 5, 03 for Bend-1 , 4, 75 for Bend-1 -sep and 4, 28 for Bend-∞ . The maximum number of Benders cuts, 14, was added for the Bend-1 -sep for an instance with n = 10, |B| = 5000, |N B| = 6667 and C = 1 4 |B|. For the 1 -norm case, in Figure 1, we can see that the Benders algorithm Bend-1 solves all the instances in approximately 200 seconds, whereas the Benders algorithm using the separability property, Bend-1 -sep, needs a bit more time. Nevertheless the performance of both methods is very similar. The same trend can be observed for the MILP formulations, the one without the separability property could solve in the end more instances within the same time limit. However, BLP-1 -2 works better for the big instances, n = 10, 20 (see Figure 3 in the Appendix).

With respect to the ∞ -norm case, Figure 2 shows that the Benders algorithm solves all the instances in less than 51 seconds, meanwhile BLP-∞ -1 only solves, in the same time, approximately one half of the instances.

More details about the Computational Results can be found in Tables 1 and2 in the Appendix. We report the average CPU times (CPU), and the numbers of problems, out of 5, solved to optimality (#OPT), for each type of instance and each formulation or algorithm.

Conclusions

This paper considers models for the location of controversial facilities. Controversial facilities must be understood as those facilities such that their placement induces a disagreement among users with different, non-aligned or opposite interests. Semiobnoxious facility location and the location and protection of critical infrastructures We model these situations by a bilevel optimization problem. The first level locates primary facilities trying to be as far away as possible from the secondary ones, which in turns, wish to be as close as possible to the primary ones. We develop mathematical programming formulations for the above mentioned bilevel programs as well as some algorithms that perform very-well in all our experiments that range from small problems on the plane (n=2) with up to |B| = 10000, possibilities for the primary facilities until dimension n = 20 and |B| = 10000.

BLP-1 -1 BLP-1 -2 Alg-1 Alg-1 -sep BLP-∞ -1 Alg-∞ n |B| |NB| C' #

  j∈B c j y j ≤ C, y j ∈ {0, 1}, j ∈ B, x ∈ arg min

Figure 1 .

 1 Figure 1. Performance profile graph of #solved instances for the different proposed models for the 1 norm.

Figure 2 .

 2 Figure 2. Performance profile graph of #solved instances for the different proposed models for the ∞ norm.

Figure 3 .

 3 Figure 3. Performance profile graph of #solved instances for the different proposed models for big instances (n = 10, 20) for the 1 norm.

INSTANCES

Table 1 .

 1 Numerical results for BLP under the l 1 and ell ∞ norm.

			OPT CPU	#OPT CPU	#OPT CPU #OPT CPU #OPT CPU #OPT CPU
	2 1000	250 3	5	10,82	5	10,36	5	0,41	5	0,50	5	2,13	5	0,56
	2 1000	250 4	5	4,80	5	5,89	5	0,62	5	0,52	5	2,22	5	1,03
	2 1000	333 3	5	13,39	5	10,07	5	0,35	5	0,41	5	2,19	5	0,65
	2 1000	333 4	5	7,57	5	4,98	5	0,55	5	0,62	5	2,47	5	1,05
	2 1000	500 3	5	11,70	5	12,76	5	0,38	5	0,43	5	8,53	5	0,81
	2 1000	500 4	5	7,86	5	6,35	5	0,54	5	0,52	5	3,02	5	0,95
	2 2000	500 3	5	36,59	5	37,27	5	0,75	5	1,44	5	6,34	5	0,79
	2 2000	500 4	5	13,36	5	11,66	5	0,65	5	0,73	5	9,15	5	1,36
	2 2000	667 3	5 217,01	5 152,79	5	0,68	5	0,85	5	7,20	5	0,92
	2 2000	667 4	5	20,58	5	14,54	5	0,45	5	0,60	5	6,14	5	1,17
	2 2000 1000 3	5	56,48	5	78,63	5	0,45	5	0,61	5 15,29	5	0,81
	2 2000 1000 4	5	29,85	5	20,60	5	0,52	5	0,65	5 12,74	5	1,56
	2 5000 1250 3	5 204,62	5 508,34	5	1,19	5	1,58	5 48,51	5	1,62
	2 5000 1250 4	5 124,68	5	97,10	5	1,49	5	2,36	5 52,24	5	1,70
	2 5000 1667 3	5 381,53	5 268,18	5	1,03	5	1,32	5 45,06	5	1,47
	2 5000 1667 4	5 135,94	5 198,89	5	2,08	5	1,66	5 37,01	5	1,96
	2 5000 2500 3	5 416,01	5 369,43	5	0,80	5	1,12	5 139,15	5	2,34
	2 5000 2500 4	5 107,17	5 405,78	5	1,23	5	1,89	5 38,32	5	1,92
	2 10000 2500 3	4 708,80	3 955,21	5	1,68	5	2,57	5 137,97	5	3,29
	2 10000 2500 4	5 325,03	5 355,51	5	2,44	5	3,75	5 81,66	5	3,72
	2 10000 3333 3	5 390,82	5 296,18	5	1,94	5	2,78	5 121,14	5	3,42
	2 10000 3333 4	5 329,09	5 439,71	5	3,09	5	3,46	5 159,98	5	3,16
	2 10000 5000 3	5 506,03	4 671,22	5	2,16	5	2,09	4 530,49	5	3,38
	2 10000 5000 4	5 516,73	5 416,86	5	2,83	5	3,27	5 141,33	5	3,77
	3 1000	250 3	5	9,67	5	19,03	5	0,21	5	0,27	5	5,95	5	0,70
	3 1000	250 4	5	58,74	5	43,20	5	0,63	5	0,71	5	3,09	5	1,18
	3 1000	333 3	5	21,51	5	50,24	5	0,36	5	0,37	5	6,99	5	0,60
	3 1000	333 4	5	29,05	5	32,52	5	0,65	5	0,90	5	4,68	5	1,21
	3 1000	500 3	5	26,14	5	28,70	5	0,34	5	0,49	5 10,17	5	0,57
	3 1000	500 4	5	33,46	5	47,28	5	0,38	5	0,52	5	6,86	5	0,73
	3 2000	500 3	5 120,41	5 144,41	5	0,38	5	0,53	5 63,77	5	1,19
	3 2000	500 4	5 389,05	5 301,49	5	0,76	5	1,22	5 21,33	5	1,41
	3 2000	667 3	5 316,33	5 232,75	5	1,04	5	0,92	5 31,18	5	1,04
	3 2000	667 4	5 129,52	5 224,55	5	0,47	5	0,56	5 21,81	5	1,48
	3 2000 1000 3	5 113,90	5 261,64	5	0,65	5	0,53	5 81,42	5	1,38
	3 2000 1000 4	5 179,31	5 249,42	5	0,62	5	0,71	5 11,05	5	0,88
	3 5000 1250 3	4 1014,05	3 1310,24	5	1,95	5	1,85	5 131,33	5	1,94
	3 5000 1250 4	3 1352,73	2 1442,08	5	2,43	5	2,24	5 106,38	5	2,12
	3 5000 1667 3	4 763,66	4 669,57	5	1,14	5	1,28	5 180,78	5	1,94
	3 5000 1667 4	3 1102,02	4 1047,36	5	1,56	5	1,91	5 213,59	5	2,08
	3 5000 2500 3	5 636,19	5 689,33	5	0,76	5	1,14	5 391,97	5	2,07
	3 5000 2500 4	5 592,53	2 1223,05	5	1,42	5	1,58	5 77,89	5	1,90
	3 10000 2500 3	1 1778,05	1 1800,00	5	2,42	5	3,38	5 899,86	5	3,09
	3 10000 2500 4	1 1727,79	0	-	5	5,59	5	5,16	5 198,48	5	4,18
	3 10000 3333 3	3 1127,45	1 1732,91	5	2,67	5	3,91	5 832,12	5	3,17
	3 10000 3333 4	3 1253,69	2 1321,56	5	4,79	5	6,15	4 797,61	5	3,91
	3 10000 5000 3	3 1379,12	1 1771,16	5	2,86	5	4,41	4 470,05	5	4,01
	3 10000 5000 4	0	-	0	-	5	4,18	5	5,79	5 325,39	5	4,87

Table 2 .

 2 Numerical results for BLP under the 1 and l ∞ norm.

Acknowledgements

Martine Labbé has been partially supported by the Fonds de la Recherche Scientique -FNRS under Grant(s) no PDR T0098.18. The last two authors have been partially supported by MINECO Spanish/FEDER grants number MTM2016-74983-C02-01.

Appendix


Sébastien Seguy
email: asebastien.seguy@insa-toulouse.fr

Gilles Dessein
email: bgilles.dessein@enit.fr

Lionel Arnaud
email: clionel.arnaud@enit.fr

Control of chatter by spindle speed variation in high-speed milling

Keywords: chatter, milling, stability, spindle speed variation, surface roughness

High-speed milling operations are often limited by regenerative vibrations. The aim of this paper is to analyze the effect of spindle speed variation on machine tool chatter in high-speed milling. The stability analysis of triangular and sinusoidal shape variations is made numerically with the semi-discretization method. Parametric studies show also the influence of the frequency and amplitude variation parameters. This modeling is validated experimentally by variable spindle speed cutting tests with a triangular shape. Stable and unstable tests are analyzed in term of amplitude vibration and surface roughness degradation. This work reveals that stability must be considered at period variation scale. It is also shown that spindle speed variation can be efficiently used to suppress chatter in the flip lobe area.

Introduction

The productivity of cutting processes is severely limited by vibrations that arise during the material removal process. One reason for such vibrations is the surface regeneration. The theory of regenerative machine tool chatter is based on the work of Tobias and Fishwick [START_REF] Tobias | Theory of regenerative machine tool chatter[END_REF]. This knowledge initially dedicated to the turning process has been adapted to milling [START_REF] Altintas | Analytical prediction of stability lobes in milling[END_REF] and led to the development of the stability lobes theory. Since then several improved model and analysis techniques have been appeared including detailed analysis of the governing delay differential equation and time domain simulations, see, for example [START_REF] Merdol | Multi frequency solution of chatter stability for low immersion milling[END_REF][START_REF] Insperger | On the higher-order semi-discretizations for periodic delayed systems[END_REF]. These models all use the so-called stability lobe diagrams, which allows to choose the spindle speed associated with the axial depth of cut with a chatter free machining. In many practical cases, the choice of the optimal speed is difficult because contradictory parameters interact with productivity [START_REF] Seguy | Toolpath dependent stability lobes for the milling of thin-walled parts[END_REF].

Another technique to reduce chatter vibrations is the spindle speed variation. As opposed to variable pitch cutters, spindle speed variation can effectively be used in a wider spindle speed range, since the frequency and the amplitude of the speed variation can easily be adjusted in CNC machines even during the machining process. In the 1970's, Takemura et al. [START_REF] Takemura | Active suppression of chatter by programmed variation of spindle speed[END_REF] presented the first simple model to study the stability of variable speed machining; they predicted significant shift of the stability lobes to higher depth of cuts, but the experimental tests showed only small improvements. Sexton and Stone [START_REF] Sexton | A stability analysis of single point machining with varying spindle speed[END_REF] developed a more realistic model, and they found some improvements in the stability properties for low spindle speeds. Moreover, they showed that the presence of transient vibrations may further reduce these gains.

The study of the stability for variable speed machining requires a special mathematical analysis to compute stability lobes. Tsao et al. [START_REF] Tsao | A new approach to stability analysis of variable speed machining systems[END_REF] have developed a model taking the angular coordinates as variables instead of time. Insperger and Stépán [START_REF] Insperger | Stability analysis of turning with periodic spindle speed modulation via semidiscretization[END_REF] showed that the semi-discretization method can effectively be used for the stability analysis of turning at variable speed. They showed that the critical depths of cut can be increased for low speeds, but for the high-speed domain, no improvement was found. Recently, Zhang et al. [START_REF] Zhang | Stability analysis on spindle speed variation method for machining chatter suppression[END_REF] presents a systematic stability analysis of spindle speed variation based on a machining chatter model of non-linear delay differential equation. Experimental results were presented by Al-Regib et al. [START_REF] Al-Regib | Programming spindle speed variation for machine tool chatter suppression[END_REF] showing gains on surface roughness for turning at low speed.

The modeling of variable speed milling is more complex than that of turning, since the speed variation frequency and the tooth passing frequency interact and the resulting system is typically quasi-periodic. Still, there are mathematical techniques to determine approximate dynamic properties. Sastry et al. [START_REF] Sastry | Floquet theory based approach for stability analysis of the variable speed face-milling process[END_REF] used Fourier expansion and applied the Floquet theory to derive stability lobe diagrams for face milling. They obtained some improvements for low spindle speeds. Recently, Zatarain et al [START_REF] Zatarain | Stability of milling processes with continuous spindle speed variation: analysis in the frequency and time domains, and experimental correlation[END_REF] presented a general method in frequency domain to the problem, and show that varying spindle speed can effectively be used to chatter suppression. Another approach is to use time domain simulation [START_REF] Altintas | In-process detection and suppression of chatter in milling[END_REF] that makes it possible to obtain more detailed information like the amplitude of vibrations, the chip thickness or the cutting forces. Simulations were sometimes made at high speed, but to the best knowledge of the authors, experimental validations are only made at low spindle speed [START_REF] Al-Regib | Programming spindle speed variation for machine tool chatter suppression[END_REF][START_REF] Sastry | Floquet theory based approach for stability analysis of the variable speed face-milling process[END_REF][START_REF] Zatarain | Stability of milling processes with continuous spindle speed variation: analysis in the frequency and time domains, and experimental correlation[END_REF][START_REF] Altintas | In-process detection and suppression of chatter in milling[END_REF][START_REF] Yang | Chatter suppression with multiple time-varying parameters in turning[END_REF].

In this paper, the stability of variable speed milling is analyzed in the high-speed domain, for spindle speeds corresponding to the first flip (period doubling). Theoretical stability predictions are obtained using the semi-discretization method based on [START_REF] Insperger | On the higher-order semi-discretizations for periodic delayed systems[END_REF], and the results are confirmed by experiments.

Milling modeling

Variation of the spindle speed. In the literature, mostly sinusoidal, triangular or square-ware modulations are considered [START_REF] Lin | The effects of variable speed cutting on vibration control in face milling[END_REF]. Here, the sinusoidal and the triangular variation shown in Fig. 1 are compared and analyzed. Assume that the spindle speed variation is periodic at period Tv with a mean value N0 and an amplitude Na, that is,

      0 va N t N t T N N S t    
, where     v S t S t T  is the shape function. For a triangular modulation the shape function is defined as:

          1 4mod , / if 0 mod , / 2 3 4mod , / if / 2 mod , t T T t T T St t T T T t T T              . ( 1
)
Here, mod(t,T) denotes the modulo function, for example, mod(12, 5) = 2. For a sinusoidal modulation  

St is defined as: (3)

  2 sin v t St T      . ( 2
) N 0 T v /2 T v 3T v /2 t [s] N 0 + N a N 0 -N a 2T v 0 N(t) [tr/min] N 0 T v /2 T v 3T v /2 t [s] N 0 + N a N 0 -N a 2T v 0 N(t) [tr/min]
RVA represents the amplitude ratio between the amplitude Na and the mean value N0, it is always less than 1. In order to have reasonable cutting conditions, its value was limited at 0.3. This represents a variation of 30% of the spindle speed. This represents also a variation of 30% of the feed by tooth due to the constant feed velocity. RVF is the ratio between the variation frequency fv and the average spindle frequency N0. It may be higher than 1. If RVF equal to 1, then the spindle makes one variation period during only one revolution of the tool.

Mechanical model.

A schematic diagram of the milling process is shown in Fig. 2. The structure is assumed to be flexible in the x direction, while the feed is parallel to the y direction.

Fig. 2: Mechanical model of the milling process

The dynamic model is given by the following equation:

        x mx t cx t kx t F t    , ( 4
)
where m is the modal mass, c is the damping, k is the stiffness and  

x Ft is the cutting force in the x direction. According to a linear cutting law, the x component of the force is given by:

      1 cos sin z x p R j T j j j F t A K K h t       , ( 5
)
where Ap is the axial depth of cut, j  is the angular position of the jth cutting edge and KT and KR are the specific tangential and radial cutting coefficients. The chip thickness is expressed by:

            sin cos j j z j j h t g t f x t x t t          , ( 6
)
where the function   j gt is a unit function, it is equal to 1 when the tooth j is cutting, otherwise it is equal to 0. Here, fz is the feed per tooth,   xt is the current position of the tool and

    x t t  
is the position at the previous cut. The regenerative delay   t  is periodic in time due to the spindle speed variation.

Stability analysis.

Stability of the milling process with variable spindle speed can be estimated theoretically via the analysis of the governing Eqs. ( 4), ( 5) and [START_REF] Takemura | Active suppression of chatter by programmed variation of spindle speed[END_REF]. According to the Floquet theory of delayed differential equations [START_REF] Insperger | On the higher-order semi-discretizations for periodic delayed systems[END_REF], an infinite dimensional operator  (monodromy operator) can be associated with the system that gives the connection between the current state of the cutter and the state one period earlier:

    x t T x t    .
(7) This period T is equal to the principal period of the system. Thus, for the milling model with spindle speed variation, the regenerative time delay  and the spindle modulation period Tv are periodic. We assume that the ratio of the modulation period Tv and the mean time delay  is a rational number, i.e., v qT p  with p and q being relative primes. Thus, the system is periodic at the principal period v qT T  consequently, the Floquet theory of periodic DDEs can be applied. Note that if the ratio of Tv and  is not rational, then the system is quasi-periodic and the Floquet theory cannot be used. Stability properties are described by the eigenvalues of the monodromy operator. If all the characteristic multipliers are in modulus less than 1, then the system is asymptotically stable. Here, the semi-discretization method [START_REF] Insperger | On the higher-order semi-discretizations for periodic delayed systems[END_REF] is used. The basic point of this method is that they provide a finite dimensional matrix approximation of the monodromy operator and the eigenvalues of this matrix can be determined numerically. This method was validated by time domain simulation [START_REF] Seguy | From the spindle speed selection to the spindle speed variation for chatter control in thin wall milling: modelling and experiments[END_REF].

Selection of the optimal parameters

Optimal area. Stability lobe can be constructed by scanning the cutting conditions (spindle speed and axial depth of cut) for a couple of (RVA, RVF) parameters and for a given shape. Fig. 3 presents the lobes for constant, triangular and sinusoidal spindle speed.

The critical depth of cut can be seen to be increased by speed variation for some ranges of spindle speeds, but for some other ranges, the critical depth of cut is less than that of the constant spindle speed. For example, a cutting process with an axial depth of cut of 1 mm and a spindle speed of 9100 rpmthat is unstable for constant spindle speedcan be stabilized by a speed variation. But globally, the critical depth of cut can be seen to be increased essentially in the area of the first flip lobe by spindle speed variation. Comparison of the shape variation. The effectiveness of the spindle speed variation is tested on a flip lobe area (N0 = 9100 rpm) with a triangular and sinusoidal shape variation. Fig. 4 shows a contour plot presenting the maximal depth of cut available without chatter for the amplitude (vertical axis) and frequency (horizontal axis) with an average spindle speed of 9100 rpm. At constant speed, the maximal depth of cut is 0.5 mm. It can be seen that the triangular shape always increase the critical depth of cut Ap for any RVA and RVF values. For some domains, even Ap = 2.4 mm can be achieved that corresponds to 380% improvements. As mentioned earlier, the choice of the frequency and the amplitude variation is limited by the spindle dynamics. Within the range of variation, the maximum acceleration of the spindle is 100 rev/s 2 , i.e., 6000 rpm/s. Considering the limits of the dynamic spindle for the triangular shape, the optimal choice is to use low frequency modulation with high amplitude. Such a point is denoted by point D or E. Fig. 4: Parametric study for triangular and sinusoidal shape (N0 = 9100 rpm and Ae = 2 mm) A similar plot is determined for a sinusoidal shape in Fig. 4. Like the triangular shape, the sinusoidal shape always increases the critical depth of cut. For some area, an improvements of 500% can be achieved (Ap = 3 mm). In fact, for a same frequency and amplitude parameters, the sinusoidal shape is more effective than the triangular shape. But if we consider the limit of the spindle, the sinusoidal shape givesfor the machinistonly an improvement of 300%. In fact, the sinusoidal shape needs more acceleration than the triangular shape (see Fig. 1). For a given spindle dynamic, the sinusoidal shape allows no more improvements. Surprisingly, a high frequency of variation, coupled with a low amplitude variation, produce a little gain. The most effective parameter is the amplitude variation.

Because triangular shape allowsfor a given dynamic spindlethe maximal gain for the machinist, only this spindle speed variation was analyzed in the experimental part.

Experimental work

The machining tests were carried out on a high-speed milling center. The average feed by tooth is 0.1 mm/tooth. The tool is an inserted mill with 3 teeth, 25 mm diameter without helix angle. The spindle speed rotation is managed by a sub-program using a synchronous function.

A removable part in aluminum alloy (2017A) is machined in contouring with a radial depth of cut (Ae) of 2 mm. Fig. 2 shows the experimental set-up of the single degree of freedom flexure, compliant in the x direction. The tool is considered infinitely rigid compared to the workpiece. The vibrations of the part were measured by a laser velocimeter.

Collected in Table 1, the dynamic characteristics of the system were determined by hammer impact test and the cutting force coefficients were determined in coherence with previous work [START_REF] Seguy | Toolpath dependent stability lobes for the milling of thin-walled parts[END_REF]. Some special cases. Firstly, cutting tests were conducted at spindle speed of 9100 rpm with RVA = 0,005, RVF = 0,15 ( fv = 22,75 Hz) and a depth of cut of 1mm. This test uses the maximum acceleration of the spindle with high frequency variation and low amplitude. The results thus obtained are presented in Fig. 5. Test (A) refers to the variable speed machining, test (B) to the constant spindle speed machining (note that these cutting parameters are in the first flip lobe, see point B in Fig. 3). Based on the theoretical predictions, these two tests are unstable.

The displacement measured with and without spindle speed variation is the same. Moreover, these two unstable machining generates exactly the same poor surface roughness on the workpiece. For this unstable machining, chatter is clearly identified.

Both machining shows also the minor importance of the variation frequency, like simulation predicts. It is not necessary to change very fast the spindle speed, even with maximum acceleration, because this is not the most important parameter.

Then, other cutting tests were performed at spindle speed of 9100 rpm and depth of cut of 1 mm using variable spindle speed with RVA = 0.08, RVF = 0.0125 ( fv = 1.9 Hz). This variation corresponds to point C in Fig. 4. Based on the theoretical predictions, the critical depth of cut is about 0.8 mm, i.e., the process is unstable.

Fig. 6 shows the recorded displacement history and the surface roughness variation during the cutting process. After about two periods of speed variation, the amplitude of the displacement increases very strongly. The envelope curve of the vibrations is illustrated by dashed lines. This behaviour corresponds to the predicted instability. However, chatter does not appear suddenlysome time is needed for it fully develops. In the case of our cutting test, this time span is about 0.7 s. If the duration of the machining process is less than this period, then large amplitude vibrations do not develop. In a sense, this case can be considered as "practically stable" machining, in spite of the fact that the process itself is globally unstable. If the duration of the machining process is longer, then chatter develops. This effect can also be observed on the surface roughness of the workpiece. The left hand side of the machined profile, where the cutting was started, has a pitch of approximately 0.3 mm, with roughness of 1.6 μm. For an ideally symmetric tool, the pitch of the machined profile is equal to the feed per tooth. However, if the tool has a runout greater than the roughness of the surface, then it leaves only one mark per revolution. The tool used in the tests had a runout of 10 μm, the feed per tooth was 0.1 mm and the tool had 3 teeth, thus the pitch of the machined profile is expected to be approximately 0.3 mm for stable machining. In fact, the pitch varies slightly around 0.3 mm, since the constant feed velocity and the variable spindle speed produce a varying feed per tooth. As the tool passes over and chatter develops, the surface roughness gets worse and worse. By the end of the workpiece, the roughness increased to 13.9 μm.

t [s] t [s] t [s] t [s] Experimental displacement [m]
Stabilization by spindle speed variation. Consider the unstable machining process with spindle speed of 9100 rpm and depth of cut of 1 mm. Spindle speed variation is applied according to point D in Fig. 4. The corresponding parameters are RVA = 0.2, RVF = 0.0046875 ( fv = 0.71 Hz). Based on the theoretical predictions, the critical depth of cut is about 2 mm, i.e., the system with variable spindle speed is predicted to be stable. In Fig. 7, test (D) refers to the variable speed machining and test (B), to the constant spindle speed machining.

During the test with variable spindle speed, no chatter was observed. The amplitude of the vibrations was less than 0.01 mm, the roughness was 1.75 μm and the pitch of the machined profile was 0.3. For a 3 fluted tool with runout of 10 μm and feed per tooth 0.1 mm, this corresponds to stable machining. These all refer to a stable cutting process. For constant spindle speed, chatter was clearly identified. The amplitude of the vibrations was about 0.07 mm, the roughness was 3.7 μm and the pitch of the machined profile was 0.6 mm that refers to the period doubling chatter (see point B in Fig. 3). In these cases, identification of chatter was unambiguous, the pitches of the machined profiles were uniform all along the workpiece.

Conclusions

Variable spindle speed machining was studied for high-speed milling at around the first flip lobes. Stability properties were predicted using the semi-discretization method. Moreover different combinations of the amplitude and the frequency of the speed modulation were analyzed in order to find the optimal technique to suppress chatter. It was found that the stability properties can always be improved (i.e., the critical depth of cut can always be increased) by spindle speed variation within the unstable domain of the first flip lobe. Amplitude was also shown to have a greater effect on the stability of the process than frequency. In this way, for a given dynamic spindle, the triangular shape allows the maximal gain for the machinist.

Cutting tests were performed for certain spindle speeds in the flip domain in order to verify the theoretical predictions. The concept of stability for variable spindle speed machining was shown to differ slightly from than of constant spindle speed, since the principal period of the system is equal to the spindle variation period (or to its integer multiple) instead of the tooth passing period. If the machining process is unstable, but the development of chatter requires more time than the duration of the process, then it can be considered as "practically stable" machining. The stabilizing effect of spindle speed variation was clearly verified experimentally, a period doubling chatter was suppressed by applying a proper spindle speed variation.

Fig. 1 :

 1 Fig. 1: Typical triangular and sinusoidal shape variation

Fig. 3 :

 3 Fig. 3: Stability diagrams for variable speed milling with RVA = 0.3 and RVF = 0.003

Fig. 5 :Fig. 6 :

 56 Fig. 5: Two machining unstable (Ap = 1 mm and N0 = 9100 rpm)

Fig. 7 :

 7 Fig. 7: Chatter suppression by spindle speed variation (Ap = 1 mm and N0 = 9100 rpm)

Table 1 :

 1 Flexure modal parameters and cutting coefficients

	m [kg]	f0 [Hz]	ξ [%]	KT [MPa]	KR [MPa]
	1.637	222.5	0.50	700	140

Acknowledgments

This work was supported in part by the European Union Interreg IIIa AEROSFIN, by the Midi-Pyrénées Region Project "Complex workpiece machining in aeronautical context", by the János Bolyai Research Scholarship of the HAS and by the HSNF under grant no. OTKA K72911. The authors would also like to thank Grégoire PEIGNÉ for his valuable suggestions.