▶ δ -(v ) -set of arcs in A incoming to v ∈ V Min a∈A c a x a S.t. a∈δ -(v ) x a = 1, v ∈ V , x a = p∈P h p a λ p , a ∈ A, p∈P λ p ≤ K , x a ∈ {0, 1}, a ∈ A, λ p ∈ {0, 1}, p ∈ P.
Route-based (extended) formulation

▶ δ -(v ) -set of arcs in A incoming to v ∈ V Min a∈A c a x a S.t. a∈δ -(v ) x a = 1, v ∈ V , x a = p∈P h p a λ p , a ∈ A, (π a ) p∈P λ p ≤ K , (µ) 0 ≤ x a ≤ 1, a ∈ A, 0 ≤ λ p ≤ 1, p ∈ P.
L ≤ 1λ ≤ U, x ∈ Z n + × R m + , λ ∈ Z k + .
A generic extended formulation ▶ x -a vector of original (natural) variables ▶ λ -a vector of extended variables (feasible sequences of decisions) ▶ H -mapping between natural and extended variables Min cx

S.t. Ax ≥ a, Bx ≥ b, ← robust 1 cuts x = Hλ, L ≤ 1λ ≤ U, x ∈ Z n + × R m + , λ ∈ Z k + .
A generic extended formulation ▶ x -a vector of original (natural) variables ▶ λ -a vector of extended variables (feasible sequences of decisions) ▶ H -mapping between natural and extended variables Min cx [START_REF] Belov | A cutting plane algorithm for the one-dimensional cutting stock problem with multiple stock lengths[END_REF]Scheithauer 2002, 2006) Gomory cuts for the 1D and 2D Cutting Stock Problems [START_REF] Jepsen | A Non-Robust Branch-And-Cut-And-Price Algorithm for the Vehicle Routing Problem with Time Windows[END_REF]Pisinger 2006, 2008) Subset-row cuts for Vehicle Routing Problems [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF] ▶ Replacing arc variables x in the set-partitioning constraints and relaxing to inequality:

S.t. Ax ≥ a, Bx ≥ b, ← robust cuts x = Hλ, Dλ ≥ d, ← non-robust 2 cuts L ≤ 1λ ≤ U, x ∈ Z n + × R m + , λ ∈ Z k + .

Contents

a∈δ -(v ) p∈P h p a λ p ≤ 1, v ∈ V . (1)
▶ Aggregating (1) for a set C ⊂ V , |C| = 3, with multiplier 1 2 :

p∈P 1 2 v ∈C a∈δ -(v ) h p a λ p ≤ 3 2 , (2)
▶ Performing Chvátal-Gomory rounding of (2):

p∈P     1 2 v ∈C a∈δ -(v ) h p a     λ p ≤ 1,
Subset-row cuts: example of violation

    v ∈C a∈δ -(v ) p v h p a     λ p ≤ v ∈C p v
All best possible multiplier vectors p for Chvátal-Gomory rounding of up to 5 constraints were found by (Pecin, Pessoa, Poggi, Uchoa, and Santos 2017):

Covering cuts

General idea

Take a subset of elements, and calculate the minimum number of "decision sequences" (columns) that should cover it.

Odd-circuit cuts (Nemhauser and Park 1991) Edge in an odd circuit in a graph should be covered by at least 3 matchings.

Strong k -path inequalities [START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF] A set C of customers should be visited by at least ⌈ i∈C d i /Q⌉ routes. Let V (p) be the set of customers visited by p, then

p∈P: V (p)∩C̸ =∅ λ p ≥ i∈C d i Q .
Strong k -path inequalities4 : illustration

∩ S = {i 1 , i 2 , i 4 , i 5 }, R 2 ∩ S = {i 2 , i 3 }, R 3 ∩ S = {i 1 , i 6 }, R 4 ∩ S = {i 5 , i 6 } and R 5 ∩ S = {i 4 , i 3 }.
Customer demands are reported in Fig. 1 next to each customer. This solution Capacitated Location-Routing Problem

▶ Depots k ∈ K with capacity W k ▶ Identical vehicles of capacity Q ▶ Clients i ∈ V with demand d i ▶ Matrix c of travelling costs Depot Client
Minimize the total travelling and depot opening costs ▶ such that every client is served ▶ each route starts and finished at a same opened depot ▶ depot and vehicle capacities are satisfied

Computational results for instances with all depots open

From [START_REF] Liguori | Non-Robust Strong Knapsack Cuts for Capacitated Location-Routing and Related Problems[END_REF] Extended formulation with linking variables ▶ x 0 -linking natural variables ▶ x 1 , x 2 -natural variables for blocks 1 and 2 ▶ λ 1 , λ 2 -extended variables representing all feasible solutions in blocks 1 and 2.

Min c 0 x 0 + c 1 x 1 + c 2 x 2 S.t. A 0 x 0 + A 1 x 1 + A 2 x 2 ≥ a, (x 0 , x 1 ) = H 1 λ 1 , (x 0 , x 2 ) = H 2 λ 2 , 1λ 1 = 1, 1λ 2 = 1, (x 0 , x 1 , x 2 ) ∈ Z n 0 +n 1 +n 2 , (λ 1 , λ 2 ) ∈ {0, 1} k 1 +k 2 .

Consistency cuts 5

▶ Let x0 be a particular partial solution for linking variables

x 0 (a pattern). ▶ Let P 1 (x 0 ) and P 1 (x 0 ) be the sets of solutions of block 1 and 2 which "match" pattern x0 .

Then the following constraint is valid

p∈P 1 (x 0 ) λ 1 p = p∈P 2 (x 0 ) λ 2 p .
Consistency cuts: computational results for the temporal knapsack instances Limited-memory Chvátal-Gomory rank-1 cuts (Pecin, Pessoa, Poggi, and Uchoa 2017) Optimal solution of instance M-n200-k16 D. Pecin et al.

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.

▶▶

  Vehicle Variable x a -arc a ∈ A is used in the solution or not ▶ Variable λ p -feasible route p ∈ P is used in the solution or not ▶ h p a = 1 if and only if path p contains arc a, otherwise 0

Fig. 1

 1 Fig. 1 Example of a set S that satisfies constraint (12) but violates constraint (13)

Figure 2 :▶

 2 Figure 2: Plot of time vs number of instances solved for the methods CP, BP, CPLEX, and CDG when all given a single thread and a 3-hour time limit. Method # of threads Time limit [h] #Optimal Avg. gap [%] Avg. time [s] CFM 1 1 1 3 8 0 . 0 7 N / A

Fig. 3

 3 Fig. 3 Optimal solution of M-n200-k16, value 1274 7.1 Handling symmetrical instances A high level of symmetry is observed in most instances from [14]. This impacts negatively the BCP performance, making the cutting and branching operations less 31 / 38

Marcus Poggi de[START_REF] Aragão | Integer program reformulation for robust branch-and-cut-and-price[END_REF]. "Integer program reformulation for robust branch-and-cut-and-price". In: Annals of Mathematical Programming in Rio. Ed. by Laurence A. Wolsey. Búzios, Brazil, pp. 56-61

[START_REF] Jepsen | A Non-Robust Branch-And-Cut-And-Price Algorithm for the Vehicle Routing Problem with Time Windows[END_REF]. A Non-Robust Branch-And-Cut-And-Price Algorithm for the Vehicle Routing Problem with Time Windows. Technical report 06/03. Dept. of Computer Science, University of Copenhagen

Mads Jepsen, Bjorn Petersen, Simon Spoorendonk, and David Pisinger (2008). "Subset-Row Inequalities Applied to the Vehicle-Routing Problem with Time Windows". In: Operations Research 56.2, pp.

497-511.

[START_REF] Baldacci | An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts[END_REF]. "An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts". In: Mathematical Programming 115, pp. 351-385.

Jens Vinther[START_REF] Clausen | An Exact Algorithm for a Rich Vehicle Routing Problem with Private Fleet and Common Carrier[END_REF]. "Consistency Cuts for Dantzig-Wolfe Reformulations". In: Operations Research Ahead of Print.

R. S. and François Vanderbeck (2021). BaPCod -a generic Branch-And-Price Code. Technical report HAL-03340548. Inria Bordeaux -Sud-Ouest.

[START_REF] Pessoa | A Generic Exact Solver for Vehicle Routing and Related Problems[END_REF]. "A Generic Exact Solver for Vehicle Routing and Related Problems". In: Mathematical Programming 183, pp. 483-523.

Arbitrary cuts of Chvátal-Gomory rank 1

Chvátal-Gomory rounding using a vector p of multipliers:

All best possible multiplier vectors p for Chvátal-Gomory rounding of up to 5 constraints were found by [START_REF] Pecin | Limited memory Rank-1 Cuts for Vehicle Routing Problems[END_REF]: Valid inequalities can be generated again by Chvátal-Gomory rounding... but not only.

The master knapsack polytope Theorem ( [START_REF] Aráoz | Cyclic group and knapsack facets[END_REF]))

The coefficient vectors ξ of the knapsack facets ξx ≤ 1 of polytope conv{t ∈ Z W + :

W q=1 qt q ≤ W }, are the extreme points of the following system of linear constraints

Road Load Knapsack Cuts (Liguori, Mahjoub, Marques, S., and Uchoa 2021)

Given a depot k ∈ K and a vector ξ ∈ R W k + satisfying (3) and ( 5), the inequality

is valid for the extended formulation of the CLRP.

Consistency cuts 5

▶ Let x0 be a particular partial solution for linking variables x 0 (a pattern). ▶ Let P 1 (x 0 ) and P 1 (x 0 ) be the sets of solutions of block 1 and 2 which "match" pattern x0 .

Then the following constraint is valid

When decomposition forms a staircase structure, there are no linking constraints, and linking variables are binary, then consistency cuts are enough to obtain an optimal integer solution (no branching is necessary).

Consistency cuts: application


Carolina Reyes
email: creyes@nt.tuwien.ac.at

Thibault Hilaire
email: thibault.hilaire@lip6.fr

Steffen Paul
email: steffen.paul@item.uni-bremen.de

Christoph F Mecklenbräuker

Evaluation of the Root Mean Square Error Performance of the PAST-Consensus Algorithm

  L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Subspace methods have become a very important subject of study in contemporary signal processing. There are several subspace estimation algorithms such as MUSIC, ESPRIT, and ESF which calculate for instance, the direction of arrivals (DOA) of plain waves impinging on a sensor array or estimate the frequencies of sinusoids that lie in a specific space [START_REF] Moonen | A singular value decomposition updating algorithm for subspace tracking[END_REF], [START_REF] Mitchley | Evaluation of selected subspace tracking algorithms for direction finding[END_REF], [START_REF] Kumaresan | Estimating the angles of arrival of multiple plane waves[END_REF], [START_REF] Rabideau | Fast rank adaptive subspace tracking and applications[END_REF], [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF].

