State Observability in Presence of Disturbances: the Analytic Solution and its Application in Robotics

Abstract-This paper presents the analytic solution of a fundamental open problem in the framework of state estimation/nonlinear observability, which is the Unknown Input Observability problem (UIO problem). The problem consists in deriving the analytic criterion that allows us to automatically obtain the state observability in presence of disturbances (or unknown inputs). In other words, the problem is to extend the well known observability rank condition to the case when the dynamics are also driven by unknown inputs. Enunciated in the seventies by the control theory community, this problem was only solved in the linear case. The solution here provided holds for nonlinear systems in presence of a single unknown input. The first part of the paper presents this analytic solution. Very surprisingly, the complexity of the overall analytic criterion is comparable to the complexity of the observability rank condition. The second part of the paper applies this analytic criterion to a robotics system when its dynamics are affected by an external disturbance (e.g., due to the presence of wind). To corroborate the results of our observability analysis we perform extensive simulations and we show that, a simple estimator based on an Extended Kalman Filter, provides results that agree with what we could expect from the observability analysis.

I. INTRODUCTION

State observability is a necessary condition that a state must satisfy to be estimated. In robotics, where the systems can be often very complex, intuitive reasoning to obtain the observable state can fail. For, performing an observability analysis prior to solving an estimation problem can be fundamental. This has been done by many authors to solve several problems (e.g., in SLAM [START_REF] Bryson | Observability Analysis and Active Control for Airbone SLAM[END_REF], [START_REF] Huang | Observability-based Rules for Designing Consistent EKF SLAM Estimators[END_REF], [START_REF] Huang | A Quadratic-Complexity Observability-Constrained Unscented Kalman Filter for SLAM[END_REF], [START_REF] Lee | On the Observability and Observability Analysis of SLAM[END_REF], [START_REF] Perera | On the linear and nonlinear observability analysis of the SLAM problem[END_REF], [START_REF] Trawny | Observabilitybased Consistent EKF Estimators for Multi-robot Cooperative Localization[END_REF], [START_REF] Vidal-Calleja | On the Observability of Bearing-only SLAM[END_REF], [START_REF] Wang | Observability analysis of SLAM using fisher information matrix[END_REF], in visual-inertial sensor fusion [START_REF] Hesch | Towards consistent vision-aided inertial navigation[END_REF], [START_REF] Hesch | Camera-IMU-based Localization: Observability Analysis and Consistency Improvement[END_REF], [START_REF] Huang | Towards consistent visual-inertial navigation[END_REF], [START_REF] Jones | Visual-inertial navigation, mapping and localization: A scalable real-time causal approach[END_REF], [START_REF] Kelly | Visual-inertial simultaneous localization, mapping and sensor-to-sensor self-calibration[END_REF], [START_REF] Kottas | On the consistency of vision-aided inertial navigation[END_REF], [START_REF] Li | Improving the accuracy of EKF-based visualinertial odometry[END_REF], [START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF], [START_REF] Mirzaei | A Kalman filter-based algorithm for IMU-camera calibration: Observability analysis and performance evaluation[END_REF], [START_REF] Panahandeh | Observability Analysis of a Vision-Aided Inertial Navigation System using Planar Features on the Ground[END_REF], in sensor calibration [START_REF] Censi | Simultaneous calibration of odometry and sensor parameters for mobile robots[END_REF], [START_REF] Chao | IMU-RGBD Camera 3D Pose Estimation and Extrinsic Calibration: Observability Analysis and Consistency Improvement[END_REF], [START_REF] Mirzaei | 3D Lidar-Camera Intrinsic and Extrinsic Calibration: Observability Analysis and Analytical Least Squares-based Initialization[END_REF]). Investigating the observability properties is very simple in the linear case. Unfortunately, real systems are very rarely characterized by linearity.

The control theory community has introduced the analytic criterion necessary to check the state observability for nonlinear systems provided that all the system inputs are known (i.e., in absence of disturbances). This is the observability rank condition introduced by Herman and Krener in 1977 [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF]. In accordance with this criterion, it is possible to obtain all the observability properties of a nonlinear system by performing automatic computation. On the other hand, in many real scenarios, one or more disturbances can dramatically impact the system dynamics. A disturbance can be considered as an unknown input (UI). Its presence can dramatically affect the observability properties of the state. This is for instance the case of a drone that operates in presence of wind (or a boat that operates in presence of stream). The wind (or the stream) is in general unknown, time-variant and acts on the system dynamics as an unknown input.

The problem of finding the analytic criterion able to determine the observability properties in presence of disturbances was enunciated and firstly investigated by the automatic control community in the seventies [START_REF] Basile | On the observability of linear, time invariant systems with unknown inputs[END_REF], [START_REF] Bhattacharyya | Observer design for linear systems with unknown inputs[END_REF], [START_REF] Guidorzi | On Wonham stabilizability condition in the synthesis of observers for unknown-input systems[END_REF], [START_REF] Wang | Dorato Observing the states of systems with unmeasurable disturbance[END_REF]. In particular, Basile and Marro provided the solution of this problem in the linear case [START_REF] Basile | On the observability of linear, time invariant systems with unknown inputs[END_REF]. In robotics, most of the systems are characterized by nonlinear dynamics, even in the very simple case of planar robotics. Additionally, the presence of disturbances cannot be ignored in many cases and can significantly alter the observability properties. For these reasons, this problem in the nonlinear case has been investigated by the robotics community and partial solutions have recently been proposed [START_REF] Belo | 3 Known Landmarks are Enough for Solving Planar Bearing SLAM and Fully Reconstruct Unknown Inputs[END_REF], [START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF]. These solutions are heuristic solutions to the UIO problem and only provide sufficient conditions for the state observability. They are based on a suitable state extension.

This paper provides the analytic solution of the UIO problem in the nonlinear case when the dynamics are driven by multiple known inputs and a single unknown input. A preliminary solution of this problem was presented in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF]. Here, this result is completed by providing all the convergence properties of the algorithm in [START_REF] Martinelli | Nonlinear Unknown Input Observability: Analytical expression of the observable codistribution in the case of a single unknown input[END_REF]. Finally, the general analytic solution that holds for any nonlinear system driven by any number of disturbances is provided in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF].

In [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF], [START_REF] Isidori | Nonlinear Control Systems[END_REF] the observability properties of a nonlinear system are obtained by computing the observable codistribution 1 . The computation of this codistribution is the core of the observability rank condition introduced in [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF]. In order to deal with the case of unknown inputs, we need to derive a new algorithm able to generate the observable codistribution. In the first part of section II we remind the reader the algorithm to compute the observable codistribution in the case without disturbances. Then, we introduce the new algorithm that generates the entire observable codistribution in presence of a single unknown input. In section II-B we provide its convergence properties, by a comparison with the case without disturbances. The solution of the UIO problem in the case of a single unknown input is summarized in section II-C. Then, the paper applies this criterion to obtain the observability properties of a robotic system. In section III we define our robotics problem. It consists of a robot that moves in 3D. This could be an aerial vehicle equipped with airspeed, gyroscopes and monocular vision sensors (but it could also be an underwater vehicle). We assume that the vehicle dynamics are affected by an external disturbance (e.g., due to the presence of wind). We apply the proposed analytic criterion to obtain the observability properties of this system. In section III-D we perform extensive simulations and we introduce a simple estimator. The results obtained by estimating the state with this estimator fully agree with our observability analysis. In section IV we provide our conclusion. We conclude this section by also mentioning new important analytic results recently obtained.

II. EXTENSION OF THE OBSERVABILITY RANK CONDITION

We refer to a nonlinear control system with m u known inputs (u

[u 1 , • • • , u mu ] T
) and a single unknown input or disturbance (w). The state is the vector x ∈ M , with M an open set of R n . We assume that the dynamics are nonlinear with respect to the state and linear with respect to the inputs (both known and unknown). Finally, for the sake of simplicity, we will refer to the case of a single output y (the extension to multiple outputs is straightforward). Our system is characterized by the following equations:

     ẋ = mu i=1 f i (x)u i + g(x)w y = h(x) (1)
