Purushothama Kolar

Laurent Naveen

Emmanuel Massoulié

Aline Baccelli

Don Carneiro Viana

Towsley

On the Interaction between Content Caching and Request Assignment in Cellular Cache Networks

determine routing mechanisms and associated cache update rules that together achieve a target trade-o↵ between cachemiss probability and bandwidth loads.

Related Work:

The problem of caching at femtocell base stations or heterogeneous caches has been recently studied in [START_REF] Dehghan | On the Complexity of Optimal Routing and Content Caching in Heterogeneous Networks[END_REF][START_REF] Golrezaei | FemtoCaching: Wireless Video Content Delivery Through Distributed Caching Helpers[END_REF][START_REF] Poularakis | Exploiting User Mobility for Wireless Content Delivery[END_REF]. Authors in [START_REF] Golrezaei | FemtoCaching: Wireless Video Content Delivery Through Distributed Caching Helpers[END_REF] address content placement at the femtocell base stations to maximize the probability of cache hit. A similar setup, but with mobile users, is studied in [START_REF] Poularakis | Exploiting User Mobility for Wireless Content Delivery[END_REF]. In both these works, request routing is purely based on content availability (so that cache hits are maximized), while ignoring the bandwidth costs at the base stations. Further, the authors focus on static content placement algorithms. For instance, in [START_REF] Golrezaei | FemtoCaching: Wireless Video Content Delivery Through Distributed Caching Helpers[END_REF] the content placement problem is formulated as that of maximizing a monotone submodular function under matroid constraints; a static greedy algorithm that is within a constant factor of the optimal is proposed.

Borst et al. in [START_REF] Borst | Distributed Caching Algorithms for Content Distribution Networks[END_REF] address the problem of content placement and routing to optimize the total bandwidth cost, but their formulation is limited to linear cost functions, in contrast to our work where we allow for general convex costs. Moreover, the focus in [START_REF] Borst | Distributed Caching Algorithms for Content Distribution Networks[END_REF] is again on designing static algorithms. Similarly, there are other works in the literature addressing the problem of static content placement so as to optimize some linear system cost [START_REF] Baev | Approximation Algorithms for Data Placement Problems[END_REF][START_REF] Dai | Collaborative Hierarchical Caching with Dynamic Request Routing for Massive Content Distribution[END_REF][START_REF] Jaho | Social Similarity Favors Cooperation: The Distributed Content Replication Case[END_REF].

There is an extensive literature on online cache management algorithms and in particular on LRU (Least Recently Used) caches [START_REF] Bianchi | Check Before Storing: What is the Performance Price of Content Integrity Verification in LRU Caching? SIGCOMM CCR[END_REF][START_REF] Dan | An Approximate Analysis of the LRU and FIFO Bu↵er Replacement Schemes[END_REF][START_REF] Fricker | A Versatile and Accurate Approximation for LRU Cache Performance[END_REF]. There is also an abundant literature on dynamic request routing [START_REF] Bonald | Insensitive Load Balancing[END_REF][START_REF] Jonckheere | Optimal Insensitive Routing and Bandwidth Sharing in Simple Data Networks[END_REF], but very few work address the design of online algorithms for both routing and cache management, and the interaction between the two. Notable exceptions are the works on cache networks [START_REF] Abedini | Content Caching and Scheduling in Wireless Networks With Elastic and Inelastic Tra c[END_REF][START_REF] Amble | Content-Aware Caching and Tra c Management in Content Distribution Networks[END_REF], which focus on stabilizing queues of pending requests by joint management of bandwidth and caches. In contrast we instead focus on loss-based models and on the objective of cost minimization rather than queue stabilization.

With a similar focus as ours, the authors in [START_REF] Jiang | Orchestrating Massively Distributed CDNs[END_REF] study the problem of cache replication and request forwarding in CDNs comprising distributed caches, e.g., household settop boxes. However, the objective is again limited to linear tra c costs. Moreover, they focus on scenarios with many servers with identical characteristics. This is in contrast with our work as we consider a system with few potentially heterogeneous caches; our limiting regime is instead obtained by scaling the arrival rates and the chunk size. There are other work in the literature, e.g., [START_REF] Leconte | Bipartite Graph Structures for E cient Balancing of Heterogeneous Loads[END_REF][START_REF] Tan | Optimal Content Placement for Peer-to-Peer Video-on-Demand Systems[END_REF], where proportional placement strategies and matching algorithms for routing requests to content available caches are studied. These di↵er again from our present work by focusing on scenarios with many identical servers.

SYSTEM MODEL

Consider a system comprising a finite collection of basestations (BSs) S and locations L where a location `2 L can be thought of as representing the set of all physical locations within the range of a subset of BSs. For instance, in Figure 1 the alternating black and grey shaded horizontal segments represents di↵erent locations. Thus, each location is connected to a subset of the BSs. These interconnections can be e↵ectively depicted as in Figure 2, where location ìs connected to BSs s and s 0 while `0 is connected to s 0 and s 00 . Let A denote the adjacency matrix representing the interconnections, i.e., A `s = 1 implies that location `is within the range of BS s; A `s = 0 otherwise. We assume a finite set C of contents in the system with the size of each content being one unit. We allow content chunking whereby a content can be completely reconstructed using k independently coded packets of the content. Each BS s 2 S includes a cache that can store up to Ms units of contents (memory constraints). We allow the BSs to store chunks of each content instead of complete copies. We assume that the chunks stored at di↵erent servers are always independent so that it only su ces to download k chunks from any subset of BSs to completely reconstruct the content.

Requests for content c arrive at location `according to a Poisson process with rate

(k) `c = k `c,
which can be thought as the aggregate rate for c from of all users stationed at `. The number of chunks of content c, available at server s, restricts the rate at which s can serve incoming requests for c (content-availability constraints). Content download requests, if accepted, result in service at unit rate. In this sense we consider a loss model (in contrast to a queuing model).

The bandwidth availability at a server s 2 S is reflected using a bandwidth-cost function, C (k) s (r), which represents the cost incurred by s for serving requests at a total rate of r. We assume that C (k) s , k 1, can be expressed as

C (k) s (r) = Cs( r k )
, where Cs, s 2 S, is strictly convex and increasing. We refer to the above system as the k-th system.

To obtain a benchmark for performance comparison, we consider a system where the caching and routing variables are relaxed to be reals. Let ✓sc 2 [0, 1] denote the fraction of content c stored at BS s, and kr s`c the rate at which s serves requests for c from `. Then, the memory constraint at BS s 2 S is given by P c ✓sc  Ms. Next, since only a fraction ✓sc of a request for c can be served by s, we have the following content availability (CA) constraints:

kr s`c  A `sk `c✓sc, s 2 S, `2 L, c 2 C
, where we have approximated the number of ongoing requests (of content c at `) by the steady state arrival rate, k `c (assuming an average service time of one). In order to obtain an unconstrained optimization formulation, we replace the above hard memory and CA constraints with strictly convex and increasing penalty functions, b Cs and C s`c , in the cost to be minimized 1 . b Cs(✓) denotes the cost of storing ✓ units of contents at BS s. For the CA constraint, denoting x := kr

s`c k `c✓sc, C s`c ( x k
) represents the cost of respecting (x negative) the CA constraint, or exceeding the constraint (x positive) by x units. Since the BSs are physically constrained to always satisfy the CA constraints, a positive x represents the content cache-miss rate from BS s for requests from `for c. Thus, C s`c ( x k ) represents the cost incurred for fetching the missed content rate from the main data center at the backend. Another possibility is to let C s`c ( x k ) = b max( x k , 0) for some b > 0, which can be interpreted as a cost of b per cache miss.

We now define the total cost incurred as,

C({✓sc}, {r s`c }) = P s Cs ⇣ P `c A `sr s`c ⌘ + P s b Cs ⇣ P c ✓sc ⌘ + P s`c A `sC s`c ⇣ r s`c `c✓sc ⌘ (1)
and propose the following optimization problem:

Total Cost (TC):

Minimize: C({✓sc}, {r s`c }) (2a)
Over:

✓sc 2 [0, 1], s 2 S, c 2 C (2b) r s`c 0, s 2 S, `2 L, c 2 C (2c)
Subject to:

X s A `sr s`c = `c, `2 L, c 2 C. (2d)
The equality constraint (2d) is used to ensure that all incoming requests are handled by the BSs. In Section 3, we propose an online algorithm to solve this optimization problem. We show that our algorithm converges to the optimal solution of (2) as the scale of the system in terms of arrival rate and number of chunks per content (i.e., k) increases.

ONLINE ALGORITHM

