Mahmoud Ali
email: mahmoud.ali@univ-orleans.fr

Mohamed Rady

Mohamed A A Attia

Emad M M Ewais

Consistent coupled optical and thermal analysis of volumetric solar receivers with honeycomb absorbers

Keywords: Concentrated solar power, volumetric solar receiver, absorptivity, porous ceramics, Monte Carlo ray tracing, CFD

In concentrating solar power plants with central towers, successful design of volumetric solar receivers requires proper understanding of the interaction between optical, heat transfer, and fluid flow phenomena occurring at the microscopic scale of receiver structure material and their effect on the overall solar to thermal efficiency. In the present article, coupled, 3D, optical, heat transfer and fluid flow numerical models have been developed for the analysis and design of honeycomb volumetric receiver modules. The optical model considers the absorptivity and micro dimensions of honeycomb absorber structure and employs a Monte Carlo ray tracing technique to calculate and analyze the absorbed solar heat flux distribution. This, in turn, is employed as a volumetric heat source term at the solid surface for consistent heat transfer and fluid flow modeling using a realistic solution domain and proper boundary conditions. The validated models have been employed to investigate the effects of different types of absorber materials, material absorptivity, and air flow rate on the performance of the solar receiver. It has been shown that positive volumetric effect and high solar-to-thermal efficiency can be obtained by controlling the absorbed radiation heat flux distribution within the honeycomb receiver using surface coating of the absorber material.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Introduction

Concentrating Solar Power (CSP) is a promising and sustainable technology for electricity generation, chemical fuels production, solar heating and cooling, and water desalination [START_REF] Behar | A review of studies on central receiver solar thermal power plants[END_REF][START_REF] Blanco | Review of feasible solar energy applications to water processes[END_REF][START_REF] Gordon | High-efficiency solar cooling[END_REF][START_REF] Nakamura | Hydrogen production from water utilizing solar heat at high temperatures[END_REF].

Many types of CSP technologies have been utilized for electricity production such as solar tower power plants (STPPs), parabolic trough, parabolic dish, and linear Fresnel solar collectors [START_REF] Barlev | Innovation in concentrated solar power[END_REF].

Among the several types of CSP technologies, STPP with open volumetric solar receiver is one of the most promising technologies [START_REF] Ho | Review of high-temperature central receiver designs for concentrating solar power[END_REF].

STPP is based on the concept of reflecting and concentrating sunlight, using sun-tracking mirrors (heliostats), on the receiver at the top of a central tower. The volumetric solar receiver (VSR) absorbs the solar radiation and converts it to high-temperature heat. Then, a heat transfer fluid (HTF) fluid flows through the receiver and removes the high-temperature heat by forced convection mechanism [START_REF] Avila-Marin | Volumetric receivers in solar thermal power plants with central receiver system technology: a review[END_REF]. This heat is finally used to generate electricity using steam or gas turbine power plants [START_REF] Kribus | The promise and challenge of solar volumetric absorbers[END_REF]. Air is commonly used as a HTF due to its availability and stability.

In a typical solar tower power plant, VSR is a key element that determines the solar-tothermal efficiency of the plant. In order to reduce the thermal losses from the VSR to the environment, the outlet temperature of the HTF should be higher than the absorber frontal surface temperature indicating what is referred to in the literature as positive volumetric effect However, all VSR tested to date show high negative values of volumetric effect [START_REF] Kribus | The promise and challenge of solar volumetric absorbers[END_REF][START_REF] Romero | An update on solar central receiver systems, projects, and technologies[END_REF]. Different volumetric solar receivers' designs have been investigated numerically and experimentally. They differ by the structure and the material of the receiver. Both ceramic as well as metallic materials have been tested [START_REF] Fricker | Proposal for a novel type of solar gas receiver[END_REF][START_REF] Heinrich | Technology Program Solar Air Receiver-2.5 MWt System Test on Volumetric Air Receiver Technology[END_REF][START_REF] Hellmuth | Performance of a wire mesh solar volumetric air receiver[END_REF][START_REF] Pitz-Paal | Evaluation of the CATREC II Receiver Test[END_REF][START_REF] Hellmuth | Modeling and optimum design of a wire mesh solar volumetric air receiver[END_REF][START_REF] Tellez | Design of" SIREC-1" Wire Mesh Open Volumetric Solar Receiver Prototype[END_REF][START_REF] Chavez | Testing of a porous ceramic absorber for a volumetric air receiver[END_REF][START_REF] Pitz-Paal | A new concept of a selective solar receiver for high temperature applications[END_REF][START_REF] Reale | Final design of a multi cavity volumetric solar receiver[END_REF][START_REF] Hoffschmidt | The development strategy of the HiTREC volumetric receiver technology up scaling from 200 kW h via 3 MW h up to 10 MWel[END_REF][START_REF] Hoffschmidt | Performance evaluation of the 200-kWth HiTRec-II open volumetric air receiver[END_REF][START_REF] Hoffschmidt | Test results of a 3 MW solar open volumetric receiver[END_REF][START_REF] Tellez | Thermal performance of SolAir 3000 kWth ceramic volumetric solar receiver[END_REF]. Different structures include extruded honeycombs with parallel channels [START_REF] Hoffschmidt | The development strategy of the HiTREC volumetric receiver technology up scaling from 200 kW h via 3 MW h up to 10 MWel[END_REF][START_REF] Hoffschmidt | Performance evaluation of the 200-kWth HiTRec-II open volumetric air receiver[END_REF], wire meshes [START_REF] Hellmuth | Modeling and optimum design of a wire mesh solar volumetric air receiver[END_REF][START_REF] Tellez | Design of" SIREC-1" Wire Mesh Open Volumetric Solar Receiver Prototype[END_REF], ceramic foams [START_REF] Wu | Coupled radiation and flow modeling in ceramic foam volumetric solar air receivers[END_REF][START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF][START_REF] Zaversky | Numerical and experimental evaluation and optimization of ceramic foam as solar absorber-Single-layer vs multi-layer configurations[END_REF][START_REF] Wang | Thermal performance analysis of porous media receiver with concentrated solar irradiation[END_REF][START_REF] Alaa | Optical study of using ceramic foams for volumetric solar receivers[END_REF], and printed structures [START_REF] Capuano | Optimized volumetric solar receiver: Thermal performance prediction and experimental validation[END_REF]. Silicon carbide (SiC) honeycomb absorber structure has been used in many STPPs such as SolAir3000 power plant [START_REF] Tellez | Thermal performance of SolAir 3000 kWth ceramic volumetric solar receiver[END_REF].

Literature survey

Several numerical and experimental studies have been carried out to investigate the optical and thermal performance of VSRs using ceramic honeycombs as an absorber. Figure 1 shows the volumetric solar receiver geometry and the components of a VSR module with honeycomb absorber. Different scales of interaction between optical and thermal phenomena can be envisaged from this figure. Previous numerical studies reported in the literature employed a

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

4 physical model at the channel scale and focused on analyzing the influence of channel dimensions, such as channel width and wall thickness [START_REF] Lee | Consistent heat transfer analysis for performance evaluation of multichannel solar absorbers[END_REF][START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF][START_REF] Sanchez | Design Parameters Influence on Flux Distribution Through Prismatic Channels of Volumetric Absorbers[END_REF][START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF] and channel cross-section shape [START_REF] Yilbas | A Solar Volumetric Receiver: Influence of Absorbing Cells Configuration on Device Thermal Performance[END_REF] on the performance of the receiver. The important findings of these studies are summarized below.

A B C D

Fig. 1. Central receiver solar plant with honeycomb volumetric absorber: (A) Plant scale [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF], (B) Receiver front area [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF], (C) Single Module absorber and cup, (D) Channel dimensions Fend et al. [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF] carried out numerical and experimental studies to analyze the influence of channel dimensions on the thermal performance of the receiver. Two SiC honeycomb samples with different geometric parameters have been investigated. For sample 1: d = 2.18 mm & t = 0.55mm and for sample 2: d = 1.43 mm & t = 0.25 mm. The radiative heat flux distribution on the absorber surface has been estimated using one dimensional Beer Lambert Law (BLL) and has been applied as a heat source in a CFD model using COMSOL software. Their results indicated that the sample with smaller d and t (sample 2) shows better performance in terms of air outlet temperature and thermal efficiency due to the increased specific surface area.

However, in both cases, the air outlet temperature was less than the frontal surface temperature of the solid absorber indicating negative values of volumetric effect.

Lee et al. [START_REF] Lee | Consistent heat transfer analysis for performance evaluation of multichannel solar absorbers[END_REF] used the 3D Monte Carlo ray tracing method (MCRT) and a 1-D CFD model to investigate the influence of channel width (0.5 -1.5 mm), wall thickness (0.15 -0.25 mm), and absorptivity (α = 0.8, 1) on the performance of the receiver. Their results indicate that increasing the absorptivity leads to a decrease in the reflection losses from the frontal surface of the receiver and enhance the receiver performance. In agreement with Fend et al. [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF], they found that honeycomb channels with smaller channel width showed favorable performance in terms of HTF outlet temperature and thermal efficiency. However, Sanchez et al. [START_REF] Sanchez | Design Parameters Influence on Flux Distribution Through Prismatic Channels of Volumetric Absorbers[END_REF] and Cagnoli et al. [START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF], using ray tracing code and CFD model, have shown that favorable honeycomb receiver performance has been obtained by increasing the channel width.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

It should be noted that the previous studies differ in using either simple one-dimensional Beer Lambert law [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF][START_REF] Yilbas | A Solar Volumetric Receiver: Influence of Absorbing Cells Configuration on Device Thermal Performance[END_REF][START_REF] Smirnova | Determination of critical thermal loads in ceramic high concentration solar receivers[END_REF][START_REF] Capuano | Numerical models of advanced ceramic absorbers for volumetric solar receivers[END_REF] or 3D Monte Carlo ray tracing [START_REF] Lee | Consistent heat transfer analysis for performance evaluation of multichannel solar absorbers[END_REF][START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF][START_REF] Gomez-Garcia | Numerical analysis of radiation propagation in a multi-layer volumetric solar absorber composed of a stack of square grids[END_REF] for performing the optical analysis of the receiver and estimating the radiative heat flux distribution on the absorber surface at the scale of single channel dimensions. The later represents the main heat source that may significantly affect the thermal performance of the receiver. Gomez-Garcia et al. [START_REF] Gomez-Garcia | Numerical analysis of radiation propagation in a multi-layer volumetric solar absorber composed of a stack of square grids[END_REF] show that BLL can be used to describe the radiative heat flux distribution within the honeycomb structure near the entrance of the channel only for depths below two pitches.

