Marcella Autran

Burlier Drummond

Luiz Gustavo
email: lgustavozelaya@gmail.com

Amaranto Zelaya Cruz

Lopes Pereira

Felix Mora-Camino

Luiz Gustavo

Felix Mora
email: felix.mora@enac.fr

M.Sc Marcella A B Drummond
email: marcella.drummond@enac.fr

Luis Gustavo

Dr.-Ing Amaranto L Pereira
email: amaranto@pet.coppe.ufrj.br

Dr Félix Mora-Camino

Planning airports within regional multimodal transportation systems

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

RESUMO

O objetivo deste trabalho é propor uma ferramenta de avaliação de impactos advindos pela instalação de um novo aeroporto em uma região onde já exista uma rede de transportes multimodal. Este problema é crucial na escolha da localização de aeroportos e no planejamento de terminais, de conexões de vôo e de frotas utilizadas pelo aeroporto. Uma das principais dificuldades desta tarefa é obter uma boa previsão da demanda futura, considerando a inserção do novo aeroporto na rede de transporte multimodal existente. Tais impactos influenciam diretamente na avaliação das condições de operação, custos e receitas do novo aeroporto. A abordagem proposta utiliza-se de dois modelos de otimização diferentes: Um modelo é dedicado à previsão da demanda de passageiros e carga, considerando a alteração da acessibilidade da rede de transporte multimodal e o outro define a oferta global de transporte de acordo com o comportamento de maximização do lucro no sistema de transporte implicado. O processo de previsão da demanda é baseado em uma abordagem de maximização de entropia com restrições de níveis gerados ou atraídos pelas atividades econômicas dos distintos centros urbanos. Isso permite determinar a intensidade e a distribuição das matrizes de origem-destino por tipo de produto e modalidade de transporte. O esquema de solução proposto é composto de um processo interativo entre a solução da previsão da demanda e o problema de otimização da oferta: o problema de distribuição da maximização da entropia gera a matriz de origem-destino, dada uma estrutura de custo/capacidade, enquanto o problema de otimização da oferta gera essa estrutura de custo/capacidade caracterizando a acessibilidade, dadas as matrizes de origem-destino. A abordagem proposta é ilustrada com um estudo de caso de inserção de um novo aeroporto de médio porte em uma área rural de rápido desenvolvimento.

PALAVRAS-CHAVE: Rede, otimização, transporte multimodal, aeroportos.

ABSTRACT

In this communication is proposed an appraisal tool for measuring the impacts of the installation of a new airport within a region already connected by a multimodal transportation network. This problem is crucial when choosing location for airports and then when planning the investments in terminals, connection links and fleets. One of the main difficulties of this task is related with the estimation of future demand over the modified transportation network which has direct influence on the evaluation of operational conditions, revenues and costs of the new airport. The proposed approach makes use of two different optimization models: One model is devoted to passenger and freight demand forecasting taking into account the modified accessibility of the multimodal transportation network, the other one defines the global transport supply according with a profit maximization behavior for the involved transport system. The demand forecasting process is based on an entropy maximization approach with distinct urban centers' economic activities generated or attracted constraints levels to determine the intensity and the distribution of origin-destination matrices by the kind of product and transport modality. The proposed solution scheme is composed of an iterative process between the current solution for demand forecasting and the supply optimization problem: the entropy maximizing distribution problem provides the origin-destination matrix given a cost/capacity structure, while the supply optimization problem provides this cost/capacity structure characterizing accessibility, given the origin-destination matrices. The proposed approach is illustrated with the case of the study of the insertion of a new airport of medium size in a rapid developing rural area.

KEY WORDS:

Network, optimization, multimodal transportation, airports.

INTRODUCTION

Transportation and economic development have always been entwined sectors [START_REF] Jain | A New Broad Composite Model of Regional Economic Development as Framework for Research[END_REF], [START_REF] Petrakos | A Generalized Model of Regional Economic Growth in European Union[END_REF]. It is admitted that, in general, for transport as a whole, the increase in the provision of transport services in a country increases the overall productivity, which improves the national GDP. However, the relationship between transportation and economic activity is quite complex [START_REF] Ishutkina | Analysis of the interaction between air transportation and Economic Activity: A worldwide Perspective[END_REF].

The main objective of this study is to provide a framework to assess the impact of the implantation or the development of a new terminal facility in a multimodal transportation network in a fast developing regional economy. This is done in two steps: first a model of transportation demand is developed for passengers and freight, then supply of transportation service is optimized leading to the definition of necessary terminal capacities.

In this study, the proposed approach is illustrated with the case of the creation of a new airport in a rural area where ground transportation networks already provide transportation services to passengers and freight. Considering air transportation, it provides a rather low investment cost solution to create access to people, markets, ideas and capital. Thus air transportation can be an efficient enabler for regional economic development.

PREDICTION OF GENERATED PASSENGERS TRANSPORTATION DEMAND

Here it is supposed that demand for air transportation at a given airport is the result of different factors: the distribution and composition of socio-economic activities, the competition with other transportation modes and the competition with other airports. This leads us to adopt respectively the following generation and attraction models for passengers demand.

A possible model for passenger generation is:

i ar i a D ar i ar i pax iar V Pop D      R r A a I i i    , , (2)
where pax iar D is the annual passenger attraction to urban center i with respect to people related with activity a, Here we consider that transportation modalities from a given urban center to and from another urban center can be distinguished, from the user point of view, by their generalized transportation costs. Here generalized transportation costs include fares, transportation delays, waiting times and comfort. Let m ij  be the generalized transportation cost from center i to center j through transportation modality m. Then, considering that an a priori distribution of trips [ arm ij T ˆ] is available for each pair (a,r) and each modality m, as well as an a priori generalized transportation cost from center j to center i with each modality m, m ji ˆ [START_REF] Drummond | An optimization approach for long planning of multimodal freight transportation systems[END_REF], (Drummond, Mora-Camino, & Pereira, 2008), we can consider the solution of the following max entropy problem [START_REF] Handou | Contribution à l'Optimisation d'un Réseau de Transport Aérien: Proposition d'un Modèle basé sur la Logique Floue et la Maximisation Entropique[END_REF] between each center and for all pair (a,r):

  ) / ( max arm ij m ij arm ij m ij arm ij m ij M m J j R r A a I i T T T Log T ij i arm ij               (3) R r A a I i D O T i pax iar pax iar arm ij M m J j ij          , , (4)


M m i I j I i A a T i arm ij       , , , 0 (5)
where arm kj T is the distribution trip from urban center k to center j, for each pair (a,r) and using the same modality m.

Observe that this is a new version of the entropy maximization model where instead of traffic flows, transportation costs have been introduced since it is considered here that they are the main drivers of travel behavior without being exclusive. This procedure will be also adopted for the other predictive models.

PREDICTION OF FREIGHT DEMAND DISTRIBUTION

Here we distinguish for a given location between freight imports and exports. It is considered that for activity sector a a i a V level of production needs an amount of i a i aa V  to be imported, while the final re-exported production is given by

i b i ba a b A b i a V V i      , , where the positive coefficients i ba  are Leontief multipliers. Let i r i ar i A a R r Pop i      
be the imported final costumer consumption at location i and let im C Im and im Ex C be respectively the importation and exportation handling and storing capacities at location i for transport mode m. An estimate of the modal and spatial distribution of imported goods by center i is given by:

  ) / ( ) / ( max Im am ij m ij am ij m ij am ij m ij M m J j A a am ji m ji am ji m ji am ji m ji M m J j A a I i P P Log P P P Log P i Ex a i i a i am ji                        (6) I i A a Pop V i i r i ar i R a R r i a aa am ji M m J j i i a             , Im Im     (7) I i A a V V i i b i ba a b A b i a am ij M m J j i i Ex a            , ,   (8)
i I i am ji M m J j I i am ij M m J j I i A a Ex a a                   Im (9) I i M m A a J j i a am ji      , , , 0 Im  (10) I i M m A a J j i Ex a am ij      , , , 0  (11)
where am ji  is the imported goods flow of sector a with the modality m from j to the urban center i,

TRANSPORTATION SUPPLY MODEL

The proposed supply model takes into account m types of flows: m vehicle flows according to the different transport modalities and freight flows using the resulting transport capacity. The fleets of vehicles and their operation generate fixed and variable costs, while freight flows are the main source of revenue for the freight transportation sector. To each transportation modality m is associated a network linking the N transportation centers such as:

R m = [G m ,[ m ij f ]] (12)
where G m is the graph (Berge, 2001) associated to the m th modality, [ m ij f ] is the flow of vehicles of the m th modality over G m . The flow of vehicles associated to each transportation mode satisfies to conservation and positivity constraints:

M m f f m ik G i j N G k i m ji m m m     ) , () , ( (13)
0  m ij f (i, j) G m , mM (14)
where N m is the number of vertices of graph G m . Here flows integrity constraints are not taken into account since no scheduling or routing problem will be formulated according to these flows which should only provide a global view of the future development of the multimodal transportation system. However for sake of realism, a fleet capacity constraint (Park & Regan, 2006) can be introduced: , the optimization problem of the supply of capacity to transport passengers and freight can be considered. To provide a tractable formulation of this problem, a global performance index for the whole transportation sector over a future period of time is adopted. In the case in which passengers and freight transportation [START_REF] Semet | Recherche Operationnelle et Transport de Marchandises[END_REF], [START_REF] Crainic | Planning models for freight transportation[END_REF], [START_REF] Lubis | Multimodal freight transport network planning[END_REF] are segregated, this will result in two independent problems. But when the fleet of some modality (air transportation, interurban bus transport, etc) can be used simultaneously to take care of the two types of transportation activities, these two problems are linked at least by common capacity constraints. Then the optimization criterion of the global problem can be given by:

m m m ij N G j i m ij F D d f m m    ) , ( M m  (
M m G j i f H m m ij m ij P m am ij aF m A a        , ) , ( 0    (19)
and under demand level constraints:

M m G j i A a am ij am ij      , ) , ( , 0   (21) M m G j i T H i A a arm ij R r m ij       , ) , ( 0 (22)
with positivity conditions:

                                 )) ( ( ) ( ) ( ) ( max 0 ) , ( 0 ) , ( , , m ji m ij i A j m T m T A i m ij Hm ij am ij am ij A a G j i m ij Fm ij m m vF m F M m G j i m ij m ij am ij am ij A a M m H f f f c c H c c f c F c c H P m m m m ij m ij m ij      CONINFRA
M m G j i A a H m ij am ij      , ) , ( , 0 0  (23)
where aF m  is a freight occupancy rate and P m  a passenger occupancy rate with respect to modality m.

When the solution of the above problem is such that for modality m:

 A a i I j I i      , ,
such as

am ij am ij    or am ij am ji .    (24)
this means that all demand for freight transport is not satisfied. Since at this stage, no terminal capacity has been imposed, this can only be the consequence of insufficient fleet size for modality m:

m G j i M m     ) , ( , such that m ij m ij P m am ij aF m A a f H        (25)
and fleet size of modality m should be increased by

m F  according to: m m m ij m ij P m am ij aF m A a N G j i arm ij A a R r m ij P m am ij aF m A a F D d H T m m              )) ( ) (( ) , (       (26)
Then the supply optimization problem can be restarted with this new fleet constraint level and new fixed costs m F c 0 . Then once a solution satisfying constraints:

M m G j i A a am ij am ij     , ) , ( ,   (27) M m G j i T H i A a arm ij R r m ij      , ) , ( (28)
is obtained, a lower bound for the necessary capacity of the different terminals can be estimated:

M m I i C f f C m i m ji i j j m ij G j i j m i m         , ) , ( ) , ( (29)
If these levels are under existing capacity at some terminal, increased capacity should be pursued for this terminal and the supply optimization problem should be restarted with an increasing terminal fixed costs m T c 0 . In the case in which the global performance of the transportation sector is under acceptable levels, transportation prices should be reviewed as well as the impacted generation and attraction levels for freight and passengers through price sensitivities such as: The implantation of a new airport in a center which had to use other transportation modalities including remote airports at other centers will have many impacts [START_REF] Caris | Planning problems in intermodal freight transport: Accomplishments and prospects[END_REF] on local development:

m ij O ar     / , m ij O ar     / , m ij i ar     / , m ij i ar P   /  , m ij i a P V   / , m ij i r P   /  (
-Transportation costs and volumes will be redeployed among the different modalities; -Sectorial activity levels will be redefined; -Import/export flows will be redistributed, -The revenue distribution among the population will be modified on the medium/long run. The above demand and supply models can be used to appraise the marginal effect created by this new situation. Then, let us consider first an initial situation in which modality m 0 is not present at location i and the volumes of activity and transportation, satisfying all the above capacity and interrelations, are given by:

  J j m M m I k A a D O T V i ij i am ji am ij pax iar pax iar arm kj i a      , , , , , , , , 0 
      
(31) with a distribution of revenues given by: R r A a Pop

i i ar i   ,  (32)
Then, introducing the new airport facility at center i, for the same level of activity in the different economic sectors, importation and exportation flows can be considered to be solutions of the following problems:

)) / ( ) / ( ( max Im am ij m ij am ij m ij am ij m ij M m J j A a am ji m ji am ji m ji am ji m ji M m J j A a I i P P Log P P P Log P i Ex a i i a i am ji                        (33) I i A a Pop V i i r i ar i R a R r i a aa am ji M m J j i i a             , 0 Im Im     (34) I i A a V V i i b i ba a b A b i a am ij M m J j i i Ex a            , , 0   (35) i I i am ji M m J j I i am ij M m J j I i A a Ex a a                   Im (36) I i M m A a J j i a am ji      , , , 0 Im  (37) I i M m A a J j i Ex a am ij      , , , 0  (38)
while travel demand will be now such as: Other economic scenarios, considering special development projects, can be considered in this same framework by adjusting the right hand side of equations 34, 35 and 40. Anyway, the net revenue of activity a will be affected by the redistribution of import and export freight flows and will be changed by:

) / ( max arm ij m ij arm ij m ij arm ij m ij M m J j R r A a I i T T T Log T ij i arm ij               (39) R r A a I i D O T i o pax iar o pax iar arm ij M m J j ij          , , (40)
     i am ji m ji J j am ij m ij J j am ji am ji m ji m M m J j am ij am ij m ij m M m J j i a A a P P P P R i a i i Ex a i ij a i ij Ex a                       ) ) ( ) ( 0 Im 0 0 0 Im 0         (42)
The variation of net revenue of transportation mode m ( m 0 ) is given by

pax im R  + fre im R  where:   I i m M m T T c R i ij J j arm kj arm kj pax ijm m ij I k J j R r A a pax im ij ij                 , ) ( ) ( 0    (43)     I i m M m c P R i ij J j am ij am ij am ji am ji fre ijm m ij J j A a fre im               0 ) ( ) ( ) (        (44)
Then, depending of the above estimated benefits/losses for each activity sector and transportation modality, a first sight to potential impacts of the implantation of the terminal facility can be obtained. When an activity sector is found to be beneficiated by the implantation, many opportunities open to it: from keeping away the profit to reinvest it in that same location or others [START_REF] Petrakos | A Generalized Model of Regional Economic Growth in European Union[END_REF]. With respect to the ground transportation modalities, even if at first sight the implantation of the new airport seems to retrieve some of their demand, the possible increase of activity sectors should be able to overcompensate the loss of this specific demand. The resulting specialization of ground transportation modalities should be also beneficial to them through increased cost rationalization and to users by improving service quality.

To pursue this analysis, scenarios built up from investment policies in activity sectors and transportation networks, terminals and fleets, should be established. In that case the proposed structure could be used in inter temporal analysis where the lower level (demand forecasting) and the upper level (activity levels) interact along the chosen horizon of time [START_REF] Handou | A two level optimization approach for the long term planning of a large air transportation network[END_REF].

ILLUSTRATION OF PROPOSED APPROACH

Here we considered the application of this framework to the implantation of a new regional airport in a fast developing rural agro-industrial area in the center of Brazil.

The area under study is crossed by a main national highway while a new freight railway is planned and has been considered in the study.

The planned airport is supposed to operate light and single aisle aircraft directly towards these regional and national centers on a weekly basis (turbojets) and on demand (light aircraft). A runway of 2.2 K m is projected with three taxiways. The airport capacity is estimated on 300 K/year for passenger and 150 K ton/year for freight. The passenger terminal building will be about 10 K m 2 and the cargo terminal building about 3 K m 2 , within the total airport area of 18 M m 2 . is the main operated aircraft, but the dimensioning aircraft will be A330. This investment project is estimated to the amount to US$ 250 M. Three main activities have been considered: agro-industry (1.5 B US$, 2009), mineral extraction (0.6 B US$, 2009) and tourism (100 M US$, 2009). The total population of the sub region is about 800 000 inhabitants with an average household annual revenue of 4000 US$ with 5% having a household annual revenue higher than 60 000 US$. Activity production costs and prices, transportation flows, costs and prices levels are available from regional and national statistics.

The resulting optimization problems and inequality sets are of small dimensions, so that standard numerical techniques (linear programming) can be applied to solve them with a limited computational effort [START_REF] Scheel | Mathematical programs with equilibrium constraints: stationarity, optimality et sensitivity[END_REF], [START_REF] Kanafani | Transportation Demand Analysis[END_REF], . This allows to analyze a whole set of differentiated scenarios as well as to perform sensitivity analysis with respect to main exogenous parameters (in general unit prices).

CONCLUSION

The goal of this paper was to propose an impact appraisal tool with respect to the installation of new airport within a region already connected by a multimodal transportation network. One of the main difficulties of this task is related with the estimation of future demand over the modified transportation network which has direct influence on the evaluation of operational conditions, revenues and costs of the new airport.

The suggested approach makes use of two different optimization models: One model is devoted to passenger and freight demand forecasting taking into account the modified accessibility of the multimodal transportation network, the other one defines the global transport supply according with a profit maximization behavior for the involved transport system. The demand forecasting process is based on an entropy maximization approach with distinct urban centers' economic activities generated or attracted constraints levels to determine the intensity and the distribution of origin-destination matrices by the kind of product and transport modality.

The presented solution scheme is composed of an iterative process between the current solution for demand forecasting and the supply optimization problem: the entropy maximizing distribution problem provides the origin-destination matrix given a cost/capacity structure, while the supply optimization problem provides this cost/capacity structure characterizing accessibility, given the origin-destination matrices.

The proposed approach is illustrated with the case of the study of the insertion of a new airport of medium size in a rapid developing rural area. This case study shows that from the points of view of statistics availability and computational complexity, the proposed approach is easy to be performed, allowing the analysis of many different scenarios.

  passenger annual demand from urban center i with respect to people related with activity a, i A a  and revenue level r, R r  . Pop i is the population of center i, i ar  is the percentage of population of center i involved in activity a with revenue level r, I being the set of urban centers considered ( N I  ) and O ar  is the corresponding trip generation factor. Then a possible passenger attraction model is given by:

  a at location i producing attraction at rate i ar  to outside people of revenue level r.

  an a priori distribution and cost that are available for the modality m. Similar notations are adopted for exported goods: am ij

  block time for travel, including departure and arrival activities, between vertices i and j , D m is the average time availability of a vehicle and m F is the fleet size with respect to modality m. SUPPLY OPTIMIZATION MODEL Given a distribution of potential passenger demand [ arm ij T ] and freight demand [ am ij  ], as well as prices for passenger and freight transportation, respectively m ij  and m ij P

  variables are, for each transportation modality and each origin-destination, the vehicle flow levels m ij f , the freight transport capacity am ij  and the passengers transport capacity m ij H . The first part of this maximization is about revenues: the am ij P are unit freight transport fares and the m ij  are unit passenger transportation fares. Then, the second part is about costs: fleet possession cost, transportation cost and terminal cost. The associated constraints are composed of flows positivity (14), flows conservation at transportation nodes (13) and fleet capacity (15), under flow transportation capacity constraints:

CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -5º TRANSPORTATION INFRASTRUCTURE CONFERENCE) August 10 th to 12 th 2011

	São Paulo -Brasil
	05-300
	ISSN 1983-3903
	CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -
	5º TRANSPORTATION INFRASTRUCTURE CONFERENCE)
	August

10 th to 12 th 2011 São Paulo -Brasil

2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -5º TRANSPORTATION INFRASTRUCTURE CONFERENCE) August 10 th to 12 th 2011

	São Paulo -Brasil
	05-300
	ISSN 1983-3903
	CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -
	5º TRANSPORTATION INFRASTRUCTURE CONFERENCE)
	August

10 th to 12 th 2011 São Paulo -Brasil