The main objective of this section is to propose an algorithm for assigning incoming requests to BSs based on their current "prices"; the key idea is that these prices reflect both the load as well as the content availability at the respective BSs (which is in contrast to the common practice of assigning requests purely based on either load or content Due to space constraints, we only briefly mention about the optimality and stability properties of our algorithm; theoretical details and proofs are available in our technical report [START_REF] Naveen | On the Interaction between Content Caching and Request Assignment in Cellular Cache Networks[END_REF].

Let ⇥sc(t) denote the number of chunks of content c stored at BS s at time t, and G slc (t) the total number of chunks of content c served (possibly for more than one request) from s to `at time t. Thus, ⇥sc(t)/k is the fraction of content c stored at s at time t, and R slc (t) := G slc (t)/k is the rate at which s serves requests for c from `. The chunk service times are i.i.d (independent and identically distributed) exponential random variables with unit rate.

Cache Updates: Updating a BS's cache for a content involves either adding new chunks or removing some existing chunks. In the k-th system, let k⌫sc denote the rate at which BS s updates the cache for content c. Suppose BS s initiates a cache update for c at time t. Then, the number of chunks that are added or removed is proportional to

⇥sc(t) = Qs(t) + X `A`s `cP s`c (t), (3)
where (using F 0 to denote the derivative of F ),

Qs(t) := b C 0 s ⇣ X c ⇥sc(t) k ⌘ (4)
P s`c (t) := C 0 s`c ⇣ R s`c (t) k `c ⇥sc(t) k ⌘ . (5)
Formally, ⇥sc(t) is updated as follows: ⇥sc(t + ) = ⇥sc(t) +[✏ ⇥sc(t)], where t + denotes the time instant just after t, ✏ > 0 is a constant, and the notation [x] represents the integer closest to x. We refer to Qs(t) and P s`c (t), respectively, as the memory price and the content availability (CA) price (w.r.t (`, c)) of BS s at time t.

Remark: For the sake of analysis, we assume that the caches are updated instantaneously at the update instants, although in practice a small time is required to download new chunks from a backend data center. We further assume that the backend data center always supplies a new batch of independently coded chunks, so that it is possible to reconstruct a content by downloading k chunks from any subset of BSs, without worrying about chunks being dependent.

Note that ⇥sc(t) is simply the negative of the gradient (w.r.t ✓sc) of the cost function in [START_REF]Cisco Visual Networking Index: Global Mobile Data Tra c Forecast Update[END_REF]. Thus, our cache update strategy is essentially the gradient descent algorithm. The novelty of our work appears next, where we propose a request assignment strategy that couples the update instants of {R s`c (t)} variables with the request arrival times. Further, incoming requests are assigned to download chunks from a (state dependent) selected subset of BSs, thus always increasing the current R s`c (t) values of the selected BSs. This di↵ers from the gradient descent algorithm where, depending on the gradient, an update could also recommend decreasing the value of some rate variables. We formally discuss our request assignment strategy next.

Request Assignments: Suppose a request for content c arrives at location `at time t. Then, the decision as to how many chunks should be downloaded from di↵erent BSs is based on the "prices" associated with the BSs connected to `. Formally, for each (`, c), we define the price of BS s (such that A `s = 1) at time t as, ⇧ s`c (t) = Ps(t)

+ P s`c (t) (6)
where, P s`c (t) is the CA price (recall (5)) and

Ps(t) := C 0 s ⇣ X `0c 0 A `0s R s`0c 0 (t) k ⌘
represents the bandwidth price of BS s.

Remark: Although the price depends on various caching and routing variables, its calculation in practice can be based on, for instance, an estimate of the current load at a BS to compute its bandwidth price; The content availability price can be simplified to reflect the cost of fetching the unavailable chunks from the backend data center. We will study the performance of such a light-weight scheme in Section 4. Now, given the BS prices, let S `c(t) denote the set of all minimum price BSs reachable from `, i.e.,

S

`c(t) = argmin

s n ⇧ s`c (t) : A `s = 1 o . (7)
Our request assignment strategy centers on downloading sufficient chunks of c from the BSs in S `c(t) to construct the content at `. However, if there are an insu cient number of chunks at these BSs, signals in the form of dummy chunk assignments are issued to increase their CA price and thus trigger future replication of content c chunks at these BSs. Formal details follow, where we have chosen to distribute the dummy chunks equally among all the BSs in S `c(t) (see (8a) and (8a)); however, note that any other dummy chunk assignment would work equally well, as long as the total number of chunks allocated (actual+dummy) is k.

Let n denote the number of BSs in S `c(t) (i.e., n = |S `c(t)|).

Now, arrange the BSs in S `c(t) in any arbitrary order, for instance, e.g., in increasing order of the BS IDs (assuming that each BS is given a unique integer ID). Suppose (s1, s2, • • • , sn) is such an arrangement, then BS si is assigned to serve

s i `c(t) = (1)
s i `c(t) + (2)
s i `c(t) chunks (ac-tual+dummy, resp.) of content c, where (using bxc to denote the largest integer smaller than x),

for i = 1, • • • , n, (1)
s i `c(t) = min n k i 1 X j=1 (1)
s j `c(t), ⇥s i c(t) o , (2)
s i `c(t) = 6 6 6 4 k P n j=1 (1)
s j `c(t) n 7 7 7 5 , i = 1, • • • , n 1 (2) sn`c (t) = k n X j=1 (1)
s j `c(t) n 1 X j=1 (2)
s j `c(t).
No chunks are downloaded from the BSs not in S `c(t), i.e., we have s`c (t) = 0 for s / 2 S `c(t).

Thus,

s i `c(t) (1)
is the actual number of chunks that si can serve, while

P n j=1 (2)
s j `c(t) represents the total number of chunks that are not present in the caches, i.e., results in cache misses. Although, physically it is not possible for si to serve (2) s i `c(t) number of additional chunks, we proceed to allocate these as dummy chunks served by si. As mentioned earlier, the idea is to increase the CA price P s i `c(t) (recall ( 5)), so that more chunks of c will be replicated by s during its subsequent cache updates (see [START_REF] Amble | Content-Aware Caching and Tra c Management in Content Distribution Networks[END_REF]). Practically, the dummy chunks are downloaded from the backend data center.

Finally, the associated rate variables are updated according to the following expression: R s`c (t + ) = R s`c (t)+ s`c (t). Theoretical Guarantees: Using mean field approximation and Lyapunov function techniques, we show that our online algorithm is optimal and stable in a limiting fluid regime that is obtained by scaling k, i.e., the arrival rates and the number of chunks per content; formal details are available in our technical report [START_REF] Naveen | On the Interaction between Content Caching and Request Assignment in Cellular Cache Networks[END_REF]). From a practical standpoint, our theoretical results imply that for a large system2 the performance of the online algorithm is expected to converge in time to a close-to-optimal performance.

SIMULATION EXPERIMENTS

We first consider a simple setting comprising 2 locations and 3 BSs; the content catalog size is |C| = 10. After demonstrating the e cacy of the joint-price based request assignment schemes for this setting, we present results for a larger system.

We assume that BSs 1 and 3 are femtocells serving locations 1 and 2, respectively, while BS 2 is a common cellular BS connected to both the locations (recall Figure 2). Thus, the requests arriving at a location can be assigned either to the respective femtocell or to the cellular BS. The total (unscaled) content request rate from each location is 1, i.e., P c `c = 1 for `= 1, 2. The content popularity at each location follows a Zipf distribution with parameter 0.8, i.e., the probability that an incoming request is for the i-th most popular content is proportional to 1/i 0.8 . In order to obtain a scenario where the same content can have di↵erent popularities at di↵erent locations, we randomly shu✏e the distributions for the two locations. The (unscaled) cache update rate of each BS for each content is 1, i.e., ⌫sc = 1.

The various costs (recall ( 1)) used in the simulations are as follows:

• Bandwidth cost: Cs(r) = exp(2(r Rs)), where R1 = R3 = 0.5 and R2 = 0.1; thus, the cellular BS is costlier than the femtocell BSs.

• Memory cost: b

Cs(✓) = exp(5(✓ Ms)) where M1 = M3 = 2 while M2 = 8; thus, the cellular BS has more memory compared with the femtocells.

• Content availability (CA) cost: C s`c (x) = exp(bx); we are particularly interested in studying the e↵ect of b on performance. A large value of b, which yields a steep cost function, results in a system where cache misses are costlier than bandwidth. Thus, b can be used to trade o↵ bandwidth cost and cache-miss probability. We refer to b as the cache-miss exponent.

Performance of the online algorithm: In Figure 3 we first plot the optimal performance curve, obtained by directly solving the problem in [START_REF] Abedini | Content Caching and Scheduling in Wireless Networks With Elastic and Inelastic Tra c[END_REF] for di↵erent values of b (curve labeled 'Optimal'). As highlighted in the figure, the cache-miss exponent b increases from left to right along all the curves. Also shown in the figure is the performance of the online algorithm, described in Section 3, for a system with k = 10 chunks per content. Although the online algorithm is proven to converge to the optimal performance as k ! 1, it is interesting to observe that its performance is comparable with the optimal already for k = 10.

Before proceeding further, let us briefly discuss the trade-o↵ curves. Note that, as b increases the cache-miss probability improves, however at the expense of an increased bandwidth cost. This is because, for a larger b, the cellular BS o↵ers a lower price for the incoming requests as it holds more contents. Thus, most of the requests are forwarded to the cellular BS, increasing its (and the overall) bandwidth cost. In contrast, when b is small, bandwidth cost is the key component in deciding the BSs' prices. Since femtocells are less expensive in this regard, more requests are now forwarded to the femtocells, resulting in more cache misses as they only hold few contents.