Moreover, BLL neglects the direction effects of the incoming solar rays as well as rays scattering inside the honeycomb channel. Therefore, using BLL may lead to improper characterization of the radiative flux distribution and overestimation of the HTF outlet temperature. Further analysis to clarify these points and in-depth analysis for the proper choice of the solution domain for consistent numerical modeling of the coupled optical and thermal phenomena inside the receiver represent one of the main objectives of the present article.

Previous studies focused on using a material with high absorptivity (SiC) to reduce the reflection losses from the frontal surface of the receiver and increase the optical efficiency [START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF].

Numerical studies carried out by Yilbas and Shuja [START_REF] Yilbas | A Solar Volumetric Receiver: Influence of Absorbing Cells Configuration on Device Thermal Performance[END_REF] and Kasaeian et al. [START_REF] Kasaeian | Performance comparison between the geometry models of multi-channel absorbers in solar volumetric receivers[END_REF] indicate that cell configuration with triangular channel cross-section results in higher air outlet temperature and higher thermal efficiency and then followed by square, rectangular, hexagonal, and circular cross-sections. However, negative values of volumetric effect have been observed for all shapes of channel cross-sections evaluated. Gomez-Garcia et al. [START_REF] Gomez-Garcia | Numerical analysis of radiation propagation in a multi-layer volumetric solar absorber composed of a stack of square grids[END_REF] performed a numerical optical study, using MCRT method, and show that an increase in the penetration depth of concentrated solar radiation can be obtained by using a novel solar receiver composed of a stack of thick square grid.

In order to achieve the required positive volumetric effect and high solar-to-thermalefficiency, radiation propagation within the honeycomb receiver should be maximized [START_REF] Kribus | The promise and challenge of solar volumetric absorbers[END_REF]. Both the architecture and the high absorptivity of SiC honeycombs limits radiation propagation inside the honeycomb channel. As a result, the maximum penetration depth of the absorbed heat flux is usually less than 3.3 times the channel width [START_REF] Lee | Consistent heat transfer analysis for performance evaluation of multichannel solar absorbers[END_REF]. This leads to higher temperature at the frontal surface of the receiver, higher emission losses, lower solar-to-thermal efficiency, and lower air outlet temperature.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6

Outline and Contribution of Present Work

The major contributions of the present article and significant differences with previous reported work can be summarized as follows:

• Consistent numerical modeling of the coupled optical and thermal phenomena inside the volumetric receiver including proper solution domain, boundary conditions, optical modelling, radiation effects, and absorbed solar radiation distribution.

• Demonstration of the basic concepts to obtain positive values of volumetric effect and high solar-to-thermal efficiency for volumetric receivers. The essential requirement is shown to increase radiation propagation within the honeycomb receiver. However, both the architecture and the high absorptivity of SiC honeycombs have been shown to limit radiation propagation inside the honeycomb channel.

• An efficient method to maximize the penetration depth within the honeycomb absorber is introduced by judicious control of the absorptivity of the honeycomb channel along its length. Of course, the solar absorptivity of the frontal surface should be as maximum as possible to reduce the reflection losses. However; the absorptivity of the internal walls can be controlled to achieve the required radiation flux distribution and obtain favorable temperature distributions. According to the author's knowledge, this method has not been yet studied in previous research work.

• A novel honeycomb receiver made of Alumina has been introduced by using judicious surface coating of the absorber material. The coating has been designed with the objective of reducing the reflection losses from the frontal surface and increasing the penetration depth of the absorbed solar heat flux inside the honeycomb receiver channel. The new introduced coated Alumina honeycomb absorbers show a favorable receiver performance in terms of air outlet temperature and solar-tothermal efficiency compared to conventional SiC honeycomb absorbers.

• The proposed method can be applied to improve the performance of honeycomb volumetric receivers.

The present article is organized as follows: section 4 describes the physical model for the optical and CFD analysis. The heat transfer and Fluid Flow Modeling are illustrated in section 5.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

7

In section 6, the results of coupled optical and thermal models are presented and discussed.

Conclusions and key findings of the present study are presented in section 8.

Physical Models for Optical and CFD Analysis

