What frequency bandwidth to run cellular network in a given country? -a downlink dimensioning problem Bartłomiej Błaszczyszyn † and Mohamed K. Karray * Abstract-We propose an analytic approach to the frequency bandwidth dimensioning problem, faced by cellular network operators who deploy/upgrade their networks in various geographical regions (countries) with an inhomogeneous urbanization. We present a model allowing one to capture fundamental relations between users' quality of service parameters (mean downlink throughput), traffic demand, the density of base station deployment, and the available frequency bandwidth. These relations depend on the applied cellular technology (3G or 4G impacting user peak bit-rate) and on the path-loss characteristics observed in different (urban, sub-urban and rural) areas. We observe that if the distance between base stations is kept inversely proportional to the distance coefficient of the path-loss function, then the performance of the typical cells of these different areas is similar when serving the same (per-cell) traffic demand. In this case, the frequency bandwidth dimensioning problem can be solved uniformly across the country applying the mean cell approach proposed in [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF]. We validate our approach by comparing the analytical results to measurements in operational networks in various geographical zones of different countries.

Index Terms-cellular network; dimensioning; bandwidth; spatial heterogeneity, QoS-homogeneous networks

I. INTRODUCTION

A systematic increase of the traffic in wireless cellular networks leads to a potential degradation of the users' quality of service (QoS). In order to prevent such degradation, network operators (besides upgrading the technology) keep adding new base stations and/or supplementary frequency bandwidth. Planning of this network dimensioning process requires the knowledge of the relations between the QoS parameters, traffic demand, density of base station deployment, and the operated frequency bandwidth. Such relations are usually developed for homogeneous network models representing, separately, a typical city or a rural area. To the best of our knowledge, a global approach to some larger geographical zone, as e.g. some given country, typically with an inhomogeneous urbanization, has not yet been proposed. This is the main goal of this paper.

Our work is based on [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF], where a typical-cell approach was proposed to the aforementioned study of the QoS, with the information-theoretic characterization of the peak-bit rate, the queueing-theoretic evaluation of the cell performance, and a stochastic-geometric approach to irregular but homogeneous network deployment. In this paper we extend this approach to an inhomogeneous network deployment, i.e., a network which is supposed to cover urban, sub-urban and rural areas, for which path-loss characteristics and the network deployment densities are different.

A key observation that we make in this regard is that if the distance between neighbouring base stations is kept inversely proportional (equivalently, the square root of the network density is proportional) to the distance coefficient of the path-loss function, then the typical cells of these different areas exhibit the same "response" to the traffic demand. Specifically, they exhibit the same relation between the mean user throughput and the per-cell traffic demand. Thus, one can say that this inhomogeneous network offers to its users homogeneous "QoS response" across all different areas. More precisely: cells in rural areas will be larger than these in urban areas, and the traffic demand per surface can be different, with urban areas generally serving more traffic. However, in the case when the mean traffic demand per cell is equal, these areas will "offer" the same user's throughput.

As a consequence, the bandwidth dimensioning problem can be addressed and solved for one particular, say urban, area and applied to the whole network, appropriately adapting the density of base station in sub-urban and rural areas to the specific value of the distance coefficient in the path-loss function.

We validate our approach by comparing the obtained results to measurements in 3G and 4G operational networks in different countries in Europe, Africa and Asia.

A. Related works

The evaluation of the performance of cellular networks is a complex problem, but crucial for network operators. It motivates a lot of engineering and research studies. The complexity of this problem made many actors develop complex and time consuming simulation tools such as those developed by the industrial contributors to 3GPP (3rd Generation Partnership Project) [START_REF] Gpp | TR 36.814-V900 Further advancements for E-UTRA -Physical Layer Aspects[END_REF]. There are many other simulation tools such as TelematicsLab LTE-Sim [START_REF] Piro | Simulating LTE cellular systems: an open source framework[END_REF], University of Vien LTE simulator [START_REF] Mehlführer | The Vienna LTE simulators -Enabling reproducibility in wireless communications research[END_REF][START_REF] Simko | Optimal pilot symbol power allocation under time-variant channels[END_REF] and LENA tool [6,[START_REF] Baldo | A new model for the simulation of the LTE-EPC data plane[END_REF] of CTTC, which are not necessarily compliant with 3GPP.

A possible analytical approach to this problem relies on the information theoretic characterization of the individual link performance; cf e.g. [START_REF] Goldsmith | Variable-rate variable-power MQAM for fading channels[END_REF][START_REF] Mogensen | LTE Capacity Compared to the Shannon Bound[END_REF], in conjunction with a queueing theoretic modeling and analysis of the user traffic; cf. e.g. [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF][START_REF] Bonald | Wireless downlink data channels: user performance and cell dimensioning[END_REF][START_REF] Hegde | Capacity of multiservice WCDMA Networks with variable GoS[END_REF][START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF][START_REF] Rong | Performance evaluation of cellular networks offering TV services[END_REF][START_REF] Karray | A queueing theoretic approach to the dimensioning of wireless cellular networks serving variable bit-rate calls[END_REF]. Recently, [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF][START_REF] Błaszczyszyn | Performance laws of large heterogeneous cellular networks[END_REF] proposes an approach combining queueing theory and stochastic geometry allowing to study the dependence of the mean user throughput on the traffic demand in a typical city. This analytical approach is validated with field measurements at the scale of a city. At this scale, not only the network average metrics but also their spatial distributions in the considered, homogeneous, zones are well predicted in [START_REF] Jovanovic | QoS and network performance estimation in heterogeneous cellular networks validated by real-field measurements[END_REF]. This prior work does not consider inhomogeneous networks, observed at the scale of a larger geographical zone or a whole country. In particular, the simulation tools are too time consuming to be applicable at this scale.

B. Paper organization

We describe our general cellular network model, comprising the geometry of base stations, propagation loss model, traffic dynamics, and service policy in Section II. In Section III we propose a model permitting to get the relation of users' QoS to the key network parameters at a country scale. This approach is validated by comparing, in Section IV, the obtained results to measurements in various operational networks in different countries. We show also how to solve the bandwidth dimensioning problem.

II. HOMOGENEOUS MODEL DESCRIPTION

In this section we briefly recall the model from [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF].

A. Network geometry and propagation

The network is composed of base stations (BS) whose locations are modelled by a stationary point process Φ = {X n } n∈Z on R 2 with intensity λ > 0. We assume that Φ is stationary and ergodic. Each base station X n emits a power denoted by P n . We assume that {P n } n∈Z are i.i.d. marks of Φ.

The propagation loss comprises a deterministic effect depending on the distance between the transmitter and receiver called path-loss, and a random effect called shadowing. The path-loss is modeled by a function

l (x) = (K |x|) β , x ∈ R 2 , (1)
where K > 0 and β > 2 are two parameters called distance coefficient and exponent respectively. The shadowing between BS X n and all the locations y ∈ R 2 is modelled by a stochastic process S n (y -X n ). We assume that the processes {S n (•)} n∈Z are i.i.d. marks of Φ.

The inverse of the power received at location y from BS X n is denoted by

L Xn (y) = l (y -X n ) P n S n (y -X n ) , y ∈ R 2 , n ∈ Z
and will be called, with a slight abuse of terminology, propagation loss between y and X n . In order to simplify the notation, we shall omit the index n of BS X n . Each BS X ∈ Φ serves the locations where the received power is the strongest among all the BS; that is

V (X) = y ∈ R 2 : L X (y) ≤ L Y (y) for all Y ∈ Φ (2)
called cell of X.

The signal-to-interference-and-noise (SINR) power ratio in the downlink for a user located at y ∈ V (X) equals

SINR (y, Φ) = 1/L X (y) N + Y ∈Φ\{X} ϕ Y /L Y (y) , ( 3
)
where N is the noise power and each ϕ Y ∈ [0, 1] is some interference factor accounting for the activity of BS Y in a way that will be made specific in Section II-C. We assume that {ϕ Y } Y ∈Φ are additional (not necessarily independent) marks of the point process Φ.

A single user served by BS X and located at y ∈ V (X) gets a bit-rate R (SINR (y, Φ)), called peak bit-rate, which is some function of the SINR given by [START_REF] Piro | Simulating LTE cellular systems: an open source framework[END_REF]. Particular form of this function depends on the actual technology used to support the wireless link.

B. Traffic dynamics

We consider variable bit-rate (VBR) traffic; i.e., users arrive to the network and require to transmit some volume of data at a bit-rate decided by the network. Each user arrives at a location uniformly distributed and requires to transmit a random volume of data of mean 1/µ. The duration between the arrivals of two successive users in each zone of surface S is an exponential random variable of parameter γ × S (i.e. there are γ arrivals per surface unit). The arrival locations, inter-arrival durations as well as the data volumes are assumed independent of each other. We assume that the users don't move during their calls. The traffic demand per surface unit is then equal to ρ = γ/µ, which may be expressed in bit/s/km 2 .

The traffic demand to a given cell equals

ρ (X) = ρ |V (X)| , X ∈ Φ, ( 4
)
where |A| is the surface of A.

We shall assume that each user in a cell gets an equal portion of time for his service. Thus when there are k users in a cell, each one gets a bit-rate equal to his peak bit-rate divided by k. More explicitly, if a base station located at X serves k users located at y 1 , y 2 , . . . , y k ∈ V (X) then the bit-rate of the user located at y j equals 1 k R (SINR (y j , Φ)), j ∈ {1, 2, . . . , k}.

C. Cell performance metrics

We consider now the stationary state of the network in the long run of the call arrivals and departures. Using queuing theory tools, it is proven in [15, Proposition 1] that:

• Each base station X ∈ Φ can serve the traffic demand within its cell if this latter doesn't exceed some critical value which is the harmonic mean of the peak bit-rate over the cell; that is

ρ c (X) := |V (X)| V (X) 1/R (SINR (y, Φ)) dy -1 , ( 5
) • The mean user throughput in cell V (X) equals r (X) = max(ρ c (X) -ρ (X) , 0). (6)
• The mean number of users in cell V (X) equals

N (X) = ρ (X) r (X) . (7)
• Moreover, we define the cell load as

θ (X) = ρ (X) ρ c (X) = ρ V (X) 1/R (SINR (y, Φ)) dy. ( 8
)
• The probability that the base station has at least one user to serve (at a given time) equals

p (X) = min (θ (X) , 1) . (9)
We assume that a BS transmits only when it serves at least one user. Then, as proposed in [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF], we take ϕ Y = p(Y ) in the SINR expression [START_REF] Piro | Simulating LTE cellular systems: an open source framework[END_REF]. Thus (8) becomes

θ (X) = ρ V (X) 1/R   1/L X (y) N + Y ∈Φ\{X} min(θ(Y ),1) LY (y)   dy, (10)
which is a system of equations with unknown cell loads {θ (X)} X∈Φ .

D. Network performance metrics

The network performance metrics are defined by averaging spatially over all the cells in the network in an appropriate way [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF]. In particular, it follows from the ergodic theorem for point processes [START_REF] Daley | An introduction to the theory of point processes[END_REF]Theorem 13.4.III] that the average cell load equals

