Cyril Roussillon

Aurélien Gonzalez

Joan Solà

Jean-Marie Codol

Nicolas Mansard

Simon Lacroix

Michel Devy

RT-SLAM: a generic and real-time visual SLAM implementation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Motivation

Progresses in image processing, the growth of available computing power and the proposition of approaches to deal with bearings-only observations have made visual SLAM very popular, particularly since the demonstration of a real-time implementation by Davison in 2003 [START_REF] Davison | Real-time simultaneous localisation and mapping with a single camera[END_REF]. The Extended Kalman Filter (EKF) is widely used to solve the SLAM estimation problem, but it has recently been challenged by global optimization methods that have shown superior precision for large scale SLAM. Yet EKF still has the advantage of simplicity and faster processing for problems of limited size [START_REF] Strasdat | Real-time monocular slam: Why filter?[END_REF], and can be combined with global optimization methods [START_REF] Estrada | Hierarchical slam: Real-time accurate mapping of large environments[END_REF] to take the best of both worlds.

Monocular EKF-SLAM reached maturity in 2006 with solutions for initializing landmarks [START_REF] Kwok | An efficient multiple hypothesis filter for bearingonly slam[END_REF] [12] [START_REF] Montiel | Unified inverse depth parametrization for monocular slam[END_REF]. Various methods for landmark parametrization have been analyzed [START_REF] Sola | Consistency of the monocular ekf-slam algorithm for three different landmark parametrizations[END_REF], and the literature now abounds with contributions to the problem.

This article presents RT-SLAM 1 , a software framework aimed at fulfilling two essential requirements. The first one is the need for a generic, efficient and flexible development tool, that allows to easily develop, evaluate and test various approaches or improvements. The second one is the need for practical solutions for live experiments on robots, for which localization and mapping require real-time execution and robustness. RT-SLAM is a C++ implementation based on Extended Kalman Filter (EKF) that allows to easily change robot, sensor and landmark models. It runs up to 60 Hz with 640 × 480 images, withstanding highly dynamic motions, which is required for instance on humanoid or high speed all terrain robots. RT-SLAM is available as open source software at http://rtslam.openrobots.org.

Section 2 details the architecture of RT-SLAM and section 3 presents some of the techniques currently integrated within, to define an efficient inertial/visual SLAM approach. Section 4 analyzes results obtained with a hand-held system assessed thanks to ground truth, and section 5 concludes the article by a presentation of prospects.

Overall Architecture

Fig. 1 presents the main objects defined in RT-SLAM. They encompass the basic concepts of a SLAM solution: the world or environment contains maps; maps contain robots and landmarks; robots have sensors; sensors make observations of landmarks. Each of these objects is abstract and can have different implementations. They can also contain other objects that may themselves be generic. Map. The maps contain an optimization or estimation engine: for now RT-SLAM uses a standard formulation of EKF-SLAM. Since this solution is very well documented in the literature [START_REF] Davison | Monoslam: Real-time single camera slam[END_REF], it is not detailed in depth here. Indirect indexing within Boost's ublas C++ structures is intensively used to exploit the sparsity of the problem and the symmetry of the covariances matrices.

Robot. Robots can be of different type according to the way their state is represented and their prediction model. The latter can be either a simple kinematic model (constant velocity, constant acceleration, . . . ) or a proprioceptive sensor (odometric, inertial, . . . ), as illustrated section 3.3. The proprioceptive sensor is an example of generic object contained in robot objects, as different hardware can provide the same function.

Sensor. Similarly to robots, sensors can also have different models (perspective camera, panoramic catadioptric camera, laser, . . . ), and contain a generic exteroceptive sensor hardware object (firewire camera, USB camera, . . . ). In addition, as sensors belong to the map, their state can be estimated: this opens the possibility for estimating other parameters such as extrinsic calibration, time delays, biases and gain errors, and the like.

Landmark. Landmarks can be of different type (points, lines, planes, . . . ), and each type can have different state parametrization (Euclidean point, inverse depth point, . . . ). Moreover the parametrization of a landmark can change over time, as explained section 3.2. A landmark also contains a descriptor used for data association, which is a dual description to the state representation.

As shown Fig. 1(a), it is worth noticing that landmark objects are common to the different sensors, all of them being able to observe the same landmark (provided they have compatible descriptors for this landmark of course). This allows to greatly improve the observability of landmarks compared to a system where the sensors are strictly independent. In the particular case of two cameras for instance, landmarks can be used even if they are only visible from one camera or if they are too far away for a stereovision process to observe their depth (this process was introduced in [START_REF] Sola | Bicamslam: Two times mono is more than stereo[END_REF] as BiCam SLAM).

Observation. In RT-SLAM, the notion of observation plays a predominant role. An observation is a real object containing both methods and data. One observation is instantiated for every sensor-landmark pair, regardless of the sensor having actually observed the landmark or not, and has the same lifetime as the associated landmark. The methods it contains are the conventional projection and back-projection models (that depend on the associated sensor and landmark models), while the stored data consist of results and intermediary variables such as Jacobian matrices, predicted and observed appearances, innovations, event counters and others, that allow to greatly simplify and optimize computations.

Managers. In order to achieve full genericity wrt landmark types, in particular to allow the concurrent use of different landmark types for one sensor, two different manager objects are added: data manager and map manager. Their implementations define a given management strategy, while their instantiations are dedicated to a certain landmark type. The data manager processes the sensors raw data, providing observations of the external landmarks. For this purpose, it exploits some raw data processors (for feature detection and tracking), and decides which observations are to be corrected and in which order, according to the quantity of information they bring and their quality. For example it can apply an active search strategy and try to eliminate outliers as described in section 3.1. The map manager keeps the map clean, with relevant information, and at a manageable size, by removing landmarks according to their quality and the given policy (e.g. visual odometry where landmarks are discarded once they are not observed, or multimap slam where maps are "closed" according to given criteria). These managers communicate together: for example, the data manager may ask the map manager if there is enough space in the map to start a new initialization process, and to allocate the appropriate space for the new landmark.

3 Inertial/visual SLAM within RT-SLAM

Active search and one-point RANSAC

The strategy currently implemented in RT-SLAM's data manager to deal with observations is an astute combination of active search [START_REF] Davison | Monoslam: Real-time single camera slam[END_REF] and outliers rejection using one-point RANSAC [START_REF] Civera | 1-point ransac for ekfbased structure from motion[END_REF].

The goal of active search is to minimize the quantity of raw data processing by constraining the search in the area where the landmarks are likely to be found. Observations outside of this 3σ observation uncertainty ellipse would be anyway considered incompatible with the filter and ignored by the gating process. In addition active search gives the possibility to decide anytime to stop matching and updating landmarks with the current available data, thus enabling hard real-time constraints. We extended the active search strategy to landmark initialization: each sensor strives to maintain features in its whole field of view using a randomly moving grid of fixed size, and feature detection is limited to empty cells of the grid. Furthermore the good repartition of features in the field of view ensures a better observability of the motions.

Outliers can come from matching errors in raw data or mobile objects in the scene. Gating is not always discriminative enough to eliminate them, particularly right after the prediction step when the observation uncertainty ellipses can be quite large -unfortunately at this time the filter is very sensitive to faulty corrections because it can mistakenly make all the following observations incompatible. To prevent faulty observations, outliers are rejected using a one-point RANSAC process. It is a modification of RANSAC, that uses the Kalman filter to obtain a whole model with less points than otherwise needed, and provides a set of strongly compatible observations that are then readily corrected. Contrary to [START_REF] Civera | 1-point ransac for ekfbased structure from motion[END_REF] where data association is assumed given when applying the algorithm, we do the data association along with the one-point RANSAC process: this allows to look for features in the very small strongly compatible area rather than the whole observation uncertainty ellipse, and to save additional time for raw data processing.

Landmark parametrizations and reparametrization

In order to solve the problem of adding to the EKF a point with unknown distance and whose uncertainty cannot be represented by a Gaussian distribution, point landmarks parametrizations and initialization strategies for monocular EKF-SLAM have been well studied [START_REF] Davison | Real-time simultaneous localisation and mapping with a single camera[END_REF] [8] [START_REF] Montiel | Unified inverse depth parametrization for monocular slam[END_REF]. The solutions now widely accepted are undelayed initialization techniques with inverse depth parametrization. Anchored Homogeneous Point [START_REF] Sola | Consistency of the monocular ekf-slam algorithm for three different landmark parametrizations[END_REF] parametrization is currently used in RT-SLAM.

The drawback of inverse depth parametrizations is that they describe a landmark by at least 6 variables in the stochastic map, compared to only 3 for an Euclidean point (x y z)

T . Memory and temporal complexity being quadratic with the map size for EKF, there is a factor of 4 to save in time and memory by reparametrizing landmarks that have converged enough [START_REF] Civera | Inverse Depth to Depth Conversion for Monocular SLAM[END_REF]. The map manager uses the linearity criterion proposed in [START_REF] Montiel | Unified inverse depth parametrization for monocular slam[END_REF] to control this process.

Motion prediction

The easiest solution for EKF-SLAM is to use a robot kinematic model to do the filter prediction, such as a constant velocity model: R = (p q v w) T where p and q are respectively the position and quaternion orientation of the robot, and v and w are its linear and angular velocities.

The advantage of such a model is that it does not require complicated hardware setup (only a simple camera), but its strong limitation is that the scale factor is not observable. A second camera with a known baseline can provide a proper scale factor, but one can also use a proprioceptive sensor for the prediction step. Furthermore it usually provides a far better prediction with smaller uncertainties than a simple kinematic model, which brings several benefits: 1. search areas for matching are smaller, so processing is faster, 2. linearization points are more accurate, so SLAM precision is better, or one can reduce the framerate to decrease CPU load for equivalent quality, 3. it allows to withstand high motion dynamics.