where f i (x), i = 1, • • • , m u , and g(x) are vector fields in M and the function h(x) is a scalar function defined on the open set M . We denote by x 0 the state in M where we are interested in deriving the observability properties. We denote by the symbol ∇ the gradient with respect to the state x. Additionally, for a given codistribution Ω and a given vector field θ (both defined on the open set M ), we denote by L θ Ω the codistribution whose covectors are the Lie derivatives along θ of the covectors in Ω. Finally, given two vector spaces V 1 and V 2 , we denote by V 1 + V 2 their sum, i.e., the span of all the generators of both V 1 and V 2 .

Note that this section directly provides the analytic results. All the analytic proofs can be found in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF].

A. Analytic computation of the observable codistribution

We start by reminding the reader the algorithm that computes the observable codistribution in absence of disturbances (namely, the core of the observability rank condition introduced by Herman and Krener [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF]). In absence of disturbances, i.e., when also the input w is a known input, the observable codistribution is generated by the following recursive algorithm (see [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF] and [START_REF] Isidori | Nonlinear Control Systems[END_REF]):

Algorithm 1 Observable codistribution when w is known 1) Ω 0 = span{∇h}; 2) Ω m = Ω m-1 + mu i=1 L f i Ω m-1 + L g Ω m-1
In presence of multiple outputs, we only need to add to the codistribution Ω 0 , the span of the gradients of the remaining outputs.

Let us consider now the case when the input w is unknown. We will denote by L 1 g = L 1 g (x) the first order Lie derivative of the function h(x) along the vector field g, i.e.:

L 1 g L g h ∇h • g(x) (2)
The analytic computation of the observable codistribution is based on the assumption that L 1 g = 0 on a given neighbourhood of x 0 . If this hypothesis is not satisfied, it is possible either to redefine the output such that this assumption holds or to prove that the presence of the disturbance does not affect the state observability 2 . Hence, we are allowed to assume that L 1 g = 0 on a given neighbourhood of x 0 . The algorithm that generates the observable codistribution is the following:

Algorithm 2 Observable codistribution when w is unknown 1) Ω 0 = span{∇h}; 2) Ω m = Ω m-1 + mu i=1 L f i Ω m-1 + L g L 1 g Ω m-1 + mu i=1 Li φm-1 ∇h where the vectors i φ m ∈ R n (i = 1, • • • , m u )
are defined by the following algorithm:

Algorithm 3 1) i φ 0 = f i ; 2) i φ m = [ i φm-1, g] L 1 g
where the parenthesis [•, •] denote the Lie brackets of vector fields, defined as follows:

[f, g] ∂g ∂x f (x) - ∂f ∂x g(x)
In presence of multiple outputs, we only need to add to the codistribution Ω 0 , the span of the gradients of the remaining outputs. Note that, in presence of multiple outputs, the function L 1 g is still a scalar function since it is still defined by using a single output. The result is independent of the chosen output (provided that L 1 g does not vanish 3 ). We conclude this subsection by remarking the differences between algorithm 1 and 2. They are:

• In the recursive step, the term L g Ω m-1 must be replaced by L g L 1 g

Ω m-1 . In other words, the vector field that corresponds to the unknown input (i.e., g) must be rescaled by dividing by L 1 g .

• The recursive step also contains the sum of the contributions mu i=1 Li φm-1 ∇h. In other words, we need to compute the Lie derivatives of the gradient of the output h(x) along the vector fields obtained through the recursive algorithm 3.

B. Convergence properties of algorithm 2

We start by summarizing the convergence properties of algorithm 1. In [START_REF] Isidori | Nonlinear Control Systems[END_REF] it is proven that this algorithm converges in an open and dense set of M . In particular, from lemmas 1.9.1, 1.9.2 and 1.9.6 in [START_REF] Isidori | Nonlinear Control Systems[END_REF], we know that the convergent codistribution is obtained at the smallest integer m for which Ω m = Ω m+1 and m ≤ n -1. This analytic result is a consequence of the fact that, all the terms that appear at the recursive step of algorithm 1, are the Lie derivative of the codistribution at the previous step, along fixed vector fields (i.e., vector fields that remain the same at each step of the algorithm).

Let us consider now the recursive step of algorithm 2. As mentioned above, the first difference with the recursive step of algorithm 1, is the term

L g L 1 g Ω m-1 , which replaces the term L g Ω m-1 .
However, this new term is still the Lie derivative of the codistribution at the previous step along the fixed vector field g L 1

g

. Hence, the presence of this term, does not alter the convergence properties with respect to the ones that hold for algorithm 1. On the other hand, the nature of the last term at the recursive step of algorithm 2 (i.e., mu i=1 Li φm-1 ∇h) significantly differs from all the other terms. Indeed, the vector fields along with we are computing the Lie derivatives (i.e., i φ m-1 ) change at each step m of algorithm 2. In other words, they are not fixed vector fields. Additionally, the Lie derivative is not computed on the codistribution at the previous step of the algorithm but directly on the gradient of the output (∇h).

Therefore, we distinguish the case when the contribution of this term is included in the other terms from the case when it is not included. In the former case, the convergence criterion of algorithm 2 would simply consist of the inspection of the equality Ω m+1 = Ω m , as for algorithm 1 and we conclude that, in this case, the convergence of algorithm 2 is attained in at most n -1 steps 4 .

The following result provides a simple criterion to check, by a finite number of operations, if a given system belongs to this former case. Without loss of generality, we consider the case when m u = 1 (we denote by φ m 1 φ m ). We have the following result:

Lemma 1 Let us denote by Λ j the distribution generated by φ 0 , φ 1 , • • • , φ j (i.e., the vectors obtained by running algorithm 3) and by m(≤ n -1) the smallest integer for which Λ m+1 = Λ m (n is the dimension of the state x). If L φj L 1 g = 0, ∀j = 0, • • • , m, the convergence of algorithm 2 occurs in at most n -1 steps and it occurs at the smallest integer j such that Ω j+1 = Ω j .

Proof: See the proof of lemma 6 in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF] Note that the existence of an integer m(≤ n -1) such that Λ m+1 = Λ m is proved in [START_REF] Isidori | Nonlinear Control Systems[END_REF]. In particular, in the first chapter of [START_REF] Isidori | Nonlinear Control Systems[END_REF], lemmas 1.8.2 and 1.8.3 analyze the convergence of Λ j with respect to j. Note that our case corresponds to the case analyzed in [START_REF] Isidori | Nonlinear Control Systems[END_REF] when ∆ is the span of f , and when τ 1 , • • • , τ q are the single vector field g. Additionally, in [START_REF] Isidori | Nonlinear Control Systems[END_REF] there is not the term L 1 g . Namely, our case corresponds to the case analyzed in [START_REF] Isidori | Nonlinear Control Systems[END_REF], when L 1 g = 1. On the other hand, it is immediate to repeat all the steps in lemmas 1.8.2 and 1.8.3 of [START_REF] Isidori | Nonlinear Control Systems[END_REF] with a general expression of L 1 g and prove the following result: The distribution Λ j converges to Λ * on an open and dense set in M and the convergence is achieved at the smallest integer for which we have Λ m+1 = Λ m . Additionally, we have

Λ m+1 = Λ m = Λ * and m cannot exceed n -1.
Let us derive an upper bound for the number of computations necessary to check if the assumptions of lemma 1 are met. We need to check that L φj L 1 g = 0, ∀j = 0, • • • , m, where m is the smallest integer for which Λ m+1 = Λ m . Hence, we need to run at most the first n -1 steps of algorithm 3 to obtain φ

0 , φ 1 , • • • , φ n-1 and compute at most n Lie derivatives of L 1 g , i.e., L φ0 L 1 g , L φ1 L 1 g , • • • , L φn-1 L 1 g
. Therefore, we have: Remark 1 Checking that the assumptions of Lemma 1 are met, only requires to perform a finite number of (simple) operations Let us consider now the case when the assumptions of Lemma 1 are not met. This is the most frequent case. The criterion needs to compute the following key quantity:

τ L 2 g h (L 1 g ) 2 (3)
When the assumptions of Lemma 1 are not satisfied, it is possible to prove that it exists m such that ∇τ ∈ Ω m (and therefore ∇τ ∈ Ω m ∀m ≥ m ). Additionally, the convergence of algorithm 2 has been reached when Ω m+1 = Ω m , and m ≥ m . Finally, in this case, the convergence is attained in at most n + 2 steps5 .