However, these traditional methods present a high complexity due to the fact that they rely on eigenvalue/singular value decomposition for estimating the signal or the noise subspace. These algorithms have to calculate the sample covariance matrix every time a new sample arrives, which leads to a high amount of internal processing.

The goal of this paper is to evaluate the performance of the PAST-Consensus algorithm developed in [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF]. This method is a distributed version of the Projection Approximation Subspace Tracking (PAST) [START_REF] Yang | Projection approximation subspace tracking[END_REF], a well-known algorithm whose major advantage is the considerably low complexity. PAST introduces a new signal subspace model interpretation and presents a solution where the signal subspace is calculated in a recursive manner at time t, depending on the previous subspace estimate at time t -1 and the new sample x(t).

Likewise, Consensus Propagation is a remarkable protocol widely used for obtaining averages over a network [START_REF] Moallemi | Consensus propagation[END_REF], [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF], [START_REF] Olfati-Saber | Consensus protocols for networks of dynamic agents[END_REF]. The work in [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF] has introduced a distributed algorithm based on Consensus Propagation for frequency estimation over a wireless sensor network.

Notation: We denominate a column vectors in underlined boldface and matrices in uppercase and boldface. The superscript H represents complex conjugation, ⊤ transposition, tr(. ) the trace operator and . the Euclidean vector norm.

Organization of the paper:

Section II presents an overview of the analytical model and constraints. The section III introduces the PAST-Consensus algorithm as implemented in [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF]. Finally, in Section IV we will present some simulation results that evaluate the root mean square error for several signal to noise ratios, as well as the signal subspace error expressed in terms of principal angle components.

II. AN OVERVIEW OF THE PROJECTION APPROXIMATION

SUBSPACE TRACKING ALGORITHM Let x(t) ∈ C N be the data vector observed at time t, composed of r narrow-band signal waves impinging on a planar array of N sensors hidden in additive noise (zero-mean with variance σ 2 ). The data vector is modelled as

x(t) = r i=1 s i (t)a(w i ) + n(t) = A(ω)s(t) + n(t) (1)
with,

A =        1 1 1 e jω1
e jω2 e jωr . . . . . . . . . e (N -1)jω1 e (N -1)jω2 e (N -1)jωr

       , s(t) =     s 1 (t) . . . s r (t)     and n(t) =     n 1 (t) . . . n r (t)     .
Here, A(ω) is a deterministic N × r signal mixing matrix depending on ω = (ω 1 , . . . , ω r ). Its columns are plane wave steering vectors a(ω i ) and the n-th element of a(ω) is defined by

[a(ω i )] n = exp j 2π λ (ξ n cos θ i + η n sin θ i ) / √ N
where ω i = cos θ i are the frequency that lies in the subspace of A and that we intend to track. The Cartesian coordinates of the n-th sensor node are (ξ n , η n ) and λ is the wavelength, s(t) is the vector of complex signal amplitudes and n(t) is uncorrelated additive noise. We allow ω to be slowly time-varying and seek an estimate for matrix W (t) whose columns span the same subspace as A(ω(t)). A tracking algorithm estimates W (t) by a function of the previously estimated matrix W (t -1) and the new observation x(t) alone. Particularly, Yang [START_REF] Yang | Projection approximation subspace tracking[END_REF] approximates the cost function to minimize

J(W (t)) = t i=1 β t-i x(i) -W (t)W H (t)x(i) 2 (2)
by

J ′ (W (t)) = t i=1 β t-i x(i) -W (t)y(i) 2 (3)
with

y(i) = W H (i -1)x(i) . (4)
Here, the vector y(t) stores the modified data vector resulting from the new incoming data x(t) and the previous estimated signal subspace W (t). Besides, the matrix W ∈ C N ×r is constrained to rank r < N and the minimization of cost function ( 3) results in a low-complexity update of the signal subspace.

III. PAST-CONSENSUS

A. Consensus Propagation

In the following, consider a wireless sensor network composed of N = 36 nodes placed in an irregular cartesian grid, as shown for instance in Figure 1. The average distance between nodes is λ/2 and the transmission range is set to 1.44λ/2. For every node n = 1, . . . , N in the network, N n denominates the set of the adjacent nodes including itself. In this setting each node in the neighborhood N n can correctly receive broadcasted messages from nodes within its neighborhood with probability 1. One can observe on Algorithm 1 that at the end of step t -1, every node n sends its own estimation of the average y n (t -1) and a weight w n to its adjacent nodes N n . Then, at the beginning of step t, every node n receives the pairs {y i (t -1), w i } i∈Nn from its neighbors and compute a new average, to be sent at the end of step t.

The aforementioned weights are constant for every node and defined as

w n = 1/ |N n | ∀1 ≤ n ≤ N. (5)
Finally, the local observation vector x n (t) is obtained with the aggregation of data observed by the adjacent nodes, denoted by {x j (t -1)} j∈Nn .

B. PAST-Consensus propagation algorithm

The distributed subspace tacking algorithm presented in [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF], has been conceived for sensor network applications. Due to its low computational complexity O(nr), scalability and robustness, it seems a suitable approach for distributed sensor network applications.

Algorithm 1: Consensus Propagation Algorithm

for n := 1, 2, . . . , N do Input: {y j (t -1), w j } j∈Nn are the pairs sent to node n in step t -1

y n (t) = j∈Nn y j (t -1)w j / j∈Nn w j (6)
Broadcast the pair {y n (t), w n } to all nodes in N n Output: y n (t) is the estimation of the average in step t at node n endfor Sensor systems are well known for a wide range of applications such as the monitoring of environmental phenomena, i.e. avalanches, floods. Constant improvement is sought in the areas of power consumption, internal processing load, scalability and robustness.

The PAST-Consensus algorithm is an alternative for tackling some of this problems. In addition, it locally averages the vector y n (t) in n with information from N n . Afterwards, n broadcasts its local observation x n (t), the locally filtered r-dimensional vector y n (t), and a weight w n . To conclude, y n (t) in equation ( 4) contains information from the updated signal subspace at t-1 as well as new observation data x n (t). The aforementioned algorithm is given by Algorithm 2.

In the following, we will show some simulation results which analyze the performance of the algorithm in terms of the root mean square error and calculate the error between the original and the estimated subspace at node n. Here, the distance between two subspaces is equal to two times the sum of the squared sines of the principal angles [START_REF] Edelman | The geometry of algorithms with orthogonality constraints[END_REF].

IV. SIMULATIONS

In order to make a consistent data performance evaluation, we propose the following simulation scenario. Consider a planar sensor network with N = 36 nodes, whose position resembles an imperfect Cartesian grid. The exact sensor array geometry is documented in [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF] and represented in Figure 1.

Algorithm 2: PAST-Consensus [START_REF] Reyes | Distributed projection approximation subspace tracking based on consensus propagation[END_REF] Input: β, P 1 (0), . . . , P N (0), W 1 (0), . . . , W N (0) for t := 1, 2, . . . do for n := 1, 2, . . . , N do Input: x(n) aggregate x n (t) = S n x(t -1) from all nodes∈ N n y n (t) = W H n (t -1)x n (t) apply Algorithm 1 for locally averaging y n (t)

h n (t) = P n (t -1)y n (t) g n (t) = h n (t)/[β + y H n (t)h n (t)] P n (t) = 1 β [P n (t -1) -g n (t)h H n (t)] e n (t) = D n (x(t) -W n (t -1)y n (t)) W n (t) = W n (t -1) + e n (t)g H n (t) broadcast {x n (t), y n (t), w n } to all nodes∈ N n endfor endfor
Here, every sensor network observes 1000 snapshots in time (t = 1, ..., 1000) and tracks the instantaneous frequency ω 1 (t) = 0.1 of a single impinging wave which is constant over time. The forgetting factor is chosen to be β = 0.97, W (0) and P (0) are initialized as identity matrices. We assume that messages sent by node n are received correctly by all nodes in its neighborhood N n . The spectral MUSIC algorithm is used for extracting the r = 1 frequency component ω1 (t) from the signal subspace matrix W (t).

In the first experiment, we consider the performance at node 17 and compare the behavoir of the PAST-Consensus algorithm (see Figure 2, curve b)) with: a) the original PAST-algorithm in [START_REF] Yang | Projection approximation subspace tracking[END_REF], where only node 17 and its neighborhood N 17 are considered (N = 6). c) the original PAST algorithm applied to the global sensor network. Here, we suppose that each node can simultaneously access all nodes in the network, and we keep using the non-distributed version of [START_REF] Yang | Projection approximation subspace tracking[END_REF]. Again, the starting point for such calculations is the node 17. d) the average RMSE when considering all nodes in the network N = 36.

In our simulations, the root mean square error formula at each node for one single frequency (r = 1) is defined as:

RMSE n = 1 901 1000 t=100 |ω 1,n (t) -ω1,n (t)| 2 , ( 7
)
where ω 1,n (t) is the true frequency value and ω1,n (t) its estimate. Since the Algorithm 2 does not require any a priori data statistics, the initial information tends to be the least reliable. Due to this fact, we neglect the first one hundred snapshots in our graphs in order to obtain a more representative error sampling.

For the different simulations, we consider the additive noise n(t) (see eq. ( 1)) with a Signal Noise Ration (SNR) from Note that this definition of RMSE n is the root of the time averaged square error, which is different from the classical definition based on the ensemble average of the square error.

In Figure 2, we observe the performance at node 17 for different signal to noise ratios. In a) the original PAST algorithm is analyzed for the neighborhood N 17 (N = 6, r = 1). Likewise, in b) we examine the root mean square error (RMSE) for the same neighborhood, but employing the PAST-Consensus based on Consensus Propagation. Both errors a) and b) behave the same way at the beginning, since the initialization parameters for the algorithms were the same. However, our distributed approach seems to outperform the centralized solution after -19dB, where it stays constant. Note that for SNR values higher than 10 dB the error factor decays faster. As expected in c), the error evaluation for n=17 considering the whole neighborhood N = 36 has the best performance, since all nodes have access to all data available in the sensor network.

In the second experiment, we consider again 1000 iterations for the case of two crossing frequencies, ω 1 =-0.5 to 0.5 and ω 2 =0.5 to -0.5. The simulation runs are carried out with the same noise realization. However, the SNR for the image a) is set to 5dB and for b) the SNR= -10dB. We observe keep simulating for the node 17.

Figure 3 shows the principal angles between the subspace spanned by the columns of the estimated signal subspace W n (t) and the original signal subspace spanned by the columns of the matrix A, as in [START_REF] Moonen | A singular value decomposition updating algorithm for subspace tracking[END_REF]. Notice that the graph a) shows a good tracking capability for a signal to noise ratio set to 5dB. However, in the graph b) the estimation accuracy diminishes for a SNR=-10dB. The principal angles are zero if the subspaces compared are the same. However, this is not case for the graphs c) and d), where the angles never reach zero. Nevertheless, one can observe that the evolution andd): Distance between two subspaces expressed in terms of its principal angles. Namely, the signal subspace matrix Wn(t) and the signal mixing matrix A of N 17 (see eq. ( 1)). Again, the calculations are for n = 17, SN R = 5dB and SN R = -5dB correspondingly of the error in the subspace seems not to be dramatically affected for different SNR values. The reason for this is that the PAST-Consensus algorithm, as previously mentioned, estimates the signal subspace in a recursive manner and does not need any knowledge of the eigenvalues. The algorithms based on eigenvalue decomposition are very sensitive to low SNR values, because the gap separating the signal from the noise subspace is reduced and is more difficult to make a correct identification of both.