In the case of an Inertial Measurement Unit (IMU), the robot state is then:

R = (p q v a b w b g) T
where a b and w b are the accelerometers and gyrometers biases, and g the 3D gravity vector. Indeed it is better for linearity reasons to estimate the direction of g rather than constraining the robot orientation to enforce g to be vertical.

A special care has to be taken for the conversion of the noise from continous time (provided by the manufacturer in the sensor's datasheet) to discrete time. As the perturbation is continuous white noise, the variance of the discrete noise grows linearly with the integration period.

Image processing

Point extraction. Point extraction is based on Harris detector with several optimizations. Some of them are approximations: a minimal derivative mask [-1, 0, 1] is used, as well as a square and constant convolution mask, in order to minimize operations. This allows the use of integral images [START_REF] Viola | Rapid object detection using a boosted cascade of simple features[END_REF] to efficiently compute the convolutions. Additional optimizations are related to active search (section 3.1): only one new feature is searched in a small region of interest, which eliminates the costly steps of thresholding and sub-maxima suppression.

Point matching. Point matching is based on Zero-mean Normalized Cross Correlation (ZNCC), also with several optimizations. Integral images are used to compute means and variances, and a hierarchical search is made (two searches at half and full resolution are sufficient). We also implemented bounded partial correlation [START_REF] Di Stefano | Zncc-based template matching using bounded partial correlation[END_REF] in order to interrupt the correlation score computation when there is no more hope to obtain a better score than the threshold or the best one up to now. To be robust to viewpoint changes and to track landmarks longer, tracking is made by comparing the initial appearance of the landmark with its current predicted appearance [START_REF] Davison | Monoslam: Real-time single camera slam[END_REF].

Results

Fig. 2 shows the hardware setup that has been used for the experiments. It is composed of a firewire camera rigidly tied to an IMU, on which four motion capture markers used to measure the ground truth are fixed. Two different sequences are used, referred to as low dynamic and high dynamic sequences. Both were acquired indoor with artificial lights only, with an image framerate of 50 Hz synchronized to the inertial data rate of 100 Hz. The motion capture markers are localized with a precision of approximately 1 mm, so the ground truth has a precision of σ xyz = 1 mm in position and σ wpr = 0.57 • in angle (baseline of 20 cm).

Movies illustrating a run of every experiment are provided at: http://rtslam.openrobots.org/Material/ICVS2011.

Constant velocity model

The inertial data is here not used, and the prediction is made according to a constant velocity model. Fig. 3 presents the estimated trajectory and the ground truth for the low dynamic sequence, and Fig. 4 shows the absolute errors for the same run. Besides the global scale factor which is manually corrected, the camera trajectory is properly estimated. The movie shows that loops are correctly closed, even after a full revolution.

The position error raises up to 10 cm after quick rotations of the camera: this is due to a slight drift of the scale factor caused by the integration of newer landmarks (when closing loops, the scale factor is restimated closer to its initial value). Also, uncertainty ellipses remain relatively large throughout the sequence: these issues are solved by the use of an IMU to predict the motions.

4.2 Inertial/visual SLAM Fig. 5 shows the trajectory estimated with the high dynamic sequence, and Fig. 6 and7 show the behavior of inertial SLAM. Here, all of the 6 degrees of freedom are successively excited, then y and yaw are excited with extreme dynamics: the yaw angular velocity goes up to 400 deg/s, the rate of change of angular velocity exceeds 3,000 deg/s 2 , and accelerations reach 3 g. It is interesting to note that the time when SLAM diverges (around t = 65 s) corresponds to motions for which the angular velocity exceeds the limit of the IMU (300 deg/s) and where its output is saturated.

The IMU now allows to observe the scale factor, and at the same time reduces the observation uncertainty ellipses and thus eases the active search procedure. Conversely, the divergence of the SLAM process at the end of the sequence illustrates what happens when vision stops stabilizing the IMU: it quickly diverges because the biases and the gravity cannot be estimated anymore.

Outlook

We have presented a complete SLAM architecture whose genericity and performance make it both a useful experimentation tool and efficient solution for robot localization. The following extensions are currently being developed: using a second camera to improve landmarks observability, a multimap approach to cope with large scale trajectories and multirobot collaborative localization, and the use of line landmarks complementarily to points to provide more meaningful maps and to ease map matching for loop closure. Fig. 6. Errors of inertial/visual SLAM estimated over 100 runs on the dynamic sequence, as a function of time (in seconds) -the difference between each run is due to the randomized landmark selection, see section 3.1. All the runs remain in the neighborhood of the 99% confidence area. The angular ground truth uncertainty is not precisely known: its theoretical value is both overestimated and predominant over SLAM precision in such a small area. Fig. 7. Inertial/visual SLAM NEES [START_REF] Bailey | Consistency of the ekf-slam algorithm[END_REF] over 100 runs on the dynamic sequence. The average NEES is quickly out of the corresponding 99% confidence area, but as the 100 runs were made on the same sequence they are not independant and the average NEES should rather be compared to the simple NEES confidence area. The filter appears to be very conservative with angles but this is due to the overestimated ground truth uncertainty as explained in Fig. 6.

The architecture of RT-SLAM allows to easily integrate such developments, and also to consider additional motion sensors: to increase the robustness of the system, it is indeed essential to consider the various sensors usually found on board a robot (odometry, gyrometers, GPS). Eventually, it would be interesting to make RT-SLAM completely generic wrt the estimation engine as well, in order to be able to use global optimization techniques in place of filtering.

  (a) Main objects in a SLAM context. Different robots Rob have several sensors Sen, providing observations Obs of landmarks Lmk. States of robots, sensors and landmarks are stored in the stochastic map Map. There is one observation per sensor-landmark pair. (b) Objects hierarchy in RT-SLAM. Each individual map M in the world W contains robots R and landmarks L. A robot has sensors S, and an observation O is created for every pair of sensor and landmark. In order to allow full genericity, map managers MM and data managers DM are introduced.

Fig. 1 .

 1 Fig. 1. The main objects in RT-SLAM

Fig. 2 .

 2 Fig. 2. The experimental setup composed of a Flea2 camera and an XSens MTi IMU, and screen captures of the 3D and 2D display of RT-SLAM.

Fig. 3 .

 3 Fig. 3. Illustration of low dynamic trajectory, constant velocity SLAM (with scale factor of 2.05 manually corrected). Estimated trajectory parameters and ground truth, as a function of time (in seconds).

Fig. 4 .

 4 Fig.4. Errors of the estimated parameters of the trajectory shown Fig.3. The 99% confidence area corresponds to 2.57σ bounds. The SLAM 99% confidence area, that does not include ground truth uncertainty, is also shown.

Fig. 5 .

 5 Fig. 5. Illustration of high dynamic trajectory, inertial/visual SLAM. Estimated trajectory parameters and ground truth, as a function of time (in seconds).

RT-SLAM stands for "Real-Time SLAM"


Jaemin Lee

Nicolas Mansard
email: nicolas.mansard@laas.fr

Jaeheung Park
email: park73@snu.ac.kr

Intermediate Desired Value Approach for Continuous Transition among Multiple Tasks of Robots

As the capability of robots is getting improved, more various tasks are expected to be performed by the robots. Complex operation of the robots can be composed of many different tasks. These tasks are executed sequentially, simultaneously, or in a combined way of both. This paper discusses the transition issue among multiple tasks on how the transition can be effectively and smoothly achieved. The proposed approach is to compose intermediate desired values to smooth the transitions rather than to modify control laws. The approach can be practically used on robotic systems without modification on their specific control algorithms. In this paper, multi-points control and joint limit avoidance are performed as applications of the proposed approach.

I. INTRODUCTION

The robots are mostly programmed to execute a specific pre-defined set of tasks in a well defined environment. However, as more complex robots are designed to operate in dynamic environments, various tasks or set of tasks should be executed sequentially or simultaneously [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF], [START_REF] Baerlocher | Task-Priority Formulations for the Kinematic Control of Highly Redundant Articulated Structures[END_REF]. Dealing with these situations, it is important to design strategies for transition among various tasks [START_REF] Dragonir | Dynamic Task-priority allocation for kinematically redundant robotic mechanisms[END_REF], [START_REF] Baerlocher | An Inverse Kinematic Architecture Enforcing an Arbitrary Number of strict priority levels[END_REF], [START_REF] Siciliano | A General Framework for Managing Multiple Tasks in Highly Redundant Robotic Systems[END_REF].

Considering general task transitions, this paper is focused on the situation when a task is added or deleted. Once established, this approach can be used for many task transition situations. Such cases can occur also in traditional applications such as dealing with joint limit where joint limit avoidance can be defined as a task in specific regions [START_REF] Mansard | A Unified approach to integrate unilateral constraints in the stack of tasks[END_REF]. Feature tracking in visual servoing is another application where some of the features are lost or added during operation [START_REF] Mansard | A Unified approach to integrate unilateral constraints in the stack of tasks[END_REF], [START_REF] Mansard | Task sequencing for sensor-based control[END_REF], [START_REF] Mansard | Continuity of Varying-Feature-Set Control laws[END_REF]. Also, Singularity avoidance in manipulation [START_REF] Chiaverini | Singularity-Robust Task-Priority Redundancy Resolution for Real-Time Kinematic Control of Robot Manipulators[END_REF] and switching behaviors of humanoid robots would be also one of the cases [START_REF] Khatib | Whole-Body Dynamic Behavior and Control of Human-like Robots[END_REF], [START_REF] Sentis | Control of Free-Floating Humanoid Robots Through Task Prioritization[END_REF], [START_REF] Yoshida | Humanoid Motion Planning for Dynamic Tasks[END_REF].

Without proper consideration of transition, there will be discrete and abrupt changes at the input of the robot, which is typically either desired joint velocity or joint torque. The goal of the paper is to provide an intelligent transition method that enables the robot to perform smooth task transition.

In this paper, the concept and validation are presented in a kinematic control approach, where joint velocities are the inputs to the robot given desired task space velocities. In the following section, the source of discontinuity is explained. One approach dealing with this discontinuity is † Jaeheung Park is the corresponding author. park73@snu.ac.kr then proposed by introducing intermediate desired values as a new task specification.

The proposed approach is compared with the continuous inverse solution in [START_REF] Mansard | Continuous Control Law from Unilateral Constraints : Application to Reactive Obstacle Avoidance in Operational Space[END_REF], which modifies the inversion process. The continuous inverse approach is to compose a control law by modifying the inverse operation with an activation matrix. It is revealed that the two approaches are equivalent in the case of equal priority tasks. However, the new approach can be directly applied to hierarchically constructed tasks, where the continuous inverse solution is limited. Also, the new approach is easily applicable to many existing platforms since it does not require to change control or inverse process.

In addition, it is briefly discussed the case when there is not enough DOF for two tasks with different priorities. In such case the lower-priority task cannot be fully performed. Its control result is discussed to provide an intuition of the priority-based control and the effectiveness of the proposed approach for prioritized multiple tasks.

The experimental validations of the proposed approach were conducted in two scenarios: multi-point control using a 3 DOF planar robot in simulation and joint limit avoidance using a 7 DOF physical robot. The smooth and effective task transitions are demonstrated in both cases.

II. DISCONTINUITY FROM TASK ADDITION OR REMOVAL

When there are two tasks or features, x 1 and x 2 , the corresponding Jacobians are defined as J 1 and J 2 .

ẋ1 = J 1 q, ẋ2 = J 2 q (1)
When only one task, x 1 , is to be controlled, the associated joint velocity, q1 , can be computed as

q1 = J + 1 ẋ1 (2)
where (.) + denotes a pseudo-inverse of the quantity. The inverse could be replaced by any of the generalized inverses. Now, if another task, x 2 , is to be controlled together with x 1 , the joint velocity can be computed as

q = J + ẋ (3)
where

J = J 1 J 2 , ẋ = ẋ1 ẋ2 . ( 4
)
The joint velocity for both tasks (Equation ( 3)) might be discontinuous from the joint velocity for only the first task (Equation ( 2)) at the instance when the second feature is inserted.

One explanation for this discontinuous behavior is that it comes from the property of pseudo-inverse. The pseudoinverse computes a solution to meet the feature specification and minimize the 2-norm of q. Obviously, when we achieve one feature or two features the solution would not be continuous.

More physically meaningful explanation is as follows. When we control the first task or feature, the computed joint velocity, q 1 , indirectly, affect the task, x 2 , although this task is not interested during the execution. Then, the problem is that we later attempt to control or specify this feature with different values all of sudden. Therefore the joint velocity will show a discrete change.

Fundamentally, the problem is that the solution, q1 , from Equation (2) generates a value for the second task as J 2 q1 since we did not specify it. Then, when we specify a desired value, ẋ2 , ẋ2 = J 2 q1 (5)

Therefore, when we add this task, x 2 , the specification for x 2 in Equation ( 4) should be continuous from the value, J 2 q1 . Otherwise, the pseudo-inverse will give a solution with a discrete change. The discrete change would occur during task removal due to the same reason.

III. INTERMEDIATE DESIRED VALUE APPROACH

A. Concept of Continuous Transition Strategy

From the observation in the previous section, a proposed approach is to design a task specification of ẋ2 smoothly from J 2 q1 , where q1 is from Equation [START_REF] Baerlocher | Task-Priority Formulations for the Kinematic Control of Highly Redundant Articulated Structures[END_REF], such that the solution would not show any abrupt change.

and

ẋi 1 = ẋ1 ẋi 2 = h 2 ẋ2 + (1 -h 2 )J 2 J + 1 ẋ1 (7)
where h 2 is an activation variable from 0 to 1 when the task, x 2 , is inserted. The superscript, i, denotes the intermediate value of the quantity. Then the joint velocity is

q = J + ẋi . ( 8
)
This solution is continuous if we change h 2 from 0 to 1 continuously.

B. Generalization

An activation parameter for the first task, x 1 , can also defined as h 1 . Then, q = J + ẋi ( 9)

ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 J + 2 h 2 ẋ2 ẋi 2 = h 2 ẋ2 + (1 -h 2 )J 2 J + 1 h 1 ẋ1 . ( 10
)
This solution is in fact the same as the one for one task with activation matrix, H, in [START_REF] Khatib | Real-time obstacle avoidance for manipulators and mobile robots[END_REF]. Refer to Appendix I for the proof of equivalence.

C. More than Two tasks

When there are tasks with more than two activation parameters, Equation (10) can be extended. For example, the solution for the tasks with three activation parameters is the following:

q = J + ẋi (11)
And

ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 q{h 2 ,h 3 } ẋi 2 = h 2 ẋ2 + (1 -h 2 )J 2 q{h 1 ,h 3 } ẋi 3 = h 3 ẋ3 + (1 -h 3 )J 3 q{h 1 ,h 2 } , ( 12
)
where q{h m ,h n } denotes the solution for the tasks m and n with activation parameters of h m and h n . The derivation of Equations ( 11) and ( 12) is straight-forward from Equations ( 9) and [START_REF] Khatib | A Unified Approach for Motion and Force Control of Robot Manipulators: The Operational Space Formulation[END_REF].

IV. INTERMEDIATE DESIRED VALUE APPROACH FOR PRIORITIZED TASKS

A new task may be added with a different priority from the existing tasks. The solution of the previous section can be re-arranged in this case. For example, when we had a task set as x 2 and the associated Jacobian, J 2 , another higher prioritized task set, x 1 , can be inserted, such as joint limits or collision avoidance. The controller for x 2 only was

q2 = J + 2 ẋ2 . ( 13
)
Now with the higher prioritized task, x 1 ,

q = J + 1 ẋ1 + N 1 q0 (14)
where

N 1 = I -J + 1 J 1 (15) q0 = (J 2 N 1 ) + ( ẋ2 -J 2 J + 1 ẋ1 ). (16)
However, transition from Equation [START_REF] Yoshida | Humanoid Motion Planning for Dynamic Tasks[END_REF] to Equation ( 14) will have discontinuity. Therefore, a similar approach to the one introduced in the previous section can be applied.

q = J + 1 ẋi 1 + N 1 qi 0 ( 17
)
where

ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 q2 q2 = J + 2 ẋ2 qi 0 = (J 2 N 1 ) + ( ẋ2 -J 2 J + 1 ẋi 1 ). (18)
The above approach provides a continuous solution for 0 ≤ h 1 ≤ 1 because q2 from Equation ( 13) is the same as q from Equation ( 17) when h 1 = 0 (See Appendix II).

However, the above solution does not have an activation parameter for x 2 . This could be implemented by introducing where

h 2 . q = J + 1 ẋi 1 + N 1 qi 0 (19)
ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 qi 2 qi 2 = J + 2 h 2 ẋ2 qi 0 = (J 2 N 1 ) + (h 2 ẋ2 + (1 -h 2 )J 2 J + 1 h 1 ẋ1 -J 2 J + 1 ẋi 1 ). (20)
The equations ( 19) and ( 20) can be written as

q = J + 1 ẋi 1 + N 1 (J 2 N 1 ) + ( ẋi 2 -J 2 J + 1 ẋi 1 ) (21)
where

ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 J + 2 h 2 ẋ2 ẋi 2 = h 2 ẋ2 + (1 -h 2 )J 2 J + 1 h 1 ẋ1 (22)
A. What is the difference from the two approaches with or without priorities?

As we can see from Equation ( 21) and ( 22), the only difference in formulations for the approaches with or without priorities is the computational structure. The intermediate goal values, ẋi , are computed in the same procedure.

If the two tasks are independent, the two approaches with or without priorities should give the same result. The difference comes only if they are conflicting.

B. More than Two Tasks

The same extension approach as the previous solution applies. For the three task sets,

q = qi 1 + qi 2 + qi 3 (23) qi 1 = J + 1 ẋi 1 qi 2 = N 1 (J 2 N 1 ) + ( ẋi 2 -J 2 qi 1 ) qi 3 = N 1 N 2|1 (J 3 N 1 N 2|1 ) + ( ẋi 3 -J 3 ( qi 1 + qi 2 )) (24)
and

N 2|1 = I -(J 2 N 1 ) + J 2 N 1 (25) ẋi 1 = h 1 ẋ1 + (1 -h 1 )J 1 q{h 2 ,h 3 ,prio} ẋi 2 = h 2 ẋ2 + (1 -h 2 )J 2 q{h 1 ,h 3 ,prio} ẋi 3 = h 3 ẋ3 + (1 -h 3 )J 3 q{h 1 ,h 2 ,prio} , (26)
where q{h m ,h n ,prio} denotes the solution for the tasks m and n with activation parameters of h m , h n , and given priorities.

V. SIMULATION AND EXPERIMENTAL RESULTS

The proposed approach for task transition has been verified in ROBOTICSLAB [17] simulation environment and physical robot. ROBOTICSLAB provides not only physics-based simulation environment but also real time control module and programmable interface for the user.

In the following sub-sections, the multiple point control with different priorities is demonstrated in simulation with a 3-DOF planar robot. Then, the end-effector control with joint limit avoidance is executed on a 7-DOF manipulator. During the experiments, the kinematic control law is implemented as

ẋ = -λ (x -x d ) (27)
where x and x d are the current and desired task, and λ is a positive gain for decreasing task error.

A. Task transition in a 3 DOF planar robot: Simulation

In this section, the proposed task transition approach is demonstrated in a 3 DOF planar robot when a new task inserted with a lower priority than that of an existing task. The purpose of the section is to demonstrate the important concepts in a relatively simpler case. Also, the effect of the priorities is explained and shown when the tasks are partially conflicting due to the lack of degrees of freedom of the robot.

The two positioning tasks are defined as positioning of the end-effector and the 2nd link (Figure 1).

ẋ1 = ẋe,x ẋe,y ẋ2 = ẋlink2,x ẋlink2,y (28)
where (.) e,(.) and (.) link2,(.) denote the tasks related to the end-effector and the end point of the 2nd link, respectively. The robot begins to control only the end-effector position, x 1 .

Later the 2nd link position control is inserted after 2 seconds (t = 4 sec) (Figure 2 and3).

In Figure 2, the control results are plotted when the second task (positioning task of the 2nd link) is abruptly inserted as a low priority task. The results of the proposed continuous approach are plotted in Figure 3. The activation parameter h 2 increased from 0.0 to 1.0 for 0.5 seconds. The plot of The results related to 2nd link the first task is omitted for both cases since their tracking performances are good and almost the same. The abrupt change in Figure 2 is successfully removed by the proposed continuous approach in Figure 3.

Note the fact that the inserted task is not executed perfectly due to the lack of degrees of freedom for both tasks. This phenomenon is explained from the global and instantaneous point of view in Figure 4 and5, respectively. At the end of the control in Figure 4, the position of the 2nd link comes to the closest point to the desired position within the workspace, which is a circle centered at the end-effector position.

In the instantaneous view, the second task execution can be decomposed into the effect of higher-priority task and its control in the task null space (Figure 5).

ẋ2 = J 2 q = ẋ2,prio + ẋ2,null (29)
where ẋprio

= J 2 J + 1 ẋ1 ẋnull = J 2 N 1 q0 . ( 30
)
The term, ẋ2,prio , is the effect from the control of the previous task with a higher priority. Then, the control in the null space, ẋ2,null , set ẋ2 to be closest to its desired value. Although the inserted task cannot be fully performed because of not enough degrees of freedom for both tasks, The results related to 2nd link it is demonstrated that the task is executed optimally under priority-based control. Using recursive control structure with the proposed approach, many tasks can be controlled partially or completely with various levels of priorities.

B. 7-DOF manipulator with elbow joint limit: Experiment

Joint limit task for a 7 DOF manipulator is controlled through the proposed approach in this section. It must be considered as a higher prioritized task than other given tasks. When the joint limit task is inserted, continuity of task transition can be ensured by the proposed approach.

1) Definition of Prioritized Tasks: In this section, two task sets are considered with priorities. One is a joint limit task and the other is consisted of position and orientation of endeffector. The priority of joint limit task is higher than that of task related to the end-effector.

