NUMERICAL SIMULATION OF LIQUID METAL INFILTRATION AND SOLIDIFICATION INSIDE A CAPILLARY TUBE

A metal foam is a porous structure whose solid matrix has a large fraction of interconnected cells. The objective of our study is to define a new manufacturing process via casting that produces homogenous open cell metal foams. This comes down to studying the infiltration and solidification process of a liquid metal inside a porous mould. The metal foams are characterized by their high porosity and permeability values; thus, the size of the mould pore is small enough to be considered as a capillary tube. Therefore, the focus of this paper is on modelling and numerical simulation at a local scale of the infiltration and solidification of liquid metal inside a capillary tube.

In our future work, we will derive a macroscopic model at the mould scale by up-scaling: the volume-averaged equations will be obtained from the local scale model using the volume averaging method [5]. Then, the numerical results will be compared to the metal foam casting tests.

The main difficulty lies in the presence of 4 phases; mould, air, liquid and solid metal: a one-domain approach is chosen in the modelling of the problem. We subsequently obtained one set of equations as a function of two volumetric fractions: the metal/air fraction in the pore and the liquid fraction in the metal. Before the solidification starts, the liquid metal infiltrates the capillary thus the liquid metal-gas interface must be tracked. This is done implicitly using the VOF (Volume Of Fluid) method [START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF], where the advection equation of the volumetric metal fraction is coupled to the Navier-Stokes equation. The solidification process of the liquid metal is accounted for using the enthalpy-porosity approach [START_REF] Brent | The Enthalpy-Porosity Technique For Modeling Convection-Diffucion Phase Change : Application To The Melting Of A Pure Metal[END_REF][START_REF] Rösler | Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments[END_REF]. In addition, to ensure a zero velocity field in the mould and in the solid metal, two penalty coefficients as a function of the phase indicators were added to the momentum equation.

The above mathematical model is implemented in a CFD tool: OpenFOAM [4]. A numerical validation is performed by comparing the numerical results with well-known solutions of test cases. Then, the numerical results of the liquid aluminium solidification are compared to the fluidity tests carried out experimentally. Subsequently, a numerical parametrical study enabled us to find the relationships between the distance flowed by the liquid before its solidification and the infiltration time as function of the pressure difference, heat transfer coefficient and metal superheat.

INTRODUCTION

The competitiveness in the metal industry result of severe products specifications in order to reduce the cost of the manufacturing processes. Thus, metal foam presents a promising material that keeps the high mechanical properties of the metal while reducing the weight up to 90%. Metal foam has variety of applications depending on whether the cells are connected or not; e.g. open-cells foam are used in heat exchanger to increase heat transfer, while closed-cells foam are employed as impact-absorbing material such in vehicle's crash box. Several patented manufacturing processes have been introduced during the last two decades. We count: Alporas by Shinko Wire Co. Ltd. [START_REF] Miyoshi | ALPORAS Aluminum Foam: Production Process, Properties, and Applications[END_REF], Alusion TM by CYMAT [7], Duocel foam by ERG aerospace cooperation [8]. However, the metallic melts are foamed or by injecting gas into liquid or by the use of blowing agents and normally these processes need a stabilizing mechanism. As a result, the products are not homogenous, thus the geometry and physical properties cannot be controlled. Therefore, the CTIF[9] presented a manufacturing process of metal foam via casting, i.e. CastFoam[10] that produce a well defined homogenous foams. However, this process needs to be elaborated before being employed in the metal industry. Thus, the infiltration and solidification of liquid metal inside the porous mould should be studied. Metal foams obtained using this process, have high permeability and porosity values; hence the size of the mould pore is small and can be modelled as a capillary tube. Subsequently, our objective is to study the infiltration and solidification of liquid metal in a capillary tube.

NOMENCLATURE

MATHEMATICAL MODEL