V. SUMMARY AND CONCLUSION

We have calculated the distance between the estimated signal subspace W (t) and the original subspace introduce by the matrix A, as in (1) in terms of the principal angles. The results are quite interesting: On one hand the error between subspaces does not lead to zero, it rather stays constant over time for values between 0.03 and 0.01 radians. On the other hand, one can observe that the signal subspace interpretation in [START_REF] Yang | Projection approximation subspace tracking[END_REF] leads to an algorithm more robust against errors, when considering too closely spaced or very weak signals.

The RMSE simulation result indicates that the PAST-Consensus algorithm has a performance in terms of the RMSE between the (centralized) PAST algorithm [START_REF] Yang | Projection approximation subspace tracking[END_REF] which uses the data of all sensors and a PAST algorithm which uses only the sensor data of the immediately neighboring nodes. In future research, we will optimize the weighting coefficients used for consensus propagation which will further reduce the RMSE of the PAST-Consensus algorithm while conserving its good tracking capability and low complexity.

Fig. 1 .

 1 Fig. 1. Sensor network with neighborhood N 17 for radius 9.

Fig. 2 .

 2 Fig. 2. Root Mean Square Error comparison for r = 1: a) PAST result for N 17 , N = 6, b) PAST-Consensus result for neighborhood N 17 , N = 6, c) PAST result for whole sensor array N = 36

Fig. 3 .

 3 Fig. 3. Figures a) and b): Two crossing frequencies (r=2, n=17) showing the tracking capabilities of the PAST-Consensus algorithm for SNR=5dB and SNR=-5dB respectively. Figures c)and d): Distance between two subspaces expressed in terms of its principal angles. Namely, the signal subspace matrix Wn(t) and the signal mixing matrix A of N 17 (see eq. (1)). Again, the calculations are for n = 17, SN R = 5dB and SN R = -5dB correspondingly

ACKNOWLEDGMENT Support was provided by the fFORTE-Women in Technology Program at Vienna Univ. of Technology, the Austrian Ministry for Science and Research, the Christian Doppler Lab Wireless Technologies for Sustainable Mobility, and the WWTF project Distributed Information Processing for Spatio-Temporal Fields in Wireless Sensor Networks. We thank the national research network on Signal and Information Processing in Science and Engineering (NFN-SISE) for many fruitful discussions.


Thibault Hilaire

On the Transfer Function Error of State-Space Filters in Fixed-Point Context

Keywords: Digital filter implementation, coefficient sensitivity, fixed-point implementation

This paper presents a new measure used for the implementation of filters/controllers in state-space form. It investigates the transfer function deviation generated by the coefficient quantization. The classical L2-sensitivity measure is extended with precise consideration on their fixed-point representation, in order to make a more valid measure. By solving the related optimal realization problem, fixed-point accurate realizations in state-space form can be found.

I. INTRODUCTION

The majority of control or signal processing systems is implemented in digital general purpose processors, DSPs 1 , FPGAs 2 , etc. Since these devices cannot compute with infinite precision and approximate real-number parameters with a finite binary representation, the numerical implementation of controllers (filters) leads to deterioration in characteristics and performance. This has two separate origins, corresponding to the quantization of the embedded coefficients and the roundoff errors occurring during the computations. They can be formalized as parametric errors and numerical noises, respectively. The focus of this paper are parametric errors, but one can refer to [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF]- [START_REF] Mullis | Synthesis of minimum roundoff noise fixed point digital filters[END_REF] for roundoff noises, where measures with fixed-point consideration already exist. It is also well known that these Finite Word Length (FWL) effects depend on the structure of the realization. In state-space form, the realization depends on the choice of the basis of the state-vector. This motivates us to investigate the coefficient sensitivity minimization problem. It has been well studied with the L 2 -measure [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF], [START_REF] Hinamoto | Analysis and minimization of L 2 -sensitivity for linear systems and two-dimensional state-space filters using general controllability and observability gramians[END_REF]. However this measure only considers how sensitive to the coefficients the transfer function is, and does not investigate the coefficients' quantization, which depends on the fixed-point representation used. In [START_REF] Hinamoto | Analysis and minimization of L 2 -sensitivity for linear systems and two-dimensional state-space filters using general controllability and observability gramians[END_REF], the transfer function error is exhibited for the first time, however, only for quantized coefficients with the same binary-point position.

This paper investigates the transfer function deviation generated by the coefficient quantization with precise consideration on their fixed-point representation. The classical L 2 -sensitivity analysis is shown in section II, whereas the new approach, based on fixed-point consideration, is presented in section III. A comparison with the L 2 -sensitivity and some scaling considerations are provided. Finally, the optimal realization problem is solved in section IV and a numerical example is exhibited before conclusion.

II. L 2 -SENSITIVITY ANALYSIS

Let (A, b, c,d) be a stable, controllable and observable linear discrete time Single Input Single Output (SISO) statespace system, i.e.

x(k

+ 1) = Ax(k)+bu(k) y(k)=cx(k)+du(k) (1)
where

A ∈ R n×n , b ∈ R n×1 , c ∈ R 1×n and d ∈ R. u(k) is the scalar input, y(k) is the scalar output and x(k) ∈ R n×1
is the state vector. Its input-output relationship is given by the scalar transfer function h : C → C defined by:

h : z → c(zI n -A) -1 b + d. ( 2
)
The quantization of the coefficients introduces some uncertainties to A, b, c and d leading to A +∆A, b +∆b, c +∆c and d +∆d respectively. It is common to consider the sensitivity of the transfer function with respect to the coefficients, based on the following definitions.

Definition 1 (Transfer function sensitivity) Consider X ∈ R m×n and f : R m×n → C differentiable with respect to all the entries of X. The sensitivity of f with respect to X is defined by the matrix

S X ∈ R m×n : ∂f ∂X S X with (S X ) i,j ∂f ∂X i,j . (3)
Applied to a scalar transfer function h where h(z) depends on a given matrix X, ∂h ∂X is a transfer function of a Multiple Inputs Multiple Outputs (MIMO) system. Definition 2 (L 2 -Norm) Let H : C → C k×l be a function of the scalar complex variable z. H 2 is the L 2 -norm of H, defined by:

H 2 1 2π 2π 0 H (e jω ) 2 F dω (4)
where . F is the Froebenius norm.

Proposition 1 If H is the MIMO state-space system (K, L, M , N ), then its L 2 -norm can be computed by

H 2 2 = tr(NN + MW c M ) (5) = tr(N N + L W o L) (6)
where W c and W o are the controllability and observability Gramians, respectively. They are solutions of the Lyapunov equations

W c = KW c K + LL , W o = K W o K + M M . ( 7
)
Proof: See [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF]. Gevers and Li [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF] have proposed the L 2 -sensitivity measure to evaluate the coefficient roundoff errors. It is defined by

M L2 ∂h ∂A 2 2 + ∂h ∂b 2 2 + ∂h ∂c 2 2 + ∂h ∂d 2 2 (8)
and can be computed by ∂h ∂A (z)=G (z)F (z), ∂h ∂b (z)= G (z), ∂h ∂c (z)=F (z) and ∂h ∂d (z)=1, with

F (z) (zI n -A) -1 b, G(z) c(zI n -A) -1 . ( 9
)
This measure is an extension of the more tractable but less natural L 1 /L 2 sensitivity measure proposed by V. Tavs ¸anoglu and L. Thiele [START_REF] Tavs | Optimal design of state-space digital filters by simultaneous minimization of sensibility and roundoff noise[END_REF] (

∂h ∂A 2 1 instead of ∂h ∂A 2 2 in (8)).
Remark 1 To simplify the expressions, it is also possible to regroup all the coefficients in one unique matrix Z:

Z Ab c d . (10)
Then, with L 2 -norm property, M L2 = ∂h ∂Z 2 2 . From (9) and the associated state-spaces, the sensitivity transfer function ∂h ∂Z can be described by the MIMO state-space system ( Ã, B, C, D) with

à Ab c 0 A , B 0 b I n 0 , C I n 0 0 c , D 00 0 1 . (11)
The proposition 1 is used to compute M L2 . See [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF] and [START_REF] Hilaire | On the compact formulation of the derivation of a transfer matrix with respect to another matrix[END_REF] for more details.

Applying a coordinate transformation, defined by x(k)

U -1 x(k) to the state-space system (A, b, c,d), leads to a new equivalent realization (U -1 AU , U -1 b, cU ,d).
Since these two realizations are equivalent in infinite precision but are no more equivalent in finite precision (fixed point arithmetic, floating-point arithmetic, etc.), the L 2 -sensitivity then depends on U , and is denoted M L2 (U ). In this case, it is natural to define the following problem: Problem 1 (optimal L 2 -sensitivity problem) Considering a state-space realization (A, b, c,d), the optimal L 2 -sensitivity problem consists of finding the coordinate transformation U opt that minimizes M L2 :

