N Sghaier

S Geoffroy

M Prat
email: mprat@imft.fr

H Eloukabi

S Ben Nasrallah

Evaporation-driven growth of large crystallized salt structures in a porous medium

Subflorescence refers to crystallized salt structures that form inside a porous medium. We report a drying experiment revealing major development of subflorescence in the dry region of the porous medium away from the liquid zone. Using a combination of image analyses and numerical computations, we show that the growth is directly correlated to the evaporation flux distribution along the boundary of the growing salt structure. This indicates that the salt is transported into the domain occupied by the salt structure in the porous medium up to the structure periphery, where salt deposition takes place. This is confirmed when a growing salt structure encounters dry subflorescence formed earlier during the drying process. The dry subflorescence is reimbibed and resumes its growth. The analysis also suggests that the solution within the growing subflorescence is in equilibrium with the crystallized salt wall. These results shed light on the growth mechanisms of subflorescence, a phenomenon that can play a fundamental role in several important issues such as carbon dioxide sequestration or salt weathering.

I. INTRODUCTION

Salt crystallization induced by evaporation is an important issue in relation with several domains, such as soil physics [1], the underground sequestration of CO 2 [2], and salt weathering, i.e., the attack of rock and building materials by salt [START_REF] Goudie | Salt Weathering Hazards[END_REF]. It is now widely accepted that the crystallization process can induce sufficient stresses on pore walls to cause major damage to porous materials, e.g., [START_REF] Scherer | [END_REF]. As discussed in [5], the crystallized salt structures forming at the surface of a porous medium are referred to as efflorescence whereas the crystallized salt structures forming inside a porous material are referred to as subflorescence. Subflorescence is potentially the more damaging phenomenon since it corresponds to crystallization inside internal pores. While the analysis of the interplay between evaporation and efflorescence development has been the subject of several studies in recent years, e.g., [5][6][7][8][9][10], the understanding of subflorescence development within a porous material is much less advanced.

A better understanding is proposed in this paper, through a simple evaporation experiment allowing the growth of subflorescence to be directly observed.

II. EXPERIMENTAL SETUP

The two-dimensional transparent model porous medium considered in this study consisted of a monolayer of 1-mmdiameter glass beads. The beads were randomly thrown onto a glass plate coated with a hydrophobic sticky film. After removal of the excess beads, a second glass plate was put on top of the bead monolayer. In order to obtain a hydrophobic model, the glass beads and the glass plates were rendered hydrophobic by silanization, a process that leads to a contact angle of about 107°-110°regardless of the NaCl concentration [11]. Three lateral edges were sealed with RTV (a hydrophobic silicone elastomer) so that the vapor could escape only through the remaining open top edge of the model. The motivation for considering a hydrophobic model was twofold: (i) to obtain a particularly simple drying pattern (an almost flat liquid-gas front) and (ii) to limit the influence of secondary capillary liquid films, e.g., [12], that would make the analysis more complicated. However, subflorescence patterns and growth mechanisms similar to those reported here have been observed in hydrophilic models [13]. Thus the main results concerning the growth of subflorescence reported in this paper are general and not specific to hydrophobic systems. The model was initially fully saturated up to level z = H (see Fig. 1) with a NaCl aqueous solution of known initial concentration C 0 = 25 g/100 g solution. As can be seen from Fig. 1, an additional hydrophobic layer similar to the main micromodel was placed on top of the model. This 10-mm layer did not contain any liquid when the experiment started. It introduced an additional mass transfer resistance (water transport occurred by vapor diffusion in this layer), which limited the evaporation rate and prevented the formation of efflorescence (see [5,10] and references therein for the factors controlling the formation of efflorescence).

The 7-cm-high and 9.5-cm-wide porous medium was placed vertically in a small, transparent Plexiglas chamber at controlled temperature (22 ± 1 °C). The relative humidity in the chamber was stabilized using a saturated LiCl solution [relative humidity (RH) = 12%]. The relative humidity and temperature near the open edge of the model were recorded during the experiments using a sensor (Kimo VT300) connected to a PC and no noticeable evolution was observed. Variation of the model weight during drying was measured with a precision (0.01 g) scale, which was connected to the PC, while the phase distribution within the porous medium was recorded by a charge-coupled device (CCD) camera set on the side of the chamber. Images were taken every 25 min.

III. RESULTS

The overall duration of the experiment was 67 days and the mass loss of the sample over this duration was 2.5 g. This gives a low mean evaporation flux of j = 4.4 × 10 -6 kg/m 2 /s. As illustrated in Fig. 1, the drying process was characterized by an almost flat traveling front separating the region invaded by the gas phase from the fully saturated liquid zone. This is the pattern expected in a hydrophobic system, e.g., [14]. Moreover, the gravity effects helped to limit the formation of front heterogeneities here [15]. Although the system was hydrophobic, the contact angle (ß107°-110°) was not sufficiently high to prevent the transient trapping of liquid bridges, possibly interconnected [16], in the contacts between beads and between the beads and the plates in the region invaded by the gas phase. The transient existence of the liquid bridges was seen from the formation of small subflorescence spots occurring away from the main front when the liquid bridges evaporated. Interestingly, a larger subflorescence, referred to as the "sleeping" subflorescence in Fig. 1, also formed in the top region of the micromodel. This subflorescence developed in the liquid bridges region without being in direct contact with the main front. This suggests that liquid bridges can occasionally form an interconnected system carrying the dissolved salt between the main front and the growing subflorescence. After a while, the "sleeping" subflorescence stopped growing and dried. This can be inferred from the change in the subflorescence color from gray when wet to much whiter when dry [see Fig. 1(b) and compare the two main subflorescences]. The explanation is that the chain of liquid bridges feeding the subflorescence no longer forms a hydraulic connection between the subflorescence and the main liquid region after a while. Then the liquid bridges dry out and the subflorescence ceases to grow when all the salt originally present in the liquid bridges has precipitated. Additional details are given in Sec. IV.

Later, when the front had traveled over about two thirds of the micromodel height H , a little protuberance formed along the front as the result of the local pinning of the front on some local defects in the micromodel. This created a local overconcentration in salt, which gave rise to the formation of a major subflorescence, referred to as the wet subflorescence in Fig. 1 and as the main subflorescence in the rest of this paper. The main subflorescence remained attached to the liquid-saturated zone during its growth because of the long pinning of the front where the main subflorescence started growing.

As illustrated in Fig. 1, the first striking observation is that the main subflorescence progressively developed over a macroscopic length, comparable to the size, H , of the micromodel. The second somewhat surprising result is that the growth was not directed toward the liquid solution but in the opposite direction, i.e., towards the open edge of the network and into the gas-invaded region. This was in contrast with experiments in single capillary tubes [17], which have shown that the crystals usually grow into the liquid-saturated region and not in the gas region of the tube.

The color of the main subflorescence was gray when growing, which indicated that the subflorescence was wet. This suggests that the region occupied by the subflorescence was, in fact, porous. As a result, ions could be transported from the solution feeding the subflorescence at its base, within the subflorescence and up to its external boundary, where the salt precipitated and formed a new deposition layer. According to this picture, evaporation at the boundary of the subflorescence is the phenomenon inducing the flow within it and therefore the ion transport.

From a mean field -continuous picture, we should have u s = j/ρ along the subflorescence boundary, where u s is the normal component of the mean velocity in the liquid at the subflorescence boundary, j the local evaporation flux, and ρ the solution density. This equation expresses that the evaporation flux is balanced by a liquid flow (due to capillarity) toward the point considered on the subflorescence boundary.

This scenario is supported by what happens when the wet subflorescence reaches the "sleeping" subflorescence. This moment is referred to as the "encounter" of the two subflorescences. When this happens, the "sleeping" subflorescence wakes up and rapidly (in approximately 2 h) becomes gray, which indicates capillary invasion of the sleeping subflorescence by the solution transported through the wet subflorescence, and grows again (see Fig. 1).

To demonstrate that subflorescence growth is indeed driven by evaporation, the evaporation flux along the boundary of the wet (also referred to as active) subflorescence during its growth was computed as follows. First we determined the contour of the active subflorescence during its growth every 500 min (=every 8 h and 20 min). This involved a series of image processing steps using the MATLAB image processing tool box. The details of the image processing will be presented elsewhere. Such a contour is presented in Fig. 2(a).