30) IMPACT EVALUATION FRAMEWORK OF A NEW REGIONAL AIRPORT CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -5º TRANSPORTATION INFRASTRUCTURE CONFERENCE) August 10 th to 12 th 2011

	São Paulo -Brasil
	05-300
	ISSN 1983-3903
	CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -
	5º TRANSPORTATION INFRASTRUCTURE CONFERENCE)
	August

10 th to 12 th 2011 São Paulo -Brasil

Embraer 195 CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -5º TRANSPORTATION INFRASTRUCTURE CONFERENCE) August 10 th to 12 th 2011

	São Paulo -Brasil
	05-300
	ISSN 1983-3903
	CONINFRA 2011 -5º CONGRESSO DE INFRAESTRUTURA DE TRANSPORTES (CONINFRA 2011 -
	5º TRANSPORTATION INFRASTRUCTURE CONFERENCE)
	August

10 th to 12 th 2011 São Paulo -Brasil


Sonia Cafieri
email: sonia.cafieri@enac.fr

Pierre Hansen
email: pierre.hansen@gerad.ca

Lucas Létocart
email: lucas.letocart@lipn.univ-paris13.fr

Leo Liberti
email: liberti@lix.polytechnique.fr

Frédéric Messine
email: frederic.messine@n7.fr

Compact relaxations for polynomial programming problems ⋆

Keywords: polynomial, nonconvex, MINLP, sBB, reformulation, convex relaxation, RLT

Reduced RLT constraints are a special class of Reformulation-Linearization Technique (RLT) constraints. They apply to nonconvex (both continuous and mixed-integer) quadratic programming problems subject to systems of linear equality constraints. We present an extension to the general case of polynomial programming problems and discuss the derived convex relaxation. We then show how to perform rRLT constraint generation so as to reduce the number of inequality constraints in the relaxation, thereby making it more compact and faster to solve. We present some computational results validating our approach.

Introduction

We target Mixed-Integer Nonlinear Programming (MINLP) problems of the form:

min x f (x) g(x) ≤ 0 Ax = b x L ≤ x ≤ x U ∀ i ∈ Z x i ∈ Z,            (1)
where x, x L , x U ∈ R n , Z ⊆ N = {1, . . . , n}, A is a full rank m × n matrix, b ∈ R m , f : R n → R and g : R n → R m ′ are polynomial functions of x.

We describe an extension to polynomial programming of an existing automatic reformulation technique [START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF][START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF][START_REF] Liberti | Compact linearization of binary quadratic problems[END_REF] called reduced Reformulation-Linearization Technique (rRLT). This technique was originally defined only for quadratic problems subject to linear equality constraints. It replaces some of the quadratic terms with suitable linear constraints, called rRLT constraints. These turn out to be a subset of the RLT constraints for quadratic programming [START_REF] Sherali | A new reformulation-linearization technique for bilinear programming problems[END_REF]. The original RLT linearizes all quadratic terms in the problem and generates valid linear equation and inequality cuts by considering multiplications of bound factors (terms like x ix L i and x U ix i ) by linear constraint factors (the left hand side of a constraint such as n j=1 a j x jb ≥ 0 or n j=1 a j x jb = 0). Since bound and constraint factors are always non-negative, so are their products: this way one can generate sets of valid problem constraints. An extension of the RLT to polynomial programming is described in [START_REF] Sherali | A global optimization algorithm for polynomial programming problems using a reformulation-linearization technique[END_REF], and to more general factorable programming problems in [START_REF] Sherali | Global optimization of nonconvex factorable programming problems[END_REF]. These results find their practical limitations in the extremely large number of adjoined constraints. Some heuristic techniques [START_REF] Sherali | A new reformulation-linearization technique for bilinear programming problems[END_REF][START_REF] Sherali | New reformulation linearization/convexification relaxations for univariate and multivariate polynomial programming problems[END_REF] were proposed to help filter out RLT constraints which are redundant. In the rRLT the presence of linear equality constraints in the original problem allows the generation of only those linear RLT constraints that are guaranteed to replace a set of quadratic terms.

We aim to improve performance of spatial Branch-and-Bound (sBB) algorithms targeted at nonconvex NLPs and MINLPs; in particular, the rRLT tightens the lower bound computed by solving a convex relaxation of (1) at each sBB node. We make two original contributions. First, we extend rRLT theory from quadratic to polynomial programs. Second, as rRLT constraint generation depends on an arbitrary choice (the basis of a certain matrix) we show how to choose this basis in such a way as to yield a more compact (i.e., fewer constraints) convex relaxation, denoted by rRLT-C; the rRLT-C relaxation may be weaker than the rRLT one, but experiments show that the loss in tightness is greatly offset by the gain in CPU time taken to solve it. We assume our polynomial programs to be dense up to Sect. 3.2 for simplicity, and deal with sparsity in Sect. [START_REF] Sherali | A new reformulation-linearization technique for bilinear programming problems[END_REF].

Notationwise, we deal sometimes with indexed symbols which are scalars and indexed symbols which are vectors; in order to avoid ambiguities, we denote with boldface all indexed symbols indicating vectors. For example, w ij is a scalar but w ij is the vector (w ij1 , . . . , w ijn ); in line with current optimization literature, we do not use boldface to indicate the name of a whole array, so w might be either a scalar or an array depending on the context.

The rest of this paper is organized as follows: Sect. 2 extends rRLT to polynomial programming. In Sect. 3 we discuss how to construct the rRLT-C compact relaxation. In Sect. 4 we address sparse polynomial programs. Sect. 5 discusses some computational experiments on randomly generated instances.

rRLT for Polynomial Programming

The results presented herein extend [START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF] to the general polynomial case. Let Q = {2, . . . , q}. For each monomial x j1 • • • x jp appearing in the original problem [START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF] where p ∈ Q, we define a finite sequence J = (j 1 , . . . , j p ) and, consistent with the notation introduced by Sherali [START_REF] Sherali | A global optimization algorithm for polynomial programming problems using a reformulation-linearization technique[END_REF], consider defining constraints of the following form:

w J = ℓ≤|J| x j ℓ (2)
(for |J| = 1, i.e. J = (j), we also define w J = x j ). For all p ∈ Q, J ∈ N p and any permutation π in the symmetric group S p we have that w J = w πJ by commutativity. We therefore define an equivalence relation ∼ on N p stating that for J, K ∈ N p , J ∼ K only if ∃π ∈ S p such that J = πK. We then consider the set of equivalence classes N p = N p /∼ to quantify over when indexing added variables w J . In practice, we choose an equivalence class representative for each J ∈ N p which we also denote by J. With a slight abuse of notation, if J ′ ∈ N p ′ and J ′′ ∈ N p ′′ such that p ′ + p ′′ = p and (J ′ , J ′′ ) is in the equivalence class represented by J ∈ N p , we write (J ′ , J ′′ ) = J. We also define, for all p ∈ Q,

M p = 1<p ′ ≤p N p ′ and M 1 p = p ′ ≤p-1 N p ′ .
We multiply the original linear constraints Ax = b by all monomials ℓ≤p-1

x j ℓ and replace them by the corresponding added variables w (J ′ ,j) , where J ′ ∈ N p-1 . This yields the following rRLTS:

∀p ∈ Q, J ′ ∈ N p-1 A w J ′ = bw J ′ , (3)
where w J ′ = (w (J ′ ,1) , . . . , w (J ′ ,n) ). We then consider the companion system:

∀p ∈ Q, J ′ ∈ N p-1 A z J ′ = 0. (4)
Since ( 4) is a linear homogeneous system, there is a matrix M such that the companion system is equivalent to M z = 0, the columns of which are indexed by sequences in M p . We let B ⊆ M p and N ⊆ M p be index sets for basic and nonbasic columns of M . We define the following sets:

C = {(x, w) | Ax = b ∧ ∀p ∈ Q, J ∈ N p (w J = ℓ≤|J| x j ℓ )} (5)
R N = {(x, w) | Ax = b ∧ ∀p ∈ Q, J ′ ∈ N p-1 (A w J ′ = bw J ′ ) ∧ ∀J ∈ N (w J = ℓ≤|J| x j ℓ )}. ( 6
)
Theorem 2.1. For each partition B, N into basic and nonbasic column indices for the companion system M z = 0, we have C = R N .

We remark that for this proof to hold, all possible nonlinear monomials must be present in the problem, which is generally not the case. We address this problem in Sect. 4. A different treatment of the essentially the same concepts, which only employs a bases of A instead of the (larger) companion system, was given in [START_REF] Sherali | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF].

Replacing C with R N for some nonbasis N effectively replaces some nonlinear monomial terms with linear constraints, and therefore contributes to simplify the problem. A convex relaxation for the reformulated problem is readily obtained by applying monomial convexification methods in the literature .

Compact convex relaxation

First, we remark that virtually no practical polynomial problem exhibits all possible nonlinear monomials. Let M = M n and M 1 = M 1 n . We introduce two sets: β ⊆ M indexing all nonlinear monomials appearing in the original problem (1) and β ′ ⊆ M indexing all nonlinear monomials appearing in [START_REF] Liberti | Linearity embedded in nonconvex programs[END_REF] and the rRLTS (3). Reduced RLT constraints are likely to give rise to compact yet tight convex relaxations if N β, in view of the fact that, by Thm. 2.1, only monomials indexed by N need appear in the formulation -so the lower and upper relaxations to monomials outside N can be dropped. Furthermore, the proof of Thm. 2.1 also implies that the number of monomials that can be replaced is equal to the rank ρ of the rRLTS. If the original problem has few monomials, N might not be a proper subset of β, and in practice this occurrence is not rare. Limited to quadratic polynomials, we address this problem in [START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF]; in Sect. 4 we propose a technique to deal with sparse polynomial programs.

In this section we discuss a choice of N whereby the monomial relaxations that are dropped define "large volumes", and are therefore more likely to be dominated by the relaxations of monomials in N . Intuitively, this should yield a compact relaxation whose bound is not far from the normal rRLT relaxation.

Convexity gap

Definition 3.1. Consider a function f : X ⊆ R n → R. Let f (x) be a convex lower bounding function for f and f (x) be a concave upper bounding function for f . Then the set

S = {(x, w) | f (x) ≤ w ≤ f (x)} is a convex relaxation of the set S = {(x, w) | w = f (x)}.
We define the convexity gap V (S) between S and S to be the volume of the set S; namely,

V (S) = x∈X ( f (x) -f (x))dx. ( 7
)
We denote the convexity gap for a quadratic term x i x j with V ij .

Convexity gap for a quadratic term