U opt = arg min U invertible M L2 (U ). (12)
In [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Probems[END_REF], it is shown that the problem has one unique solution.

Hence, for example, a gradient method can be used to solve it.

III. TRANSFER FUNCTION ERROR

A. Fixed-point implementation

In this paper, the notation (β, γ) is used for the fixed-point representation of a variable or coefficient (2's complement scheme), according to Figure 1. β is the total wordlength of the representation in bits, whereas γ is the wordlength of the fractional part (it determines the position of the binary-point). They are fixed for each variable (input, states, output) and each coefficient, and implicit (unlike the floating-point representation). β and γ will be suffixed by the variable/coefficient they refer to. These parameters could be scalars, vectors or matrices, according to the variables they refer to.

± 2 1 2 0 2 -1 ... ... 2 β-γ-2 β -γ -1 β γ 2 -γ integer part fractional part s Fig. 1. Fixed-point representation
Let us suppose that the wordlength of the coefficients β Z is given 3 . Then, the coefficients Z ij are represented in fixedpoint by (β Zij ,γ Zij ) with:

γ Zij = β Zij -2 - log 2 |Z ij | . ( 13
)
where the a operation rounds a to the nearest integer less or equal to a (for positive numbers a is the integer part).

Remark 2

The binary point position is not defined for null coefficients, however this is no problem because these coefficients will not be represented in the final algorithm (the null multiplications are removed). So, in order to consider coefficients that will be quantized without error, we introduced a weighting matrix δ Z such that

(δ Z ) ij 0 if Z ij is exactly implemented 1 otherwise. ( 14
)
The exactly implemented coefficients are 0, ±1 and also positive and negative coefficients of power of 2.

Remark 3 In some specific computational cases the fixedpoint representation chosen for the coefficients is not always the best one as defined in [START_REF] Parhi | VLSI Digital Signal Processing Systems: Design and Implementation of Digital Controllers[END_REF]. For example, in the Roundoff Before Multiplication scheme, some extra quantizations are added to the coefficients, in order to avoid shift operations after multiplications [START_REF] Hilaire | Bit accurate roundoff noise analysis of fixed-point linear controllers[END_REF]. Only the classical case (corresponding to the Roundoff After Multiplication) is considered here, as defined by equation [START_REF] Parhi | VLSI Digital Signal Processing Systems: Design and Implementation of Digital Controllers[END_REF].

Remark 4 It is also possible to choose the same fixedpoint representation for all the coefficients (determined by the coefficient with the highest magnitude). But in that case, the lowest coefficients (in magnitude) do not have an appropriate representation. They are coded with less meaningful bits and have a higher relative error. When the ratio between the greatest and lowest magnitude is too high, then underflows occurs for the lowest coefficients that cannot be represented. For example, this is common for the Direct Form realizations with high (or low) L 2 -gain.

During the quantization process, the coefficients are changed from Z into Z † Z + ∆Z. For a best-roundoff quantization, the {∆Z i,j } are independent centered random variables uniformly distributed [START_REF] Sripad | A necessary and sufficient condition for quantization error to be uniform and white[END_REF] within the ranges

-2 -γ Z ij -1 Z i,j < 2 -γ Z ij -1
, so their second-order moments are given by

σ 2 Zij E (∆Z ij ) 2 (15) = 2 -2γ Z i j 12 δ Zij , (16)
where E{.} is the mean operator.

B. Transfer function error

Due to the quantization of the coefficients, the transfer function is changed from h to h † h +∆h. This degradation can be evaluated in a statistical way with the following definition.

Definition 3 (Transfer function error) A measure of the transfer function error can be statistically defined by [START_REF] Hinamoto | Analysis and minimization of L 2 -sensitivity for linear systems and two-dimensional state-space filters using general controllability and observability gramians[END_REF] σ 2 ∆h

1 2π 2π 0 E ∆h e jω 2 dω. (17)
The transfer function error is a tractable measure that can be evaluated with the following propostion.

Proposition 2

The transfer function error is given by:

σ 2 ∆h = ∂h ∂Z × Ξ Z 2 2 ( 18
)
where × is the Schur product, Ξ Z ∈ R (n+1)×(n+1) defined by:

(Ξ Z ) ij 2 -β Z ij +1 √ 3 Z ij 2 (δ Z ) ij if Z ij =0 0 if Z ij =0 (19)
and x 2 is nearest power of 2 lower than |x|:

x 2 2 log 2 |x| . ( 20
)
Proof: A first order approximation gives

∆h(z)= i,j ∂h ∂Z ij (z)∆Z ij , ∀z ∈ C. (21)
Hence

E ∆h e jω 2 = i,j ∂h(e jω ) ∂Z ij 2 σ 2 Zij ( 22
)
because the random variables ∆Z ij are independent. Considering ( 13) and ( 16) for non-null coefficients, we get

σ 2 Zij = 4 3 2 -2β Z ij Z ij 2 2 (δ Z ) ij (23)
and

σ 2 ∆h = i,j (Ξ) ij ∂h ∂Z ij 2 2 (24)
Then, with ∂h ∂Z ij = ∂h ∂Zij and (4), eq. ( 18) holds.

Remark 5 In the classical case where the wordlength of the coefficients are all the same (equal to β), we can define a normalized transfer error σ2 ∆h defined by:

σ2 ∆h 3σ 2 ∆h 2 -2β+2 . ( 25
)
This measure is now independent of the wordlength, and can be used for some comparisons.

C. Comparison with the classical M L2 measure

It is of interest to remark the relationship with the classical M L2 measure. In [START_REF] Hinamoto | Analysis and minimization of L 2 -sensitivity for linear systems and two-dimensional state-space filters using general controllability and observability gramians[END_REF] where the transfer function error appears for the first time, the coefficients are supposed to have the same fixed-point representation, so their second-order moment (σ 2 Zij ) are all equal and denoted σ 2 0 . The M L2 satisfies then

M L2 = σ 2 ∆h σ 2 0 . (26)
Here, the transfer function error σ 2 ∆h can be seen as an extension of the M L2 measure with fixed-point considerations. The sensitivity is weighted according to the variance of the quantization noise of each coefficient.

The M L2 measure considers each coefficient with the same weight, even if the quantization of one particular coefficient induces a small or big modification of this coefficient. To avoid this problem, a normalization has been created. Since a coordinate transformation U = αI n does not change A, but only multiplies the coefficients of b by 1 α , and those of c by α, it was of interest to set a condition on b or c for normalization.

The L 2 -dynamic-range-scaling constraints have been introduced by Jackson in [START_REF] Jackson | Roundoff-noise analysis for fixed-point digital filters realized in cascade or parallel form[END_REF] and Hwang in [START_REF] Hwang | Dynamic range constraint in state-space digital filtering[END_REF]. It consists in scaling the state-variable vector so as to prevent overflows or underflows during its evaluation and it imposes a condition on b (it avoids big values for c and small for b, and vice versa). (A,b,c,d) is said to be L 2 -scaled if the transfer functions from the input to each state have unitary L 2 -norms, i.e.:

Definition 4 (L 2 -scaling) A state-space realization

e i (zI n -A) -1 b) 2 =1, ∀1 i n (27)
where e i is the column vector of appropriate dimension and with all elements being 0 except for the i th element which is 1.

With proposition 1 applied to the system (A, b, e i , 0), the L 2 -scaling constraints (27) can be expressed as:

(W c ) i,i =1, ∀1 i n (28)
where W c is the controllability Gramian of the state-space system (A, b, c,d).

In addition, it has been shown in [START_REF] Hilaire | Low parametric sensitivity realizations with relaxed l 2dynamic-range-scaling constraints[END_REF] that the L 2 -scaling is not necessary to implement a realization without overflows. Two equivalent choices are possible for the implementation:

• define a binary-point position for each state (for example, the same as the input), and apply a scaling to them in order to adapt the peak values of each state to the chosen binary-point position. This scaling can also be based on a relaxed L 2 -scaling, as in [START_REF] Hilaire | Low parametric sensitivity realizations with relaxed l 2dynamic-range-scaling constraints[END_REF]. • or set the binary-point position for each state, according to (29), to make sure that the fixed-point representation of the states avoids state-overflows:

γ xi = β xi -2 - log 2 up x i . ( 29
)
where up

x i is an upper-bound of the i th state, given by a L 1 -norm:

up x i = e i (zI n -A) -1 b 1 max u, (30)
or estimated by a L 2 -norm [START_REF] Hilaire | Low parametric sensitivity realizations with relaxed l 2dynamic-range-scaling constraints[END_REF], [START_REF] Parhi | VLSI Digital Signal Processing Systems: Design and Implementation of Digital Controllers[END_REF]:

up x i δ e i (zI n -A) -1 b 2 max u. (31)
The L 2 -scaling constraints have to be applied in the first case, whereas no other constraint has to be applied in the second case.. Contrary to the L 2 -sensitivity measure M L2 , the transfer function error σ 2 ∆h can be applied in a general cases and be meaningful. In the very particular case where all the coefficients have the same fixed-point representation, these two measure are equivalent. But for other cases, only the σ 2 ∆h measure is a meaningful measure of the degradation of the transfer function due to the quantization of the coefficients.

D. Scaling

Let us consider a scaling of the states: x(k) is changed in U -1 x(k) with U an invertible diagonal matrix. The realization Z 0 is changed into Z 1 :

Z 1 = T -1 Z 0 T , with T = U -1 I n (32)
Proposition 3 (Invariance to scaling) A scaling with powers of 2 (U diagonal with U ii =2 pi ,p i ∈ Z, 1 i n) does not change the transfer function error σ 2 ∆h .

Proof: Let F 2 (x) denotes the fractional value of log 2 |x|:

F 2 (x) log 2 |x|-log 2 |x| (33)
Then the operator . 2 satisfies

ab 2 = a 2 b 2 2 F2(a)+F2(b) (34)
and hence

(Z 1 ) ij 2 = T -1 ii 2 (Z 0 ) ij 2 T jj 2 Φ ij (35) with Φ ij 2 j F2(T -1 ii )+F2(T jj )+F2((Z0) ij ) k . So, Ξ Z Z1 is deduced from Ξ Z Z0 by Ξ Z1 ij = Ξ Z0 ij T -1 ii 2 T jj 2 Φ ij . (36)
The similarity on Z 0 changes the sensitivity such that

∂h ∂Z Z1 = T ∂h ∂Z Z0 T -1 , (37)
then

∂h ∂Z ij Ξ ij Z1 = ∂h ∂Z ij Ξ ij Z0 T -1 ii 2 T -1 ii T jj 2 T jj Φ ij .
(38) Now we can remark that Φ ij ∈{ 1, 2, 4} and Φ ij =1if the power of 2 are used for the scaling. Also a 2 a =1if a is a power of 2.

IV. OPTIMAL REALIZATIONS

It is now possible to consider the optimal transfer function error problem. It consists of finding the optimal coordinate transformation U opt :

U opt = arg min U invertible σ 2 ∆h (U ). (39)
Since σ 2 ∆h is invariant to power-of-2 scaling, this optimization problem has an infinite number of solutions. So it could be of interest to normalize all the coordinate transforms with regard to an extra consideration. For example, this could be a L 2scaling constraints, even if it is not necessary here.

One possible normalization is to apply on a realization the power-of-2 scaling V:

V ii = (W c ) ii 2 , 1 i n (40)
where W c is the controlability Gramian of this realization. Then, the normalized realization has the following property:

Proposition 4 A normalized realization satisfies the relaxed-L 2 -scaling constraints, i.e.:

1 (W c ) ii < 4, 1 i n. ( 41
)
This relaxed-constraints were proposed in [START_REF] Hilaire | Low parametric sensitivity realizations with relaxed l 2dynamic-range-scaling constraints[END_REF] as an extension of the strict-L 2 -scaling constraints (28) that still prevents the implementation from overflows.

Proof: After the normalization, the new realization will have a controlability Gramian changed into W

c = V -1 W c V -, hence (W c ) ii =   (W c ) ii (W c ) ii 2   2 . ( 42
)
We can remark that ∀x, 1 x x 2 < 2, so the normalized realization satisfies (41). Of course, any other normalization is possible, but this one allows us to use the existing L 2 -scaling or relaxed L 2 -scaling methods [START_REF] Hilaire | Low parametric sensitivity realizations with relaxed l 2dynamic-range-scaling constraints[END_REF], [START_REF] Hinamoto | Minimization of l2 sensitivity of one-and two dimensional state-space digital filters subject to l2-dynamic-range-scaling constraints[END_REF].

Then, the optimal problem can be defined as Problem 2 (normalized σ 2 ∆h -optimal realization) Considering a state-space realization (A, b, c,d), the σ 2 ∆h -optimal realization can be found by solving the following optimization problem:

U opt = arg min U invertible σ 2 ∆h (UV), ( 43
)
where V is a diagonal scaling matrix such that

V ii = (U -1 W c U -) ii 2 (44)
Proof: V is built in such way that the coordinate transformation T = UV applied on the realization performes the normalization. Since the σ 2 ∆h measure is non smooth, this optimization problem can be solved with a global optimization method such as the Adaptive Simulated Algorithm (ASA) [START_REF] Ingber | Adaptive simulated annealing (ASA): Lessons learned[END_REF], [START_REF] Chen | Adaptive Simulated Annealing for optimization in signal processing applications[END_REF]. A gradient-base method such as the quasi-Newton algorithm leads to local optima, which in practice are however not too far from the global optimum. The FWR Toolbox 4 was used for the numerical examples, and few minutes of computation were here required on a desktop computer.

V. N UMERICAL EXAMPLE

Let us consider the filter with coefficients given by the Matlab command butter(4,0.05), and some equivalent (in infinite precision) realizations: Z 1 : the Direct Form II, Z 2 : the balanced realization, Z 3 : the normalized σ2 ∆h -optimal realization (obtained with ASA and proposition 2). The following table gives the σ2 ∆h measure of these different realizations. This could be compared to the a posteriori difference between the transfer function h and the transfer function with quantized coefficients (denoted h † ). The wordlengths used are 16, 14 and 10 bits (but 10 bits are not enough to preserve the stability of the Direct Form II, that cannot be used with so few bits): The realization Z 3 has of course the lowest transfer function error. Even with few bits, the degradation of realization Z 3 remain low, compared to Z 1 and Z 2 . Even if it is non-optimal, Z 2 performs quite well with 16-bit implementation, but is less accurate than Z 3 with less bits. As we can remark, the coefficients of these three realizations have different magnitudes, and that legitimates the use of σ

realization σ2 ∆h ' ' h -h † ' ' 2 16

VI. CONCLUSION

This paper has presented the optimal realization problem in fixed-point context and exhibited a new measure to evaluate the coefficients roundoff errors. Compared to the classical L 2sensitivity measure, the transfer function error is a meaningful measure for every realization, even for non-L 2 -scaled ones. A normalization procedure used to solve the optimal realization problem has been presented with a numerical example, but there is still further work to be done, specially to develop an ad hoc optimization algorithm. This measure could be easily adapted to the Specialized Implicit Framework [START_REF] Hilaire | A unifying framework for finite wordlength realizations[END_REF] that allows to encompass all the existing structures (direct forms, cascade/parallel decompositions, lattices, δ or ρ-operator based realizations, etc.). Some other measures, like the pole-sensitivity measure, could be extended to consider fixed-point coefficients.

+29648• 2

 2 2 ∆h .Here are the 16-bit fixed-point coefficients of Z † 3 (each one has a different binary-point position):-15 +27141•2 -18 +20820•2 -20 -30467•2 -19 -32227•2 -19 +24569•2 -20 +29679•2 -15 +22295•2 -17 -31725•2 -20 +19083•2 -22 -31503•2 -20 -31152•2 -19 +29148•2 -15 +30424•2 -22 -32633•2 -15 +22733•2 -17 +21076•2 -20 -32727•2 -21 +29154•2 -15 -26416•2 -31 +28776•2 -24 -32739•2 -22 -25371•2 -26 -32767•2 -18 +16771•2 -29 at http://fwrtoolbox.gforge.inria.fr

  This work has been partially funded by the NFN SISE project (National Research Network "Signal and Information Processing in Science and Engineering") and the Institute of Communications and Radio-Frequency Engineering of Vienna University of Technology.France. T. Hilaire is now with the Laboratory of Computer Science (LIP6) of the University Pierre & Marie Curie of Paris, France.

1

Digital Signal Processors 2 Field Programmable Gate-Array

In FPGA or ASIC, it is of interest to consider the wordlength as optimization variables, in order to find hardware realizations that minimize hardware criteria like power consumption or surface, under certain numerical accuracy constraints, like L 2 -sensitivity ones[START_REF] Rocher | Fixed-point configurable hardware components[END_REF]. This is not considered in this paper.


Carolina Reyes
email: creyes@nt.tuwien.ac.at

Thibault Hilaire
email: thibault.hilaire@lip6.fr

Christoph F Mecklenbräuker

Distributed Projection Approximation Subspace Tracking Based on Consensus Propagation

We develop and investigate a distributed algorithm for signal subspace tracking with a wireless sensor network without the need for a fusion center, in order to improve the robustness and scalability. We assume that all sensor nodes may broadcast messages to sensors in their neighborhood defined by a finite (small) communication radius.

To this aim, we start from Projection Approximation Subspace Tracking (PAST) which is a well-investigated algorithm suitable for implementation in a fusion center. We arrive at a distributed approximation of the PAST algorithm by letting each sensor broadcast its local observation variable xn(t) and a filtered observation vector y n (t) to its neighborhood. Vice versa, the received messages at sensor node n from its neighborhood are fused by employing consensus propagation.

Finally, we investigate the proposed distributed algorithm in simulation runs.

I. INTRODUCTION

Wireless sensor networks consist of a number of nodes which sense environmental changes and report to other nodes inside the network. There are many applications ranging from industrial, building, home system automation to monitoring of chemicals in hydrology, agriculture, and pollution control, as well as the prediction of avalanches and land slides.

In many applications, aggregate functions of the sensor data are more important than individual node data. When the scale of the wireless sensor network is large, the average value collected from the whole network is more important than the value collected by one single node, as stated in [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF], [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]. The relevant subspace in which the sensor data are confined is an important aggregate statistic which is a pre-requisite for many types of data compression techniques, signal detection, classification, and localization.

Sensor nodes are adequate for deployment in harsh environments or over large geographical areas, however most of the state-of-the-art sensor network architectures are rather centralized. The resulting drawbacks motivate us to the current research:

• Since the sensor nodes might be located in remote areas, the use of batteries for powering is necessary. However, when the lifetime of these batteries is used up, then it is very difficult to replace them. Due to this fact, it is a major concern to decrease the power consumption as much as possible.

• The communication channel constrains the data rate for message transmission. We concentrate our efforts on distributing the Projection Approximation Subspace Tracking (PAST) algorithm [START_REF] Yang | Projection approximation subspace tracking[END_REF] with a low amount of message passing among sensor nodes. The latter is motivated by the high power consumption for transmitting messages.

Notation: We denote a column vectors in underlined boldface and matrices in uppercase and boldface. The superscript H represents complex conjugation, transposition, tr(. ) the trace operator, E[. ] the expectation and . the Euclidean vector norm.

Organization of the paper: Section II introduces the PAST algorithm, the analytical model and corresponding constraints. The section III presents a survey on distributed averaging as an alternative for decentralization. The fusion of both PAST and consensus propagation methods in order to obtain a distributed version of the original algorithm is presented in Section IV. Finally, Section V shows and discusses specific simulation results which highlight the tracking capabilities of our distributed approach.

II. PROJECTION APPROXIMATION SUBSPACE TRACKING

In the recent years, subspace methods have become a very researched and remarked method in the realm of modern signal processing. Some of its applications encompasses blind communication channel identification [START_REF] Akhtman | Decision directed channel estimation employing projection approximation subspace tracking[END_REF], estimation of direction of arrivals from signals impinging on an antenna array and image compression, to mention some. Methods based on the singular value decomposition of the sample covariance matrix for every time step lead to a high computational burden, cf. [START_REF] Kumaresan | Estimating the angles of arrival of multiple plane waves[END_REF], [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF]. The PAST algorithm does not require a sample covariance estimate. It estimates the signal subspace at time t recursively, depending on the previous subspace estimate at time t -1 and the new observation x(t). Yang [START_REF] Yang | Projection approximation subspace tracking[END_REF] introduces a new signal subspace model interpretation, considering an unconstrained minimization function. The subspace tracking is implemented by employing the recursive least squares algorithm, based on an appropriate projection approximation.

A. Mathematical model

Let x(t) ∈ C N be the data vector observed at time t, composed of r narrow-band signal waves impinging on a planar array of sensors hidden in additive noise (zero-mean with variance σ 2 ). The data vector is modelled as

x(t)=A(ω)s(t)+v(t) , (1)
A(ω)= a(ω 1 ) a(ω 2 ) ... a(ω r ) , (2)
s(t)=     s 1 (t) . . . s r (t)     , v(t)=     v 1 (t) . . . v N (t)     . (3)
A(ω) is a deterministic N × r matrix depending on ω = (ω 1 ,...,ω r ) whose columns are plane wave steering vectors a(ω i ). The n-th element of a(ω) is

[a(ω i )] n = exp j 2π λ (ξ n cos θ i + η n sin θ i ) / √ N
where frequency ω i = cos θ i . The Cartesian coordinates of the n-th sensor node are (ξ n ,η n ) and λ is the wavelength. N is the number of sensors and r the number of impinging waves. s(t) is the signal vector and v(t) is uncorrelated additive noise.

We allow ω to be slowly time-varying, i.e. ω = ω(t) and seek an estimate for matrix W (t) whose columns span the same subspace as A(ω(t)).Atracking algorithm estimates W (t) by a function of W (t -1) and x(t) alone.

B. Brief review of PAST

Yang [START_REF] Yang | Projection approximation subspace tracking[END_REF] approximates the cost function to minimize

J(W (t)) = t i=1 β t-i x(i) -W (t)W H (t)x(i) 2 (4)
by

J (W (t)) = t i=1 β t-i x(i) -W (t)y(i) 2 (5)
with

y(i)=W H (i -1)x(i) . (6)
Here, the matrix W ∈ C N ×r is constrained to rank r < N . The minimization of cost function ( 5) results in a lowcomplexity update of the signal subspace. The computational complexity is of order 3Nr+ O(r 2 ) operations per time step, where N is the number of sensors and r the number of tracked signals.

In Algorithm 1, y(t) stores the modified data vector resulting from the multiplication between the old signal subspace W H (t -1) and the new data vector x(t) observed at the sensor array. The following computational steps in Algorithm 1 are derived from the Recursive Least Square algorithm [START_REF] Haykin | Adaptive Filter Theory[END_REF] to track the signal subspace matrix W (t).

III. DISTRIBUTED AVERAGING ALGORITHMS

The PAST algorithm is structured in a centralized way, and we aim to distribute it for sensor network applications. Namely, we wish to achieve two important goals for large scale distributed systems: (i) Fault tolerance: consider large scale systems that are prone to communication disruptions due to link/node failures (Robustness). (ii) Simplicity and Algorithm 1: PAST algorithm by Yang [START_REF] Yang | Projection approximation subspace tracking[END_REF] Input: β, P (0), W (0) for t := 1, 2, ... do Input:

x(t) y(t)=W H (t -1)x(t) h(t)=P (t -1)y(t) g(t)=h(t)/[β + y H (t)h(t)] P (t)= 1 β [P (t -1) -g(t)h H (t)] e(t)=x(t) -W (t -1)y(t) W (t)=W (t -1) + e(t)g H (t)
Output: W (t) endfor scalability: pursuit of a simple topology network able to adapt itself according to the volatility of the system. The use of a centralized network limits the scalability of the system.

In recent work, distributed adaptive algorithms have been proposed to address the issue of estimation over distributed networks. These new algorithms outperform the classical centralized solution, but are based on specific network topologies which lead to scalability constrains. For example, [START_REF] Kleinberg | Algorithm Design[END_REF] considers the design of distributed architectures based on randomized graph models.

Different strategies are used for distributed fusion and average of the information. Some algorithms are based on interaction between a node and all its adjacent nodes in the network, so as to reach a consensus, namely the Consensus Propagation algorithms, introduced by Olfati-Saber and Murray in [START_REF] Olfati-Saber | Consensus protocols for networks of dynamic agents[END_REF], [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF].

Another strategy is based on pairwise communications: the gossip-based algorithms such as the push-sum protocol introduced by Kempe in [START_REF] Kempe | Gossip-based computation of aggregate information[END_REF] are good alternatives when low communication overhead is desired.

In this paper, we are using a vector version of the Consensus Propagation, with constant weights (see [START_REF] Moallemi | Consensus propagation[END_REF] for the general case, with non-constant weights), as presented in Algorithm 2.

For every node n =1 , . . . , N in the network, N n denotes the set of the adjacent nodes, including itself. We assume that every node in the neighborhood N n can correctly receive broadcasted messages from node n with probability 1, cf. Figure 1.

At the end of step t -1, every node n sends to its adjacent nodes N n its own estimation of the average y n (t -1) and a weight w n . Then, at the beginning of step t, every node n receives the pairs {y i (t -1),w i } i∈Nn from its neighbors and compute a new average, to be sent at the end of step t. The constant weights are chosen as

w n =1/ |N n | . (7)
Algorithm 2 can be rewritten as follows. Let Y (t) ∈ C N ×r be a matrix aggregating all y n (t) Matrix Y (t) is evolving according to Y (t + 1) = ΩY (t)

Y (t) y 1 (t) ... y N (t) . (9
where

Ω = diag Aw -1 A diag(w), (10)
w = w 1 ,...,w N and A is the adjacency matrix of the undirected graph (A i,j =1if nodes i and j can communicate). The speed of the convergence is given by the second largest eigenvalue of Ω [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF].

IV. PAST-CONSENSUS PROPAGATION ALGORITHM Recall the centralized Algorithm 1, that is employed within each node of the system for distributing the computations, with minor differences that are explained below. First, considera wireless sensor network with a planar square array structure and a number of nodes N . Every node n observes (environmental) data x n (t), assumed to be scalar, that are shared with all nodes j ∈N n inside a specific broadcasting neighborhood. This means that each node n receives the observations from its adjacent nodes x j (t -1), that are used for aggregating the local observation vector x n (t). We define the local observation vector x n (t) as the aggregation of {x j (t -1)} j∈Nn , i.e.

x n (t)=S n x(t -1) [START_REF] Moallemi | Consensus propagation[END_REF] where the |N n |× N selection matrix S n is defined by

(S n ) ij = 1 if j = i-th node ∈N n 0 , (12)
and x(t) is the data vector observed in the whole network. This information exchange permits that all nodes involved in the communication process have knowledge of the neighboring observations. We propose to locally average the vector y n (t) in node n by fusing information aggregated in its associated neighborhood N n . The local averaging is described by Eq.( 8) within Algorithm 2. Here, we calculate the internal y n (t) taking into account the weighted y j (t -1). This is subsequently divided by the sum of the accumulated weight of w n and the neighboring weights w j . Thus, a distributed average version of y n (t) is assured by means of the Consensus Propagation algorithm.

After defining in which part of the algorithm the distribution actually takes place, we introduce the Algorithm 3. The local vector y n (t) is calculated at every node of the system using the Algorithm 2 and the following calculations follow the original Algorithm 1 except for the diagonal weighting matrix D n = S n diag(w)S n . The reason why we decided to distribute this particular variable is because according to (6), this vector contains information from the updated signal subspace at t -1 as well as new arriving observation data x n and x j . As a result, every sensor in the system has indirect knowledge of these parameters.

Algorithm 3: PAST-Consensus

Input: β, P 1 (0),...,P N (0), W 1 (0), . . . , W N (0)

for t := 1,

V. SIMULATIONS

We simulate a planar sensor network with N = 36 nodes, cf. Figure 1. The locations of the nodes are chosen as an imperfect Cartesian grid. The average distance between nodes is λ/2 and the transmission range is 1.44λ/2.

The sensor network observes 1000 snapshots in time. The sensor network tracks the instantaneous frequencies of r =2 impinging waves which are linearly time-varying. The true instantaneous frequencies are shown in the figures below. The forgetting factor β =0.97 is chosen. We initialize W (0) and P (0) to identity. We assume that messages sent by node n are received correctly by all nodes in neighborhood N n . For displaying the subspace tracking behavior vs. time, we apply the spectral MUSIC algorithm [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] for extracting r frequency components from the signal subspace matrices W (t) over time t. All simulation runs are carried out with one and the same noise realisation.

In Figure 2, we observe the tracking of the original PAST Algorithm 1. The true instantaneous signal frequencies are represented by blue lines. The green dots correspond to the tracked frequencies. Likewise, in Figure 3, we plot the tracking behavior of the original PAST algorithm for the subarray defined by the neighborhood N 17 which contains 6 sensors. Finally, in Figure 4, we illustrate the tracking behavior of Algorithm 3. These preliminary results indicate that its behavior is closer to Figure 2 than Figure 3.

VI. SUMMARY AND CONCLUSION

PAST-Consensus (Algorithm 3) implements a distributed variant of the PAST algorithm. Every sensor node n employs the PAST algorithm (Algorithm 1) to locally track the signal subspace and uses Algorithm 2 to couple its own observation and internal state vector with its neighborhood N n of sensor nodes. To this aim, every sensor node n broadcasts its local observation x n (t), its locally filtered r-dimensional vector y n (t), and a weight w n to its neighborhood N n of sensor nodes. Thus, every node broadcasts 2(r + 1) + 1 real-valued variables per time step. From these preliminary results, we believe that a well performing low-cost implementation of a distributed PAST algorithm with low communication overhead is a feasible goal.

9 Fig. 1 .

 91 Fig. 1. Sensor network with neighborhood N 17 for radius 9.

Fig. 2 .Fig. 3 .

 23 Fig. 2. PAST result for whole sensor array (N = 36, r =2)

Fig. 4 .

 4 Fig. 4. Algorithm 3 result for sensor No. 17 (r =2)

ACKNOWLEDGMENT

  Support was provided by the fFORTE-Women in Technology Program at Vienna Univ. of Technology, the Austrian Ministry for Science and Research, the Christian Doppler Lab Wireless Technologies for Sustainable Mobility, and the WWTF project Distributed Information Processing for Spatio-Temporal Fields in Wireless Sensor Networks. We thank Erwin Riegler (ftw) and the national research network on Signal and Information Processing in Science and Engineering (NFN-SISE) for many fruitful discussions.

  Broadcast the pair {y n (t),w n } to all nodes in N n Output: y n (t) is the estimation of the average in step t at node n endfor

	w j	(8)
	j∈Nn

) Algorithm 2: Consensus Propagation Algorithm for n := 1, 2, . . . , N do Input: {y j (t -1),w j } j∈Nn are the pairs sent to node n in step t -1 y n (t)= j∈Nn y j (t -1)w j /


Hristu Culetu
email: hculetu@yahoo.com

On a curved geometry for a rotating frame

Keywords: rotating frame, imperfect stress tensor, energy conditions, curved geometry, flux density

We made a slight modification of the flat metric for a rotating observer, to keep track of the energy spent to generate the rotational motion. The source of curvature energy-momentum tensor represents an imperfect fluid and has constant components, depending on a parameter b which, for physical reasons, is taken equal to the angular velocity ω. The energy density and all pressures are non-negative in the first version of the model from Sec.4. The improved version (Sec.5) gives more realistic expressions of the main quantities of the physical system (though the dominant energy condition is not satisfied), with the help of a positive cosmological constant included in Einstein's equations.

Introduction

In a paper on a relativistic rotating disk [1], Weber presents the dominant view that one needs simply apply traditional relativistic concepts directly and all problems and paradoxes disappear. However, Klauber [2] considers that the issue is far from settled, with a lot of inconsistencies unresolved by the standard solution (see also [3]). Klauber emphasizes that the disk surface experiences no tidal stresses and this is true for any observer, an opinion which is at odds with the traditional point of view.

The Sagnac experiment [4,5] is a Michelson-Morley type experiment for rotating motion and it showed that the local speed of light in a circumferential direction on rotating frames is not invariant [3]. In Klauber's view [2], we have to start with Sagnac's results, not those of Michelson-Morley and rederive relativity theory for rotating frames, following the same steps Einstein did for rectilinear motion. In that case there is no Lorentz contraction along the disk rim, the disk surface is Minkowskian and the clock syncronization turns out to be exactly what one finds by using light rays from a clock located at the disk center.

Nikolic [6] remarked that each point of a rotating disk belongs to a different noninertial frame, although there is no relative motion among different points on the disk. This is similar with the situation faced by a static Rindler observer located at different points on a rotating platform. The author of [6] showed further that, contrary to the prediction of the standard approach, one finds that an observer on the rim of the disk will see equal lengths of other different moving objects as an inertial observer whose instantaneous position and velocity are equal to that of the observer on the rim. Moreover, he also note that the general accepted expression for the space line-element in a non-time-orthogonal frame is found sometimes inappropriate. One also remarks that, according to Debski et al. [7], only motion with constant velocity is relative: acceleration is absolute (and, in particular, rotation) and the principle of relativity does not apply to accelerating motion.

It is our purpose in this work to take advantage of this property of acceleration and look for a different (i.e. curved) geometry in a rotating reference system. That means the new spacetime cannot be obtained from the flat one by a coordinate transformation. In Sec.2 we present the curved metric proposed to be valid for a rotating observer. Sec.3 is devoted to the kinematical quantities associated to a static observer including the covariant acceleration, shear and vorticity tensors. Sec.4 investigates the properties of the source of curvature stress tensor, including the energy conditions. A new version with a cosmological constant gives more realistic results in Sec.5. The Sagnac effect is briefly analysed in Sec.6. The paper concludes in Sec.7 with few final remarks.

Throughout the paper we use geometrical units c = G = 1, unless otherwise specified.

Rotating curved metric

Let us consider the line element of Minkowski spacetime in cylindrical coordinates

ds 2 = -dT 2 + dR 2 + dZ 2 + R 2 dΦ 2 . (2.1)
The motion of a reference system, rotating with constant angular velocity ω is obtained by the coordinate transformation [5]

t = T, r = R, z = Z, ϕ -ωt = Φ, (2.2)
which yields

ds 2 = -(1 -ω 2 r 2 )dt 2 -2ωr 2 dtdϕ + dr 2 + dz 2 + r 2 dϕ 2 . (2.3)
The line-element (2.3) is the well-known Born metric. The time-time component of the metric g tt = -1 + ω 2 r 2 is vanishing at r = 1/ω, which is the radius of the light cylinder. ωr being a velocity, we have to take v = ωr < 1. In addition, one thus avoids the signature flip when r > 1/ω. The light cylinder is not a horizon because [8,9]

g tt - g 2 tϕ g ϕϕ = -1 ̸ = 0. (2.4)
Our main conjecture is that any accelerating observer finds in a curved spacetime, generated by the agent who creates the acceleration. That should be valid both for linear and for rotating systems. In other words, the geometry in a rotating frame cannot be obtained by a coordinate transformation from the flat metric. Therefore, we propose the following spacetime, supposed to be valid for a rotating observer

ds 2 = -(1 -ω 2 r 2 )dt 2 -2ωr 2 dt dϕ + dr 2 1 -b 2 r 2 + dz 2 + r 2 dϕ 2 , ( 2.5)
where b a positive constant (an inverse length) and r < 1/b. We take again r < 1/ω, to get a positive -g tt . It is worth noting that for ωr << 1, b • r << 1 the metric (2.5) acquires a Minkowskian form, in cylindrical coordinates. Same is valid if ω = b = 0 (no rotation).

As far as the curvature invariants (the scalar curvature and the Kretschmann scalar) are concerned, the geometry (2.5) yields

R a a = 2b 2 , K = 4b 4 , ( 2.6)
G t t = G z z = -b 2 , G ϕ t = -b 2 ω. (2.7)
As for the curvature invariants, we see that all the above components are constants, being zero for b = 0.

Kinematical quantities

We now consider a static observer in the geometry (2.5), having a velocity vector field, irrespective of the value of b

u b = ( 1/ √ 1 -ω 2 r 2 , 0, 0, 0 ) , u b = ( - √ 1 -ω 2 r 2 , 0, 0, - ωr 2 √ 1 -ω 2 r 2 ) , (3.1)
where u b u b = -1 and the Latin indices take the values (t, r, z, ϕ). The only nonzero component of the covariant acceleration

a b = u a ∇ a u b is a r = - 1 -b 2 r 2 1 -ω 2 r 2 • ω 2 r, (3.2)
which leads to the proper acceleration

A ≡ √ a b a b = √ g rr |a r | = (1-b 2 r 2 )ω 2 r/(1- ω 2 r 2 ) 3/2
). Because a r is always negative, an observer at rest at r = const. should push oneself to the axis of rotation for to maintain the static position, an effect of the centrifugal force. We know from experiments that, for a rotating observer with ωr << 1, (2.4) gives us a r = -ω 2 r. Therefore, we will consider b = ω from now on, such that Eq.(2.5) becomes