ẋ1 = qlimit , ẋ2 = v e ω e , (31)
where v e and ω e are the linear and angular velocities of the end-effector.

When the joint position is not in joint limit region, the end-effector positioning task is controlled by Now when the joint limit task with high priority is inserted without intermediate desired value method, the kinematic control with the priorities provides the solution of

q = J + 2 ẋ2 . ( 32
)
q = J + 1 ẋ1 + N 1 (J 2 N 1 ) + ( ẋ2 -J 2 J + 1 ẋ1 ). (33)
The result of this kinematic control through an abrupt change produces oscillatory behavior. The oscillation comes from rapid task transition at the activation border (Figure 8). It means that pushing toward and backward at the joint limit are iterated in sequence.

In the experiments, the end-effector of the robot is controlled to follow a line to reach a desired goal position. While the end-effector task is being executed, the Elbow joint reaches a joint limit buffer at right before 5 seconds (Figure 7, 8, and 9). Figure 8 shows instability due to the iterative abrupt changes in control due to the conflicting joint limit and end-effector control. The results of the proposed approach are plotted in Figure 9, which is explained in detail in the following sub-sections.

2) Activation Buffer: For the joint limit task execution, the activation parameter is defined as a function of the corresponding joint angle. In the activation buffer, the activation parameter increases from 0.0 at the beginning of the buffer to 1.0 at the limit. For inserting high-priority task, the activation parameter, h 1 , is defined as