lim |A|→∞ X∈Φ∩A θ(X) Φ(A) = E 0 [θ (0)] = ρ λ E 1 R (SINR (0, Φ))
,

where A is a ball centered at the origin having radius increasing to infinity and E 0 is the expectation with respect to the Palm probability associated to Φ. The second equality in the above equation is proved in [1, Proposition 3]. Moreover, it is observed in [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF] that the network performance metrics are well approximated by the mean cell model. Specifically, we define a virtual cell having traffic demand ρ := E 0 [ρ (0)] and load θ := E 0 [θ (0)]. For other performance metrics, the mean cell mimics the behavior of the true cells given in Section II-C; that is it has critical traffic demand deduced from [START_REF] Goldsmith | Variable-rate variable-power MQAM for fading channels[END_REF] ρc := ρ θ mean user's throughput deduced from ( 6)

r := max (ρ c -ρ, 0) ,
and the mean number of users deduced from ( 7)

N := ρ r .
The load equations [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF] become for the mean cell

θ = ρ λ BS E 1/R 1/L X * (0) N + θ Y ∈Φ\{X} 1/L Y (0)
, where X * is the location of the BS whose cell covers the origin.

III. FROM SCALING EQUATIONS TO INHOMOGENEOUS

NETWORKS

We begin by studying some scaling laws observed in homogeneous networks. Then we will introduce some inhomogeneous networks which are able to offer homogeneous QoS response to the traffic demand.

A. Scaling laws for homogeneous networks

Consider a homogeneous network model descried in Section II. For α > 0 consider a network obtained from this original one by scaling the base station locations Φ ′ = {X ′ = αX} X∈Φ , the traffic demand intensity ρ ′ = ρ/α 2 , distance coefficient K ′ = K/α and shadowing processes S ′ n (y) = S n y α , while preserving the original (arbitrary) marks (interference factors) ϕ ′ X ′ = ϕ X , X ∈ Φ and powers P ′ n = P n . For the rescaled network consider the cells V ′ (X ′ ) given by ( 2) and their characteristics ρ 4), ( 5), ( 6), ( 7), ( 8) and ( 9), respectively.

′ (X ′ ), ρ ′ c (X ′ ), r ′ (X ′ ), N ′ (X ′ ), θ ′ (X ′ ), p ′ (X ′ ) calculated as in (
Proposition 1: For any

X ′ ∈ Φ ′ , we have V ′ (αX n ) = αV (X n ) while ρ ′ (X ′ ) = ρ(X), ρ ′ c (X ′ ) = ρ c (X), r ′ (X ′ ) = r(X), N ′ (X ′ ) = N (X), θ ′ (X ′ ) = θ(X).
Proof. Observe that

L ′ αXn (y) = P n S ′ n (y -αX n ) (K ′ |y -αX n |) β = P n S n y α -X n K y α -X n β = L Xn y α
and by (2) one gets the required dilation relation for the individual cells. Moreover, using (3) one gets by simple algebra SINR (y, Φ ′ ) = SINR (y/α, Φ). Starting from ( 5) and making the change of variable z = y/α, it follows that

ρ ′ c (X ′ ) = |V ′ (X ′ )| V ′ (X ′ ) 1/R (SINR (y, Φ ′ )) dy -1 = α 2 |V (X)| αV (X)
1/R (SINR (y/α, Φ)) dy

-1 = ρ c (X) .
The remaining desired equalities then follow from the fact that

ρ ′ (X ′ ) = ρ α 2 |αV (X)| = ρ |V (X)| = ρ (X) .

Corollary 1:

Consider the rescaled network as in Proposition 1 with the interference factors taken ϕ X = min (θ (X) , 1). Then the load equations [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF] are the same for the two networks Φ and Φ ′ . Therefore, the load factors solving these equation are the same θ ′ (X ′ ) = θ (X) and consequently,

ρ ′ c (X ′ ) = ρ c (X), r ′ (X ′ ) = r(X), N ′ (X ′ ) = N (X), θ ′ (X ′ ) = θ(X).
From the above observations we can deduce now that the considered scaling of the network parameters preserves the mean cell characteristics defined in Section II-D. Denote by E ′0 the expectation with respect to the Palm probability associated to Φ ′ .

Corollary 2: Consider the rescaled network as in Proposition 1 with arbitrary interference factors, possibly satisfying the load equations [START_REF] Borst | User-level performance of channel-aware scheduling algorithms in wireless data networks[END_REF].

Then E ′0 [ρ ′ (0)] = E 0 [ρ (0)] and E ′0 [θ ′ (0)] = E 0 [θ (0)].
Consequently, the mean cells characteristics associated to Φ and Φ ′ as described in Section II-D are identical.

Proof. It follows from the inverse formula of Palm calculus [START_REF] Baccelli | Stochastic Geometry and Wireless Networks, Volume I -Theory, ser. Foundations and Trends in Networking[END_REF]Theorem 4

.2.1] that E 0 [ρ (0)] = ρ λ . Similarly, E ′0 [ρ ′ (0)] = ρ ′ λ ′ = ρ
λ which proves the first desired equality.