ds 2 = -(1 -ω 2 r 2 )dt 2 -2ωr 2 dt dϕ + dr 2 1 -ω 2 r 2 + dz 2 + r 2 dϕ 2 , ( 3.3)
Consequently, the source of curvature comes directly from the angular velocity.

In other words, the agent who generates the rotational motion creates the energy leading to curvature. It is worth noting that the value of a r given by (2.5) is valid for any ωr < 1, in contrast with a r = -ω 2 r given by (2.4), that is valid for ωr << 1.

The velocity vector given by (3.1) leads to the vanishing of the scalar expansion Θ = ∇ a u a = 0 and of the shear tensor of the particles worldlines

σ a b = 1 2 (h ac ∇ c u b + h c b ∇ c u a ) - 1 3 Θh a b + 1 2 (a a u b + a b u a ), (3.4)
where h a b = g a b + u a u b is the projection tensor onto the direction perpendicular to u a and σ a b expresses the distorsion of the worldlines in shape without change in volume. As far as the vorticity tensor

Ω a b = 1 2 (h ac ∇ c u b -h c b ∇ c u a ) - 1 2 (a a u b -a b u a ) (3.5)
is concerned, the only nonzero components are

Ω t r = ω 2 r (1 -ω 2 r 2 ) 3/2 , Ω r ϕ = ωr √ 1 -ω 2 r 2 , Ω ϕ r = - ω r √ 1 -ω 2 r 2 , ( 3.6)
with the modulus

