Marcella Autran
email: marcella@pet.coppe.ufrj.br

Burlier Drummond

Catherine Mancel
email: catherine.mancel@enac.fr

Amaranto Lopes Pereira
email: amaranto@pet.coppe.ufrj.br

Felix Mora-Camino
email: felix.mora@enac.fr

Félix Mora-Camino

An optimization approach for long planning of multimodal freight transportation systems

Keywords: Planning, Network Flows, Bi-level Programming

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

Introduction

In this communication is considered the problem of long term forecasting of the freight traffic growth in a large multimodal transportation system. This problem is crucial when planning the necessary investments in terminals, networks and fleet resources. One of the main difficulty of this task is related with the estimation of future demand over the network which has direct influence on the operational conditions and costs. The proposed approach makes use of two different optimization models: One model is devoted to freight demand forecasting, the other one defines the transport supply according with a profit maximization behavior for the whole sector operating in this area. The supply optimization model considers simultaneously two classes of flows: vhicle flows providing transportation capacity and passengers flows generating revenues to the operators. Hence, no classical flows in networks optimization technique is available to solve this problem while two level solution techniques considering vehicle flows at the first level and freight flows at the other can be considered. Then a global optimization problem is associated to demand and operations costs. Each of the two optimization problems, taken separately is convex, however, the whole problem, through the definition of revenue and of global cost constraints is non convex.

A proposed solution schemes is composed of an iterative process between the current solutions of the demand and the supply optimization problems: the entropy maximizing problem provides the freight origin-destination matrix given a fare structure, while the supply optimization problem provides the fare structure given a freight origin-destination matrix. To enforce convergence while maintaining convexity of the two problems, the effective-potential demand level constraints of the supply optimization problem are modified.

Distribution of Demand

To perform the prevision of demand, it is supposed that an a priori trip distribution is available. This information may be the result of a prediction of trip distribution for a previous time period or the result of local forecasting studies. It is also supposed [START_REF] Assad | Multicommodity network flows-a survey[END_REF] that predictions of generation and attraction potential levels

{ } N i O i ..., , 1 , ∈ and { } N j D j ..., , 1 , ∈
, are available.

Local cost elasticities