Then we determined an indicator of local growth rates of the subflorescence as follows. Consider two successive The last step consisted in determining the evaporation flux along the contours. The method is illustrated in Fig. 3. Since the growth of the subflorescence was slow, the diffusion problem governing the water vapor concentration c within the gas phase in the micromodel could be considered as quasisteady. The problem to be solved thus read ∇ • (D∇c) = 0 in g where g is the domain occupied by the gas in the system formed by the micromodel and the additional top layer. On top of the domain (at z = H + h) we imposed c = 0, with a zero flux condition on the lateral edges and c = 1 on the subflorescence contour and the liquid-gas main front. For this type of computation, the exact value of D is not important. We took D = 1 in the region occupied by the gas and D = 0.01 where there were dry subflorescences in g since the pores were partially blocked by the salt crystals in these regions.

This problem was solved using the commercial simulation software COMSOL MULTIPHYSICS. The numerical solution allowed the local flux along the contour of active subflorescence to be evaluated as j = D∇c • n where n is the unit vector normal to the contour. This flux was computed where there were markers along contour number n + 1, which gave j k . Then the correlation between the subflorescence local growth rates r k and local evaporation fluxes j k was evaluated from the computation of the correlation coefficient,

χ (t) = N 1 [d k (t) -dk (t)][j k (t) -jk (t)] N 1 [d k (t) -dk (t)] 2 N 1 [j k (t) -jk (t)] 2 ,
where t is the time, and the overbar, ─, indicates the arithmetic average.

As can be seen from Fig. 4, the correlation coefficient was very high, i.e., close to 1, as expected, demonstrating that the growth of the subflorescence was directly driven by the evaporation process.

IV. DISCUSSION

While the general picture of subflorescence development can be understood from the results reported in the previous sections, two points need further discussion: (i) why the subflorescence growth is not stable, i.e., does not lead to a smooth interface but a ramified structure, and (ii) how ion transport is possible within the subflorescence without pore blocking due to salt precipitation or noticeable dissolution of the subflorescence away from the growing regions. Point (i) can be explained as follows. Due to random heterogeneities, local protrusions form at the evaporating interface, which lead to higher evaporation fluxes, because of the higher local vapor flux (due to the Laplacian). The accompanying salt precipitation creates a fine (water-wet) porous medium, which induces a capillarity-driven flow of water, and the creation of a positive feedback loop that leads to instability, hence the growth of subflorescence. This is essentially the same explanation as given in Ref. [18], where further details can be found, to explain the irregular shape of salt trees growing from a salt solution placed in a Petri dish. Point (ii) is perhaps more puzzling and requires more development.

Since the subflorescence is a porous structure, one can begin with the Darcy's scale equation governing the ion transport within the subflorescence, i.e. [19],

∂ρ εSC ∂t + ∇ • (ρ εU s SC) = ∇ • (ρ SεD * s ∇C) + f, ( 1
)
where ρ is the density of the solution, ε is the porous medium porosity, S is the liquid saturation in the region occupied by the subflorescence, D * s is the effective diffusion coefficient for the dissolved salt, C is the ion mass fraction, U s is the average interstitial velocity, and f is the dissolution-precipitation rate per unit volume of porous medium, which is expressed as

f = -ερ cr ∂S cr ∂t , ( 2
)
where ρ cr is the density of salt crystal and S cr the crystal saturation (volume fraction of pore space occupied by crystals within the subflorescence). Following [19], f is then expressed as

ερ cr ∂S cr ∂t = K cr (σ -1) g , ( 3
)
where σ is defined by

σ = m m sat . ( 4
)
In Eq. ( 4), m is the molality; m sat is the molality at saturation (m sat = 6.1 molar for NaCl). Expressed in term of mass fraction C, σ reads

σ = m m sat = C(1 -C sat ) (1 -C)C sat , ( 5
)
where the solubility C sat = 26.4% for NaCl. According to the data reported in Ref. [19], g = 1, K cr ≈ 0.41 kg/m 3 /s (precipitation) or 0.0128 kg/m 3 /s (dissolution) in the case of halite (halite is the crystallized form of NaCl), while the supersaturation σ is at most of the order of 1.6, according to [20]. The physical picture is clear. Evaporation at the periphery of efflorescence induces the flow of the solution within the subflorescence. This flow carries the ions toward the evaporative surface where the salt precipitates. As in any reactive flow, an important parameter is the Damköhler number comparing the reaction time scale and the (convective) transport time scale.

Here the Damköhler number can be expressed as

Da ≈ K cr [σ -1]L jC sat , ( 6
)
where j is the mean evaporation flux at the periphery of subflorescence and L is a characteristic length that can be taken as the mean height of subflorescence. Substituting numerical values representative of our experiment in Eq. ( 6) leads to Daâ1. The result is that the reaction is fast compared to transport. Under these circumstances, it can be assumed that the ion mass fraction in the subflorescence is essentially equal to the solubility mass fraction C sat . If the solution was not in equilibrium with the crystallized salt wall, precipitation should lead to pore blockage within the subflorescence or dissolution to shrinking of subflorescence. The fact that the precipitationdissolution phenomena within the subflorescence are negligible is consistent with the observations since the shape of subflorescence away from the active evaporation region does not change discernibly during the subflorescence growth. Unlike the classical case of ion transport in a porous medium, e.g., [8,21], where the key parameter controlling the ion distribution is the Peclet number, which characterizes the competition between convective and diffusive transport, ion transport in the porous subflorescence is thus controlled only by the convective transport (C ≈ C sat within the pores inside the subflorescence, which implies ∇C ≈ 0 and thus negligible diffusion transport).

At the bottom of the main subflorescence, the ion supply rate is

Q m = J s C sat , ( 7
)
where J s is the evaporation rate from subflorescence. Let V be the volume of porous domain occupied by the solution feeding the subflorescence at its base. A simple mass balance gives

ρ ε dV dt = -J s -J f , ( 8
)
where J f is the evaporation rate from the fraction of the boundary of this region not in contact with subflorescence.

Since the ion leaves this shrinking region only through the subflorescence, the ion mass balance for this region is

ρ ε d CV dt = -C sat J s , (9)
where C denotes the average ion mass fraction in this region.

From the visual inspection of the main subflorescence growth, one can make the rough approximation that the subflorescence growth is rapid compared to the change in the volume V, and so V cst during the subflorescence growth. Then Eq. ( 9) yields

ρ εV d C dt ≈ -C sat J s , (10)
which indicates that the average ion mass fraction in volume V decreases during the subflorescence growth. Thus the subflorescence growth tends to desalinate the solution feeding it at its base. This is consistent with the observations that show growth of new subflorescence structures [visible in the lower region of the micromodel in Fig. 1(d)] from shrinking volume V, starting only some time after the saturated region at the bottom of our micromodel has detached from the main subflorescence.

This also indicates that the main subflorescence stops growing not because of precipitation phenomena blocking the pores but because it stops being fed by the solution at its base when the subflorescence ceases to be in hydraulic contact with the solution as the result of invasion front receding due to evaporation.

The same lines of argument can be developed to explain how the subflorescence referred to as the "sleeping" subflorescence "goes to sleep." When the subflorescence ceases to be hydraulically connected to the solution, it forms a liquid-saturated porous structure in equilibrium with the solution contained in it (C ≈ C sat throughout the pore space). Then this subflorescence dries as, once again, the ions are transported to the evaporation active surface of the subflorescence where the salt precipitates to form new porous salt layers. This scenario is completely in line with observations which show preferential growth of this subflorescence at its upper part where evaporation fluxes are obviously greater and there is no pore blockage since the subflorescence resumes its growth when rewetted later as the result of the development of the main subflorescence.