Ω ≡ √ Ω a b Ω b a = √ 2 ω √ 1 -ω 2 r 2 . ( 3.7)
In contrast, the vorticity vector ω a has only one nonzero component:

ω z = Ω/ √ 2.

Stress tensor properties

We look now for the source of curvature of the geometry (3.3), namely we need the components of the stress tensor to be inserted on the r.h.s. of Einstein's equations for to get (3.3) as an exact solution. From (2.7) we immediately get the nonzero components

8πT t t = 8πT z z = -ω 2 , 8πT ϕ t = -ω 3 . (4.1)
It is clear that we are dealing with an imperfect fluid, with different transversal pressures. Therefore, we start with the general form [10,11,12]

T a b = (p r + ρ)u a u b + p r g a b + (p ϕ -p r )n a n b + (p z -p r )s a s b + u a q b + u b q a , (4.2)
where ρ is the energy density of the fluid, p r , p ϕ , p z , q a are, respectively, the radial pressure, the pressure along the ϕ-direction, the pressure along the zdirection and the flux density vector. In addition,

s a = (0, 0, 1, 0), n a n a = s a s a = 1, u a n a = u a s a = n a s a = 0. (4.3)
From the conditions u a n a = 0, n a n a = 1 and the velocity vector from (3.1), one obtains