x 2 i . The convex envelope of the set ξ = {(x i , w ii ) | w ii = x 2 i , x L i ≤ x i ≤ x U i } (where i ≤ n) consists of

the area between the function x 2

i and the chord. The convexity gap of ξ is:

V ii = x U i x L i (x L i ) 2 + (x U i ) 2 -(x L i ) 2 x U i -x L i (x i -x L i ) -x 2 i dx = 1 6 (x U i -x L i ) 3 . ( 8
)
Convexity gap for a bilinear term x i x j . For all i < j ≤ n, the convex envelope of the set

{(x i , x j , w ij ) | w ij = x i x j , x L i ≤ x i ≤ x U i , x L j ≤ x j ≤ x U j } is a tetrahedron ∆ in R 3 given by the McCormick inequalities [9, 10]. The vertices of ∆ are: (x U i , x U j , x U i x U j ), (x U i , x L j , x U i x L j ), (x L i , x U j , x L i x U j ), (x L i , x L j , x L i x L j ). Let µ = x L i -x U i ν = x L j -x U j ξ = x L i x L j -x U i x U j ζ = x L i x U j -x U i x L j D = 0 1 1 1 1 1 0 ν 2 (1 + a 2 ) µ 2 (1 + c 2 ) µ 2 + ν 2 + ξ 2 1 ν 2 (1 + a 2 ) 0 µ 2 + ν 2 + ζ 2 µ 2 (1 + d 2 ) 1 µ 2 (1 + c 2 ) µ 2 + ν 2 + ζ 2 0 ν 2 (1 + b 2 ) 1 µ 2 + ν 2 + ξ 2 µ 2 (1 + d 2 ) ν 2 (1 + b 2 ) 0 .
The volume V ij of ∆ can be computed using the Cayley-Menger formula in 3 dimensions [START_REF] Jiao | Geometrical ambiguity of pair statistics I. point configurations[END_REF], i.e.

V ij = (x L i , x U i , x L j , x U j ) = √ 2 D 24 .
Convexity gap for a multilinear monomial. By [START_REF] Rikun | A convex envelope formula for multilinear functions[END_REF], the convex under-and over-approximating envelopes of a multilinear monomial x 1 • • • x p of degree p, are polyhedral. Therefore, the facet defining inequalities of the enveloping polytope can be computed using the 2 p polytope vertices v ℓ = (x ℓ1 1 , . . . , x ℓp p , j≤p x ℓj j ) in R n+1 , for every p-sequence ℓ ∈ {L, U } p . We carry out this computation using the Porta software [START_REF] Christof | The porta manual page[END_REF]: although its worst-case complexity is exponential in p, it is practically efficient for low values of p. Since our computational experiments only address problems up to degree 4, this methodological choice is appropriate. Standard methods to compute volumes of polytopes exist [START_REF] Avis | User's Guide for lrs[END_REF]. Our preliminary implementation uses the volume of the smallest bounding box -this will be changed later.

If J is the (ordered) sequence of p variable indices appearing in a multilinear monomial µ(x), we let V J denote the convexity gap for µ(x).

Convexity gap for a general monomial. By associativity, for any p ∈ Q and a given sequence J ∈ N p it is always possible to express the monomial µ(x) = ℓ≤p x j ℓ as a product of multilinear factors (

x j1 • • • x j2 ) • • • (x j3 • • • x j4 ).
After replacement by the appropriate added variable, the monomial is reduced to w J1 • • • w J2 , where J 1 = {j 1 , . . . , j 2 } and J 2 = {j 3 , . . . , j 4 }. Associativity can then be re-applied recursively. This allows us to use the results of the preceding sections to derive a convexity gap for µ(x). This approach, albeit simple, is similar, up to commutativity, to the standard reformulation exploited by sBB implementations [START_REF] Smith | A symbolic reformulation/spatial branch-and-bound algorithm for the global optimisation of nonconvex MINLPs[END_REF][START_REF] Adjiman | A global optimization method, αBB, for general twice-differentiable constrained NLPs: I. Theoretical advances[END_REF][START_REF] Sahinidis | BARON 7.2.5: Global Optimization of Mixed-Integer Nonlinear Programs[END_REF][START_REF] Liberti | Writing global optimization software[END_REF][START_REF] Belotti | Branching and bounds tightening techniques for Optimization[END_REF] in view to obtain the convex relaxation of general monomials, and more specifically to the approach followed in [START_REF] Cafieri | On convex relaxations of quadrilinear terms[END_REF] for quadrilinear monomials. Explicit monomial convex envelopes are also known for trilinear terms [START_REF] Meyer | Trilinear monomials with positive or negative domains: Facet s of the convex and concave envelopes[END_REF][START_REF] Meyer | Trilinear monomials with mixed sign domains: Facets of the convex and concave envelopes[END_REF] and univariate terms of odd degree [START_REF] Liberti | Convex envelopes of monomials of odd degree[END_REF].

Since we already treated the multilinear case separately, we assume that µ(x) is not multilinear. In such cases, the recursive techniques described above yield nonlinear convex relaxations. Because our computational results only refer to linear relaxations, however, we take a simpler approach and use interval arithmetic [START_REF] Moore | Introduction to Interval Analysis[END_REF] to compute an interval range [µ L , µ U ] such that µ L ≤ µ(x) ≤ µ U for all x ∈ [x L , x U ]. If J is the (ordered) sequence of variable indices appearing in the monomial µ(x), we compute V J as the volume of the bounding box

[x L , x U ] × [µ L , µ U ].

Choosing a good basis for the companion system

Let B, N be the basic/nonbasic sets of column indices of the companion system (4), which we write in this section as M z = 0, or, equivalently, as

M B z B + M N z N = 0. As shown in Sect. 2, the elements of B, N are sequences J ∈ M. For S ⊆ M and p ∈ Q we define V S,p = J∈S |J|=p V J and V S = p∈Q V S,p . If, for all p ∈ Q, V N,p < V β,p
then the total convexity gap of R N is smaller than that of C. Thus, we aim to find N such that V N,p is minimized, or equivalently, to find B such that V B,p is maximized for all p ∈ Q. This yields the multi-objective problem:

∀p ∈ Q max V B,p M B is a basis of (4) (9)
Next, we show that ( 9) is equivalent to a single-objective problem. Consider a block diagonal m × n matrix Ā with r blocks A s (each a rectangular, full-rank m s × n s matrix for each s ≤ r) having a basis indexed by the set B ⊆ {1, . . . , n}. For all s ≤ r let α s be the set of column indices corresponding to the submatrix A s of Ā. For a matrix T let span(T ) be the space spanned by the columns of T . Lemma 3.2. Let i, j ≤ n be such that i ∈ B, j ∈ B and π be the swap (i, j). If i ∈ α q and j ∈ α t with q = t then the columns of Ā indexed by π B do not form a basis.

We remark that M is a block-diagonal matrix. Instead of showing a formal proof of this fact, which would be long and tedious, we exhibit an example for the case of polynomials of degree at most 3. Example 3.3 does not exploit any specific property of the given matrix, and therefore appropriately illustrates what happens in the general case.

Example 3.3. The companion system (4) M z = 0 derived from the system Ax =

b with A = a 1 a 2 a 3 a 4 a 5 a 6 , x ∈ R 3 with Q = {2, 3} is:                                 a 1 a
     = 0
The number of monomials of n variables of degree exactly p is given in [START_REF] Katzman | Counting monomials[END_REF] as

p + n -1 p . From this, it follows that M has m p∈Q p + n -2 p -1 rows and p∈Q p + n -1 p columns.
Theorem 3.4. Any solution B of (9) maximizing V B also maximizes V B,p for all p ∈ Q.

The single-objective problem max{v B | M B a basis of (4)} has a matroidal structure and can therefore be solved using a greedy algorithm. Evidently, the technique based on compact rRLT constraints is not significant whenever the bounds are the same across all decision variables, because in this case all the V J 's are equal. This, however, is rarely true if bound tightening techniques [START_REF] Messine | Enclosure methods for multivariate differentiable functions and application to global optimization[END_REF][START_REF] Liberti | Writing global optimization software[END_REF][START_REF] Belotti | Branching and bounds tightening techniques for Optimization[END_REF][START_REF] Belotti | Feasibility-based bounds tightening via fixed points[END_REF] are used as a preprocessing step; and it is never true during sBB with rectangular partitioning schemes, as the variable ranges are partitioned at each node. We also emphasize that the remark given in [START_REF] Sherali | Computability of global solutions to factorable nonconvex programs: Part I -Convex underestimating problems[END_REF] p. 11 is also valid in our setting: for polynomial degrees ≥ 3, the proposed convex relaxation might not be monotonically increasing w.r.t. following branches of the sBB tree down from the root, thus preventing the sBB from converging. This can be fixed at each node by adjoining those bound factor inequalities (derived from the multiplication of different bound factors) that involve the branching variable.

Dealing with sparsity

Polynomial problems are rarely dense; this might prevent the set N , introduced in Sect. 2 to index the added variables corresponding to nonbasic columns of the companion system (4), from being a subset of β, the set indexing the added variables corresponding to all monomials appearing in (1). We deal with this possibility by looking for a subset ρ of rows of Ax = b to be multiplied by a subset σ of added variables indexed by M 1 . The rRLTS (3) and its companion system (4) are derived in the same way as for dense polynomial problems.

Theoretically, we should require N β: as remarked in Sect. 3, the rRLTS is then likely to give rise to compact yet tight convex relaxations. In practice, in sparse problems, this requirement will often force ρ or σ to be empty, which means that the reformulation does not take place. Following the principle that every equation corresponds to one fewer degree of freedom, we aim to find ρ, σ such that the size of the sparse rRLTS obtained by multiplying rows in ρ by added variables in σ exceeds the number of new monomials (i.e. monomials not in β) generated by these multiplications. This flexibility in the choice of ρ, σ has a trade-off in terms of relaxation compactness. In order for Thm. 2.1 to hold, the quantifications J ∈ N p in (5) and J ′ ∈ N p-1 in ( 6) should be replaced by J ∈ τ , J ′ ∈ σ respectively, where τ = {(J ′ , j)

| J ′ ∈ σ (|J ′ | = p) ∧ j ∈ N }. Since σ ⊆ M 1 ,
fewer monomials are replaced by the rRLTS, yielding a relaxation which might not be much more compact than the standard RLT relaxation.

We formalize the problem of finding suitable ρ, σ by considering a bipartite graph that represents the incidence of monomials indexed by M in products of rows