Drawbacks of the online algorithm:

The implementation of the online algorithm however requires the BSs to maintain routing variables in the form of number of chunks of each content downloaded by requests from each location. This would require identifying each chunk with its request's location, which is not possible in practice when the number of ongoing requests are large. Further, estimates of the request arrival rates, `c, are required by the BSs to compute their content availability prices. Finally, the cache update in the online algorithm is server-driven requiring the BSs to update their caches regularly. This would result in unnecessary updates whenever the request rate is less.

A light-weight signaling scheme: In view of the above drawbacks, we modify the online algorithm to obtain a lightweight signaling scheme for request assignments. We refer to it as the joint-price scheme since, as in the online algorithm, the request assignment here is based jointly on the content and bandwidth availability. The bandwidth price, as before, remains to be the same convex increase function of the current load on the BS; thus, by using the congestion signals from the BS, it is possible to estimate the bandwidth cost at a location. However, for the content availability price, we now simply use the cost of fetching unavailable chunks if an insu cient number of chunks of the requested content is available at the BS. Formally, the content availability price of BS s for request c (from any location) at time t is given by b(k ⇥sc(t))/k, where ⇥sc(t) is the number of chunks of content c at BS s at time t and b now actually represents the cost of fetching the unavailable fraction of content c. Thus, as with the online algorithm, a larger b implies that the cache misses are costly.

We use the above joint price signaling scheme in conjunction with the LRU (Least Recently Used) cache eviction algorithm [START_REF] Bianchi | Check Before Storing: What is the Performance Price of Content Integrity Verification in LRU Caching? SIGCOMM CCR[END_REF][START_REF] Dan | An Approximate Analysis of the LRU and FIFO Bu↵er Replacement Schemes[END_REF], which is a request driven cache management strategy. Under LRU, cache replacements are initiated by incoming requests that do not find the desired content in the cache (a cache-miss event). As a consequence, the least recently used content is removed from the cache to free memory for storing the new content. Since we deal with contents at the chunk level, we have implemented a variation of LRU, referred to as the CLRU (Chunk LRU) strategy. In CLRU, the incoming request at a cache for m chunks of a content triggers the relocation of such m chunks to the head of the cache queue to be served. The placement of extra cached chunks at the cache (i.e., if more than m chunks are cached) remain unchanged. In the opposite case, if chunks are missing (i.e., if less than the requested m chunks are cached) and the cache is full, new chunks of the content are downloaded and placed at the head of the queue by removing least recently used chunks at the tail of the queue.

Performance of joint-price scheme: In Figure 3 we have depicted the performance of the joint-price signaling scheme. Also shown in the figure are the performance of two heuristic policies (commonly used in practice) where the request assignment is purely based on bandwidth availability or content availability signals (points 'Load-only Signaling' and 'Content-only Signaling', respectively). The joint-price scheme provides a favorable trade o↵ between bandwidth cost and cache-miss probability, encompassing the heuristic policies as its extreme end points.

Although the joint-price scheme su↵ers some degradation in performance when compared with the online algorithm, it is useful to emphasize that the simplicity of the former scheme would render it more favorite than the latter when it comes to considerations for practical implementation. Moreover, the performance degradation is less for low to moderate values of cache-miss cost b, which is typically the case in practice.

Large network: Motivated by the observations made in the previous section, we proceed to study the performance of the joint-price scheme for a large network. We now consider a content catalog of size 1000. The number of locations are increased to 10. It will be useful to think of a location `2 {1, 2, • • • , 10} as the segment [0, `] on the real line (recall Figure 1 for illustration). A cellular base station is located at the center (i.e., at 5) so that, as in the earlier case, it can cover all the locations. However, we now increase the number of (small-size) femtocells to 9, which are located at positions 1, 2, • • • , 9; each small-size femtocell has a range of 1. Thus, each location (except the end locations) is served by two adjacent femtocells. To model the heterogeneity in the femtocell base stations, we introduce 5 moderate-size femtocells. The moderate femtocells have a larger range of 2 and are located at 3, 4, • • • , 7. Note that, the inner locations are covered by more femtocells (small+moderate) than the outer ones.

Small-size femtocells have smaller caches (of size 10) compared to moderate-size femtocells (whose cache capacity is 100), and the latter's capacity is smaller than the cellular BS which can store up to Ms = 800 contents. Bandwidth cost increases as we move from the small to moderate femtocell and to cellular BS; recalling the cost functions from the previous section, the respective values of Rs are 5, 2.5 and 1.

In Figure 4 we first plot the performance achieved by the joint-price signaling scheme along with the performance of the two heuristic schemes. Again, we observe that the jointprice scheme can achieve a large range of performance values, which is in contrast to the heuristic schemes whose perfor-mance is fixed at either optimizing bandwidth cost or cachemiss probability, without regarding for the other metric.