The second equality follows e.g. by the ergodic argument

E ′0 [θ ′ (0)] = lim |A|→∞ 1 Φ ′ (A) X ′ ∈Φ ′ ∩A θ ′ (X ′ ) = lim |A|→∞ 1 Φ(A/α) X∈Φ∩A/α θ(X) = lim |A|→∞ 1 Φ(A) X∈Φ∩A θ(X) = E 0 [θ (0)] .
where the second equality is due to Proposition 1.

B. Inhomogeneous networks with homogeneous QoS response to the traffic demand

Consider a geographic region (say a country), which is composed of urban, suburban and rural areas. The parameters K and β of the distance loss model (1) depend on the type of the zone. For example the COST-Hata model [START_REF]Evolution of land mobile radio (including personal) communications[END_REF] gives the distance loss function in the form 10 log 10 (l (x)) = A + B log 10 (|x|) where A and B are given in Table I and the distance |x| is in km. The corresponding value of the parameter K = 10 A/B as well as the ratio K urban /K are also given in this table.

The density of base stations also depends on the type of the zone: it is usually much higher in urban than in rural areas. Indeed, the distance D between two neighboring base stations in rural areas may be up to 10 times larger than in dense urbane zones (where a typical value of D is about 1km).

Consider now a situation where the average distance between neighbouring base stations is kept inversely proportional to the distance coefficient of the path-loss function: D × K = const, or, in other words,

K i / √ λ i = const, (11)
where λ is the network density in the given zone (assumed homogeneous) and K i is the path-loss distance factor observed in this zone. In our example shown in Table I we should thus have D = 1, 5 and 8 km respectively for urban, suburban and rural zones. 1 Then the scaling laws proved in Section III-A say that locally, for each homogeneous area of this inhomogeneous network, one will observe the same relation between the mean user throughput and the (per-cell) traffic demand. Cells in rural areas will be larger than these in urban areas, however in the case when they are charged in the same way by the traffic demand (which, per unit surface, needs to be correspondingly larger in urban areas) then they offer the same QoS. In other words, one relation is enough to capture the key dependence between the QoS, network density, bandwidth and the traffic demand for different areas of this network. It can be used for the network dimensioning. 2 1 Indeed, when operators deploys networks, they firstly aim to assure some coverage condition which has the form (D × K) β = const. Moreover, in urban zones, networks have to be densified not only for capacity constraints, but also to assure coverage for indoor users. 2 The approximation of an inhomogeneous network by a piecewise homogeneous network can be made more precise e.g. following the ideas presented in [START_REF] Błaszczyszyn | Approximate decomposition of some modulated-Poisson Voronoi tessellations[END_REF][START_REF] Błaszczyszyn | Approximations of functionals of some modulated-Poisson Voronoi tessellations with applications to modeling of communication networks[END_REF]. Note also that in our approach we ignored the fact that the pathloss exponent β is not the same for urban and the two other zones considered in Table I. This is clearly an approximation. Its quality can be studied using a network equivalence approach [START_REF] Błaszczyszyn | Equivalence and comparison of heterogeneous cellular networks[END_REF] allowing for different path-loss exponents.

Environment

A B K = 10

IV. NUMERICAL RESULTS

The relations of the network performance metrics to the traffic demand described in Section II-D were already validated in [START_REF] Błaszczyszyn | How user throughput depends on the traffic demand in large cellular networks[END_REF] comparing them to operational network measurements in some typical cities in Europe. Thus we shall focus on the validation of its extension to more large areas, and ultimately to a whole country as proposed in Section III. Once this validation carried, we shall solve numerically the bandwidth dimensioning problem.

A. Real-field measurements

We describe now the real-field measurements. The raw data are collected using a specialized tool which is used by operational engineers for network maintenance. This tool measures several parameters for every base station 24 hours a day. In particular, one can get the traffic demand, cell load and number of users for each cell in each hour. Then we estimate the network performance metrics for each hour averaging over all considered cells. The mean user throughput is calculated as the ratio of the mean traffic demand to the mean number of users. The mean traffic demand is used as the input of our analytical model.

B. A reference country

We consider a reference country, called in what follows Country 1, divided into 5 regions. Each of these regions comprises a large mix of urban, suburban and rural zones. In what follows we consider 3G and 4G networks deployed in Country 1.

1) 3G network: a) Numerical setting: Country 1 is covered by a 3G network at carrier frequency f 0 = 2.1GHz with frequency bandwidth W = 5MHz, base station power P = 60dBm and noise power N = -96dBm. The peak bit-rate equals to 30% of the ergodic capacity of the flat fading channel; i.e., The shadowing random variable is lognormal with unit mean and logarithmic-standard deviation σ = 9.6dB and mean spatial correlation distance 50m.

R (SINR) = 0.3 × W E log 2 1 + |H| 2 SINR
Each BS comprises three antennas having each a threedimensional radiation pattern specified in [START_REF] Gpp | TR 36.814-V900 Further advancements for E-UTRA -Physical Layer Aspects[END_REF]. The BS and mobile antenna heights equal 30m and 1.5m respectively. The pilot channel power is taken equal to 10% of the total base station power.

The theoretical relations of the network performance metrics to the traffic demand are calculated for such typical urban zone and compared, to measurements in large areas, and ultimately the whole country, applying the approach developed Section III. c) Network performance: Figure 1(a) shows the mean cell load θ versus mean traffic demand per cell ρ = ρ/λ calculated analytically and obtained from measurements for the 5 considered regions. The mean number of users N and user's throughput r versus mean traffic demand per cell ρ are plotted in Figures 1(b) and 1(c) respectively. Observe that the analytical curves fits well to the measurements.

In order to analyze the effect of the size of the regions over which we average the cell performance metrics, we decompose Country 1 into a regular grid composed of squares (meshes) of equal side. Figure 2(a) shows the mean cell load versus traffic demand for mesh sizes 3, 10, 30 and 100km respectively. Observe that the larger is the mesh size, the closer are measurements to the analytical (mean) model.

We aim now to check the hypothesis made in Section III-B, namely that zones with different density of base stations observe the same dependence of the mean cell characteristics on the mean per-cell traffic demand. Figure 2(b) shows the mean cell load versus traffic demand for mesh size 30km where we distinguish the meshes according to the distance between neighboring BS; distances within [0, 2] km (resp. [6,[START_REF] Goldsmith | Variable-rate variable-power MQAM for fading channels[END_REF] km) corresonding to urban (resp. rural) zones. Observe that there is no apparent difference between urban, sub-urban and rural meshes.

2) 4G network: The peak bit-rate is calculated by R (SINR) = b×W log 2 (1 + SINR/a) where a = 3, b = 1.12 are estimated from fitting to the results of a link simulation tool.

a) Numerical setting for carrier frequency 2.6GHz: We consider a 4G network at carrier frequency f = 2.6GHz with frequency bandwidth W = 20MHz, base station power P = 63dBm and noise power N = -90dBm.

The typical urban zone is the same as described in Section IV-B1b except for the distance loss parameter K. Following the COST-Hata model [START_REF]Evolution of land mobile radio (including personal) communications[END_REF] this latter is calculated by assuming that it is proportional to f 2/β ; which gives