h 1 =            1.0
if q elbow ≥ q upperLimit f (q elbow ) if q upperLimit -β < q elbow < q upperLimit 0.0 if q lowerLimit + β ≤ q elbow ≤ q upperLimit -β g(q elbow ) if q lowerLimit < q elbow < q lowerLimit + β 1.0 if q elbow ≤ q lowerLimit (34) where β is the buffer length (β = 0.2(rad)), q upperLimit = 1.6(rad), and q lowerLimit = -1.6(rad) in this example. The velocity with respecto to x velocity with respect to y Null space Fig. 5. Instantaneous task execution of the low-priority task, x 2 , in the 3 DOF planar robot. The execution of the inserted task with the lower priority is consisted of the effect from the high-priority task and its control in the task null space. However, when there is not enough DOF for both tasks, the control of the low-priority task in the null space partially achieve its goal. The term ẋd is the desired value for ẋ2 and ẋa is the actual command for ẋ2 due to the lack of DOF and dependency between task sets. 6. Activation parameter of task transition, h 1 , for the joint limit task in the 7 DOF manipulator. It is parameterized with the joint angle in the buffer zone. activation parameter is plotted in Figure 6. Equation (35) are activation functions varying in activation buffers.

0 0.2 0.4 0.6 0.8 1 T T ! $ F W L Y D W L R Q E X I I H U $ F W L Y D W L R Q E X I I H U -R L Q W O L P L W -R L Q W O L P L W Fig.
f (q elbow ) = 0.5 + 0.5sin( π β (q elbow -(q upperLimit -β )) - π 2 )
g(q elbow ) = 0.5 + 0.5sin( π β (q elbowq lowerLimit ) + π 2 )

(35)
3) Intermediate Desired Value: Discontinuity and instability of task execution occur when tasks are transitioned abruptly. For continuous task transition, the intermediate desired value for the joint limit task for the elbow is set to be ẋi

1 = h 1 ẋ1 + (1 -h 1 )J 1 J + 2 ẋ2 (36)
where h 1 is the activation function defined in Equation (34).

When the elbow joint is in the activation buffer zone, the joint Oscillation and instability due to the discrete task transition. When the activation border is reached, joint limit task is inserted with a higher priority, the oscillation starts due to the discrete task transition. The experiment is terminated at around 9.4 seconds due to the instability of the robot. (a) Joint velocity. (b) Elbow joint angle near the activation border. (c) Elbow joint angle plotted in a smaller scale (zoom-in view). The joint angle oscillates around activation border. limit task has to be inserted without discontinuity and instability. In Figure 9, there are no discontinuity and instability in the case which is controlled by the proposed approach. The joint position of elbow has been successfully controlled at the border of activation buffer and the activation parameter has varied continuously with respect to the position of the elbow joint. In this case, the elbow joint does not reach its limit. In the cases when more joints enter into the joint limit regions, the same approach for more tasks can be simply used.

VI. CONCLUSION

The intermediate desired value approach for task transition is presented in this paper. The proposed approach is to provide continuity during transition using existing control frameworks with minimal performance compromise. The approach without priorities are shown to be equivalent to continuous inverse method [START_REF] Mansard | Continuous Control Law from Unilateral Constraints : Application to Reactive Obstacle Avoidance in Operational Space[END_REF]. Then, it is further extended to the cases of tasks with priorities.

When task transition is executed by the proposed approach, the transitioned tasks can be executed smoothly with various levels of priorities including high priority and low priority. When there is not enough DOF of the manipulator, the task with a lower priority may not be fully performed but partially. The simulation result of such a case is shown to demonstrate that the proposed approach can be used successfully with the priority-based control framework.

We are currently working on extending the proposed approach for the dynamic control framework. And the issue of computation cost of the proposed approach will be investigated in the future.

APPENDIX I EQUIVALENCE WITH CONTINUOUS INVERSE METHOD

The proposed approach is to modify the desired values for x 2 to smooth the solution out. This solution can be re-arranged as follows: q = J + ẋ1 h 2 ẋ2 + (1h 2 )J 2 J + 1 ẋ1 (I.1)

q = J + I 0 (1 -h 2 )J 2 J + 1 h 2 I ẋ1 ẋ2 (I.2)
where

J + I 0 (1 -h 2 )J 2 J + 1 h 2 I = J + I 0 J 2 J + 1 0 + h 2 0 0 -J 2 J + 1 I = J + I 0 J 2 J + 1 0 + h 2 I - I 0 J 2 J + 1 0 = J + 1 0 + h 2 J + -J + 1 0 = J ⊕H | h 1 =1 (I.3)
During the above derivation the following property has been used:

JJ + 1 = J 1 J 2 J + 1 = I J 2 J + 1 (I.4) → J + JJ + 1 = J + I J 2 J + 1 J + 1 = J + I J 2 J + 1 (I.5)
The operator from ẋ to q in Equation (I.2) is equivalent to the continuous inverse J ⊕H in [START_REF] Mansard | Continuous Control Law from Unilateral Constraints : Application to Reactive Obstacle Avoidance in Operational Space[END_REF] when h 1 = 1. In a more general case when both h 1 and h 2 are not equal to zero or one, Equation (I.2) becomes q = J + h 1 I (1h 1 )h 2 J 1 J +

2

(1h 2 )h 1 J 2 J + 1 h 2 I ẋ1 ẋ2 (I.6) The equivalence with the continuous pseudo inverse can be easily derived.

J ⊕H = J + h 1 I 0 J 2 J + 1 0 + h 2 0 J 1 J + 2 0 I -h 1 h 2 0 J 1 J + 2 J 2 J + 1 0 (I.7)
APPENDIX II CONTINUITY IN PRIORITY BASED APPROACH Equation (17), when h 1 = 0, is