C. Summary of the analytic criterion

We conclude this section by outlining the steps to investigate the weak local observability at a given point x 0 of a nonlinear system driven by a single disturbance and several known inputs. In other words, to investigate the weak local observability of a system defined by a state that satisfies (1). Basically, these steps are the steps necessary to compute the observable codistribution (i.e., the steps of algorithms 2 and 3) and to prove that the gradient of a given state component belongs to this codistribution.

1) For the chosen x 0 , compute L 1 g (= L 1 g h). In the case when L 1 g = 0 redefine the output, as explained in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF].

2) Compute the codistribution Ω 0 and Ω 1 (at x 0 ) by using algorithm 2. 3) Compute the vector fields i φ m (i = 1, • • • , m u ) by using algorithm 3, starting from m = 0, to check if the considered system is in the case dealt by Lemma 1. In this case, set m = 0, use the recursive step of algorithm 2 to build the codistribution Ω m for m ≥ 2, and skip to step 6.

4) Compute τ = L 2 g h (L 1 g ) 2
j , j = 1, • • • , n) belongs to Ω * (namely if ∇x j ∈ Ω * )
on a given neighbourhood of x 0 , then x j is weakly locally observable at x 0 . If this holds for all the state components, the state x is weakly locally observable at x 0 . Finally, if the dimension of Ω * is smaller than n on a given neighbourhood of x 0 , then the state is not weakly locally observable at x 0 .

III. ROBOTIC APPLICATION

A. The system

We consider a vehicle that moves in a 3D-environment. We assume that the dynamics of the vehicle are affected by the presence of a disturbance (e.g., this could be an aerial vehicle in presence of wind). We assume that the direction of the disturbance is constant in time and a priori known. Conversely, the disturbance magnitude is unknown and time dependent. The vehicle is equipped with speed sensors (e.g., airspeed sensors in the case of an aerial vehicle), gyroscopes and a bearing sensor (e.g., monocular camera). We assume that all the sensors share the same frame (in other they are extrinsically calibrated). Without loss of generality, we define the vehicle local frame as this common frame. The airspeed sensors measure the vehicle speed with respect to the air in the local frame. The gyroscopes provide the angular speed in the local frame. Finally, the bearing sensor provides the bearing angles of the features in the environment expressed in its own local frame. We consider the extreme case of a single point feature and, without loss of generality, we set the origin of the global frame at this point feature (see figure 1 for an illustration). Additionally, we assume that the z-axis of the global frame is aligned with the direction of the disturbance. Our system can be characterized by the following state: X [x, y, z, q t , q x , q y , q z ] T [START_REF] Bryson | Observability Analysis and Active Control for Airbone SLAM[END_REF] where r = [x, y, z] is the position of the vehicle in the global frame and q = q t + q x i + q y j + q z k is the unit quaternion that describes the transformation change between the global and the local frame. The dynamics are The z-axis of the latter is aligned with the direction of the disturbance (assumed to be known and constant in time). The speed V is the vehicle speed with respect to the air, which differs from the ground speed because of the disturbance (w).

affected by the presence of the disturbance. The disturbance is characterized by the following vector (in the global frame):

w = w   0 0 1   ( 5
)
where w is its unknown magnitude.

In the sequel, for each vector defined in the 3D space, the pedix q will be adopted to denote the corresponding imaginary quaternion. For instance, regarding the vehicle position, we have: r q = 0 + x i + y j + z k. Additionally, we denote by V and Ω the following physical quantities:

• V = [V x , V y , V z ]
is the vehicle speed with respect to the air expressed in the local frame (hence, w k +qV q q * is the vehicle speed with respect to the ground expressed in the global frame).

• Ω [Ω x Ω y Ω z ] is the angular speed (and Ω q = 0 + Ω x i + Ω y j + Ω z k). The dynamics of the state are:

  ṙq = w k + qV q q * q = 1 2 qΩ q (6)
The monocular camera provides the position of the feature in the local frame (F q = -q * r q q) up to a scale. Hence, it provides the ratios of the components of F :

h cam (X) [h u , h v ] T =
(q * r q q) x (q * r q q) z , (q * r q q) y (q * r q q) z T [START_REF] Chao | IMU-RGBD Camera 3D Pose Estimation and Extrinsic Calibration: Observability Analysis and Consistency Improvement[END_REF] where the pedices x, y and z indicate respectively the i, j and k component of the corresponding quaternion. We have also to consider the constraint q * q = 1. This provides the further observation:

h const (X) q * q (8)
Our system is characterized by the state in (4), the dynamics in (6) and the three outputs h u , h v and h const in ( 7) and (8).

B. Observability in absence of disturbance

Our system is characterized by the state in (4), the dynamics in [START_REF] Guidorzi | On Wonham stabilizability condition in the synthesis of observers for unknown-input systems[END_REF] with w = 0 and the three outputs h u , h v and h const in [START_REF] Chao | IMU-RGBD Camera 3D Pose Estimation and Extrinsic Calibration: Observability Analysis and Consistency Improvement[END_REF] and [START_REF] Hermann | Nonlinear Controllability and Observability[END_REF].

By comparing ( 6) with ( 1) we obtain that our system is characterized by six known inputs (m u = 6) that are:

u 1 = Ω x , u 2 = Ω y , u 3 = Ω z , u 4 = V x , u 5 = V y and u 6 = V z .
Additionally, we obtain:

f 1 = 1 2           0 0 0 -q x q t q z -q y           , f 2 = 1 2           0 0 0 -q y -q z q t q x           , f 3 = 1 2           0 0 0 -q z q y -q x q t           , f 4 =           q 2 t + q 2
x -q 2 y -q 2 z 2q t q z + 2q x q y 2q x q z -2q t q y 0 0 0 0

          , f 5 =           2q x q y -2q t q z q 2 t -q 2 x + q 2 y -q 2 z 2q t q x + 2q y q z 0 0 0 0           , f 6 =           2q t q y + 2q x q z 2q y q z -2q t q x q 2 t -q 2 x -q 2 y + q 2 z 0 0 0 0          
Finally, in absence of disturbance we have: g = [0, 0, 0, 0, 0, 0, 0]

T In this case we can apply the observability rank condition, i.e., algorithm 1, to obtain the observable codistribution. We compute the codistribution Ω 0 by computing the gradients of the three functions h u , h v and h const . We obtain that this codistribution has dimension equal to 3. We use algorithm 1 to compute Ω 1 . We obtain that its dimension is 4. In particular, the additional covector is obtained by the gradient of the following Lie derivative:

L f 4 h u
In other words:

Ω 1 = span ∇h u , ∇h v , ∇h const , ∇L f 4 h u ,
All the remaining first order Lie derivatives have gradient that is in the above codistribution. Additionally, by an explicit computation, it is easy to realize that Ω 2 = Ω 1 . This means that algorithm 1 has converged and the observable codistribution is Ω 1 . By an explicit computation, it is possible to check that the gradients of the components of the vector F belong to Ω 1 . This means that all the observable modes are the components of F , i.e., the position of the feature in the local frame (obviously, the fourth observable mode is the norm of the quaternion). In particular, no component of the vehicle orientation is observable.

C. Observability in presence of the disturbance

We now consider the case when the dynamics are affected by the presence of the disturbance. By comparing ( 6) with ( 1) we obtain that the vector fields that characterize the dynamics are the same that characterize the dynamics in absence of disturbance with the exception of the last one, which becomes: g = [0, 0, 1, 0, 0, 0, 0] T To derive the observability properties we apply the proposed analytic tool, by following the seven steps provided in section II-C. We have:

First Step

We start by computing the Lie derivatives of the outputs h u , h v and h const along the vector field g. We find that the result differs from zero for the first two outputs. Hence, we use the first output (h u ) to define L 1 g (we could choose also the second output h v ). In particular, we obtain: L 1 g L g h u = {-y(2q t q z -2q x q y ) -x(q 2 t -q 2 x + q 2 y -q 2 z )}/{[z(q 2 t -q 2 xq 2 y + q 2 z ) + 2x(q t q y + q x q z ) + 2y(q y q z -q t q x )] 2 }

Second Step

We compute the codistribution Ω 0 by computing the gradients of the three functions h u , h v and h const . This coincides with the case without disturbance, and we obtain that this codistribution has dimension equal to 3.