K = K 0 × (f /f 0 ) 2/β = 7964km -1 .
b) Numerical setting for carrier frequency 800MHz: We consider also a 4G network at carrier frequency f = 800MHz with frequency bandwidth W = 10MHz, base station power P = 60dBm and noise power N = -93dBm.

The typical urban zone is the same as described in Section IV-B1b except for the average distance between two neighboring BS which is now D u ≃ 1.5km and the distance loss parameter K = K 0 × (f /f 0 )

2/β = 4283km -1 . Figures 4(a) and 4(b) show the mean user's throughput versus mean traffic demand per cell for the 4G network at carrier frequencies 2.6GHz and 800MHz respectively. Even if currently we do not observe traffic demand larger than 1000kbit/s per cell, we show the predicted curves in their whole domains as they will serve for the network dimensioning, cf. Section IV-D.

C. Two other countries

We consider now two 3G networks in two other countries, called Country 2 and Country 3. The numerical setting is the same as in Section IV-B1a. Since we don't have detailed information on the propagation characteristics in these countries, we make the hypothesis that the typical urban zone has the same characteristics as for Country 1; i.e. those described in Section IV-B1b.

Figures 5(a), 5(b) and 5(c) show the mean cell load, users number and user's throughput versus mean traffic demand per cell for the two considered countries. Note that for Country 3 the user's throughput stagnates for traffic demand exceeding 800kbit/s/cell. This is due to congestion; some users being blocked at access as may be seen in Figure 5(b).

D. Bandwidth dimensioning

We come back now to the question stated in the title of this article: What frequency bandwidth is required to run cellular network guaranteeing sufficient QoS for its users? This is the frequency dimension problem often faced by the network operators.

Here we consider the mean user throughput in the downlink as the QoS metric. We have remarked at the end of Section III-B that (conjecturing) the relation [START_REF] Bonald | Wireless downlink data channels: user performance and cell dimensioning[END_REF] allows one to capture the key dependence between the mean user throughput and the per cell traffic demand (given the network density and the frequency bandwidth) for an inhomogeneous network focusing only on one homogeneous type of network area. We use for this purpose the urban area. Examples of the corresponding user-throughput versus traffic demand curves have been presented on Figures 1(c), 4(a) and 4(b). Note that these curves depend on (increase with) the frequency bandwidth. The solution of the frequency dimensioning problem (for a given technology and a reference, urban, inter-BS distance) consists in finding, for each value of the traffic demand per cell, the minimal frequency bandwidth such that the predicted mean user throughput reaches a given target value.

Figure 6 plots the solutions for this problem in several considered cases. More precisely, it gives the predicted frequency bandwidth allowing one to assure the target r = 5 Mbit/s mean user's throughput in the network, for a given (per cell) traffic demand. We show solutions for 3G and 4G networks with As a possible application of the above dimensioning strategy let us mention the following actual problem. In Country 1, the configurations currently deployed for 4G are f = 2.6GHz with D u = 1km and f = 800MHz with D u = 1.5km. The operator is wondering whether to densify the network for f = 800MHz decreasing D u to 1km. Figure 6 shows how much bandwidth can be economized in this case.

V. CONCLUSION

We propose a model permitting to calculate the performance of 3G and 4G wireless cellular networks at the scale of a whole country comprising urban, suburban and rural zones. This model relies on an observation that the distance coefficient in the propagation loss function depends on the type of zone in such a way that its product to the distance between neighboring base stations remains approximately constant. It is then shown that the network performance, in terms of the relations of the mean cell load, number of users and user's throughput to mean traffic demand, are the same for the different types of zones.

This theoretical result is validated by field measurements in 3G and 4G cellular networks in various countries. Then we solve the bandwidth dimensioning problem for 3G and 4G networks at different carrier frequencies; i.e. we plot the frequency bandwidth as function of the traffic demand per cell to assure a user's throughput of 5Mbit/s. This curve is crucial for operators to predict the frequency bandwidth required to serve the continuing increase of traffic in the next decades.

We shall attempt in future work to extend the present approach to the uplink; in particular to account for power control effect and the specificity of the locations of the interferers in this case.