From the above analysis, it is expected that the development of subflorescence must have an impact on the overall drying rate. The evaporation rate is roughly inversely proportional to the distance between the menisci inside the porous medium and the inlet of the porous medium. The result of subflorescence growth is to bring menisci much closer to the micromodel inlet and thus to reduce the mean distance between the menisci present in the system and the inlet. The conclusion is that the evaporation rate must increase with the development of subflorescence. This can be clearly seen from Fig. 5, which shows the evolution of the evaporated mass of water with time during the growth of the main subflorescence. The evaporation rate increases markedly during the growth of the main subflorescence and its subsequent encounter with the sleeping subflorescence, which causes the rewetting of the latter and subsequent significant subflorescence growth in the inlet region of the micromodel. Then the evaporation rate decreases consistently and significantly as the system formed by the union of the main subflorescence and the former sleeping subflorescence dries out (which corresponds to the plateau on the right in Fig. 5). The dry subflorescence forms a region of lower porosity and thus of less conductance for the diffusive transport of the vapor across the dry region of the model porous medium. This additional resistance combined with the fact that the gas-liquid front has receded into the system during the growth of the subflorescence explains the reduced evaporation rate after the drying of the main subflorescence. The final increase of the evaporation rate that can be discerned on the far right of the mass versus time curve in Fig. 5 is attributed to the growth of the bottom subflorescence clearly visible in the lower region of the micromodel in Fig. 1(d).

V. CONCLUSION

This study shows that development of subflorescence driven by evaporation can occur over macroscopic distances within a porous medium. The analysis indicates that most of the salt structure is in equilibrium with the solution traveling through it, which explains why the structure can continue to grow without significant pore blockage by precipitation or dissolution of the structure away from the evaporative surface. Here the precipitation kinetics does not play any role in the subflorescence growth kinetics, which is controlled by the transport of ions to the subflorescence evaporative surface and thus the evaporation flux distribution and its evolution at the periphery of the subflorescence. This is so because the precipitation time scale is much smaller than the ion transport time scale.

This study also suggests that a major development of subflorescence in a porous medium could be detected from the variation of the evaporation rate, which tends to increase as the result of subflorescence development. This is in contrast with classical drying kinetics under stationary external drying conditions, which never show an increase in the evaporation rate.

This study is of broader interest in relation with the widely studied domain of Laplacian growth [START_REF] Meakin | Fractals, Scaling and Growth Far from Equilibrium[END_REF]. Here, the growth results form an interesting coupling between evaporation (determined from the solution of the Laplace equation), capillarity (which induces the liquid flow in the subflorescence and thus the ion transport), ion transport, and crystallization. The interplay between these phenomena is likely to lead to various patterns and this clearly deserves further studies. This is supported by our recent study on efflorescence [6,10]. From the results reported in Refs. [6,10], we surmise that major subflorescence developments such as the ones reported in the present paper can be expected for sufficiently slow evaporation rates and systems with relatively large pores. When these conditions are not met, subflorescence development could be inhibited by limited solution transport capacity toward its boundary due to viscous effects and/or possibly to salt precipitation (when the evaporation rate is sufficiently high for the ion transport time scale not to be considered large compared with the precipitation time scale; it can also be hypothesized that the probability of blocking pores would increase with decreasing pore size), leading to possibly different subflorescence patterns. This is supported by the results reported in Refs. [START_REF] Goncalves | [END_REF]24] showing decreasing rates of drying for salt-laden real materials. The next step is therefore to investigate the influence of evaporation rate and pore size on subflorescence growth.

FIG. 1 .

 1 FIG. 1. Images of model porous medium with subflorescence at different stages of development. The white dashed line approximately marks the separation line between the liquid-saturated zone and the gas-invaded region. Panel (d) shows the system near the end of the drying process. Note the white color of the main subflorescence, indicating that it is now dry, and the major development of new subflorescence at the bottom of the system.

FIG. 2 .

 2 FIG. 2. (Color online) (a) Example of contour of active subflorescence extracted from the raw image after a series of image processing steps; (b) sketch of the method used for determining the local growth rate of subflorescence; (c) visualization of the successive salt deposition layers. Each layer corresponds to the growth occurring during 500 min. The contour of the former "sleeping" dry subflorescence is clearly visible in the upper region.

FIG. 3 .

 3 FIG. 3. Two-dimensional finite element computation of evaporation flux along the contour of subflorescence. (a) Finite element mesh; (b) computed water vapor concentration isolines.

FIG. 4 .

 4 FIG. 4. Correlation coefficient χ between subflorescence local growth rates and local evaporation fluxes at the subflorescence boundary. The drop in the correlation coefficient on the right-hand side of the plot corresponds to the encounter between the two subflorescences (this is so notably because we wrongly imposed c = 1 along the contour of the former dry subflorescence even when this subflorescence was not yet reimbibed). The insets show two examples of rescaled distribution of growth rate (figure on the left in each pair of figures) and evaporation flux (figure on the right) along the contour, illustrating the excellent correlation between local growth rates and local fluxes.

FIG. 5 .

 5 FIG. 5. (Color online) Evaporated mass of water versus time (inset). The main figure shows the variation of the evaporated mass of water during the period of growth of the main subflorescence. The dashed lines are guides to the eye materializing the change in the evaporation rate (slope of lines) when the main subflorescence develops and dries out (after union with the sleeping subflorescence).


Ali El

Olivier Brun
email: brun@laas.fr

Offline Re-Optimization of VNF Placement Decisions for Existing Network Slices -A Game-Theoretic Algorithm

Keywords: Network Function Virtualization, Virtual Network Function -Forwarding Graph, Slice as a Service, Offline Allocation, Game Theory, Optimization

Network Function Virtualization (NFV) simplifies the deployment of network services by leveraging virtualization technologies to make the management of network functions more flexible and cost efficient. The deployment of these services requires the allocation of Virtual Network Function -Forwarding Graph (VNF-FG), which implies placing and chaining VNFs according to the requests of VNF-FGs. In this paper, we consider the offline allocation of VNF-FG problem to improve resource utilization and reduce total costs. We focus on how VNF-FG demands are routed so as to optimize resource utilization without adding capacity to the infrastructure. Given a non-linear cost function associated to each network resource, we formulate the problem as a non-linear single-path routing problem in an extended graph. Then, we propose to adapt a single-path routing heuristic algorithm inspired from game theory to solve it. We show that this algorithm converges and establishes its approximation ratio in a number of cases. Experimental results obtained for different network topologies and different cost functions show that this algorithm provides very good quality solutions with a substantially lower computing times compared to the optimal solution.

Introduction

Network Function Virtualization (NFV) opens up the possibility to replace hardwired middleboxes implementing specific network functions (e.g., firewalls or Web proxies) by Virtual Network Functions (VNFs) running in datacentres and operated as cloud services. In addition to reducing the cost and complexity of the network, this approach enables to deploy and compose VNFs so as to create network slices on demand. Network slices are isolated logical networks coexisting simultaneously on the same Physical Substrate Network (PSN), each one being tailored to the needs of a specific application. The VNFs composing the network slice as well as the virtual links established between them are described by a VNF Forwarding Graph (VNF-FG). A VNF-FG, according to ETSI, is a directed graph of an ordered set of VNFs that compose an end-to-end service [START_REF] Nfv | Network Functions Virtualisation (NFV): Architectural framework[END_REF].

The logical service graph described by a VNF-FG is decoupled from the network and has to be mapped to existing PSN capabilities [START_REF] Herrera | Resource allocation in nfv: A comprehensive survey[END_REF]. This implies not only the fulfillment of slice's requirements in terms of Quality of Service (QoS), but also considering the constraints of the underlying infrastructure. Two types of decisions have to be made: (a) where to run the VNFs, and (b) how to interconnect them in the physical network, taking into account specific VNF ordering requirements. They are usually made so as to optimize a certain performance objective, which may be related for instance to the resource utilization of the PSN, while satisfying a number of technical constraints (e.g., some VNFs should be placed in the vicinity of end users). This problem, which is known as the VNF Forwarding Graph Embedding (VNF-FGE), is an extension of the Virtual Network Embedding problem [START_REF] Fischer | Virtual network embedding: A survey[END_REF][START_REF] Laghrissi | A survey on the placement of virtual resources and virtual network functions[END_REF] which is known to be NP-hard [START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF]. The VNF-FGE problem has been studied in the online context, where requests to create network slices are submitted one after the other [START_REF] Soualah | A green vnf-fg embedding algorithm[END_REF][START_REF] Mechtri | Vnf placement and chaining in distributed cloud[END_REF][START_REF] Carpio | Replication of virtual network functions: Optimizing link utilization and resource costs[END_REF][START_REF] Quang | Single and multi-domain adaptive allocation algorithms for vnf forwarding graph embedding[END_REF][START_REF] Amine | Shortening the deployment time of SFC by adaptively querying resource providers[END_REF], and in the offline context where multiple VNF-FG have to be embedded in the physical network at the same time [START_REF] Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF][START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF][START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in nfv-enabled inter-dc network[END_REF][START_REF] Guo | Cost-aware placement and chaining of service function chain with vnf instance sharing[END_REF].