Actually, the physical model for solar tower power plant consists of sun-tracking mirrors (heliostats) reflecting and concentrating sunlight on the receiver. However, in laboratory testing, solar furnace facilities are used to investigate the thermal and optical performance of VSRs due to high concentration ratios and high temperatures that can be achieved by these systems [START_REF] Levêque | High-flux optical systems for solar thermochemistry[END_REF]. As shown in Fig. 2, the solar furnace used in the present study is similar to the facility used in [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF] and composed of a reflector, a parabolic concentrator, and a VSR target similar to the one used in SolAir 3000. The total surface area of the reflector and the concentrator are 27 m 2 and 13.5 m 2 , respectively. The focal distance of the concentrator is 3.8 m. The solar furnace system redirects and focuses the incoming sunlight onto the receiver module surface. The total depth of the absorber (Z) is 50 mm with channels of square crosssection. The honeycomb absorber receives the radiative flux and converts it to hightemperature heat. Then, air flows through the receiver and removes the thermal energy by forced convection. The solar furnace optical model for the system shown in Fig. 2 has been developed using an experimentally validated object-oriented Monte Carlo based ray tracing software (Tonatiuh © ) [START_REF] Blanco | The Tonatiuh Software Development Project: An open source approach to the simulation of solar concentrating systems[END_REF][START_REF] Blanco | Experimental validation of Tonatiuh using the Plataforma Solar de Almería secondary concentrator test campaign data[END_REF][START_REF] Mutuberria | Results of the empirical validation of Tonatiuh at Mini-Pegase CNRS-PROMES facility[END_REF]. In a typical Tonatiuh ray-tracing simulation, the solar power is discretized into

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 8
photons. Each photon has the same amount of energy that represents the power per photon ( PPP ), which is calculated using:

) 1 ( . Sol t P PPP N =

Where

. Sol P is the solar power and t N is the total number of photons. The sun shape has been considered as pillbox and the value of direct normal irradiation (DNI) used in the present analysis is 1000 W/m 2 . Since the surface area of the receiver is smaller than the solar image reflected from the concentrator, it can be assumed that the frontal surface of the receiver (A) has a constant flux distribution. The area-averaged absorbed solar heat flux "

. Sol q that will be used as a source term for heat transfer analysis has been computed using:

) 2 ( " . . abs Sol PPP q A × Ν =

Where

. abs Ν is the number of absorbed photons. A sensitivity analysis has shown that a minimum number of 10 8 photons is required to reach an independent value of absorbed heat flux at the frontal area of the honeycomb receiver.

The ray-tracing algorithm is used to trace and collect the absorption and scattering data of the simulated photons through the system, starting from the light source (the sun) going towards the reflector, the concentrator and finally reaching the receiver. Ray tracing is also used for tracking the paths of the simulated photons through the receiver until they exit it or completely absorbed. This allows to determine the absorbed heat flux at each surface of the honeycomb channels (upper, lower, right, and left surfaces). Complete description of the optical model and analysis of results have been reported by the present authors elsewhere [START_REF] Alaa | Optical study of using ceramic foams for volumetric solar receivers[END_REF].

Tonatiuh has been designed to simulate either specular or diffuse surface behavior. The reflection mechanism of a real ceramic material depends on the manufacturing process. It is neither purely specular nor purely diffuse, but it is a mixture of the two. The effects of specular and diffuse reflectivity on the accuracy of the optical model have been addressed before by Cagnoli et al. [START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF]. The results show that the assumption of specular reflection is accurate enough to reproduce the actual receiver behavior. In the present study, specular directional reflectivity has been implemented in Tonatiuh.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 9

The optical analysis showed that the absorptivity of the material plays a dominant role in determining the radiative flux distribution on the honeycomb receiver surfaces. Also, the radiative flux distribution of honeycomb absorbers is highly non-uniform with prominent peaks and surface dependent. The latter finding is very important for consistent coupling of optical and CFD models. This variation of solar heat flux distribution should be properly included in the choice of physical model for CFD analysis. In the previous numerical studies [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF][START_REF] Capuano | Numerical models of advanced ceramic absorbers for volumetric solar receivers[END_REF], only onequarter of the honeycomb structure cells is used as a solution domain. This choice is valid only when the absorbed heat flux distributions on all honeycomb structure surfaces are identical.

Therefore, in this study, an appropriate solution domain has been used as depicted in Fig. 3.

The proposed solution domain allows the authors to set the absorbed heat flux of each surface as a standalone one in consistency with the optical analysis. Further details on the governing equations of CFD analysis, boundary conditions implementations and mesh generation are outlined in Section 3.

Heat transfer and Fluid Flow Modeling

Governing Equations

The complete model of volumetric solar receiver couples the fluid dynamics and the conjugate heat transfer physics. In the present study, ANSYS © software has been used for numerical modeling of the CFD phenomena inside a single channel of the numerical solution domain shown in Fig. 3. The solution domain comprising both fluid and solid domains has been discretized using ANSYS © meshing software. The grid consists of tetrahedron elements as shown in Fig. 4. In comparison with hexahedron elements, tetrahedron elements are more

M A N U S C R I P T A C C E P T E D

ACCEPTED MANUSCRIPT 10 efficient, require less computation demand, but are more complex to generate. In addition, five inflation layers placed in the fluid domain at the interface with the solid domain in order to solve the boundary layer where large velocity and temperature gradients usually exist. As a measure of the mesh quality, the Skewness control has been adopted. The Skewness is a measure of the relative distortion of an element compared to its ideal shape and is scaled from zero (Excellent) to one (unacceptable). For all meshes created in the present work, the Skewness value has never exceeded 0.7 indicating a good quality mesh [START_REF] Fluent | ANSYS fluent theory guide 15[END_REF]. The air flow inside the absorber channel is considered to be steady, incompressible and laminar. Under these assumptions, the 3D governing continuity, momentum, and energy equations are written as follows [START_REF] Pletcher | Computational fluid mechanics and heat transfer[END_REF]:

) 3 ( ( ) 0 i i u x ρ ∂ = ∂ ) 4 ( ( ) i j i j i j j u dp u u x x x x ρ µ   ∂ ∂ ∂ = - +     ∂ ∂ ∂ ∂   (5)
( )

f i p f f i i i T u c T K x x x ρ ∂   ∂ ∂ =   ∂ ∂ ∂  
The energy equation for the solid phase is written as [START_REF] Pletcher | Computational fluid mechanics and heat transfer[END_REF]:

(6)

0 s s v i i T K S x x   ∂ ∂ + =   ∂ ∂  
In the above equations, the subscript s refers to the solid and f refers to the fluid. The radiative heat transfer between the inner walls of the channel and between the inner walls and the ambient has been taken into account using the surface to surface radiation model (S2S) [START_REF] Fluent | ANSYS fluent theory guide 15[END_REF]. In the S2S radiation model, the surfaces are assumed to be gray and diffuse (absorptivity is equal to the emissivity and is independent of the wavelength). Each internal surface of the channel is divided into a number of discrete surfaces ( k = 1: N). Then view factors, which are a function of orientation, size and distance of the k surface, are computed [START_REF] Fluent | ANSYS fluent theory guide 15[END_REF]. The radiative heat flux leaving from one surface k to another surface is calculated according to:

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT ( ) 4 , , 1 out k k k k in k q T q = ε σ + -ε (8)
Where:

, out k q Radiative heat flux leaving surface k (W/m 2 ) k ε Absorptivity of surface k σ Stefan-Boltzmann constant (W/m 2 K 4 ) k T Temperature of surface k (K) , in k q
Energy flux incident on the surface k (W/m 2 )

In order to determine the overall performance of volumetric solar receiver, the solar to thermal efficiency (

= × × - = = & η η η (9)
Where . opt η , . Th η , air P , air in T and air out T are the optical efficiency, thermal efficiency, rate of heat transfer to air, and air inlet and outlet temperatures, respectively.

Boundary Conditions

A constant velocity inlet boundary condition with a static temperature of 300 K has been . . , 0 ( ) The heat transfer coefficient has been taken from [START_REF] Ho | Review of high-temperature central receiver designs for concentrating solar power[END_REF] and the ambient temperature is assumed to be 300 K.

conv s amb s z Q hA T T → = - (10)

Absorber Materials

Ceramics are the most suitable materials for VSR due to their high operating temperature and thermal stability. In the present study, two ceramic materials have been selected and studied. These include Alumina (Al 2 O 3 ) and SiC. SiC is widely employed in high temperature volumetric solar receiver [START_REF] Fend | Two novel high-porosity materials as volumetric receivers for concentrated solar radiation[END_REF][START_REF] Agrafiotis | Evaluation of porous silicon carbide monolithic honeycombs as volumetric receivers/collectors of concentrated solar radiation[END_REF], due to its high absorptivity (about 0.8) [START_REF] Sciti | Suitability of ultra-refractory diboride ceramics as absorbers for solar energy applications[END_REF][START_REF] Sciti | Tantalum diboride-based ceramics for bulk solar absorbers[END_REF], high operating temperature and good mechanical and physical properties. Alumina has been introduced in the present study for volumetric solar receiver. Alumina is a low cost well-balanced material offering high wear resistance, high mechanical strength, and thermal stability. In addition, as compared to SiC, Al 2 O 3 has high oxidation resistance at high temperature. However, the solar absorptivity of Al 2 O 3 is low (around 0.2) [START_REF] Sani | Optical properties of black and white ZrO2 for solar receiver applications[END_REF]. Optical analysis shall help in understanding the impact of its low absorptivity on the radiative flux distribution, the resulting optical efficiency and penetration depth of solar radiation in comparison with conventionally used materials such as SiC. Table 1 summarize the thermo-physical properties of SiC and Al 2 O 3 used in the present study. It should be noted that Silicon Carbide and Alumina used in volumetric solar receivers and implemented in the present study are opaque materials. Pure silicon carbide is colorless and transparent. However, the production of transparent, translucent, and opaque ceramics is obtained by controlling its microstructure and inclusion of impurities during the manufacturing process [START_REF] Kim | Microstructure comparison of transparent and opaque CVD SiC[END_REF][START_REF] Kim | Microstructure and optical properties of transparent alumina[END_REF]. The optical properties of Al 2 O 3 and SiC used in the present study have been reported in recent references [START_REF] Sciti | Suitability of ultra-refractory diboride ceramics as absorbers for solar energy applications[END_REF][START_REF] Sciti | Tantalum diboride-based ceramics for bulk solar absorbers[END_REF][START_REF] Sani | Optical properties of black and white ZrO2 for solar receiver applications[END_REF][START_REF] Auerkari | Mechanical and physical properties of engineering alumina ceramics[END_REF] dealing with characterization of materials for solar receivers. The emissivity of Alumina has been reported by P. Auerkari [START_REF] Auerkari | Mechanical and physical properties of engineering alumina ceramics[END_REF]. For wavelength up to 2.5 µm (wavelength of interest), the average value of emissivity is around 0.2 and it decreases with increasing the temperature. The temperature dependence of thermo-physical properties of SiC, Al 2 O 3 and air are taken into account as listed in table 1. = + [START_REF] Munro | Material properties of a sintered α-SiC[END_REF] Specific Heat ( p C ) (J. Kg -1 . K -1 )

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

0.003 1110 0.15 425

T P C T e - = + - [52]
Absorptivity ( α )(-)/Emissivity ( ε )(-) 0.8 [START_REF] Sciti | Suitability of ultra-refractory diboride ceramics as absorbers for solar energy applications[END_REF][START_REF] Sciti | Tantalum diboride-based ceramics for bulk solar absorbers[END_REF]] Alumina

Density ( ρ ) (kg. m -3 ) 3750 [53] Thermal Conductivity ( s K ) (W. m -1. K -1 ) 12227 4.5536 s K T = - + [54]
Specific Heat ( p C ) (J. Kg -1 . K -1 ) 197620 1429.4

P C T = - [54]
Absorptivity ( α )(-)/Emissivity ( ε )(-) 0.2 [START_REF] Barreto | Three-dimensional modelling and analysis of solar radiation absorption in porous volumetric receivers[END_REF][START_REF] Sani | Optical properties of black and white ZrO2 for solar receiver applications[END_REF] Air Viscosity ( µ ) (kg/m s)

Sutherland law [START_REF] Sutherland | The viscosity of gases and molecular force[END_REF] Thermal conductivity (

f K ) (W. m -1. K -1 )
f f f f K T T T - - - - = × - × + × - × [23]
Specific heat ( p C ) (J. Kg -1 . K -1 )

P f f f f C T T T T - - - - = × -× + × - × + × [23]

Numerical Solution Methods and Model Validation

A commercial software ANSYS fluent © has been used to solve the 3-D governing equations.

The SIMPLE algorithm and second-order upwind scheme have been used for pressure-velocity coupling and field variables interpolation from cell centers to the faces of the control volume.

A grid independence study has been carried out by analyzing the air outlet temperature ( air out T

) and the frontal surface temperature of the absorber ( solid in T ) for different grid sizes. The results are shown in table 2. It is clear that a grid with about 2,884,612 cells is fine enough.

Further refinements of the grid have no effect on the results. The results have been also checked as a function of the convergence criteria and independent results have been obtained

for a set value of 10 -6 for all discrete conservation equations. In order to validate the current optical and CFD models, a comparative study has been carried out with experimental studies available in the literature. Mey-Cloutier et al. [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF] reported experimental data for SiC honeycomb volumetric solar receiver. The experimental work has been carried out using solar furnace facility at PROMES-CNRS Laboratory, France. The Table 3 shows a comparison of air outlet temperature and absorbed solar heat flux on the frontal surface predicted by the current model and experimental data of Mey-Cloutier et al. [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF]. It can be observed that air outlet temperature and concentrated heat flux predicted from the present numerical model and those obtained by Mey-Cloutier et al. [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF] are in a good agreement.

Table 3. Comparisons of radiative heat flux and mean air outlet temperature from the honeycomb solar receiver with the experimental data of Mey-Cloutier et al. [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF].

Experimental data [START_REF] Mey-Cloutier | Experimental study of ceramic foams used as high temperature volumetric solar absorber[END_REF] Simulation results Deviation Concentrated solar heat flux (kW/m 2 ) 795 ± 3% 760.25 4.371%

Air outlet temperature (K) 1170 ± 10 K 1105.27 5.53%

Results and Discussions

The previous results have demonstrated the capability of the present coupled numerical optical and CFD models for the analysis and design of volumetric solar receivers. In the next sections, the models shall be employed to investigate the effects of different types of absorber materials material absorptivity, micro dimensions of the honeycomb structure, and air flow rate on the performance of the solar receiver.

Analysis of Absorbed Solar Heat Flux Distribution

Figure 5 shows a schematic of honeycomb channel and incident solar rays. As shown in Fig. 5 (a), the incidence angle of solar flux (β) is the angle between the normal direction to the frontal surface of the receiver and the line connecting the center of the heliostat module to the aiming point. In real commercial solar thermal power plant (actual and large receiver), thousands of heliostat modules are used to concentrate the incident solar radiation on the receiver surface and one volumetric solar receiver module may receive solar rays from multiple heliostats. As a result, there will be a range of incidence angle values and a dominant value is usually used in the analysis [START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF][START_REF] Barreto | Three-dimensional modelling and analysis of solar radiation absorption in porous volumetric receivers[END_REF][START_REF] Chen | Thermal and chemical analysis of methane dry reforming in a volumetric reactor under highly concentrated solar radiation[END_REF][START_REF] Nakakura | Efficiency and heat loss analysis of honeycomb receiver varying air mass flow rate and beam width[END_REF][START_REF] Nakakura | Conjugate radiation-convectionconduction simulation of volumetric solar receivers with cut-back inlets[END_REF]. This dominant value can be controlled by changing the tilt angle of the receiver (θ R ). As can be seen in Fig. 5(b), the value of the incidence angle (β) should be higher than zero to ensure that incident solar rays are not parallel to the axis of the honeycomb channel and exit the receiver without being absorbed. Previous studies show that an incidence angle less than 10 degrees will enhance the radiation propagation within the honeycomb receiver, however, some rays will exit the receiver without being absorbed [START_REF] Cagnoli | Coupled optical and CFD parametric analysis of an open volumetric air receiver of honeycomb type for central tower CSP plants[END_REF]. Nakamura et al. [START_REF] Nakakura | Efficiency and heat loss analysis of honeycomb receiver varying air mass flow rate and beam width[END_REF] concluded that the effect of incidence angle from 10 O to 45 O on the receiver efficiency is negligible. This can be attributed to very small values of channel height as compared to channel length. Figure 6 Figure 7 shows a detailed absorbed heat flux distribution on different surfaces of the honeycomb channel as a function of the axial depth (Z). For a higher absorptivity material (SiC of α = 0.8), the majority of concentrated solar radiation is intercepted by the upper surface and the flux peak is located at the entrance of the Honeycomb channel (

0 Z mm =
). The reflected radiative power from the upper surface is absorbed by the lower surface. Thus, the radiative flux peak is located on deeper sections ( Z around 7 mm). In addition, the maximum penetration depth of the radiative flux within the honeycomb absorber is less than 10 mm.

Such a flux distribution with small penetration depth may lead to higher temperature at the frontal surface of the receiver as compared to the air outlet temperature from the receiver.

This results in relatively higher values of thermal losses and lower values of solar to thermal efficiency.

In case of material with low absorptivity (Al 2 O 3 ), the percentage of the reflected radiative solar power from the upper surface to the lower surface is higher when compared with the percentage share of lower surface for a higher absorptivity material (SiC). However, the total absorbed radiative power by the frontal surface (see Fig. 8) and internal walls is very small compared to SiC Optical analysis shows a multi-reflection effect in the absorber channel of For the purpose of increasing the penetration depth of the radiative heat flux and the optical efficiency of Al 2 O 3 honeycomb absorber, the frontal and internal surfaces of the honeycomb channel can be coated with a high absorptivity selective coating such as Pyromark 2500 © hightemperature paint [START_REF] Gray | High temperature performance of highefficiency, multi-layer solar selective coatings for tower applications[END_REF]. The solar absorptivity and emissivity of the coating are 0.95 and 0.9, respectively [START_REF] Ho | Characterization of Pyromark 2500 paint for high-temperature solar receivers[END_REF]. Since most of the optical losses are due to reflection losses from the frontal surface and photons escaping from the exit of the receiver, it is proposed to apply a selective coating (Pyromark 2500) on the frontal surface (

0 Z mm =
) and the last 25 mm depth from the exit (

50 Z mm =
) of Al 2 O 3 honeycomb absorber (see Fig. 9).

It should be noted that Pyromark 2500 high-temperature paint has been used on central receivers, including the Solar One Central Receiver Pilot Plant. As reported by the indicate that the solar absorptance of the coating can decrease with time at elevated temperatures, and repainting may be required every few years [START_REF] Ho | Characterization of Pyromark 2500 paint for high-temperature solar receivers[END_REF]. Pyromark 2500 was implemented in the analysis performed in the present study in order to employ realistic and reasonable values of solar absorptivity and emissivity of high temperature coating materials.

The development of solar selective coatings for high-temperature solar receivers is currently a topic of research [START_REF] Gray | High temperature performance of highefficiency, multi-layer solar selective coatings for tower applications[END_REF]. Other superior coating materials are most likely expected to be produced. However, the results of the present analysis shall be applicable to other coating materials of similar values of optical properties.

Employment of the present suggested enhancement method, shows that the absorbed heat flux on the frontal surface is increased to reach a value higher than that of SiC (See Fig. 8,C-

Al 2 O 3 )
. Moreover, a maximum peak of the radiative heat flux exists in all surfaces after 25 mm from the entrance of the honeycomb channel (see Fig. 10). The second peak in absorbed solar flux occurs, as expected, at the beginning of the suggested coating layer. Since the optical efficiency of AL 2 O 3 is below 20% as can be seen in Fig. 8. This means that about 80% of the incident solar radiation is still available for absorption in the coated part of the channel. The absorption of this huge number of photons results in a second peak in absorbed solar heat flux that may exceed the first peak occurring at the channel inlet section. Using this technique results in preventing photons from escaping from the exit of the receiver, increasing the penetration depth of the radiative heat flux, decreasing the reflection losses from the frontal surface of the absorber and consequently increasing the optical efficiency to reach values of about 87.56%. For the range of in u studied, both air outlet temperature and solar-to-thermal efficiency of C-Al 2 O 3 honeycomb absorber proposed in the present work are higher than those of SiC and Al 2 O 3 absorbers. This can be attributed as due to the higher penetration depth of the solar heat flux and the lower convective and radiative heat losses from the frontal surface. This performance is advantageous when one considers the low cost and high-temperature oxidation resistance of Al 2 O 3 as compared to conventional SiC honeycomb absorbers. In a first evaluation of the effectiveness of the current proposed method, the authors suggested applying high absorptivity coating at the surface of exit half of absorber channel made of low absorptivity material (AL 2 O 3 ). The results show an increase in the optical efficiency, penetration depth of absorber heat flux, positive volumetric effect, and high solar to thermal efficiency. Based on the present analysis, the authors believe that the suggested location of coated section should be always applied away from the channel entrance. Further parametric analysis for the optimization of this coated length shall be carried in the future.

In previous studies [START_REF] Fend | Numerical investigation of flow and heat transfer in a volumetric solar receiver[END_REF], the effects of radiative exchange between the inner walls of the honeycomb channel have been neglected. Figure 15 shows the impact of radiative emission on the solar to thermal efficiency of honeycomb structure receiver made of different materials.

Overall, it has been observed that neglecting the radiative exchange will lead to overestimation of solar to thermal efficiency. The difference between • Optical, CFD models for the analysis and design of volumetric solar receivers.

• Effects of absorber material and air flow rate on the performance are analyzed.

• A novel honeycomb receiver made of Alumina has been introduced.

• Coated alumina honeycomb absorbers show a favorable receiver performance.

Fig. 2 .

 2 Fig. 2. Schematic diagram of the solar furnace optical model with volumetric solar receiver module.

Figure 3 .

 3 Figure 3. CFD solution domain.

Figure 4 .

 4 Figure 4. Numerical solution domain with generated mesh

kcS

  Thermal conductivity (W/m K) p Specific heat (J/kg K) v S Volumetric heat source (W/m 3 ) In the present model, energy equations for the solid and fluid phases are coupled in a conjugate heat transfer problem. The radiative heat flux distribution obtained from the optical simulation has been implemented as a volumetric source in the solid phase energy equation using a User Defined Function (UDF). This source function represents the inter-linkage between the optical and CFD models. The present volumetric source term in implemented only in the surface control volumes of the solid phase. This is a well-known practice in modeling of surface heat flux by an equivalent volumetric source term [43]. The solar radiation source term has been calculated using: Volumetric heat source (W/m 3 ) PPP Power per photon (W) . abs Ν Number of absorbed photons (-) s V Volume of solid surface control-volume at the solid surface (m 3 )

  used at the fluid domain entrance and a zero-gauge pressure outlet boundary condition has been used at the domain outlet. All contact walls between the fluid and solid domains are set to no-slip non-penetrating walls. All other surfaces have been set to symmetrical boundary conditions. At the frontal surface of the receiver, the convection and radiation losses represent the boundary conditions for energy equation. The convection losses and radiation losses

  following parameters have been collected and used as model parameters: DNI = 1000 W/m 2 , reflector area = 27 m 2 , concentrator area = 13.45 m 2 , sample thickness = 6 cm, sample diameter = 5 cm, atmospheric pressure = 84.5 kPa and mass flow rate of air ( m & ) = 1 g/s

  Section 4.1 is devoted to the analysis of absorbed solar heat flux distribution on different surfaces of the honeycomb solar receiver as obtained from the optical modeling. Two different materials with high and low absorptivity, namely SiC and Al 2 O 3 are employed. Optical analysis shall help in understanding the influence of radiative flux distribution, the resulting optical efficiency and penetration depth of solar radiation. New enhancement methods to improve the performance of honeycomb solar receivers are suggested. Sections 4.2 and 4.3 are devoted to the analysis of heat transfer and fluid flow inside the receiver and their effects on the overall receiver performance.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Schematic of honeycomb channel and incident solar rays: (a) Solar incidence angle, (b) Incident solar ray tracing, (c) Channel surfaces

Al 2 O 3 . 3 Fig. 7 .

 2337 Fig. 7. Absorbed solar heat flux distribution on internal surfaces of honeycomb absorber channel for: (a) SiC; (b) Al 2 O 3 .

Fig. 8 .

 8 Fig.8. Effect of material absorptivity on the absorbed heat flux on frontal surface of honeycomb absorber and optical efficiency.

Fig. 9 . 21 Fig. 10 . 6 . 2 Flow

 9211062 Fig.9. Schematic of present suggested enhancement method using Alumina honeycomb absorber coated with Pyromark 2500 © .

Fig. 11 . 22 Figure 12

 112212 Fig. 11. Velocity (a) and temperature fields (b) on the axial middle plane of SiC honeycomb absorber.

3 Fig. 12 .

 312 Fig. 12. Temperature field on the axial middle plane of a honeycomb absorber made of (a) Al 2 O 3 and (b) C-Al 2 O 3 .

Figure 13

 13 Figure 13 depicts the mean temperature distribution of the solid and fluid phases of honeycomb structure receiver made of different materials as a function of the absorber depth using 0.4741 / in u m s = . As mentioned before, in case of honeycomb receiver made of Alumina (low absorptivity material) the majority of the concentrated sunlight on the frontal surface of the receiver is reflected back to the ambient, and the photons can run away from the exit of the receiver due to the low absorptivity of the material. Therefore, as can be seen from Fig. 13, air outlet temperature is very low (around 690 K) in comparison with SiC and C-Al 2 O 3 . Moreover, the solid and fluid phases do not reach the thermal equilibrium. There is a temperature difference between the two phases even at the exit of the receiver. This leads to unfavorable receiver performance in terms of air outlet temperature and thermal efficiency. On the other hand, coated alumina (C-Al 2 O 3 ) honeycomb absorbers show favorable receiver performance in

Fig. 13 . 6 . 3

 1363 Fig.13. Temperature distribution along the centerline of a honeycomb absorber made of different materials with u in = 0.4741 m/s.6.3 Evaluation of Optical and Thermal Performance of the Honeycomb Receiver

m

  & ) leads to a decrease in air outT , and a decrease of the average air temperature within the honeycomb channel. This leads to a decrease in the temperature of the internal walls of the honeycomb channel as well as the frontal surface temperature. Moreover, decreasing the frontal surface temperature with the increase of in u leads to a reduction in the convective and the ambient. Therefore, it can be seen that the solarto-thermal efficiency of the receiver ( ) . s Th → η increases with the increase of in u .