The liquid metal infiltration and solidification inside a capillary tube is handled by solving the governing equations presented in this section, using a one-domain approach. The fluid flow is governed by the Navier-Stokes equations for an incompressible fluid. Therefore, the continuity and the momentum equations are given by: 0

) .(      v   t (1) v v v vv v B A K p t                   ) ( ) .( (2)
With               K
, curvature of the metal-gas interface

In two-phase flow model the momentum equation contains an additional source term which is the pressure jump at the liquid-gas interface arisen from the surface tension between the two fluids. Using the volume-of-fluid (VOF) method [START_REF] Rusche | Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions[END_REF] this term that is the third term in the r.h.s. of equation ( 2) is evaluated by resolving the transport equation of the metal volumetric fraction :

  0 ) 1 ( . ) .(             r v v t (3) With g l r v v v  
, the relative velocity

The surface curvature is function of volumetric fraction and is required for the determination of surface tension force. Therefore, the calculation of the volumetric fraction distribution must be accurate especially in case of high physical properties ratios i.e. density, thermal conductivity..., where small errors in evaluating the volumetric fraction leads to significant errors in calculating the effective physical properties. Thus, equation ( 3) contains an additional convective term known as the compression term active only within the interface region ) 1 (    ; its role is to get a sharper liquid-gas interface. As we can see, the momentum equation contains also two supplementary artificial terms, which will be discussed later on after presenting the solidification model.

The modelling of convection-diffusion phase change of the metal is done using a modified enthalpy-porosity method [START_REF] Brent | The Enthalpy-Porosity Technique For Modeling Convection-Diffucion Phase Change : Application To The Melting Of A Pure Metal[END_REF]. The advantage of this method is that it respects the one domain approach. The method is modified by writing the energy equation as function of temperature instead of enthalpy and is extended by a continuous liquid fraction function, i.e. error function [START_REF] Rösler | Shell-and-tube type latent heat thermal energy storage: numerical analysis and comparison with experiments[END_REF]. The energy conservation equation is written as follow:

) ( ) .( T k H t H        v (4)
The latter contains two unknowns, fortunately the enthalpy can be written as function of the temperature and the volumetric fractions:  the metal/air fraction, g l the liquid fraction in the metal and ε the fluid phase indicator.

                                    T T m m T T g g T T s s l l l l T T l l ref ref ref ref d cp d cp cp g L g d cp g H             ) 1 ( ) 1 ( ) 1 ( (5)
Replacing the enthalpy with its expressions in equation ( 4), and neglecting the density change between the liquid and solid phases of the metal, we get the following energy equation:

                               l l T T s l metal g t g L d cp cp T k T cp t T cp ref v v     ) ( ) .( (6)
As we mentioned earlier, the liquid fraction is considered to be a continuous function of temperature given by:

2 1 ) ( 8 2 1            fusion elt m l T T T erf g (7)
With elt m T and fusion T being reciprocally the arithmetic mean and the difference between liquidus and solidus temperatures of the metal in question. After replacing the liquid volumetric fraction g l with its expression in equation ( 6), the energy conservation equation is function of one unknown, the temperature.

In a one domain approach, the four phases; mould, air, liquid and solid metal are considered as one effective fluid throughout the domain having physical properties weighted as function of the phase fractions.

g metal     ) 1 (    (8) g l     ) 1 (    (9)       m m g g s l l l metal p cp cp cp g cp g c        ) 1 ( ) 1 ( ) 1 (        (10)
These effective properties are obtained by supposing a linear variation between the physical properties of each phase at the different interfaces. However, this simple relation cannot be used when calculating the effective conductivity because of the high ratio between the metal and air conductivity values. Therefore, using Patankar [START_REF] Suhas | Numerical Heat Transfer and Fluid Flow[END_REF] formulation, the effective conductivity is given by:

    m g s l l l k k k g k g k ) 1 ( ) 1 ( ) 1 ( 1                     (11)
Now, reviewing the momentum equation, two sources terms were added, they were derived from the Carman-Kozeny equations for flow in a porous media. Therefore, the condition of zero velocity fields in the solidified metal is verified using the penalty coefficient A that acts only in metal region and is function of the liquid metal volumetric fraction:

b g g C A l l    3 2 ) 1 (  (12)
When the calculated temperature T is higher than the melting temperature, the liquid fraction is equal to 1 and the penalty coefficient vanishes, thus the momentum equation returns to its normal form. However, the calculated temperature is lower than the melting temperature, solidification process occurs, and the liquid fraction is equal to zero. Thus, the penalty coefficient has a large value that multiplies the velocity vector which result in a zero velocity field. The constant b has a small value; its role is to avoid division by zero.

The second additional source term in the momentum equation, is to ensure a zero velocity field in the mould region. Thus, the penalty coefficient is function of the fluid phase indicator:

D B ) 1 (    (13)
In the fluid region ε equal 1 and the penalty coefficient turns to zero. However, in the mould region ε equal 0 and the penalty coefficient has large value multiplying the velocity vector therefore insuring a zero velocity solution.

NUMERICAL VALIDATION

The mathematical model presented earlier was implemented in a CFD open source software OpenFOAM [4]. The latter is written in C++ for linux operating systems, and contains essential codes to solve flow and heat transfer problems. It uses the finite volume method where the space is divided into discrete control volumes, and time is split into a number of time-steps. The solution procedure of the transition simulations can be summarized as follow: first the mesh is generated, and the variable fields are initialized. After that, the time loop starts, and the time step is calculated in order to verify the Courant-Friedrichs-Lewy (CFL) condition. Then, the transport equation of the metal volumetric fraction ( 3) is solved which enables us to calculate the surface tension force. Subsequently, the metalgas interface is constructed and the physical properties are updated as function of volumetric fractions (equations (8, 9, 10 & 11). Next, the pressure-velocity coupling is done using the PIMPLE (PISO-SIMPLE) scheme of OpenFOAM where the continuity and the momentum equations (1, 2) are solved. Afterwards, the convergence of the solution is monitored by calculating the continuity error. At this stage, the pressure and velocity fields as well as the metal volumetric fraction distribution are known for the current time step. Using the velocity field solution, the energy equation ( 6) can now be solved in order to obtain the temperature field and the metal liquid volumetric fraction g l , which will be used to update the penalty source term of the momentum equation. Now all the variable fields are known, if the final time is not reached yet, the time loop continues with the calculation of the next time step.

The numerical validation of the code is done by comparing the numerical solutions with analytical solutions and correlations of test cases available in the literary. In this section, the mould is not taken into consideration. Therefore, the capillary tube is modelled as a cylindrical tube of 1.2 mm radius, and 120 mm length. The mesh is structured 2D radial quadrangle extruded to 3D using Salome software[12], and is shown in the figure below.

Figure 1 Global view of the geometry and mesh of the tube

One phase flow

The first step of our code validation is by simulating an isothermal one phase flow and comparing it to the Poiseuille solution. This is done by fixing the right thermal conditions. Thus, the tube is initially filled with liquid aluminium at a temperature higher than its fusion temperature and a zero thermal flux is set at the tube wall. A Dirichlet boundary condition of constant uniform velocity profile is imposed at the inlet and an absolute pressure value of 0 Pa is fixed at the outlet. The adimensional number of Reynolds is taken to have a value of 600, which corresponds to an average velocity of 0.1 m/s. The entrance length after which the flow is fully developed is calculated as follow:

D L m Re 0576 . 0  (14)
The radial velocity profile at an axial section farther from the entrance length (x=100mm) is compared to the Poiseuille solution and shown in Figure 2. Both solutions are in good agreement, and the relative error at the tube axe is less than 1%.

Figure 2 axial velocity profile at x=100mm

For an internal forced convection, there are two possibilities of thermal conditions at the tube wall: fixed temperature or constant thermal flux, accordingly to which the thermal entrance length and the correlated dimensionless number of Nusselt are given. Using the aluminium properties, the thermal entrance length for both cases is negligible in comparison with the mechanical one. In addition, the high thermal conductivity of the aluminium yields to a uniform radial thermal profile in the case of fixing a constant temperature at the wall. Therefore, the validation of the thermal behaviour is done by fixing a constant flux at the wall. In order, to identify the value of the flux that must be used as a boundary condition; the energy that should be extracted from the metal in order to lower its temperature from in T to out T is calculated. For a constant heat flux at the wall, we have [START_REF] Padet | Convection thermique et massique Nombre de Nusselt: partie 1[END_REF]: In the previous thermal simulation the heat flux imposed at the tube wall, was sufficient enough to lower the liquid aluminium temperature in order to calculate the Nusselt number once the mechanical and thermal regime are established but not high enough to induct the solidification of the metal liquid. Now, in order to observe the solidification phenomena of the liquid flowing inside the tube, the heat flux imposed at the tube wall is much higher. The tube is initially filled with liquid aluminium at a temperature of 1073 K, driven by a pressure difference of 300 Pa between the inlet and outlet. The evolution of the liquid fraction is shown in the pictures below: The solidification process of the liquid aluminium cannot be validated by an analytical solution, but the numerical results show a correct physical behaviour.

      

Two phase flow

At this stage, the solidification code should be tested in simulating the infiltration and solidification process of liquid metal into the capillary tube. For an isothermal two-phase flow, special attention is given to the interface velocity. In fact, after a certain time the interface shape remains the same. This means, that the interface can be considered as a solid part where all the particles move at the same velocity, which corresponds necessarily to the fluid mean velocity; as Dussan states in her article [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF] "the speed of the contact line 2 a Q  ". In order to understand the motion of the two fluids and the interface dynamic we refer to the article mentioned earlier [START_REF] Dussan | On the spreading of liquids on solid surfaces: static and dynamic contact lines[END_REF]; where it is demonstrated that the displacing fluid undergoes a rolling motion and that the displaced fluid has a similar but more complex behaviour. To visualize this, the relative velocity vectors in the interface frame of reference are shown below: For an internal convection two-phase flow, the tube wall is considered as adiabatic, and the two fluids have different temperatures that are higher than the solidification temperature of the aluminium. The liquid aluminium infiltrating the tube cools down by losing its energy to the air that has a lower temperature.

Figure 6 Axial temperature profile

As mentioned earlier, the solidification cases cannot be validated by analytical comparison. However, it is of our interest to observe the solidification process of the liquid aluminium infiltrating the tube. That is why, in this section we are going to show the solidification behaviour of the liquid aluminium under an imposed negative flux boundary condition. In aim to eliminate the air temperature effect on the solidification kinetics of the liquid aluminium, the latter is taking to be equal to the liquid aluminium temperature T=1073K.

RESULTS AND DISCUSSIONS

  profile obtained at (x=100mm) is shown in Figure3and gives the following Nusselt number: the correlated Nusselt number yields to a relative error of 1.45%.

Figure 3

 3 Figure 3 Temperature radial profile at x=100mm

Figure 4

 4 Figure 4 solidification process of liquid aluminium (in red liquid phase, in blue solid phase); a-at start time the tube is completely filled by liquid aluminium, b-the aluminium solidification starts at the tube wall at a distance from the entry, c-the solidified aluminium progress in the radial direction, dthe solidified aluminium close the tube at a section further from the entrance and stops the fluid flow, e-the solidification progress at the entry direction, f-at the end the aluminium is completely solidified

Figure 5

 5 Figure 5 Relative velocity vectors in the region of the liquidgas interface Figure 5 shows that the displacing fluid flows forward at the tube centre, and then rolls at the interface vicinity to flows backward near the wall. At the other hand, the displaced fluid flows backward near the wall, and rolls at the interface to flows forward at the tube centre.For an internal convection two-phase flow, the tube wall is considered as adiabatic, and the two fluids have different temperatures that are higher than the solidification temperature of the aluminium. The liquid aluminium infiltrating the tube cools down by losing its energy to the air that has a lower temperature.

Figure 7

 7 Figure 7 liquid aluminium solidification progress; a-the liquid aluminium infiltrates the tube, b-the solidification starts at the wall behind the interface, c-the solidified aluminium closes the flow passage, d-the solidification propagates in both directions, the interface oscillates, e-at the end the aluminium is completely solidified


Siim-Toomas Marran

Artjom Lind

Amnir Hadachi

Introducing Cellular Network Layer into SUMO for Simulating Vehicular Mobile Devices' Interactions in Urban Environment

Keywords: Simulation, Mobility Big Data, Call Detail Records, Cellular Networking, Urban Road Traffic, GPS Devices, Mobile Devices, SUMO

During the last decade researchers have been demonstrating the importance of mobile data or CDR data in depicting the human mobility patterns. However, this type of data is not easy to get access to from mobile operators. Besides, in order to make this type of data available and enable their usage for the scientific communities the process can face many constraints that can constitute obstacle. From this perspective, this paper introduces a way to produce realistic real-life mobility logs through the traffic simulation tool SUMO, which has been enhanced with a cellular network layer to mimic cellular networking behavior.

INTRODUCTION

The mobility data is a broad connotation that refers to a various types of data, which describes people's movement and activities. Some of the common types of the mobility data are Call Detail Records (CDR) [START_REF] Hadachi | Cell phone subscribers mobility prediction using enhanced markov chain algorithm[END_REF], Global Positioning System (GPS) logs [START_REF] Hadachi | Approach to estimate travel time using sparsely sampled gps data in urban networks[END_REF], Vehicular ad hoc network (VANET) logs [START_REF] Karnadi | Rapid generation of realistic mobility models for vanet[END_REF], Radiofrequency identification (RFID) system logs [START_REF] Finkenzeller | Rfid handbook: fundamentals and applications in contactless smart cards, radio frequency identification and near-field communication[END_REF], WiFi access points related data [START_REF] Bulut | Wifi access point deployment for efficient mobile data offloading[END_REF], etc. Moreover, with all the advancement in information and telecommunication technologies (ICT). Mobile data has a great potential in sensing the urban mobility dynamics. From this perspective, it is clear that mobile data or CDR data are valuable. In general, CDR are data records, produced inside of the cellular network between the user equipment and the base station, documenting various telecommunication transactions. Besides, it is estimated according to (Statista, 2017) that 4.77 billions mobile users produce daily unimaginable amount of data, with restricted access. The mobility big data has been used for developing many different algorithms and applications such as triangulation and trilateration, Kalman filter for mobile positioning [START_REF] Lind | A new approach for mobile positioning using the cdr data of cellular networks[END_REF] , handling the crowd management [START_REF] Pan | Crowd sensing of traffic anomalies based on human mobility and social media[END_REF] , deploy real-time incident detection solutions [START_REF] Zaldivar | Providing accident detection in vehicular networks through obd-ii devices and androidbased smartphones[END_REF] , manage and avoid congestions in the traffic ,develop Vehicle-To-Infrastructure (V2I) applications , bind the mobil-ity with various Internet of Vehicles (IoV) and Geographic Information Systems (GIS) to achieve new systems upgraded with the human activity related functionality. As a consequence, there occur moments when it is not feasibly easy to acquire real data for specific purpose, we rely on simulating and generating synthetic data based on mathematical models. For example, we have open source SUMO simulator (DLR, 2017) that provides the ability to simulate real-life urban traffic at microscopic level. This open source simulation package provides realistic simulation due to the OpenStreetMap (OSM) (OSM, 2017) map import functionality and in-built configurable behavior models. The routable entities can be equipped with various devices. In this paper, we are integrating we are integrating SUMO with mobile and GPS devices. Hence, we introduced also a cellular network layer with all its functionalities linking the routable devices and the network. This allows us to generate from the microscopic traffic simulator the associated CDR and GPS logs, which can be later used for the training and testing various machine learning models in the real life context.

RELATED WORK

SUMO is a suite of applications to prepare and to simulate various real-life traffic scenarios. Therefore, it has been used for many intelligent transport systems related research topics. For example, the project TrafficOnline is a good illustration of the use of SUMO for simulating the usage of mobile phone data, where the real-world GSM data was used to develop a telephony model (consisting of the following properties: call start, call duration), to determine the travel times [START_REF] Krajzewicz | Recent development and applications of sumo-simulation of urban mobility[END_REF]. However, the dynamic properties, e.g. cell size variations, cell changing mechanics, of the GSM network were not considered. Another related research topic is VANET simulations. The researcher tried to create an integration between the traffic simulator SUMO and third party network simulators, such as OMNet++ [START_REF] Omnet++ | OMNeT++ Discrete Event Simulator -Home[END_REF], where they have created bidirectionally coupled hybrid simulators and Viens [START_REF] Sommer | On the need for bidirectional coupling of road traffic microsimulation and network simulation[END_REF]. Besides, the list is long for example there are also NS-2 (ns 2, 2017) and NS-3 (ns 3, 2017) that have been used for VANET system's for detailed and realistic performance evaluation. Similar to Veins, TraNS (Traffic and Network Simulation Environment) based on NS-2, was introduced to combine two disjointly developed research simulators into one in order to deal with VANET performance evaluation [START_REF] Piorkowski | Trans: realistic joint traffic and network simulator for vanets[END_REF]. The NS-3 was introduced as replacement to NS-2 and it was used in the VANET crossroad scenario to evaluate the performance of HWMP, OLSR and DD routing protocols, while using IEEE 802.11p standard and TwoRayGround Propagation Loss Model, to send multiple Constant Bit Rate (CBR) flows over UDP between 20 source-destination pairs. [START_REF] Kolici | Performance evaluation of a vanet simulation system using ns-3 and sumo[END_REF] Furthermore, we have the emergence of the vehicular networking concept, where the vehicle is connect to everything. Therefore, the birth of many type of communications related to VANET such as Vehicle-To-Everything (V2X), vehicle-to-vehicle (V2V), and vehicle-to-infrastructure (V2I). The main purpose of introducing the connectivity into the vehicle to solve issues related to traffic control and management via the share or collection of reliable information from the vehicles. For this purpose, many simulators have been created such as VSimRTI [START_REF] Wedel | V2x-based traffic congestion recognition and avoidance[END_REF]. Obviously VANET solutions aren't only viable solutions to tackle the road traffic estimation problem. Due the widespread of the mobile communications we can estimate the traffic density on the bigger and more important roads with the help of cellular network, since many commuters are carrying mobile devices with them. Such approach has proved to be accurate and very capable to detect and quantify state changes. [START_REF] Bolla | Road traffic estimation from location tracking data in the mobile cellular network[END_REF] 3 SYSTEM DESIGN AND ARCHITECTURE

The purpose of this simulation is to produce CDRlike activity logs with the help of microscopic road traffic simulator SUMO, where we added an extra layer -cellular networking. This integration of a new layer introduces to SUMO's core and support packages many changes and new features.

Models

The most important step in the process of integrating the cellular network into SUMO is to define the models of each new component.

Core mechanics of the simulation

The core mechanics of the road traffic simulation SUMO has been gathered into one class which contains the road network and orchestrates the simulation. This class is initialized during the start of the program and is given to the network builder as a reference to create the infrastructure for the simulation.

During the road network building, the cellular network is being loaded into the simulation network.

Then, when the simulation performer has all the vital dependencies, then the simulation cycle starts. The core of the simulation contains of multiple event controllers, which are being ran after every tick. The simulation have been divided into four components according to their event controller action types and intermediate smoothing and they are as follows:

Begin of time step events is the original event controller in SUMO, which, as an example, takes care of traffic light events, trigger type events(lane speed, rerouting), routing device, and pedestrian move events. Intermediate actions, depends on the previous event control cycle and smooths the network for the next event control cycle. During this process, the system is detecting regularly for collisions, checking the traffic lights, checking if the edges are active, plans the vehicles movement, executing the vehicles movement, run lane changes, load new routes and mobile events for the cellular simulation, insert new vehicles to the network, insert new events, etc. During the execution of vehicles movement, there is a case when the vehicle has arrived to the destination and, immediately afterwards, the vehicle shall be removed from the simulation. This function has been altered to have checks for existing mobility devices (GPS, mobile) to run their last events immediately and then invalidating those events to stop them from recurring. This entire code alternation prevents losing event data due the premature removal of the vehicle. End of time step events is an event running cycle which runs the events dependent on the updated location of the routable elements. There are also step by step logger events, vehicular devices specific events (Bluetooth, GPS, mobility management related state update) and vehicle flow calibrator events. Mobility events is a newly added event running cycle entirely for the purpose of mobility event simulation. The events are created by the mobile devices from the data, which has been attached to the vehicle after it has been loaded into the simulation during the intermediate actions.

Data import and export

The simulator depends on the data which has been prepared beforehand for the simulation and therefore imported into the application through Extensible Markup Language (XML) files. SUMO has in its package multiple programs which help to produce the network and other simulation related data in XML document format to run the simulation properly. For the cellular networking layer, there were created multiple XML data import handlers and their respective XML Schema Definition (XSD) files. First, cellular network related data like antennas and their respective transmitters with coverage areas. Second, mobile events data, which were described in the Section 3.1.

For the export there were designed output XML formats for the GPS and mobile devices. When the spe-cific routable entity finished the navigation in the network then all its devices data were exported into the resulting XML file.

SIMULATOR IMPLEMENTATION

The implementation of cellular network log generator and it's integration into SUMO base code was done based on introducing two prerequisite applications: Mobility Event Simulation Generator (MES-GEN) and Cell Coverage Area Generator (HEXA-GEN).

MESGEN

Mobility Event Simulation Generator (MESGEN) is a supplementary application required by the main simulator application. MESGEN is responsible for generating the cellular events-time-line for each vehicle (cellular device inside vehicle). SUMO simulator uses the cellular events-time-line to trigger the events in certain cell towers during the simulation when vehicle traverses the corresponding coverage area. Several cellular activity profiles have been defined and are illustrated in Table 1. Each profile is specified by three attributes: frequencies of calls, SMS or and data usage activities. The events generated by MES- Listing 1 The appearance of the MESGEN output file <vehicles> <vehicle id="0" depart="0.00"> <mEvents> <mEvent begin="114.00" end="137.00" type="CALL_REC"/> <mEvent begin="363.00" end="450.00" type="WEB_COMM"/> <mEvent begin="490.00" end="516.00" type="CALL_INIT"/> <mEvent begin="627.00" end="669.00" type="CALL_INIT"/> <mEvent begin="731.00" end="742.00" type="CALL_REC"/> <mEvent begin="829.00" end="856.00" type="CALL_INIT"/> </mEvents> </vehicle> </vehicles>

Minimum length and maximum length of different events can be configured through constant variables. It is also possible to define time buffers between occurring events. The values of low and medium frequency events thresholds and random number seed can also be set through the configuration file.

HEXAGEN

HexagonGen is an accessory level script, written in Python, to generate 120 °sectorized cells for the cellular network behaviour simulation. It takes an input of latitude and longitude coordinate values, size (length of hexagons side), hexagon network width, and height values. The given coordinate pair is the center of the drawn network (tower position). The result of the HexagonGen is an XML file which contains the cell coverage elements and can be fed directly into SUMO.

Algorithm 1 Mobile's events generator

Precondition: V ∈ N set of vehicles, P ∈ N 3 set of profiles, T ∈ N vehicle departure times, l call maximal call duration, lsms maximal interval between sending SMS, l data maximal web surfing duration, E set of event types as was specified in formed. Figure 1 illustrates how mobile device component is integrated into SUMO and how does it interact with other modules.

Mobility management

During the simulation the state of mobile device is managed by MobilityManagement. The state is updated at the end of each simulation step event. Vehicles are initialized based on the parameters file, which specifies what devices are equipped (mobile, GPS, etc). Mobility management entity is only created in case of vehicle being equipped with a mobile device. Mobil-ityManagementUpdate is an extension of SUMO Command base abstract event class and is responsible for triggering the state update. The new state is assigned based on previous one, the state transitions are illustrated in the Figure 2. Mobility management has a reference to the mobile device to call out mobility related procedures from the lower layers. Then, there are mobility management state related variables; the T312 timer is for reminding a mobile device to make periodic location updates after every user defined steps (by default 10 steps). Some state related extra scenarios, define whether the mobile has been turned on, the SIM card inserted, the IMSI attached, and define whether the mobile has PLMN related information. In the Figure 2 a mobile mobility life cycle is described visually for better understanding. At the start of the mobility management update, the creation of the current time step's mobility events for the simulation for this specific device takes place. During the mobility management update event, there will be a call out to the management cores update function which increments instantly the periodic location update with respect to the previous state that goes through state machine to determine the next steps: a) If the mobility management state is classified as NULL the core update will change mobile to the state as it was just MM CONNECTION ACTIVE starts the cell selection process for the handover purposes, if needed.After cell selection, the work-flow will move from the device back into the mobility management, where the state of MM will be updated and location update initiated. If the state machine has done its work, then there is a successive check to determine if the periodic update is needed and, if it is, then it will be executed.

The rule with the periodic location update is that if a mobile has not had any location updates during some certain amount of simulation steps, then the mobility management has to step in and remind the transmitter that it still exists and has not been lost during the commute.

Wave propagation

During the cell selection, the mobile device has multiple challenges. It has to determine the nearest polygons, determine in which one of them it resides, and determine the signal strength from the base stations transmitter. The signal strength has been envisioned by the free space wave propagation model:

P r (d) = P t G t G r λ 2 (4π) 2 d 2 L , (1)
Where P r (d) is the received power, P t is the transmitter output power, G t is transmitter gain G r is antenna gain, d is the distance between the MS and the tower, and L is the system loss factor. [START_REF] Nishith Tripathi | Cellular Communications: A Comprehensive and Practical Guide[END_REF]. This latter indicates the power received by the antenna under ideal conditions. in addition, the free space model predicts the powers decay to be the negative square root of the distance. In our case, we have not included the system loss factor. Besides, the distance in the equation 1 is calculated between the tower and the vehicle based on Haversine distance:

d = 2r Arcsin sin 2 φ 2 -φ 1 2 + f (φ 1 , φ 2 ) f (φ 1 , φ 2 ) = cos(φ 1 )cos(φ 2 )sin 2 λ 2 -λ 1 2 , ( 2
)
where φ is the latitude, λ the longitude and r is the Earth's radius in meters (6 371 000 meters).

During the signal strength calculation, we use the configuration supplied tower transmitter output strength (dBm) and the transmitter's frequency. From the frequency, we calculate the wavelength. Wavelength equals the speed of light divided by the transmitter emitted frequency.

Location update, Cell selection and reselection, Handover

The procedures are key component of keeping track of the mobility management level state, this can be seen also in the Figure 2.

Cell selection is a mobile device level method for measuring the current signal strength value of the serving antenna and manage handover procedure.

During the cell selection, the first step is to find sur- [START_REF] Guttman | R-trees: A dynamic index structure for spatial searching[END_REF] After searching the nearest polygons to the device, we check through the vector of cells and determine whether the vehicle is in the cell or not using the winding number algorithm. The vehicle might be in multiple cells; therefore, we are calculating the best signal strength and accordingly, we pick the transmitter and cell to camp on. In essence, the winding number algorithm is an algorithm which counts the number of times the polygon winds around the point of interest. If the point is not inside of the polygon, then the resulting winding number is 0. (Kai Hormann, 2001) During the cell selection, the registration to the tower or the deregistration from the tower will be determined. This is emulating the registering the mobile devices location area into the visitor location register. The final step is the status update of the mobility management.

Location update is the mobility management level procedure for updating the devices location in the network. Location update resets the periodic update timer (sets the timer back to 0) and updates the mobility management state.

Handover procedure is triggered by the state machine at the level of MM CONNECTION ACTIVE when there is a better quality cell flagged true.

Mobility events and the simulation cycle

We have covered the mobility management tied event object called MobilityManagementUpdate. In this section, we will discuss the mobile device related event.

MobilityEvent extends the SUMO Command abstract event base class. The mobility events have their own MSEventControl container-event queue and the events are being executed at the very end of the simulation step, after the mobility management update events.

Each mobility event has a type, which is declared in the enum of MobilityEventType in the class of MSMo-bilityEventData. The event has the mobile device reference to let the mobile know about the radio link failure and trigger the log creation about the occurred event. Other attributes of the mobility event are the events start time and end time. Calls and the web communication events depend on the timeframe variables.

Mobility CDR-like logs

Call detail record (CDR) is the information about incoming and outgoing mobile activities, e.g., calls or SMS messages. The data in the CDR is about the event originated and the terminated parties (both sides phone numbers), time of connection through the starting time and the call duration, call event type, unique generated id for the record, etc. [START_REF] Horak | Telecommunications and Data Communications Handbook[END_REF] Mo-bilityLog class is a data structure that, emulates the essence of CDR. The attributes of the mobile logs are illustrated in Table 3 (also visible in the Figure 1):

Results

The results are illustrated through the outcome of the logs generated by the simulator. The figure 4 is reflecting the interactions between the mobile devices, GPS devices and the cellular network.

GPS logs

MSDevice GPS related logs can be created in two ways: through main SUMO configuration file or through SUMO route file. Listing 2 illustrates the GPS signal log, it contains timestamp, longitude, and latitude.

Listing 2 GPS coordinates (longitude and latitude) with a timestamp <gps-output> <vehicle id="1"> <gps-device id="gps_1"> <signal time="1.00"> <coordinates latitude="58.29626476" longitude="26.44804043"/> </signal> <signal time="3.00"> <coordinates latitude="58.29628805" longitude="26.44810251"/> </signal> <signal time="5.00"> <coordinates latitude="58.29622327" longitude="26.44831755"/> </signal> </vehicle> </gps-output>

Cellular network simulation logs

The MSDevice Mobile class generates logs in case the corresponding mobile device has configuration values set to true, and the mobile data events have been generated by MESGEN. In addition, the cellular network coverage file (XML) should be specified. Otherwise there would be no connectivity for the devices and the events cannot occur.

The configuration of the mobile devices is similar to the GPS devices. One must define the mobile device into the route file or one can set all the vehicles to carry a mobile device in the SUMO configuration file.

In the results, the cellular network behaviour simulation written into XML file are illustrated in Listing 3. The logging entity is defined in the Figure 1 (class MobilityLog) and is explained in the Section 5.6.

Listing 3 CDR-like mobility logs <mobile-output> <vehicle id="1"> <event timestamp="92.00" eventType="KEEP_ALIVE" imsi="1" cgi="t-438-1" signalStregnth="-21.01" latitude="58.29" longitude="26.46292491" speed="13.68143750"/> <event timestamp="94.00" eventType="SMS_REC" imsi="1" cgi="t-438-1" signalStregnth="-21.70533381" latitude="58.29089506" longitude="26.46336419" speed="12.64468285"/> </vehicle> <vehicle id="27"> <event timestamp="27.00" eventType="CALL_INIT" imsi="27" cgi="t-852-1" signalStregnth="-21.11894073" latitude="58.25022813" longitude="26.43376626" speed="0.00000000"/> <event timestamp="28.00" eventType="CALL_INIT" imsi="27" cgi="t-852-1" signalStregnth="-21.21370266" latitude="58.25021502" longitude="26.43379995" speed="2.45816080"/> </vehicle> </mobile-output>

CONCLUSION

In this article, we are introducing a new layer into SUMO simulator. Our layer simulates the cellular network behavior and generates the logs of the interactions between the mobile devices and the mobile network. In addition, in order to make the integration with SUMO we created two main components: Mobility Event Simulation generator and mobile cell coverage generator. The results of our simulator is a set of mobile logs similar to CDR data that reflects the communication and the interactions between the vehicular mobile devices in the SUMO's vehicles and the cellular network. In general, the outcome is very encouraging and there is many enhancements that can be added; especially with in regards to the mobile network model for reflecting the mobile network protocols and behavior.

Figure 1 :

 1 Figure 1: The class diagram of MSDevice Mobile and its sub-entities.

Figure 2 :

 2 Figure 2: The simple work flow of the mobility management.

  Figure 3: Overview of R-tree processrounding cellular network cells. Those cells (polygons) are structured into the R-tree entity during the creation of the mobility tower control unit, where they were added for a faster search, as illustrated in the Figures3(a) and 3(b), R-tree is a tree data structure to handle spatial data efficiently by indexing shapes for the future access.[START_REF] Guttman | R-trees: A dynamic index structure for spatial searching[END_REF] After searching the nearest polygons to the device, we check through the vector of cells and determine whether the vehicle is in the cell or not using the winding number algorithm. The vehicle might be in multiple cells; therefore, we are calculating the best signal strength and

Figure 4 :

 4 Figure 4: 5G standard proposes use of micro-cells, therefore an example of HexagonGen generated micro-cells in SUMO.

Table 1 :

 1 User cellular activity profiles (counting amount of events per hour). The numbers in the table correspond to

	events frequency as follows: none (0), low (1), medium (2)
	and high (3).
	Profile Calls SMS Data
	Casual 1	1	1
	Only Keep Alive 0	0	0
	Business 3	1	1
	Teenager 1	3	2
	Student 1	1	2
	Talkative 2	1	1
	Media streamer 1	1	3
	GEN are similar to those of cellular event data model
	(original model was simplified by removing not used
	events). The MESGEN event types are illustrated in
	Table 1. The simplified flowchart of the cellular mo-
	bility management entity is illustrated in Figure 2. For
	example MESGEN events SMS SEND and SMS REC
	occur mimics the sequence of cellular network equip-
	ment activities that result in delivering an SMS mes-

sage (send or receive). MESGEN generates a history of events for each vehicle participating in simulation

Table 2 :

 2 MESGEN event type descriptions The result events are associated to a vehicle and have time specified as an offset of SUMO simulation life-time. An example of MESGEN output is illustrated in the Listing 1 having vehicles with associated cellular events. The pseudo code in the Algorithm 1 illustrates the algorithm used in MESGEN generator.

	Message Type Description
	SMS REC	SMS Received
	SMS SEND	SMS Sent
	CALL INIT	Calling out
	CALL REC	Receiving a call
	WEB COM	Uses mobile data
	using the simulation life-time, vehicle activity pro-
	files, road network information and cell coverage ar-
	eas.

Table 1

 1

	1: function TRIGGEREVENT(e,t, l)
	2:	Triggers event of type e ∈ E at t, lasting l sec.
	3: end function
	4: function RAND(n 0 ,n 1 )
	5:	n ∼ U([n 0 , n 1 ])	random number from (n 0 , n 1 )
	6:	return n
	7: end function
	8: function MESGENVEHICLE(v,tstart , t end )
	9:	n calls ,nsms,n data ← Pv	mobile profile of v
	10:	n all = ∑ Pv
	11:	t ← Tv	departure time of v
	12:	while t end -t ≤ 0 do
	13:	n = RAND(0, n all )
	14:	if n ≤ n calls then
	15:	l = RAND(0, l call )	Call duration
	16:	if RAND(0, 1) > 0.5 then	Call direction
	17:	TRIGGEREVENT(CALL INIT,t, l)
	18:	else
	19:	TRIGGEREVENT(CALL REC,t, l)
	20:	end if
	21:	t ← t + l
	22:	else if n ≤ n calls + nsms then
	23:	if RAND(0, 1) > 0.5 then	SMS direction
	24:	TRIGGEREVENT(SMS SEND,t,0)
	25:	else
	26:	TRIGGEREVENT(SMS REC,t,0)
	27:	end if
	28:	t ← t + RAND(0, lsms)
	29:	else
	30:	l = RAND(0, l data )	Web Surfing Duration
	31:	TRIGGEREVENT(WEB COM,t, l)
	32:	t ← t + l
	33:	end if
	34:	end while
	35: end function
	36: function MESGENALL(tstart , t end )
	37:	for all v ∈ V do MESGENVEHICLE(v,tstart , t end )
	38:	end for
	39: end function
	5 MOBILITY BEHAVIOUR AND
		MANAGEMENT
	5.1 Microscopic simulation devices
	SUMO provides basic interfaces to communicate with
	the vehicle and abstract methods, which need to be
	implemented to develop a working device.
	GPS device (MSDevice GPS class in the code) is a
	test device, it was created with the objective to study
	SUMO simulation cycle for the upcoming mobil-
	ity device. MSDevice GPS has two subclasses called
	RouteGPSInfo and GPSSignalUpdate. RouteGPSInfo is a
	data structure that records geographical coordinates
	and the timestamp. GPSSignalUpdate represents GPS
	signal, it is an extension of SUMO's Command class
	(base SUMO microsim event class). During every
	simulation step, if the vehicle has a device, the event
	is triggered the vehicles movement.
	Mobile device (MSDevice Mobile class in th code) is
	emulating multiple functions, which in real GSM

stack are performed by several logical units. Next section will explain how mobile device simulation is per-

Table 3 :

 3 Mobility Log Record Attributes

	Attribute Description
	timestamp Amount of seconds since
	simulation started
	event type CALL INIT, CALL REC,
	SMS SEND, SMS REC,
	WEB COM, KEEP ALIVE
	imsi subscriber ID
	(vehicle ID)
	cgi Common Gateway Identifier (CGI)
	the transmitter id
	(or simply a cell id)
	signal strength radio signal strength
	during the event
	latitude latitude coordinate
	(vehicle actual location)
	longitude longitude coordinate
	(vehicle actual location)
	speed vehicle speed
	during the event (in m/s)

Acknowledgement

This research work was supported by IUT34-4 "Data Science Methods and Applications" (DSMA) project.


Jilles S Dibangoye

Olivier Buffet

Learning to Act in Decentralized Partially Observable MDPs

Keywords: Decentralized Partially Observable Stochastic Control, Multi-Agent Reinforcement Learning Processes décisionnels de de Markov partiallement observables et décentralisés, Apprentissage par Renforcement Multi-Agents

We address a long-standing open problem of reinforcement learning in decentralized partially observable Markov decision processes. Previous attempts focussed on different forms of generalized policy iteration, which at best led to local optima. In this paper, we restrict attention to plans, which are simpler to store and update than policies. We derive, under certain conditions, the first near-optimal cooperative multi-agent reinforcement learning algorithm. To achieve significant scalability gains, we replace the greedy maximization by mixed-integer linear programming. Experiments show our approach can learn to act near-optimally in many finite domains from the literature.

Introduction

Decentralized partially observable Markov decision processes (Dec-POMDPs) emerged as the standard framework for sequential decision making by a team of collaborative agents [START_REF] Bernstein | The Complexity of Decentralized Control of Markov Decision Processes[END_REF]. A key assumption of Dec-POMDPs is that agents can neither see the actual state of the system nor explicitly communicate their noisy observations with each other due to communication cost, latency or noise, hence providing a partial explanation of the double exponential growth at every control interval of the required memory in optimal algorithms [START_REF] Hansen | Dynamic Programming for Partially Observable Stochastic Games[END_REF], [START_REF] Szer | MAA*: A Heuristic Search Algorithm for Solving Decentralized POMDPs[END_REF], [START_REF] Oliehoek | Optimal and Approximate Q-value Functions for Decentralized POMDPs[END_REF], [START_REF] Amato | Incremental Policy Generation for Finite-Horizon DEC-POMDPs[END_REF], [START_REF] Oliehoek | Incremental Clustering and Expansion for Faster Optimal Planning in Dec-POMDPs[END_REF], [START_REF] Dibangoye | Structural Results for Cooperative Decentralized Control Models[END_REF][START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. While planning methods for finite Dec-POMDPs made substantial progress in recent years, the formal treatment of the corresponding reinforcement learning problems received little attention so far. The literature of multi-agent reinforcement learning (MARL) can be divided into two main categories: concurrent and team approaches [START_REF] Tan | Multi-agent Reinforcement Learning: Independent vs. Cooperative Agents[END_REF], [START_REF] Panait | Cooperative multi-agent learning: The state of the art[END_REF].

Perhaps the dominant paradigm in MARL is the concurrent approach, which involves multiple simultaneous learners: typically, each agent has its learning process. Self-interested learners, for example, determine their best-response behaviors considering their opponents are part of the environment, often resulting in local optima [START_REF] Brown | Iterative Solutions of Games by Fictitious Play[END_REF], [START_REF] Hu | Multiagent Reinforcement Learning: Theoretical Framework and an Algorithm[END_REF], [START_REF] Littman | Markov games as a framework for multi-agent reinforcement learning[END_REF]. While concurrent learning can apply in Dec-POMDPs, a local optimum may lead to severely suboptimal performances [START_REF] Peshkin | Learning to Cooperate via Policy Search[END_REF], [START_REF] Zhang | Coordinated Multi-Agent Reinforcement Learning in Networked Distributed POMDPs[END_REF], [START_REF] Kraemer | Multi-agent reinforcement learning as a rehearsal for decentralized planning[END_REF], [START_REF] Nguyen | Policy gradient with value function approximation for collective multiagent planning[END_REF]. Also, methods of this family face two conceptual issues that limit their applicability. The primary concern is that of the co-adaptation dilemma, which arises when each attempt to modify an agent behavior can ruin learned behaviors of its teammates. Another major problem is that of the multi-agent credit assignment, that is, how to split the collective reward signal among independent learners.

Alternatively, the team approach involves a single learner acting on behalf of all agents to discover a collective solution [START_REF] Salustowicz | Learning Team Strategies: Soccer Case Studies[END_REF], [START_REF] Inria | When Evolving Populations is Better Than Coevolving Individuals: The Blind Mice Problem[END_REF]. Interestingly, this approach circumvents the difficulties arising from both the co-adaptation and the multi-agent credit assignment. Coordinated agents, for example, simultaneously learn their control choices and the other agent strategies assuming instantaneous and free explicit communications [START_REF] Guestrin | Coordinated Reinforcement Learning[END_REF], [START_REF] Kok | Sparse Cooperative Q-learning[END_REF]. While methods of this family inherit from standard single-agent techniques, they need to circumvent two significant drawbacks: the explosion in the state space size; and the centralization of all learning resources in a single place. Recently, team algorithms ranging from Q-learning to policy-search have been introduced for finite Dec-POMDPs, but with no guaranteed global optimality [START_REF] Wu | Monte-Carlo Expectation Maximization for Decentralized POMDPs[END_REF], [START_REF] Liu | Stick-breaking policy learning in Dec-POMDPs[END_REF][START_REF] Liu | Learning for Decentralized Control of Multiagent Systems in Large, Partially-Observable Stochastic Environments[END_REF], [START_REF] Kraemer | Multi-agent reinforcement learning as a rehearsal for decentralized planning[END_REF]. So, it seems one can either compute local optima with arbitrary bad performances or calculate optimal solutions but assuming noise-free, instantaneous and explicit communications.

A recent approach to optimally solving Dec-POMDPs suggests recasting them into occupancy-state MDPs (oMDPs) and then applying (PO)MDP solution methods [START_REF] Dibangoye | Optimally Solving Dec-POMDPs As Continuousstate MDPs[END_REF][START_REF] Dibangoye | Optimally solving Dec-POMDPs as Continuous-State MDPs: Theory and Algorithms[END_REF]Dibangoye et al. [ , 2015a[START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. In these oMDPs, the states called occupancy states are distributions over hidden states and joint histories of the original problem, and actions called decision rules are mappings from joint histories to controls [START_REF] Nayyar | Optimal Control Strategies in Delayed Sharing Information Structures[END_REF], [START_REF] Oliehoek | Sufficient Plan-Time Statistics for Decentralized POMDPs[END_REF], [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. This approach achieves scalability gains by exploiting the piece-wise linearity and convexity of the optimal value function. Since this methodology was successfully applied for planning in Dec-POMDPs, it is natural to wonder which benefits it could bring to the corresponding MARL problem. Unfortunately, a straightforward application of standard RL methods to oMDPs will face three severe limitations. First, occupancy states are unknown, and hence must be estimated. Second, they lie in a continuum making tabular RL methods inapplicable. Finally, the greedy maximization is computationally demanding in decentralized stochastic control problems [START_REF] Radner | Team Decision Problems[END_REF], [START_REF] Dibangoye | Point-based incremental pruning heuristic for solving finite-horizon DEC-POMDPs[END_REF], [START_REF] Kumar | Constraint-based dynamic programming for decentralized POMDPs with structured interactions[END_REF], [START_REF] Oliehoek | Heuristic search for identical payoff Bayesian games[END_REF].

This paper extends the methodology of Dibangoye et al. to MARL, focussing on the three major issues that limit its applicability. Our primary result is the proof that, by restricting attention to plans instead Inria of policies, a linear function over occupancy states and joint decision rules, which is simple to store and update, can capture the optimal performance for Dec-POMDPs. We further use plans instead of policies in a policy iteration algorithm, with the plan always being improved with respect to a linear function and a linear function always being driven toward the linear function for the plan. Under accurate estimation of the occupancy states, the resulting algorithm, called occupancy-state SARSA (oSARSA) [START_REF] Rummery | On-line Q-learning using connectionist systems[END_REF], is guaranteed to converge with probability one to a near-optimal plan for any finite Dec-POMDP. To extend its applicability to higher-dimensional domains, oSARSA replaces the greedy (or soft) maximization by a mixed-integer linear program for finite settings. Altogether, we obtain a MARL algorithm that can apply to finite Dec-POMDPs. Experiments show our approach can learn to act near-optimally in many finite domains from the literature.

We organize the remainder of this paper as follows. Section 2 extends a recent planning theory, starting with a formal definition of finite Dec-POMDPs. We proceed with the introduction of a framework for centralized MARL in Dec-POMDPs in Section 3. Also, we discuss our solutions to the three limitations mentioned above. We further present the resulting algorithm oSARSA along with convergence guarantees in Section 4. Finally, we conduct experiments in Section 5, demonstrating our approach can learn to act optimally in many finite domains from the literature.

2 Planning in Dec-POMDPs as oMDPs

Finite Dec-POMDPs

A finite Dec-POMDP is a tuple M . " pn, X, tU i u, tZ i u, p, r, , λ, b 0 q, where n denotes the number of agents involved in the decentralized stochastic control process; X is a finite set of hidden world states, denoted x or y; U i is a finite private control set of agent i P v1; nw, where U " U 1 ˆ¨¨¨ˆU n specifies the set of controls u " pu 1 , . . . , u n q; Z i is a finite private observation set of agent i, where Z " Z 1 ˆ¨¨¨ˆZ n specifies the set of observations z " pz 1 , . . . , z n q; p describes a transition probability kernel with conditional density p u,z px, yq; r is a reward model with immediate reward rpx, uq, we assume rewards are two-side bounded, i.e., for some c P R `, @x P X, u P U : |rpx, uq| ď c; is the planning horizon; γ P r0, 1s denotes the discount factor; and b 0 is the initial belief state with density b 0 px 0 q. We shall restrict attention to finite planning horizon ă 8 since an infinite planning horizon solution is within a small scalar ą 0 of a finite horizon solution where " rlog λ pp1 ´λq {cqs.

Because we are interested in MARL, we assume an incomplete knowledge about M , i.e., p and r are either unavailable or only through a generative model. Hence, the goal of solving M is to find a plan, i.e., a tuple of individual decision rules, one for each agent and time step: ρ . " pa 1 0: , . . . , a n 0: q. A tth individual decision rule a i t : O i t Þ Ñ PpU i q of agent i prescribes private controls based on the whole information available to the agent up to the tth time step, i.e., history of controls and observations o i t " pu i 0:t´1 , z .

" γ 1 r t `¨¨¨`γ r , where γ t denotes the time-step dependent weighting factors, often set to γ t " λ t for discounted problems or γ t " 1 for the average reward case. For any control interval t, joint plans a 0:t of interest are those that achieve the highest performance measure Jpa 0:t q

. " E a0:t tR 0 | b 0 u starting at b 0 , where E a0:t t¨u denotes the expectation with respect to the probability distribution over state-action pairs joint plan a 0:t induces, in particular Jpρq . " Jpa 0: ´1q for ρ . " a 0: ´1. One can show that, in Dec-POMDPs, there always exists a deterministic plan that is as good as any stochastic plan [see Puterman, 1994, Lemma 4.3.1]. Unfortunately, there is no direct way to apply the theory developed for Markov decision processes [START_REF] Bellman | Dynamic Programming[END_REF], [START_REF] Puterman | Markov Decision Processes, Discrete Stochastic Dynamic Programming[END_REF] to Dec-POMDPs, including: the Bellman optimality equation; or the policy improvement theorem. To overcome these limitations, we rely on a recent theory by Dibangoye et al. that recasts M into an MDP, thereby allowing knowledge transfer from the MDP setting to Dec-POMDPs.

Occupancy-State MDPs

To overcome the fact that agents can neither see the actual state of the system nor explicitly communicate their noisy observations with each other, [START_REF] Szer | MAA*: A Heuristic Search Algorithm for Solving Decentralized POMDPs[END_REF] and later on [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF] suggest formalizing M from the perspective of a centralized algorithm. A centralized algorithm acts on behalf of the agents by selecting a joint decision rule to be executed at each control interval based on all data available about the system, namely the information state. The information state at the end of control interval t, denoted ι t`1

.

" pb 0 , a 0:t q, is a sequence of joint decision rules the centralized algorithm selected starting at the initial belief state. Hence, the information state satisfies the following recursion: ι 0 . " pb 0 q and ι t`1

. " pι t , a t q for all control interval t, resulting in an ever-growing sequence. To generalize the value from one information state to another one, Dibangoye et al. introduced The process the occupancy states describe is known as the occupancy-state Markov decision process (oMDP), and denoted M 1 .

" pS, A, R, T, , γ, s 0 q. Similarly to POMDPs, it was proven that Dec-POMDPs can be recasted into MDPs, called oMDPs, and an optimal solution of the resulting oMDP is also an optimal solution for the original Dec-POMDP [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. M 1 is an -steps deterministic and continuous MDP with respect to M , where S . " Y tPv0; ´1w S t is the set of occupancy states up to control interval ´1; A . " Y tPv0; ´1w A t is the set of joint decision rules up to control interval ´1; R is the reward model; and T is the transition rule; s 0 is the initial occupancy state, which is essentially the initial belief in M ; γ and are as in M . It is worth noticing that there is no need to construct explicitly M 1 ; instead we use M (when available) as a generative model for the occupancy states T ps t , a t q and rewards Rps t , a t q, for all control intervals t.

To better understand why we use plans instead of policies and how they relate, consider the MDP case. The solution of any finite MDP called a policy π : S Þ Ñ A can be represented as a decision tree, where nodes are labelled with actions and arcs are labelled with states. Since an oMDP is also an MDP, its policies can also be represented as decision trees, except that actions are decision rules and states are occupancy states. In contrast to standard MDPs, oMDPs are deterministic. This means that only a single branch in the decision-tree representation-i.e., a sequence of actions-is necessary to act optimally in oMDPs. A single branch of a decision tree is called a plan. Hence policies are more general than plans, but in deterministic MDPs, both can be employed while achieving optimal performance (plans inducing an open-loop approach, and policies a closed-loop approach). We shall restrict attention to plans because they are more concise than policies. Below, we review a closed-loop approach based on the dynamic programming theory [START_REF] Bellman | Dynamic Programming[END_REF].

Inria

For any finite M , the Bellman equation is written as follows: for all occupancy state s t P S t , and some fixed policy π,

V π t ps t q . " Rps t , πps t qq `γ1 V π t`1 pT ps t , πps t qqq (1)
with boundary condition V π p¨q . " 0, describes the return of a particular occupancy state s t when taking decision rule a t " πps t q prescribed by π. The equation for an optimal policy π ˚is referred to as the Bellman optimality equation: for any control interval t, and occupancy state s t ,

V t ps t q . " max atPA Rps t , a t q `γ1 V t`1 pT ps t , a t qq (2)
with boundary condition V ˚ p¨q . " 0. Unfortunately, occupancy states lie in a continuum, which makes exact dynamic programming methods infeasible. Interestingly, when optimized exactly, the value function solution of (2) along with the boundary condition is always piece-wise linear and convex in the occupancystate space [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF].

Lemma 1. For any arbitrary M 1 , the solution V 0: of ( 2) is convex in the occupancy-state space. If we restrict attention to deterministic policies and finite M (and corresponding M 1 ), the solution of ( 2) is piece-wise linear and convex in the occupancy-state space. Hence, the optimal value at any occupancy state s t is as follows:

V t ps t q . " max αtPΓt xs t , α t y, (3)
where xs t , α t y is used to express the expectation of a linear function α t (also called α-vectorin the probability space defined by sample space X ˆO, the σ-algebra X ˆO and the probability distribution s t ; and Γ t is the set of all tth α-vectors.

Proof. Consider the αvector induced when agents follow plan a t: ´1 from control interval t onwards, denoted α a t: ´1 and given by α a t: ´1 px, oq . " EtR t |x t " x, o t " o, a t: ´1u. Hence, the optimal value starting at any occupancy state s t is given by taking the maximum over values of all possible plans from control interval t onwards: V t ps t q " max a t: ´1 xs t , α a t: ´1 y. In addition, the linearity of the expectation also implies that α a t: ´1 is linear in the occupancy-state space. The proof directly follows from [Rockafellar, 1970, Theorem 5.5].

Lemma 1 shows that for any arbitrary M and corresponding M 1 , the solution of (2), represented by sets Γ 0: , is convex in the occupancy-state space. Each α-vector defines the value function over a bounded region of the occupancy-state space. In addition, it is associated with a plan, defining the optimal plan for a bounded region of the occupancy-state space. Sets Γ 0: are iteratively improved by adding a new α-vector that dominates current ones over certain regions of the occupancy-state space. The α-vector to be added is computed using point-based Bellman backup operator H:

rHΓ t`1 sps t q "
arg max

α a t : aPAt,αt`1PΓt`1 xs t , α a t y,
where α a t px, oq .

" Etrpx, uq `γ1 α t`1 py, po, u, zqq|au, for each hidden state x P X, and joint history o P O. To keep the number of α-vectors manageable, one can prune those that are dominated over the entire occupancy-state space. All in all, the oMDP reformulation permits us to solve finite M by means of M 1 using near-optimal planning methods leveraging on the special structure of the optimal value function [START_REF] Shani | A survey of point-based POMDP solvers[END_REF]. This methodology results in the current state-of-the-art algorithm to optimally solving finite Dec-POMDPs [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. So it seems natural to wonder if the same methodology can also succeed when applied to the corresponding reinforcement-learning problem. In other words, how can a centralized algorithm learn to coordinate a team of agents with possibly contradicting perceptual information?

3 Learning in Dec-POMDPs as oMDPs Using the oMDP reformulation, a natural approach to achieve centralized RL for decentralized stochastic control suggests applying exact RL methods. In the Q-learning algorithm [START_REF] Watkins | Q-Learning[END_REF], for example, one would learn directly the Q-value function when following a fixed policy π: for any control interval t P v0; ´1w,

Q π t ps t , a t q . " Rps t , a t q `γ1 V π t`1 pT ps t , a t qq (4)
with boundary condition Q π p¨, ¨q " 0. The policy improvement theorem provides a procedure to change a sub-optimal policy π into an improved one π Howard [1960]: for any control interval t P v0; ´1w,

πps t q . " arg max atPAt Q π t ps t , a t q.
(5)

Unfortunately, this approach has three severe limitations. First, the occupancy states are unknown and must be estimated. Second, even if we assume a complete knowledge of the occupancy states, they lie in a continuum, which precludes exact RL methods to accurately predict α-vectors even in the limit of infinite time and data. Finally, the greedy maximization required to improve the value function proved to be NP-hard in finite settings [START_REF] Radner | Team Decision Problems[END_REF], [START_REF] Dibangoye | Point-based incremental pruning heuristic for solving finite-horizon DEC-POMDPs[END_REF], [START_REF] Kumar | Constraint-based dynamic programming for decentralized POMDPs with structured interactions[END_REF], [START_REF] Oliehoek | Heuristic search for identical payoff Bayesian games[END_REF].

Addressing Estimation Issues

Although mappings T and R in M 1 are unknown to either agents or a centralized algorithm, one can instead estimate on the fly both T ps 0 , a 0:t´1 q and RpT ps 0 , a 0:t´1 q, a t q for some fixed plan ρ . " a 0: ´1 through successive interactions of agents with the environment. To this end, we shall distinguish between two settings. The first one assumes a generative model is available during the centralized learning phase, e.g. a black box simulator; and the second does not. In both cases, we build on the concept of replay pool [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF], except that we extend it from stationary single-agent domains to non-stationary multi-agent domains.

If a generative model is available during the learning phase, then a Monte Carlo method can approximate T ps 0 , a 0:t´1 q and RpT ps 0 , a 0:t´1 q, a t q arbitrary closely. To this end, the generative model allows the agents to sample experiences generated from M . An -steps experience is a 4-tuple ξ . " px 0: ´1, u 0: ´1, r 0: ´1, z 1: q, where x 0: ´1 are sampled hidden states, u 0: ´1 are controls made, r 0: ´1 are reward signals drawn from the reward model, and z 1: are the resulting observations, drawn from the dynamics model. If we let D ρ .

" tξ ris u iPv1:Kw be the replay pool of K i.i.d random samples created through successive interactions with the generative model, then empirical occupancy state ŝt « T ps 0 , a 0:t´1 q and reward Rt « RpT ps 0 , a 0:t´1 q, a t q corresponding to the current D ρ are given by: for any control interval t P v0 : ´1w, ŝt px, oq

.

" 1 K ř K i"1 δ x px ris t q ¨δo pu ris 0:t , z ris
1:t q (6) and Rt . " 1 K ř K i"1 r ris t , (7)
where δ x p¨q and δ o p¨q denote the delta-Dirac mass located in hidden state and joint history pair, respectively. By the law of large numbers the sequence of averages of these estimates converges to their expected values, and the standard-deviation of its error falls as 1{ ? K [Sutton and Barto, 1998, chapter 5]. The error introduced by Monte Carlo when estimating T pŝ t´1 , a t´1 q instead of T ps 0 , a 0:t´1 q is upper bounded by 2 { ? K. The proof follows from the performance guarantee of the policy-search algorithm by [START_REF] Bagnell | Policy Search by Dynamic Programming[END_REF]. Hence, to ensure the learned value function is within ą 0 of the optimal one, one should set the replay-pool size to K " Θp4 When no generative model is available, the best we can do is to store samples agents collected during the learning phase into replay pools D ρ , one experience for each episode within the limit size of K. We maintain only the K recent experiences, and may discard1 hidden states since they are unnecessary for the updates of future replay pools and the performance measure. The rationale behind this approach is that it achieves the same performances as a Monte Carlo method for the task of approximating T ps 0 , a 0:t´1 q and RpT ps 0 , a 0:t´1 q, a t q given a fixed plan ρ . " a 0: ´1. In fact, if we let D ρ be a replay pool of K i.i.d. samples generated according to ρ, the empirical occupancy state ŝt « T ps 0 , a 0:t´1 q and reward Rt « RpT ps 0 , a 0:t´1 q, a t q corresponding to D ρ are given by ( 6) and ( 7), respectively. One can further show this approach preserves performance guarantees similar to those obtained when using a generative model.

Addressing Prediction Issues

The key issue with large spaces of occupancy states and decision rules is that of generalization, that is, how experiences with a limited subset of occupancy states and decision rules can produce a good approximation over a much larger space. Fortunately, a fundamental property of oMDPs is the convexity of the optimal value function over the occupancy-state space, see Lemma 1. Building on this property, we demonstrate a simple yet important preliminary result before stating the main result of this section.

Lemma 2. For any arbitrary M 1 (resp. M ), the optimal Q-value function is the upper envelope of sets Ω 0: of α-vectors over occupancy states and joint decision rules: for any control interval t, Q t ps t , a t q " max qtPΩ t xs τ d a τ , q t y, where q t P Ω t are appropriate α-vectors, and s τ d a τ denotes the Hadamard product2 .

Proof. We proceed by induction to prove this property. In the following we assume that all operations (e.g. integrals) are well-defined in the corresponding spaces. For control interval t " ´1, we only have to take into account the immediate reward and, thus, we have that Q ˚ ´1ps ´1, a ´1q " Rps ´1, a ´1q. Therefore, if we define the set Ω ˚ ´1 " tq ´1u, where q ´1px, o, uq .

" rpx, uq, the property holds at control interval t " ´1. We now assume the property holds for control interval τ `1 and we show that it also holds for control interval τ . Using (2) and (4), we have that, Q τ ps τ , a τ q " Rps τ , a τ q γ1 max aτ`1 Q τ `1pT ps τ , a τ q, a τ `1q, and by the induction hypothesis, let s τ `1 .

" T ps τ , a τ q:

Q τ `1ps τ `1, a τ `1q " max qPΩ τ `1 ÿ

x,o,u s τ px, oqa τ pu|oq ÿ y,z,u 1 p u,z px, yqa τ `1pu 1 |o, u, zqqpy, po, u, zq, u 1 q.

With the above,

Q τ ps τ , a τ q " max aPAτ`1,qPΩ τ `1 ÿ x,o,u s τ px, oqa τ pu|oqrrpx, uq `γ1 ÿ y,z,u 1
p u,z px, yqapu 1 |o, u, zqqpy, po, u, zq, u 1 qs.

At this point, we can define the bracketed quantity as q aτ`1 px, o, uq . " rpx, uq `γ1 ÿ y,z,u 1 p u,z px, yqa τ `1pu 1 |o, u, zqqpy, po, u, zq, u 1 q.

Note that α-vector q aτ`1 is independent of occupancy state s τ and decision rule a τ for which we are computing Q τ . With this, we have that Q τ ps τ , a τ q " max q a : aPAτ`1,qPΩ τ `1 xs τ d a τ , q a y and, thus the lemma holds.

Lemma 2 generalizes the convexity property demonstrated in Lemma 1 from optimal value functions over occupancy states to optimal value functions over occupancy states and decision rules. As a consequence, finite sets Ω 0: ´1 of α-vectors can produce solutions arbitrarily close to the optimal Q-value function Q 0: ´1. Though Q-value function Q 0: ´1 generalizes from a pair of occupancy state and decision rule to another one, storing and updating a convex hull is non trivial. Instead of learning the optimal Q-value function over all occupancy states and decision rules, we explore a simpler yet tractable alternative, which will prove sufficient to preserve ability to eventually find an optimal plan starting at initial occupancy state s 0 .

Theorem 1. For any arbitrary M 1 (resp. M ), the Q-value function Q ρ 0 : ´1 under an optimal plan ρ ˚. " a 0: ´1 starting at initial occupancy state s 0 is linear in occupancy states and decision rules: Q ρ t ps t , a t q " xs t d a t , q ρ t y where q ρ t .

" arg max qtPΩ t xT ps 0 , a 0:t´1 q d a t , q t y .

Proof. The proof derives directly from Lemma 2. First, notice that any arbitrary non-dominated joint plan ρ induces a sequence of α-vectors q ρ 0: ´1 stored in Ω 0: ´1, which proves the Q-value function under a fixed plan is linear over occupancy states and joint decision rules. In addition, each α-vector q ρ t P Ω t describes the expected returns from t P v0; ´1w onward, when agents follow non-dominated joint plan ρ. If we let ρ ˚be a greedy joint plan with respect to Q 0: ´1, then q ρ 0: ´1 is maximal along tT ps 0 , a 0:t´1 qu tPv0; ´1w .

Theorem 1 proves that the Q-function for a given optimal joint plan achieves performance at the initial occupancy state s 0 as good as the Q-value function for an optimal joint policy. Standard policy iteration algorithms search for an optimal joint policy, which requires a finite set of α-vectors to approximate V ˚/Q ˚, hence the resulting PWLC approximator is tight almost everywhere. Building upon Theorem 1, we search for an optimal ρ, which requires only a single α-vector to approximate V ρ {Q ρ , thus the resulting linear approximator is loose everywhere except in the neighborhood of a few points. The former approach may require less iterations before convergence to an optimal joint policy, but the computational cost of each iteration shall increase with the number of α-vectors maintained. The latter approach may require much more iterations, but all iteration shares the same computational cost.

Addressing Plan Improvement Issues

A fundamental theorem in many RL algorithms is the policy improvement theorem, which helps improving policies over time until convergence. This section introduces a procedure to improve a plan starting with a sub-optimal one. Suppose we have determined the value function V ρ 0: ´1 for any arbitrary ρ . " a 0: ´1. For some control interval t P v0; ´1w, we would like to know whether or not we should change decision rules a 0:t to choose ā0:t ‰ a 0:t . We know how good it is to follow the current plan from control interval t onward-that is V ρ t -but would it be better or worse to change to the new plan? One way to answer this question is to consider selecting ā0:t at control interval t and thereafter following decision rules a t`1: ´1 of the existing ρ. The value of the resulting joint plan is given by Jpā 0:t´1 q `γ1 V ρ t`1 pT ps 0 , ā0:t´1 qq. The key criterion is whether this quantity is greater or less than Jpρq. Next, we state the plan improvement theorem for oMDPs.

Theorem 2. Let ρ . " a 0: ´1 and ρ .

" ā0: ´1 be any pair of plans and J 0: be a sequence of α-vectors such that, for all t, J t px t , o t q . " Etα 0 r 0 `. . . `αt r t |b 0 , x t , o t , a 0:t´1 u. Let st . " T ps 0 , ā0:t´1 q and s t

. " T ps 0 , a 0:t´1 q be occupancy states at any control interval t P v0; ´1w under ρ and ρ, respectively. Then, xā 0:t ˚´1 , a t ˚: ´1y such that t ˚" arg max tPv0; ´1w xs t ´st , J t ´γ1 V ρ t y is as good as, or better than, ρ.

Inria

Proof. The proof follows from the difference between the performance measure of ρ . " a 0: ´1 and ρ .

" ā0: ´1. Let ς t pρ, ρq be the advantage of taking plan xā 0:t´1 , a t: ´1y instead of ρ: for any control interval t P v0; ´1w, ς t pρ, ρq " Jpā 0:t´1 q `γ1 V ρ t pT ps 0 , ā0:t´1 qq ´Jpρq " Jpā 0:t´1 q ´Jpa 0:t´1 q `γ1 pV ρ t ps t q ´V ρ t ps t qq " xs t ´st , J t ´γ1 V ρ t y.

If we let t ˚. " arg max t"0,1,..., ´1 ς t pρ, ρq, then plan xā 0:t ˚´1 , a t ˚: ´1y achieves the highest advantage among plan set txā 0:t´1 , a t: ´1yu tPv0; ´1w constructed based on ρ. If t ˚" 0, then xā 0:t ˚´1 , a t ˚: ´1y " ρ, and no improved plans were found from plan set generated from ρ. Otherwise, new xā 0:t ˚´1 , a t ˚: ´1y must be better than ρ.

Theorem 2 plays the same role in the plan space as does the policy improvement theorem in the policy space. More precisely, after sampling a plan, i.e., a sequence of decision rules, it tells us which of these decision rules will improve the current plan. More specifically, it shows how, given ρ . " a 0: ´1 and α-vector q ρ 0: ´1, we can easily evaluate a change in ρ at any control interval to a particular (possibly improved) plan. To ease exploration towards promising plans, we investigate the -greedy maximization (or soft-maximization). At each control interval t and occupancy state s t , it randomly selects ât with probability ; otherwise, it greedily selects ât w.r.t. the current Q-value function

Q ρ t : ât . " arg max at : a 1 t PA 1 t ,...,a n t PA n t Q ρ t ps t , a t q,
where ρ . " â0: ´1. Unfortunately, this operation proved to be NP-hard for finite M Radner [1962], [START_REF] Dibangoye | Point-based incremental pruning heuristic for solving finite-horizon DEC-POMDPs[END_REF], [START_REF] Kumar | Constraint-based dynamic programming for decentralized POMDPs with structured interactions[END_REF], [START_REF] Oliehoek | Heuristic search for identical payoff Bayesian games[END_REF]. Searching for the best decision rule requires enumerating all of them, which is not possible in large planning horizons. Instead, we present a mixed-integer linear programming (MILP) method, which successfully performs the greedy maximization for finite M . Though MILP is NP-hard in the worst case, the solution of its LP relaxation, which is in P, is often integral in our experiments. In other words, the solution of the LP relaxation is already a solution of the MILP. A similar observation was done before by MacDermed and Isbell. Mixed-Integer Linear Program 1 builds on MacDermed and Isbell [2013], which introduced an integer program for the greedy maximization in finite M . We also exploit the occupancy state estimation, in which ŝt replaces s t , and the current α-vector q ρ t . Mixed-Integer Linear Program 1 (For finite M ).

Maximize

at,a 1 t ,...,a n t ř x ř o ŝt px, oq ÿ u a t pu|oq ¨qρ t po, uq (8)
s.t.: ř u j a t pu j , u i |oq " a i t pu i |o i q, @i, u i , o (9) ř u a t pu|oq " 1, @o (10)
where ta t pu|oqu and ta i t pu i |o i qu are positive and boolean variables, respectively.

Mixed-Integer Linear program 1 optimizes positive variables ta t pu|oqu uPU,oPOt , one positive variable for each control-history pair. More precisely, each variable represents the probability a t pu|oq of control u being taken given that agents experienced joint history o. Constraints must be imposed on these variables to ensure they form proper probability distributions (10), and that they result from the product of independent probability distributions (9), one independent probability distribution for each agent. In order to make the description of the conditional independence,

a t pu|oq " a 1 t pu 1 |o 1 q ˆ¨¨¨ˆa n t pu n |o n q, (11)
Algorithm 1 The oSARSA Algorithm Initialize ḡ " ´8, ρ and q 0: ´1 arbitrary, and D ρ. while q 0: ´1 has not converged do Select -greedily ρ w.r.t. q 0: ´1 and D ρ.

Compose D ρ with N trajectories tξ rτ s u N τ "1 . Estimate pg, ςq from r ř ´1 t"0 Rt |D ρ , ŝ0 " s 0 s. If g ´ς ě ḡ then pρ, ḡ, D ρq " pρ, g `ς, D ρ q. Update α-vectors q 0: ´1 as described in ( 12). end while we use additional variables

ta i t pu i |o i qu iPv1;nw,u i PU i ,o i PO i t .
Marginalizing out both sides of ( 11) over all control-history pairs of all agents except agent i, denoted ´i, leads to (9). That is not sufficient to ensure conditional independence in general. If we further constrain ta i t pu i |o i qu to be boolean, then system of equations ( 9) implies (11). Given ( 9) and (10), agent variables ta i t pu i |o i qu u i PU i ,o i PO i t describe a proper probability distribution, so we omit corresponding constraints. Our greedy maximization approach is fundamentally different from previous ones, including the integer program by [START_REF] Macdermed | Point Based Value Iteration with Optimal Belief Compression for Dec-POMDPs[END_REF] and the constraint optimization program by [START_REF] Kumar | Constraint-based dynamic programming for decentralized POMDPs with structured interactions[END_REF], [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF]. First, while previous approaches made use of boolean variables, we use both positive and boolean variables instead. Next, prior approaches optimize a value function represented as a convex hull; we optimize an α-vector instead.

The oSARSA Algorithm

This section presents the oSARSA algorithm with tabular representations and function approximations (using either linear functions or deep neural networks) along with convergence guarantees. oSARSA algorithms are specializations of Policy Iteration, except that we use plans instead of policies. For the sake of conciseness, we describe a generic algorithm, which can fit to either tabular or approximate representations.

In Dec-POMDPs, the goal of oSARSA is to learn q 0: ´1, a sequence of α-vectors of an optimal plan ρ ˚.

In particular, we must estimate q t px, o, uq for the current plan ρ and for all reachable state x, joint history o, control u, and any control interval t. At the same time, the algorithm changes ρ towards improved plans according to the plan improvement theorem. The improved plans are constructed by exploring the occupancy-state space according to -greedy plans (see Section 3.3). To provide good estimations, we store all experiences in data set D ρ , from which we estimate the occupancy states and returns under ρ for any control interval (see Section 3.1). Upon estimating occupancy state ŝ and selecting joint decision rule a, we update parametrized α-vector q t with parameter θ t using q t`1 , D ρ and a To establish the convergence of oSARSA, we introduce the following assumptions.

Inria

Theorem 3. Consider assumptions: (1) The stepsizes tβ τ u τ "1,2,... satisfy Robbins and Monro's conditions;

(2) The occupancy states ŝ0: ´1 and immediate returns R0: ´1 are accurately estimated; and (3) Every pair of reachable occupancy state and joint decision rule is visited infinitely often. Under these assumptions, the sequence q rτ s 0: ´1 generated by oSARSA converges with probability 1 to q 0: ´1. Proof. Under these assumptions, we define H ρ that maps a sequence of α-vectors q 0: ´1 to a new sequence of α-vectors H ρ q 0: ´1 according to the formula: for all hidden state x, joint history o and control u, at control interval t, pH ρ q 0: ´1qpx, o, uq " rpx, uq `γ1 12), we converge with probability 1 to q ρ 0: ´1. It now remains to be verified that the plan improvement step of the oSARSA algorithm changes the current plan for an improved one. Initially, ḡ is arbitrarily bad, so any new plan is an improved one. Then, ḡ " Jpρq for the current best plan ρ since occupancy state and return are accurately estimated. Hence, when ever g ě ḡ, we know that the new plan ρ yields a performance measure Jpρq superior to Jpρq, thus ρ improves ρ. We conclude the proof noticing that in finite M , the number of deterministic plans is finite. As a consequence, by visiting infinitely often every pair of occupancy state and decision rule we are guaranteed to visit all deterministic plans, hence an optimal one.

It is now important to observe that we meet assumption (2) in Theorem 3 only when M is available. Otherwise, we rely on confidence bounds rg ´ς, g `ςs, e.g. Hoeffding's inequality, on estimate g « Jpρq.

In particular, we use lower-bounds g ´ς on sample means instead of the sample means g themselves, to limit situations where g is overestimated. Small data sets often lead to suboptimal solutions, but as the number of experiences in data set D ρ increases, sample means and corresponding lower bounds get close to the mean, i.e., ς tends to 0. We require an accurate estimation of a plan's performance to know for sure its performance is above that of any other plans we may encounter. However, an estimation of a sample plan's performance can be refined over time until it becomes accurate. For example, if the algorithm samples a promising plan-i.e., its confidence bounds suggest it might achieve a better performance than that of the current plan-the algorithm can progressively refine it until it becomes accurate. Of course, having a good initial estimate can significantly speed up the convergence. In the extreme case, the initial estimate is the true value. It is also worth noticing that the memory complexity of the oSARSA algorithms is linear with the size of an α-vector, i.e., Θp|D ρ |q; and its time complexity is linear with the episodes.

Experiments

We ran the oSARSA algorithm on a Mac OSX machine with 3.8GHz Core i5 and 8GB of available RAM. We solved the MILPs using ILOG CPLEX Optimization Studio. We define features to use sequences of K last joint observations instead of joint histories, hence the dimension of the parameter vector θ is |X|p|U ||Z|q K for finite M .

We evaluate our algorithm on multiple 2-agent benchmarks from the literature all available at masplan.org: Mabc, Recycling, Gridsmall, Grid3x3corners, Boxpushing, and Tiger. These are the largest and the most challenging benchmarks from the Dec-POMDP literature. For each of them, we compare our algorithm to the state-of-the-art algorithms based on either a complete or a generative model: FB-HSVI [START_REF] Dibangoye | Optimally Solving Dec-POMDPs as Continuous-State MDPs[END_REF], RLaR [START_REF] Kraemer | Multi-agent reinforcement learning as a rehearsal for decentralized planning[END_REF], and MCEM [START_REF] Wu | Monte-Carlo Expectation Maximization for Decentralized POMDPs[END_REF]. We also reported results of the state-of-the-art model-free solver: (distributed) REINFORCE Peshkin et al. [2000]. For REINFORCE and oSARSA, we used hyper-parameters and β ranging from 1 to 10 ´3 with a decaying factor of 10 4 , sample size |D| " 10 4 . We use maximum episodes and time limit 10 5 and 5 hours, respectively, as our stopping criteria.

Surprisingly, REINFORCE performs very well on domains that consist of weakly coupled agents, see Figure 1. However, for domains with strongly coupled agents, e.g., Tiger or BoxPushing, it often gets stuck at some local optima. In contrast, oSARSA converges to near-optimal solutions when enough resources are available over all domains, see Figure 1 and Table 1. Regarding the most challenging benchmarks, which require more resources, oSARSA stops before the convergence to a near-optimal solution; yet, it often outperforms the other RL algorithms. RLaR can achieve near-optimal result for small domains and short planning horizon ( ď 5q, assuming there exists a unique optimal plan. As for MCEM, it can solve infinite horizon problems, but similarly to REINFORCE may get stuck in local optima; this is essentially as they both use a form of gradient descent in a parametrized policy space. Table 1: Comparing V ρ ps 0 q of all solvers when available, where the "-" sign mean "out of memory" and/or "out of time".

Discussion

This paper extends a recent but growing (deep) MARL paradigm [START_REF] Szer | MAA*: A Heuristic Search Algorithm for Solving Decentralized POMDPs[END_REF], Dibangoye et al.

[2016], [START_REF] Kraemer | Multi-agent reinforcement learning as a rehearsal for decentralized planning[END_REF], [START_REF] Mordatch | Emergence of Grounded Compositional Language in Multi-Agent Populations[END_REF], [START_REF] Foerster | Counterfactual Multi-Agent Policy Gradients[END_REF], namely RL for decentralized control, from model-based to model-free settings. This paradigm allows a centralized algorithm to learn on behalf of all agents how to select an optimal joint decision rule to be executed at each control interval based on all data available about the system during a learning phase, while still preserving ability for each agent to act based solely on its private histories at the execution phase. In particular, we introduced tabular and approximate oSARSA algorithms, which demonstrated promising results often outperforming state-of-the-art MARL approaches for Dec-POMDPs. To do so, oSARSA learns a value function that maps pairs of occupancy state and joint decision rule to reals. To ease the generalization in such high-dimensional continuous spaces, we restrict attention to plans rather than policies, which in turn restricts value functions of interest to linear functions. To speed up the greedy maximization, we used a MILP for finite settings-we shall use a gradient approach instead of a MILP for continuous settings in Inria Figure 1: Comparing V ρ ps 0 q of solvers with " 8 and γ " 0.9, where x-axis denotes the number of episodes during training.

future works. Finally, we present a proof of optimality for a MARL algorithm when the estimation error is neglected. We shall investigate an approach to relax this somewhat restrictive assumption, perhaps within the probably approximately correct learning framework.

The RL for decentralized control paradigm is significantly different from the standard RL paradigm, in which agents have the same amount of information during both the learning and the execution phases. Another major difference lies in the fact that learned value functions in standard (deep) RL algorithms are mapping from histories (or states) to reals. In contrast, oSARSA learns a value function that maps occupancy-state/decision-rule pairs to reals-spaces of occupancy states and joint decision rules are multiple orders of magnitude larger than history or state spaces. As a consequence, standard (MA)RL methods, e.g. REINFORCE and MCEM, may converge towards a local optimum faster than oSARSA, but the latter often converges towards a near-optimal solution. oSARSA uses occupancy states instead of joint histories mainly because occupancy states are (so far minimal) sufficient statistics for optimal decision-making in Dec-POMDPs-using joint histories instead of occupancy states may lead to suboptimal solutions except in quite restrictive settings. For example, RLaR learns value functions mapping history/action pairs to reals, but convergence towards an optimal solution is guaranteed only for domains that admit a unique optimal joint plan-which essentially restricts to POMDPs [START_REF] Kraemer | Multi-agent reinforcement learning as a rehearsal for decentralized planning[END_REF].

  the concept of occupancy states. The occupancy state at control interval t, denoted s t . " Ppx t , o t |ι t q, is a distribution over hidden states and joint histories conditional on information state ι t at control interval t. Interestingly, the occupancy state has many important properties. First, it is a sufficient statistic of the information state when estimating the (current and future) reward to be gained by executing a joint decision rule:In addition, it describes a deterministic and fully observable Markov decision process, where the next occupancy state depends only on the current occupancy state and next joint decision rule, for all y P X, o P

	Rps t , a t q	. "	ÿ	ÿ	s t px, oq	ÿ	a t pu|oq ¨rpx, uq.
			x	o		u

O, u P U, z P Z: T ps t , a t q . " s t`1 s t`1 py, po, u, zqq . " a t pu|oq ÿ x s t px, oq ¨pu,z px, yq.

  where β τ is a step size, and quantity ∇q t px, o, uq denotes the gradient of q t at px, o, uq w.r.t. some parameter θ t . Using tabular representations (e.g., finite/small M ), θ t " q t and thus ∇q t px, o, uq is a unit vector e x,o,u whose value at px, o, uq is one and zero otherwise. Using linear function approximations (e.g., continuous/large M ), q t px, o, uq . " φ t px, o, uq J θ t , where ∇q t px, o, uq " φ t px, o, uq is the feature vector at px, o, uq. Algorithm 1 shows the pseudocode of oSARSA.

						t`1 by means of temporal
	difference learning: for all px, o, uq,
	θ	rτ `1s t	. " θ	rτ s t `βτ E ŝ,a,D,at`1 tδ t ∇q t px, o, uqu rτ s	(12)
		δ t " r `γ1 q t`1 py, o 1 , u 1 q rτ s	´qrτs t px, o, uq,

  Etv t`1 py, o ' pu, zqqu, where v t px, oq . " q t px, o, ρpoqq and ρpoq is the control prescribed by ρ at joint history o. Then, the plan evaluation step of the oSARSA algorithm is of the form

	q t rτ `1s	px, o, uq " p1 ´βt qq	rτ s t px, o, uq `βt κ t px, o, uq, rτ s
	κ rτ s t px, o, uq " pH ρ q 0: ´1qpx, o, uq `wt px, o, uq, rτ s
	where w t px, o, uq " rpx, uq `γ1 v	rτ s t`1 py, o ' pu, zqq ´pH ρ q	rτ s 0: ´1qpx, o, uq is a zero mean noise term.
	Using this temporal-difference update-rule, see (

RR n°9179

Dibangoye & Buffet

2 q.Inria

Note that one should keep hidden states when available since they often speed up the convergence.

@px, o, uq : rsτ d aτ spx, o, uq . " sτ px, oq ¨aτ pu|oq.RR n°9179