q = J + 1 J 1 J + 2 ẋ2 + N 1 (J 2 N 1 ) + ( ẋ2 -J 2 J + 1 J 1 J + 2 ẋ2 ) = {J + 1 J 1 + N 1 (J 2 N 1 ) + (J 2 -J 2 J + 1 J 1 )}J + 2 ẋ2 = {J + 1 J 1 + N 1 (J 2 N 1 ) + (J 2 N 1 )}J + 2 ẋ2 = {J + 1 J 1 + N 1 -N 1 + N 1 (J 2 N 1 ) + (J 2 N 1 )}J + 2 ẋ2 = {(J + 1 J 1 + N 1 ) -N 1 (I -(J 2 N 1 ) + (J 2 N 1 ))}J + 2 ẋ2 = {I -N 1 N 2|1 }J + 2 ẋ2 = J + 2 ẋ2
(II.1)

When h 1 = 0, q = J + 2 ẋ2 + h 1 {J + 1 -N 1 (J 2 N 1 ) + J 2 J + 1 }( ẋ1 -J 1 J + 2 ẋ2 ) (II.2)

  x1

Fig. 1 .

 1 Fig. 1. 3 DOF planar robot simulation: there are two task sets, x 1 and x 2 . (a) initial position of 3-DOF planar robot (b) positioning task of the end effector with high priority, x 1 , is being executed without the other task, x 2 . (c) Positioning task of 2nd link, x 2 , is added with a lower priority. (d) x 1 and x 2 are executed with priorities.

Fig. 2 .

 2 Fig. 2. Example of discontinuity of the control law when a task is inserted abruptly for the 3 DOF planar robot. (a) Position of 2nd link of 3 DOF planar robot during the control. (b)The inserted task set, x 2 , is plotted over time. The desired position of 2nd link is out of reach due to the high prioritized task, x 1 . The final position is the closest reachable position, as high prioritized task, x 1 , is executed. As soon as x 2 is inserted, undesirable jerky motion occurs at the end-effector and 2nd link paths.

Fig. 3 .

 3 Fig. 3. Example of task transition executed by proposed control approach for the 3 DOF planar robot. (a) Position of the 2nd link of 3 DOF planar robot during the control. (b)The inserted task set, x 2 , is plotted over time. The parameters of the simulation, initial position, desired position, trajectories, operating time and etc., are the same as the case of discontinuity (Figure2). Compared to the result of the abrupt transition, the paths of the end-effector and 2nd link are smooth and there is no undesirable motion.

3 link 3 Fig. 4 .

 334 Fig.[START_REF] Baerlocher | An Inverse Kinematic Architecture Enforcing an Arbitrary Number of strict priority levels[END_REF]. Task execution of the two tasks: the end-effector, x 1 , and the 2nd link, x 2 , of the 3 DOF planar robot with different priorities. The final position of the 2nd link, which has the lower priority, does not reach its goal position but to the closest position from the desired position in the reachable positions (the circle centered at the final position of the end-effector).

Fig. 7 .

 7 Fig. 7. Joint limit task execution of the 7 DOF manipulator with positioning and orientation tasks. (a) Initial position of 7-DOF manipulator (b) Positioning and orientation task, x 2 , execution without joint limit (c) x 2 execution with joint limit (d) Return positioning task, x 2 , execution as elbow joint is extricated from joint limit (e) Final position of 7-DOF manipulator

Fig. 9 .

 9 Fig. 9. Execution of joint limit task by the proposed continuous control law. The joint limit task is being executed successfully without any instability with the end-effector task simultaneously. (a) Joint velocity (b) Elbow joint angle: The joint angle converges within activation buffer.(c) Elbow joint angle plotted in a smaller scale (zoom-in view)

VII. ACKNOWLEDGMENTS

This work was supported by the National Research Foundation of Korea (No.2010-0005799) and Seoul R&BD Program (No.JP100106).


Antoine Beylot
email: a.beylot@brgm.fr

G Hernandez-Rodriguez

J Villeneuve

A numerical approach for compiling full Physical Supply-Use Tables (PSUTs) under conflicting information

Keywords: Constrained optimization technique, Mass balance, Physical Input-Output Tables (PIOTs), Physical Supply-Use Tables(PSUTs)

. Accordingly there is a strong need to limit the cost and time required for PSUTs and PIOTs compilation while at the same improving their reliability. This work proposes a numerical approach to balance full PSUTs under conflicting information using the notation of constrained optimization. The mass balancing identities of PSUTs, in terms of products and activities, are applied in a mathematical technique which fulfills all requirements of constrained optimization techniques. Following the theoretical framework as defined in a first part, a tentative PSUT for the Netherlands 2006 is presented as a case study (being validated). Such an approach could constitute a major advance for the practice of compiling PSUTs before deriving PIOTs, since it aims to remove the necessity of manually tracing conflicting information. In addition, this article gives some guidelines for future research (e.g. a full PSUT compiled for France for the year 2006) which could help to make PIOTs more relevant and cost-effective.

INTRODUCTION

Physical Supply-Use and Input-Output Tables (PSUTs and PIOTs) provide a comprehensive picture of anthropogenic material flows within the economy and in interaction with the natural environment. They complement the corresponding monetary tables (Monetary Supply-Use Tables, MSUTs) by registering flows of physical products, extraction of materials from nature, supply and use of wastes, emissions to nature and stock changes [START_REF] Giljum | Conceptual foundations and applications of physical input-output tables[END_REF][START_REF] Gravgård Pedersen | The system of environmental and economic accounts-2003 and the economic relevance of physical flow accounting[END_REF][START_REF] Hoekstra | Constructing physical input-output tables for environmental modeling and accounting: Framework and illustrations[END_REF]. There are three main advantages that can be identified by using PSUTs and PIOTs: integration of physical data, improvement of monetary statistics and environmental-economic modelling [START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF]. Of course these advantages come at a price. The data requirements and cost of producing PSUTs and PIOTs are large, and increase when the tables become more complex. We propose a mathematical model based on which the PSUTs and PIOTs can be produced more costefficiently. Such a technique is currently being developed by the writers of the present paper. It includes a cost-saving strategy consisting ofdata reconciliation &error estimation.

Data reconciliation is a technique to improve the quality of measured data. These measurements are inherently inaccurate or subject to failures. Using erroneous data for accounting analysis and decision-making may yield distorted conclusions and result in improper decisions [START_REF] Narasimhan | Data reconciliation & gross error detection,Data reconciliation & gross error detection: an intelligent use of process data[END_REF]. Accurate data is therefore essential for compiling PSUTs and subsequently analyzing material flows. One measure of data inaccuracy is the consistency with regard to the mathematical models describing the accounting system. Among the more classical models used for describing a functioning are the balance relationships (mass, component, species, enthalpy, etc.). If all of these models are structurally perfectly known, some of them depend on parameters which are difficult to assess. Therefore, it becomes very hazardous and mathematically not correct to reconcile operation data with regards to an uncertain model without taking this last fact into account. We propose to take this information into account in the reconciliation procedure, assuming some knowledge about the precision of the values of the input parameters.

This paper provides an overview of the literature on PSUTs and PIOTs. It also proposesan optimization model (currently under validation) which could help to make the PSUTs and PIOTs more relevant and cost-effective.Rudimentary PSUTs for the Netherlands for the year 2006 are only presented as an illustration for future applications.In a more prospective view, this paper proposes another approach for evaluating uncertainty related to model input parameters: the possibilistic or fuzzy approach [START_REF] Guyonnet | Comparing two methods for addressing uncertainty in risk assessments[END_REF].

SUPPLY-USE AND INPUT-OUTPUT TABLES

General framework

The framework of PSUTs and PIOTs, including their accounting identities, has been widely detailed in the literature, in particular by [START_REF] Giljum | Conceptual foundations and applications of physical input-output tables[END_REF][START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF][START_REF] Hoekstra | Constructing physical input-output tables for environmental modeling and accounting: Framework and illustrations[END_REF][START_REF] Schmidt | Documentation of the final model used for the scenario analyses[END_REF]. We consider the latter's description and notations in what follows. Every table mentioned below corresponds to a specified period and geographical area (e.g. France for the year 2006). The Supply matrix, V, of dimensions activities by products, reports the supply of products per human activity. The Use matrix, U, of dimensions products by activities, details the intermediary consumptions of products per human activity. These two tables are completed by the Import and Export vectors, N and E, of dimensions products by one, which report the exchanges of products with the rest of the world. Finally the vector Y, of dimensions products by one, stands for the final consumption vector. This first set of tables is traditionally accounted for in monetary terms (i.e. MSUTs) and is correspondingly reported in physical units in PSUTs.

PSUTs additionally include the environment as a source of raw materials (matrix R of dimensions resources by activities) and as a sink for residuals and emissions (respectively matrices W V and B of dimensions products by activities and emissions by activities). Finally, W U and ∆S, of dimensions products by activities, respectively represent the use of residuals and the addition to stocks of products and residuals. The accounting identities that structure the PSUTs are based on the material balance principle (Figure 1). On the one hand, on a product perspective:

𝑉 + 𝑁 = 𝑈 + 𝑌 + 𝐸 (1)
And on the other hand, on an activity perspective:

𝑉 + 𝑊 𝑣 + 𝐵 + ∆𝑆 = 𝑈 + 𝑊 𝑢 + 𝑅 (2)
The balanced PSUTs may finally be converted into PIOTs by using one of the following assumptions: the product technology assumption, the industry technology assumption, and the assumption of fixed industry sales structure or the assumption of fixed product sales structure. This step is widely detailed and discussed in the literature [START_REF]Eurostat, Eurostat Manual of Supply, Use and Input-Output Tables[END_REF]

Applications and existing case studies