We consider that the demand for a given origin destination i-j is such that there exists a positive constant ij λ with :

) ( 0 0 ) ( ) ( ij ij ij e T T ij ij ij ij π π λ π π - - = (1)
where ij π is a mean cost between i and j. Considering that the potential origin and destination levels

∑ ≠ = = = N i j j i ij N à i O T , 1 1 ) 0 ( (2)
∑ ≠ = = = N j i i j ij N à j D T , 1 1 ) 0 ( (3)
we have the following constraints:

∑ ≠ = = = N i j j i ij N à i O e T ij , 1 1 π λ (4) ∑ ≠ = = = N j i i j ij N à j D e T ij , 1 1 π λ (5)
Observe that a more general relation for (1) could be:

) ( ) ( ) ( ) ( 0 0 ij ij ij ij ij ij ij ij f T f T π π π π = (6)
where ij f is a continuous positive monotonous decreasing function such as:

1 ) 0 ( = ij f and 0 ) ( lim = +∞ → x f ij x (7)

Estimation of demand distribution

The demand distribution estimation problem is taken here as a constrained entropy maximization problem. Then to a choice of an

instance I={[ ij T ˆ], { } N i O i ..., , 1 , ∈ , { } N j D j ..., , 1 , ∈ , ]
[ ij λ } is associated the following maximization problem, Problem I: [START_REF] Brotcorne | A bilevel model for toll optimization on a multicommodity transportation network[END_REF] under constraints (4), [START_REF] Alou | A multilevel modelling approach for air transportation system[END_REF] and

( ) ∑∑ =≠ = - N i N i j j ij ij ij T T Ln T 1, 1 / max
j i N à j N à i T ij ≠ = = ≥ , 1 , 1 , 0 (9)
The adopted optimization criteria, a conditional entropy function, is representative of the global distortion between the a priori and the predicted demand distributions. The instance is said to be consistent if the following condition is fulfilled :

∑∑ == = = N i N j j i T D O 11 (10)

Analysis of the solution

The above optimization problem being convex it is useful to introduce the Lagrangian associated to this problem :

( ) ∑∑ ∑∑ ∑∑ =≠ = =≠ = =≠ = - + - + - = N j N j i i j ij j N i N i j j i ij i N i N i j j ij ij ij D e T O e T T T Ln T L ij ij 1, 1 1, 1 1, 1 ) ( ) ( / π λ π λ β α (11)
The first order optimality conditions are such as :

0 = ∂ ∂ ij T L j i N à j N à i ≠ = = , 1 , 1 (12)
or :

j i N à j N à i ij e j ij ij e i ij T ij T Ln ≠ = = = + + - - , 1 , 1 0 ) / ( 1 π λ β π λ α (13)
Then the solution is such as :

j i N à j N à i ij ij j i e ij ij e T T ≠ = = + + - = , 1 , 1 ) ) ( 1 ( * ˆπ λ β α (14)
In the case that ij ij π λ is very small, this expression can be approximated by:

j i N à j N à i ij ij j i e T T ij ij ≠ = = + + + - = , 1 , 1 )) 1 )( ( 1 (
* ˆπ λ β α (15)
The dual variables are such as:

N to i i O j i ij p N i j j ij 1 ) ( , 1 = = + ∑ ≠ = β α θ ( 16
) N to j j D j i ij p N j i i ij 1 ) ( , 1 = = + ∑ ≠ = β α θ (17) where ) 1 ( ˆ- = ij ij e ij T ij π λ θ and ij ij e ij e p π λ = (18)
however, the numerical solution of an instance of Problem I will be obtained using linear -convex algorithm such as the Simplex-Convex or Frank-Wolfe algorithm [START_REF] Assad | Multicommodity network flows-a survey[END_REF]. Another interesting direction for numerical resolution is through the solution of the geometric primal associated to this problem [START_REF] Mora-Camino | Introduction à la Programmation Géométrique[END_REF].

Freight Supply Model

The proposed model takes into account two types of flows and their corresponding constraints [START_REF] Berge | The theory of graphs[END_REF], [START_REF] Ford | Flows in networks[END_REF]: vehicle flows along the network according to the available fleets and freight flows using the resulting transport capacity [START_REF] Dempe | Foundations of bilevel programming[END_REF]. The fleet of vehicle and its operation generate fixed and variable costs, while freight flows are the main source of revenue for the freight transportation sector.

Vehicle Flow Network Model

Considering that freight fleet is composed of M different vehicles classes. To each freight terminal i, i∈A, and to each freight modality is associated the set of freight terminals directly reachable with this modality : A C (i). Then M freight transportation networks can be introduced:

R m = [G m ,[ m ij f ]] { } M m , , 1 L ∈ (19)
where G m is the graph [A, Γ m ] where the successor function Γ m is such as :

Γ m (i) =A m (i) i∈ A (20) [ m ij f
] is the flow of freight vehicles over G m . It satisfies to conservation and positive ness constraints:

m ik N i j j N i k k m ji f f ∑∑ ≠ =≠ = = , 1, 1 { } M m , , 1 L ∈ (21) 0 ≥ m ij f , i, j ∈A { } M m , , 1 L ∈ (22)
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 air transportation network. However for sake of realism, fleet capacity constraints are introduced:

m m m ij N ii A j m ij F D d f m ≤ ∑∑ =∈ 1) ( (23)
where m ij d is the block time for a freight vehicle of class m to go from terminals i and j , D m is the medium time availability of a vehicles of class m. It can be also of interest to introduce terminal capacities such as:

A i K f f i i j m ji dm i j m ji am m m ∈ ≤ + ∑ ∑ Γ ∈ Γ ∈ - ) ( ) ( 1 ω ω (24)
where i K is the capacity at terminal i j. Here am ω and dm ω are capacity parameters attached to the different types of operations (arrival or departure) and fleets.

Freight Flows Modeling

The M freight networks provide the physical support for the freight flows network defined as :

R f = [G f ,[ ij φ ]] (25)
where

G f = [A, Γ f ] and Γ f is such as: Γ f (i) = { } ) ( , , 1 i A m M m L U ∈ i∈A (26)
[ ij φ ]is the freight flows over this network. It obeys to the following capacity constraints :

A j i f C m ij m M m ij ∈ ≤ ≤ ∑ = , 0 1 φ (27)
Within graph G f , for each origin-destination pair, a set of concurrent paths composed of a succession of trips is retained according to directness criteria:

Ch ij ={ Ch n ij , n =1 to n max } i, j ∈ A ( 28
)
Here it is considered that freight firms assign demand for a given O-D pair between these different possible paths according with their available capacity. This approach is acceptable when available capacity is close to demand levels.

Let [α(i,j,n,k,l) ] be the incidence matrix between path Ch n ij and arc (k,l) of

G pax : α(i,j,n,k,l) = 1 if (k,l)∈ Ch n ij α(i,j,n,k,l) = 0 else.
Then the freight flow between terminals i and j is given by : ( )

j i A j i j i n l k k N k N k l l N n ij ij ≠ ∈ = ∑∑ ∑ =≠ == , , ) , , , , ( ln 1, 11 θ α φ (29)
where ijn θ is the flow of freight between terminals i and j using the n th path between them. The mean fare for the i-j origin-destination pair is then given by :

A j i ijn ij N n ijn ij N n n ij ij ∈ ∑ = ∑ = = , ) 1 /( ) 1 ( θ θ π π (30)

Supply Optimization Model

Given a distribution of potential demand

] [ * ij T , associated to [ ij T ] { } N i O i ..., , 1 , ∈ , { } N j D j ..., , 1 , ∈
, and ] [ ij λ , as well as to mean air transportation fares ] [ ij π , the optimization of the supply (capacities and fares) by the freight firms can be considered:

Problem II is concerned with the optimization of the global economic performance over a future period of time of the firms operating the freight transportation network.

Here the decision variables are the vehicle flows (

{ } M m f m ij , , 1 , L ∈
) between the different terminals and the fares ( n ij π ) applied to each selected path between the terminals. To solve this problem it is also necessary to introduce the effective flows of freight ( ijn θ ) associated to each selected path between the terminals.

The optimization criterion of Problem II is given by:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + - + ∑∑ ∑ ∑∑∑ =Γ ∈ = =Γ ∈= ) ( ) 1 /( ) ( max 1 1 11 , , N ij m ij m ij m mv mf M m r N ij ijn n ij N n f c pax ij ij n ij ij f c F c c τ θ π θ π (31)
where τ is a return rate and the mv c are coefficients related with fixed and variable fleet and vehicle flows costs.

Problem II must satisfy the following constraints:

m ik N i j j N i k k m ji f f ∑∑ ≠ =≠ = = , 1, 1 { } M m , , 1 L ∈ i, j ∈A (32) m m m ij N ii A j m ij F D d f m ≤ ∑∑ =∈ 1) ( (33) m ij m M m k N k N k l l N n f C j i n l k kl ∑ ∑∑ ∑ = =≠ == ≤ 1 ln 1, 11 ) , , , , ( θ α i, j ∈A (34) ∈ ≤ ∑ = j i T ij N n ijn ij , * 1 θ A , i ≠ j (35)
with the positiveness conditions: Each of the two optimization problems, taken separately is a convex programming problem, however, the whole problem when integrating Problem I into Problem II through constraint ( 35) is non convex. Then it seems interesting to solve numerically these problems separately and design an interactive process towards equilibrium. In this case, some questions are of interest: -the easiness to solve numerically each of the optimization problems, -the guarantee of convergence of the iterative process, -the speed of convergence, -the quality of the limit solution. These two linked problems constitute a non standard bilevel programming problem [START_REF] Alou | Pricing in air transport systems : a multilevel approach[END_REF] where Problem II is the leader's problem while Problem I is the follower's problem. This results from the adoption of a deregulated point of view in which firms operate the freight transport network according to their direct economic interest without considering any social surplus. Constraint (35) plays a central role in the articulation of the two problems and the convergence of their solutions towards a common global solution. Since this constraint transmits to the revenue optimization problem the reaction of demand with respect to changes in mean fares between the different origins and destinations, it is useful to make apparent this effect so that limited fare values will be provided by the solutions of Problem II. However, the relation between origin-destinations flows ] [ * ij T and fare levels ] [ ij π is quite complex and at least non linear. So, to maintain the convexity of Problem II, this constraint is replaced by its first order approximation where the reference values are the solutions of the two problems at the previous iteration:

0 ln ≥ k θ k, l ∈ A, n ∈ N kl (36) { } M m f m ij , , 1 0 L ∈ ≥ i, j ∈ A (37) ∈ ∈ ≥ j i N n ij n ij , , 0 π A (38)

Global Solution Scheme

j i A j i ≠ ∈ , ,
The partial derivatives present in relation (39) can be computed by differentiating relations (15), ( 16),(17) and (18) at the solutions of Problem I and II at iteration k-1 and solving linear systems of equations. The following inequalities are guaranteed:

0 1 * ≤ - ∂ ∂ k n ij ij T π , 0 1 * ≥ - ∂ ∂ k n lm ij T π 1 * 1 * - ∂ ∂ << - ∂ ∂ k n ij ij T k n lm ij T π π (40) j i A j i ≠ ∈ , , , m l A m l ≠ ∈ , , , ) , ( ) , ( j i m l ≠
So it appears that starting from low fares, an increase of fares on trips linking origin i and destination j at solution of Problem II at iteration k-1 will imply a decrease of potential demand on the same origin-destination pair at the same iteration of Problem I and then (40)

Problem II ] [ * ij T ] [ ij π [ m ij f ] [ ln k θ ] Problem I θ π θ π c - , max under (42)} is : ⎩ ⎨ ⎧ - = + - + = ) exp( 0 ) 2 /( ) 1 1 ( k k c k k π λ θ θ λ λ π λ π (43)
It is clear that the limit when k →+ ∞ of (π k , θ k ) is ((1+c λ)/λ, θ 0 e -(1+c λ ). Here, the convergence rate is such as:

) 2 1 1 ( ) 2 / 1 ( - - - - = - k k k k π π π π (44)
Other bilevel schemes have been considered in [START_REF] Alou | Pricing in air transport systems : a multilevel approach[END_REF] and [START_REF] Brotcorne | A bilevel model for toll optimization on a multicommodity transportation network[END_REF] for different transportation problems while numerical convergence conditions have been discussed in [START_REF] Scheel | Mathematical programs with equilibrium constraints: stationarity, optimality and sensitivity[END_REF].

Conclusion

This communication has considered the problem of long term forecasting of the flows and traffic growth in a freight transportation network. The proposed approach has introduced two different optimization models: One model devoted to demand forecasting and the other describing a profit maximization supply behavior by freight firms. An entropy maximization approach is used to determine origindestination matrices. The formulation of this problem introduces in a new way elasticity of demand with respect to fares. The supply optimization model considers simultaneously two classes of flows: vehicle flows providing freight transportation capacity and goods flows generating revenues to the freight operators. Then a global optimization problem is associated to each scenario with respect to demand and operations costs. The proposed solution scheme is composed of an iterative process between the current solutions of the demand and the supply optimization problems: the entropy maximizing problem provides the passengers origin-destination matrix given a fare structure, while the supply optimization problem provides the fare structure given a passengers origin-destination matrix. Convergence conditions are discussed for this iterative process between two problems which can be seen as inverse of each other.

  It appears that Problem I and Problem II are strongly interdependent: while Problem I provides to Problem II potential levels of demand ] [ * ij T (constraint (35)), Problem II provides mean fare values ] [ ij π (relation (30)) to Problem I.

Figure 1 :

 1 Figure 1: Interaction between problems


Maxime Maria
email: maxime.maria@univ-poitiers.fr

Sébastien Horna
email: sebastien.horna@univ-poitiers.fr

Lilian Aveneau
email: lilian.aveneau@univ-poitiers.fr

Topological Space Partition for Fast Ray Tracing in Architectural Models

Keywords: Interactive Ray Tracing, Beam Tracing, Acceleration Structure, Architecture, Topology, Cells-and-Portals

come

INTRODUCTION

Ray-tracing is a rendering technique allowing to compute high-quality realistic images. It allows to render all kinds of visual phenomena such as reflection, refraction, direct or global illumination. Its main disadvantage lies in its high computational cost. Indeed, for each pixel to render, at least one ray is cast through a virtual scene. This simulates the light transport through the scene, and returns the color of the pixel. For a given ray, the problem consists in finding the nearest intersection, among the geometry of the scene. From Whitted's work [START_REF] Whitted | An improved illumination model for shaded display[END_REF] to nowadays, many methods have been proposed to improve ray tracing efficiency. Usually, a well-suited acceleration structure is used to reduce the number of useless ray intersection tests.

As part of architectural environments, lighting simulation can be useful, for instance to visualise a building before its construction. General acceleration structures such as kd-trees [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF], BSPtrees [START_REF] Fuchs | On visible surface generation by a priori tree structures[END_REF], regular grids [START_REF] Fujimoto | Arts: accelerated ray-tracing system[END_REF] or bounding volume hierarchy [START_REF] Rubin | A 3-dimensional representation for fast rendering of complex scenes[END_REF][START_REF] Kay | Ray tracing complex scenes[END_REF] are not fully adapted to architectural scenes. While they can be used for a broad range of applications, they can produce bad rendering results, particularly with architectural scenes.

Usually, buildings are rendered using a specific acceleration structure, called cells-and-portals. In previous works [START_REF] Airey | Towards image realism with interactive update rates in complex virtual building environments[END_REF][START_REF] Teller | Partitioning and ordering large radiosity computations[END_REF][START_REF] Meneveaux | A new partitioning method for architectural environments[END_REF][START_REF] Fradin | Out-ofcore photon-mapping for large buldings[END_REF], the efficiency of cells-and-portals have been demonstrated using only neighbourhood information between volumes.

In this article, we propose a new acceleration structure for ray tracing in architectural environments. This structure consists in a topological model corresponding to a 3 dimensional space partition. It benefits from all topological information available, i.e. neighbourhood relations between volumes (as in previous works), and also between faces, edges and vertices. This new structure is well-suited for ray casting purpose, ray tracing with reflection, refraction and direct or global illumination. It comes with an efficient beam traversal algorithm, that uses all the topological relations to schedule and speed up ray intersection tests. Using this topological information, a simple calculation associates each cell with a list of potentially visible light sources. Thus, our structure traversal leads to interactive frame rates using ray tracing, even for scenes containing many light sources (cf. Figure 1). This paper is organised as follows: Section 2 reviews previous works on acceleration structure for ray tracing in general and then specifically focuses on those dedicated to architectural environments. Section 3 introduces our topological acceleration structure. Section 4 presents our topological strategy used to accelerate the traversal of one ray through the scene. Section 5 extends the traversal of our structure for reflection and lighting, as examples of its capabilities. Section 6 gives and analyses some ray-tracing results. At last, Section 7 concludes and gives some future work lines.

RELATED WORKS

Since 1968 and Appel's works [START_REF] Appel | Some techniques for shading machine renderings of solids[END_REF], ray tracing has been used to compute high quality realistic images [START_REF] Whitted | An improved illumination model for shaded display[END_REF][START_REF] Cook | Distributed ray tracing[END_REF][START_REF] Glassner | An introduction to ray tracing[END_REF] including reflection and refraction phenomena or global illumination account. While the first interactive ray tracer was proposed by [START_REF] Parker | Interactive ray tracing[END_REF] using a large shared memory supercomputer, interactive frame rates rendering using ray tracing still remains a non trivial task. Most of the time, an acceleration structure is built over the scene in order to improve ray tracing performances. In this section, we briefly present the common acceleration structures. Then, we focus on structures dedicated to the architectural environments.

General Acceleration Structures

A naive approach of tracing a ray through a scene would consist in testing the intersection between rays and all the polygons making up the scene. For interactive rendering, such an approach is obviously prohibited. An acceleration structure is designed to decrease this linear complexity. It relies on a special organisation of the geometry which aims to speed up a random access to a polygon to a logarithmic time. As a consequence, the number of ray intersection tests performed is highly reduced and the rendering is faster. Many acceleration structures have been proposed in the literature such as kd-trees [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF], bounding volume hierarchies (BVH) [START_REF] Rubin | A 3-dimensional representation for fast rendering of complex scenes[END_REF][START_REF] Kay | Ray tracing complex scenes[END_REF] or regular grids [START_REF] Fujimoto | Arts: accelerated ray-tracing system[END_REF]. A survey comparing all these ac-celeration structures on CPU can be found in [START_REF] Havran | Heuristic Ray Shooting Algorithms[END_REF].

These structures are general, in that they work with any kind of scenes. Nevertheless, in some cases their traversal cost can become very high. That is the case for architectural environments. First, buildings can contain little but highly detailed objects inside big empty spaces. This kind of configuration is well known as the "teapot in the stadium" problem. It slows down traversal performances of regular structures like regular grids. A second problem appears with architectural environments with walls in general configurations. With axis-aligned partitioning methods, such as kd-tree or BVH, the cutting splits do not respect the wall positions and produce a complex subdivision scheme. These two problems explain for the most part why general acceleration structures are not favoured to render architectural scenes.

Architectural Acceleration Structures

Architectural scenes are made up of large occlusive surfaces such as floors, ceilings or walls. From a given point of view only a small part of the scene is visible. The structure called cells-and-portals takes advantage of this particular structural organisation. In such a structure, a cell represents a volume in which light can be propagated, such as a room or an opening. A portal corresponds to a face incident to two cells (cf. Figure 2). A cells-and-portals is stored in a graph, where a node represents a cell and an edge corresponds to a portal. The principle of traversal consists in scanning such a graph. From a given cell C, a ray r can exist only by a face incident to C. If this outgoing face is a portal, then r is propagated in the neighbouring cell; else it is stopped (cf. Figure 2(d)). Due to the local complexity of this process, it remains efficient even for large architectural scenes.

Room

Several methods have been proposed to construct such a structure starting from a list of polygons. Airey et al. [START_REF] Airey | Towards image realism with interactive update rates in complex virtual building environments[END_REF] described a method for architectural scenes with axis-aligned walls. The subdivision is performed using a kd-tree. It has been extended for general convex cells (non necessarily axis-aligned) using a BSP-tree by Teller et al. [START_REF] Teller | Visibility preprocessing for interactive walkthroughs[END_REF]. Meneveaux et al. [START_REF] Meneveaux | A new partitioning method for architectural environments[END_REF]) used a rule-based system which consists in finding the openings (i.e. the portals), in order to use them to construct cells that fit better with walls. All these works aim at extracting neighbourhood information between volumes. The problem of these methods is that the cells are not optimally determined so that a room could be uselessly split several times.

In [START_REF] Fradin | Out-ofcore photon-mapping for large buldings[END_REF], cells-and-portals are directly extracted from architectural environments which are built by using a dedicated topological modeller. Resulting scenes are adapted for lighting simulation in large buildings as shown in [START_REF] Fradin | Out-ofcore photon-mapping for large buldings[END_REF][START_REF] Fradin | A hierarchical topology-based model for handling complex indoor scenes[END_REF]. Neighbourhood information between volumes are directly known from the topology. This work produces an optimal cells-andportals structure in which a room is represented by an unique cell. Nevertheless, it has two main drawbacks. Firstly, the full potential of the topological model is not exploited: Indeed, the neighbourhood relations between faces, edges and vertices are known but unused. Secondly, the construction is entirely manual so that designing a scene is a long and tedious task. In addition, scenes can not be validated when no architectural constraints or editing operations are provided.

In this paper, we use the model introduced by Horna et al. in [START_REF] Horna | Consistency constraints and 3d building reconstruction[END_REF]. It is also a topological model, but it can produce automatically the piece of information needed to cell-and-portals traversal. This model is described in more details in Section 3. The associated modeller allows to automatically generate the 3 dimensional structure from an architectural plan [START_REF] Horna | Building 3d indoor scenes topology from 2d architectural plans[END_REF]. Moreover, contrary to Fradin et al. we propose to exploit all its topological properties, leading to a very efficient traversal.

TOPOLOGICAL STRUCTURE

In this section, we present our topological acceleration structure for architectural environments.

Generalized Maps

An architectural environment can be considered as a set of volumes. Each volume corresponds to a unique element in the scene (room, wall, door, window, etc.). Therefore, an adjacency graph is sufficient to store such a structure.

Our acceleration structure is based on the topological model introduced by Horna et al. [START_REF] Horna | Consistency constraints and 3d building reconstruction[END_REF]. It uses the notion of topological dimension, which is associated with geometrical objects: Dimension 3 with a volume, dimension 2 with a face, dimension 1 with an edge and dimension 0 with a vertex. This model contains topological relations between all these geometric elements, from dimension 0 to 3. It is based on generalized maps, or G-maps [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasi-manifolds[END_REF], extended by the formal definition of a set of topological, geometrical and semantic constraints specifically dedicated to architectural modelling. These constraints ensure the model validity during all the modelling process. In an architectural environment, all elements are real world objects. In such a scene, two elements can not occupy the same space, and the whole space is filled, authorizing a space partition. Consequently, topologically speaking, a building should be defined as a subdivision of space into volumes, faces, edges and vertices. In this article, we only recall the main principles of generalized maps applied to buildings, i.e. the neighbourhood information. All the definitions and properties of n-dimensional generalized maps can be found in [START_REF] Lienhardt | N-dimensional generalized combinatorial maps and cellular quasi-manifolds[END_REF]. The Figure 3 shows how topological information is organised (for simplification, links of dimension 0, binding two vertices to make an edge, are not represented).

E 2 E 1 1D link Face (a) F 2 F 1 2D link (b) V 1 V 2 3D link (c)

Architectural Topological Model

Architectural scenes are represented by a 3 dimensional oriented and closed space subdivision composed of elements with a significant thickness. Each From this model, we design a cells-and-portals structure in the following way: In the cells-andportals philosophy, a cell is a volume where light can be propagated. This corresponds to a given subset of semantic: ROOM, DOOR, EXTERIOR. Then, with our architectural model, the cells can be deduced automatically from the volumes.

In the same manner, a portal is a face which can be crossed by the light. With our model, portals are automatically deduced from semantics. Indeed, the only way to stop the light is to find a non cell volume. Then, using the semantic information stored into our structure, a face is a portal if and only if it links two cells (cf. Figure 4(b)).

In fact, our model is also used to find automatically the semantics of each face, even the non portal one. This allows to deduce the reflecting properties of non portal faces, and is exploited during the rendering. For example, a face incident to a ROOM and a WALL is identified as a WALLPAPER, while one inci-dent to a DOOR and a WALL becomes a GATE-POST.

Our Acceleration Structure

Architectural environments are mostly composed of rooms made with large vertical planes for the walls, and horizontal planes for floor and ceiling. Hence, we represent each volume by two horizontal faces (upper and lower ones) and a set of vertical bounding faces.

Each vertical face can be divided into several vertical parts, eg. to represent a door or a window (cf. Figure 5).

Our structure results in a space partition, representing a cells-and-portals structure in which a cell corresponds to an unique room. Nevertheless, for rendering efficiency our structure traversal is optimised for convex cells. Then, each concave cell is automatically subdivided into a set of convex cells (sharing portals). Thus, finally, a room is made of several cells.

Our acceleration structure is stored into a compact graph structure. A node corresponds to a cell and the graph links are the topological links of dimension 3 (the neighbourhood between cells i.e. the portals). Each node (so cell) contains its 2 dimensional lower face made of vertices. It is ordered counter clockwise (i.e. from right to left from inside the volume). Each vertex contains a list of horizontal edges representing the vertical faces, with their semantics. It is ordered from the bottom up. At least, a node contains a ceiling height. This particular data structure (organised according to topological elements) allows to schedule and then to speed up the traversal process. Next section describes the traversal algorithm for one ray.

STRUCTURE TRAVERSAL

The traversal of our structure requires displacement in a graph. In fact, it consists in a classical cells-

GRAPP2014-InternationalConferenceonComputerGraphicsTheoryandApplications

and-portals traversal but improved by using available topological links in all dimensions.

In this section, we firstly present the principle of a basic ray casting in our topological structure. Then, we recall the Plücker coordinates, which are used for intersection tests. Finally, we explain our topological strategy use to improve the outgoing face search.

Basic Ray Casting

Without considering the optimisation brought by our topological strategy for the outgoing face search, our acceleration structure is a kind of cells-and-portals. Then, a basic ray casting is a straightforward problem. All the rays are propagated through the scene starting from the volume containing the camera. Nevertheless, this algorithm hides two key points.

The first on resides in the determination of the initial volume (or camera volume). This volume is calculated once and only once. Then, when the camera is moved by the user, the camera volume is updated according to the translation vector defining the movement; this costs nothing since it simply consists in casting a finite and short ray. Thus, the updated camera volume corresponds to the ray arrival volume.

V 1 : ROOM 1 r 1 V 2 : ROOM 2 V 3 : ROOM 3 V 4 : EXTERIOR WALLPAPER r 2 r 3 r 4 r 5 GROUND PORTAL FACADE
Figure 6: Scene traversal: Rays start from volume V 1 and are propagated when they hit a face PORTAL, or are stopped when they hit a boundary face (WALLPAPER, GROUND, etc.). As an example, ray r 1 hits a face in V 1 with semantics WALLPAPER, while ray r 4 goes through three portals before being propagated inside volume V 4 of semantics EX-TERIOR.

The second key point lies in finding the outgoing face of a given ray in a given volume. If this face has a semantics PORTAL, then the ray is propagated through its neighbouring volume by using topological relation of dimension 3. Else, the hit face is a boundary face (semantics WALLPAPER, GROUND, etc.) and the ray is stopped (cf. Figure 6).

As mentioned in Section 2.2, the global number of faces making up the scene does not really matter due to the local complexity of the traversal. Given that the propagation of a ray from a volume to its neighbour is constant by using topological relation of dimension 3 (i.e. volume neighbouring), the efficiency of such an algorithm lies in finding quickly the outgoing face of a given ray inside a volume. To improve that and thus the global traversal, we propose an outgoing face search scheduling based on the topology of the scene. While previous works only used neighbouring relation between volumes (of dimension 3), we also use relation between faces (dimension 2), edges (dimension 1) and vertexes (dimension 0), in order to speed up the outgoing face search. Section 4.3 explains in details this process.

Pl ücker Coordinates

Plücker coordinates was firstly introduced in [START_REF] Plücker | On a new geometry of space[END_REF]. Generally speaking, they are used to represent linear subspaces of dimension k in a projective space of dimension n. In computer graphics, they represent oriented lines of 3 dimensional geometrical space. They are particular point in the projective 5 dimensional space P 5 . More precisely, an oriented line l defined by two distinct points A and B is represented by a sextuplet of coordinates Π l = {u = b-a : v = a×b}, where the 3 dimensional vectors a and b are respectively the position of the points A and B. Hence, the vector u corresponds to the direction of the line, while v is its position (or its mechanical moment).

As a example of Plücker coordinates usefulness, they are used in visibility problems as shown in [START_REF] Teller | Global visibility algorithms for illumination computations[END_REF][START_REF] Charneau | Exact, robust and efficient full visibility computation in the Plücker space[END_REF][START_REF] Fang | Mesh-based Monte Carlo method using fast ray-tracing in Plücker coordinates[END_REF][START_REF] Mora | Lazy visibility evaluation for exact soft shadows[END_REF]. Here, we use them to determinate the relative orientation between two oriented lines via the operator side [START_REF] Shoemake | Plücker coordinate tutorial[END_REF]. This operator is defined for two lines l and l ′ as side(Π l , Π l ′ ) = u.v ′ + u ′ .v. If the result is negative, l turns clockwise around l ′ ; if it is positive, l turns in the counterclockwise around l ′ ; else, when side is null, the two lines have a common point (in fact they are coplanar, and so they meet at least at infinity).

Outgoing Face Search

Using the topological relations available in our structure, we develop a simple and efficient algorithm for searching the outgoing face of a given ray into a volume (cf. Algorithm 1). It is divided into two main steps.

Let us recall that each volume is convex. Then a ray can only exit through an unique face. First, we consider the volume as an infinite vertical polyhedron, ignoring the upper and lower faces. Therefore, we start to search the outgoing infinite vertical section (Algorithm 1 -line 3-7).

An infinite vertical section is bounded by two vertical lines oriented from the bottom up. Using Plücker lines, it is trivial to test if a ray exits through a given vertical section, by testing its side product with the two bounding lines. Indeed, a ray intersects such a section when it goes through the left (resp. right) of the right (resp. left) line. This search is done linearly, by scanning the vertices of the lower face. It uses both topological relations of dimension 1 and 0 (i.e. links between edges and vertices). Considering the low number of vertical sections per volume (in average 4.37) a more complex scan (like dichotomy) is useless.

Algorithm 1: Outgoing face search for a given ray r into a volume V . Require: V : volume; r: ray;

1: {First step: search the outgoing vertical section} 2: i ⇐ 0; 3: while (sideVertical(r,

V .vSec[i].left) ≤ 0 � sideVertical(r, V .vSec[i].right) > 0) do 4:
i ⇐ i + 1; 5: end while 6: {Second step: search the outgoing face in the vertical section i} It should be noticed that this first step is done strictly in dimension 2 (see Figure 7). If x and y are the 2D coordinates of a line bounding the wall (with direction u = z), and for a ray defined with Plücker coordinates Π l = {π 0 , π 1 , π 2 , π 3 , π 4 , π 5 }, then the side operator is reduced to: y × π 0 -x × π 1 + π 5 .

7: h ⇐ V .vSec[i].hSec[]; 8: if sideHorizontal(r, h[0]) < 0 then
Once the outgoing infinite vertical section is found, we search for the outgoing face. We consider the outgoing vertical section as a set of infinite faces bounded by two horizontal lines (see Figure 8). Let us recall that these lines are oriented from right to left in our acceleration structure. Then we can use Plücker side product to test the ray intersection. We start from the lower edge (Algorithm 1 -line 9-12): If the result is negative, the ray turns clockwise around the edge (the ray goes below the edge); then the floor is hit and can be returned. Otherwise, we repeat this process with upper edges by scanning them (by using topological relations of dimension 2, i.e. links between faces) until finding the outgoing face (Algorithm 1line 13-17). If none face is found, then the ceiling is returned (Algorithm 1 -line 18). Once more, a nonlinear strategy is useless in front of the low number of faces (in average 2.11). Considering that bounding lines are only horizontal, Plücker side operator is reduced to four products and five additions.

To sum up, outgoing face search is optimised by using topological relations of all dimensions:

• Dimensions 0 (vertex) and 1 (edge) to scan the lower face of a volume and to find the outgoing vertical section.

• Dimension 2 (face) to change of face within a vertical section and determine the outgoing face.

• Dimension 3 (volume) to move from a volume to its neighbour when propagating a ray.

GRAPP2014-InternationalConferenceonComputerGraphicsTheoryandApplications

EXTENDED USES

This section discusses about the integration of our model into a ray tracing application to render architectural scenes, using wide ray packets and secondary rays (reflection, refraction, and light rays).

Ray Packet

The basic ray casting algorithm is modified in order to group the rays in packets. In practice, we use can very wide packet: depending on the number t h of threads, we divide the image into 2 × t h packets. In our experiments, this leads to 16 packets for 1024 × 1024 pixels.

A packet is defined by its 4 corners, corresponding to 4 pixels in the ray casting (and idem for secondary rays). From the initial volume, the outgoing face of this 4 rays is searched according to Algorithm 1. When the results are coherent (the 4 rays goes through the same face), the algorithm computes successfully the outgoing face for all internal rays.

Else, when rays are incoherent, the algorithm can not continue for the given packet. For each corner rays, all the current data is stored: i.e. the identifiers of the volume, of the vertical section and of the face of this vertical section. Then, the packet is split in 2 or 4 parts, depending on the incoherence and each part continue in the same way, from the same initial position.

In practice, this strategy is very efficient. It uses SIMD float point instruction as, for instance, in the kd-tree packet traversal proposed in [START_REF] Wald | Interactive rendering with coherent ray tracing[END_REF]. Since the new packets continue from where they are split we do not spend time in a new structure scanning.

Reflection and Refraction

Our model allows to associate any kind of behaviour to a face thanks to the face semantics. This allows to add some reflection or refraction properties to the surfaces, as an example for panes.

When a ray hit a reflective or refractive surface, its trajectory is modified. For a both reflective and refractive surface, it is also subdivided. Since our traversal algorithm have a local scope, the cost of a reflective or refractive ray stays about the same that for primary rays. For example, with a specular surface, a ray can be reflected a given number of times. For each bounce, a new ray have to be cast, with new Plücker coordinates. In fact, specularity account is not a difficult task. The new coordinates of the specular ray are computed according to a traditional scheme (cf. Section 4.2) and the ray is propagated by using the Algorithm 1, or its packet version. The main difference resides in the starting cell: The volume containing the reflection point is not the same that the camera one. Fortunately, this volume do not have to be computed: It corresponds to the arrival volume of the primary ray.

Direct Illumination

Computing the direct illumination is a good way to evaluate the robustness of our structure in case of spatially incoherent rays. That implies new rays (light rays) to be cast through the scene, to determine the point-to-point visibility between the primary ray intersections and the point light sources. Each new ray is simply cast, as a primary ray (see Section 4). Then, lighting is computed according to the visibility and the chromatic intensity of each visible source. The main difficulty for direct illumination account in architectural environments lies in the huge number of light sources which highly increase the number of rays to be cast, leading to low frame rates. Fortunately, given that buildings are composed of large occlusive surfaces, too many useless tests are performed (cf. Figure 9(a)). Cells-and-portals structures are used in [START_REF] Airey | Towards image realism with interactive update rates in complex virtual building environments[END_REF], [START_REF] Teller | Visibility preprocessing for interactive walkthroughs[END_REF], [START_REF] Luebke | Portals and mirrors: simple, fast evaluation of potentially visible sets[END_REF] or [START_REF] Meneveaux | A new partitioning method for architectural environments[END_REF] to compute the potentially visible set (or PVS) in architectural environments. From the same perspective and for efficiency purpose, we choose to compute the set of potential visible lights for each cell making up our scenes. We proceed by launching continuous beams from a light source, one through each opening of the volume containing the light source(cf. Figure 9(b)). Parts of the beam which hit a face POR-TAL are recursively propagated through neighbouring cells until being stopped by a non-portal face. As usual, the propagation is performed by using topological relation of dimension 3. Each cell crossed by the beam is potentially illuminated by the light source. Instead of having to cast a light ray for each light source, computing PVS allows to highly reduce the number of lights to a local list of potentially visible sources. As an example, in one of our test scene which contains 199 light sources, the average number of light rays per image is reduced from 208.67 millions to 10.14 millions (-95.14%).

RESULTS

We evaluate the performances of our acceleration structure on a Intel R � Core TM i7 CPU 960 @ 3.20GHz with 12GB of RAM. The traversal algorithm runs in parallel using OpenMP. The packet version also uses 4-wide SIMD floating point instructions (Intel R � SSE). Performances are given in terms of:

• FPS: frames calculated per second.

• Mrays/s: millions of rays cast per second. The rendering is done for image of 1024 × 1024 pixels with one ray per pixel (without anti-aliasing). In order to have relevant statistics, we use a way-point system to keep the same points of view for each algorithm tested. These points of view result from a walk-through into each scene; then the number of points of view can differ according to the size of the scene.

Notice that the rendering time include the shading, using a Lambertian brdf. With ray-casting, a virtual point light is positioned at the camera location, and used for the calculation of the cosine of the incident direction with respect to the surface normal.

Scenes

We present four of our test scenes which are more or less complex in terms of size, number of polygons and light sources. Their main characteristics are summarised into Table 1. The smallest scene corresponds to a house made up of 8 rooms (cf. Figure 12). It has 1.166 polygons and contains 8 light sources. The scenes BUILD1 and BUILD2 represent two distinct floors of a large administration building (cf. Figures 13,1 and 11). They are respectively made up of 46 and 135 rooms, and contains 8.557 and 29.937 polygons with 55 and 199 light sources. The last scene, BUILD3, is a very large building made up of 13 floors with 1755 rooms, with 369 thousands of polygons and 2557 point light sources (cf. Figure 14). These four scenes result directly from our topological modeller.

GRAPP2014-InternationalConferenceonComputerGraphicsTheoryandApplications

Results and Discussion

Primary Ray Casting

We propose two different algorithms for primary ray casting. While the first one works pixel per pixel, the second corresponds to the packet traversal described in Section 5.1. These two algorithms have been evaluated with each of our test scenes. We can notice that primary ray casting using our topological structure is very efficient (see Table 2). The Figure 11(a) shows the rendering for the second building test scene, using a simple ray casting. The Figure 11(b) proposes the same point of view, but using and emphasizing ray packets. Clearly, the packet are split only near cell boundaries. For large wall, ground or ceiling parts, a wide ray packet can be efficiently traced with our acceleration structure.

Table 3 compares the performances of these two methods with those measured by using the ray tracing kernel Embree [START_REF] Ernst | Embree: Photo-realistic ray tracing kernels[END_REF]. We can observe both of them are faster than Embree. Our ray casting method is more than 4 times faster and the packet version is more than 30 times faster. Since our algorithm have a local complexity, the size of the scene does not matter since the data could fit easily on the available memory. This complexity takes into account the number of volumes crossed by the ray and the number of polygons making up these volumes (which impact on the outgoing face search). Indeed, performances are reduced in case of a cell split a lot of times or in case of volumes containing too many polygons (i.e. not enough split). Thus, the scene must be convexified using a well suited strategy.

Analogous to the well-known surface area heuristic (SAH) [START_REF] Goldsmith | Automatic creation of object hierarchies for ray tracing[END_REF][START_REF] Macdonald | Heuristics for ray tracing using space subdivision[END_REF], a concave cell have to be split considering its area and the expected cost of tracing an infinite random ray through that cell. In fact, we search for a split scheme which have the lowest SAH cost possible. The main problem with convexification occurs in case of curved walls as illustrated on Figure 10. The Figure 10(a) shows a room convexified with the method presented in [START_REF] Fernández | A practical algorithm for decomposing polygonal domains into convex polygons by diagonals[END_REF]. Because of curved walls, many splits are performed, creating a dense partition. Then, rendering an image requires to cross a lot of volumes so that performances are highly reduced. Our method is illustrated in Figure 10(b): Convexification is extended by constructing a bounding box around curved walls to restraint all the splits in a small space, leading to a better SAH cost. The use of such a strategy offers a good compromise and leads to better performances. Notice that convexification process is automatically performed with our modeller and takes from seconds to few minutes according to the size of the scene.

Direct Lighting

Taking into account direct lighting means casting a light ray for each light source. For one light source this should twice the number of rays. Then, for x light sources this can dramatically reduce the number of frames per second. As explained in Section 5.3, we compute the PVS for each cell in our structure to reduce the number of light rays.

Table 2 shows that direct lighting makes the number of rays per second drop linearly. Nevertheless, this reduction is not linear with respect to the number of light sources. Some rendering examples are proposed in the Figures 1,12, 13 and 14. This fall in terms of FPS is brought by two causes. First, new Plücker coordinates have to be computed for each light ray to be cast. Secondly, it is due to spatial ray incoherence which leads to incoherent memory accesses. The robustness of our topological structure in terms of number of rays cast per second allows to think that it could be used efficiently to compute global illumination in architectural environments.

CONCLUSIONS

This article proposes a new acceleration structure for ray tracing dedicated to architectural environments. It is generated automatically from a topological modeller, without human interaction. It takes advantage of the architectural building constraint: it corresponds to a 3 dimensional oriented and closed space subdivision composed of elements with a significant thickness.

Our structure is a kind of advanced cells-andportals, with two main improvements: Firstly, thanks to our topological model, cells-and-portals cutting planes fit perfectly to the building topology. Secondly, while previous cells-and-portals traversal algorithms used only neighbouring relations between volumes, we propose to take advantage of all the topological information available, to schedule the outgoing face search, and to accelerate the whole process.

The traversal of our acceleration structure is implemented on CPU, both with single rays or with ray packets. Our experimentation shows that our structure traversal is very efficient. Our rendering tool computes hundreds of millions of rays per second, using only a single 4-core processor. Consequently, it allows to render images interactively, for ray casting or ray tracing, with direct lighting or not.

We think that our structure and its traversal algorithm can run onto a GPU, taking advantage of the high number of cores. It should be interesting to study the results of such a structure on GPU, especially when a large number of light sources are used.

We also think that our structure can be used to include piece of furniture quite easily. Indeed, using bounding boxes surrounding high detailed objects, we should be able to include them into the topological model. During the rendering step, the rendering could be done by combining our structure with a more classical one. Thus, our topological structure traversal performances would highly depend on the acceleration structure used to render the object.

Then, a logical next step will be to adapt our topological partition to general objects in order to include it directly in our architectural structure.

The main disadvantage of our model relies on the topological modeller. It could not be generated with only a geometric description of scene, without any topological information. As a future work, we plan to adapt our ideas to general objects described only by their geometry. We aim to use generalized maps to connect geometric elements themselves, in order to generate automatically a topological acceleration structure from scratch, for architectural scenes or not. This step would allow to directly include furniture in our architectural structure and thus to avoid the dependence on an auxiliary acceleration structure.

Figure 1 :

 1 Figure 1: Examples of rendering in a building made up of 30k polygons, for 1024 × 1024 pixels without anti-aliasing. Left corresponds to a packet ray casting, right illustrates direct lighting from the 199 ponctual light sources in that scene.

Figure 2 :

 2 Figure 2: 2D cells and portals representation and usage: (a) building representation; (b) rooms and doors correspond to cells; (c) portals are shared by 2 cells; (d) rays hitting a portal are propagated through neighbouring cells.

Figure 3 :

 3 Figure 3: Topological links in each dimension: (a) links of dimension 1 bind edges to make faces; (b) links of dimension 2 tie together faces to generate volumes; (c) volumes are tied with links of dimension 3.

Figure 4 :

 4 Figure 4: Topology and semantic: (a) 2 dimensional scene complying the topological and semantic properties; each face (in 3D, each volume) corresponds to a unique element (room, wall, etc.); topology represents neighbouring relations. (b) An edge (in 3D, a face) incident to two cells becomes a PORTAL while a non-portal edge has a specific semantics (WALLPAPER, FACADE, GATE-POST).volume is identified by a unique semantic according to its nature: ROOM, DOOR, WALL, GROUND, CEIL-ING or EXTERIOR. The Figure4(a) illustrates a scene complying to this property. Note that all descriptions made in dimension 2 are extensible to dimension 3, since G-maps are homogeneous in all dimensions.From this model, we design a cells-and-portals structure in the following way: In the cells-andportals philosophy, a cell is a volume where light can be propagated. This corresponds to a given subset of semantic: ROOM, DOOR, EXTERIOR. Then, with our architectural model, the cells can be deduced automatically from the volumes.In the same manner, a portal is a face which can be crossed by the light. With our model, portals are automatically deduced from semantics. Indeed, the only way to stop the light is to find a non cell volume. Then, using the semantic information stored into our structure, a face is a portal if and only if it links two cells (cf. Figure4(b)).In fact, our model is also used to find automatically the semantics of each face, even the non portal one. This allows to deduce the reflecting properties of non portal faces, and is exploited during the rendering. For example, a face incident to a ROOM and a WALL is identified as a WALLPAPER, while one inci-

Figure 5 :

 5 Figure 5: Volume characteristics in building: They are composed of 2 horizontal faces and several vertical faces corresponding to volume border.

  ; {go through the floor} 10: end if 11: for k = 1 to V .walls[i].nb hSec-1 do 12:if sideHorizontal(r, e[k]) for 16: {Not a regular face: outgoing face is the ceiling!} 17: return V .ceiling;

Figure 7 :Figure 8 :

 78 Figure 7: First step of outgoing face search: (a) the cell is considered as an infinite vertical polyhedron; a vertical section is bounded by two lines oriented from the bottom up; outgoing vertical section is 6; (b) the process is strictly done in 2D.E

Figure 9 :

 9 Figure 9: Direct illumination: (a) for the ray r, 11 light rays are launched while only 4 are visible (in blue); (b) cells crossed by light beams are potentially illuminated, allowing to decrease the number of light rays.

  (a) Convexification using[START_REF] Fernández | A practical algorithm for decomposing polygonal domains into convex polygons by diagonals[END_REF]'s method; Rays cross 9.139.670 volumes in all to render the image (c). (b) Convexification extended with bounding boxes and topology; Rays cross 3.422.570 volumes in all to render the image (c). (c) Associated rendered image; 314,31 FPS with convexification (a) and 377,28 FPS with (b).

Figure 10 :

 10 Figure 10: Comparison of the convexification of a room.

  (a) Ray casting (48 FPS). (b) Ray packets (356 FPS).

Figure 11 :

 11 Figure 11: Illustration of primary ray casing with scene BUILD2.

Figure 12 :

 12 Figure 12: HOUSE without (a) or with (b) lighting from 8 light sources.

Figure 13 :

 13 Figure 13: BUILD1 without (a) or with (b) direct lighting. The scene contains 55 light sources.

  (a) 155 FPS. (b) 0.24 FPS.

Figure 14 :

 14 Figure 14: BUILD3 without (a) or with (b) direct lighting from 2557 lights. Most of the lights are visible outside, leading to a small FPS.

Table 1 :

 1 Characteristics of our test scenes: number of rooms, cells, polygons, lights and memory size in megabytes.

	Scenes	# rooms # cells # polygons # lights memory (Mb)
	HOUSE BUILD1 BUILD2 BUILD3	8 46 135 1755	509 2.333 8.258 95.932	1.166 8.557 29.937 369.521	8 55 199 2557	0.46 2.03 6.82 80.13

Table 2 :

 2 Results: for each scene we give the average number of frame calculated per second (FPS), and the average number of rays cast per second in millions (Mrays/s).

	Scenes Results	HOUSE FPS Mrays/s	BUILD1 FPS Mrays/s	BUILD2 FPS Mrays/s	BUILD3 FPS Mrays/s
	Ray casting Packet casting Packet Lighting	65.09 388.16 128.71	68.26 407.01 419.75	51.28 384.13 66.64	53.78 370.68 367.05	51.35 375.13 60.48	53.85 393.35 352.09	46.86 341.66 34.49	49.14 358.25 246.91

Table 3 :

 3 Comparison with Embree ray tracing software[START_REF] Ernst | Embree: Photo-realistic ray tracing kernels[END_REF].

	Scenes	Embree Ray casting Packet casting
	HOUSE BUILD1 BUILD2 BUILD3	18 12 10 8	65.09 51.28 51.35 46.86	388.16 384.13 375.13 358.25

GRAPP2014-InternationalConferenceonComputerGraphicsTheoryandApplications

TopologicalSpacePartitionforFastRayTracinginArchitecturalModels

ACKNOWLEDGEMENTS

Authors thanks the Région Poitou-Charentes for their funding support.


Sonia Cafieri
email: sonia.cafieri@enac.fr

Leo Liberti
email: liberti@lix.polytechnique.fr

F Rédéric Messine

Bertrand Nogarede
email: nogarede@n7.fr

Discussion about Formulations and Resolution Techniques of Electrical Machine Design Problems

Keywords: analytical model, formulation, local optimization, inverse problem, design, electrical machine, mixed integer nonlinear programming

In this paper, we discuss, through a simple example, the impact of two different but equivalent formulations used by standard local optimization solvers (here fmincon of MatLab). We show that even if the two formulations are equivalent in a mathematical sense (no loss of global optima) it is not completely true in a numerical way; using 1000 starting points, we show that it is quite difficult for the designer to find a starting point yielding a convergence of the algorithm to a local minimum and to find better points yielding the global solution (previously found using a global optimization algorithm). Furthermore, we discuss how to deal with the insertion of the integer variable p, representing the number of pole pairs of the machine, inside the problem of design which uses a standard local continuous optimization code to be solved.

I. INTRODUCTION

The interest of the electromagnetical actuators design combining optimization algorithms and analytical models has already been shown in a lot of works, [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF]- [START_REF] Messine | Deterministic Global Optimization using Interval Constraint Propagation Techniques[END_REF], [START_REF] Nogarede | Machines Tournantes: principes et constitutions[END_REF], [START_REF] Nogarede | Optimal Design of Permanent-Magnet Machines using Analytical Field Modeling[END_REF], [START_REF] Slemon | Modeling and design optimization of permanent magnet motors[END_REF], [START_REF] Wurtz | A Methodology and a Tool for the Computer-Aided Design with Constraints of Electrical Devices[END_REF]. The design of electromechanical actuators is understood as an inverse problem, i.e. from the characteristic values given by the schedule of conditions (for example the torque), obtain the structure, the dimensions and the material compositions of the actuator constitutive parts, [START_REF] Fitan | A General Analytical Model of Electrical Permanent Magnet Machine dedicated to Optimal Design[END_REF], [START_REF] Fitan | The Electromagnetical Actuators Design Problem: a General and Rational Approach[END_REF].

In [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF], we have shown that even for optimal dimensioning problems which can be formulated as continuous-constrained optimization problems, it was really difficult to find the global optima using classical optimization algorithms, like Lagrangian Augmented Techniques [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF], [START_REF] Wurtz | A Methodology and a Tool for the Computer-Aided Design with Constraints of Electrical Devices[END_REF]. Moreover, we have shown that exact global optimization based on interval analysis can be perfectly adapted to solve these kinds of problem [START_REF] Fitan | The Electromagnetical Actuators Design Problem: a General and Rational Approach[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems[END_REF]- [START_REF] Messine | Optimal design of electromechanical actuators: a new method based on global optimization[END_REF].

Unfortunately, the exact global algorithm named IBBA (for Interval Branch and Bound Algorithm [START_REF] Hansen | Global Optimization using Interval Analysis[END_REF], [START_REF] Keafott | Rigorous Global Search: Continuous Problems[END_REF], [START_REF] Messine | A Deterministic Global Optimization Algorithm for Design Problems[END_REF], [START_REF] Ratschek | New Computer Methods for Global Optimization[END_REF]), is difficult to implement on a computer (about 10000 lines of a Fortran90 code), and nowadays, it does not exist an industrial version of IBBA or of a similar one. Furthermore, its application to solve such problems is not so easy and it often needs to develop or to adapt some parts of the code.

Concerning stochastic algorithms, such as genetic algorithms or simulated annealing, they are difficult to be efficiently used to solve this kind of problems due to the equality and inequality constraints which define many small realizable domains.

In fact, the inverse problems of the electromechanical actuators design are more general than optimal dimensioning problems, see [START_REF] Fitan | The Electromagnetical Actuators Design Problem: a General and Rational Approach[END_REF]. Thus, these problems must be formulated as mixed-constrained optimization problems, see [START_REF] Fitan | A General Analytical Model of Electrical Permanent Magnet Machine dedicated to Optimal Design[END_REF], [START_REF] Fitan | The Electromagnetical Actuators Design Problem: a General and Rational Approach[END_REF] for more complete formulations.

In this paper, we discuss about the impact of different formulations on the efficient use of a deterministic local optimization solver (fmincon of MatLab), through a simple example of a slotless electrical rotating machine with magnetic effects. This example was first presented in [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF] and also studied in a lot of papers such as [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF], [START_REF] Messine | Deterministic Global Optimization using Interval Constraint Propagation Techniques[END_REF]- [START_REF] Nogarede | Optimal Design of Permanent-Magnet Machines using Analytical Field Modeling[END_REF], [START_REF] Wurtz | A Methodology and a Tool for the Computer-Aided Design with Constraints of Electrical Devices[END_REF].

The design problem can be formulated as follows:

                     min πβl a D λ (D -2e -l a ) s.t. Γ em = π 2λ (1 -K f ) √ k r βE ch ED 2 (D + E)B e E ch = AJ cu = k r EJ 2 cu K f =1.5pβ e+E D B e = 2laP D ln( D+2E D-2(la +e) ) C = πβBe 4pBiron D p = πD ∆p (1)
where D(m) is the bore diameter, λ the diameter over length ratio, l a (m) the thickness of the permanent magnets, E(m) the winding thickness, C(m) the thickness of yoke, β the polar arc factor, B e (T ) the magnetic field in the air-gap, J cu (A/m 2 ) the current areal density, K f a semi-empiric magnetic leakage coefficient (established by numerical simulations), e(m) the thickness of the mechanical air gap, p the number of pole pairs, k r a coefficient of occupation, B iron the magnetic field in the iron, Γ em the electromagnetical torque, P the magnetic polarization and ∆ p the polar step. For this study, we fix Γ em =1 0 N.m, P =0 .9T , k r =0 .7, B iron =1 .5T , E ch =1 0 11 A/m and ∆ p =0 .1m the polar step as in [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF]. The other parameters can vary inside the following intervals: This problem is mixed (the parameter p is an integer number) and also non-homogeneous (some parameters, such as J cu , have high values and the others are quite small, such as D, e, ...). The designer must take care about this point.

D(m) ∈ [0.01, 0.5], λ ∈ [1, 2.5], l a (m) ∈ [0.003, 0.05], E(m) ∈ [0.001, 0.05], C(m) ∈ [0.001, 0.05], β ∈ [0.8, 1], B e (T ) ∈ [0.1, 1], J cu (A/m 2 ) ∈ [10 5 , 10 7 ], K f ∈ [0.01, 0.
In Section II, we present two different but equivalent formulations of the design problem and we discuss about their efficiency when the optimization problem is solved using the function 'fmincon' of MatLab 7 . In Section III, we discuss about the fact that the parameter p representing the number of pole pairs can become an integer variable of the optimization code. This yields a mixed-integer non-linear optimization problem which is difficult to solve using classical continuous local search algorithms such as fmincon of MatLab. Some concluding remarks are given in Section IV.

II. FORMULATIONS AND STARTING POINTS

One of the main difficulty in the use of a deterministic local optimization algorithm is the introduction of a starting point to initialize it. Then, we have three possibilities: (i) the algorithm converges to a local solution; (ii) the algorithm works very slowly and, after a while, it reaches the maximum number of iterations or of function evaluations (here, both are fixed to 30000); (iii) the algorithm does not converge and fails.

If we are on stage (ii), we can restart the algorithm by increasing the maximum number of iterations and of function evaluations but it could not work (case (iii) or still case (ii)).

In this section, p is fixed to 5 and we discuss about the two following formulations:

               Formulation 1:
-D, λ, l a ,E,C,β,J cu ,K f ,e: 9 variables.

-B e is a function depending on D, λ, ... -Equations in (1) yield 5 equalities constraints and 2 inequality ones: 0.1 ≤ B e (D, λ, ...) ≤ 1.

-Add constraints to the bounds of the variables.

(2)

       Formulation 2:
-D, λ, l a ,E,C,β,J cu ,K f ,e,B e : 10 variables.

-Equations in (1) yield 6 equalities constraints.

-Add constraints to the bounds of the variables.

(

These two explicit formulations are:

• Formulation 1:

                             min (D,λ,••• )∈D⊂I R 9 V m (D, λ, ...)=πβl a D λ (D -2e -l a ) Γ em = π 2λ (1 -K f ) √ k r βE ch ED 2 × (D + E)B e (D, λ, ...) E ch = AJ cu = k r EJ 2 cu K f =1.5pβ e+E D C = πβBe 4pBiron D p = πD ∆p Be(D, λ, ...) ≥ 0.1 Be(D, λ, ...) ≤ 1 (4)
Where B e (D, λ, ...) is defined as an auxiliary function and returns

2laP D ln( D+2E D-2(la +e) ) . • Formulation 2:                            min (D,λ,••• )∈D⊂I R 10 V m (D, λ, ...)=πβl a D λ (D -2e -l a ) Γ em = π 2λ (1 -K f ) √ k r βE ch ED 2 × (D + E)B e E ch = AJ cu = k r EJ 2 cu K f =1.5pβ e+E D B e = 2laP D ln( D+2E D-2(la +e) ) C = πβBe 4pBiron D p = πD ∆p (5)
D define the initial domain where the solution is searched for; it defines lower and upper bounds for all the variables.

Mathematically this two formulations are equivalent because the only difference is the addition of a new variable B e ∈ [0.1, 1] in Formulation 2 and the transformation of 2 inequality constraints into an equality one. Therefore, this reformulation involves nothing in the mathematical point of view: any solution (global minimum) of Formulation 1 is also a solution for Formulation 2 and conversely.

Table I presents numerical results obtained for the two described formulations, using the same starting points. These starting points are randomly generated between the bounds of the variables by using the function 'rand' of MatLab (for Formulation 2, one random number is necessary for B e ). Some attentions must be paid to discard the fact that using some random numbers, some evaluations of negative logarithms could be performed. We performed 1000 iterations of 'fmincon' for each formulation using 1000 randomly generated starting points. In the first line of Table I we give the percentage of starting points yielding a local minimum (the percentage of failure is quite high). In the second line, we show the percentage yielding the best solution found during the 1000 iterations (comparison of the values at 10 -7 ). This percentage is very low, confirming that a lot of care must be paid in the choice of a good starting point. In the following of the table, the best and worst values of local minima and the best and worst CPU-times used for one iteration are reported. The last line indicates that taking the starting point in the middle of the hypercube defined by the bounds of the variables, Formulation 1 provides a bad solution and Formulation 2 does not converge; this proves the non-equivalence of the two formulations in a numerical sense. However, we note that the two formulations provide similar results. The percentage of time when the two formulations provide equivalent local values (within a tolerance of 10 -7 ) at the same iteration is about 3.3%. We remark that Formulation 2 is more efficient providing best convergence to the best local minimum (which is the global one); it seems to be more efficient to use this formulation with the introduction of a new variable B e because the expressions of the equations to be satisfied become less complicated.

The best found local solutions for Formulation 1 and Formulation 2 are the following: Thus, the two formulations provides two similar best local solutions.

Moreover, we note that the number of pole pairs fixed to p =4yields a best local minimum equal to 8.4272e -5 for Formulation 1 and 8.4097e -5 for Formulation 2, which are about 15% larger than the value corresponding to the case p =5.

In the next section, we discuss how to deal with p as an integer variable of the optimization problem yielding a difficult mixed integer non-linear program.

III. INTEGER PARAMETER

In such a dimensioning problem, the number of pole pairs p was previously considered as a fixed parameter, p =5. In fact, Quasi-Newton algorithms such as fmincon of MatLab needs continuous formulations and twice differentiable functions. As a consequence, integer parameters must be dealt with attention.

In the literature, the integer parameter p is usually considered as a continuous one. The user has to convert the obtained real value of p into an integer one, [START_REF] Kone | Le Dimensionnement des Actionneurs Electriques: un Problème de Programmation Non-Linéaire[END_REF], [START_REF] Nogarede | Optimal Design of Permanent-Magnet Machines using Analytical Field Modeling[END_REF]. Here, by considering p ∈ [1, ••• , 10], we performed a multistart method with 1000 starting points (randomly generated) and the results are presented in Table II, where p * corresponds to the best found local minimum and p mid to the local minimum found using the middle of the the domain as a starting point.

The best found local minima for Formulation 1 and Formu- We can note that these two solutions are different even if their corresponding volume values are very close to each other.

From results in Table II, we remark that Formulation 2 gives again the best results in terms of the percentage of the found local minima and of the best found local minimum. The real values of the number of pole pairs p in both cases are very close and larger than 4.86, thus confirming that the integer value p =5(used for experiments discussed in Section I) is a good choice.

Note that only 31 starting points over 1000 provide close optimal solutions where the values of the volume differ at most of 10 -7 .

If we fix p * =5, which is the integer value the most close to p * , we obtain the following table containing the values of the equalities constraints. They must be close to 0 (less than 10 -8 ) to be considered as satisfied constraints:

Functions Formulation 1 Formulation 2 Va = 6.9894e -5 6.9851e -5 eq1 = 9.7964e -5 3.7897e -4 eq2 = -4.7982e -5 -8.2906e -5 eq3 = -4.3537e -3 -3.9579e -3 eq4 = 1.3943e -4 -1.2990e -5 eq5 = 1.3826e -1 1.4987e -4 eq6 = - 1.3974e -1
Of course, the equality constraints are now not satisfied and therefore one has to change some values of the solution previously found. In the following, we discuss some of these possibilities to deal with the integer parameter.

Another possibility to deal with the fact that p is integer, is to insert a new equality constraint which is satisfied when

p ∈{1, ••• , 10}: (p -1) × (p -2) ו••×(p -10) = 0.
Numerical results for both formulations are shown in Table III, where again p * corresponds to the best found local minimum. As in the results previously discussed, the two solutions are different even it close values are found for the volume of magnet.

Note that in this case only 49 starting points over 1000 provide close optimal solutions where the values of the volume differ at most of 10 -7 .

Another strategy to deal with p consists in starting from the real solutions given in Table II, and changing the last added equality constraint to a lower degree polynomial expression as follows:

(p -([p * ] -1)) × (p -[p * ]) × (p -([p * ]+1))=0,
where [p * ] represents the most close integer value of p * .I n our experiments, we used for p * the value 4.86 that was found by the two formulations for the continuous problem. Results are reported in Table IV. We note that the percentage to find a local minimum is a little lower than in the previous case, while the percentage to find the best local minimum is significantly improved. Therefore, this seems to show that it is beneficial to reduce the zone of research, concentrating the efforts around p =5 , in order to improve the chances to find the global minimum. Again the two solutions differ, dispite of equivalent values of the volume of the magnet. We remark that only 30 starting points over 1000 provide close optimal solutions, where the values of the volume differ at most of 10 -7 .

In Figures 1 and2 we compare the two formulations in terms of the percentage of found local minima and the percentage of best found local minima respectively, considering the different choices of p discussed above. So, these figures summarize the results presented in Tables I to IV. From Figure 1, one can see that there is not a clearly most efficient method to find the largest number of local solutions. This actually depends on the considered problem. We can note that considering the parameter p as a variable, the second formulation seems to be more efficient.

In Figure 2, it is most clear that, using the second formulation, a higher number of starting points allows to obtain the best found local solution. Thus, as the values of the percentage of best found local solution are low, we can have more confidence on the local solutions found by the second formulation in order to provide the global solution. We have presented two reformulations of a design problem, showing that the designer must take care about the formulation of the problem in order to make more efficient the use of a local standard algorithm such as that proposed by MatLab in fmincon.

First, we have discussed the solutions found using the two formulations, where the number of pole pairs is fixed to 5 and 4. Formulation 2 provides the best percentage of local solutions found.

Second, we have discussed about the fact that p has to be an integer variable, providing a mixed integer non-linear problem. Using the two formulations and some different strategies to deal with the integer variable, we have remarked that Formulation 2 is again most efficient in providing the best local solution. However, the two best local solutions, found using one or the other formulation, are different but give the same value of the volume of the magnet. Furthermore, using the middle of the domain as a starting point, we found that Formulation 1 provides a solution which is in general far from the best one, while using Formulation 2 the process never converges.

  3], e(m) ∈ [0.001, 0.005] and p ∈{1, ••• , 10}. Furthermore, the function to be minimized corresponds to the magnet volume (V m = πβl a D λ (D -2e -l a )).

XIX

  International Conference on Electrical Machines -ICEM 2010, Rome 978-1-4244-4175-4/10/$25.00 ©2010 IEEE

Fig. 1 .

 1 Fig. 1. Percentages of found local minima

Fig. 2 .

 2 Fig. 2. Percentages of best found local minima

TABLE II COMPARISON

 II BETWEEN FORMULATIONS 1 AND 2 WITH p ∈[START_REF] Fitan | A General Analytical Model of Electrical Permanent Magnet Machine dedicated to Optimal Design[END_REF][START_REF] Nogarede | Machines Tournantes: principes et constitutions[END_REF]

		Formulation 1	Formulation 2
	% local min	51.7%	58.1%
	% best min	1.4%	1.8%
	best min value	6.9894e -5	6.9851e -5
	p *	4.8617	4.8603
	worst local min	1.3261e -3	1.3740e -3
	best CPU-time	0.00s	0.02s
	worst CPU-time	15.02s	20.47s
	x0=mid	8.8957e-4	2.2991e-4
	p mid	3.4728	3.7814
	lation 2 are:
	Best local min Formulation 1 Formulation 2
	D * =	0.154754	0.154707
	λ * =	2.5	2.5
	l * a =	0.003	0.003
	E * =	0.003066	0.002659
	C * =	0.00506	0.005356
	β * =	0.8	0.8
	B * e =	0.37953	0.401686
	J * cu =	6825814.56	7329490.86
	* f =	0.153290	0.137949
	e * =	0.001	0.001
	p * =	4.861748	4.8603


Marcella Autran

Burlier Drummond

Gustavo Luis
email: lgustavozelaya@gmail.com

Amaranto Lopes Cruz

Felix Pereira

Mora-Camino

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

Luis Gustavo

Zelaya Cruz

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

Felix Mora-Camino
email: felix.mora@enac.fr

A bi-level scheme for assessing the impact of air transportation on local development

Keywords: Transportation networks, airports, air transportation, regional development

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

INTRODUCTION

Transportation and economic development have always been entwined sectors. 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 air transportation and economic activity is quite complex.

With respect to regional economy, air transportation 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.

The main objective of this study is to provide a framework to assess the impact of the implantation of a new airport in a rural area where some other airports as well as ground transportation networks already provide transportation services to passengers and freight. It is supposed that each airport, pre-existing or new is attached to a particular urban center.

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.

The generation model: Here we consider that transportation modalities from a given local urban center to and from a major urban center outside the considered region 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 local center i to interregional center j through generalized transportation modality m. Figure 1 gives a sketch of the two retained level space analysis. Fig. 1 Hierarchical view of intra and inter regional transport These generalized modalities include multimodal transportation corresponding for instance to the use of concurrent transportation facilities (airports) from other local centers. The modalities available for local center i to regional center j are given by the set ij M , with:

  R r A a I i Pop O i O ar i ar i pax iar     , 0 ,   
A ij A ij ij M M M   (2)
(3) local regional centers

j 1 j 2 j 3 j 4 j 5 i 1 i 2 i 3 i 4 i 5 51st ERSA 2011

SPECIAL SESSION Air Transportation and Local Development (ZL SS)

where A ij M is the set of transportation paths using a local airport, starting from the airport of local center i when it exists, A ij M is the set of other transportation paths between local center i and interregional center j.

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 interregional center j to local center i with each modality m [ m ji ˆ], we can consider the solution of the following max entropy problem for each local center and each 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               R r A a I i D O T i pax iar pax iar arm ij M m J j ij          , , I i C T pax m i arm kj I k J j R r A a ij i           0  arm ij T where ij arm kj T
is the distribution trip from the local regional center k to the interregional center j,

for each pair (a,r) and using the same modality m from the intermediate city i to arrive at j.

pax m i C
is the infrastructure capacity of the local center i, considering the passenger modality m used.

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.

Then the total passenger demand for airport i is estimated by:

ij i arm kj I k J j R r A a i T          
In the case where:

pax i arm kj I k J j R r A a i C T ij            (6) (9) (4) (5) (8) (7)
the no-capacity saturation assumption will lead to install an added capacity given by:

pax i i i C C    

PREDICTION OF FREIGHT DEMAND DISTRIBUTION

Here we distinguish 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 reexported production is given by

i b i ba a b A b i a V V i      , ,
           i i r i ar i R a R r i a aa am ji M m J j A a Pop V i i a                Im Im 0  am ji  i i a M m A a J j    , , Im
Where am ji  is the imported goods flow of sector a with the modality m from j to the regional center i, m ji P is its generalized transportation price and An estimate of the modal and spatial distribution of exportations from local center i is given by:

i ij a i M m im am ji J j A a C       Im Im  (10) (11) (13) (12) (14) ) / ( max am ij m ij am ij m ij am ij m ij M m J j A a P P Log P i Ex a i am ij            i i b i ba a b A b i a am ij M m J j A a V V i i Ex a             , i im Ex am ij J j A a M m C Ex a i        i i Ex a am ij M m A a J j     , , 0  4. LOCAL DEVELOPMENT PERFORMANCE INDICATORS
Following the above model of regional economics and transportation, we get different and, in general, competing performance indicators:

The net revenue of activity sector a at local center i is given by: 0

) ( ) ( Im ia am ji m ji J j am ij m ij J j M m i a i a a i a P P V p R a Ex a i                i A a 
where a p is the generalized price of the goods, i a  is the average variable cost and 0 ia  is the fixed cost of production of goods in activity sector a of the local center i.

The net revenue of transportation mode m is given by

pax im R + fre im R
where : Since m 0 is the index assigned to air transportation activity, here part of the cost of modality m 0 at local center i is the main component of airport revenue which is given by: (15) Its net revenue is given by:

ij J j pax m i arm kj pax ijm m ij I k J j R r A a pax im M m c T c R ij               0 ) ( I i  I i M m c c P R ij J j fre m i am ij am ji fre ijm m ij J j A a fre im             0 ) ( ) (  
(17) (16) (18) (19) (21) (20)     ) 1 ( ) ( ) )( 1 ( ) 1 ( ) 1 ( 0 0 0 0 0 0 0 fre v am ji am ij fre ijm arm kj R r I k pax ijm J j pax v fre f fre m i pax f pax m i aerop i c T c c c R ij                     
)) ( ) ( ( 0 0 am ji am ij fre iv arm kj R r I k J j pax iv aerop if aerop i c T c c R ij            
where it is supposed that the fixed cost of airport i is such as: Then finally, the distribution of revenue among the population can be related to activity levels by a relation such as:

) ( 0 0 0 Im im Ex im fre i m i pax i aerop if C C C c     
R r A a I i V Pop i i a i ar i ar i ar i      , , 0   

IMPACT EVALUATION FRAMEWORK OF A NEW REGIONAL AIRPORT

The implantation of a new airport in a local center which had to use other transportation modalities including remote airports at other regional centers will have many impacts on local development:

-Transportation costs and volumes will be redeployed among the different modalities;

-Sectorial activity levels will be redefined; [START_REF] Jain | A New Broad Composite Model of Regional Economic Development as Framework for Research[END_REF] (24) [START_REF] Ishutkina | Analysis of the interaction between air transportation and Economic Activity: A worldwide Perspective[END_REF] (25) -Import/export flows will be redistributed, -The revenue distribution among the population will be modified.

The above transportation and sectorial 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        
with a distribution of revenues given by:

R r A a Pop i i ar i   ,



Then, introducing the new airport facility at regional center i, to get the same level of activity, importation and exportation flows can be considered to be solutions of the following problems:

) / ( max Im am ji m ji am ji m ji am ji m ji M m J j A a P P Log P i a i am ji            i i r i ar i R a R r i a aa am ji M m J j A a Pop V i i a                 Im Im 0  am ji  i i a M m A a J j    , , Im and ) / ( max am ij m ij am ij m ij am ij m ij M m J j A a P P Log P i Ex a i am ij            i im Ex am ij J j A a M m C Ex a i        i i Ex a am ij M m A a J j     , , 0
 while travel demand will be now such as:

i ij a i M m im am ji J j A a C       Im Im  i i b i ba a b A b i a am ij M m J j A a V V i i Ex a               , ( 27
) (26) (28) (29) (31) (32) (33) (35) (34) (36) (30) ) / ( 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               R r A a I i D O T i pax iar pax iar arm ij M m J j ij          , ,   I i C T pax m i arm kj I k J j R r A a ij i           0  arm ij T
Then, the net revenue of activity a will be changed by:

     ) ) ( ) ( 0 Im 0 0 0 Im 0 i a i i Ex a i ij a i ij Ex a 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 P P P P R                              i A a 
The variation of net revenue of transportation mode m ( m 0 ) is given by

pax im R  + fre im R  where:   0 ) ( ) ( i ij J j arm kj arm kj pax ijm m ij I k J j R r A a pax im m M m T T c R ij ij                   I i      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 ) ( ) ( ) (       
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 new airport 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. When the activity levels saturate some capacity constraints at solution of problems (28-31), (32-35) and (36-39) it will be interesting to restart the analysis with increased capacity levels.

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 an inter temporal analysis where the lower level (demand forecasting) and the upper level (activity levels) interact along the chosen horizon of time.

ILLUSTRATION OF PROPOSED APPROACH

Here we considered the application of this framework to the implantation of a medium regional airport in a fast developing rural agro-industrial area in the center of Brazil (Western Bahia State). The considered region is shown in map 1:

Map 1. National position of study area

The area under study consists of a set of six local centers located near a main highway (BR242) and the planned LESTE-OESTE freight railway has been considered in the study.

Only one short runway airport is already operating in the sub region with turbo-powered light aircraft to connect this area with two larger regional airports in Bahia state where turbo jets link them to major regional and national centers (Salvador, Ilhéus, Brasilia, Belo Horizonte and Recife). The planned airport is supposed to operate light aircraft and turbo jets directly towards these regional and national centers on a weekly basis (turbojets) and on demand (light aircraft).

Three main activities have been considered: agro-industry (1.5 B US$, 2009), mineral extraction (0.6 B US$, 2009) and tourism (100 MM US$, 2009). The total population of the sub region is about 300 000 inhabitant with an average household annual revenue of 6000 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 dimension, so that many numerical techniques can be applied to solve them with a limited computational effort allowing the study of differentiated scenarios.

CONCLUSION

This communication has considered the problem of long term forecasting of passengers and freight transportation in a large network and its impact when a new air transport infrastructure is created or improved. The proposed approach has introduced two different optimization models: A model devoted to demand forecasting and another one defining the economic performance of local activity sectors and transportation modalities. The prediction of generated demand of passengers and freight process is based on an entropy maximization approach with flexible origin-destination levels to determine the intensity and the distribution of new origin-destination vectors.

Then, a bi-level approach considering passengers and freight flows at the first level and activity flows at a second level has been developed with the objective of providing a tool for the appraisal of different regional development scenarios related with air transportation investment.

DV

  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 local center i, i ar  is the percentage of population of local center i involved in activity a with revenue level r (Index 0  a for i ar  represent residual generation of trips from people of revenue level r, I i  , I being the set of regional urban centers considered) and O ar  is the corresponding trip generation factor.and the attraction model: is the annual passenger attraction to urban center i with respect to people related with activity a, is level of activity a at a location i producing attraction at a i ar  rate, to outside people of revenue level r.

  where the positive coefficients i ba 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 local center i is given by:



  ˆ is an a priori distribution and cost that are available for the modality m. Appling the same argument, am ij ˆare used for the exported goods.

  respectively the fixed cost of transportation of passengers and goods which are used in the local center i with the modality m.



  are the proportions of airlines fixed costs corresponding to airport taxes while pax v  and fre v  are the proportions of airlines average variable costs corresponding to airport taxes, either for passengers or goods air transportation.

  imports and exports.

The case study considered of rural a area in the center of Brazil, characterized by agro-industrial activities, shows that from the points of view of statistics availability and computational complexity, the proposed approach is feasible.