a i x = b i of Ax = b multiplied by monomials indexed by M 1 . Let U = {(i, J ′ ) | i ≤ m ∧ J ′ ∈ M 1 } and define a set E such that, for (i, J ′ ) ∈ U and J ∈ M β, {(i, J ′ ), J} ∈ E if there is j ∈ N such that a ij = 0 and (J ′ , j) = J. Consider the bipartite graph G = (U, M β, E): we want to find an induced subgraph G ′ = (U ′ , V ′ , E ′ ) of G, with U ′ ⊆ U and V ′ ⊆ M β, such that: |U ′ | is maximum, |U ′ | > |V ′ | and V ′ = N E (U ′ ) = {v ∈ M β | ∃u ∈ U ′ ({u, v} ∈ E)}.
We define this problem using a Mathematical Programming (MP) formulation: we employ binary variables u i,J ′ for all (i, J ′ ) ∈ U and v J for all J ∈ M β:

max (i,J ′ )∈U u i,J ′ (i,J ′ )∈U u i,J ′ ≥ J ∈β v J + 1 ∀{(i, J ′ ), J} ∈ E v J ≥ u i,J ′ u ∈ {0, 1} |U | v ∈ {0, 1} |M β| .                  (10)
We can then define

ρ = {i ≤ m | ∃J ′ ∈ M 1 (u i,J ′ = 1)} and σ = {J ′ ∈ M 1 | ∃i ≤ m (u i,J ′ = 1
)}. This is a Binary Linear Program (BLP), so it cannot be solved in polynomial time with standard MILP technology (i.e. using a Branch-and-Bound algorithm).

Bipartite matching based algorithm for (10)

We propose an efficient algorithm based on bipartite matching for solving [START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF].

With respect to a matching M of G, a vertex of G is exposed if it is incident to an edge which is not in M . For a subset of edges

F ⊆ E, let U(F ) = {u ∈ U | ∃e ∈ F (u ∈ e)} and V(F ) = {v ∈ M β | ∃e ∈ F (v ∈ e)}.
If M is a maximum matching of G such that there exists an exposed vertex u ∈ U U(M ), the search for an augmenting path from u which is alternating with respect to M will fail by maximality of M [START_REF] Berge | Two theorems in graph theory[END_REF]. A dilation rooted in u with respect to a maximum matching M is a maximal simple alternating path p u (seen as a sequence of edges in E) in G, from u to a vertex u ′ ∈ U(M ), whose even-indexed edges are in M and odd-indexed edges are in E M . Dilations are the certificates used in cardinality bipartite matching algorithms to prove that the current matching is optimal; informally, they certify the failure to find an augmenting alternating path to increase the cardinality of the current matching in the classical bipartite matching algorithm (see Fig. 10-3 in [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF]). In order to deal with the case of isolated vertices, if u ∈ U is isolated, then it is exposed w.r.t. the empty matching, and any empty path p u rooted at u is a dilation; in this case, with a slight abuse of notation, we define U(p u ) = {u}. A dilation set is the set of edges in all dilations rooted at u with respect to M ; dilation sets can be found in polynomial time using breadth first search (bfs) from u (see [START_REF] Papadimitriou | Combinatorial Optimization: Algorithms and Complexity[END_REF], Sect. 10.2). By Lemma 4.2 above, dilation sets P are such that

|U(P )| > |V(P )|.
Algorithm 1 Matching-based algorithm for solving [START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF].

Require: A bipartite graph G = (U, M β, E) Ensure: A nontrivial subgraph G ′ = (U ′ , V ′ ,
G ′ = (U ′ , V ′ , E ′ ) of G such that |U ′ | is maximum, |U ′ | > |V ′ | and V ′ = N E (U ′ )
, or determines that no such subgraph exists.

Computational results

The best practical indication of the importance of rRLT techniques for linear equality constrained polynomial programming is best observed by the computational results given in [START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF], where the presence or absence of rRLT yields differences of up to five orders of magnitude on a class of pooling problems (these are sparse polynomial problems of degree 2) from the oil industry. Those results are limited to quadratic programs and employ a restricted version of the results for sparsity given in Sect. 4. We are currently in the process of extending our code to deal with polynomial programming, so in this paper we only show empirically that the compact rRLT-C formulation generally takes less time to solve and yields bounds that are not much than those given by the rRLT formulation.

We generate two sets of random polynomial programming instances of degree 2,3 and 4 with varying numbers of variables and linear equality constraints (all the monomials are in the objective function, weighted by random scalars). Set 1 consists of 10 instances of degrees 2 and 3 with random variable ranges whose widths follow a Gaussian distribution. Set 2 consists of 8 instances of degrees 3 and 4 with random variable ranges whose widths follow a superposition of two Gaussian distributions (n/2 ranges have width of order 1, the other n/2 have width of order 10); set 2 is designed to simulate the typical sBB node after a few levels of branching, where some of the variable ranges have become small whereas others are still at their original bounds.

For all these instances we construct and solve a linear relaxation with no rRLT constraints (column labeled "simple" in Table 1), the rRLT linear relaxation (column labeled "rRLT") and the rRLT-C linear relaxation (column labeled "rRLT-C"). We recall that rRLT-C is like rRLT without the constraints relaxing monomial terms corresponding of basic columns of the companion system (4). We record bound value and CPU time. Notice some of the generated instances are infeasible: this is consistent with the fact that in a typical sBB search tree some of the nodes represent infeasible subproblems. Since the infeasibility is determined by the linear relaxation, performance on infeasible LPs is also an important factor. These results were obtained using CPLEX 11 [START_REF]ILOG: ILOG CPLEX 11.0 User's Manual[END_REF] on a Pentium Xeon 2.4GHz CPU with 8 GB of RAM running Linux. The results in Table 1 show that rRLT brings considerable benefits to bound tightness within polynomial optimization, and that the (significant) CPU time reduction yielded by rRLT-C is not offset by an excessive loss in bound quality with respect to rRLT: the cumulative bound worsening is 0.07% against a time improvement of nearly 40%. The CPU time taken by the simple relaxation is of course much lower than those of the rRLT relaxations, but already for the quadratic case it was shown in [START_REF] Liberti | An exact reformulation algorithm for large nonconvex NLPs involving bilinear terms[END_REF] that this time difference is not sufficient to offset the benefits of the bound improvement -hence the corresponding values do not appear in Table 1. We remark that rRLT-C does not always yield better CPU time results. This is simply because the relationship between CPU time and number of constraints in solving an LP is neither regular nor monotonic. Conversely, for all p ∈ Q and J ∈ N p we define:

A Appendix: Proofs

z J = w J - ℓ≤p x j ℓ . (11)
The definition of R N implies (3). We consider the companion system (4): since B, N form a partition of the columns of M z = 0 we can write the latter as M B z B + M N z N = 0. By [START_REF] Jiao | Geometrical ambiguity of pair statistics I. point configurations[END_REF] and the fact that the definition of R N imposes the condition ∀J ∈ N (w J = ℓ≤|J| x j ℓ ), it follows that ∀J ∈ N (z J = 0), and hence M N z N = 0. Since B is a basic index set, we have M B z B = 0; since M B is invertible, M z = 0 has the unique solution z = 0. This, in turn, again by [START_REF] Jiao | Geometrical ambiguity of pair statistics I. point configurations[END_REF], implies ∀J ∈ N p (w J = ℓ≤|J| x j ℓ ) as in the definition of C. Hence, the conditions defining R N imply the conditions defining C. ⊓ ⊔ Lemma A.2 (Lemma 3.2). Let i, j ≤ n be such that i ∈ B, j ∈ B and π be the swap (i, j). If i ∈ α q and j ∈ α t with q = t then the columns of Ā indexed by π B do not form a basis.

Proof. Let A ′ be the (square) m × m invertible matrix corresponding to the basis indexed by B. Since span( Ā) = span(A ′ ) and B ⊆ {1, . . . , n}, A ′ is a block-diagonal matrix with r full-rank, square m s × m s blocks A ′ s such that span(A s ) = span(A ′ s ) for all s ≤ r, and such that span(A ′ ) = s≤r span(A ′ s ). Suppose π = (i, j) where i ∈ α q , j ∈ α t for some q = t ≤ r; let A ′′ q be the matrix A ′ q without column i, and let A ′′ t be the matrix A ′ t with column j adjoined. Since A ′ q , A ′ t are square, A ′′ q , A ′′ t are rectangular with span(A ′′ q ) span(A ′ q ) and span(A ′′ t ) = span(A ′ t ). Thus, letting A ′′ s = A ′ s for all s ∈ B {q, t}, we have s≤r span(A ′′ s ) s≤r span(A ′ s ) = span( Ā), which implies that π B does not index a basis of Ā. This concludes the proof.

⊓ ⊔ Theorem A.3 (Thm. 3.4). Any solution B of (9) maximizing V B also maximizes V B,p for all p ∈ Q.

Proof. Let M be the matrix of the companion system (4) and B ′ be a solution of (9) for which, for a given p ∈ Q, we have V B ′ ,p > V B,p . We remark that B ′ can be obtained from B by applying a sequence of t pivots. These correspond to a column index permutation π that can be expressed as the product of transpositions π = k≤t (I k , J k ), where I k indexes an exiting column and J k an entering column of M for all k ≤ t: thus,

B ′ = B {I k | k ≤ t} ∪ {J k | k ≤ t}.
Since M is block diagonal, then as in the proof of Lemma 3.2 we have that M B , M B ′ are block diagonal. Furthermore, by Lemma 3.2, for all k ≤ t we have that the column indices I k , J k are in the same block of M . Furthermore, each block corresponds to the subsystem of (4) given by ∀J ′ ∈ N p-1 A z J ′ = 0 for a fixed p > 1: in other words, each block corresponds to a given p > 1. Thus, for a given k ≤ t, in order for B = (I k , J k )B to correspond to a basis of M , then I k , J k must index columns of M indexed by monomials with the same degree p: this implies that V B,p ′ = V B,p ′ for all p ′ = p. Hence, if V B,p > V

G ′ = (U ′ , V ′ , E ′ ) of G such that |U ′ | is maximum, |U ′ | > |V ′ | and V ′ = N E (U ′ ),
or determines that no such subgraph exists.

Proof. Assume G ′ = ∅. This means that U(M ) = U at the first iteration, which implies that N E (U ) = M β, therefore making G an infeasible instance. Now assume G ′ = ∅; then it is the union of different dilation sets found in different iterations of the while loop, and by Lemma 4.2 is such that |U ′ | > |V ′ |. We have to show that V ′ = N E (U ′ ), so let {u, v} ∈ E with u ∈ U ′ . Suppose v ∈ (M β) V ′ ; then {u, v} was never part of a dilation set at any iteration, and the only way v ∈ V ′ might occur is because u was removed from U at a certain iteration i whilst v was never removed. Consider the iteration (call it i) where u is removed: v is not removed, so v is not incident to any edge in the dilation set P found at itn. i. Let u ′ be the root node of P ; since u is incident to P , there is a dilation p ⊆ P from u ′ to which u is incident (i.e. there is an alternating path p ′ from u ′ to u), which implies that u is incident to M . Since p ⊆ P , v is not in p, which means that {v, u} ∈ M . Now, either v is incident to M or not. If it is, then there is an edge {v, u ′′ } ∈ M , which means that p ′ ∪ {u ′′ , v} is an alternating path longer than p ′ , which can be extended to a dilation using bfs, against the assumption that v is not adjacent to edges in P . If it is not, then v is an exposed vertex; but then p ′ can be extended to an augmenting alternating path from We remark that a more careful running time analysis of Alg. 1 would consider that |U |, |M β| and |E| decrease at each iteration, and would therefore yield considerably lower worst-case complexity order.

Lemma 4 . 1 .

 41 For a dilation p u from u in G w.r.t. a maximum matching M , we have |U(p u )| = |V(p u )| + 1.Lemma 4.2. If p, p ′ are different dilations in E, then |U(p ∪ p ′ )| > |V(p ∪ p ′ )|.

Theorem A. 1 (

 1 Thm. 2.1). For each partition B, N into basic and nonbasic column indices for the companion system M z = 0, we have C = R N . Proof. Since (3) is implied by Ax = b and N ⊆ M p , it follows that C ⊆ R N .

  u ′ to v by simply adding {u, v} to p ′ as the last edge, contradicting the maximality of M . Maximality of |U ′ | follows because all dilations (including isolated vertices) are considered. The worst-case running time is O(|U | min(|U |, |M β|)|E|): a factor O(|U |) for the while loop, a factor O(min(|U |, |M β|)) for updating the data structures (the cardinality of M is bounded above by min(|U |, |M β|)), and a factor O(|E|) in order to find maximum matchings and dilation sets (which are in the data structures used to find the matching in Step 3).

  Proposition 4.3. In polynomial time, Alg. 1 finds a subgraph

	3:	Let M ⊆ E be a maximum matching in G
	4:	if |U(M )| < |U | then
	5:	Find u ∈ U exposed and a corresponding dilation set Pu
	6:	Let H = (U(Pu), V(Pu), Pu)
	7:	Update G ′ ← G ′ ∪ H
	8:	Update G ← G H
	9:	else
	10:	break
	11:	end if
	12: end while
	13: return G ′

E) solving