Whereas PSUTs are better suited in an accounting perspective, they most often need to be converted into PIOTs to address environmental issues [START_REF] Hoekstra | Constructing physical input-output tables for environmental modeling and accounting: Framework and illustrations[END_REF]:

 PIOTs can primarily be intended to derive environmental information: environmental pressure indicators, composition of products, element cycles in the economy, dematerialization indicators and physical trade balance indicators;

 PIOTs can be further used for environmental modelling purposes, either to analyze the impact of a certain change in final demand on output (impact analysis) or to impute requirements in raw materials and emissions to a specific final demand (imputation to final demand).

In their literature review [START_REF] Giljum | Conceptual foundations and applications of physical input-output tables[END_REF][START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF] list several compiled PIOTs among which: an Austrian PIOTs, for the year 1983, which was the first attempt to calculate a PIOTs [START_REF] Kratena | Eineokologischevolkswirtschaftlichegesamtrechnung fur osterreich. Die Umwelt Input Output Tabelle 1983. Wien: Institutfürsozial-, wirtschafts-und umweltpolitischeForschung[END_REF]; full PIOTs for Germany (for the years 1990 and 1995) [START_REF] Stahmer | Physische Input-Output Tabellen[END_REF][START_REF]StatistischesBundesamt. A Physical Input-Output Table for Germany 1995[END_REF], for Denmark (1990 updated to 2002) [START_REF] Gravgård Pedersen | Physical Input-Output Tables for Denmark[END_REF][START_REF] Mulalic | Material Flows and Physical Input Output Tables -PIOT for Denmark 2002 based on MFA[END_REF] and for New Zealand (1997/98) [START_REF] Mcdonald | Development of a New Zealand Physical Input-Output Table[END_REF]; an aggregated PIOTs for Italy (1995) [START_REF] Nebbia | Contabilità monetaria e contabilità ambientale[END_REF] and a detailed PIOTs for Finland (1995) [START_REF] Mäenppä | Physical input-output tables of Finland 1995solutions to some basic methodological problems[END_REF]. This list can be additionally expanded to the full PIOTs compiled for United Kingdom (covering the period 1997-2004) [START_REF] Wiedmann | National and Regional Physical Accounts (Material Flows) for the United Kingdom[END_REF] and for 22 countries of the European Union (2003) [START_REF]Overall mapping of physical flows and stocks of resources to forecast waste quantities in Europe and identify life-cycle environmental stakes of waste prevention and recycling[END_REF].

Despite converting PSUTs to PIOTs is generally necessary, PSUTs may also be directly used for environmental modelling. In particular, PSUTs may be used to forecast future waste quantities, environmental impacts and benefits related to changes in economic activities and policies [START_REF] Schmidt | 25-year forecasts of the cumulated physical stocks, waste generation, and environmental impacts for each scenario for EU-27 and for the case study countries[END_REF].

Constructing PSUTs involves compiling data which to a larger or smallerextent are inconsistent.A research problem, which has not yet found its final solution, is how to reconcile various sources of information in balancing consistent PSUTs, taking into account all information in themost efficientmanner possible. The purpose of the following sections is to present some directions for balancing PSUTs.

Compiling PSUTs: data inventory, uncertainty and inconsistency

In a literature review, [START_REF] Giljum | Conceptual foundations and applications of physical input-output tables[END_REF] report four main methodological differences between existing tables: the level of sector aggregation (from 27 activities in the Danish table to 59 in the German one), the system boundaries, including or not plants and forests, the inclusion or exclusion of different material categories such as water and air, and the base year.

In order to complement this list of methodological differences between existing tables, we performed the review of five studies for which the compilation procedure was sufficiently documented (respectively [START_REF]Overall mapping of physical flows and stocks of resources to forecast waste quantities in Europe and identify life-cycle environmental stakes of waste prevention and recycling[END_REF][START_REF] Mcdonald | Development of a New Zealand Physical Input-Output Table[END_REF][START_REF] Mulalic | Material Flows and Physical Input Output Tables -PIOT for Denmark 2002 based on MFA[END_REF][START_REF]StatistischesBundesamt. A Physical Input-Output Table for Germany 1995[END_REF][START_REF] Wiedmann | National and Regional Physical Accounts (Material Flows) for the United Kingdom[END_REF]). From this review, both the data inventory and the treatment of inconsistencies appear as being treated differently from one study to another whereas they are of core importance in the compilation.

PSUTs compilation is mainly driven by the availability of statistical data. The latter generally originate from different sources and are in some cases obtained from rough assumptions in the absence of more accurate information. Focusing on the Supply of products table (V) and on the Use of products table (U), without considering imports/exports, emissions, waste, stocks and resources in a first approach, four distinct kinds of data are observed to be usually implemented (Table 1):

1. Sectorial data on the physical supply and use of products, expressed in mass units. These data are directly extracted from national statistical databases. This is the "ideal" case, in the sense that these data can supposedly be directly implemented in the PSUT as such, without any conversion. 2. Sectorial data on the physical supply and use of products, expressed in other units than mass, e.g. in volume or number of items. These data are extracted from national statistical databases and are converted into masses by use of adequate factors. 3. Coefficients of the monetary supply and use tables, as for example annually reported by Eurostat, those need to be converted into physical terms by use of product prices. Import/export commodity prices per net weight may be used as surrogates for their domestic supply/use. 4. Process-specific data, extracted from Life Cycle Inventories and expressing the amount of inputs per unit of output of a specific product. These need to be up-scaled before their implementation into the Use of products table (U). Yes Upscaling [START_REF]Overall mapping of physical flows and stocks of resources to forecast waste quantities in Europe and identify life-cycle environmental stakes of waste prevention and recycling[END_REF] PSUTs compilation may require combining several types of data. This is in particular the case for the compilation of the PSUTs of the 22 countries of the European Union for the year 2003 [START_REF]Overall mapping of physical flows and stocks of resources to forecast waste quantities in Europe and identify life-cycle environmental stakes of waste prevention and recycling[END_REF]. These PSUTs were compiled by primarily using mass data for products supplied and used, and required as a complement to convert monetary data, process-specific data and number of items into masses in order to cope with missing data. Finally, it is worth reminding that this short review on data types focused on tables V and U (Supply and Use), but that similarly, data from multiple sources are also necessary to compile the tables of imports/exports, emissions, resources, waste generation and use, and stocks.

All the data aforementioned are intrinsically uncertain and convey errors. Uncertainty originates from many sources. For each type of data used in the compilation of Supply and Use tables (V and U), a few causes of uncertainty are reported in the following:

 Physical supply and use data per sector, in mass units, may for example be inaccurate due to errors in reports from enterprises or due to errors in data aggregation;

 Similarly to mass data on products supplied and used per sector, statistical sectorial data on the number of items produced and used may be inaccurate. In addition, factors to convert the number of items into masses are in most cases rough estimates. As an example, converting the number of pairs of shoes annually produced at the scale of a country into a mass (e.g. in tons) requires setting an average mass per pair of shoes.This conversion factor is necessarily inaccurate at the scale of a country production;

 The price per mass of a product category used to convert monetary tables into physical tables may be representative for supplies, but generally fails in representing the different uses of the product (in distinct sectors of the economy). The assumption of homogenous sectorial prices is not valid, as highlighted by discrepancies between monetary and physical input-output model outcomes [START_REF] Weisz | Physical and monetary input-output analysis: What makes the difference?[END_REF];

 Process-specific data extracted from Life Cycle Inventories databases are generally average data from a small number of plants and processes, for a given country. Thesedata may therefore not be representative at the scale of the production of a country.

The compilation of PSUTs leads to inconsistencies in the mass balancing, as obviously suggested by the large number of uncertainties associated with data and as usually observed in case studies. However, data uncertainty is generally not addressed, whereas inconsistences are handled manually by modifying coefficients of the Supply-Use Tables (as e.g. performed by [START_REF]Overall mapping of physical flows and stocks of resources to forecast waste quantities in Europe and identify life-cycle environmental stakes of waste prevention and recycling[END_REF]).Balancing PSUTs by manual correction are often rather costly to maintain and not easy to document even with the help of modern electronic data processing. Consequently a balancing tool appears necessaryto allow the user to build balanced-PSUTs with considering data uncertainties, and therefore making PSUTs more relevant and cost-effective. Such a tool is currently being developed by the writers of the present paper.

BALANCING PSUTS

At the start of balancing an estimate is available for every entry of the PSUTs. In spite of all efforts on compiling real-preliminary estimates, it has to be expected that inconsistencies in the estimates remain. How can inconsistencies be detected and how can they be solved in order to get balanced PSUTs?

In our knowledge, no general theory or useful mathematical programs are available. However, in balancing it is very important to follow a systematic approach to solve the problems. The balancing process is particularly important in the case of detecting and correcting many weaknesses of primary statistics. Moreover, balanced PSUTscan be used for many other purposes than just balancing the national accounts (as mentioned previously). Some experiences [START_REF]Eurostat, Eurostat Manual of Supply, Use and Input-Output Tables[END_REF] show thatcombination of manual and automatic statistical techniques and proceduresis the best workable solution to establish a supply and use system.

Preliminary proposed-strategy: data reconciliation & error estimation

Data reconciliation is a technique that has been developed to improve the accuracy of measurements by reducing the effect of random errors in the data. The principal difference between data reconciliation and other filtering techniques is that data reconciliation explicitly makes use of mass balance identities and obtains estimates of the variables by adjusting measurements so that the estimates satisfy the mass balance constraints [START_REF] Narasimhan | Data reconciliation & gross error detection,Data reconciliation & gross error detection: an intelligent use of process data[END_REF].Thereby, data reconciliation improves the accuracy of sectorial national statistical data by adjusting the measured data so that they satisfy the material balance identities.

In general, data reconciliation can be formulated by the following constrained least-squares optimization problem:

𝑀𝑖𝑛 𝑦 𝑖 -𝑥 𝑖 2 𝑛 𝑖=1 (3)
Subject to