27 Fig. 15 .

 2715 Fig. 15. Impact of radiative emission on the solar to thermal efficiency of honeycomb structure receiver made of different materials.

1 (Figure 2 :

 12 Fig.14. Variation of: (A) air outlet temperature; (B) temperature of the absorber at the entrance; (C) solar-to-thermal efficiency with inlet velocity.

14 Table 1 .

 141 Thermo-physical properties of SiC, Al 2 O 3 , and air.

	Silicon Carbide
	Property Density ( ρ ) (kg. m -3 )	Value 3100				Ref. [52]
	Thermal Conductivity ( s K ) (W. m -1. K -1 )	K	s	52000 e T	1.24 10 437	5	T

--×

Table 2 .

 2 Grid independence study for CFD analysis.

	Grid size	air out T	(K)	solid in T	(K)


Marcos Vinícius

Santos Issa

Americo Cunha

Francisco Soeiro

Cross-Entropy Method

Cross-entropy Method In Structural Optimization with Dynamic Constraints

come   L'archive ouverte pluridisciplinaire

 Gradient based methods are not possible in some cases.

 Metaheuristics may be an alternative, but may have high computational cost or may be prohibitive.

 Cross-entropy method (CE) has been used successful in combinatorial optimization and estimation of rare events in the last two decades  Due to complex geometric configurations, structural optimization can be extremely non-linear, requiring the use of very efficient optimization algorithms.

 Propose a Cross-Entropy framework for structural optimization and investigate its accuracy and efficiency.

 Compare optimal values found by the CE with other results obtained by Sequential quadratic programming (SQP) and Genetic algorithm (GA). [2] R. T. Haftka and Z. Gürdal, Element of Structural Optimization, 3rd ed., Kluwer Publishers, Dordrecht, Netherlands, 1992. [3] M.V.S Issa, A. Cunha Jr, F. J. C. P. Soeiro, A. Pereira, Structural optimization using the Cross-entropy method, In XXXVIII Congresso Nacional de Matemática Aplicada e Computacional, Campinas-SP, Brazil, 2018. https://dx.doi.org/10.5540/03.2018.006.02.0443  Only global search algorithms should be used to obtain optimal solutions for these optimization problems with dynamic constraints.

METHOD

 In general, the SQP obtains better results than the GA and the CE. However, in this case, the metaheuristic, GA, obtain better results than those of SQP, which is a gradient-based optimization method, and CE has a satisfactory result.

BENCHMARK TEST Nucleus of Modeling and Experimentation with Computers (NUMERICO)



  Structural optimization has the typical objectives as mass reduction (weight), change natural frequencies (avoid resonance), reduce internal stress and REDUCE COST.

FIGURE 1 -

 1 FIGURE 1 -Benchmark test.

FIGURE 2 -

 2 FIGURE 2 -Structural optimization.

FIGURE 3 -

 3 FIGURE 3 -CE algorithm.

TABLE 1 -

 1 Comparison between the results obtained with different techniques optimization the area of each bar.

Source: HAFTKA; GÜRDAL , 1992

[1] R. Y. Rubinstein and D. P. Kroese, Simulation and the Monte Carlo Method, Wiley, 3rd edition, 2017.

TABLE 2 -

 2 Illustration of the areas obtained by the three optimization methods considering natural frequencies constraint.


Michel Tosin

Americo Cunha

Flávio Codeço Coelho

Flávio Codeço

Coelho Sensitivity Analysis

In Zika

Sensitivity Analysis in Zika Virus Dynamics and a Model Discrepancy Approach

 In future works, the authors intend explore this new framework to quantify model discrepancy and then improve its predictions;

  INTRODUCTION São Paulo, July 29th -August 9th, 2019 Suppose a data set of measures of the . The -th observation is given by  Inference problem become use data information to update the prior probability density func-tion(PDF), defined for . The solution corresponds posterior PDF;  From Bayes' rule, To define a good point of start, the Maximum Entropy Principle is applied to construct the most informative prior distribution. FINAL REMARKS  Comparative results of global Sobol' Indices show how the lack of some parameters can change the sensibility effect of the others;  A framework for statistical inference exploring Polynomial Chaos to measure the model discrepancies was presented;

[ 1 ]FIGURE

 1 FIGURE 1 -Zika transmission representation.

FIGURE 4 -

 4 FIGURE 4 -Sobol' indices of the model. FIGURE 5 -Sobol' indices with constant.

FIGURE

  FIGURE 7 -Sobol' indices without .


Victor Bucarey
email: vbucarey@ulb.ac.be

Martine Labbé
email: mlabbe@ulb.ac.be

Discussion of Fairness and Implementability in Stackelberg Security Games

Keywords: Fairness, Implementability, Stackelberg Security Games

In this article we discuss the impact of fairness constraints in Stackelberg Security Games. Fairness constraints can be used to avoid discrimination at the moment of implementing police patrolling. We present two ways of modelling fairness constraints, one with a detailed description of the population and the other with labels. We discuss the implementability of these constraints. In the case that the constraints are not implementable we present models to retrieve pure strategies in a way that they are the closest in average to the set of fairness constraints.

Introduction

In the last years, Stackelberg Security Games have been applied in several real domains such as airport security [START_REF] Pita | Deployed armor protection: the application of a game theoretic model for security at the los angeles international airport[END_REF], IRIS for security of flights [START_REF] Jain | Security games with arbitrary schedules: A branch and price approach[END_REF], ports [START_REF] Shieh | Protect: A deployed game theoretic system to protect the ports of the united states[END_REF] and border [START_REF] Bucarey | Building real stackelberg security games for border patrols[END_REF] patrolling and fare evasion [START_REF] Correa | Fare evasion in transit networks[END_REF], among others. In these games, the leader, also called defender must protect a set of targets with limited resources available from a possible attack performed by the follower or attacker. The payoff structure depends only on whether the target attacked is being protected or not [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF]. Each defender resource has a set of possible schedules, that is, the possible subset of targets they can protect in one strategy. If all the resources have the same possible set of schedules we name this game with homogeneous resources. Otherwise, it is called with heterogeneous resources. Several mixed integer formulations, for both, general Bayesian Stackelberg games and Bayesian Stackelberg Security games are presented in [START_REF] Paruchuri | Playing games for security: An efficient exact algorithm for solving bayesian stackelberg games[END_REF], [START_REF] Casorrán | A study of general and security stackelberg game formulations[END_REF].

One of the key points of the scalability of SSG is the representation of the set of strategies of the defender. Instead of taking into account every single pure strategy, this set is represented through the frequency in which each target is protected. These frequencies are called coverage distribution.

In some applications, payoffs matrices are generated in a black-box. In others, they are built by real data. In both cases these payoffs can lead to discriminative outcomes. In the first case, this is due to the fact that there is not a known methodology. In the second case, this is due to the manner in which data is gathered. For instance, according to [START_REF] Pierson | A large-scale analysis of racial disparities in police stops across the united states[END_REF], officers generally stop black drivers at higher rates than white drivers, and stop Hispanic drivers at similar or lower rates than white drivers. Another example came from the United Kingdom. According to the Governamental institution Ethnicity in the UK 3 , in 2016/17, there were 4 stop and searches for every 1,000 White people, compared with the 29 stop and searches for every 1,000 Black people. Among specific ethnic groups, the Chinese and Mixed White/Asian groups consistently had the lowest rates of stop and search since 2009/10 [START_REF] Hargreaves | Police Powers and Procedures, England and Wales, Year Ending 31[END_REF].

The relationship between algorithms and discrimination is a major concern for the European Union Agency for Fundamental Rights -FRA [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF]. They state:

When algorithms are used for decision making, there is potential for discrimination against individuals. The principle of non-discrimination, as enshrined in Article 21 of the Charter of Fundamental Rights of the European Union (EU), needs to be taken into account when applying algorithms to everyday life.

And later,

Big data claims to be neutral. It isn't . . . machine learning depends upon data that has been collected from society, and to the extent that society contains inequality, exclusion or other traces of discrimination, so too will the data. Consequently, unthinking reliance on data mining can deny members of vulnerable groups full participation in society.

While there are several real world problems where fairness is important, there is no a specific way to measure fairness. For instance, for classifications problems, it is used as fairness constraints the following considerations: disparate treatment which implies that the probability of any classifier output does not change after observing a sensitive feature; disparate impact which implies that the probability of classifying with a positive value; and disparate treatment which implies that the probability of misclassification does not change with some sensitive feature [START_REF] Zafar | Fairness constraints: A flexible approach for fair classification[END_REF].

In the SSG context, bias in the data can be translated into allocate more/less surveillance focused in race, wealth/poverty or any other type of discrimination.

In this work we study different ways to include constraints in order to avoid discrimination SSGs from a tactical point of view. We are interested in studying how these considerations could be implemented and how much we lose by adding these considerations in terms of expected utilities.

SSG can be seen as the problem of allocating resources to targets in a random fashion. A relevant work that came from the problem of designing randomized allocations is presented in [START_REF] Budish | Designing random allocation mechanisms: Theory and applications[END_REF]. In this context, authors define the concept implementability of a set of constraints when any random allocation under this set of constraints can be decomposed in a convex combination of pure allocations satisfying this set of constraints. Sufficient and necessary conditions for the implementability of a set of constraints are given, based on the bihierarchy structure of the constraints.

Sometimes, there are constraints that are necessary in real applications but they are not implementable in general. In this article we show that including fairness constraints in a detailed description of the population inside each target might not be implementable. That means, that in average solutions satisfy those constraints but the pure strategies that implement this solution could not satisfy them. We will show some models and algorithms to retrieve pure strategies minimizing the violations of such constraints.

The questions that we aim to answer are the following:

-Is it possible to model coverage distributions including fairness considerations? -Are they implementable? If not, how can we include those considerations in practical settings? -How much does the defender lose, in terms of expected utility, by including these considerations?

Our first contribution is to model fairness constraints in the coverage probabilities in SSGs. We present two models, one focused on a detailed description of the population, the second one based on labels on the targets. Our second contribution is to show that the model based on labels is implementable. Our third contribution is to provide a methodology for implementing schedules allocating low probability to strategies that violate more the set of non-implementable constraints.

The rest of the paper is as follows. In Section 2 we introduce the main concepts related to SSGs, Implementability and Random Allocations. We give an introductory example for this problem. In Section 3, we provide the models that are discussed in this work. In Section 4 we provide a discussion about the implementability of the coverage distribution returned by these models. Also, some extensions are presented. In Section 5 computational results are shown. Our models are tested in a realistic instance presented in Section 6. Finally, our conclusions are presented in Section 7.

Problem Statement and Notation

Stackelberg Security Games and compact formulations

In SSGs, the leader, named in this context defender, must protect from the followers, named attackers, a set of targets J. The payoffs for the players only depend on whether a target is protected or not. In consequence, several strategies have identical payoffs. Thus, we denote by D k (j|p) the utility of the defender when an attacker of type k ∈ K attacks a covered target j ∈ J and by D k (j|u) the utility of the defender when an attacker of type k ∈ K attacks an unprotected target j ∈ J. Similarly, the utility of an attacker of type k ∈ K when successfully attacking an unprotected target j ∈ J is denoted by A k (j|u) and that attacker's utility when attacking a covered target j ∈ J is denoted by A k (j|p). We denote π k the probability of the defender facing attacker k.

In the heterogeneous resources setting there is a set Ω of resources, |Ω| = m, in which each one can be allocated to a possible subset of targets J ω . If J = J ω for each ω ∈ Ω, we call it homogeneous resources. A pure strategy i ∈ I for the leader is an allocation of resources to targets. That is

I =    {a ω ∈ {0, 1} |Jω| } ω∈Ω : j∈Jω a ωj ≤ 1 ∀ ω ∈ Ω   
in the case of heterogeneous resources, or

I =    a ∈ {0, 1} |J| : j∈J a j ≤ m   
in the case of homogeneous resources. In this context, in the homogeneous case |I| = m k=1 n k . Authors in [START_REF] Kiekintveld | Computing optimal randomized resource allocations for massive security games[END_REF] provide a compact formulation using the transformation

c j = a∈I:aj =1 x a , (1)
where c j represents the frequency coverage of target j and x a represents the probability of playing strategy a ∈ I. The formulation for computing a SSE in the homogeneous case (HOM) is stated as follows:

(HOM) max c,q,f,s

k∈K π k f k (2) s.t. f k ≤ D k (j|p)cj + D k (j|u)(1 -cj) + M (1 -q k j ) j ∈ J, k ∈ K(3) 0 ≤ cj ≤ 1 j ∈ J (4)
j∈J cj = m (5)

j∈J q k j = 1 j ∈ J, k ∈ K(6) 0 ≤ s k -A k (j|p)cj -A k (j|u)(1 -cj) ≤ M (1 -q k j ) j ∈ J, k ∈ K(7) f k , s k ∈ R (8)
q k j ∈ {0, 1} j ∈ J, k ∈ K(9)
where c j represents the frequency with which target j is protected. Variables f k represent the expected utility for the defender of facing attacker k and s k the expected utility of attacker k. The objective function (2) to be maximised represents the expected utility of the defender. Expression (3) states an upper bound for f k which is tight for the target selected by each attacker k. Expression (4) and ( 5) define the coverage probabilities c j . Expression [START_REF] Hargreaves | Police Powers and Procedures, England and Wales, Year Ending 31[END_REF] states that each attacker selects a pure strategy. For each type of attacker, expression [START_REF] Jain | Security games with arbitrary schedules: A branch and price approach[END_REF] states that the target attacked maximizes their expected utility. An extension to the heterogeneous case is stated by introducing variables c ωj satisfying

c j = ω∈Ω:j∈Jω c ωj j ∈ J, (10)
and adding constraints limiting the amount of coverage for every single resource:

j∈Jω c ωj ≤ 1 ω ∈ Ω. (11)
By solving these optimization problems, a coverage distribution is obtained. In order to obtain a mixed strategy x in the original space I that fits with (1) the following method described in Algorithm 1 could be applied. An example is presented in Figure 1.

Algorithm 1 Box Method

Require: c ∈ R |J| feasible coverage.

Step 1: For each resource ω ∈ Ω, consider a column of height 1.

Step 2: For each resource ω, fill up the column with the values of cj, with j ∈ Jω

Step 3: Define x by extending each rectangle line into a horizontal line crossing all columns. The area between two horizontal lines represents a defender strategy. This area identifies a set of targets protected, at most one for each resource ω. The height of this area represents the probability of the corresponding strategy. x {1,4} = 0.1

x {1,3} = 0.4 x {2,4} = 0.2 x {2,5} = 0.2 x {3,5} = 0.1 (c)

Population and unfair allocations

The optimal solution found through the optimisation problem stated above relies on how the payoffs matrices are computed. Anyway, discrimination issues can be performed by the security agency, putting less (or more) resources over some groups. We study two models to avoid such situation.

The first model considers that in each target j ∈ J there is a population p jt of type t between a set of possible population types T . Examples of T can include different races, religions, income level, among others. In order to avoid discrimination issues, it could be desirable that the total coverage allocated to each population type is proportional to the amount of inhabitants.

A second model can be developed with a slightly different description of the population. Instead of considering a fraction of population in each target, consider that each target has a label L. Each label ∈ L denotes the most representative population in the target. For example, the Latin, Asian and African areas in the main cities in Europe. Also, in some cities there is a division in High income, Medium class and Low income areas.

In the following example, we show that even in small instances, the problem of unfairness may occur. We will use this example along the article to introduce the main concepts of this work.

Example 1. Consider the following instance with five targets, J = {j 1 , j 2 , j 3 , j 4 , j 5 }, three attacker types, K = {k 1 , k 2 , k 3 } and the total population divided in three types, T = {t 1 , t 2 , t 3 }. Payoffs of the game and the description of the population in each target is represented in Table 1. We consider π = (π k1 , π k2 , π k3 ) = (0.5, 0.3, 0.2) the distribution of probability over the set of the attackers. We aim to allocate two homogeneus resources, i.e., m = 2. In the last row of Table 1.(b) a label for each target is stated. In this example, a target is labeled by l i if the population type t i is the most representative. We label each target in L = { 1 , 2 , 3 } corresponding with the most representative population type. In other words, j = i if i = arg ti∈T max p jti . It would be desirable that the coverage that each part of the population receives is proportional to the total population that they represent.

k1 k2 k3 j D k 1 A k 1 D k 2 A k 2 D k 3 A k 3 j1 0 -
Table 2.(a) shows the coverage given by model (HOM) and Table 2.(b) shows how this coverage is allocated to each part of the population. The proportional . Population type t 1 receives proportionally 97.143% more coverage than the population that they represent while population type t 3 receives a -42.4% less. Note that 1 represents the 20% of the targets, and 2 and 3 the 40% each one. In Table 2.(c), we show what is the proportion of resources allocated to each label. In this instance, both evaluations give an unbalanced allocation of the resources.

Random Allocations, Implementability and Bi-hierarchies.

We now introduce some mathematical notions that allow us to develop the models and results of this work. The problem of finding an SSE can be seen as finding a specific set of random allocations between resources Ω and a set of targets J under some considerations. Then, the coverage vector c ωj can be seen as a random allocation as presented in [START_REF] Budish | Designing random allocation mechanisms: Theory and applications[END_REF]. They study constraints of the type:

q L S ≤ (ω,j)∈S c ωj ≤ q U S S ⊆ Ω × J,
where S is called constraint set and q L S , q U S are positive integers named quotas on S. The full of set of constraints is named constraint structure H.

A random allocation c is implementable under the constraint structure H and quotas q = {(q L S , q U S )} S∈H if it can be written as a convex combination of pure allocations feasible under H and quotas q. Constraint structure H is universally implementable if, for any integer quotas q, every random allocation satisfying constraints in H is implementable under q.

A constraint structure is a hierarchy if for any pair S, S either S ⊆ S or S ⊆ S or S ∩ S = ∅. A constraint structure H is a bihierarchy if it can be partitioned in two hierarchies H 1 and H 2 , that is,

H 1 ∩H 2 = ∅ and H 1 ∪H 2 = H.
Necessary and sufficient conditions over the constraint structure and implementability are given by [START_REF] Budish | Designing random allocation mechanisms: Theory and applications[END_REF] through the following two theorems: Theorem 1. (Sufficiency) If a constraint structure is a bihierarchy and quotas q are integers, then it is universally implementable. Theorem 2. (Necessity) If a constraint structure contains all the rows and columns constraints and is not a bihierarchy, then it is not universally implementable.

In the SSG context, the implementability of a set of constraints means that each coverage distribution feasible can be decomposed in a convex combination of pure strategies, or pure patrols, all of them satisfying the constraints. Note that if we only restrict to allocate resources to targets, the set of constraints forms a bihierarchy and all the coverage distributions c are implementable.

Models

In this section we discuss about how to restrict the coverage distributions taking into account the issues described in Section 2.2. We describe two ways of modelling fairness constraints: First we model constraints with a detailed description of the population; Then, we restrict the possible coverage considering aggregated information in terms of labels in each target. By doing this, we generate coverage probabilities that are not significantly correlated with sensitive features, as race, income level, etc.

Focus on the population

In this setting we assume a description of the population in each target j, given by the percentage of population of type t ∈ T denoted by p jt . In order to avoid discrimination issues, we might consider to restrict the amount of coverage that each population receives.

That is, the coverage vector c should satisfy

q L t ≤ j∈J c j pjt ≤ q U t t ∈ T. (12)
where q = (q L t , q U t ) are the quotas of the total coverage performed over the population t ∈ T .We denote pjt = pjt t ∈T p jt the fraction of population type t in target j. In this work, we use quotas of the form:

q L t = (1 -α)m j∈J pjt q U t = (1 + α)m j∈J pjt ( 13
)
where α is the maximum acceptable percentage of deviation from the total fraction of the population t multiplied by the number of the resources available (e.g., 10%).

In constraints [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF], we assume that the total coverage inside a target is distributed proportionally to each population type. This assumption can be relaxed, through introducing some nonlinear relationship between the coverage and the probability of being covered given the composition of the population inside target j. Anyway, this topic is out of the scope of this paper, but it is still an interesting research question.

As we mentioned before, constraints [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF] are not universally implementable. This is mainly because these constraints do not induce integer extreme points. In consequence, the vector coverage c under its form can not be decomposed in terms of pure strategies satisfying constraints [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF]. We will discuss how to deal with this issue in Section 4.

Label focus

A different approach can be stated as follows. There exists for each target a label ∈ L, representing a type of population representative in that target. In that case, we should think about protecting targets with an amount of coverage proportional to the percentage of the population they represent. Information about population is aggregated in each target, and in consequence, it is a more relaxed way of modelling fairness constraints.

Formally, define J the targets labeled by . Note that {J } ∈L defines a partition of J. For each label , there is a minimum and a maximum number of resources assigned to protect zones in J , denoted again by q L and q U respectively. Then, constraints on the amount of coverage can be stated as:

q L ≤ j∈J c j ≤ q U ∈ L. ( 14
)
where the quotas q are given by:

q L = (1 -α)m |J | |J| q U = (1 + α)m |J | |J| (15)
In this case, we use integer quotas to establish the implementability result in Theorem 3. In the model focused on the population, even using integer quotas the implementability result does not hold. Example 1 continued Consider α = 25%. In Table 3, we show the bounds for each type of population and labels for both models respectively. Figure 2(a) shows the optimal coverage in each target by including the fairness constraints in both models and the coverage given by (HOM). For a fixed α, the model on labels is less restrictive than the one focused on the population. This explains the difference between the optimal coverage given by (HOM). By the same reason, the difference in terms of Defender Expected utility for different values of α for the model of population is greater than the model with labels, as is shown in Figure 2(b).

The optimal coverage, in both cases should be implemented by sampling pure strategies. We showed in Table 4, one possible decomposition of these solutions Type q L t q U t Label q L q U t1 0.21 0.35 1 0 1 t2 0.54 0.9 using Algorithm 1. Note that strategies in model focused on population are not implementable. In particular, the third strategy in Table 4.(a) that covers targets j 2 and j 4 allocates 0 resources to the population type t 3 . 4. Decomposition for the optimal coverage in (a) model focused on population and (b) model focused on labels, using Algorithm 1. In both cases we use α = 25%.

Implementability and extensions

In this section we discuss the implementability of both models. First, we show that if quotas used in the model of labels are integers, then every coverage satisfying fairness constraints are implementable. Then, we discuss how to find pure strategies that are the closest to the set of fairness constraints in the model of population and that fits with the coverage distribution. Finally, we discuss about some extensions of the models and when they preserve this property.

Labels are implementable, Population constraints do not

Now we show that the model focused on labels is implementable. Theorem 3. If for each ∈ L, q L , q U are integers, then conditions ( 14) are universally implementable.

Proof. We prove that the set of constraints forms a bihierarchy in the problem of allocation targets to resources. In particular, any vector coverage satisfies:

0 ≤ c ωj ≤ 1 ω ∈ Ω, j ∈ A(ω) (16) 0 ≤ c ωj ≤ 0 ω ∈ Ω, j / ∈ A(ω) (17) 0 ≤ j∈J c ωj ≤ 1 ω ∈ Ω (18) 0 ≤ ω∈Ω c ωj ≤ 1 j ∈ J (19) q L ≤ ω∈Ω j∈J c ωj ≤ q U ∈ L. (20)
Conditions ( 16) and ( 17) represent singletons {(ω, j)} ∈ Ω × J. Condition (18) can be represented as the sets {ω} × J. Condition ( 19) is represented by Ω × {j}. Finally, conditions (20) can be represented as the sets {(ω, j) : j ∈ J , ω ∈ Ω}. We show that by grouping conditions (16), ( 17) and (18), and grouping conditions (19) and (20), they form a bihierarchy. Formally, we define the following sets:

H 1 = {{(ω, j)} ∈ Ω × J} ∪ {{ω} × J : ω ∈ Ω} H 2 = {Ω × {j} : j ∈ J} ∪ {{(ω, j) : j ∈ J , ω ∈ Ω} : ∈ L}
First, we show H 1 is a hierarchy: Clearly each pair of singleton are disjoint. The same occurs with each pair of elements of {{ω} × J : ω ∈ Ω}. On the other hand, each singleton (ω, j) is, either included in the set {ω}×J and it has empty intersection with any other set {ω } × J, with ω = ω. Then, H 1 is a hierarchy. Now we prove that H 2 is a hierarchy. Each pair of elements in {Ω × {j} : j ∈ J} are disjoint. The same occurs with each pair of elements in {{(ω, j) : j ∈ J , ω ∈ Ω} : ∈ L} because {J } induces a partition over J. Now, if we take a pair of element of each group of conditions, for index j and , there are two cases: Either j ∈ J , in that case Ω × {j} ⊆ {(ω, j) : j ∈ J , ω ∈ Ω}; or j / ∈ J , in that case Ω × {j} ∩ {(ω, j) : j ∈ J , ω ∈ Ω} = ∅. Then, H 2 is also a hierarchy.

Clearly H 1 ∩ H 2 = ∅. Then, the set of conditions B = H 1 ∪ H 2 , forms a bihierarchy, and then under any integer quotas, the expected allocation is implementable and the result follows. A graphical representation of the proof is shown in Figure 3

                                      (a)
                                                        (b)

Approximating the Implementability

Coverage frequencies in the model focused on population are not implementable in general. That means that they cannot be decomposed in pure strategies satisfying the fairness constraints [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF]. In any case, they can be decomposed in pure strategies in the original set of strategies, that is, pure strategies covering at most m targets. This can be preformed by Algorithm 1. On one hand, this algorithm generates different decompositions by considering different orders in which the set of targets are included. For instance, in Figure 1, targets were included in lexicographic order. If the algorithm considers order j 2 , j 3 , j 1 , j 4 , j 5 , the output probabilities x will be different. On the other hand, Algorithm 1 does not take into account if the strategies produced satisfy conditions such as fairness, or are close to satisfy them. We would like to get a decomposition such that is the fairest as possible, allocating low probability to strategies that are unfair and high probabilities to the fairest ones. Formally, we have a polyhedron P 1 as the convex hull of the binary encoding of the set I of pure strategies. Let P 2 be the polyhedron of all the coverage vectors satisfying the fairness constraint in the model focused on population, that is P 2 = {c ∈ P 1 : c satisfies (12)}. We want to find a convex decomposition of a point in P 2 in terms of vertices of P 1 such that the weighted sum of the violations of constraints ( 12) of each strategy is minimised. The weights used in this optimisation problem come from the convex decomposition. By doing this, we aim to achieve a set of strategies implementing the optimal fair coverage, but at the same time, allocating low probability to the strategies that are unfair. Now we present some models to find such decomposition. We formulate the following non-linear model in order to get a decomposition in pure strategies where each pure strategy minimises the violation of the constraint [START_REF] Rights | Bigdata: Discrimination in data-supported decision making[END_REF]. Let M = {1, . . . , U } be a set of indices, where U is an upper bound on the number of strategies needed to decompose c in terms of pures strategies in I. In our first model, we create a vector of variables a i as the binary encoding of strategy i ∈ M, where a ij = 1 if the target j is covered by strategy i. The formulation is as follow: min a,λ, i∈M t∈T λ i it (21)

i∈M λ i = 1 (22) i∈M λ i a ij = c j j ∈ J (23) j∈J a ij ≤ m i ∈ M (24)
-it + q L t ≤ j∈J a ij pjt ≤ q U t + it t ∈ T, i ∈ M (25) a ∈ {0, 1} |J||M| , λ ∈ R |M| (26) ≥ 0 (27)
where λ i is the weight in the decomposition and t measures the violation of the fairness constraint for each strategy. The objective function (21) minimizes the wheigted violation of constraints (25). Equation ( 22) states that the weights must sum to 1 and Equation ( 23) ensures that the convex combination fits with c. constraints ( 24) and (26) define the pure strategy i. Expression (25) and the fact that t ≥ 0 defines the maximum deviation for each population type. This formulation is a non-convex mixed integer non-linear problem untractable even for small instances. We linearize this model by introducing variables γ ij λ i a ij and µ it = λ i it , and re-scaling constraints (25) by λ i . This mixed integer linear problem (MILP) is:

(DEC) min a,γ,µ i∈M t∈T µ it (28)
constraints ( 22), ( 24), ( 26)

i∈M γ ij = c j j ∈ J (29)
-µ it + λ i q L t ≤ j∈J γ ij pjt ≤ λ i q U t + µ it t ∈ T, i ∈ M (30) γ ij ≤ λ i j ∈ J, i ∈ M (31) γ ij ≤ a ij j ∈ J, i ∈ M (32) a ij + λ i -1 ≤ γ ij j ∈ J, i ∈ M (33) µ ≥ 0 ( 34
)
where constraints (31), ( 32) and (33) defines the product λ i a ij . We name this MILP (DEC). This formulation has two main drawbacks. First, we have to know a priori an upper-bound of the number of strategies to achieve the best decomposition. Secondly, the linear relaxation has always optimal value equal to zero. This means the formulation is a weak formulation and in consequence algorithms for solving MILP implemented in commercial optimization software might perform very poorly. We will show this issue in the next section.

In order to solve this problem in an efficient way, we propose the following column generation algorithm. Consider a set of feasible strategies. For a given strategy, it is straightforward to compute the violation of fairness constraints. We denote the violation of strategy a i by v(a i ). Then, we state the following linear optimisation problem (MP):

(MP) min λ i∈M λ i v(a i ) (35) s.t. i∈M λ i a i = c (36)
i∈M λ i = 1 (37)

λ i ≥ 0 ( 38
)
with the difference, that here a i and v(a i ) are parameters known for the problem. We denote α ∈ R |J| and β ∈ R the dual variables associated to constraints (36) and (37) respectively. The column generation algorithm works as follow: First choose a set of feasible set of pure strategies. They can be retrieved using Algorithm 1. Solve (MP) and get the dual variables. Then, compute the most negative reduced cost of a possible new strategy. This optimization problem, called column generator (CG), is stated as follow:

(CG) min

a t∈T v t (a) -j∈J α j a j -β (39) s
.t. j∈J a j ≤ m (40)
-v t (a) + q L t ≤ j∈J a j pjt ≤ q U t + v t (a) t ∈ T, i ∈ M (41) a ∈ {0, 1} (42)
where the objective function minimizes the reduced cost of the new strategy generated. If the new strategy generated by (CG) has a positive reduced cost, then the algorithm stops and the optimal solution of the master is the optimal solution of the whole problem. Otherwise, a new strategy is added to M with cost v(a) = t∈T v t (a). Example 1 continued. Now we decompose the coverage c pop using formulation (DEC). The decomposition is showed in Table 5. The distance decreases from 0.284 to 0.258. Note that the strategy that covers targets j 2 and j 4 (allocating 0 resources to the population type t 3 ) is not present anymore. 5. Decomposition of c pop using (DEC).

Some Extensions

Here we present two extensions to the discussion in this topic:

Multiple Labels: We consider the setting where each target has multiple labels representing different dimensions of analysis. It could be the case, where each target is characterised by race, religion, wealth, etc. Consider the set of attributes labels A = {L 1 , L 2 , . . . , L |A| }, , with its corresponding set of labels. Similarly to [START_REF] Zafar | Fairness constraints: A flexible approach for fair classification[END_REF], quotas (q L k , q U k ) on each label are stated, for each label k in the k-th attribute.

We would like to state an extension of Theorem 3 for this setting. If the labels satisfy that J k ⊆ J k or J k ∩ J k = ∅ for each k ∈ L k and k ∈ L k , then constraints forms a hierarchy. Then, Theorem 3 applies directly.

If it is not the case, we define L = L 1 × L 2 × . . . × L |A| which clearly induces a partition over J. If we set quotas as in ( 15) the result will be implementable under labels L, but is a relaxation of the original problem. We would need to find quotas (q L , q U ) for each , that produces the same constraints that quotas (q L k , q U k ). If we are able to find integer quotas satysfying this, then Theorem 3 applies again.

Penalizing violations with different weights: As in the Example 1, maybe it could be the case that decision makers would not prefer strategies that do not cover one type of population at all. The model (DEC) and the column generation, penalizes in the same way strategies that cover more than the quotas and the ones that cover less. This symmetry can be broken by introducing different variables to measure the violation, -and + , and replacing in the objective function (21) and constraints (25) by: min a,λ, -, + i∈M t∈T λ i (κ t

- it + + it ) (43)
-- it + q L t ≤ j∈J a ij pjt ≤ q U t + + it t ∈ T, i ∈ M (44)
where κ t > 1 is a parameter assigning more weight to under allocate protection to population type t. Linearization techniques and column generation can be straight applied as before.

Computational Experiments

In our computational results we investigate three questions: First, the impact in the defender expected utility by including fairness considerations in SSGs. Secondly, we test the computational performance of (DEC) and the column generation approach in order to get strategies close to be implementable. Finally, we test how much we win by decomposing coverage probabilities in almostimplementable pure strategies with our method instead with the Box method presented in Algorithm 1.

Experimental Setting: We test our methods in randomly generated games. For each n ∈ {20, 30, 40}, m ∈ {5, 10} and {1, 3} we generate 10 instances. Payoff matrices were generated uniformly such that D(j|p), A(j|u) ∼ U (0, 100) and D(j|u), A(j|p) ∼ U (-100, 0). For each game generated, an amount of 1000 "inhabitants" were allocated among the targets. Finally, we divided that population in |T | = 3 or |T | = 7 types, by running random partitions that sum the total population allocated in each target. With this setting, we aim to obtain defender's expected rewards comparables when population is divided in 3 or 7 types.

All experiments have been carried out using CPLEX 12.8 and Python 3.6, in a single thread on a server with a 3.40Ghz Intel i7 processor and 64 GB of memory.

Defender's Expected Utility: In order to measure the impact, we run the models of Population and Labels with different values of α ∈ 5%, 10%, 25%, 50%. In Figure 4, we show the average defender's expected utility in function of α, and separated by T . The model that consider labels in the targets is less restrictive, thus for a fixed α returns a bigger defender expected utility. Logically, both of them are upper bounded by the SSE returned by (HOM). Also, if we consider a more detailed description of the population (i.e. more types), the defender expected utility decreases. Finally, as we increase α, the defender expected utility for the three models converges to the SSE value. Efficiency of Decomposition Methods: For α = 25%, we solve all the instances with the model focused on population. We test the MILP (DEC) and the column generation. For the formulation (DEC), we use an upper bound on the number of strategies necessaries to decompose the coverages equal to U = |J| and U = 1.5|J|. We compare them with the output of Algorithm 1, using a lexicographic order.

We use a time limit of 600 seconds for the algorithms. Figure 5.a shows the average runtime in logarithmic scale. (DEC) considering U = |J|, respectively U = 1.5|J| hits the time limit in 97%, resp. 99%, of the instances. Algorithm 1 takes less than 1 mili-second in being performed. The Column Generation takes between 0.75 seconds and 4 seconds.

In figure 5.b, we compare the minimimum weighted violation and the weighted violation of the solution generated by Algorithm 1. Algorithm 1 returns consistently a solution that violates the more compared to the solution returned by the Column Generation. We note that the minimum weighted violation decreases as the size of the instance increase. From the practical point of view, the models that miminizes the weighted violation generates fairer allocations in reasonable time.

Case Study

In this section we study the performance of our models in a realistic instance, based on data retrieved from Ñuñoa, a municipality in Santiago of Chile. This instance consists in 1266 census blocks where the data considered is level of income (medium high, medium, low), the demand of policial resources DEM , the population in each census block and a measure of the criminal activity (amount of reported crime RC). The demand of policial resources is computed by Carabineros de Chile, the national police of Chile, as the amount of resources necessary to do general deployments, court orders, monitoring establishments and extraordinary services as in [START_REF] Bucarey | Shape and balance in police districting[END_REF]. To consider one target for each block is expensive to solve in terms of computational times, so we apply the clustering algorithm integrated in QGIS [START_REF] Team | Qgis geographic information system[END_REF], and we aggregate the data to reduce the size of the problem. We finally consider 250 targets represented in Figure 6.(a).

In this municipality there are three types of population in terms of level of education and income (strongly correlated). People with low income, medium income and high income, who receive 432US$, 750US$ and 1400$ per month. The geographical distribution of the income level is represented in Figure 6.(b). We consider two types of attackers: the first one, named k 1 , related to the criminal activity and the second one, named k 2 , related to the general activities that Carabineros perform. Both of them will appear with the same frequency, so π k1 = π k2 = 0.5. Payoffs are built proportionally to DEM and RC. The penalties for the attacker where used as a constant for each type of attacker (the punishment will be equal no matter where the attacker is caught). The specific parameters are shown in Table 6.

D k 1 (j|p) D k 1 (j|u) A k 1 (j|p) A k 1 (j|u) D k 2 (j|p) D k 2 (j|u) A k 2 (j|p) A k 2 (j|u) 0 -DEMj DEMj -100 0 -RCj RCj -300 Table 6. Payoff matrices built in the case study.

We test the model (HOM) (without any fairness consideration) and the models focused on population and labels. We deploy 120 homogoneus policial resources and we use the fairness parameter α = 0.1. All models took less than 1 minute to return the optimal coverage distribution. Results are shown in Figure 7. The model focused on the population allocates resources where the model without fairness consideration does not. Both, high and low income areas which, in model HOM, were not covered now present patrol assigments. The model focused on labels, on the other hand, presents a result with similar behaviour as HOM. Even when results are comparable, these models become a useful tool to ensure a fair distribution of police resources without having significant impact on its optimal allocation. This could be important for both police decision making and governmental policies.

Conclusions

In this paper, we have studied the impact of fairness constraints in SSG. We have presented two models: one imposing fairness constraints in a detailed description of the population, the other, imposing constraints in targets that are labeled. These models aim to allocate fair distribution of resources amongst the population, avoiding discrimination issues from a tactical point of view. We study our model on a realistic instance.

We have shown that imposing constraints with labels on the targets is implementable, meaning that each coverage distribution satisfying these constraints can be decomposed in pure strategies, all of them satisfying these constraints. This is not the case with the model with the detailed description of the population. In this case, we propose a MILP formulation to find the decomposition that is

Fig. 1 .

 1 Fig. 1. Example of Algorithm 1 with 5 targets and two homogeneous resources. (a) Coverage probability (b) Step 1 and 2. (c) Step 3.

Table 1 .

 1 6 34 -30 42 -9 15 -44 32 -29 49 -33 j2 11 -26 9 -16 47 -3 12 0 37 -39 16 -48 j3 2 -4 11 -39 9 -11 15 -16 25 -47 26 -5 j4 0 -35 3 -11 37 -6 22 -32 29 -23 21 -48 j5 9 -28 20 -48 47 -33 7 -14 42 -25 30 -40 Description of Example 1. a) Payoffs: Each row represent a target j. For each attacker type it is described D k = (D k (j|p), D k (j|u)) and D k = (A k (j|p), A k (j|u)). b) Description of the population detailed by types and labels.