We use algorithm 2 to compute Ω 1 . We obtain that its dimension is 5. In particular, the additional two independent covectors are obtained by the gradients of the following two Lie derivatives:

L f 4 h u , L g L 1 g h v
In other words:

Ω 1 = span ∇h u , ∇h v , ∇h const , ∇L f 4 h u , ∇L g L 1 g h v
All the remaining first order Lie derivatives have gradient that is in the above codistribution.

Third Step

We compute 1 φ 1 , 2 φ 1 , 3 φ 1 , 4 φ 1 , 5 φ 1 and 6 φ 1 by using algorithm 3. We obtain that all these vectors vanish. As a result, all the subsequent steps of algorithm 3 provide null vectors. Therefore, the assumptions of Lemma 1 are trivially met. We set m = 0 and we skip to the sixth step.

Sixth Step

We use algorithm 2 to compute Ω 2 and we obtain:

Ω 2 = Ω 1 + span ∇L f 4 L g L 1 g h v
Hence, its dimension is 6. Finally, by using again algorithm 2 it is possible to compute Ω 3 and to check that Ω 3 = Ω 2 . This means that the algorithm has converged and the observable codistribution is Ω * = Ω 2 .

Seventh Step

By computing the distribution orthogonal to the codistribution Ω * we can find the continuous symmetry that characterizes the unobservable space [START_REF] Martinelli | State Estimation Based on the Concept of Continuous Symmetry and Observability Analysis: the Case of Calibration[END_REF]. By an explicit computation we obtain the following vector:

-y, x, 0, - q z 2 , - q y 2 , q x 2 , q t 2 T
This symmetry corresponds to an invariance with respect to a rotation around the z-axis of the global frame. This means that we have a single unobservable mode that is the yaw in the global frame 6 . We conclude by remarking that the presence of the disturbance, even if its magnitude is unknown and is not constant, makes observable the roll and the pitch angles. This result is similar to the result that we obtain in the case of visual and inertial sensor fusion in presence of gravity. The presence of gravity makes observable the roll and the pitch angles, even if its magnitude is unknown [START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF].

What it is non intuitive in the case now investigated, is that, not only the magnitude of the disturbance is unknown, but it is also time dependent.

D. Simulations Simulated trajectories and robot sensors

The trajectories are simulated as follows. The equations in (6) are discretized with a time step of 0.01 s. Each trial lasts 200 s. The initial vehicle speed is set to zero. The initial position is set equal to [0.5 0.5 0.5]m. The vehicle motion is randomly generated. The angular speed, i.e. Ω, is Gaussian. Specifically, its value at each step follows a zero mean Gaussian distribution with covariance matrix equal to (1 deg) 2 I 3 , where I 3 is the identity 3×3 matrix. At each time step, the vehicle speed is incremented by adding a random vector with zero mean Gaussian distribution. In particular, the covariance matrix of this distribution is set equal to σ 2 I 3 , with σ = 10 -4 m. Finally, we simulate the disturbance as a vector along the z-axis, whose magnitude is generated as a random Gaussian variable, with mean value 0.2 ms -1 and variance 0.05 2 m 2 s -2 . Typical trajectories, obtained with this setting, are displayed in figures 2 and 3.

The vehicle is equipped with proprioceptive sensors able to measure at each time step the speed with respect to the air and the angular speed. These measurements are affected by errors. Specifically, each measurement is generated at every time step of 0.01 s by adding to the true value a random error that follows a Gaussian distribution. The mean value of this error is zero and the standard deviation is 0.01 ms -1 for the airspeed and 1 deg s -1 for the gyroscope. Regarding the exteroceptive measurements, they are generated at a lower frequency. Specifically, the measurements are generated each 0.1 s. Also these measurements are affected by errors. Specifically, each measurement is generated by adding to the true value a random error that follows a zero mean Gaussian distribution, with variance 1 deg 2 . We simulate two point features, one at the origin and the other one at the position [1 0 1]m (note that the result of the observability analysis is independent of the number of point features).

Estimation results

We adopt an EKF that estimates an extended state that includes the state in (4) and the unknown input (w) together with its first order time derivative (w (1) dw dt ). In other words, the state is: x, y, z, q t , q x , q y , q z , w, w (1) T . Its dynamics are given in [START_REF] Guidorzi | On Wonham stabilizability condition in the synthesis of observers for unknown-input systems[END_REF] with the two further equations ẇ = w (1) and ẇ(1) = w (2) . In order to implement the prediction phase of our EKF we have to provide the value of w (2) ( d 2 w dt 2 ), which is unknown. We set this quantity to zero. Note that the simulated trajectory does not satisfy this hypothesis since the disturbance is randomly generated.

However, the estimator is able to provide good performance as it is shown in figures 4-7. Figure 2 displays the estimated trajectory in absence of the disturbance. The red line is the trajectory obtained by only using the proprioceptive measurements, while the black line is the trajectory obtained by the proposed EKF. The true trajectory is in blue. Both the estimated trajectories diverge. However, the divergence is slower for the trajectory estimated by the EKF. This divergence is consistent with our observability analysis. The vehicle position is not observable (both with and without disturbance). Figure 3 displays the same of figure 2 for a trial in presence of the disturbance. Note that, in this case, the z-axis has now a different scale (much larger). This is to show the trajectory obtained by only using the proprioceptive measurements, whose divergence is much faster due to the presence of the disturbance along the z-axis. In the remaining figures, we show the performance of the proposed EKF in estimating the distances of the two features and the roll, the pitch and the yaw angles. In accordance with our observability analysis, the distances of the features are observable both with and without the disturbance. In contrast, the roll and the pitch angles are observable only in presence of the disturbance. Finally, the yaw angle is always unobservable. All these results are clearly confirmed by our simulation results. Figure 4 shows the error on the two distances, in absence and in presence of the disturbance, respectively. We note that, the error does not increase with time in both cases. Figures 5, 6 and7 show the error on the roll, pitch and yaw angles, respectively. The left side is in absence of the disturbance while the right side is in presence of the disturbance. In absence of the disturbance, the error increases with time for all the three angles. In presence of the disturbance, it does not increase for the roll and the pitch angle. We conclude this section by remarking that our simulations fully agree with the observability analysis provided in section III-C. Also note that we showed the result of a single trial. However, we found similar results by running many trials.

IV. CONCLUSION

In this paper we provided the analytic criterion to perform the observability analysis of a system in presence of disturbances together with an application of this criterion in order to obtain the observability properties of a robotic system. The analytic criterion is the solution of a fundamental open problem in nonlinear observability/estimation theory, which is the nonlinear unknown input observability problem. The first part of the paper provided this solution in the case of nonlinear systems driven by a single unknown input and multiple known inputs. In the second part of the paper we analyzed the observability properties of a system that consists of a robot that moves in 3D and when its dynamics are affected by an external disturbance (e.g., due to the presence of wind). For this robotic system the analytic criterion provided automatically the observable states. To corroborate the results of this observability analysis, we performed extensive simulations. We showed that, a simple estimator based on an Extended Kalman Filter, provides results that fully agree with what we could expect from the observability analysis.

Note that more robotic applications of our analytic solution can be found in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF]. In particular, in section 5.1, we analyze the observability of the well known unicycle when one of the two inputs (which are the linear and the angular speed) is unknown. When the linear speed is unknown and when the observations consist of angle measurements (e.g., the ones provided by a monocular camera), all the measurements only consist of angle measurements (the angular speed and the camera measurements). In this case, it is obvious that, with respect to the case when the linear speed is a known input, the observability space loses one degree, which corresponds to the absolute scale. It is amazing that, our simple analytic solution, automatically detects this unobservable direction (see more details in section 5.1 in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF]).

Finally, in this paper we provided the analytic solution of the UIO problem for driftless systems and in the case of a single unknown input. Recently, the general case that includes a nonlinear drift and multiple unknown inputs was investigated and the analytic criterion to check the state observability of these more general systems has been derived [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]. To this regard, we remark two important aspects:

1) The derivation of the analytic condition in the general case builds on the solution provided in this paper. The main difference comes from the fact that many scalar quantities become tensors (e.g., L 1 g becomes a two index tensor of type (1, 1); the scalar τ in (3) becomes a three index tensor of type (2, 1)).