[START_REF] Al-Khayyal | Jointly constrained biconvex programming[END_REF]

, or ∅ if none exists 1: Let G ′ = ∅ 2: while |U | > 0 do

Table 1 .

 1 Comparing bound strength and CPU time for linear relaxations.

	instance	simple	rRLT	rRLT-C
	name p n m q bound bound time bound time
	1 1 20 10 2 2864.8 2982.63 0.20 2982.63 0.17
	2 1 30 10 2 5286.23 5517.85 2.56 5517.85 1.29
	3 1 10 3 3 48.6115 478.184 0.25 478.184 0.30
	4 1 10 3 3 infeas	infeas 0.23	infeas	0.22
	5 1 20 3 3 infeas	infeas 14.81 infeas 12.92
	6 1 10 5 3 infeas	infeas 0.11	infeas	0.10
	7 1 10 5 3 infeas	infeas 0.16	infeas	0.16
	8 1 20 5 3 infeas	infeas 29.01 infeas 31.41
	9 1 10 7 3 infeas	infeas 0.19	infeas	0.15
	10 1 20 7 3 infeas	infeas 65.22 infeas 46.98
	11 2 10 3 3 130.693 1546.4 0.20 1542.39 0.13
	12 2 10 3 3 18.8459 772.417 0.09 772.417 0.09
	13 2 15 5 3 17.3797 701.723 10.13 701.588 11.03
	14 2 16 8 3 infeas	infeas 488.05 infeas 287.38
	15 2 7 2 4 infeas	infeas 0.17	infeas	0.17
	16 2 8 3 4 47.4445 3468.56 1.26	3458.4 0.85
	17 2 10 3 4 26.2698 4038.69 250.88 4038.68 131.95
	18 2 12 3 4 56.9232 13127.5 166.62 13118.6 109.95

  B,p then V B > V B , contradicting the maximality of B. Induction on t concludes the proof. ⊓ ⊔ Lemma A.4 (Lemma 4.1). For a dilation p u from u in G w.r.t. a maximum matching M , we have |U(p u )| = |V(p u )| + 1.Proof. This follows because a dilation p u is alternating w.r.t. M and starts and ends in U .Lemma A.5 (Lemma 4.2). If p, p ′ are different dilations in E, then |U(p ∪ p ′ )| > |V(p ∪ p ′ )|.Proof. Since all dilations are alternating w.r.t. M , p ∪ p ′ is incident to as at least many vertices of U as of M β; but because p, p ′ are also incident to their exposed root vertices, which are not themselves incident to M , the result follows by Lemma 4.1.

Proposition A.6 (Prop. 4.3). In polynomial time, Alg. 1 finds a subgraph

to Dr. Tatjana Davidović and Christoph Dürr for useful discussions. Financial support by grants: Digiteo Chair 2009-14D "RMNCCO", Digiteo Emergence 2009-55D "ARM" is gratefully acknowledged.


Finding and Proving the Optimum: Cooperative Stochastic and Deterministic Search

Jean-Marc Alliot 1 and Nicolas Durand 1,2 and David Gianazza 1,2 and Jean-Baptiste Gotteland 1,2 Abstract. In this article, we introduce a global cooperative approach between an Interval Branch and Bound Algorithm and an Evolutionary Algorithm, that takes advantage of both methods to optimize a function for which an inclusion function can be expressed. The Branch and Bound algorithm deletes whole blocks of the search space whereas the Evolutionary Algorithm looks for the optimum in the remaining space and sends to the IBBA the best evaluation found in order to improve its Bound. The two algorithms run independently and update common information through shared memory.

The cooperative algorithm prevents premature and local convergence of the evolutionary algorithm, while speeding up the convergence of the branch and bound algorithm. Moreover, the result found is the proved global optimum.

In part 1, a short background is introduced. Part 2.1 describes the basic Interval Branch and Bound Algorithm and part 2.2 the Evolutionary Algorithm. Part 3 introduces the cooperative algorithm and part 4 gives the results of the algorithms on benchmark functions. The last part concludes and gives suggestions of avenues of further research.

Background

Evolutionary Algorithms (EAs) appeared in the 60s with Holland [START_REF] Holland | Adaptation in Natural and Artificial Systems[END_REF] and became popular in the late 80s with Goldberg [START_REF] Goldberg | Genetic Algorithms[END_REF]. They can be very efficient to solve large dimension problems but are difficult to handle (many parameters need to be chosen and are very often problem dependant). They often get trapped in local optima (premature convergence).

Interval Branch and Bound Algorithms (IBBAs) were first introduced by Hansen [START_REF] Hansen | Global optimization using interval analysis[END_REF] in the 90s and combined interval analysis with a Branch and Bound algorithm to reduce the size of the domain containing the optimum. They are able to prove the optimality of the solution but can rarely handle large dimension problems.

According to Alander [START_REF] Jarmo | An indexed bibliography of genetic algorithms: Years 1957-1993[END_REF] who studied the bibliography on genetic algorithms from the 50s to 93, very few articles were related to intervals and none of them dealt with IBBA-EA cooperation. In [START_REF] Jourdan | Hybridizing exact methods and metaheuristics: A taxonomy[END_REF], Jourdan, Basseur and Talbi proposed in 2009 a taxonomy of exact methods and metaheuristics hybridizations. It appears that most of the hybridization between metaheuristics and exact methods concern discrete or combinatorial optimization.

IBBA and EA hybridizations were introduced by Sotiropoulos [START_REF] Sotiropoulos | A new hybrid genetic algorithm for global optimization[END_REF] in 1997 and used by Zhang [START_REF] Zhang | A new interval-genetic algorithm[END_REF] in 2007. Both approaches are integrative combinations, as described by Puchinger and Raidl [START_REF] Puchinger | Combining metaheuristics and exact algorithms in combinatorial: A survey and classification[END_REF]. In Sotropoulos' article, the first step of the algorithm uses a branch and bound to reduce the size of domain to a list of boxes (with a size 1 Institut de Recherche en Informatique de Toulouse, name.surname@irit.fr 2 Laboratoire "Mathématiques Appliqués et Informatique" de l'ENAC smaller than ǫ). Then a genetic algorithm initializes its population in every box and updates the upper bound of the minimum searched. A shrinking box is used to improve the lower bound of the minimum searched. A new population is generated after updating the bounds and the corresponding box list. Zhang incorporates a genetic algorithm in the Interval Branch and Bound algorithm to improve the bounds and the remaining intervals list order.

Our approach is different as the IBBA and the EA cooperate but run independently. They share and update common information that helps both of them to accelerate their convergence.

Standard algorithms

Interval branch and bound

The Interval Branch and Bound Algorithm (IBBA) is basically a Branch and Bound algorithm operating in a search space of intervals. It requires to re-code the function using interval arithmetic [START_REF] Moore | Interval Analysis[END_REF].

Let us consider

I = {[a, b]|a ≤ b, (a, b) ∈ R 2
} the set of compact intervals in R, and I(R) n the set of n-dimensional interval vectors (or boxes). The basic operations of interval arithmetic are defined as follows:

[a, b] + [c, d] = [a + c, b + d] (1a) [a, b] -[c, d] = [a -d, b -c] (1b) [a, b] * [c, d] = [min{ac, ad, bc, bd}, max{ac, ad, bc, bd}] (1c) [a, b]/[c, d] = [a, b].[1/d, 1/c] if 0 / ∈ [c, d] (1d)
The usual real-valued functions (cos, sin, log, and so on) can also be extended to interval arithmetic. There are now a large number of interval arithmetic implementations with various bindings to different languages, such as MPFI [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the MPFI library[END_REF] for C, C++, or SUN interval arithmetic implementation for Fortran 95 or C++ [START_REF]C++ Interval Arithmetic programming manual[END_REF][START_REF]Fortran95 Interval Arithmetic programming manual[END_REF].

In the rest of this document, we shall denote x = (x1, . . . , xn) the real vectors, and X = (X1, . . . , Xn) the interval vectors (boxes). An interval function F : I n → I is said to be an interval extension of the real-valued function f :

R n → R if f (x) ∈ F (X) whenever x ∈ X. An interval function F is said to be inclusion monotonic if X ⊂ Y implies F (X) ⊂ F (Y ).
A result due to Moore ( [START_REF] Moore | Interval Analysis[END_REF], [START_REF] Moore | Methods and applications of interval analysis[END_REF]) states that if F is an inclusion monotonic interval extension of f (or more shortly, an inclusion function), then F (X) contains the range of f (x), for all x ∈ X.

Interval methods for solving optimization problems rely on the above result and the use of deterministic branch and bound techniques to find the optima of f . An initial domain X0 is split into smaller boxes (branching) evaluated using the inclusion function F (bounding). The subregions that do not contain a global minimizer of f are discarded. The basic cut-off test allowing the algorithm to Algorithm 1 Interval branch and bound (maximization)

1: Q ← X0 2: while Q not empty do 3:

Extract X with highest priority pX from Q 4:

if upperbound(F (X)) < f best or width(X) ≤ ǫx or width(F (X)) ≤ ǫ f then 5:

Go to step 3 6:

else 7:
Split X in two sub-boxes X1 and X2

8:

for i ∈ {1, 2} do 9:
ei ← f (midpoint(X i ))

10:

if ei > f best then 11:

f best ← ei 12:
X best ← X i 13:

end if

14:
Insert X i into Q, with priority ei 15:

end for 16:

end if 17: end while 18: return (X best , f best ) discard a subregion X consists in comparing the bounds of F (X) to the best estimator of the optimum found so far. Boxes that are not discarded are divided again in smaller boxes until the desired precision for F (X) is reached (or when X becomes too small). Valid boxes are inserted in a priority queue, which is sorted according to an estimator of the quality of the box. In this simple Branch-and-Bound algorithm, the estimator is just the image of the midpoint of the box. The search stops when the whole domain has been explored (the priority queue is empty).

Several refinements of this basic algorithm exist: monotonicity test when the derivatives of f are available 3 , concavity test, local search procedures for enhancing the best estimator, etc. These procedures may (or may not) be efficient, depending on the problem. In order to keep this article as simple and clear as possible, we opted to use the basic algorithm described above, discarding all these refinements. The interval branch and bound algorithm (IBBA) is described in algorithm 1

Evolutionary algorithm

Evolutionary algorithms, popularized by David Goldberg ([4]) and Michalewicz [START_REF] Michalewicz | Genetic algorithms+data structures=evolution programs[END_REF], are inspired by Darwin's theory of evolution. A population of individuals (points of the search space) is selected according to its fitness, and recombined using crossover and mutation operators. The process is repeated until a termination criterion is met, as described in algorithm 2.

Several refinements have been introduced in this evolution scheme, (among others elitism, scaling, sharing). The reader may refer to [START_REF] Eiben | Introduction to Evolutionary Computing[END_REF] for a description of genetic algorithms and other evolutionary algorithms also inspired from evolutionary theory.

In this article, we have used a real-coded genetic algorithm, where the population is made of N real values vectors. The population is randomly initialized, with uniform probability, within the bounds of the search space. Before selecting the pool of parents, a sigma truncation [START_REF] Goldberg | Genetic Algorithms[END_REF] scaling is applied to the fitness values, followed by a clusterized sharing (step 4). The selection/reproduction itself is made us-Algorithm 2 Evolutionary algorithm (EA)

1: Initialize population 2: while termination criterion is not met do 3:

Evaluate raw fitness of population elements 4:

Apply scaling and sharing operations on raw fitness 5:

Create new population according to new fitness criterion 6:

Replace some elements by mutation and crossover 7: end while 8: Return best elements of population ing the stochastic remainder without replacement [START_REF] Goldberg | Genetic Algorithms[END_REF] principle (step 5).

The crossover and mutation operators are then applied with respective probabilities Pc and Pm (Pc + Pm < 1) to the pool of parents, as follows:

• crossover: two different elements p1 and p2 are randomly drawn from the parents' pool and recombined into two children using an arithmetic crossover. Each child is defined by αp1

+ (1 -α)p2
where α is a real value randomly chosen in a given interval. The process is repeated ⌊ N.Pc 2 ⌋ times to create ⌊N.Pc⌋ children. • mutation: ⌊N.Pm⌋ elements are drawn from the pool of parents.

For each drawn vector, a number k of values is randomly selected, and a Gaussian noise is added to the selected values, thus providing the mutated vector. Assuming the vectors are of dimension n, k is randomly chosen so that k ≤ n. This creates ⌊N.Pm⌋ children.

At the end of the crossover/mutation process, the parents are replaced by their respective children and the new generation of N population elements replaces the previous one. The process is repeated until a termination criterion -maximum time here -is met. We could have chosen other evolutionary algorithms such as Particule Swarm Optimization [START_REF] Kennedy | Particle swarm optimization[END_REF], Differential Evolution [START_REF]Differential Evolution A Practical Approach to Global Optimization[END_REF] or CMA-ES [START_REF] Hansen | Evaluating the cma evolution strategy on multimodal test functions[END_REF]. These algorithms might (or might not) have been more efficient than a real-coded EA. However, the goal of this article is not to find the fastest or most efficient algorithm, but to show how the two approaches (stochastic and deterministic) cooperate. We therefore chose the algorithm we were the most comfortable with.

Parallel cooperative algorithm

When hybridizing the genetic and interval branch and bound algorithms, we adopted the following cooperation scheme. The two algorithms run in parallel. Shared memory is used to exchange information between the two programs. A third thread is used to perform some common operations on elements of both threads.

IBBA thread

The Branch and bound thread is very similar to the Branch and bound algorithm described in section 2.1. The thread is described in algorithm 3.

The main differences between the IBBA algorithm and IBBA thread of the cooperative algorithm are outlined below:

• Shared memory is used to retrieve the best evaluation found by the evolutionary algorithm (step 4). This best evaluation is used to update the bounding value of the IBBA thread, thus speeding up the process of cutting intervals.

f bestag ← GetF romSharedM em(f bestag ) 5: f best ← max(f best , f bestag ) 6:
Extract X with best priority pX from Q

7: if upperbound(F (X)) < f best or width(X) ≤ ǫx or width(F (X)) ≤ ǫ f then 8:
Go to step 6 9:

else 10:
Split X in two sub-boxes X1 and X2

11:

for i ∈ {1, 2} do 12:
ei ← f (midpoint(X i ))

13:

if ei > f best then 14:

f best ← ei 15: bestbb ← midpoint(X i ) 16:
X bestbb ← X i 17:

P utT oSharedM em(bestbb) Other operations are performed on the priority queue of the IBBA thread by the UPDATE thread at the synchronization point. They are described in section 3.3.

EA thread

The evolutionary algorithm thread is also very similar to the evolutionary algorithm described in section 2.2. This thread is described in algorithm 4. P utT oSharedM em(f bestag )

6:
bestbb ← GetF romSharedM em(bestbb)

7:
Replace worst population element by bestbb 8:

Evaluate bestbb raw fitness 9:

Apply scaling and sharing operations on raw fitness 10:

Create new population according to new fitness criterion 11:

Replace some elements by mutation and crossover 12: end while 13: Return best element of population The main differences are outlined below:

• The EA thread puts in shared memory the best evaluation found so far (step 5), which will be retrieved by the IBBA thread. • The EA thread gets from the shared memory the best element found so far by the IBBA thread (step 6) and then replaces its worst population element by this element.

Other operations are performed by the UPDATE thread on the EA population at the synchronization point (step 3). These operations are described in section 3.3.

UPDATE thread

The UPDATE thread is triggered every t seconds. It is described in algorithm 5. Sleep for duration t 3:

Wait for and then Suspend EA thread and IBBA thread 4:

for i = 1 to N do 5: dmin ← +∞ 6: N Q ← Q 7:
while N Q not empty and dmin = 0 do 8:

Extract (X, pX ) from N Q

9:
if upperbound(F (X)) < f best then 10:

Suppress X from Q 11:

else 12:
if elt(i) ∈ X then if distance(elt(i), X) < dmin then 16:

dmin ← distance(elt(i), X) Resume EA thread and IBBA thread 31: end loop

The thread first waits for the IBBA and the EA thread to reach their synchronization point, and suspends them before performing any operation.

The thread then examines in turn the N elements of the population of the EA thread. For each element elt(i), it performs a lookup in the priority queue Q of the IBBA thread. This queue contains all the interval vectors (boxes) of search space that are still valid. For each element elt(i), the thread finds the minimal distance dmin of this element to the closest box Xc in queue Q (in the process the thread also suppresses from Q boxes whose upper-bounds are lower than the current best evaluation in step 10). Then:

• if dmin is equal to zero, then we have found a box X that contains elt(i) and elt(i) is in an admissible zone of search space. Thus elt(i) is kept inside the EA population. If f (elt(i)) is better than the current priority p X of box X that contains elt(i) then we have found a better estimator for the maximum in box X, and the priority of box X in queue Q is updated to f (elt(i)). • if dmin is not zero then elt(i) is outside the admissible search space. Then we project elt(i) on the closest box Xc and replace in the EA population elt(i) by this projection.

The projection algorithm is simple and described in algorithm 6.

Algorithm 6 Projection algorithm (step 27 of algorithm 5)) 1: for j = 1 to n do if elt(i)(j) ∈ Xc(j) then if upperbound(Xc(j)) < elt(i)(j) then 4:

elt(i)(j) ← upperbound(Xc(j))

5:

else 6:
elt(i)(j) ← lowerbound(Xc(j))

7:

end if

8:
end if 9: end for elt(i) is a real vector in R n , while Xc is an interval real vector in I(R) n . For each dimension j we check if elt(i)(j) is inside interval Xc(j). If elt(i)(j) is not inside the interval then we replace elt(i)(j) by the closest element of interval Xc(j), which is either the upper bound or the lower bound of Xc(j).

The UPDATE thread has two main goals:

1. Put all the population elements of the EA thread back into the admissible search space. This will increase the speed of convergence of the EA, and will also take the EA out of local minima as soon as these minima have been ruled out by the IBBA thread. In fact, on some examples developed in section 4 we will see that even the best element of the EA thread can be suppressed and projected elsewhere by the UPDATE thread when this element is inside a local optimum. 2. Re-sort the IBBA priority queue, thus focusing the search in the IBBA thread on the "interesting" part of the search space, and increasing the IBBA convergence speed.

The UPDATE thread is a costly one, especially when there are many boxes in the priority queue Q. Thus, it should not be triggered too often, but often enough to fulfil its two goals. For simplicity's sake, we have only presented here a simple strategy (timer interval) for triggering this thread. But other, more efficient strategies can be used to trigger it, based on the size of the priority queue, the evolution of the population in the EA thread. Moreover, some implementation tricks can be used to accelerate it. However, again for simplicity's sake, we present in the following results the simple, basic algorithm.

Understanding the algorithm

In this section we are going to graphically present a few examples in order to understand how the cooperative algorithm works. Statistical tests and results will be presented in section 4.

We will first consider the Griewank function in dimension 6. Griewank is a classical example [START_REF] Pohlheim | Example of objective functions, documentation of the Matlab Genetic and evolutionary algorithm toolbox[END_REF], even if not a very good one regarding global optimization, as Locatelli has shown in [START_REF] Locatelli | A note on the Griewank test function[END_REF] that the function becomes easier to optimize for large dimensions with stochastic algorithms. Moreover the Griewank function is partially separable, which makes convergence of both EA and IBBA algorithms extremely fast.

It is now customary to use a variant of the Griewank function, the rotated Griewank function [START_REF] Suganthan | Problem definitions and evaluation criteria for the CEC 2005 special session on real-parameter optimization[END_REF].

f (x) = D i=1 zi 4000 - D i=1 cos( zi √ i ) + 1 with z = M (x -o)
where M is a random rotation matrix and o a random vector. To have results easy to read we maximize here the function g(x) =

1 1+f (x)
The rotated Griewank function is not separable . The nonseparability of the variables turns the inclusion function of the IBBA, which is very efficient for the regular Griewank function, into a very inefficient one. It is currently impossible to find the optimum of the R-Griewank function with a simple IBBA algorithm as soon as D is larger than 7. Thus 6 is a good value to see how the cooperative algorithm works, and to compare the convergence of all three algorithms.

On Figure 1, we first compare the cooperative algorithm with the standard Evolutionary Algorithm and with the Branch and Bound algorithm 4 . These results are only an indication of the general be- haviour of the cooperative algorithm, and statistical results will be presented in the next section. It is already clear however that the cooperative algorithm is much faster than both the EA and the IBBA algorithms, while proving the result, as the IBBA does.

On Figure 2, we see how the cooperative algorithm finds and proves the optimum in 25s. The red line is the value of the internal best evaluation found by the EA thread. The green line is the internal value of the best evaluation found by the IBBA thread. Deep blue crosses are the times when the EA thread sends to the IBBA thread a better evaluation than the one that the IBBA has. The pink crosses are the times when the IBBA thread sends to the EA thread a better element than the one the EA thread has. The light blue crosses are the times when the UPDATE thread destroys the best element of the EA thread because it is outside the searchable domain (the EA thread is stuck in a local optimum). We can see on this figure that the algorithms collaborate in an extremely efficient way. All mechanisms are used during the run.

With 8 variable, the IBBA algorithm can never find a solution in a reasonable amount of time, while the cooperative algorithm can. The EA algorithm performance on the 8 variable function depends on luck as it can sometimes get stuck in a local optimum (the cooperative algorithm never does). On Figure 3, we have a detail of the convergence of the cooperative algorithm on the 8 variable rotated Griewank. We choose one of the example with the longest running time, in order to clearly see what happens here. The EA algorithm is usually dragging up the IBBA (deep blue crosses), at least at the beginning of the search. However, from 40s and up to 45s, it is the IBBA which is taking the EA algorithm out of local minima. After 45s, the algorithm is already in the vicinity of the optimum. The IBBA is more efficient than the EA in performing a local optimization (all the pink crosses at the top of the graphic). The light blue crosses at the top are simply the result of the IBBA thread "killing" search space at a very fast pace. The global optimum is found with the required precision at the last pink cross (85s). Thereafter, the algorithm is just completing the proof by searching and cutting the remaining search space (the last clear blue crosses).

As a last example we will discuss the Michalewicz function [START_REF] Pohlheim | Example of objective functions, documentation of the Matlab Genetic and evolutionary algorithm toolbox[END_REF]:

f (x) = D i=1 sin(xi) sin( ix 2 i π ) 20
This function is difficult to optimize because of the steepness of the curve (the 20-th power), and is interesting because there are very few results available for large D. In [START_REF] Pohlheim | Example of objective functions, documentation of the Matlab Genetic and evolutionary algorithm toolbox[END_REF], the highest D for which the optimum is presented is D = 10, and the same goes for [START_REF] Molga | Test functions for optimization needs[END_REF] and [START_REF] Idoumghar | A novel hybrid evolutionary algorithm for multi-modal function optimization and engineering applications[END_REF].

Of course, the optimum is never proved as it is found by stochastic algorithms. It was thus a challenge to find and prove the optimum of the Michalewicz function for D = 12 variables. The function optimized is g(x) = f (x) + D, in order to keep g positive. On Figure 4 we see how the cooperative algorithm behaves. The cooperation is present all the way up to the optimum. On Figure 5, we have displayed the times when the IBBA thread had to kill the best element of the EA thread which was stuck in a local optimum. This happens often because the function has a lot of local optima and because the vicinity of the optimum is extremely small due to the steepness of the function.

Statistical tests and results

In Table 1, statistical results on the rotated Griewank function are presented on 100 runs using the EA alone, the IBBA alone, and the cooperative algorithm. The search space is [-400, 600] n with n ∈ {6, 7, 8, 9, 10}. The time was limited to 1800 seconds. We used a 3.40GHz Intel Xeon E3-1270.

The Evolutionnary Algorithm (described in part 2.2) parameters are: 1000 chromosomes, Pc = 0.5 and Pm = 0.3. An arithmetic crossover is used with α ∈ [-0.5, 1.5]. The mutation operator adds a random noise in the [-0.5, 0.5] interval to each variable of the function. We set σ = 2 for sigma truncation scaling and used the clusterized sharing described by Yin and Germay [START_REF] Yin | A fast genetic algorithm with sharing scheme using cluster analysis methods in multimodal function optimization[END_REF]. The algorithm stops when the distance between the current best element and the optimum (1 in this case) is less than 10 -4 , or when the alloted time (1800s) is over.

For the Interval Branch and Bound Algorithm, ǫx and ǫ f (see algorithm 1) were set to 10 -2 and 10 -4 . The algorithm stops when the Q list is empty, or when the alloted time is over.

The same parameters and stopping criteria are used for the cooperative algorithm. For each algorithm, Table 1 gives the number of runs that found the optimum in less than 1800 seconds, the mean time duration, and the corresponding standard deviation in seconds.

Results show that the IBBA can only deal with small dimensions (≤ 6) in a reasonable time. The EA approach is sensitive to dimension as well. The EA would certainly give much better results if its parameters and operators were optimized for the Griewank function but we did not concentrate on this issue. The Cooperative Algorithm always gives much better results than the IBBA and EA.

Conclusion

In this article, we have presented a cooperative algorithm that combines the advantages of the global stochastic optimization techniques and the global deterministic techniques, and we have shown that this algorithm is able to speed up the convergence of the stochastic algorithm. But the most important result is that this algorithm is able to prove the optimality of the result for very difficult functions such as R-Griewank or Michalewicz, up to 12 variables, while the best available result was, as far as we know, currently limited to 6 to 8 variables.

We have also focused on presenting the algorithm as clearly as possible, using only a standard evolutionary algorithm and a standard interval branch and bound algorithm, leaving out all the acceleration, modifications and implementation tricks. The results presented are thus easily reproducible with off-the-shelf algorithms. For simplicity's sake and lack of space, we have also limited our presentation to two functions but we have run similar tests on many more functions (Rastrigin, Schwefel, etc. . . ) with similar excellent results.

We think that this cooperative algorithm is currently the best algorithm available for proving the optimality of the result for complex and deceptive functions up to a number of variables which had, to our knowledge, never been reached.

Our next paper will present the modifications of the EA and the IBBA algorithms, along with the implementation optimizations that we have developed. These improvements speed up the cooperative algorithm by a factor of up to 50 and enable for example to find and prove the optimum of the Michalewicz function with 12 variables in less than 150 seconds on a "standard" dual core processor.

Algorithm 4

 4 Cooperative algorithm, EA thread 1: Initialize population 2: while (termination criterion not met) or (no signal from IBBA thread)

Algorithm 5

 5 Cooperative algorithm, UPDATE thread 1: loop 2:

Figure 1 .

 1 Figure 1. Comparison Cooperative/AG/BB (6 variable Griewank)

Figure 2 .

 2 Figure 2. Understanding cooperative (6 variable Griewank)

Figure 3 .

 3 Figure 3. Understanding cooperative (8 variable Griewank)

Figure 4 .

 4 Figure 4. Understanding cooperative (12 variable Michalewicz)

Figure 5 .

 5 Figure 5. Killing local minima (12 variable Michalewicz)

  size

  When the IBBA thread finds a better overall element, it updates the shared memory, and makes this element available for the EA thread (step 17).• When the IBBA thread ends, we are sure that we have found a global optimum and the IBBA thread sends a signal to the EA thread and then terminates (step 23).

	18:	end if
	19:	Insert X i into Q, with priority ei
	20:	end for
	21:	end if
	22: end while
	23: Signal EA thread and stop
	•

Table 1 .

 1 Rotated Griewank function, statistical results on 100 runs

		6	7	8	9	10
	EA Found 100	94	92	83	15
	Mean 204 864 972 1340 1678
	Sigma	92 356 389	430	34
	IBBA Found	71	0	0	0	0
	Mean 284
	Sigma 192
	Cooperative Found 100 100 100	100	100
	Mean	50	62 156	215	267
	Sigma	18	47	85	317	105

This can be done by hand for simple functions. or using automatic differentiation[START_REF] Bücker | Automatic Differentiation: Applications, Theory, and Implementations[END_REF] for complex programs.

In all figures, AG+BB=cooperative algorithm, AG=evolutionary algorithm, BB=branch and bound. The x-axis is in seconds, the y-axis is the fitness.