n a = ( - ωr √ 1 -ω 2 r 2 , 0, 0, √ 1 -ω 2 r 2 r ) , n a = ( 0, 0, 0, r √ 1 -ω 2 r 2 , ) , (4.4)
Once we acquired the above informations of the kinematical quantities we proceed solving the Einstein equations for to obtain the energy density, pressures and the flux density. By means of T a b from (4.1) and (4.4), we have

8πρ = 8πT a b u b u a = ω 2 1 -ω 2 r 2 , p r = T r r = 0, 8πp z = 8πT z z = -ω 2 . (4.5)
The pressure on the ϕ-direction is obtained from the trace 8πT

a a = -R a a = -2ω 2 8πT a a = 8π(-ρ + p r + p z + p ϕ ) = -2ω 2 , ( 4.6)
whence we have

8πp ϕ = ω 2 r 2 8πρ = ω 4 r 2 1 -ω 2 r 2 .
(4.7)

The flux density q a (with q a u a = 0) results from a multiplication of (4.2) with the velocity vector

q a = -T a b u b -ρu a , ( 4.8)
and we find that

q a = ( - ω 4 r 2 (1 -ω 2 r 2 ) 3/2 , 0, 0, ω 3 √ 1 -ω 2 r 2 ) , q a = ( 0, 0, 0, ω 3 r 2 (1 -ω 2 r 2 ) 3/2
) , (4.9) with q ≡ √ q a q a = ω 3 r/(1 -ω 2 r 2 ). One could check that all the other Einstein equations (when the Einstein tensor is vanishing) are satisfied identically.

As long as the energy conditions for T a b are concerned, one observes that the weak energy condition (WEC) (ρ > 0, ρ + p r > 0, ρ + p t > 0), the null energy condition (NEC) ( ρ + p r > 0, ρ + p t > 0), the strong energy condition (SEC) (ρ+p r > 0, ρ+p t > 0, ρ+p r +2p t > 0) and the dominant energy condition (DEC) (ρ > |p r |, ρ > |p t |) are satisfied. That is very clear from the expressions (4.5) and (4.7) for the energy density and pressures, all being non-negative excepting p z .

Let us make few general observations on the above results. The negative (and very high) value of p z = -ω 2 seems to be unrealistic. In addition, the expression of ρ from (4.5) leads to huge values for ρ. The fact that p r = 0 is also not believable. For those reasons, we look for a more convincing version.

Positive cosmological constant

To get rid of the above inconsistencies, we assume a cosmological constant is generated by rotation, besides the above stress tensor. In this case, the Einstein equations are given by

G a b + λδ a b = 8πT a b , λ > 0. ( 5.1)
The new expressions of the energy density and pressures appears as

8πρ = ω 2 1 -ω 2 r 2 -λ, 8πp r = λ, 8πp z = -ω 2 +λ, 8πp ϕ = ω 4 r 2 1 -ω 2 r 2 +λ (5.2)
(Noting that the expression of the flux density does not depend on λ). We impose a value of the cosmological constant such that p z = 0, a condition leading to λ = ω 2 . In that situation the new values from (5.2) acquire the form

8πρ = ω 4 r 2 1 -ω 2 r 2 , 8πp r = ω 2 , p z = 0, 8πp ϕ = ω 2 1 -ω 2 r 2 (5.3)
We see that 8πρ equals the proper acceleration of our static observer (in the rotating frame). A comparison of the energy conditions validity (with the previous ones from Sec.4) shows that all the energy conditions for the stress tensor are obeyed, excepting the DEC, because in this case ρ < |p ϕ |. Let us consider a numerical example. Take a cylinder with height ∆z and radius r = 1cm, rotating around its central axis with angular velocity ω = 100rad/s. Taking care of all fundamental constants and consider ωr << c, we have, for a rotating observer inside the cylinder at r = 1cm from the axis, ρ ≈ (c 4 /8πG) • ω 4 r 2 /c 4 = ω 4 r 2 /8πG = (a r ) 2 /8πG = 6 • 10 15 ergs/cm 3 , which is a reasonable value. The fact that ρ ∝ (a r ) 2 is a well-known property of a Newtonian gravitational field [13] or of the electrostatic field. The above value of ρ corresponds to a radial acceleration a r ≈ 100 • g, g being the Earth gravity on its surface.

Sagnac effect

Let us consider now two light beams counter-propagating on a circular trajectory (using a system of mirrors) in a rotating reference frame, in flat spacetime [5,14,15]. Due to the rotation of the source of light, we know that a phase shift (with respect to the interference pattern when the device is at rest) and, from here, a time delay is measured (Sagnac effect) between the two counterpropagating beams. Using an analogy with the expressions for the Coriolis and Lorentz forces, Rizzi and Ruggiero [14] showed that the effect is also valid for matter beams : Cooper pairs, charged particles, etc.

Classically, he predicted the following phase shift, for monochromatic light waves in vacuum, counter-propagating along a closed circular path in a rotating interferometer, with radius r 0 and z = const..

∆Φ = 4ωS λc , ( 6.1)
where ω is the angular velocity of the turntable, S = πr 2 0 is the area enclosed by the light path, and λ is the wavelength of light in vacuum. The corresponding time difference will be given by ∆t = λ∆Φ = 4ωS. (6.2)