2) In contrast with the surprising simplicity and compact architecture of the general solution (see sections 3.3 and 4.2 in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF]), its derivation is very complex. Specifically, it is based on a subtle analogy with the theory of General Relativity and, consequently, on the use of Ricci calculus. In this analogy, the presence of the drift term corresponds to the presence of a time dimension in relativity and the presence of unknown inputs corresponds to the space dimension in relativity (i.e., in the aforementioned analogy, the space dimension of relativity equals the number of unknown inputs). In this sense, the case analyzed in this paper corresponds to the trivial case of a space-time frozen with respect to time and with a single spatial dimension (this is the reason why the derivations of the results provided in this paper (available in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF]) do not require the use of Ricci calculus). The general solution provided in [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF] allows us to automatically derive the observability properties of any nonlinear system driven by disturbances. Note that, in sections 5.4-5.7 of [START_REF] Martinelli | Nonlinear Unknown Input Observability: The General Analytic Solution[END_REF], we revisit the visual-inertial sensor fusion problem. The results include the ones previously obtained [START_REF] Martinelli | Observability Properties and Deterministic Algorithms in Visual-Inertial Structure from Motion[END_REF], [START_REF] Martinelli | Minimalistic sensor design in visual-inertial structure from motion[END_REF]. However, with respect to the former derivation, the new one was obtained automatically and the results are more general.

and ∇τ . 5 )

 5 Use the recursive step of algorithm 2 to build the codistribution Ω m for m ≥ 2, and, for each m, check if ∇τ ∈ Ω m . Denote by m the smallest m such that ∇τ ∈ Ω m . 6) For each m ≥ m , check if Ω m+1 = Ω m and denote by Ω * = Ω m * , where m * is the smallest integer such that m * ≥ m and Ω m * +1 = Ω m * (note that m * ≤ n + 2). 7) If the gradient of a given state component (x

Fig. 1 :

 1 Fig. 1: Local and global frame for the considered problem.The z-axis of the latter is aligned with the direction of the disturbance (assumed to be known and constant in time). The speed V is the vehicle speed with respect to the air, which differs from the ground speed because of the disturbance (w).

Fig. 2 :

 2 Fig. 2: A typical simulated trajectory in absence of disturbance. In blue the true trajectory, in red the trajectory estimated by only integrating the gyroscope and the airspeed measurements and in black the trajectory estimated by the proposed EKF.

Fig. 3 :

 3 Fig. 3: A typical simulated trajectory in presence of disturbance along the z-axis. In blue the true trajectory, in red the trajectory estimated by only integrating the gyroscope and the airspeed measurements and in black the trajectory estimated by the proposed EKF.

Fig. 4 :

 4 Fig. 4: Estimation error on the distance of the two observed point features with and without disturbance (right and left side, respectively).

Fig. 5 :

 5 Fig. 5: As in fig. 4 but for the roll angle.

Fig. 6 :

 6 Fig. 6: As in fig. 4 but for the pitch angle.

Fig. 7 :

 7 Fig. 7: As in fig. 5 but for the yaw angle.

The reader non-familiar with the concept of distribution, as it is used in[START_REF] Isidori | Nonlinear Control Systems[END_REF], should not be afraid by the term distribution and the term codistribution. Very simply speaking, a distribution is a vector space defined on M (our set in R n where the system is defined). In particular, this vector space changes by moving on M . This vector space is in fact the span of a set of vector functions (vector fields) defined on M . A codistribution is the dual of a distribution. Very simply speaking (and this is enough to understand the theory of nonlinear observability) a distribution is generated by a set of column vectors. A codistribution is generated by a set of line vectors. All these vectors are vector functions (i.e., they depend on the point x ∈ M ) and they have the same dimension of x.

The reader is addressed to[START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF] to see the definition of the new output.

The case when L 1 g = 0 for all the outputs is dealt in[START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF].

In other words, it is possible to repeat all the steps of lemmas 1.9.1, 1.9.2 and 1.9.6 in[START_REF] Isidori | Nonlinear Control Systems[END_REF], for algorithm 2 without the term mu i=1 Li φ m-1 ∇h.

The reader is addressed to[START_REF] Martinelli | Nonlinear Unknown Input Observability: The Analytic Solution in the case of a Single Unknown Input[END_REF] where the proofs of all these statements are provided.

Note that the chosen global frame is aligned with he direction of the disturbance (fig.1). Hence, what is unobservable is a rotation around the direction of the disturbance.

*This work was supported by the French National Research Agency ANR 2014 through the project VIMAD A. Martinelli is with INRIA Rhone Alpes, Grenoble,


Anna Dupleix

Tien Van Nguyen

Katja Vahtikari

Mark Hughes
email: mark.hughes@aalto.fi

The anisotropic temperature rise on wood surfaces during adsorption measured by thermal imaging

The aim of the work reported herein was to investigate anisotropy in the temperature rise occurring on the surface of Scots pine (Pinus sylvestris L.) during adsorption. The temperature increase on the transverse, radial and tangential surfaces of small wood specimens during adsorption was measured by thermal imaging. The experiments were conducted in a purpose-built humidity chamber facilitating the accurate control of the internal relative humidity (RH). It was found that the temperature rise on initially oven-dry wood surfaces exposed to a RH of 95% was anisotropic, with a maximum temperature increase of around 4 °C occurring on the transverse surface and around 1 °C on both the radial and tangential surfaces. Additionally, it was observed that the technique reported can be used to accurately measure the temperature change arising from sorption, thus providing reliable information about the behaviour of wood. This information is vital if the effects of the enthalpy change occurring during water vapour sorption are to be correctly considered in building physics models and the potential to reduce energy consumption in buildings fully realised.

Introduction Background and objectives

Wood is a plentiful and renewable material with low embodied energy that is already used extensively as a building construction material. In addition to its excellent structural properties, it can act as moisture and thermal buffer [START_REF] Kraniotis | Thermography measurements and latent heat documentation of Norwegian spruce (Picea abies) exposed to dynamic indoor climate[END_REF], and this ability can lead to reduced energy demand by helping to passively control the internal relative humidity (RH) of a living space [START_REF] Osanyintola | Moisture buffering capacity of hygroscopic building materials: experimental facilities and energy impact[END_REF][START_REF] Holcroft | Heat of sorption and moisture buffering properties of building insulation materials[END_REF]. This would offer the potential to reduce the overall energy consumption of buildings both directly, through lower space heating/cooling needs, and indirectly by reducing the requirements for heating, ventilating and air-conditioning (HVAC) equipment.

When wood adsorbs water vapour from the surrounding air, it is accompanied by a release of heat-the heat of sorption-during the phase change from vapour in the air to water adsorbed in the cell wall [START_REF] Skaar | Thermal conductivity and diffusivity of wood[END_REF]; conversely, when desorption occurs, heat is required. This change in enthalpy can result in a change in surface temperature. [START_REF] Brueckner | Investigating latent heat exchange of untreated panels[END_REF], for example, investigated the surface temperature change during adsorption in untreated spruce wood panels, finding that a surface temperature increase of 2 °C occurred when the ambient relative humidity (RH) rose from 20 to 90%. Recently, the question of whether this enthalpy change could be used to reduce the energy consumption of buildings, has become the subject of speculation and research [START_REF] Kraniotis | Thermography measurements and latent heat documentation of Norwegian spruce (Picea abies) exposed to dynamic indoor climate[END_REF]. Nevertheless, for it to be correctly considered in building physics models, accurate data about the temperature changes in wood arising from sorption are required, and this forms the subject matter of the study reported herein.

Many of the properties of wood are known to be highly anisotropic, and this includes moisture transport. For example, the moisture diffusion coefficient of dry pine is approximately 60 times greater in the longitudinal direction than the transverse direction [START_REF] Siau | Transport processes in wood[END_REF]. This implies that heat may be evolved more rapidly as moisture is adsorbed on the transverse surface of wood, which might also result in a more significant temperature rise on this surface. However, it should be noted that thermal conduction is around 1.5 times greater in the axial direction than in the transverse direction [START_REF] Skaar | Thermal conductivity and diffusivity of wood[END_REF], so this must be factored in as well. In most uses, such as internal wall panelling, the tangential or radial surfaces of wood are exposed to the surroundings; it is just possible that for mediating the internal climate of a room, it may be advantageous to expose the transverse surface. The aim of this study was therefore to investigate the degree of anisotropy in the temperature rise occurring on the transverse, radial and tangential surfaces of wood undergoing sorption.

Physio-chemistry of sorption in wood

According to [START_REF] Kollmann | [END_REF], the total heat of sorption (Q S ) of wood is the overall heat released when a unit mass of wood at a particular moisture content adsorbs moisture up to the fibre saturation point (FSP), expressed as the unit mass of water in wood (kJ/kg water in wood, either dry or not ). Other authors subtract the latent heat of condensation (or vaporisation) of free water at the same temperature, DH cond , taken to be 2443.9 kJ/kg water at 25 °C [START_REF] Marsh | Recommended reference materials for the realization of physicochemical properties[END_REF], to obtain the net heat of sorption, which can be calculated for dry wood (h S ) or for wood at any higher moisture content (q S ). Thus,

• At 0% MC, the heat released (kJ) is given by Q S = h S ? DH cond [START_REF] Simo ´n | Thermodynamic analysis of water vapour sorption behaviour of juvenile and mature wood of Abies alba Mill[END_REF][START_REF] Kraniotis | Thermography measurements and latent heat documentation of Norwegian spruce (Picea abies) exposed to dynamic indoor climate[END_REF]) expressed as the unit mass of water in dry wood (kJ/ kg water in dry wood ); • At higher moisture content (MC), the heat released (kJ) is given by Q S = q S ? DH cond with q S = f(m) where m is the mass of water in wood [START_REF] Skaar | Thermal conductivity and diffusivity of wood[END_REF]) expressed as the unit mass of water in wood (kJ/kg water in wood ).