𝑔 𝑘 𝑥 𝑖 = 0 𝑘 = 1, … , 𝑚 (4)
The objective function (Equation 3) defines the total sum square of adjustment made to measurements; where y i is the measurement and 𝑥 𝑖 is the reconciled estimate for variablei. Equation 4 defines the set of model constraints (e.g. material mass balance). The deterministic natural laws of conservation of mass (or energy for process engineering) are typically used as constraints for data reconciliation because they are usually known. These types of constraints that are imposed in reconciliation depend on the scope of the reconciliation problem. Furthermore, the complexity of the solution techniques used depends strongly on the constraints imposed. For example, if we are interested in reconciling only the mass flow rates, then the material balances constraints are linear in the mass flow variables and a linear data reconciliation problem results. On the other hand, if we wish to reconcile process data (e.g. temperature or pressure measurements along with flows),then a nonlinear data reconciliation problem occurs. Note that the preliminary proposed-strategy is based on the assumption that only random errors are present in the account measurements which follow a normal (Gaussian) distribution, with zero mean and a known variance-covariance as described in what follows: a tentative PSUTs for the Netherlands.

CASE STUDY

In order to obtain a good understanding of the issues in data reconciliation for future more realistic problems (e.g. full PSUTs compiled for France for the year 2006), a simple case study is introduced here, in order to highlight the assumptions to estimate PSUTs considering data uncertainties.The question is: whatis the most efficient way to achieve our objective?

Let us consider the reconciliation of tentative PSUTs. Initially, all mass flow rates are assumed to be known: tables of supply (V), use (U), imports (N), exports (E), stocks changes (∆S), needs fulfillment (y), supply and use of residuals (W V and W U ), emissions (B) and resources (R). The flow measurements contain unknown random errors. Note that the preliminary proposedstrategy is based on the assumption that only random errors are present in the account measurements which follow a normal (Gaussian) distribution, with zero mean and a known variance-covariance. For that reason, the material input and output do not balance. The aim of reconciliation is to make minor adjustments to the measurements in order to make them consistent with the material balances. The adjusted measurements, which are referred to as estimates, are expected to be more accurate than the measurements.

Tentative PSUTs for the Netherlands 2006

As an illustration of physical accounting we have used the tentative PSUTs for the Netherlands for 2006 compiled by [START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF]. According to the author, this is very much a quick-dirty effort to show what the numbers look like for a Western country. It is however also aimed at regaining inhouse experience for future research.

Tables 2 &3 show the PSUTs. The economy has been split into four parts (agriculture, mining, industry and services) which are relevant for material flows. Note that the PSUTs source data is available at about 50-60 industries and many subcategories of wastes, natural resources, etc. The imports and exports of good are part of the Mass Flow Analysis (MFA) statistics produced by the department of environmental accounts. Other data is derived from several accounts such as: air emission, waste, energy and water. The other components, for which physical data is not available, are estimated using the monetary values and appropriate prices from the import and export data [START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF]. It is assumed that the flows of all the PSUTs are known and that these measurements contain random errors (they follow a normal distribution, with zero mean and a known variancecovariance). If we denote the true value of the flow rate 𝑖 by the variable 𝑥 𝑖 and the corresponding measured value by 𝑦 𝑖 , then we can relate them by the following equation:

𝑦 𝑖 = 𝑥 𝑖 + 𝜀 𝑖 𝑖 = 1 … 𝑛 (5)
Where𝜀 𝑖 is the random error in measurement 𝑦 𝑖 .Flows must fulfill the balancing identities presented in Equations 1-2. This means that commodity (q) and industry totals (g) have to be equal in both tables (Figure 1). Obviously the measured values do not satisfy the above equations, since they contain random errors (Tables 2 and3). It is desired to derive estimates of the flows that satisfy the above flow balances. Intuitively, we can impose the condition that the differences between the measured and estimates flows should be as small as possible (Equation 3). Moreover, we can assume that the error variances for all the measurements are known. Thus, the reconciliation problem is a typically constrained optimization problem with the objective function given by Equation 3 and the constraints given by Equations 1-2. The solution of this optimization problem can be obtained mathematically by means of the preliminary proposed-strategy (being evaluated).It could enhance the benefits of applying data reconciliation techniques. Indeed, by taking into account all the available information about the data, it can prevent from erroneous decisions. A current study aims to validate numerically these assumptions and to extend the strategy to more complex PSUTs (e.g. full PSUTs compiled for France for the year 2006).

CONCLUSIONS AND RECOMMENDATIONS

Physical supply-use and input-output tables offer a detailed description of material flows within an economy and in interaction with the environment. These tables can be used either to directly derive environmental information or to perform environmental modelling. However, despite PIOTs and their source PSUTs have clearly appeared as popular in the last decades, as shown by the increasing number of publications and case studies on this issue; their compilation still requires large scale efforts. In addition, despite data estimates generally implemented in PSUTs come from multiple sources and are more or less accurate, errors related to these estimates are generally not accounted for in the PSUTs compilation. Accordingly there is a strong need to limit the cost and time required for PSUTs and PIOTs compilation while at the same improving their reliability. The strategy of data reconciliation &error estimation proposed to reconcile tables would benefit greatly by addressing these two core issues. First it would significantly reduce the time necessary to produce PSUTs. Secondly, it would also improve the accuracy of data by adjusting the measured values so that they satisfy the process constraints.

  and will not be considered further in this paper. The derived PIOTs reportin columns the intermediary consumptions of products, emissions, resource consumptions, stocks changes, waste generation and use associated with the production of one additional unit (e.g. 1 ton) of the corresponding product or activity.

	Balanced PSUT	Activities Import	Needs Fulfilment	Export Total
	Products	V'	N			Q
	Total	g'
	Products	U		Y	E	q
	Stock changes	-∆S
	Supply of residuals	-Wv
	Use of residuals	Wu
	Resources	R
	Emissions	-B
	Total	g'
		Figure 1: BalancedPSUTs [17].

Table 1 :

 1 Literature review of data types and sources in V and U compilation

	Input data for PSUTs (V and U) compilation	Data source	Need for data conversion	Conversion factor	Example of study
	Statistical annual data of production and use, per sector. In mass units	National Statistical Institutes	No	N.A.	[1, 14, 20, 22,]
	Statistical annual data of use and production, per sector. In units other than mass	National Statistical Institutes	Yes	Mass per unit	[1, 20, 22]
	Monetary Supply Use Table	National Statistical Institutes	Yes	Mass per unit monetary	[1, 13]
		Life Cycle Inventories
	Life Cycle Inventories	databases (e.g.
		ecoinvent)

Table 2 :

 2 Physical supply table for the Netherlands 2006-millions tons[START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF]

				Industries		Imports Cons. Total
			Agr.	Min.	Ind.	Serv.
		Agr.	39	0	0	0	24		63
	Commodities	Min. Ind.	0 0	113 0	4 218	4 6	157 144		277 368
		Serv.	0	0	0	1	0		1
	Supply resid.		4	0	48	7	11	9	69
	Emissions		10	3	47	107	0	37	203
	Total		53	115	317	124	337	46	981

Table 3 :

 3 Physical use table for the Netherlands 2006-millions tons[START_REF] Hoekstra | Physical input-output tables: Developments and future[END_REF]

				Industries		Final demand		Total
			Agr.	Min.	Ind.	Serv.	Cons.	Exp.	Inv.
		Agr.	2	0	30	1	6	16	0	56
	Commodities	Min. Ind.	2 12	8 0	210 127	11 28	1 47	80 183	0 4	312 400
		Serv.	0	0	0	0	0	0	0	0
	Use of resid.	Waste	1	0	53	12	0	13	0	67
	Raw	Ores/fuel		175
	materials	Water	208	5	3652 11179	729	0	0	15773
	Total		226	189	4072 11230	784	293	4	16608


Camille Fauth
email: camille.fauth@gmail.com

Béatrice Vaxelaire

Jean-François Rodier

Pierre-Philippe Volkmar

Fayssal Bouarourou

Fabrice Hirsch

Rudolph Sock

A spatiotemporal prospective study of speech in patients with or without recurrent laryngeal nerve paralysis after thyroid surgery

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Background

The current preliminary study is part of programmatic research work that is being conducted in our laboratory. This specific research is being carried out in collaboration with the Department of Surgical Oncology of the Paul Strauss Centre in Strasbourg. This centre offers a privileged advantage for our research as it receives, every year, between 650 and 750 patients suffering from some pathology of the thyroid gland. Consequently, it is possible to carry out systematic analyses of the diverse consequences of thyroid surgery on the voice and speech of the patients, and also to conduct longitudinal studies on a large number of speakers.

Let us recall that the thyroid gland, located at the base of the neck, is responsible for secretion of hormones indispensable for important functions of the organism such as cardiac rhythm, intestine activity, body temperature, transformation of fats into sugar, appetite, humour, sleep, body weight equilibrium and growth in the child.

During ablation of this gland, the surgeon takes great care to preserve the two recurrent nerves, responsible for mobility of the vocal folds. However, after surgery, at the laryngeal level, the voice and speech of the patient may be altered with the following consequences: a breathy or hoarse voice, irregular laryngeal vibrations or diplophonia, modification of voice periodicity, augmentation of jitter and shimmer, and of noise relative to harmonics, a reduction of intensity, effects of voice fatigue, alteration of airflow, etc. In some cases, when one recurrent nerve is inadvertently injured during surgery, unilateral paralysis may appear.

The aim of this preliminary study is to analyse some spectral and temporal characteristics of the speech of patients with recurrent paralyses and also of patients without diagnosed paralyses but with severe alteration of their voice and speech after thyroid surgery. This longitudinal study falls within the perturbation and readjustment paradigm preoccupied mainly with evaluating the flexibility of the speech production and perception system, determining the range of linguistically tolerated deviations from speech "targets", and ultimately elaborating viability constraints in articulatory and acoustic terms.