Noting that the Eq.( 6.2) could be derived from (2.3) where we place r = r 0 = const., z = const. and ds 2 = 0. Solving the equation obtained [16], we get

dϕ dt = ω ± 1 r 0 , ( 6.3)
where (±) refers to the two counter-propagating beams. Integrating the above equation for a complete rotation, (

= 4ωπr 2 0 1 -ω 2 r 2 0 , ( 6.3) yields ∆t
which is (6.2) with the approximation used by Sagnac (ωr << 1). We observe that our metric (2.5) is different from (2.3) only at the term dr 2 /(1 -ω 2 r 2 ). This term does not contribute in the Sagnac experiment, where one considers r = r 0 = const.. Therefore, our geometry (2.5) is not in contradiction with the results of the Sagnac experiment.

Conclusions

We take into account in this short paper the energy spent to generate the rotational motion, on the geometry of the spacetime. As a reward of that energy consumption, the metric becomes curved, with constant curvature invariants and constant components of the energy-momentum tensor. Because of an out of diagonal component, the selected static observer is not comoving with the imperfect fluid.

In the 1st version of the model from Sec.4, an unrealistic value of the energy density of the fluid is obtained, together with a negative pressure along the z-direction. Therefore, a better version is proposed in Sec.5, with reasonable values of the main physical quantities, by means of a cosmological constant (related to the angular velocity of rotation), introduced in Einstein equations. That would be a proof of a local role played by the cosmological constant.

  which have positive constant values. Of course, when b = 0, they are vanishing because the metric becomes flat. With the help of the Maple software package one may compute the components of the Einstein tensor from Einstein's equation G ab = 8πT ab . The only nonzero components of the Einstein tensor look as


Reviewed by Eric Bouteiller (PhD candidate EPHE-PSL, CRCAO)

This timely book shines a new light on archaeology of the Han period (206 BCE-220 CE), collecting almost all the recently published data concerning aboveground funerary features of tombs left by Han era social elites. It includes numerous highquality maps, drawings, and photographs from published excavation reports. The focus of the book is an area often neglected by archeologists in China, who are usually fixated on underground architecture -and the "treasures" therein -in isolation from above-ground facilities, thus de-contextualizing the tombs. Han period graves are very complex structures with often lavishly decorated underground chambers, and features above ground that include all sorts of funerary monuments, such as memorial halls and offering shrines, but also a wider range of facilities built within extensive cemetery gardens.

Liu Zunzhi received his PhD in archeology from Zhengzhou University in 2007 and has authored or co-authored over 90 archaeological reports and research articles, mostly focusing on Han Dynasty burials. He started his archeological career at the Xuzhou Museum (Jiangsu). Xuzhou, the birthplace of the founder of the Han Dynasty Liu Bang, is one of the richest cities for Han Dynasty remains, at the same level of the two Han capitals (Chang'an and Luoyang). Liu Zunzhi is now a professor of archeology at Nankai University in Tianjin. This new center is developing capabilities in archeology and hopes to catch up with heavyweights in archeology from Beijing, Xi'an, Jilin, and Nanjing.

The remarkable growth of the Chinese economy has led to a "golden age of archeology". [START_REF] Rowan | It's a golden age for Chinese archaeology -and the West is ignoring it[END_REF] This is due both to the boom in infrastructure development and real estate construction resulting in extensive excavation activity and an increased awareness of, and new resources provided for, archeological work in recent years. There is more to read today regarding archaeology in China than for any other part of the world, with an "ocean of findings and excavation reports". [START_REF] Sebillaud | L'archéologie en Chine: histoire, fonctionnement et perspectives[END_REF] While there are no official statistics, it has been estimated that more than 50,000 Han tombs have been discovered since 1949. [START_REF] Qingzhu | 劉慶柱 estimates that 30,000 Han tombs have been discovered from 1949 to 1999[END_REF] However, only a fraction of related reports have been published. In the introduction, Liu explains that he spent nearly a decade collecting data for this study. He estimates, thanks to abundant historical records (p. 212), that during the entire Han period there were 1688 nobles (liehou 列侯) who were granted territories and buried in large tombs.

However, only 140 of these tombs have been discovered so far. Thousands of other tombs, the resting places of those from the lower ranks of society, have been excavated in each province. [4] Facing this "ocean of data," there is an urgent need for synthesis and theorization. Chinese archeologists have already done some very promising studies in this regard, but these are mostly limited to specific regions and the findings are often hidden away in unpublished academic theses. Liu Zunzhi's book is an essential step in synthesizing the numerous excavation reports to get a fuller overview of the four centuries of the Han period.

Liu does not include those unexcavated imperial mausoleums with significant above-ground features which have previously been studied in detail. His focus is on the kings (zhuhou wang 諸 侯王) and the aristocracy (liehou) whose members had the social status and economic resources to build mausoleums composed of vast underground multi-chamber tombs and above-ground cemetery gardens. He also presents his findings on the social class below the aristocracy, with "small and medium tombs." The analysis of this data constitutes only one tenth of the book (75 pages), but it is very important for understanding the burial practices of those affluent social classes not belonging to the aristocracy.

Beyond the simple collection of data, this book is groundbreaking in its focus on above-ground architectural features. Chinese archeologists are usually concerned only with rescuing artifacts and generally disregard anthropological research questions. Objectively, Han tombs are very rich in art treasures, and the official policy has been to retrieve these artifacts. As one archaeologist told me, "We only get on stage after the bulldozers have hit a tomb, so nothing remains on the surface." However, the above ground context provides numerous indications about what lies underground. A recent publication on the discovery of Baling 霸陵, the tomb of Emperor Wen of the Han 漢文帝 (202-157 BCE), is a good example of how the quality of information obtained when actively looking at above-ground features can improve archeologists' capabilities in deducing the location of a mausoleum and the identity of the owner. [START_REF] Yongying | Han Wendi Baling yizhi yanjiu 漢文帝霸陵遺址研究[END_REF] At first glance, Liu's study seems limited in scope, but it collects all the details that enable us to understand the "system" built around a tomb and go beyond mere treasure hunting. These details make it possible to re-contextualize a tomb. The author describes many types of remains that can be shortlisted as follows:

 enclosures around cemetery gardens, such as walls, fences or ditches;

 entrances of the cemeteries with que 闕 (stone pillar gates), roads, and stone sculptures;

 pits and trenches accompanying the tombs, filled with offerings, such as cattle or yong 俑 (figurine servants) serving the tomb owner after his death;

 memorial halls and offering shrines used during burials and later rituals;  hydraulic systems to collect and drain the runoff water inside the tombs or the funerary enclosures. The largest mausoleums often include ponds and pools used to store water and maybe to decorate the gardens. There were also wells, which provided drinking water for the people working on-site and perhaps also served ritual purposes.  cemetery garden, with tumuli, the natural topography, artificial landscapes and vegetation.  workshops, kilns, and stone quarries used to produce building materials. Sometimes, tombs and graveyards are enormous and require a large amount of pottery, tiles, bricks or stones which had to be locally produced.  human settlements, including housing for people guarding and maintaining the site.

These are difficult to analyze as the remains are often only tiles and pillar bases. In some cases, temporary installations for the workers who built the tomb are also recorded.

The author's list of above-ground features is quite exhaustive, but one dimension is missing. That is the tombs' localization and the topographical context which can help answer the key question of how the burial site was selected. From a traditional Chinese perspective, it would be interesting to compare the location with Chinese geomantic practices (fengshui 風水). The question is not a simple one. In the same graveyard, there can be tombs of almost the same period which are oriented in exactly opposite directions. No unique model applies. Perhaps some of the variations are due to the families' choices. From an urban planning point of view, the location of tombs could also provide an important clue about the social practices of the Han period. Tombs were located outside of cities, but how far away? Were they in "forbidden" areas, or, on the contrary, placed in open and visible locations? Many more questions could be raised about, for example, the organization of graveyards. At some sites, hundreds of tombs from the Han era can be found in one place, but the logic of their spatial organization remains to be studied.

The structure of the book is a compromise among several perspectives. Liu starts with the tombs of the kings. Due to the abundance of materials available, this part is split into two different chapters focusing on the Western Han period (202 BCE-8 CE) and the Eastern Han period (23-220 CE). Chapter three presents the tombs of the nobility. Chapter four focuses on the expansion of sacrificial infrastructures, including offering shrines and memorial halls for the elite as well as small and medium tombs throughout the period under consideration. According to Liu, these developments are probably due to the increase of economic resources and the evolution of practices and beliefs. Chapter five describes above-ground remains of small and medium tombs. It is very difficult to tell who the owners of these tombs were (as opposed to those tombs described in the previous chapters). Liu Zunzhi is aware of the importance of social context in archeological work and research. His PhD thesis focused on the relationship between tombs and society during the Han period. However, on this question, Liu is constrained by the available sources. Most excavation reports have limited information on tomb owners because the archeologists did not spend enough time on this question.

Chapter six explains the "system" behind the structures found above ground. However, this attempt stops at the descriptive stage. Liu merely states that there were many different practices, which varied according to the social origin of the tomb owner, the location, and the time of the burial. He does not attempt to identify "models" of tombs and their above-ground features. Liu could have used quantitative data and analysis to discern patterns, for example, in the frequency of each type of remains. The book would have also benefitted from an index to enable the reader to cross search data by function, geographical location, and period. This book does not fully explain the role of above-ground features in funerary behaviors during the Han period. This is not the fault of the author but primarily due to limitations in the data available (from excavation reports) and the methodology used to obtain that data. This shows the current limits of Chinese archaeology as a discipline, which is dominated by relic rescue projects and the saving of artefacts jeopardized by the steamroller of economic development. It is to be hoped that, in the future, Chinese archeologists will be able to spend more of their time and resources to follow stratigraphy and fully explore tomb surroundings to capture the burial contexts better.

In summary, this book offers a substantial and complete database on the context of Han tombs from the social elites that will facilitate to reaching a more complete understanding of funerary practices in the Han period. The implicit message of the author is also a loud call to the Chinese archeology community to "slow down" and consider the above ground features of tombs.

[

  Review] 劉尊志《漢代墓外設施研究:以王 侯墓葬与中小型墓葬為参考》(Liu Zunzhi: Auxiliary Surrounding Structures of Han Dynasty Tombs. Case Studies on Tombs of the Nobility and Small-and Medium-Sized Tombs) Liu Zunzhi 劉尊志, Handai mu wai sheshi yanjiu: Yi wanghou muzang yu zhong-xiaoxing muzang wei cankao 漢代墓外設施研究:以王侯墓葬与中小型墓葬為参考 (Auxiliary Surrounding Structures of Han Dynasty Tombs. Case Studies on Tombs of the Nobility and Small-and Medium-Sized Tombs), Beijing: Kexue Chubanshe 科學出版社, 2021. ISBN 978-7-03-068149-2; 718 pages, 278 RMB.