According to the physio-chemical theory of water adsorption in wood [START_REF] Avramidis | The basics of sorption[END_REF], water molecules are strongly attracted to the various hydrophilic chemical groups, principally hydroxyl (OH) groups, naturally present in the wood cell wall constituents. This attraction leads the water molecules to hydrogen bond (H-bond) to sorption sites in the non-crystalline cell wall regions (Salme ´n 1997; [START_REF] Engelund | A critical discussion of the physics of woodwater interactions[END_REF][START_REF] Perre | Physique et me ´canique du bois (Physics and mechanics of wood)[END_REF]. Since the H-bonds formed between water molecules and the OH groups of the cell wall polymers are always stronger than the H-bonds of the liquid water molecules between them (Perre ´2015), water molecules preferentially 'bind' to the wood rather than simply condense and form a liquid film over the wood substrate. If this was the case, the heat of sorption would equal the latent heat of vaporisation of water, DH cond (Fig. 1). However, being more stable, bound water has a lower energy level than atmospheric water, H bound water \ H atmospheric water [START_REF] Time | Hygroscopic moisture transport in wood[END_REF]. Adsorption-in which atmospheric water vapour molecules become bound to the wood-is therefore exothermic (negative enthalpy), releasing energy in Fig. 1 In the absence of wood, the phase change of water from the gaseous to solid states releases an amount of energy equal to DH cond ? DH freez , whilst in the presence of wood, the phase change of water from vapour to bound water states releases DH cond ? q s (or h S for dry wood). Adapted from [START_REF] Time | Hygroscopic moisture transport in wood[END_REF] and [START_REF] Leuk | Heat of sorption: a comparison between isotherm models and calorimeter measurements of wood pulp[END_REF] the form of heat: DH sorption =

P

hydrogen bonds H bound water À P H atmospheric water \ 0. For [START_REF] Skaar | Thermal conductivity and diffusivity of wood[END_REF], the heat of sorption is the excess binding energy of water molecules 'bound' to the wood solid substrate over that of liquid water molecules between them. It is a direct measure of the amount of water 'bonded' to the cell wall and, according to the physio-chemical theory of water adsorption in wood [START_REF] Avramidis | The basics of sorption[END_REF], it is an indirect measure of the number of sorption sites accessible in wood.

From Skaar's equations, heat of sorption decreases exponentially as wood MC increases, varying from circa 900 kJ/kg water in wood for oven-dry wood to zero when wood reaches FSP. For [START_REF] Avramidis | The basics of sorption[END_REF], this decrease is supported by the physio-chemical adsorption theory in which the binding energy of water to wood lowers when the moisture in wood increases. However, work done by [START_REF] Hofstetter | Moisture uptake in native cellulose-the roles of different hydrogen bonds: a dynamic FT-IR study using deuterium exchange[END_REF] on native cellulose does not confirm this to be the case, calling into question whether the heat of sorption can be explained by binding energies alone [START_REF] Engelund | A critical discussion of the physics of woodwater interactions[END_REF]. Regardless of the underlying mechanisms behind the heat of sorption phenomenon, the sorption of water in wood, resulting from a variation in RH, leads to a change in the temperature which might be beneficially used to affect the overall energy balance of a building.

Measuring surface temperature changes in wood

Measuring the temperature rise on wood surfaces poses several problems. First of all, sorption takes place rapidly, and for this reason when measuring temperature change as a result of a change in RH, the change in RH should be more or less instantaneous. Experimentally this is difficult to achieve, even in dedicated systems such as dynamic vapour sorption (DVS) apparatus, and leads to difficulties in data analysis [START_REF] Hill | The water vapor sorption behavior of flax fibers-Analysis using the parallel exponential kinetics model and determination of the activation energies of sorption[END_REF]. For this reason, a bespoke apparatus was designed and built for this work (see 'Experimental' section for details), which allows for stable internal temperature to be maintained whilst ensuring a rapid, though not instantaneous, change in the internal RH. A second problem relates to the measurement of surface temperature. Traditional ways of measuring temperature rely on contact methods, such as thermocouples. These may either interfere with the sorption process if they are attached to the surface or may not measure the actual surface temperature if they are embedded in the wood material and, as such, they are not ideal for measuring surface temperature changes. Thermography on the other hand has reached a state of development where it has been used to measure the temperature on a variety of materials.

The quantitative use of thermal imaging requires knowledge of the emissivity of the material constituting the target in the spectral range of the infrared camera used to perform the image acquisition [START_REF] Avdelidis | Emissivity considerations in building thermography[END_REF][START_REF] Barreira | Evaluation of building materials using infrared thermography[END_REF]. This physical quantity characterises the ability of a material to emit thermal radiation at particular wavelengths. This explains why literature values for wood emissivity have only been gathered from 7.5 to 13 lm (769 to 1333 cm -1 ) which is within the spectral range of the infrared FLIR e60 camera used in this present work (FLIR 2013). Few literature values carefully specify the wood species, temperature and moisture content used during emissivity measurements (Table 1).

However, the overall trends in the influence of temperature and moisture content can be summarised as follows: (1) dry wood behaves as a dielectric material; its emissivity decreases with temperature [START_REF] Lopez | Determination of the emissivity of wood for inspection by infrared thermography[END_REF][START_REF] Gaussorgues | La thermographie infrarouge: principes-technologie-applications (Infrared thermography: principles-technology-applications[END_REF]), ( 2) the emissivity of wood increases with its moisture content up to FSP at which point it reaches the emissivity of water (Kollmann and Co ˆte ´1968). Work done by [START_REF] Dupleix | Mid-infrared absorption properties of green wood[END_REF] confirms that, in the spectral range of interest 7.5-13 lm, wood emissivity does not vary with moisture content above FSP.

Materials and methods

In the design of the experiments reported herein, special care was taken to maximise the precisions of the measurements (camera settings) and minimise uncertainties (temperature and humidity conditions inside the chamber, measurement protocol) in order to obtain as accurate values as possible in these delicate measurements.

Experimental set-up

Specimen preparation

Boards of Scots pine (Pinus sylvestris L.) obtained from a tree grown in southern Finland were used throughout. The specimens were cut consecutively from the same board, and they contained both heartwood and sapwood. The following denomination for the specimens is used hereinafter: 'radial' specimens refer to specimens in which the exposed surface is radial, i.e. they are cut in the longitudinal-tangential plane (L-T in Fig. 2b), 'tangential' specimens are specimens in which the exposed surface is tangential, i.e. the longitudinal-radial plane (L-R in Fig. 2b) and in 'transverse' specimens, the transverse surface is exposed, i.e. the tangential-radial plane (T-R in Fig. 2b). In order to limit structural differences, the 10 replicates tested were five pairs of 'mirror' samples. However, due to limitations on the amount of material, the transverse samples were not mirror samples but came from locations close to one another. The sample size of 50 9 50 9 25 mm 3

COV (%) 2.7
Values in bold typeface derive from calculation (length 9 width 9 height) was chosen so as to be as large as possible in order to avoid any edge effects affecting the region of interest of the camera display. This sample size also represented a compromise between the space available inside the humidity chamber (see 'Humidity chamber' section), the focal length of the IR camera (195 mm) and the desire to maintain a constant surface-to-volume ratio irrespective of the grain orientation.

Humidity chamber

A purpose-built humidity chamber (Fig. 2a) was constructed. It featured two fans, a mist maker and a water tank enabling a homogeneous and constant RH of up to 95% to be reached. The fans and mist maker were controlled by an Arduino Uno microcontroller board, which automatically adjusted the conditions until the set point value of the data logger placed inside the chamber was reached. All devices were controlled, and their measurements recorded as a function of time, through a LabVIEW Ò interface. One fan was a basic model which could only be switched on or off according to the needs of the mist maker to reach the humidity set point value; the other fan, from Antec Tricool, revolved at constant speed (1600 rpm) during the experiments in order to ensure homogeneous conditions within the chamber. The dimensions of the chamber were 300 9 300 mm 2 (width 9 length), and the height was 600 mm with the tank and 400 mm without the tank. The front door was closed by magnets and, in addition, a thick transparent plastic film was used to cover the doorway to minimise disruption of the interior conditions when handling the samples placed inside. The humidity chamber was equipped with one Hygrofox Mini Ò data logger placed at the sample level-to record RH and temperature (accuracy: ± 3% and ± 0.4 °C, respectively). The plastic cover of the data logger was removed to ensure better contact of the inner sensors with their surroundings. The chamber itself was located in a conditioning room stabilised at 65% and 23 °C.

Camera

A FLIR infrared (IR) camera e60 (uncooled microbolometer detector) was chosen for its compact format so that it would fit in the humidity chamber. The position of the thermal camera was adjusted to minimise changes to the camera angle of view and the distance between the camera and the measured surface. The distance between the camera and the specimen was chosen in order to minimise the (moist) air layer in between so as to reduce possible disturbances. The angle between the camera angle of view and the horizontal plane was set at around 25°in order to reduce IR reflection from the camera lens and respect the limiting angle (45°/50°) advised by the camera manufacturers FLIR, above which the emissivity of materials such as wood starts to change (personal communication). The pertinent camera characteristics are summarised in Table 2, and further detailed information can be found in FLIR e-series datasheet (FLIR 2013).

Measurement protocol

Specimens

The specimens were covered on five faces with an Advance Ò aluminium tape previously tested for its ability to prevent any uptake of water vapour. Only the top face of the specimen was freely exposed to humidity when placed inside the chamber.

Time

The humidity chamber set point of 95% RH is reached within a couple of minutes, but in the experiments, the chamber was stabilised for at least 1 h prior to commencing the measurements in order to minimise effects from any artificial temperature sources such as the fans and the camera motors. Allowing for the temperature to equilibrate was found to be extremely important since the surface temperature rise resulting from heat of sorption was of the order of a couple of degree Celsius.

Reflective ambient temperature (RAT)

RAT takes into account the parasitic radiation from the surroundings reflected by the measured surface. It is determined by recording the mean temperature of a mirror reflecting all incident radiation when the IR camera emissivity is set to 1 (a mirror surface was obtained by covering the top face of a specimen with crumpled Advance Ò aluminium tape).

Surface temperature rise

The specimens were dried to 0% MC in a fan-assisted oven at 103 °C for 24 h. Upon removal from the oven, aluminium tape was quickly applied to 5 surfaces and a strip of black electrical tape (reference emissivity) applied to the exposed surface. The specimens were then immediately transferred to a desiccator (0% RH) where they were cooled to room temperature over silica gel for an additional 24 h. The open surface of the specimens was then covered with a 'lid' made out of aluminium foil and placed inside the humidity chamber about 1 h before the experiment began to ensure that its temperature equilibrated with that of the humidity chamber. The sample was placed inside the humidity chamber so that it was in the field of view of the IR camera. The camera recording began when the aluminium foil 'lid' was removed. This protocol ensures that the surface temperature changes on the wood specimen are recorded immediately the (dry) wood surface comes into contact with the high RH of the chamber; since the temperature and RH remain constant inside the chamber (Fig. 3b). The temperature and RH conditions inside the humidity chamber were considered to be ostensibly constant at 22.7 ± 0.4 °C and 95% ± 3% according to the records from the data logger plotted in Fig. 3b for one experiment, but being representative of all others. Throughout the test there was no 'lag' during which the conditions stabilise. This is important when it comes to detecting the rate of heat evolution and temperature rise and to the authors' knowledge this level of control has not been employed previously.

Results

Wood density

The oven-dry densities of the wood samples, calculated from the weight and dimensions of the specimens, are presented in Table 3, with the average and standard deviations shown. The low variation in the density values, given by the low COV, highlights the homogeneity of the sample densities, indicating that the mean value of 430 kg/m 3 is a representative value for all the specimens.

Surface temperature rise

A typical thermal image captured during an experiment is shown in Fig. 3a. The image shows the transverse surface of a specimen, and anatomical features, such as the growth rings, are clearly visible on the right-hand side of the image; the distinction suggests variations in temperature throughout the wood material. The right-hand side of the image shows the surface covered with black PVC tape which was used later as an impermeable coating to serve as a reference surface to illustrate the stability of the temperature conditions inside the humidity chamber during the experiments (see Fig. 4c).

Figure 4 shows that the surface temperature-measured as the average temperature over the region of interest, including early wood and late wood-increases as a function of time, DT(t), when wood, initially in the oven-dried state (i.e. 0% MC) is placed inside the humidity chamber and is suddenly exposed to 95% RH. Temperature rises on the radial, tangential and transverse surfaces are shown in Fig. 4a-c, respectively. The surface temperature increase of wood is given by DT(t) = T(t)-T i (0), where T(t) is the surface temperature recorded by the camera and T i (0) is the initial temperature of the sample surface recorded by the camera at t = 0. The emissivity values for the wood sample are taken from the literature (Table 2) and set to 0.95 (see 'Discussion' section).

The black curve is the average of the 10 replicates tested, whilst the grey curve above is the average ? the standard deviation and the grey curve below is the average-the standard deviation. Table 4 summarises, for each of the three surfaces measured, the maximum temperature increase, DT max , and the time needed to attain this value, t max .

Discussion

Anisotropic surface temperature change

The form of the DT(t) curves (Fig. 4) was found to be consistent for each of the ten replicates tested on each of the three wood surfaces (i.e. radial, tangential and transverse). However, the results from the transverse surfaces were noted to differ significantly from those of the radial and tangential surfaces (Fig. 4, Table 4). As may be seen, the surface temperature on the radial and tangential surfaces rises more rapidly (mean t max = 48 s), but reaches a lower maximum temperature (mean DT max of 1.4-1.5 °C) than the transverse surface for which the rise is higher (mean DT max of 3.6 °C) and takes longer to reach its peak value (mean t max = 282 s). The results of [START_REF] Kraniotis | Thermography measurements and latent heat documentation of Norwegian spruce (Picea abies) exposed to dynamic indoor climate[END_REF] on the tangential surface of spruce heartwood are of the same order of magnitude as the current findings, albeit the behaviour was observed to be closer to that seen on the transverse surfaces of the specimens tested in this present work (t max = 500 s/DT max = 2.1 °C). This difference is most probably attributable to the different experimental procedure adopted by Kraniotis and co-workers, who measured the surface temperature change on wood that had a nonzero initial moisture content and under conditions of slowly increasing relative humidity, changing from 25 to 90% over the course of 12 h. Further, the work of [START_REF] Kraniotis | Thermography measurements and latent heat documentation of Norwegian spruce (Picea abies) exposed to dynamic indoor climate[END_REF] was on spruce heartwood, which would doubtlessly behave differently to the Scots pine used in this study that contained both sapwood and heartwood, whilst the specimens also differed in size which would have undoubtedly affected the mass and heat transfer characteristics. These findings suggest some interesting possibilities for using wood to help passively lower the heating (or cooling) energy used in buildings. Typically solid wood panelling, for instance, is installed with the tangential surface (or a combination of tangential and radial surfaces) exposed to the room in which it is situated. The wood is also generally coated, which may prevent the exchange of moisture between wood and the surroundings. These results suggest that by exposing the 'end grain' of the wood, i.e. the transverse surface, together with the use of coatings that do not hinder the exchange of moisture [START_REF] Lozhechnikova | Towards energy efficiency through an optimized use of wood: the development of natural hydrophobic coatings that retain moisturebuffering ability[END_REF], the potential for utilising the heat of sorption phenomenon could be maximised.

Influence of wood structure

For the interpretation of the results, the following assumption is made; the main direction of propagation of moisture and heat is normal to the exposed surface.

Moisture and heat transfer are thus mainly in the longitudinal direction in specimens with the transverse surface exposed (and naturally either in the radial or tangential directions in the others) (Fig. 2b). Anatomical knowledge of the structural shape and arrangement of cells in wood [START_REF] Comstock | Directional permeability of softwoods[END_REF]) helps us to understand the higher temperature rise on the transverse surface. In the longitudinal direction, water vapour can flow more readily and for greater distances into the wood structure in the gaseous phase before adsorbing onto the surfaces of the cell wall polymers. Softwood tracheids can be several millimetres in length, so that even when cut across, the lumens may be more than a millimetre in length, allowing for rapid adsorption on to a relatively large surface area. In contrast, there is more resistance to moisture movement in the radial and tangential directions, since flow is impeded by the cell wall (Perre ´2015) and the necessity for vapour in the cavities to flow through the connecting pits. Moisture diffusion coefficients are consequently higher in the longitudinal direction (3.10 -5 cm 2 /s at 25 °C) than in the transverse (5.10 -7 cm 2 /s at 25 °C) (values given by Siau 1984 for oven-dry pine of similar 500 kg/m 3 density to that used in the present study, Table 3). Therefore, the more significant temperature rise on the transverse surface is a logical consequence of the fact that water is adsorbed, flows and forms H-bonds more easily in the longitudinal direction where the sorption sites are anatomically more accessible in wood.