Method

Patients

Today, 21 speakers (19 women and 3 men) have been recorded. This corresponds to 7 hours of audio recordings. All patients are native French speakers, between 52 and 60 years old. A short biography of each subject is kept in the medical file, with information pertaining also to pre-and post-surgery. For this specific work, 4 patients have been analysed, 2 patients with unilateral paralysis and 2 patients without unilateral paralysis, all the four patients have undergone thyroidectomy. When a patient had only been recorded in post-surgery phases, data are acquired from a control speaker, matched with age and gender, with no voice pathology and no surgery in the vocal tract area.

Corpus and recordings

The corpus consisted of sustained vowels, VCV nonsense words and a read text. The sustained vowels were the 3 extreme vowels /i, a, u/, which allow exploring the maximum vowel space of each speaker. As concerns the VCV sequences, when V1 was vowel /i/, V2 was vowel /a/, and when V1 was vowel /a/, V2 was vowel /i/. The consonant was one of the 6 plosives /p, t, k, b, d, g/. The combination of these sounds allowed creating 12 nonsense words that were inserted in the carrier sentence "Cet________ça". Vowels /i/ and /a/ were specifically chosen in order to observe the acoustic consequences of larynx trajectory from a low position for /i/ to a high position for /a/, and vice versa. Consonants /b, d, g, p, t, k/ were retained in order to observe for effects of place of articulation on intra-oral pressure in relation to timing between burstrelease and onset of vocal fold activity (VOT as defined by Lisker & Abramson, 1964), or timing between burst-release and onset of a clear formant structure (VOT as specified by Klatt, 1975). These six consonants were also chosen as they allow comparing entirely voiced sequences V1/b, d, g/V2 with partially voiced sequences V1/p, t, t/V2. The latter sequences should impose more constraints on vocal fold activity since they require alternating closed and open configurations of the glottis. Reading a text allows evaluating the effects of vocal fatigue on the speaker's production.

During recordings, the vowel to be pronounced was presented to the subject on a cardboard in a graphic form /i, a, ou/ which the experimenter was sure would elicit the desired pronunciation. The three vowels were presented in random order, and the subject was required to pronounce each vowel, sustained for about 3 seconds, 10 times per vowel. The same method was also used to record the VCV nonsense words. The short text was read at a self-selected normal speaking rate. The data were obtained in 3 phases: 1) a pre-surgery phase, the day before surgery, which gives us the patient's unaltered reference voice; 2) a first post-surgery phase, the day after surgery, when the patient's voice is altered in varying degrees, due to paralysis or not of the vocal folds; 3) a second post-surgery phase, 15 days after surgery, which allows observing probable voice and speech recuperation. If unilateral paralysis is diagnosed, the patient will be recorded once a month during his or her entire speech therapy.

Measures

Measurements obtained were: VOTs (defined supra), VTT or Voice Termination Time, i.e. voicing decay after closure within the voiceless consonant hold in VCV sequences; vowel (V1 and V2) and consonantal durations; F0, F1 and F2 values in sustained vowels.

Hypotheses

It is hypothesised that in the first post-surgery phase, due to potential paresis or paralysis of the vocal folds: 1) VOT and VTT values, which respectively require precise offset and onset of laryngeal activity, relative to articulatory closure and release of the consonant, would be significantly longer; 2) Perturbation of the timing of VOT and VTT, as intra-segmental consonantal measures, would consequently have an impact on the duration of the consonantal duration, which should be longer, at least in the VOT transitional phase towards V2; 3) Difficulty in controlling voicing would also affect vowel durations; 4) Altered laryngeal activity would modify coupling between the larynx and the vocal tract, such that perturbation of voicing at the source would also have an impact on supraglottal resonances, i.e. on F1 and F2 values, and may affect the size and shape of vowel spaces; 5) Difficulty in controlling voicing would expectedly affect F0 values; 6) With speech therapy, the abovementioned parameters would be less modified, and may begin to resemble reference values observed for the patient in the pre-surgery phase, or those of the control speaker when no pre-surgery recordings are available.

Results

Preliminary results, based on mean values and standard deviations when patients had difficulties producing sufficient items, and based on ANOVAs when patients succeeded in carrying out the tasks required of them, show that: 1) VOT and VTT values had a tendency to be longer in post-surgery phases, and especially in the first post-surgery one; variability as highlighted by standard deviations, was significantly more pronounced in post-surgery productions (typical e.g. of Klatt's VOT mean value for /iba/: 17ms for the control speaker (CS), std=4ms; 31ms in second post-surgery phase, std=15 ms; 24 ms in third post-surgery phase, std=7ms; 19ms in the fourth post-surgery phase, std 5ms).

Figure 1. VOT values for [iba] for patient PLPH with unilateral paralysis

Figure 2. VTT values for [ika] for patient PLPH with unilateral paralysis

This finding is also true for V1 and V2 durations. See Figure 3 for a typical example of vowel durations for a patient (PLPH) with unilateral paralysis (example of vowel durations for /iba/: V1 /i/ 67ms for the control speaker (CS), std=15ms; 138ms in second post-surgery phase, std=18ms; 122ms in third post-surgery phase, std=28ms; 92ms in the fourth post-surgery phase, std 17ms and V2 /a/ 68ms for the control speaker (CS), std=11ms; 80ms in second post-surgery phase, std=18ms; 70ms in third postsurgery phase, std=14ms; 73ms in the fourth post-surgery phase, std 8ms).

Figure 3. V1 and V2 durations for [iba] for patient PLPH with unilateral paralysis

Vowel spaces are significantly modified (p<.01) in the first post-surgery phase compared with the pre-surgery one; they look more conventional in the second postsurgery phase for patient (NPPG) without unilateral paralysis (see Figure 4). Figure 5 shows vowel space areas for subject CS (control speaker) and subject PLPH (patient with unilateral paralysis). Here also, differences are statistically significant (p<.01). F0 values diminish significantly (p<.01) in the first post-surgery phase, compared with the pre-surgery phase (see Figure 5 for an illustration of this fact). Figure 6 shows that vowel /i/ varies from 219 Hz (std=7) in the pre-surgery phase to 187 Hz (std=15) in the first post-surgery one. Reduction of F0 values is also observed for /a/ and /u/, whose values go from 187 Hz (std=17) to 157 Hz (std=16) and from 219 Hz (std=13) to 157 Hz (std=11), respectively. In the second post-surgery phase, all F0 values increase; F0 is at 209 Hz for /i/ (std=4), 193 Hz for /a/ (std=8) and 212 Hz (std=12) for /u/.

Figure 6. F0 for patient NPPH without unilateral paralysis

Reduction in F0 values is clear-cut for patients with unilateral paralysis in the post-surgery phase, compared with those of their respective control speakers (see Figure 7 for a typical example). It can be seen in the figure that F0 for CS is at 219 Hz (std=7) for /i/, 187 Hz (std=30) for /a/ and 219 Hz (std=31) for /u/. The patient with paralysis, in the second post-surgery phase, has lower F0 values, compared with those of the corresponding control speaker (CS): /i/ is at 172 Hz (std=12), /a/ at 155 Hz (std=17) and /u/ at 160 Hz (std=13). Results that are being analysed for the third post-surgery phase (130 Hz (std=50) for /i/, 130 Hz (std=40) for /a/ and 126 Hz (std=49) for /u/) and the fourth post-surgery phase (174 Hz (std=30) for /i/, 174 Hz (std=16) for /a/ and 177 Hz (std=33) for /u/) seem to indicate increase of F0 values, with time and speech therapy.

Conclusions and discussion

As concerns hypotheses 1 and 2, VOT and VTT values tend to be longer in post-surgery phases. However, such results shoud be taken with caution as they rely mostly on nonsense words, and also reveal high inter-speaker variability. Systematic analyses of these parameters will be pursued in order to be able to evaluate their robustness in the current investigation.

Vowel duration has proven to be a relevant parameter in the current study (hypothesis 3). It seems that difficulty in controlling vocal fold vibrations has had an impact on the length of vocalic segments.

Vowel space area is reduced in post-surgery phases, in conformity with hypothesis 4. This is probably due more to alteration of the voice source than to potential perturbations in tongue body mobility. It is also plausible that the patient, in a postsurgery phase, is using gestural economy strategies. However, the latter conjectures must be verified with the help of articulatory data.

. Hypothesis 5 is also confirmed, as F0 values diminish in post-surgery phases, with spontaneous recuperation in patients without paralysis, in the second post-surgery phase. Recuperation is slower in patients with unilateral paralysis, however speech theray and time seem to have a positif effect on F0 values.

Hypothesis 6 was related to the benefits, in time, of speech therapy. Indeed, therapy has a positive effect on altered parameters in post-surgery phases. Patients without vocal fold paralysis recover, with therapy, as from the second post-surgery phase. Patients with unilateral paralysis, show slow recuperation as from the beginning of speech therapy; which is continually beneficial to the patient. This investigation is being pursued, extending data collection to a larger set of subjects, with or without paralysis. Articulatory data will also be collected in order to confirm some of our hypotheses.

Figure 4 .

 4 Figure 4. Vowel space areas for patient NPPG without unilateral paralysis

Figure 5 .

 5 Figure 5. Vowel space areas for patient PLPH with unilateral paralysis

Figure 7 .

 7 Figure 7. F0 values for patient PLPW with unilateral paralysis

Acknowledgments

This work was supported by the program « Perturbations et Réajustements : parole normale vs parole pathologique » of the Maison Inter-Universtaire des Sciences de l'Homme Alsace (MISHA), 2008-2012, and also by an ANR Grant "DOCVACIM" to the Phonetics Institut of Strasbourg / U.R. LiLPa, E.R. Parole et Cognition.