In this work, we consider a physical infrastructure in which transmission and computing resources are available. These resources have been used to set up a number of logical networks. As the requests for new logical networks arrive one after the other in time, without knowing neither the future requests nor the logical networks that will be removed in the future, the individual placement and routing decisions made are necessarily sub-optimal in hindsight, which may result in poor resource utilisation. It may therefore be necessary to re-optimise the mapping of existing logical networks onto the physical infrastructure. It should allow for a better load distribution on existing resources, and thus increase the capacity of the physical infrastructure to accommodate new logical networks without having to add additional capacity.

Despite the previous efforts on the VNF-FG embedding problem, little considered the re-optimization of resources utilization on the physical infrastructure with non-linear cost functions. In [START_REF] Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF], the authors aimed at minimizing the number of VNF instances mapped on the infrastructure to improve resource utilization. However, the work ignored the routing cost from connecting the VNFs of the VNF-FG request. Inspired by the previous work to share the same VNF by multiple VNF-FGs, the work in [START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF] and [START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in nfv-enabled inter-dc network[END_REF] studied the offline VNF-FGE problem to reduce both deployment and forwarding traffic costs. In [START_REF] Tastevin | A graph approach to placement of service functions chains[END_REF], the authors proposed a graph-based heuristic inspired from the Viterbi algorithm to reduce the deployment cost of requested VNF-FG by reducing the number of PoPs hosting the VNF instances. However, this solution results in an unbalanced infrastructure resource utilization. In [START_REF] Zhong | Cost-aware service function chaining with reliability guarantees in nfv-enabled inter-dc network[END_REF], the authors focused at minimizing the cost of routing and deploying the requested VNFs, while guaranteeing a certain level of reliability. The work however, does not consider an efficient use of resources for future requests. The closest to our work is [START_REF] Guo | Cost-aware placement and chaining of service function chain with vnf instance sharing[END_REF] that focuses on minimizing both resource utilization and total placement costs. Despite the similarities with our work, the model considered is different from ours. First, the authors considered an infrastructure with only PoP nodes that are used to host VNFs. Hence, their work do not consider the routing problem. Second, the authors considered the possibility to deploy only one type of VNF in a PoP. Moreover, none of the previous works considered non-linear cost functions which are essential if one wishes to express that the cost per unit capacity of a resource grows with its utilization rate.

In the setting that we consider, we are given a physical network whose links have known transmission capacities and in which some nodes provide known computing resources. A non-linear cost function is associated to each network resource and allows to measure its congestion level. A set of traffic demands flow through the network, each traffic demand being characterised by its source, destination, volume and the sequence of VNFs it must flow through. The volume of a traffic demand may change after each VNF, due to, for example, encryption or data aggregation operations. An initial solution specifies, for each traffic demand, on which compute nodes the VNFs are executed and the path taken by the demand in the network between these VNFs. We show that the VNF placement and chaining problem can be formulated as a non-linear single-path routing problem in an extended graph. To solve this problem, we adapt a single-path routing algorithm inspired from Game Theory which was proposed in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. We show that this algorithm converges and establishes its approximation ratio in a number of cases. Experimental results obtained for different network topologies and different cost functions show that this algorithm provides very good quality solutions with a rather modest computation time.

The paper is organized as follows. We describe the problem addressed in this paper and introduce our main notations in Section 2. In Section 2.1, we introduce the concept of expanded network. The problem is then formulated as a non-linear single path routing problem in the expanded network in Section 2.2. We present the proposed algorithm in Section 3. Section 4 is devoted to the performance evaluation of this algorithm using diverse network topologies and cost functions. Finally, some conclusions are drawn in Section 5.

Problem Statement

We are given a physical network represented by a directed graph G = (V, E), where V is the set of nodes and E is the set of links. A subset D ⊆ V of the nodes have the required compute and storage resources to host VNF. These nodes will be called the function nodes in the following. They are used to execute one or more VNFs out of a set F of K VNFs, where function f ∈ F is available at function nodes v ∈ D(f ) ⊆ D. A function node corresponds to a logical entity representing reserved computing and storage capacity in a data centre. Multiple VNFs may be deployed at the same function node, and a VNF may be replicated at different function nodes. The transmission links of the network and the function nodes are the critical resources in our problem and in the following we shall denote by R = E ∪ D the set of these resources.

We are also given a set W of traffic demands. Each traffic demand w ∈ W is characterized by its source node sw, its destination node tw and the sequence Fw = f w 1 , f w 2 , . . . , f w Kw of network functions it must flow through in that prescribed order. The sequence is called the Service Chain associated to the traffic demand. The processing path of demand w may therefore be decomposed in different stages, where stage 0 corresponds to the transmission of the original traffic volume from the source to VNF f w 1 (excluded), stage 1 represents the processing at VNF f w 1 followed by the transmission to VNF f w 2 , etc. The volume of a traffic demand may change at each stage due to, e.g., packet encryption/decryption or to data aggregation Capacity required by demand w from resource r at stage k.

operations. Moreover, the processing capacity required per unit of traffic may be different for one VNF to another. We therefore define b w k,r as the capacity required from resource r for processing demand w at stage k. If r is a transmission link e ∈ E, then this capacity just represents the volume of traffic send from VNF f w k to VNF f w k+1 for 0 < k < K(w), and the volume of traffic from sw (resp. VNF f w K(w) ) to VNF f w 1 (resp. tw) for k = 0 (resp. = K(w)). This volume of traffic should be the same for all links e and is expressed in bit/s. Similarly, if r is a function node v ∈ D(f w k ), the capacity b w k,r represents the (fractional) number of cores required for processing the traffic of demand w by VNF f w k . Our notations are summarized in Table 1.

Figure 1 provides a simple example of the setting considered in the paper, in which there are two traffic demands. The service chain for the traffic demand s1 -t1 is F1 = (f1, f2), whereas the service chain for the traffic demand s2 -t2 is F2 = (f1). The volume of traffic demand s1 -t1 is b 1 0 = 20 units of traffic and the volume of traffic demand s2 -t2 is 10 units. To simplify, we assume that the volumes do not change along the paths from the sources to the destinations. We also assume that VNF f1 requires 1.5 cores per unit of traffic, whereas VNF f2 requires 1.0 core per unit of traffic. The capacity of all links are taken to be ce = 60 units of traffic, and the processing capacity of the function nodes are cD = 70 and cE = 40. The figure illustrates two possible deployment scenarios. In the first scenario (see Fig. 1a), all instances of all VNFs are deployed at node D, whereas in the second scenario (see Fig. 1b) the instances of VNF 1 required for processing the traffic demand s1 -t1 are executed at node E. We analyze below each scenario, assuming that the cost function for resource r is yr cr -yr :

(a) Scenario 1: In this scenario, the utilization rates of the network links are as follows: 1 3 for one link, 1 6 for two links and 1 2 for four links. This yields a total link cost of 1× 1/3 1-1/3 +2× 1/6 1-1/6 +4× 1/2 1-1/2 = 4.9. However, the total capacity required at node D is 65 cores ((20 + 10) × 1.5 for VNF f1, and 20 × 1 for VNF f2), yielding a total cost for function nodes equals to 65 70-65 = 13. Hence the total network cost for this deployment scenario is 17.9.

(a) Scenario 2: In this scenario, the utilization rates of the network links are as follows: 1 3 for seven links, 1 6 for three links and 1 2 for three other links. The total link cost is therefore 7

× 1/3 1-1/3 + 3 × 1/6 1-1/6 + 3 × 1/2 1-1/2 = 7.1.
The computing capacity required at node D is 35 cores: 10 × 1.5 for VNF f1, and 20 × 1 for VNF f2. The computing capacity required at node E is 20 × 1.5 = 30 cores. The total computing cost of function nodes is hence 35 70-35 + 30 40-30 = 4. We conclude that the total network cost in this scenario is 11.1, which shows that it achieves a better resource utilization.

The problem at hands amounts to choosing a single path through the network for each traffic demand such that the function nodes are visited in the prescribed order and so that a certain network cost function is minimized. We shall shortly discuss the network cost function to be optimized, but we would first like to emphasize that a first difficulty is related to the modelling of precedence constraints on the order in which the function nodes must be visited by each traffic request. The key observation is that the overall path from sw to tw for traffic demand w can be viewed as a collection of segments, the first segment connecting sw to one of the function nodes hosting f w 1 through some intermediate nodes, the second segment connecting the end node of the first segment with one of the function nodes hosting f w 2 , etc, until the destination node tw is reached. These segments are connected by decision points where we choose which nodes will compute the functions in Fw. As observed in [START_REF] Nguyen | Optimizing resource allocation in infrastructure networks based on network function virtualization[END_REF], it is possible to construct an expanded network to model both the segment construction and the function selection at the nodes. We shall discuss the construction of this expanded network in Section 2.1. In Section 2.2, we formulate the single-path routing problem considered in this paper using the concept of expanded network.

s 2 b 1 0 = 20 b 2 0 = 10 f 1 s 1 t 1 C B D E f 1 , f 2 t 2 A
(a) Deployment scenario 1 in which all VNFs are executed at node D.

s 2 b 1 0 = 20 b 2 0 = 10 f 1 s 1 t 1 C B D E f 1 , f 2 t 2 A
(b) Deployment scenario 2 in which instances of VNF 1 are executed at node E for the traffic demand s1 -t1.

Figure 1: A simple example network , in which we assume that VNF f 1 may be executed at nodes D and E, whereas VNF f 2 is only available at node D.

f1 f1 f1, f2 s (0) 1 s (0) 2 t (0) 1 t (0) 2 A (0) C (0) B (0) D (0) E (0) s (1) 1 s (1) 2 t (1) 1 t (1) 2 A (1) C (1)
B ( 1)

D (1) E (1) s (2) 1 s (2) 2 t (2) 1 t (2) 2 A (2) C (2) B (2) D (2) E (2) f1, f2

Expanded Network

The expanded network G * is made up of K layers. Each layer G (k) = (V (k) , E (k) ), k = 0, . . . , K, is a copy of of the original network G = (V, E). We shall use the notation v (k) (resp. e (k) ) to refer to node v ∈ V (resp. link e ∈ E) at layer k. We use each layer to find one segment of the total path for each traffic demand w. Consecutive layers are connected by artificial edges which are used to model the function selection at the nodes. More precisely, for each function node v ∈ D and each layer k = 0, . . . , K -1, there is an artificial directed edge from node v (k) to node v (k+1) . The interpretation is that a path for traffic demand w goes through this artificial link if and only if VNF f w k is executed at node v ∈ D(f w k ). In the following, we shall use the notation ℓ (k) v to refer to the artificial edge (v (k) , v (k+1) ) and we shall denote the set of artificial edges between layers k = 0, . . . , K -1 and k + 1 by

L (k) = ℓ (k) v : v ∈ D . In summary, the set of nodes of the expanded network G * is V * = K k=0 V (k) ,
whereas its set of edges is

E * = K k=0 E (k) ∪ K-1 k=0 L (k) .
The construction of the expanded network is shown in Figure 2 for the example network of Figure 1. In the following, we say that an artificial edge ℓ

(k) v ∈ L (k) is feasible for demand w if and only if v ∈ D(f w k )
, that is, it connects two consecutive copies of a node v hosting function f w k . A feasible path πw ⊂ E * for the traffic demand w ∈ W is then defined as a path from node s (0)

w to node t (Kw ) w
in the expanded network such that all artificial edges in the path are feasible for demand w. Figure 2 shows two feasible paths for the demand from node s1 to node t1. Note that the artificial edges appearing in a path specify the function nodes at which the VNFs are executed, while the other edges in the path specify the path segments used to reach these function nodes. The feasibility of a path ensures that the network functions are executed in the prescribed order.

Given a set π ⊂ E * of edges in the expanded network and a link e ∈ E of the original network, we define the constant δ e,k π as 1 if e (k) ∈ π, and 0 otherwise. Similarly, for any set π ⊆ E * and any function node v ∈ D, we define the constant δ v,k π as 1 if the artificial edge ℓ (k) v ∈ π, and 0 otherwise. In the following, the notation δ r,k π will refer to one of the above-defined constants, depending on whether resource r ∈ R is a transmission link e or a function node v. We shall say for short that resource r appears in path π at layer k of the expanded network and write (r, k) ∈ π when δ r,k π = 1. Similarly, we write r ∈ π if δ r,k π = 1 for some k.

Single-Path Routing Problem

In the following, we shall assume that we are given a set Πw of candidate feasible paths through the expanded network for each demand w ∈ W. A feasible solution for the joint routing and VNF placement problem is then defined as a vector π = (πw) ∈ Π, where πw is the path assigned to traffic demand w and Π = w∈W Πw. Given such a feasible solution π, the traffic sent by demand w on resource r ∈ R of the original network is:

y w r (πw) = K k=0 δ r,k πw b w k,r , (1)
from which it follows that the traffic flowing on resource r is yr(π) = w∈W y w r (πw). We assume that to each resource r ∈ R of the physical network is associated a non-decreasing function ϕr : R+ → R+ and that the cost of this resource has the form yrϕr(yr). The function ϕr() may be interpreted as the cost per unit of capacity of resource r, and it may depend on the total traffic flowing on that resource. A typical example is the Kleinrock function ϕr(yr) = 1/(cr -yr), where cr is a constant representing the (transmission or processing) capacity of resource r. This function assumes that the total traffic on a resource is smaller than its capacity and that the cost per unit of capacity grows unboundedly as the former approaches the latter. Other examples are the linear function ϕr(yr) = yr/cr or the quadratic function ϕr(yr) = y 2 r /c 2 r , which are valid even for yr ≥ cr. In the above examples, the cost per unit capacity grows as the traffic yr on the resource increases, in contrast with linear models in which ϕr is assumed to be a constant function. We emphasize that the functions ϕr() may differ for different types of resources or even from one resource to the other.

The goal is to find a feasible solution π ∈ Π that minimizes the network cost F (π) = r∈R yr(π) ϕr (yr(π)). Formally, the problem is as follows: minimize F (π) = r∈R yr(π) ϕr (yr(π)) (OPT) subject to:

πw ∈ Πw, w ∈ W. (2)
It is worth mentioning that problem (OPT) can be cast as a 0-1 mathematical programming problem by introducing the following binary variables:

xw,π = 1 if demand w is routed along path π in G * , 0 otherwise, (3)
so that it becomes equivalent to: minimize r∈R yr ϕr(yr)

subject to:

yr = w∈W π∈Πw K k=0 δ r,k π b w k,r xw,π, for all r ∈ R, (4)
π∈Πw xw,π = 1, for all w ∈ W, (5)
xw,π ∈ {0, 1}, for all π ∈ Πw and w ∈ W.

Except when the functions ϕr() are constant functions, the above problem belongs to the class of nonlinear mathematical programs with integer variables, a class of mathematical programs which are known to be extremely hard to solve. This motivates the development of an efficient approximation method for solving problem (OPT).

Penalized Best-Response Algorithm

As discussed in Section 2.2, the joint routing and VNF placement problem can be cast as a non-linear single path routing problem in the expanded network. The problem differs however significantly from traditional single-path routing problems as on one hand the same resource appears at different layers of the expanded network and on the other hand the volume of a traffic demand may change from one layer to the other. Nevertheless the heuristic algorithm that we propose in this section is directly inspired from an algorithm proposed in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF] for solving such problems. The idea of the algorithm is to view the traffic demands as the players of a non-cooperative game in which each each player independently optimizes its own objective function. Starting from an arbitrary initial feasible solution, the algorithm then mimics the best-response dynamics of the game, that is, the players take turns in some order to adapt their strategy based on the most recent known strategy of the others. The player objective functions are designed in such a way that (a) the best-response dynamics converges to a Nash equilibrium of the game in a finite number of steps, (b) all optimal solutions of problem (OPT) are Nash equilibria of the game, and (c) an upper bound on the approximation ratio of the algorithm can be proven in some cases.