As observed by [START_REF] Comstock | Directional permeability of softwoods[END_REF], the higher variability of the results on the radial and tangential surfaces might be due to practical sampling matters, in that no perfectly radial or tangential surfaces can be achieved on cubic specimens. In such cases, the different grain angles can vary from one sample to another adding heterogeneous features to the material structure, which is not the case for transverse samples where the slope of the grain can be truly longitudinal, especially for such small surface dimensions (50 9 50 mm).

Assessment of experimental errors

The accurate and reliable measurement of surface temperature, under known conditions of relative humidity and temperature and defined material parameters, such as moisture content and exposed anatomical surface, is essential if the potential of the latent heat of sorption is to be fully exploited in reducing energy demands in buildings and to be incorporated into building physics models.

Stability of conditions within humidity chamber

The surface temperature recorded on the reference sample, in which adsorption was prevented by a covering of impermeable tape (Fig. 4c), provides confirmation that the internal conditions within the humidity chamber remained constant throughout the experiment. This strongly suggests that the results obtained on the wood sample surfaces represented real temperature changes, despite the limitations of the thermographic method, which are discussed below.

Influence of camera accuracy

By using relative values (DT), the accuracy of the camera (± 2 °C given by the manufacturers, FLIR) is considered here to be offset, given the stability of the measurements obtained from one point to another.