3 .Fig. 2 .

 32 Fig. 2. (a) Comparison of the optimal coverage without any fairness consideration and the models focused on Labels and Population. (b) Defender expected utility in function of α.

Fig. 3 .

 3 Fig. 3. Representation of the bihierarchy. (a) H1 consisting in constraints (18) and singletons. (b) H2 consisting in constraints (19) and (20).

Fig. 4 .

 4 Fig. 4. Average defender's expected utility in function of α. (a) |T | = 3 (b) |T | = 7

Fig. 5 .

 5 Fig. 5. (a) Solution time for different methods to decompose coverage c in implementable strategies x. (b) Weighted violation of fairness constraints for each method.

Fig. 6 .

 6 Fig. 6. Data from Ñuñoa

Fig. 7 .

 7 Fig. 7. Coverage distribution in Ñuñoa.

Table 2 .

 2 Description of the solution in Example 1 in terms of population: (a) Coverage in the Strong Stackelberg equilibrium. (b) Coverage received by each type of population and the deviation from the fraction of population that they represent. (c) Coverage received by each labeled target. coverage on population type t is computed as c t = j∈J c j pjt t ∈T p jt

	Target j1	j2	j3	j4	j5	Type ct % resources Deviation
	cj 0.416 0.675 0.358 0.335 0.216		t1 0.552 27.6 %	97.143%
							t2 0.872	43.6%	21.111%
							t3 0.576	28.8%	-42.4%
		(a)					(b)
			Label c l % resources Deviation
			1	0.335 16.75 %	-16.25%
			2	1.033 51.65%	29.125%
			3	0.632	28.8%	-21.0%
						(c)

  .

		j1	j2	. . .	jn
	ω1	cω 1 j1 cω 1 j2	. . .	cω 1 jn
	ω2	cω 2 j1 cω 2 j2	. . .	cω 2 jn
	. . .	. . .	. . .	. . .	. . .
	ωm	cω m j1 cω m j2	. . .	cω m jn

http://www.ethnicity-facts-figures.service.gov.uk/

Authors have been partially supported by the Fonds de la Recherche Scientifique -FNRS under Grant(s) no PDR T0098.18.

closest to the set of fairness constraints. Also we propose a column generation method to solve this problem efficiently. Computational tests have shown that the column generation approach finds efficiently the best decomposition.