Let us think of the traffic demands as the players of the game. The strategy of player w ∈ W is the path πw it chooses in the set Πw, and a strategy profile of the game is a feasible solution to problem (OPT), that is, a vector π ∈ Π. Given the strategy of the other players π-w = (πu) u̸ =w , we shall assume that the player w seeks to solve the following problem:

minimizeπ∈Π w cw(π, π-w) = fw(π, π-w) + pw(π, π-w), (OPT-w)
where the value fw(π, π-w) associated to path π by player w reflects the cost of this path, whereas the term pw(π, π-w) is a penalty term measuring the impact of player w's choice on other players. More precisely, we shall assume that the selfish objective function of player w is:

fw(π, π-w) = r∈π y w r (π) ϕr y -w r + y w r (π) , (7)
where y -w r = u̸ =w y u r (πu) is the total traffic sent over resource r by all the other players u ̸ = w. The term fw(π, π-w) thus represents the cost of path π for player w, given the fixed strategies π-w of the other players.

In order to define the penalty term pw(π, π-w), assume that player w chooses path π ∈ Πw and consider a resource r appearing in the path π at some layer. Then, for all players u ̸ = w, the cost of this resource increases by:

y u r (πu) ϕr y -w r + y w r (π) -ϕr y -w r .
As a consequence, we define the penalty term pw(π, π-w) as follows:

pw(π, π-w) = With ( 7) and ( 8), the objective function c w (π, π -w ) = f w (π, π -w ) + p w (π, π -w ) of player w is perfectly defined. Given the strategies π -w of the other players, the path π minimizing c w (π, π -w ) in problem (OPT-w) is known as the best response of player w. The pseudocode of our heuristic is given in Algorithm 1. The algorithm starts from an initial feasible solution π (0) . At each iteration n, the players update their strategies in a given order by computing their best responses to the strategies of the others (lines 3 -5). Note that a player w deviates from its strategy π (n) w at iteration n to a new strategy π ′ if and only if c w π ′ , π (n) -w < c w π (n) . The algorithm stops when no player can decrease its cost by unilaterally deviating from its strategy, that is, π (n+1) = π (n) , which means that a Nash equilibrium has been reached (line 7).

Algorithm 1 Penalized best-response

Require:

π (0) 1: n ← 0 2: repeat 3:
for w ∈ W do end for 6:

n ← n + 1 7: until π (n+1) ̸ = π (n) . 8: return π (n)
Theorem 1 below states the main properties of Algorithm 1. The proofs of properties (a), (b) and (c) are straightforward adaptations of those of similar results in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. The proof of property (d) is notably more involved, though it follows the same lines as the proof of Theorem 3 in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]. Note that property (c) follows from property (d) by taking d = 0. (c) It computes an optimal solution of problem (OPT) for linear resource costs, that is, when the functions ϕ r () are constant functions ϕ r (y) = a r .

(d) Its approximation ratio is 2

1 d+1 -1 -(d+1)
when the functions ϕ r () are polynomials of degree

d, that is, ϕ r (y) = d n=0 a r,n y n . Proof. See Appendix A.
Algorithm 1 is optimal for linear resource costs and returns a solution whose cost is at most 5.83× the cost of an optimal solution for quadratic resource costs. Unfortunately, the approximation ratio of the algorithm for polynomial resource costs degrades quickly with the degree of the polynomial (it is already equals to 56.9 for d = 2). Our numerical results suggest however that the heuristic provides close-to-optimal solutions in most cases.

The proof of Theorem 1 shows that the game that we consider is a potential game, which implies that the worst-case complexity of Algorithm 1 is exponential in the number of traffic demands [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF]. However, as we will see in Section 4 below, its convergence in practice is much faster than suggested by this worst-case result.

Performance Evaluation

In this section, we experimentally evaluate the performance of the penalized best-response algorithm. We first describe the different types of objective functions used to solve the penalized best-response algorithm and the procedure for generating random instances in sections 4.1 and 4.2. Then, we present in section 4.3 the numerical results for the different topologies evaluated under the different types of objective functions.

Cost Functions

Given the cost per-unit of capacity of resource r function ϕ r , we denote by Φ r = y r ϕ r (y r ) the cost function of resource r. We consider four types of cost functions: linear, piece-wise linear, quadratic and the Kleinrock Function (M/M/1).

1. Linear cost function: the per-unit cost does not depend on the utilization rate of the resource r (ϕ r (y r ) = a, where a is a constant). Here, the objective function is linear with the amount of traffic which is solved using Gurobi 9.0 as the ILP solver. We note that in this case the solution may over utilize the resources capacity. We choose a = 5, and we have:

Φ r = 5 × y r (9)
2. Piece-Wise Linear cost function: we consider the following increasing per-unit cost function:

ϕ r (y r ) =      3 if 0 ≤ yr cr ≤ 1 4 , 5 -1 2 cr yr if 1 4 < yr cr ≤ 3 4 , 10 -17 4 cr yr if 3 4 < yr cr . ( 10
)
where c r is the capacity of resource r. This function expresses that it is cheap to utilize a resource with a small utilization rate, whereas, as the load approaches c r , it becomes more expensive. Problem (OTP) can then be formulated as a Mixed-Integer Linear Problem and solved using Gurobi 9.0 as the MILP solver:

minimize r∈R Φ r
subject to:

y r = w∈W π∈Πw K k=0 δ r,k π b w k,r
x w,π , for all r ∈ R, (

) π∈Πw x w,π = 1, for all w ∈ W, 11
Φ r ≥ ay r -bc r , for all r ∈ R and (a, b) ∈ C,

x w,π ∈ {0, 1}, for all π ∈ Π w and w ∈ W. ( 14)

y r ≥ 0, Φ r ≥ 0, for all r ∈ R, (15)
where C = {(3, 0), (5, 1 2 ), (10, 17 4 )}. 3. Quadratic cost function: The per-unit cost function is linear (ϕ r (y r ) = ay r ). This results in a quadratic objective function which is solved using Gurobi 9.0 as the Quadratic Integer Program solver. Similar to the previous function, the more is the utilization rate of a resource, the more expensive it becomes. Φ r = y 2 r /c 2 r (16)

4. Kleinrock cost function (M/M/1): with this function, it becomes hard to find the optimal solution. However, it guarantees that the utilization of a resource r does not exceed its capacity c r . Φ r = y r c r -y r [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF] We formulate problem (OPT) as a Bilinear Program, and we solve it with Gurobi 9.0 Bilinear solver.

y r = w∈W π∈Πw K k=0 δ r,k π b w k,r x w,π , for all r ∈ R, ( 18
) π∈Πw x w,π = 1, for all w ∈ W, (19)
y r -Φ r c r + Φ r y = 0, for all r ∈ R, (20)
x w,π ∈ {0, 1}, for all π ∈ Π w and w ∈ W.

(21)

y r ≥ 0, Φ r ≥ 0, for all r ∈ R, (22)

Generation of Random Instances

In order to evaluate the performance of the penalized best-response algorithm, we use five network topologies (see Table 2) collected from the IEEE literature and from The Internet Zoo Topology [START_REF] Knight | The internet topology zoo[END_REF]. As illustrated in Figure 3 for the NSFNET topology, we choose three source nodes (s 1 , s 2 and s 3 ) and three destination nodes (t 1 , t 2 and t 3 ) in each topology. The original topology is also expanded by adding three function nodes v 1 , v 2 and v 3 and connecting them to the other nodes. We consider three VNFs f 1 , f 2 and f 3 and assume that each one requires one core per unit of traffic. Function node v 1 can host VNFs f 1 and f 2 , whereas function node v 2 (resp. v 3 ) can host VNFs f 1 and f 3 (resp. f 2 and f 3 ). As shown in Figure 3, each VNF is thus available at two different locations.