Influence of emissivity

The choice of wood emissivity was carefully considered during the experiments in order to ensure that the temperature rise recorded was not an artefact due to any change of emissivity with the moisture content in wood during the experiments. The median value from the literature for softwood species in the corresponding wavelength region (emissivity = 0.95, [START_REF] Dupleix | Mid-infrared absorption properties of green wood[END_REF]) was chosen given the fact that wood emissivity increases with the amount of moisture in wood, and the assumption is made that wood material near the surface might equilibrate almost immediately with the surrounding RH (90%) to reach a moisture content near to the saturated state (neglecting the surface resistance, [START_REF] Siau | Transport processes in wood[END_REF]. The latter assumption also allows us to neglect the influence of changes in the wood moisture on emissivity in the same way as the influence of changes in the surface temperature of wood on emissivity is neglected because of the small increase recorded (a couple of degrees Celsius). The validity of these assumptions is confirmed in Fig. 5, which shows the sensitivity of wood surface temperature to changes in emissivity value. As may be observed a change in emissivity of ± 0.02, a variation, which might be expected to occur due to changes in moisture content, for example, results in only very slight changes (less than 0.01 °C) in the measured temperature.

Conclusion

The adsorption of water in oven-dried wood generates a surface temperature increase of several degrees Celsius, which is recordable by thermal imaging. This phenomenon-the heat of sorption-arises because of the phase change of water from vapour to 'bound' water at the sorption sites in the wood cell wall. The experimental results presented here highlight the anisotropy of the surface temperature rise on the surface of wood due to heat of sorption. The transverse surface of Scots pine wood can undergo a temperature increase of between 3 and 4 °C when exposed to high humidity (95%) conditions, whereas the temperature rise on the radial and tangential surfaces was more modest at around 1 to 2 °C, though the maximum temperature was reached rather quickly, and quicker than on the transverse surface. These results validate the ability of the experimental protocol developed here to provide source data useful for calculating the heat of sorption (in terms of unit of energy per mass of wood released) and modelling the energetic impact of wooden surfaces in building interiors. Additional work is ongoing in this respect. Further work studying the heat of sorption on wood coated with different materials would be beneficial in the practical use of this phenomenon in interior design, to benefit the future development of the application of wood in buildings.

Fig. 2 a

 2 Fig. 2 a Experimental set-up showing (b) the three directions of wood (L for longitudinal, R for radial and T for tangential) and the corresponding surfaces (radial L-R, tangential L-T and transverse T-R)

Fig. 3 a

 3 Fig. 3 a Snapshot of a radial surface taken with the IR camera as recorded during the experiments (right side covered with Advance Ò black PVC tape and left side, uncovered) and b conditions of temperature (°C) and RH (%) inside the humidity chamber during the experiments recorded with the data logger

Fig. 4

 4 Fig. 4 DT(t) recorded for (a) radial, (b) tangential and (c) transverse samples for the first 10 min of experiments and T(t) for longer time period (d). A reference sample fully covered with Advance Ò PVC black tape is also presented in (c)

Table 1

 1 Comparison of experimental emissivity values reported for pine and other wood species (when known) in the literature in the spectral range 7.5-13 lm

	Spectral range	Temp.	MC (%) Wood species	Emissivity Authors	Median
	(lm)	(°C)					values
	7.5-13	22	N/A	Pine	0.92	Lopez et al. (2013)	0.92
	7.5-13	20-25 20-22	Spruce, Douglas-fir 0.93-0.97 Dupleix et al. (2013) 0.95
	7.5-13	20-25 64 68	Beech, birch	0.95-0.98 Dupleix et al. (2013) 0.97
	Full spectrum 20	N/A	N/A	0.90	Gaussorgues (1999)	0.90
	Mean						0.93

Table 3

 3 Mean, standard deviation and COV values for the samples densities (10 parallel samples for each direction)

		Transverse	Radial	Tangential	All samples
	Mean (kg/m 3 )	428.1	432.2	429.7	430.0
	SD	12.4	22.6	16.4	17.4
	COV (%)	2.9	5.2	3.8	4.0

Table 4

 4 Mean, standard deviation and COV values for DT max and t max (10 replicates for each direction)

		DT max (°C)			t max (s)
		Transverse	Radial	Tangential	Transverse	Radial	Tangential
	Mean	3.6	1.4	1.5	282	48	48
	SD	0.2	0.3	0.3	92	4	6
	COV (%)	4.7	19.5	22.0	33	9

Acknowledgements This work was supported by the Wood Life project [2013][2014][2015][2016][2017] funded by Aalto University within the framework of the Aalto Energy Efficiency Programme (AEF). Fig. 5 Effect of emissivity on temperature simulated with FLIR Tools software. Example of a transverse sample