where |H| 2 1 .

 21 is a unit mean exponential random variable representing the fading and E [•] is the expectation with respect to H. b) Typical urban zone: We consider a typical urban zone (a city) in the considered Country 1. Knowing all BS coordinates and the surface of the deployment zone we deduce the BS density λ = 1.15 stations per km 2 which corresponds to an average distance between two neighboring BS of D u ≃ 1/ √ λ ≃ 1km. For the analytical model, the locations of BS is modelled by a homogeneous Poisson point process with intensity λ. The parameters of the distance loss model (1) are estimated from the COST Walfisch-Ikegami model [20] which gives β = 3.8, K 0 = 7117km -1 . Mean cell characteristics versus traffic demand calculated analytically and estimated from measurements for a 3G network in Country 1. cell [kbit/s] Network 3G, Carrier frequency 2.1GHz Analytical model Inter-BS distance in [6,8[ km [4,6[ km [2,4[ km [0,2[ km (b) Decomposition into different density zones Fig. 2. Mean cell load versus traffic demand in a 3G network at 2.1GHz in Country 1, with two different decompositions.

Fig. 3 .Fig. 4 .

 34 Mean cell load versus traffic demand for 4G in Country 1 at two different frequency bandwidths. Mean user throughput versus traffic demand for 4G in Country 1 at two different frequency bandwidths. c) Network performance: Figures 3(a) and 3(b) show the mean cell load versus mean traffic demand per cell calculated analytically and obtained from measurements for the 4G network at carrier frequencies 2.6GHz and 800MHz respectively. Observe again that the analytical curve fits well with the measurements.

5 .

 5 Mean cell characteristics versus traffic demand calculated analytically and estimated from measurements for a 3G network in Country 2 and 3.different frequency carriers f and different inter-BS distances D u relative to the urban area.

We thank our colleagues from Orange for fruitful exchanges related to the present work.


MATHEMATICS, UNIV. MONTPELLIER B Mohammadi

J Tuomela

UNCERTAINTY QUANTIFICATION IN THE NUMERICAL SOLUTION OF COUPLED SYSTEMS BY INVOLUTIVE COMPLETION

Keywords: Mathematics Subject Classification. Primary 58J05; Secondary 35J40, 35S15 overdetermined PDE, Uncertainty Quantification, sensitivity analysis, chromatographic adsorption, constrained coupled systems . 1

We address the issue of epistemic uncertainty quantification in the context of constrained differential systems. To illustrate our approach we have chosen a certain chromatographic adsorption model which is a coupled system of partial differential, ordinary differential and algebraic equations. The difficulty in solving this type of a system is that typically certain unknowns lack a natural time evolution equation. The standard approach in such cases is to devise specific numerical schemes which somehow try to take into account the implicit structure of the system. In our approach, we complete the system by finding the appropriate missing evolution equations. This makes the system overdetermined and more complicated in some way but, on the other hand, the completed system provides extra information useful for error estimation and uncertainty quantification. We will also show that reducing the epistemic uncertainties also leads to better estimations of aleatory uncertainties.

Introduction

In the field of Uncertainty Quantification (UQ) the uncertanties are classified in epistemic (reducible) and aleatory (irreducible) categories [START_REF] Aiaa | Guide for the verification and validation of computational fluid dynamics simulations[END_REF][START_REF] Ghanem | On the construction and analysis of stochastic models: characterization and propagation of the errors associated with limited data[END_REF][START_REF] Iaccarino | Quantification of uncertainty in flow simulations using probabilistic methods[END_REF][START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF][START_REF] Smith | Uncertainty quantification[END_REF]. There are two main approaches to estimate aleatory uncertainties. One can use either statistical techniques directly or derive new models including the stochastic features of the parameters. The methods in the latter class are usually intrusive [START_REF] Ghanem | Stochastic finite elements: A spectral approach[END_REF][START_REF] Xiu | Numerical methods for stochastic computations: A spectral method approach[END_REF]. An important issue with all these methods is their high computational complexity. To overcome this limitation one usually introduces low-order models [START_REF] Obinata | Model reduction for control system design[END_REF][START_REF] Qu | Model order reduction techniques with applications in finite element analysis[END_REF][START_REF] Schilders | Model order reduction: Theory, research aspects and applications[END_REF] which are then used for UQ instead of the initial high-fidelity models.

Epistemic uncertainties, on the other hand, can be reduced by improving the modeling or solution procedures and it is one of our objectives to show how this can be done in the context of constrained coupled differential systems. Unlike classical numerical errors, epistemic uncertainties concern both the continuous and discrete levels. They can be reduced improving either the models or their discretizations. This paper concerns both issues as an initial model is first augmented before discretization. More precisely, we aim at addressing the uncertainties related to the numerical solution of differential systems which have constraints or in other words which are overdetermined. We propose a systematic and constructive approach which produces new equations whose solutions quantify the level of the uncertainties on the variables. It is based on a reformulation of the initial system by completing the system, i.e. adding appropriate differential consequences to the system. This means that we work on the continuous level before any discretization.

We expect the completed system to be more suitable for discretization because the constraints are then better taken into account. Hence typically some standard discretization method which gives poor results when applied to the original system will be accurate and robust when applied to the completed system. We can thus easily use generic commercial tools to solve our problem. Moreover our approach makes it easier to reach higher accuracy by increasing the accuracy of the discretization for a particular variable without being limited by possible compatibility issues of discretizations of different variables. For instance, we previously showed that the inf-sup condition (see e.g. [START_REF] Ern | Theory and practice of finite elements[END_REF]) for the Stokes problem disappears after our reformulation: the system stays stable even with higher order discretizations for the pressure than velocity [START_REF] Mohammadi | Simplifying numerical solution of constrained PDE systems through involutive completion[END_REF].

In general our aim is to incorporate the new features obtained by completing the given system as upgrades of existing solvers. In this way the considerable effort required in developing a new code from scratch is avoided [START_REF]Involutive upgrade of Navier-Stokes solvers[END_REF]. These techniques are not very well known in the numerical simulation communities but hopefully the present paper will help to popularize them further.

We will illustrate our approach by analyzing a well known model of packedbed chromatography. The physical idea of this interesting model is that there are several substances in the ambient fluid which are adsorbed on the boundary of the domain of the flow. This leads to a system where there are convection diffusion equations, ordinary differential equations and nonlinear algebraic equations.

The challenge of solving this kind of the system comes from the algebraic part of the system. Previous studies have treated the algebraic part of the problem as constraints and devised special methods such that the constraints would be as well satisfied as possible, see [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF][START_REF] Lim | Computer-aided model analysis for ionic strengthdependent effective charge of protein in ion-exchange chromatography[END_REF] and references therein. The drawback is that different systems may well require different methods. To address this difficulty we use our general methodology and differentiate the constraints to obtain new dynamic equations for the "remaining" or "algebraic" variables. Remarkably this will always give a (nonlinear) convection diffusion equation for any constraints. Hence we can write time evolution equations for all variables and use any standard method to solve them. Since we can monitor the constraints we can then correct the solution such that it stays on the constraint manifold.

We think that this is useful for epistemic uncertainty reduction of constrained systems because, unlike for the initial systems, all relevant constraints of the problem are explicitly included in the completed system. In this way the errors related to constraints can effectively be controlled and the computations thus become more stable because one does not need to worry about possible instabilities due to the violation of constraints.

From the completed system we derive a related system, called the augmented system, which contains some new variables. They play the role of local uncertainty estimators and require only a very low extra computational cost. These variables can also be estimated for the initial system, as well as for different discretizations and so they provide a natural framework to measure the evolution of the uncertainties in various cases.

To implement our approach we start with an existing program where the constraints are imposed by using the Newton's method. Typically there are problems of convergence. Then we modify the program taking into account the completed and augmented system that we derived from the initial system. In particular we use the dynamics of the new equations to initialize the Newton iterations. This improves greatly the convergence properties of the method.

We will also analyze the impact of our reformulation on aleatory uncertainties when the functioning physical parameters of the model are not known exactly, but their probability density function (PDF) is assumed to be known. In particular, resolving the constraints accurately has a clear impact on the prediction of the first and second order moments of the variables.

In order to quantify uncertanties it is very important to analyze how uncertanties propagate. We introduce small perturbations at the inlet boundary and study how they propagate through the domain. It is seen that, unlike the initial model, the augmented model damps high frequency components very fast, as theoretical considerations predicted.

In section 2 we recall the general background information about PDE systems which are needed in our formulation. In section 3 we describe the general form of the model and study its structural properties. In section 4 we introduce a particular model of 2 substances which was given in [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF][START_REF] Lim | Computer-aided model analysis for ionic strengthdependent effective charge of protein in ion-exchange chromatography[END_REF] and formulate it according to our approach. In section 5 we give the numerical results and discuss the impact of our formulation on the epistemic and aleatory uncertainties. In particular, we show that it is not always possible to separate these two sources of uncertainties.

Involutive, completed and augmented systems

Let us first briefly outline our general approach to constrained systems. For more details we refer to [START_REF] Dudnikov | Linear overdetermined systems of partial differential equations, Partial Differential Equations VIII[END_REF][START_REF] Krupchyk | Overdetermined elliptic systems[END_REF][START_REF] Krupchyk | Completion of overdetermined parabolic PDEs[END_REF][START_REF] Mohammadi | Simplifying numerical solution of constrained PDE systems through involutive completion[END_REF][START_REF]Involutive upgrade of Navier-Stokes solvers[END_REF][START_REF] Pommaret | Systems of partial differential equations and Lie pseudogroups[END_REF][START_REF] Seiler | Involution -the formal theory of differential equations and its applications in computer algebra[END_REF][START_REF] Spencer | Overdetermined systems of linear partial differential equations[END_REF]. Let us consider a system:

Ãu = f . (1)
This is the physical model we start with. This initial model contains in our case partial differential, ordinary differential and algebraic equations. Now differentiating these equations we obtain new equations which we call differential consequences or integrability conditions of the initial system. Evidently if u is a smooth solution then it will also satisfy all the differential consequences of the system. One may say that by computing differential consequences one makes explicit some information of the system which is only implicit in the original system. In particular one can thus analyze some structural properties of the system which is not possible using the initial system. Now it is a deep theorem in formal theory of PDE [START_REF] Dudnikov | Linear overdetermined systems of partial differential equations, Partial Differential Equations VIII[END_REF][START_REF] Pommaret | Systems of partial differential equations and Lie pseudogroups[END_REF][START_REF] Seiler | Involution -the formal theory of differential equations and its applications in computer algebra[END_REF][START_REF] Spencer | Overdetermined systems of linear partial differential equations[END_REF] that (under certain relatively mild conditions) it is sufficient to differentiate the system a finite number of times to obtain "full" structural information the system. The system thus obtained is called the involutive form of the given system.

Note in particular that in spite of the highly abstract nature of the theorem the implications of the theorem are quite constructive. To look for the involutive form one needs only to differentiate and eliminate variables which even in nonlinear case is fairly feasible using Gröbner basis techniques which are readily available in many computer algebra systems. It turns out that for the purposes of numerical computation it is sometimes convenient to use not the full involutive form of the system, but to add just some of the integrability conditions to the original system. Hence we will use the term completed system to indicate that we may not use the full involutive system.

In all cases completed systems have more equations than unknowns. One possible way to deal with this is to solve the system in the least square sense. Another approach which we advocate here is to make the system square again by introducing new variables which will also play the role of error and uncertainty estimators.

Let us denote by A 0 the operator obtained form à by adding some appropriate differential consequences to it and consider the completed system:

A 0 u = f. (2)
The right hand side f is obtained from f in the process of computing the differential consequences. However, since A 0 is overdetermined and has more equations than unknowns there are usually no solutions if f is chosen arbitrarily. Hence our problem (2) has no solutions unless f satisfies some compatibility conditions. These conditions are given by a compatibility operator A 1 such that A 1 A 0 = 0 and (2) has a solution only if A 1 f = 0.

Since A 1 contains important information about our system we will use it in the numerical solution. Let us now introduce some function spaces V i such that A i : V i → V i+1 and consider the following complex:

0 / / V 0 A 0 / / V 1 A 1 / / V 2 / / 0.
In the ideal situation this complex is exact; i.e. A 0 is injective, A 1 is surjective and image(A 0 ) = ker(A 1 ) . This suggests that we could decompose V 1 as:

image(A 0 ) ⊕ image(A T 1 ) = V 1 ,
where A T 1 is the formal transpose of A 1 . This decomposition means that it is indeed possible to find some functional framework in which the combined operator (A 0 , A T 1 ) is bijective. So instead of trying to solve the system (2) in some least square sense, we introduce an auxiliary variable z and solve

A 0 y + A T 1 z = f. (3)
We call this system the augmented system. The augmented system is square.

In previous works, we have shown how to proceed with this approach for quite general coupled differential systems (see [START_REF] Mohammadi | Simplifying numerical solution of constrained PDE systems through involutive completion[END_REF][START_REF]Involutive upgrade of Navier-Stokes solvers[END_REF][START_REF]Involutive formulation and simulation for electroneutral microfluids[END_REF] for examples of applications on constrained physical systems such as the Stokes and Navier-Stokes equations and coupled systems for microfluidic ionic flows). In this work we will illustrate the approach with a chromatography model.

Application to a coupled chromatography model

We consider a fluid flow where there are several substances whose concentrations interest us. The concentrations of the substances are supposed to be so small that the flow can be considered to be independent of them. Further it is realistic to suppose that the flow is one dimensional and even constant. We denote by x the coordinate of the domain along the flow and by v the velocity field.

Let C j be the concentration of the substance j in the fluid phase, S j concentration in the solid phase and E j the equilibrium concentration. We suppose that there are substances and denote C = C 1 , . . . , C and similarly for other variables.

First we have identical differential equations for the components of C and S.

C t -d C xx + v C x + αk(E -S) = 0 S t -k(E -S) = 0 (4)
Here d, α and k are some positive constants. In what follows we assume that v is constant. Hence we have simple convection diffusion equations for C and ODE for S. To close the system we have in addition an algebraic equation

g(E, C) = 0 , g : R 2 → R (5)
To get more information about the time evolution of E we differentiate the above equation. In the following ∂ t denotes the total derivative with respect to t and ∂ E g is the differential of g with respect to variables E. The second differentials are denoted similarly. This gives

∂ t g = ∂ E g E t + ∂ C g C t = 0 ∂ x g = ∂ E g E x + ∂ C g C x = 0 ∂ xx g = ∂ x ∂ E g E x + ∂ C g C xx + (∂ EC g)(E x , C x ) + (∂ CC g)(C x , C x ) = 0
Using these and the original equations ( 4) we obtain an evolution equation for E:

∂ E g E t -d ∂ x ∂ E g E x + v ∂ E g E x -d (∂ EC g)(E x , C x ) -d (∂ CC g)(C x , C x ) -αk∂ C g E -S = 0 (6)
Note that in addition to smoothness we did not need to assume any particular property of g to derive this equation.

Of course this is highly nonlinear in general. However, we have the following necessary condition for the equation ( 6) to be reasonable. Lemma 1. Suppose that C and S are known and ∂ E g is invertible. Then the linearization of (6) is parabolic.

Proof. With the above hypothesis the linearized equation is of the form

E t -d E xx + lower order terms = 0
The hypothesis that ∂ E g is invertible is perfectly reasonable since this guarantees by implicit function theorem that the variables E could actually be solved in terms of C using the constraint equations [START_REF] Ghanem | Stochastic finite elements: A spectral approach[END_REF]. In other words if this hypothesis were not satisfied then that would cast some doubt on the validity of the model itself.

It is somewhat curious that otherwise the lemma is valid whatever the actual constraint map g is. Even the diffusion constant is always the same for both E and C.

Our new system can now be written as

C t -d C xx + v C x + αk(E -S) = 0 S t -k(E -S) = 0 ∂ E g E t -d ∂ x ∂ E g E x + v ∂ E g E x -d (∂ EC g)(E x , C x ) -d (∂ CC g)(C x , C x ) -αk∂ C g E -S = 0 g(E, C) = 0 (7)
Let us finally denote u = (C, S, E) and write the above system simply as A 0 u = 0. Let us then introduce the operator

A 1 = -∂ C g, 0, -1, ∂ t -d ∂ xx + v ∂ x
This operator codes the calculations by which we found the new equation ( 6) and hence A 1 A 0 = 0. We call A 1 the compatibility operator of A 0 . Now using the same idea as in our previous papers [START_REF] Mohammadi | Simplifying numerical solution of constrained PDE systems through involutive completion[END_REF][START_REF]Involutive formulation and simulation for electroneutral microfluids[END_REF] we introduce the augmented system A 0 u + A T 1 z = 0 where z are new variables. Writing this out in full gives

C t -d C xx + v C x + αk(E -S) -(∂ C g) T z = 0 S t -k(E -S) = 0 ∂ E g E t -d ∂ x ∂ E g E x + v ∂ E g E x -d (∂ EC g)(E x , C x ) -d (∂ CC g)(C x , C x ) -αk∂ C g E -S -z = 0 z t -d z xx -v z x + g(E, C) = 0 (8)
Note in particular that convection for z is "backwards". For the exact solution z would be identically zero so we can use it to quantify and analyze the numerical error. Note that formally z is a solution of the equation

A 1 A T 1 z = z tt -2dz xxt + d 2 z xxxx -v 2 z xx + (I + ∂ C g(∂ C g) T )z = 0 Now I + ∂ C g(∂ C g) T is
symmetric and positive definite so let us perform the linear stability analysis in the scalar case and consider

z tt -2dz xxt + d 2 z xxxx -v 2 z xx + b z = 0
where b > 0. Looking for the plane waves one finds that there are solutions of the form z(x, t)

= exp i(kx - √ b + v 2 k 2 t) -dk 2 t (9)
Hence the high frequency components are strongly damped but low frequency components can propagate. We might therefore expect similar behavior also for the augmented model. This would also imply that the physical model itself is quite robust with respect to small high frequency perturbations of the concentrations.

two substances

Let us now consider a specific model with two substances (salt and a certain protein) [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF][START_REF] Lim | Computer-aided model analysis for ionic strengthdependent effective charge of protein in ion-exchange chromatography[END_REF]. Hence there are two algebraic equations g 1 = g 2 = 0 which we now describe. In the following all β i are some positive constants. First we have a simple equation

g 2 = E 2 + β 1 E 1 -β 0 = 0
The map g 1 is more complicated. In general we can write it as follows:

g 1 = K 12 γ 2 C 2 γ1 E 1 β 2 -γ2 β 2 E 2 γ 1 C 1 β 2 = 0
Here K 12 is the equilibrium constant, γ j are liquid activity coefficients and γj are surface activity coefficients. Now we somewhat simplify the model given in [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF][START_REF] Lim | Computer-aided model analysis for ionic strengthdependent effective charge of protein in ion-exchange chromatography[END_REF] by taking γ 2 = γ1 = γ2 = 1 and γ 1 = exp(β 4 C 1 ). This gives

g 1 = K 12 C 2 (E 1 ) β 2 -β 2 E 2 (C 1 ) β 2 exp(β 3 C 1 )
where β 3 = β 2 β 4 . Now computing ∂ E g we can check that det(∂ E g) > 0 as long as values of E and C remain strictly positive. Since g 2 is simply a linear equation we can eliminate the variable E 2 from the system. Let us now write the equations in this case and for simplicity of notation we denote E = E 1 and note that now z is scalar. Let us further denote ψ(E, S) = E -S 1 , β 0 -β 1 E -S 2 . Then the initial system we can be written as

C t -d C xx + v C x + αk ψ = 0 S t -k ψ = 0 f (E, C) =K 12 C 2 E β 2 -β 2 (β 0 -β 1 E)(C 1 ) β 2 exp(β 3 C 1 ) = 0 (10)
Here f : R 3 → R and we denote M = f -1 (0). Hence the solution should be such that it stays on the surface/manifold M . Note that M is smooth for positive values of variables but it is singular when E = C 1 = 0. Thinking M as a graph over (C 1 , C 2 ) plane for C j positive it is seen that 0 < E < β 0 /β 1 , see Figure 1.

One can say that M has 4 distinct regimes:

(i) the low plateau M 0 where E ≈ 0 (ii) the high plateau M 1 where E ≈ β 0 /β 1 (iii) the slope M 2 where one rises steeply from M 0 to M 1 along C 1 direction.

(iv) the wall M 3 where C 2 ≈ 0

The curvature of the surface is very big when one moves from one regime to another. One may summarize this by saying that the nonlinearities of M are very strong but geometrically the shape of M is quite simple. When we specialize the equation ( 6) to the case where we have only one constraint equation it is helpful to introduce the following notations

f E = ∂f ∂E , ∇ C f = ∂f ∂C 1 , ∂f ∂C 2 , f ij = ∂ 2 f ∂C i ∂C j
Note that f 12 = f 22 = 0. Then we can write the new equation as

f E E t -d ∂ x f E E x + f E vE x -d f 11 (C 1 x ) 2 -αk ∇ C f, ψ = 0 where v = v - d ∇ C f E , C x f E Figure 1.
The constraint surface of the model [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF].

Hence the completed system is

C t -d C xx + v C x + αk ψ = 0 S t -k ψ = 0 f E E t -d ∂ x f E E x + f E vE x -d f 11 (C 1 x ) 2 -αk ∇ C f, ψ = 0 f (E, C) = 0 (11)
Let us again denote the above system simply as A 0 u = 0. The compatibility operator is thus

A 1 = -∇ C f, 0, -1, ∂ t -d ∂ xx + v ∂ x
This gives the augmented system

C t -d C xx + v C x + αk ψ -∇ C f z = 0 S t -k ψ = 0 f E E t -d ∂ x f E E x + f E vE x -d f 11 (C 1 x ) 2 -αk ∇ C f, ψ -z = 0 z t -d z xx -v z x + f (E, C) = 0 (12)

Numerical results

The aim of the simulations that follow is to explore three issues. First, we illustrate how to use the extra information contained in the augmented system in a given simulation program. Moreover we show how this in turn gives a much more accurate solution than the one obtained with the original program. Then we discuss the impact of our new formulation on epistemic uncertainties (related both to the modeling and to the numerical scheme). Finally we give an example of the impact on aleatory uncertainties in a situation where a parameter is uniformly distributed over an interval.

In the description of the model problem all quantities related to length are in centimeters, time is in seconds and concentrations are moles per litre. We consider a one dimensional domain of length 1. The constants of the model are set as in [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF][START_REF] Lim | Computer-aided model analysis for ionic strengthdependent effective charge of protein in ion-exchange chromatography[END_REF]:

α = 1.38 k = 6.0 β 0 = 12.0 β 1 = 0.0334 β 2 = 4.91 β 3 = 2.51 K 12 = 1.05 • 10 10 .
5.1. Basic setting of computations. Suppose we have a program solving the system [START_REF] Lim | A novel partial differential algebraic equation (PDAE) solver: iterative spacetime conservation element/solution element (CE/SE) method[END_REF]. We denote by h the space discretization parameter and the time step is given by δt = 0.2h 2 /(d+h|v|). The program uses explicit Euler scheme for the time discretization and standard finite differences for space discretization. To handle transport related stability issues the advection term has been discretized using a second order upwind scheme. This introduces a numerical viscosity which is bounded by 0.5h |v|. We have done the simulations with two different meshes corresponding to h = 10 -2 and h = 5 • 10 -3 . The velocity and the diffusion coefficients were chosen as v = 1 and d = 0.1. Thus for the coarser mesh the numerical viscosity is less than 0.005. Despite the fact that this is quite small compared to d it appears to be necessary to use an upwind scheme to stabilize the calculations. The initial concentrations are set to zero. Since v > 0 we must specify the inlet values at the left boundary for C 1 and C 2 . These were chosen as C 1 in = 1 and C 2 in = 0.3. Since z should be as small as possible it has been initially set to zero and z in = 0 and (z x ) out = 0. 5.2. Solution of the initial system. We compute z even for the initial model, without otherwise using its value when computing the values of other variables. In this way z can be seen as an error indicator as well as an uncertainty indicator even for the initial system (see section 5.4). This leads to the following discrete form of the system (10):

1 δt (C n+1 h -C n h ) -(d + 0.5h|v|) C n xx,h + v C n x,h + αk ψ n h = 0, 1 δt (S n+1 h -S n h ) -k ψ n h = 0, E n+1 h ← f (E n+1 h , C n+1 h ) = 0, 1 δt (z n+1 h -z n h ) -d z n xx,h -v z n x,h + f (E n+1 h , C n+1 h ) = 0, (13)
Since E variable has no natural time evolution equation in the initial system we must use the algebraic equation f (E, C) = 0 in [START_REF]Involutive upgrade of Navier-Stokes solvers[END_REF] to compute E at the next time instant. Note that efficient and accurate solution of the the algebraic equation is the main issue with model ( 13) since the time evolution of the differential equations is standard. Our aim is to show how to use our augmented model to remove this difficulty. The original program uses Newton iterations at each time and for each mesh point to compute E:

E n,p+1 h = E n,p h -f E (E n,p h , C n+1 h ) -1 f (E n,p h , C n+1 h ), E n+1 h = E n,∞ h , f (E n,p h , C n+1 h ) = K 12 C 2,n+1 h (E n,p h ) β 2 -β 2 (β 0 -β 1 E n,p h )(C 1,n+1 h ) β 2 exp(β 3 C 1,n+1 h )
As stopping criteria it uses the value of the constraint and a bound on the number of iterations in p.

5.3.

Using the augmented system. The basic idea is to take advantage of the dynamical equation for E given in [START_REF] Mohammadi | Simplifying numerical solution of constrained PDE systems through involutive completion[END_REF] and use it to provide a better initialization for the Newton iterations. Hence given C n+1 , E n and S n+1 we first compute an intermediate value Ẽn which we then use as an initial value for Newton iteration if |f | > tol:

Ẽn h = E n h +δt f n+1/2 E,h -1 d ∂ x f n+1/2 E,h E n x,h -f n+1/2 E,h vn+1/2 h E n x,h +d f 11 (C 1,n+1 x,h ) 2 + αk ∇ C f n+1/2 h , ψ n+1/2 h + z n h
, where

f n+1/2 h =f (E n h , C n+1 h ) ψ n+1/2 h = E n h -S 1,n+1 h , β 0 -β 1 E n h -S 2,n+1 h vn+1/2 h =v -d ∇ C f n+1/2 E,h , C n+1 x,h /f n+1/2 E,h . ( 14
)
But equation ( 14) involves several nonlinear expressions and brings new difficulties. Recall that z h can be considered as a local measure of the error. Hence neglecting the nonlinear terms would typically increase greatly the value of |z h |. However, ( 14) is only used as an initialization to the Newton's method so we expect that |z h | is much reduced if the iteration converges successfully. Hence it is conceivable that one can simplify ( 14) considerably, and yet maintain the good convergence properties. We propose the following formula as a candidate for such a compromise:

Ẽn h = E n h + δt d E n xx,h -(f n+1/2 E,h ) -1 z n h . ( 15
)
Indeed it turns out that with this choice the Newton iterations always succeed.

In Figure 2 we show the ratio of the number of failures to the number of the discretization points at each time step for two meshes. Here by failure we mean that after the prescribed number of iterations the convergence criterion |f | ≤ tol was not fulfilled. To pursue the computation in these cases we used the final iteration of the Newton procedure. One can therefore deduce that in our discrete solution of the initial model the constraint f (E, C) = 0 is far from being adequately satisfied. It is also interesting to investigate how much the Newton iteration changes the computed solution. In Figure 3 we compare the solutions obtained by using ( 15) without projection to the case where ( 15) is used to initialize the Newton iterations. It is seen that the projection is necessary in the regions where the solution changes rapidly. One remarks that the error decreases with mesh refinement only in regions where the solution is smooth. Mesh refinement will not therefore remove the need for this final projection. Actually the need for projection increases for finer meshes in area of rapid variations.

In Figure 4 the computed solutions are plotted both for the initial and augmented model. The variables C 1 , S 1 , E and z are shown at times t = 0.3 and t = 0.6 over the domain for the two meshes. Improving the solution of the constraint equation by an initialization with our new evolution equation for E reduces the mesh dependency of the solution.

5.4. New variable z as a measure of epistemic uncertainties. The variable z is an indicator of the epistemic uncertainty both at the continuous and the discrete levels. It can be used, for instance, to measure the impact of mesh refinements on the solution when a reference solution is not available. In Figure 4 one can see that some of the deviations between the solutions of augmented and initial systems decrease with mesh refinement. However, some qualitatively incorrect features remain present in the solution obtained by the initial model such as jumps and non-monotonic behaviors. Hence, another interest of the approach is that mesh independence of the solution is achieved much faster when the constraint is better solved. The pertinence of these solution is confirmed by z distributions. 5.5. Impact on aleatory uncertainties. Now let us see the impact of our reformulation on aleatory uncertainties. These are irreducible and related to uncertain data. Let us first consider an interval-based simulation where some of the physical parameters are given as probability distributions on an interval and not as a single value. For instance, we can consider that the transport velocity v is a random variable which is uniformly distributed on the interval [0.8, 1.2].

The initial and the augmented models have been simulated on uniform 20 and 40 points sampling of this interval on the 100 points mesh. This coincides with a classical first order numerical integration over a uniform discretization of interval [0.8, 1.2].

We are interested in mean values and standard deviations of the variables at different times over the domain. These quantities are of importance in robust engineering design which is often based on the control of the first and second order moments.

To see the impact of the sampling size on the numerical evaluation of the moments we consider two different samplings of size 20 and 40. In Figure 5 these cases are compared on the 100 points mesh. One sees that a 20 points discretization of the interval is nearly sufficient for an accurate numerical integration and evaluation of the two moments. We see no difference on the mean and slight overestimation of the standard deviation with the 20 points sampling. This is less problematic in robust design than an underestimation.

In Figure 6 we use 40 points sampling to compare the expected values IE(C 1 ) and standard deviations IE((C 1 ) 2 )-IE(C 1 ) 2 1/2 of C 1 at times t = 0.3 and t = 0.6 for the two models for the 100 and 200 points meshes. One sees that imperfect resolution of the constraint leads to a large overestimation of the second order moment of the variable. Most importantly, the error reaches its maximum nearly where the first moments by the two methods coincide. This overestimation of the second order moment is even more visible for E and S in Figures 7 and8. One also remarks that the distributions of the moments are nearly the same for the two meshes when the constraint is well resolved. This therefore gives more confidence on these moments since they are less sensitive to the quality of the discretization than the results obtained with the initial model.

Prediction of higher order moments and extreme values is of major importance for most engineering applications and one sees here an example of how reducing the epistemic uncertainties leads to better estimations of aleatory uncertainties. 5.6. Propagation of perturbations from inlet boundaries. Another classical source of aleatory uncertainty is the imprecise nature of the boundary conditions. For instance, instead of being deterministic as above, the inlet values (at the left boundary of the domain since v > 0) for C 1 and C 2 can have a time dependent random component. As an example we consider the case where the inlet values are normally distributed:

C 1 in = N (1, σ 1 ) and C 2 in = N (0.3, σ 2 )
. This modeling of uncertain inlets is reasonable as the signal to noise ratio is usually known in practice. We consider here a signal to noise ratio (SNR) of 5%. This introduces a high frequency component in the solution from the inlet boundary. The plane waves computed in (9) make us expect that the high frequency components of the augmented system are heavily damped. Therefore the uncertainties introduced through the boundary conditions would disappear when propagating in the domain. Figure 9 shows the distributions for C 1 , E and S 1 in the domain at times 0.3 and 0.6 for 20 independent simulations with the randomly perturbed boundary conditions for C 1 and C 2 . We can see indeed that the perturbations do not propagate and remain very much located near the inlet left boundary. This is particularly clear for the variable C 1 . A nice indication of this behavior is through z profiles which, in particular, show drops where the variables E and S 1 have jumps. In comparison the indicator z, also estimated for the initial model (see lower curves in Figure 4), does not present any such behavior.

Conclusion

We have discussed the quantification of epistemic uncertainties in the context of numerical solution of constrained coupled systems. We have shown how our constructive algorithm for the derivation of the involutive and augmented systems provides a framework for a systematic estimation of these variabilities. The approach has been illustrated on a coupled set of partial differential, ordinary differential and algebraic equations modeling adsorption chromatography. The key idea has been to complete the system so that all variables have natural evolution equations. In particular, it has been shown that the approach naturally introduces extra variables which also play the role of epistemic uncertainty estimators.

Working with existing software, our aim has been to introduce minimal changes in the code. It has been shown how a suitable initialization for Newton iterations used for the solution of the constraint and a control of the error through the additional variable of the augmented model improves both the robustness of the numerical scheme and the accuracy of the solution obtained.

The impact of the approach on the quantification of aleatory uncertainties has also been discussed by estimating the first and second order moments in a situation when only the probability distribution of some parameters are known. It has been shown that the prediction of the standard deviation can be quite inaccurate when the constraint in the model is not properly solved. This is of major importance in engineering applications as a robust design of a system is often based on the minimization of the mean value of a functional under a constraint on its second order moment.

Finally, the impact of uncertainties in boundary conditions has been discussed. Numerical results recover the properties predicted by the theory for the augmented model where one expects fast damping of the high frequency components of the solution introduced through small perturbations at the boundaries. This behavior is well illustrated by the indicator z which clearly behaves differently for the initial and augmented systems.

From a practical point of view, this shows that chromatographic adsorption techniques are robust with respect to small and high frequency perturbations of the species concentration and therefore cannot be used to discriminate between two samples if those are only high frequency perturbations of each other. x (cm) 'DIFF_F_200' 'DIFF_F_100'

Figure 3. On the left: differences between Ẽn h given by ( 15) and its projection at times t = 0.3 and t = 0.6 on the meshes with 100 and 200 points. On the right: corresponding variations in f .

Figure 2 .

 2 Figure 2. Proportion of discretization points where the Newton algorithm fails for the initial system (10) for the 100 and 200 (lower curve) points meshes. Iterations succeed when initialized with E from (12).

  f(E_invol) -f(E_invol_project) | )

  0

Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 56789 Figure 5. Impact of the sampling size (continuous: 40, dashed: 20) on the mean value (left column) and the standard deviation (right column) for C 1 , E and S at times t = 0.3 and t = 0.6 on the 100 points mesh for the involutive model.