We then generate 100 random instances for each network topology. A random instance is composed of 25 traffic demands, which are randomly generated as follows. The source and destination nodes of a traffic demand are chosen randomly in the sets {s 1 , s 2 , s 3 } and {t 1 , t 2 , t 3 }, respectively, according to uniform distributions. The volume of a traffic demand is drawn from a uniform distribution in the interval [START_REF] Nfv | Network Functions Virtualisation (NFV): Architectural framework[END_REF][START_REF] Amaldi | On the computational complexity of the virtual network embedding problem[END_REF], and, to simplify, we assume that it does not change along the path from the source to the destination. Finally, the processing path of a traffic demand in chosen uniformly at random in the set

{(f 1 , f 2 ), (f 1 , f 3 ), (f 2 , f 3 ), (f 1 , f 2 , f 3 )}.
Each traffic demand must thus go through two or three VNFs, which are available at multiple locations. The possible paths for a traffic demand w are generated as follows. We first compute two possible paths between the source node s w and each location where the first VNF f w 1 is available by solving a 2-shortest path problem (assuming unit weights for the edges of the expanded network). We then apply the same procedure for computing two possible paths between each possible location of the VNF f w 1 and each possible location of f w 2 , etc. The final set of candidate paths Π w is obtained by connecting the different segments forming the overall path from s w to t w .

The initial solution is obtained by choosing randomly a path for each traffic demand. We assume that all links and all function nodes have the same capacity. This capacity is computed for each random instance so as to obtain a network congestion rate of 0.83.

We report in Table 3 the number of candidate paths for a traffic demand |Π w | calculated for each studied topology. The max (min, resp.) value is obtained when the size of the processing path is three (two, resp.) resulting in |Π w | = 2 7 (|Π w | = 2 5 , resp.). We remind the reader that we have k = 2 candidate paths for each VNF option. The average value is obtained over 100 problem instances where traffic demands have a probability of 3/4 (1/4, resp.) to choose a processing path size of two (three, resp.). In the case of CESNET topology and under all the generated instances, we could not find two candidate paths for each VNF option, hence the reduced size of Π w .

The final set of candidate paths calculated for each traffic demand is used to find the optimal solution with Gurobi solver. However, when running the penalized best-response algorithm, we consider only the k ′ -shortest paths from the previously obtained set. For instance, a traffic demand in the AARNET topology has on average 52 different paths that will be used to find the optimal solution. For k ′ = 10, only the 10 shortest paths will be considered in the penalized best-response algorithm.

Numerical Results

We report in Tables 4 and5 the performance of the penalized Best Response (BR) algorithm for different cost functions under distinct topologies. We set k ′ = 10, and we evaluate the algorithm performance by calculating the relative gap between the penalized BR solution and the optimal solution from Gurobi. Following the obtained results, the penalized BR algorithm provides very good quality solutions for the different cost functions with impressive computation time. Its average relative gap to the optimal solution is under 5%, and it is almost negligible for small topologies such as CESNET that has a result value of at most 0.07%. The average execution time of the algorithm depends on the size of the topology and on the cost function used. It is also directly proportional to the number of paths considered for each demand. Nevertheless except for the PWL cost function where the MILP solver is known to be efficient and quick at solving, the solver is significantly slower for the other more complex cost functions. In the case of the quadratic cost function, the maximum relative error of the BR algorithm over the 100 instances is 3.725%. However, the if we look at the computing times in Table 5, we observe that the BR algorithm is at most 16 times faster and at worst 6.3 times faster than the solver. It is interesting to note that for the M/M/1 cost function, the solver has reached its time limit that has been set to 5 minutes, while the BR algorithm takes less than 10 seconds to execute with a solution as close as 1.43% to the solver. The box plot in Figure 4 compares the distribution of the 100 instances generated for each scenario. The line that divides the box into two parts marks the median (Q2) of the data. The upper quartile (Q3) value is represented by the upper box line. 75% of the values fall below Q3. The lower quartile (Q1) value is the lower box line. 25% of the values fall below Q1. Hence, the values inside the box represent 50% of the data. The end of the lower and upper whiskers represents the minimum and maximum values, respectively. We first note that there are no outliers in the data samples which means that all the values of the data fall within the normal range. We observe that the values obtained with the quadratic cost function are more dispersed than the values obtained from the other cost functions.

In order to evaluate the choice of k ′ , we report in Table 6 the average relative gap and the computing times over different values of k ′ . We evaluate the results with the quadratic cost function (similar results were observed for the other cost functions). As expected, the optimality gap decreases when more routing paths are added to the set of candidate paths. However, the decrease of this gap is insignificant for k ′ > 10. In the case of the CESNET topology and since on average there are only 5 possible paths for a given demand (refer to Table 3), the gap optimality remains constant for k ′ ≥ 5. On the other hand, increasing the routing paths increases the running time of the BR algorithm. In this work, we found that the best trade-off between quality solution and running time is for k ′ = 10.

Conclusion

In this paper, we have considered the offline VNF placement and chaining problem, assuming that a non-linear cost function is associated to each network resources. It is also assumed that the volume of a traffic demand can change after each VNF. With respect to previous works, the main originality of the considered model is that the cost per unit capacity of a resource is not constant, but instead that it grows with its utilization rate, which is an essential feature to achieve a better load distribution. We have formulated the problem as a single-path routing problem in an extended network and adapted an existing game-theoretic algorithm to solve it. Our numerical results suggest that the algorithm provides near-optimal solutions in a substantially lower computing times than the solver, in particular for highly non-linear cost functions.

A Proof of Theorem 1

Before proving Theorem 1, we remark that if Algorithm 1 converges, it necessarily converges to a Nash equilibrium of the game. By definition, the game is in a Nash equilibrium if and only if no player has anything to gain by unilaterally deviating from its current strategy, that is, the current strategy of each player is its best response to the strategy of the others. Formally,

π (n) = π (n) w w∈W
is a Nash equilibrium of our game if and only if:

c w (π (n) w , π (n) -w ) ≤ c w (π, π (n)
-w ), ∀π ∈ Π w , ∀w ∈ W, which is precisely the condition under which Algorithm 1 stops (see line 7).

In order to prove property (d), we will need the following technical lemma which is an extension of Lemma 2 in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF].

-x d+1 = d k=0 x d-k N j=1 a j (x + b j ) k -b j (x -a j ) k , (23)
where

x = x + N j=1 (a j -b j ) Proof. We have xd+1 -x d+1 =   N j=1 a j - N j=1 b j   d k=0 xk x d-k , = d k=0 x d-k N j=1 a j xk - d k=0 x d-k N j=1 b j xk , = d k=0 x d-k N j=1 a j (x + b j -b j ) k - d k=0 x d-k N j=1 b j (x + a j -a j ) k , = d k=0 x d-k N j=1 a j k n=0 k n (x + b j ) n (-b j ) k-n - d k=0 x d-k N j=1 b j k n=0 k n (x -a j ) n (a j ) k-n , = d k=0 x d-k N j=1 a j (x + b j ) k - d k=0 x d-k N j=1 b j (x -a j ) k ,
where the last equality follows from a j b j = 0.

We will also need Lemma 2 below (which is basically a simple application of Hölder's inequality).

Proof. See the proof of Lemma 3 in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF].

We are now in position to prove Theorem 1.

F (π) = r∈R a r y r (π) d+1 . (30)
Given an instance of the problem, let π be the solution returned by Algorithm 1 and π * be any optimal solution. Consider a traffic demand w ∈ W and consider the solution obtained from π by deviating this demand from path π w to path π * w , with the paths of the other demands unchanged. As π is a Nash equilibrium, we have c w (π * w , π -w ) ≥ c w (π w , π -w ), and with (26) it follows that

F (π * w , π -w ) ≥ F (π). ( 31
)
Consider now the impact of this deviation on the traffic flowing on the edges of the expanded network. Let y w r,k (π w ) and y w r,k (π * w ) be the traffic sent by demand w on the edge (r, k) of the expanded network before and after the deviation, respectively. We have