In the thus far implemented joint-price schemes, we have been assigning requests to download chunks strictly from the set of min-priced servers. We now relax this strict consideration by regarding any BS, whose price is within a tolerance value of ⌧ from the minimum price, as eligible to serve chunks for an incoming request (provided as before the eligible BSs are connected to the request's location); chunks are however downloaded by sorting the eligible BSs in the increasing order of their prices. In Figure 4 we have depicted the performance of the joint-price scheme for di↵erent values of ⌧ , where the earlier strict min-price approach now corresponds to ⌧ = 0; the case where all connected BSs (irrespective of their prices) are eligible corresponds to ⌧ = 1.

We observe that as ⌧ increases the trade-o↵ curves shift downwards, implying that it is possible to achieve a lower cache-miss probability for a given target bandwidth cost. However, the range of trade-o↵ that is possible reduces as ⌧ increases. For instance, using ⌧ = 100 it is not possible to reduce the bandwidth cost to less than 40; In fact, ⌧ = 1 case achieves the lowest cache-miss probability, however at the expense of being constrained to operate at a higher range of bandwidth cost. This shift in performance is because, as ⌧ increases, more BSs become eligible to serve a request so that the cache-miss probability decreases; however, more BSs are now loaded, resulting in a bandwidth scarce system. Thus, further trade-o↵ in performance can be achieved by suitably choosing the value of ⌧ .

In summary, the joint-price scheme, in conjunction with the CLRU policy for cache management, can be a practical candidate for request assignments in futuristic cellular networks where a dense deployment of heterogeneous femtocell base-stations are expected. The cache-miss cost b, along with the tolerance value ⌧ , can be used to serve as "tunable knobs" to trade-o↵ one metric for the other.

CONCLUSION

We proposed an optimization framework to study the problem of cost minimization in cellular cache networks. Towards this direction, we proposed an online algorithm for caching and request assignments which is shown to be optimal and stable in a limiting regime that is obtained by scaling the arrival rates and the content chunking. Based on the online algorithm, we proposed a light-weight joint-price scheme for request assignment that can be used in conjunction with the LRU cache management strategy. Through simulations we found that our joint-price based request assignment strategy outperforms the common practices of routing purely based on either load or content availability. Our proposed jointprice routing mechanisms are thus an appealing candidate combining sound theoretical guarantees with good experimental performance while being simple to implement.

Figure 2 :

 2 Figure 2: Illustration of the interconnections between BSs and locations. the range of BS s; A`s = 0 otherwise. We assume a finite set C of contents in the system with the size of each content being one unit. We allow content chunking whereby a content can be completely reconstructed using k independently coded packets of the content. Each BS s 2 S includes a cache that can store up to Ms units of contents (memory constraints). We allow the BSs to store chunks of each content instead of complete copies. We assume that the chunks stored at di↵erent servers are always independent so that it only su ces to download k chunks from any subset of BSs to completely reconstruct the content.Requests for content c arrive at location `according to a Poisson process with rate

Figure 3 :

 3 Figure 3: Performance curves.(unscaled) content request rate from each location is 1, i.e., P c `c = 1 for `= 1, 2. The content popularity at each location follows a Zipf distribution with parameter 0.8, i.e., the probability that an incoming request is for the i-th most popular content is proportional to 1/i 0.8 . In order to obtain a scenario where the same content can have di↵erent popularities at di↵erent locations, we randomly shu✏e the distributions for the two locations. The (unscaled) cache update rate of each BS for each content is 1, i.e., ⌫sc = 1.The various costs (recall (1)) used in the simulations are as follows:• Bandwidth cost: Cs(r) = exp(2(r Rs)), where R1 =

Figure 4 :

 4 Figure 4: Performance of the light-weight schemes for a large network.cost. In contrast, when b is small, bandwidth cost is the key component in deciding the BSs' prices. Since femtocells are less expensive in this regard, more requests are now forwarded to the femtocells, resulting in more cache misses as they only hold few contents.

Note that large here refers to a system with large arrival rates and content chunking. It does not imply a large network in terms of number of BSs and locations; all our analysis holds for a network of any finite size.


F Ribaud

M Ait-Ighil

S Rougerie

J Lemorton

O Julien

F Pérez-Fontan

Reduced multipath channel modelling preserving representative GNSS receiver testing

The increasing need for accuracy in GNSS positioning applications requires a thin knowledge of the impact of the urban propagation channel on the receiver's tracking loops. Therefore, new receivers' performance assessment requires the study of their behavior in a multipath environment. In that context, propagation channel synthesizers have to be able to provide multipath parameters that are representative of a real environment. Then, the receiver can be put to test by using these multipaths parameters to quantify its resulting positioning error.

In order to simulate a realistic environment, the channel synthesizer usually generates a high number of multipaths. This characteristic can prevent an optimal receiver performance evaluation, especially when using a hardware channel emulator, whose multipath production capacity does not normally exceed 10 units per satellite (even if the total number of channels can reach 36 or 72 in the most powerful emulator). Therefore, this study aims at elaborating reduction methods of the channel impulse response and comparing their capacity to preserve the receiver's positioning error induced by the multipath channel, in order to advocate the use of one of these methods to answer the multipath channel reduction problem.

I. Elaboration of 3 reduction methods

The first part of the study consists in the description of the considered reduction methods. They represent 3 different approaches to the problem:

A. Aggregation methods

First, the channel impulse response can be reduced by coherently summing (amplitude and phase) groups of multipaths. To do so, a clustering method has been implemented in order to create clusters of multipaths according to Delay and Doppler proximity criteria. The multipaths of the reduced channel are the centers of the clusters, located at the barycenter of the multipaths of the clusters. Such a reduction technique has been chosen for ease of application to any impulse response and its preservation of the multipaths power-delay spread and Doppler spectrum. Figure 1 illustrates an example of channel impulse response reduction using the Clustering technique, with 5 multipaths in the reduced channel.

B. Application of parametric methods

A second approach consists in optimizing the parameters of the reduced channel multipaths (delay, Doppler and amplitude) to minimize a cost function. Knowing that the receiver uses the correlation of the incoming GNSS signal with its local replica to compute the tracking loop discriminator, the cost function has been chosen to be the difference between this theoretical correlation function convoluted with the complete original multipath channel and the reduced one. The SAGE (Space Alternating Generalized Expectation maximization) iterative algorithm has been used to perform this optimization, like in [START_REF] Rougerie | Algorithmes de diversité d'antennes appliqués à la reception des signaux GNSS en environnement urbain et sur terminal mobile[END_REF]. The multipaths parameters output of the algorithm constitute the reduced channel impulse response. Figure 2 illustrates an example of channel impulse response reduction using the parametric technique, with 5 multipaths in the reduced channel. An example of this optimization realization is displayed on fig. 3, which represents the autocorrelation functions of the original complete multipath channel and reduced one. The Line Of Sight (LOS) contribution has not been represented to highlight the importance of the multipaths parameters.

C. Statistical approach

The last approach considered in this study computes the multipaths of the reduced channel as the realization of a stochastic process. The temporal evolution of the reduced channel multipaths delays and Doppler shifts follow a discrete state and time Markov channel, like in [START_REF] Schubert | Components Count Decrease of Channel Responses Using Markov Process[END_REF]. In other terms, the multipaths parameters depend at a given time on the present delay and Doppler statistic distribution and on the one at the precedent epoch to calculate the transition probabilities. The complex amplitudes of those contributions are drawn from the complete channel impulse response, knowing their delay. The multipaths of the reduced channel are realizations of this stochastic process. It can be remarked that contrary to the two previous reduction methods, this one is not deterministic and its performances can variate from realization to realization. Figure 4 illustrates an example of channel impulse response reduction using the statistical technique, with 5 multipaths in the reduced channel. Figure 5 illustrates the stochastic evolution of the delay of the multipaths of the reduced channel.

II. Study of the pseudo-range error preservation

The SCHUN (Simplified Channel for Urban Navigation) channel synthesizer [START_REF] Ait-Ighil | Enhanced Physical-Statistical Simulator of the Land Mobile Satellite for Multipath Modelling Applied to a Satellite Navigation System[END_REF] has been chosen for its channel accuracy/computation time ratio to generate a reference channel impulse response in a virtual environment in which the receiver moves (multipaths delay, Doppler and complex amplitude). This reference scenario is represented on fig. 6.The receiver moves along the blue line at a speed of 5 m.s -1 . A single satellite emitter has been considered located at azimuth 30° and elevation 35 °. The reduction methods have been applied to this reference scenario along the receiver's trajectory to compare their pseudo-range error preservation performances. The number of multipath has been reduced by the 3 reduction methods from several thousands in the reference channel to 5 units in one case and to 10 units in a parallel study. The performance of the different multipath reduction methods can then be compared according to several criteria.

A. Discriminator error preservation

First, the impact of the channel reduction process has been observed through its impact on the tracking loop's discriminator output. The principle of this criterion is to compute the deviation of the discriminator due to the multipath channel considering a given discriminator, Chip Spacing (CS) and signal modulation. This study offers a first look at the multipath channel impact on the receiver by looking at its impact on the tracking loop stable lock point (instead of the actual tracking loops). The resulting has been computed for a BPSK modulated signal along the receiver's trajectory for the reference impulse response and the different reduced ones, using an Early Minus Late (EML) Discriminator with 1 chip CS.

The comparison of the difference of static pseudo-range error between the reference channel and the different reduced channels has shown that the best reduction method according to this criterion is the optimization method (B). The clustering technique (A) preserves discriminator output error almost as good as the optimization method does. Indeed the difference between the reference channel discriminator error and the clustering reduced channel is 12 % higher than the difference between the reference channel and the optimization technique reduced one.

On the contrary, the statistic method (C) has the worst performances considering this criterion. Indeed the difference between the reference channel discriminator error and the statistical method reduced channel is 3 times higher than the difference between the reference channel and the optimization technique (B) reduced one. Those statements remain the same considering a reduced channel impulse response with 5 multipaths and 10 multipaths. Figure 7 represents the cumulative distribution of the difference between the discriminator error of the reduced channels and the original one along the receiver's path. Left plot corresponds to the 5-multipath reduced channels study and the right plot the 10-multipaths study.

B. Dynamic pseudo-range error preservation

To study the impact of the reduction method on a real receiver, a software receiver has been used to compute the receiver's pseudo-range error along its trajectory. In this section, the dynamic tracking of the GNSS signal code and carrier phase are investigated. The signal is BPSK modulated and an EML discriminator with CS=1chip has been used for the DLL. This subsection confirms the conclusion of the previous one concerning the hierarchy of the methods according to the preservation of their impact on the receiver: the optimization technique (B) has the best preservation performances, close to the Clustering (A) performances and far better than the statistical method (C) performances. This hierarchy among the channel reduction methods according to this static pseudo-range error criterion is the same with 5 or 10 multipaths in the reduced channels. Moreover, it is observed that the difference between the Clustering technique and the optimization technique performances decreases with the increase of the number of multipath in the reduced channel, tending to converge. Figure 8 represents the cumulative distribution of the difference between the dynamic pseudo-range error of the reduced channels and the original one along the receiver's path. Left plot corresponds to the 5-multipath reduced channels study and the right plot the 10-multipaths study.

III. Conclusion

This study has presented three different reduction methods, elaborated from three different approaches (multipath aggregation, channel autocorrelation function preservation and stochastic evolution of the multipaths parameters), aiming at covering a large field of possible solutions of the channel impulse response reduction problematic. The goal of these methods is to reduce the number of multipaths from several thousands to less than 10 units impacting the receiver's positioning performances as little as possible. Therefore, their validity has been investigated through the preservation of the receiver's behavior. The results of this comparison show a big difference between the Clusteringparametric methods (A-B) and the statistical method (C), which corrupts the pseudo-range error much more than the two others. As a consequence, this method can be eliminated from the possible approaches to answer the multipath channel reduction problem.

The performances of the Clustering (A) and optimization (B) techniques show to be very close for 5 multipaths and show disappear in the 10 multipaths study. It is however important to mention that the parametric technique (B) has a large computation time, whereas the Clustering technique (A) computation time is negligible as compared to the receiver's multipath processing time. Therefore, the use of the Clustering technique (A) to reduce the multipath channel can be recommended for most of the applications.

Finally, this study has used both static discriminator error and dynamic receiver's tracking error to assess the performance of the reduction methods. Both show the same conclusions. As a conclusion, the discriminator error, easy to compute, could be sufficient to evaluate the performances of a given reduction method in a given scenario. This tool can be a complement to the reduction method in itself, quantifying its error on the receiver's pseudo-range error preservation, playing the role of a representativeness index.

Figure 1 .

 1 Figure 1. Example of impulse response reduction with the Clustering technique

Figure 2 .Figure 3 .

 23 Figure 2. Example of impulse response reduction with the parametric technique

Figure 4 .

 4 Figure 4. Example of impulse response reduction with the statistic technique

Figure 6 .

 6 Figure 6. Reference scenario

Figure 7 .

 7 Figure 7. Cumulative distribution of the absolute difference between the discriminator error caused by the original channel and the reduced ones in case of 5 and 10 multipaths.

Figure 8 .

 8 Figure 8. Cumulative distribution of the absolute difference between the discriminator error caused by the original channel and the reduced ones in case of 5 and 10 multipaths.


S Rougerie

M Ait-Ighil

Vincent Fabbro

Impact of GNSS receiver tuning on the estimation of scintillation index

When crossing the ionosphere layer, GNSS signals may be impacted by scintillation, a dynamic effect which causes very fast variations of the amplitude and the phase. The GNSS receivers may lose lock, which lead to unavailability of the satellite link as well as the reduced reliability of the available information. In the present state of the art, the understanding of scintillation phenomena is still insufficient, and the available models are still not adequate for GNSS system design. Thus, many studies are carried out in order to model scintillation [1] [2] [3] and evaluate the effect on the GNSS receivers and augmentation system [4] [5]. These studies monitor the scintillation thanks to specific scintillation index. However, the estimation of these scintillation indices may be affected by the receiver itself and by the estimator. The goal of this paper is to understand how receiver or estimator tuning may impact the estimation of scintillation index. To conclude, this paper gives recommendation on how to tune GNSS receiver or estimator in order to minimize their impact.

I.

INTRODUCTION

When crossing the ionosphere layer, GNSS signals are impacted by two kinds of effects: mean effects and dynamic effects. Mean effects encompass the group delay and the signal phase advance occurring when the signal is going through the ionosphere. This delay depends on the TEC (Total Electron Content)of the ionosphere layer. It leads to a positioning error of several meters to tens of meters. Most of this delay can be corrected thanks to a dual frequency measurement, or using correction models (like Klobuchar, NeQuick or TEC grids in the case of augmentation systems) for single frequency measurements. Most of the dynamic effects are coming from the scintillation phenomena. These events generally occur over the equatorial and polar regions, but they may also occur at middle latitudes when the solar activity increases as for solar maximum. These events impact the signal and cause very fast variations of the amplitude and the phase. Therefore, the signal fades and GNSS receivers may lose lock. As a consequence, the main impact of the ionosphere dynamic effects on GNSS is the unavailability of the satellite link as well as the reduced reliability of the available information. In the present state of the art, the understanding of scintillation phenomena is still insufficient, and the available models are still not adequate for GNSS system design.

Thus, many studies are carried out in order to model scintillation [START_REF] Boscher | Modelling ionospheric effects for L band GNSS receivers at high latitudes[END_REF] [2] [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF] and evaluate the effect on the GNSS receivers and augmentation system [START_REF] Secretan | SAGAIE A GNSS Network for Investigating Ionospheric Behavior in Sub-Saharan Region[END_REF] [START_REF] Béniguel | MONITOR IONOSPHERIC MONITORING SYSTEM: GNSS PERFORMANCE ESTIMATION[END_REF]. These studies monitor the scintillation thanks to specific indices such as ROTI (Rate Of TEC Index), S4, SigmaPhi or the slope of the scintillation spectrum. However, the estimation of these scintillation indices may be affected by the receiver itself. Indeed, ROT (Rate Of TEC) spectrum and phase spectrum are estimated thanks to the PLL and DLL outputs, which are directly affected by the loop bandwidth and loop order. Thus, the scintillation indices values, which are the results of the spectrum integration (for example ROTI is obtained by integrating ROT spectrum), may depend of the receiver tuning, and also of the integration parameters (i.e. integration time and sampling frequency).

The goal of this paper is to understand how receiver / estimator tuning may impact the estimation of scintillation indices. The approach is the following: based on a theoretical model of scintillation [START_REF] Report Itu-R P | Transionospheric radio propagation, The Global Ionospheric Scintillation Model (GISM)[END_REF] [3], phase and amplitude time series are generated. These times series feed a GNSS receiver simulator, which at the end generate the traditional GNSS observable: pseudo-range and phase estimations. Based on these observables, conventional scintillation estimators are applied. The simulated and the estimated scintillation spectrum (and thus the scintillation indices by integration) are finally compared in order to evaluate the GNSS receiver influence. The paper is therefore organized as follow: will be presented in section II, a short review of the scintillation model, the GNSS simulator and the scintillation algorithms estimator. In the section III, the entire tests will be presented in order to evaluate the impact of the receiver on the phase and amplitude scintillation spectrum and indices. Finally, the section IV summarise all the receiver and estimator effects, give recommendation to minimise it, and present our conclusions and future work.

II.

Scintillation simulator

For the purpose of this study, a scintillation simulator has been developed. This simulator is divided in three main steps: a scintillation model, a GNSS receiver simulator, and a scintillation estimator.

A. Scintillation model

The electronic density Ne , , with the vector path and t the time, is modelled as a stochastic process, where the mean effect 〈 , 〉 may be given by a traditional ionosphere model, such as NeQuick or Klobuchar, and the dynamic one, reproducing the scintillation due to turbulence in the electronic density, given by a centred random process Δ , .

, 〈 , 〉 Δ , (1)
The stochastic process Δ , is model thanks to the Schkarofsky spectrum. The link between this spectrum and the amplitude / phase fluctuations of an electromagnetic waves propagating through such turbulent environment is given in [START_REF] Report Itu-R P | Transionospheric radio propagation, The Global Ionospheric Scintillation Model (GISM)[END_REF] [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF]. At the end, we define the generic scintillation spectrum of the log amplitude and phase as follow:

(

(3)

4 parameters can define these spectrums: , , , which are respectively the amplitude of the log-amplitude spectrum, the slope, the cut off frequency of the logamplitude spectrum and the cut off frequency of the phase spectrum. Please note that amplitude of the phase spectrum is given by [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF]:

(4)
Then, by filtering a random Gaussian noise by such spectrum, scintillation time series are generated. An example of FFT (Fast Fourier Transform) processed logamplitude time series is presented Figure 1.

B. GNSS receiver model

The GNSS receiver model is based on a post-correlation approach reproducing a GPS C/A signal over 10MHz bandwidth integrated over 5ms. The implemented DLL / PLL / FLL are described in [START_REF] Kaplan | Understanding GPS, Principles and Applications[END_REF]. The delay and the phase are modeled considering the geometric range between satellite and receiver, and including the phase ionospheric time series generated by the scintillation model. All the other bias (e.g. tropospheric, clock …) are not modeled as they are out of the scope of the present study. Regarding the amplitude, the log-amplitude time series generated by the scintillation model is used. Finally, the GNSS receiver model is able to provide pseudo range, phase and Doppler estimations during both scintillation and non-scintillation events. As the signal is integrated other 5ms, the outputs rate is 200Hz, which allows a 100Hz as maximum span on the spectrum.

C. Scintillation estimator

The goal here is to estimate the log-amplitude and the phase scintillation spectrum based on the GNSS observables. To do so, the following procedure has been applied [START_REF] Van Dierendonck | Ionospheric Scintillation Monitoring Using Commercial Single Frequency C/A Code Receivers[END_REF]. For the phase processing, the first step is to remove cycle slips by checking the coherency between the phase and the Doppler estimate. If the difference between the estimated Doppler and the phase derivative estimated by the PLL loop is higher than a tolerance of 10Hz, a cycle slip is assumed, and this portion of the data is removed. The second step is to remove the satellite dynamic from the spectrum. To do so, a Butterworth high pass filter of 6 th order is used. The cut off frequencies of the filter are usually between 0.1Hz and 0.3 Hz. As a third step, the first second of the post-processed time series of observables is removed in order to avoid filter transition effects. Finally, the spectrum is estimated by the mean of a FFT. For the amplitude processing, the first step is to transform the post correlated I/Q samples into log-amplitude time series. Then, a "detrending" method is applied which consist in filtering the log-amplitude time series with a Butterworth low pass filter of 6 th order with, typical cut off frequencies between 0.1Hz and 0.3 Hz, in order to extract the signal "template". This template is subtracted to the original log-amplitude time series removing its slow variations. Finally, as for the phase processing, the first second of the post-processed time series of observables is removed. The spectrum is finally estimated by FFT. The reader is reminded that the scintillation indices ! and S4 are the results of the integration of the previously presented spectrums. Thus, spectrum estimation is fully equivalent to scintillation indices estimation.

"4 4 $ % (5) ! &$ ' % (6)

III. Simulation results

For each simulation, 300 seconds of GNSS signal have been simulated. For the first 180s, no scintillation occurs. During this time, we are able to estimate the frequency response of the receiver loops. Then, scintillation occurs between 180 and 300s. The estimation of the scintillation spectrums is done during this second time interval. Figure 2 illustrates this scintillation timing through the observation of the PLL error. The acquisition and phase recovery steps of the PLL can be observed during the first second. Then the PLL reaches it steady state before being perturbed by a phase scintillation event after t=180s. Table 1 gives typical scintillation parameters extracted from [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF] which were used for the next presented simulations. In Figure 3, an example of spectrum estimation for an unrealistic case of C/N0 = 200dBHz is presented. The idea is to observe the influence of the loops only on the scintillation spectrum estimation. Figure 3, spectrum named input is estimated directly from the output of the scintillation model, without passing through the GNSS receiver nor the scintillation estimator. It is considered as the reference spectrum and therefore the target of the estimation. The input (with filter) spectrum is estimated from the output of the scintillation model, which goes through the scintillation estimator only. This spectrum is therefore independent from the GNSS receiver behavior, and depends only of the scintillation estimator parametrization. The two other spectrums, no iono and iono, are based on the GNSS receiver outputs which go through the scintillation estimator. The no iono spectrum curve can be interpreted as the frequency response of the receiver, estimated when no scintillation occur (from 2s to 180s). The iono spectrum curve is estimated during a single scintillation event.

Finally, the asymptotic response of the PLL, and the theoretical pattern of the scintillation are added in dashed line.

In Figure 3, the impact of the PLL bandwidth can be highlighted. Here, a PLL bandwidth of 5Hz was used, which affects the slope of the phase spectrum scintillation. The slope of the iono spectrum suddenly changes around 5Hz, and no longer follows the input spectrum for frequencies higher than 5Hz. Another important observation is the influence of the high pass filter used in the scintillation estimator which cuts the low frequency part of the spectrum. Before the cut off frequency (fc = 0.2Hz in this case), all the spectrum going through the scintillation estimator, no iono, iono and input (with filter) curves, are cut.

Figure 4 presents spectrum estimations with a C/N0 = 45dbHz. In order to minimize the loop bandwidth influence, 10Hz bandwidth for the PLL was used. With such a C/N0, a double slopes effect can be seen, with a floor between these two slopes. The first slope is clearly the scintillation spectrum, as the input and iono spectrums perfectly match between the cut of frequency (fc = 0.1Hz) and the end of the first slope (f=0.3Hz). The 10 -4 floor observed on the spectrum is the noise floor. It can be clearly identified thanks to the no iono spectrum, which is flat over the PLL bandwidth as expected. From this frequency, the scintillation spectrum cannot be estimated anymore as it is below the noise. The last slope is the PLL response, which is not linked to the any atmospheric effect. The consequences on the scintillation indices are the following:

The phase slope should be estimated only in the interval where the effective scintillation can be measured. Thus, the first recommendation is to estimate the slope between the cut off frequency of the high pass filter "fc" and the PLL bandwidth. This recommendation leads to use an fc as small as possible (0.1 Hz) and a PLL bandwidth as big as possible (BPLL > 15 Hz). The reader is reminded that the frequency range is limited by the sampling rate of the output. Consequently, it is not required to use a PLL bandwidth higher than the sampling rate. In the case of low C/N0, the frequency of the noise floor (note F noise ) should be estimated. Then, the phase slope should be estimated between fc and F noise .

For the σ φ index, the cut off frequency of the high pass filter "fc" will have a big influence. Indeed, a major part of the energy of the phase spectrum is at low frequency as the cut off frequency of the phase spectrum are usually very small (~ 0.01 Hz) [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF]. In Figure 5 and Figure 6, simulated and estimated σ φ are compared considering a cut off frequency for the high pass filter respectively equal to 0.2 and 0.3Hz. The scintillation parameters are summarizes in table 1, and a second order PLL with a 10 Hz bandwidth is used. The C/N0 is equal to 50 dBHz. In these two figures, the following labels are used: the reference σ φ (labeled ref) is based on the output of the scintillation model, the estimated σ φ (labeled estimated) is based on the GNSS outputs, the filtered reference σ φ (labeled ref with filter) is based on the high pass filtered time series from the scintillation model. For the last, note that the GNSS receiver is bypassed, the spectrum estimator influence only is then evaluated. Comparing these two figures, the higher the fc, the larger the estimation error is. We also observe that the estimated and the ref with filter curves are very close, independent of fc. The conclusion is that the observed bias is not coming from the receiver tuning, but from the estimation strategy itself with cuts a major part of the scintillation energy. To better highlight the influence of the cut off frequency on the estimated σ φ index, raw data from a Septentrio GPS receiver placed in polar region (Tromso, Norway) were processed with different fc values. As we can see on Figure 7, σ φ values are strongly affected by the chosen cut off frequency. The PLL bandwidth influences the high part of the spectrum. However, this part of the spectrum has lower scintillation energy as the phase spectrum is decreasing. Thus, for the σ φ estimation, which is the integration of the phase spectrum, as stated in Eq.( 6), PLL bandwidth has almost no influence. The Figure 8 shows simulations results where different PLL bandwidths were tested. The scintillation parameters are presented in table 1, the C/N0 is equal to 50 dBHz and fc = 0.2 Hz. As visible, the σ φ estimations are close to the reference, with a difference of about 0.05rad, independent of the PLL bandwidth.

B. Log amplitude spectrum estimation

The log-amplitude spectrum is based on the I/Q time series, which are almost independent of the DLL and PLL parametrization. As an example, Figure 9 presents reference and estimated spectrums using the same naming convention as in Figure 3. Here, a 1Hz bandwidth for both DLL and PLL is used. The filtering effect of both loops is not predominant on the estimated spectrum. As a conclusion, the GNSS receiver tuning seems to have almost no impact on the log-amplitude spectrum estimation. Figure 10 and Figure 11, respectively compare the log-amplitude spectrum slope and the S4 index for both estimated and simulated scenarios with several PLL bandwidths. The ionospheric parameters are the ones in table 1, and we use fc = 0.1 Hz for the "detrending step". The C/N0 is equal to 50dBHz, and the slopes were estimated from fc to 50 Hz. As visible, the bandwidth does not affect the estimation of the log-amplitude parameters.

The "detrending step" in the log-amplitude spectrum estimation process is similar to a high pass filter as it removes the low frequency part of the spectrum. This similar behavior compared to phase spectrums can be observed Figure 9. However, the cut off frequency of the theoretical log-amplitude spectrum are much higher than the cut off frequency of the theoretical phase spectrum (~ 0.1Hz) [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF]. Thus, the fc parameters will less affect the S4 estimation, as less energy is cut. A theoretical comparison, between phase and log-amplitude spectrum, is propose in the annex in order to evaluate the power losses introduced by the high pass filter. It is demonstrated that for fc = 0.1Hz, a minimum of 10 dB on the total phase scintillation power is lost by the high pass filter. For the log-amplitude spectrum, a maximum of 6dB is lost. Regarding S4index estimation, only noise can affect the log-amplitude spectrum. It will create a floor, which may bias the slope estimation (same as the phase spectrum) and therefore the S4 estimation. Removing the noise influence in S4 estimation is a well-known procedure [START_REF] Van Dierendonck | Ionospheric Scintillation Monitoring Using Commercial Single Frequency C/A Code Receivers[END_REF]. However, contrary to the phase spectrum which is affected by loop effects at high frequency, the noise is not cut at high frequency in the log-amplitude spectrum. Thus, the frequency rate, or in other words the frequency range on which to estimate the spectrum, will affect the S4 estimation as the noise contribution will increase with the bandwidth. An example is given in Figure 13, where raw data from a Septentrio GPS receiver placed in polar region (Tromso, Norway) were processed with different sampling rates. The data were collected originally at 100 Hz, and several down sampling were performed to test the S4 estimator. It can be observed that the S4 index depends on the sampling rate. Regarding the σ φ index, the low frequency component is important as the theoretical phase spectrum shape is decreasing very fast. Thus, the high pass filter will remove a major part of the scintillation energy from a spectral point of view, which finally affects the σ φ value.

The users are advised that the σ φ index is strongly linked to the cut off frequency of the high pass filter. On the other hand, the GNSS receiver tuning will almost not affect the σ φ estimation. The log-amplitude spectrum is based on the I/Q samples, which are independent from the DLL and PLL parametrization. Finally, the GNSS receiver tuning will have almost no effect on the estimation of S4 and logamplitude slope. The "detreding step" of the scintillation estimator will cut the low frequency component of the log-amplitude spectrum inducing power losses without major consequences.

In the other hand, as the log-amplitude spectrum is not affected by loops, no filtering effects append at high frequency. Thus, the higher the data rate, the higher the integration frequency range in the spectrum which finally amplifies the noise contribution in the S4 estimation. As a conclusion, the users should be aware that S4 values will depend on the sampling rate.

Future work is still to be done in order to provide indices independent from the GNSS receiver tuning and scintillation estimator parametrisation. For example, satellite motion may be suppressed from the phase spectrum with the help of the geometrically computed Doppler. This solution may be substituted to the low pass filter applied during the detrending step, which cuts a large part of the phase spectrum.

Figure 1 :

 1 Figure 1: FFT of generated amplitude time series

Figure 2 :

 2 Figure 2: PLL error without (t=0:180s) and with scintillation (t=180:300)

Figure 3 :Figure 4 :

 34 Figure 3: example of phase spectrum estimation, C/N0 = 200dBHz

Figure 5 :Figure 6 :

 56 Figure 5: estimated vs simulated σ φ with fc = 0.1Hz

Figure 7 :

 7 Figure 7: Influence of fc based on raw data from Tromso GNSS receiver, 2013/06/01, PRN 20

Figure 8 :

 8 Figure 8: SigmaPhi estimation vs PLL bandwidth

Figure 9 :

 9 Figure 9: log-amplitude spectrums

Figure 12 :

 12 Figure 12: S4 estimation for different fc parameters

Figure 13 :

 13 Figure 13: S4 estimation for different sampling rate (denoted fe in the figure)

Figure 15 :

 15 Figure 15: Power lost due to high pass filter for logamplitude spectrum

Table 1 : Ionospheric scintillation parameters

 1

			(Hz)	(Hz)
	20	3.5	0.1	0.01
	A. Phase spectrum estimation

ACKNOWLEDGEMENT

The authors would like to thank NMA and NSC for providing data in polar region.

ANNEX

In this annex, we propose to evaluate to power cut by the high pass filter from a theoretical point of view. To that, we will compute the phase power between 0 and the cut off frequency

. Here, we use the equation ( 2), [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF].

Please note p > 3 [START_REF] Galiègue | Modélisation des effets des scintillations ionosphériques sur la propagation des ondes éléctromagnétiques en bande L aux latitudes polaires[END_REF]. Then let's note the equation as:

The total power is given by:

In order to evaluate the power cut by the high pass filter, we plot 10 * log + ? @A@ , ) * @A@ Bas a function of fc and for different slop value. As we can see, a minimum of 10 dB is cut due to the high pass filter (Figure 14) for fc = 0.1Hz. The same exercise has been done for the logamplitude spectrum. Here, a maximum of 6dB may be lost if we use fc = 0.1Hz.


Des ordinateurs aux robots : les machines en informatique 1

Pierre--Yves Oudeyer Qu'est--ce--qu'un ordinateur ? Un ordinateur est un machine physique capable de manipuler automatiquement et de mémoriser des informations (en particulier réaliser des suites d'opérations arithmétiques et logiques), et que l'on peut programmer. Bien que le mot « ordinateur » soit souvent associé à l'image de nos ordinateurs personnels, avec écran et clavier, les ordinateurs sont aujourd'hui présents dans de très nombreux objets quotidiens sous la forme de processeurs « embarqués » dans ces objets : téléphones, tablettes, montres électroniques, voitures, avions, feux de circulation, caméras, appareils photos, télévisions, appareils ménagers, jouets interactifs, appareils auditifs, thermostats, et tous les robots dont nous parlerons plus bas.

Des ordinateurs mécaniques aux microprocesseurs

Aujourd'hui la plupart des ordinateurs sont des machines électroniques (c'est à dire des machines dans lesquelles le traitement et la mémorisation des informations se fait uniquement par le biais de l'activation électrique de circuits). Cependant, les premiers ordinateurs ont d'abord été mécaniques. Charles Babbage conceptualisa en particulier le premier ordinateur au début du 19 ème siècle, sous la forme d'une machine composée d'un jeu de rouages mécaniques complexes et à laquelle il était possible de fournir des instructions de calculs générales sous la forme de cartes perforées (selon un principe similaire à la programmation des machines à tisser). Cette machine était en particulier capable de réaliser des boucles et des branchements conditionnels (c'est--à--dire de décider quelle suite d'instructions réaliser en fonction du résultat d'un calcul). Plus tard, pendant la seconde guerre mondiale, les premiers ordinateurs électro-mécaniques sont apparus, utilisés en particulier pour faire des calculs balistiques pour les missiles. Ceux--ci utilisaient l'électricité pour activer et désactiver des relais mécaniques afin de réaliser des calculs. Assez rapidement, des ordinateurs entièrement életroniques apparurent, basés sur l'utilisation de tubes à vide, comme le Colossus en Grande--Bretagne qui servit à décoder les messages secrets échangés par l'armée Allemande. Les tubes à vide furent ensuite remplacés dans les années 1950 par les transistors, dont la miniaturisation, la faible consommation d'énergie et la robustesse permirent leur accumulation sur ces circuits intégrés de très petite taille, ce qui donna naissance aux microprocesseurs à la fin des années 1950. Aujourd'hui des laboratoires de recherche imaginent des alternatives aux circuits électroniques pour faire réaliser à des machines des opérations de calculs et de manipulation d'informations. En particulier, ils étudient comment utiliser certaines propriétés de la matière comme les états quantiques afin d'accélérer certains types de

calculs (avec des « ordinateurs quantiques »), ou comment utiliser certaines molécules complexes comme l'ADN (en anglais on parle de « DNA computing »).

Les composants des ordinateurs

La plupart des ordinateurs existants aujourd'hui ont quatre grands types de composants : l'unité logique d'arithmétique (ALU), l'unité de contrôle, la mémoire, et des dispositifs d'entrées--sorties. Ces parties sont interconnectées avec des câbles que l'on appelle des « bus », et qui transmettent des informations suivant un codage que l'on appelle un « protocole de communication ». L'unité de contrôle est le chef d'orchestre qui contrôle l'ensemble des autres composants : elle lit et décode les instructions et les données de la mémoire, les transforme en signaux qui vont activer d'autres composants (par exemple envoyer à l'unité logique arithmétique un calcul a effectuer, récupérer le résultat et le ranger à un emplacement particulier dans la mémoire). On appelle souvent « CPU » (Central Processing Unit en anglais) l'ensemble composé de l'unité de contrôle, de l'unité arithmétique et logique, et de cellule mémoire spéciales appelées « registres » et utilisées pour mémoriser quel est l'emplacement en mémoire de la prochaine instruction. La mémoire d'un ordinateur peut être vue comme un ensemble de cases dans lesquelles ont peut stocker et écrire des nombres (encodés selon le codage binaire : un circuit élémentaire allumé code la valeur « 1 », et un circuit élémentaire éteint code la valeur « 0 »). Chaque case à une adresse qui permet aux autres composants de la localiser. Par exemple, des instructions peuvent demander à l'unité de contrôle de « mettre le nombre 345 dans la case numéro 256 » ou d' « additionner le nombre contenu dans la case 1356 et le nombre contenu dans la case 9595 ». Le programmeur qui conçoit les instructions données à l'ordinateur peut choisir d'utiliser ces nombres pour représenter et manipuler des informations qui ne sont pas des nombres : par exemple des lettres, des mots, des phrases, des formes, des couleurs, des sons, etc. Les ordinateurs sont souvent associés à des dispositifs d'entrée--sortie qui leurs permettent d'échanger des informations avec l'extérieur. Ces dispositifs sont appelés « périphériques ». Au delà des claviers, souris, disques durs, écrans, haut--parleurs ou imprimantes, les périphériques de la plupart des ordinateurs embarqués dans notre quotidien sont des capteurs (par exemple de lumière, de mouvement, de chaleur, de position GPS, …), des actionneurs (par exemple des moteurs, des valves, des radiateurs, …), et surtout d'autres ordinateurs qui peuvent être à l'intérieur du même objet (par exemple un ordinateur personnel ou un téléphone contiennent plusieurs ordinateurs dédiés au traitement du son ou de l'image et qui communiquent les uns avec les autres), ou à l'extérieur (les autres ordinateurs connectés au réseau internet).

Les robots

Parmi les objets qui embarquent des ordinateurs, les « robots » ont aujourd'hui une importance scientifique, sociétale et économique grandissante. Les robots sont partout : dans les usines et dans les champs, au fond des mers et dans l'espace, dans les jardins et les salons. Ils ont une importance économique grandissante, et d'aucuns prédisent qu'ils seront au XXIe siècle ce que la voiture fut au XXe siècle. En outre, ils n'ont pas seulement pénétré le monde industriel, ils ont aussi pénétré notre culture et certains d'entre eux participent au renouvellement de la vision que nous avons de nous--même. D'un point de vue technique, un robot est une machine dotée de capteurs qui lui permettent de percevoir son environnement (par exemple capteur de contacts, de distances, de couleurs, de force, …), de moteurs l'autorisant à bouger et à agir sur cet environnement, et d'un système électronique ou informatique qui contrôle, au moins en partie, ce qu'effectue le robot en fonction de ce qu'il perçoit. Une caractéristique fondamentale des robots est en particulier cette rétroaction entre la perception et l'action : ce que fait le robot ne dépend pas entièrement de ses instructions, mais aussi des propriétés physiques de son environnement perçues avec ses capteurs. Ainsi, selon cette définition, les automates (tels que ceux de Jacques de Vaucanson ou de Pierre et Henri--Louis Jaquet--Droz au XVIIIe siècle) ne sont pas des robots car leurs mouvements ne dépendent pas de ce qui se passe autour du robot (leurs propriétés et leurs enchainements sont entièrement pré--déterminés par le programme et il n'y a pas de capteurs). En pratique, cette définition recouvre une très vaste diversité de machines. Par exemple, les bras articulés et programmables dans les usines automobiles, les voitures (mode d'assistance à la conduite) et les avions qui sont aujourd'hui largement automatisés, les aspirateurs qui font le ménage tout seul, certains robots ludiques et objets électroniques des magasins de jouets, ou les robots bio--mimétiques à la forme de singes ou de poissons que l'on rencontre parfois dans les laboratoires de recherche. Cette diversité n'est pas seulement une diversité de formes, c'est aussi une diversité de fond. Les mécanismes de fonctionnement interne diffèrent beaucoup d'un robot à l'autre. Deux axes importants qui les distinguent sont l'autonomie et les capacités d'adaptation et d'apprentissage. L'autonomie : il existe des robots qui agissent sans besoin qu'un humain ne les guide et d'autres dont le comportement est, soit influencé par un humain soit, presque totalement contrôlé par un humain. Par exemple, dans une usine, les robots qui tra-vaillent à la chaîne et répètent toujours le même geste le font souvent de manière autonome. En revanche, les robots utilisés dans les centrales nucléaires pour opérer dans les zones à forte radiation, sont typiquement téléguidés par un humain qui leur dit où aller et quoi faire après chaque action ; L'adaptation et l'apprentissage : il y a des robots dont le comportement est figé au départ et une fois pour toute par le programmeur, alors que d'autres robots sont capables d'acquérir de nouveaux comportements et de nouvelles connaissances par leurs expériences : leur comportement évolue en fonction de l'histoire de leurs interactions avec l'environnement. Ainsi, certains robots sont capables d'apprendre à construire une carte de leur environnement et de la réutiliser plus tard afin de pouvoir y naviguer efficacement. D'autres sont capables d'apprendre à reconnaître des objets dans des images, des actions réalisées par des humains, ou à marcher tout seuls en expérimentant eux--mêmes différentes stratégies, les évaluant à chaque test et explorant des variations des meilleures. Ces mécanismes d'adaptation sont réalisés grâce à des « algorithmes d'apprentissage », qui reposent sur la détection automatique de régularités dans des flux de données captées par le robot, et sur des méthodes d' « optimisation » permettant de raffiner de manière progressive et itérative les paramètres des solutions à un problème. Grâce à ces algorithmes d'apprentissage, certains robots sont capables d'inventer des solutions et des comportements non prédits par leur concepteur, et même de sélectionner pour eux--mêmes des objectifs qui ne sont pas pré--programmés au départ. Il est par exemple possible de programmer un robot en lui donnant des instructions lui demandant de rechercher des situations nouvelles afin d'augmenter ses connaissances sur le monde qui l'entoure. Ainsi, certains algorithmes permettent de doter ces machines de formes d'apprentissage et de créativité. Néanmoins, les capacités et les performances de ces algorithmes sont aujourd'hui, et encore pour très longtemps probablement, très faibles en comparaison des capacités d'adaptation et de raisonnement de nombreux animaux, en en particulier de l'humain. Au delà des différences techniques entre robots, on trouve également des différences de finalité : il existe en effet une grande variété de raisons, et donc de fonctions, pour lesquelles les robots sont construits et utilisés. On peut en particulier considérer trois familles de fonctions : travailler et explorer, assister l'humain, modéliser les mécanismes cognitifs et comportementaux du vivant pour mieux les comprendre. Travailler et explorer. Parmi les robots qui sont en service dans le monde, la plupart sont industriels (il y a aujourd'hui environ 9 millions de robots industriels dans le monde). Très tôt, les entreprises se sont intéressées à ces machines pour plusieurs raisons : -- D'abord, les robots peuvent être utilisés pour remplacer les travailleurs humains dans les tâches répétitives, pénibles et nécessitant de faibles compétences, comme les montages, la peinture ou les soudures de pièces ; -- En outre, ces machines sont typiquement capables de réaliser des travaux à la chaîne bien plus rapidement et efficacement que les humains. Le premier robot de ce type est apparu en 1961. Il s'appelait Unimate et était installé dans une usine automobile de General Motors. C'était un bras articulé destiné à manipu-ler et à déplacer de lourdes pièces de fonderie. Dans les années 1970, l'usage des robots dans l'industrie a pris son envol. Aujourd'hui, les robots ont pénétré toutes les branches de l'industrie et ne sont plus restreints au seul secteur de l'automobile. Par exemple, dans l'agriculture et l'agro--alimentaire, les robots vont dans les champs cueillir fruits et légumes, certains coupent, pressent et mettent en bouteille ; d'autres trient et mettent en cartons ; d'autres encore groupent en palettes. Dans les aéroports, des flottes de robots s'occupent de transporter les bagages et de les charger dans les soutes. Les robots ne sont pas seulement utiles dans l'industrie pour les tâches simples et répétitives, ils sont aussi utilisés pour travailler dans des environnements dangereux pour l'homme. L'industrie nucléaire est un exemple typique. Qu'ils soient autonomes ou partiellement téléguidés, les robots des centrales nucléaires peuvent déambuler dans les enceintes confinées et radioactives, ils peuvent manipuler les substances dangereuses et s'occuper de la maintenance des autres machines. Un autre exemple est l'industrie pétrolière : les robots sous--marins sont par exemple utilisés pour contrôler l'état de la coque des navires afin de détecter les risques d'accidents et d'alerter sur les risques de navires poubelles. Les robots sont enfin cruciaux pour l'exploration des endroits où l'homme ne peut pas aller, au premier rang desquels on trouve l'espace et les planètes du Système solaire. C'est en 1966 que le premier robot mobile arrive sur la Lune, embarqué dans la sonde Surveyor. Suivront le robot soviétique Lunokhod, puis toute la série des robots Mariner (américains). En 1997, un robot atterrit sur la planète Mars : Sojourner est propulsé par l'énergie qu'il capte grâce à ses panneaux solaires ; il envoie à la Terre des milliers de clichés et provoque un engouement du grand public. Sojourner navigue en partie de manière autonome car, étant donné la distance avec la Terre, il est très difficile de le téléguider. En 2004, une nouvelle mission robotisée concentre l'attention du monde entier : Spirit et Opportunity, équipés de spectromètres et d'un bras qui leur permet de creuser sous la surface, font la preuve que des fleuves d'eau ont coulé sur Mars. Assister l'humain dans son quotidien. Si le XXe siècle a vu l'avènement des robots travailleurs et explorateurs, à l'aube du XXIe siècle, une autre grande famille de robots prend son essor : les robots d'assistance à la personne. Après avoir pénétré les usines, le fond des mers et les planètes lointaines, ils arrivent maintenant dans nos maisons, dans les magasins, sur nos lieux de travail, dans les hôpitaux. Au delà des robots assistants aux travaux ménagers, des robots apparaissent pour accompagner les personnes qui ont des difficultés physiques ou cognitives, par exemple pour les aider à se lever et s'asseoir, pour les stimuler cognitivement quand ils ont des problèmes de mémoire, ou pour jouer le rôle de facilitateur des contacts avec la famille ou l'entourage médical. Les robots d'assistance chirurgicale sont quotidiennement utilisés dans nombre d'hôpitaux depuis une quinzaine d'années. Plus récemment, d'autres robots, beaucoup plus petits, ont commencé a pousser la porte des blocs opératoires: des chercheurs mettent au point des capsules endoscopiques robotisées et miniaturisées capable d'explorer conduits intestinaux, artères ou veines pour aider le chirurgien dans son diagnostique. Apparaissent aussi depuis quelques années des prothèses de mains ou de bras entiers qui permettent à des personnes qui ont perdu le bras d'utiliser un bras robotique pour le remplacer.

Modéliser le vivant et la cognition

Les robots sont aussi utilisés par les chercheurs comme outils pour penser et modéliser le vivant. En particulier, ils commencent à être utilisés dans les laboratoires universitaires qui essaient de comprendre comment les êtres vivants s'adaptent et se comportent si efficacement dans leur environnement naturel. En effet, ce comportement est le résultat de la dynamique complexe des interactions entre un cerveau, un corps physique et un environnement. Et cette dynamique se redéfinit en permanence car le cerveau est modifié à chaque nouvelle interaction. Or, nous disposons de très peu d'outils scientifiques pour comprendre les dynamiques des systèmes complexes, dans lesquelles de nombreux composants s'influencent les uns les autres. Les systèmes complexes ne se retrouvent d'ailleurs pas seulement dans le fonctionnement des êtres vivants, ils sont partout dans la nature : météorologie, systèmes planétaires, formation des cristaux, apparition des dunes dans les déserts... De même que les ordinateurs sont au centre de l'étude et de la simulation des systèmes météorologiques, les robots sont ainsi devenus essentiels pour développer notre compréhension des phénomènes comportementaux et cognitifs. En effet, ils permettent d'implémenter des modèles qui prennent en compte au moins une partie de la complexité de l'interaction cerveau--corps--environnement, et il est possible de les mettre à profit pour faire des expérimentations à volonté, comme par exemple « éteindre » une partie de leur cerveau artificiel pour voir ce que cela modifie en terme de comportement. Dans le cadre de ces projets, les interactions avec les neurosciences, la biologie, la psychologie ou encore l'éthologie jouent un rôle central. Ainsi, certains laboratoires étudient par exemple les mécanismes de contrôle moteur et de perception visuelle qui permettent à aux singes ou aux humains de se mouvoir avec quatre ou deux jambes, les mécanismes grâce auxquels les animaux parviennent à se repérer dans l'espace, ou encore les mécanismes de l'apprentissage et e l'évolution de la parole et du langage chez l'humain.

Ce texte est sous licence Creative Commons CC--BY. Il a été également publié dans l'ouvrage « 1,

2, 3Codez » coordonné par la Fondation Main à la Pâte (http://www.fondation--lamap.org/node/34547), aux éditions Le Pommier.