Figure 2 :

 2 Figure 2: Expanded network for the example of Figure 1. The artificial edges between consecutive layers are shown with dotted lines. Two feasible paths for the demand from node s 1 to node t 1 are shown, with blue-coloured edges for the first one and with red-coloured edges for the other one (common edges are purple-coloured).

  r∈R u̸ =w y u r (πu) ϕr y -w r + y w r (π) -ϕr y -w r = r∈π y -w r ϕr y -w r + y w r (π) -ϕr y -w r (8)

Theorem 1 .

 1 Algorithm 1 has the following properties: (a) It converges in a finite number of steps. (b) If it ever reaches a global optimum of problem (OPT), it returns this optimal solution.

Figure 3 :

 3 Figure 3: NSFNET topology which has been expanded with function nodes v 1 , v 2 and v 3 . The nodes s 1 , s 2 and s 3 are the source nodes and the nodes t 1 , t 2 and t 3 are the destination nodes of the traffic demands.

Lemma 1 .

 1 Let d ∈ N and x ≥ 0. For any non-negative numbers a 1 , a 2 , . . . , a N and b 1 , b 2 , . . . , b N such that a i b i = 0 for all i = 1, 2, . . . , N , it holds that xd+1

Lemma 2 .

 2 For non-negative numbers a e , x e z e , and non-negative integer m, e a e (x e + z e ) m+1 ≤

  , = δ r,k πw b w k,r + (δ r,k π * w -δ r,k πw )b w k,r , = y w r,k (π w ) + (δ r,k π * w -δ r,k πw )b w k,r , π * w , π -w ) = y w r (π * w ) + y -w r , = y w r (π w ) + y -w r, k) ∈ π * w \ π w , 0 otherwise, and α - r,k (w) = 1 if (r, k) ∈ π w \ π * w , 0 otherwise, The constant α + r,k(w) is 1 if the deviation increases the traffic of demand w on the edge (r, k) of the expanded network, and similarly, α - r,k (w) is 1 if the deviation decreases the traffic of demand w on the edge (r, k) of the expanded network. Note that α+ r,k (w) α - r,k (w) = 0. As δ r,k π * w -δ r,k πw = α + r,k (w) -α - r,k(w), we can rewrite (35) as follows r∈R a r   y r (π) + K k=0 (α + r,k (w) -α - r,k (w))b w k,r d+1 -y r (π) d+1   ≥ 0. (36) Applying Lemma 23 to (36) with a k = α + r,k (w)b w k,r and b k = α - r,k (w)b w k,r , we obtain r,k (w) -α - r,k (w))b w k,r , (w)b w k,r (y r (π) + y r (π * )) r,k (w) -α - r,k (w))b w k,r = w∈W K k=0 (δ r,k π * w -δ r,k πw )b w k,r = y r (π * ) -y r (π).

Table 1 :

 1 Notations.

	G = (V, E)	Physical network, where V is the set of nodes and E is the set of links.
	D	Set of function nodes.
	R	Set of network resources (transmission links and function nodes).
	F (resp. F w ) Set of K network functions (resp. service chain for demand w).
	D(f )	Set of function nodes where VNF f can be executed.
	W	Set of traffic demands.
	s w (resp. t w ) Source (resp. destination) node of demand w.
	K w	Number of VNFs in the service chain F w .
	f w k b w k,r	k th VNF in the service chain F w .

Table 2 :

 2 Topologies : number of nodes and links.

	Topology	# Nodes # Links
	AARNET	19	46
	ARPANET	25	56
	NSFNET	13	30
	IBM	18	48
	CESNET	10	18
	s2	s1	t1
				f2, f3
	v2		t2	v3
	f1, f3
			s3
	t3		v1 f1, f2

Table 3 :

 3 Number of candidate paths Π w .

	Topology Max Min Average
	AARNET	128	32	52
	ARPANET	128	32	47
	NSFNET	128	32	48
	IBM	128	32	51
	CESNET	8	4	5

Table 4 :

 4 Relative gap (% ± std) to the optimal solution with k ′ = 10. Figure 4: Box plot of the relative gap to the optimal solution.

	Topology	PWL		Quadratic	M/M/1
	AARNET	0.71 ±0.65	1.81 ±0.85	1.43 ±1.76
	ARPANET	1.15 ±0.5		3.725 ±1.55	1.5 ±0.79
	NSFNET	0.42 ±0.49	1.67 ±1.12	0.31 ±0.32
	IBM		1.7 ±0.8		2.85 ±0.98	0.53 ±0.71
	CESNET	0.051 ±0.003 0.047 ±0.008 0.070 ±0.056
		8 10					PWL function Quadratic function M/M/1
	Relative gap (%)	4 6
		2
		0	AARNET	ARPANET	NSFNET	IBM	CESNET

Table 5 :

 5 Average computing times (seconds±std) for the 100 problem instances with k ′ = 10.

	Function Costs	AARNET BR Gurobi	ARPANET BR Gurobi	NSFNET BR Gurobi	BR	IBM Gurobi	CESNET BR Gurobi
	PWL	7.1±0	1.2	8.09±0	1.72	4.93±0	0.72	6.83±0	0.9	1.96±0	0.07
	Quadratic	7.89±0 126±0 8.89±0 124±0	6.3±0	40±0	8.27±0 109±0 2.33±0	0.3
	M/M/1	8.17±0	300	9.84±0	300	7.52±0	300	9.17		300	2.79±0	300

Table 6 :

 6 Average relative gap and computing times for the 100 problem instances with the quadratic cost function.

	Topology	gap	k ′ = 5 runtime	k ′ = 10 gap runtime	k ′ = 20 gap runtime	k ′ = 32 gap runtime
	AARNET	5.2%	3.8 sec	1.8%	7.9 sec	1.4%	17 sec	1.4% 26.7 sec
	ARPANET 5.5%	4 sec	3.7%	8.9 sec	3%	19 sec	3%	30 sec
	NSFNET	1.8%	2.9 sec	1.7%	6.3 sec	1.6%	12 sec	1.6%	18 sec
	IBM	5.88% 3.5 sec 2.85% 8.2 sec 2.77% 12.8 sec 2.75%	30 sec
	CESNET	0.05% 1.9 sec 0.05% 2.3 sec 0.05% 2.4 sec 0.05% 2.4 sec

Proof of Theorem 1. We first prove property (a). We observe that the cost function c w (π w , π -w ) of player w can be rewritten as follows c w (π, π -w ) = f w (π, π -w ) + p w (π, π -w ),

where the penultimate equality is obtained by noting that y r (π, π -w ) = y -w r for all resources r ̸ ∈ π. Introducing F (π -w ) = r∈R y -w r ϕ r (y -w r ), which represents the network cost that would be obtained if demand w was removed from the network, (25) yields

which means that the network cost function F () is a potential function of the game, that is, any best-response move of a player w decreases the total network cost. As Algorithm 1 computes the best response of all players at each iteration and as the network cost F () can take only a finite number of values, this implies the convergence of the algorithm in a finite number of steps. We now prove property (b). Assume that Algorithm 1 has reached an optimal solution π * of problem (OPT). The strategy of each player w is therefore π * w . For any player w and any path π ∈ Π w , it follows from (26) that

which implies that π * is Nash equilibrium of our game. Hence, Algorithm 1 returns this solution.

Consider now property (c). Note that it directly follows from property (d). Nevertheless, we provide a direct proof below. Assume that ϕ r (y) = a r for all y ≥ 0 and all resources r ∈ R. The network cost F (π) can then be written as follows

from which it follows that it is optimal to route each demand w on a path

Note from (8) that p w (π, π -w ) = 0 in this case, from which it follows that

Therefore, the best-response strategy of a player is an optimal path, independently of the strategies of the other players. This implies that Algorithm 1 returns a globally optimal solution for linear resource costs.

Finally, let us prove property (d). We focus on the case of monomial functions ϕ r (y) = a r y d , as the extension to polynomials is straightforward (see the proof of Theorem 3 in [START_REF] Brun | A penalized best-response algorithm for non-linear singlepath routing problems[END_REF]). In that case, we have Hence, summing inequalities (38) over all w ∈ W, we get r∈R a r y r (π * )

which gives

from which we obtain

which concludes the proof.