Jeffrey E Wieselthier
email: wieselthier@itd.nrl.navy.mil

Gam D Nguyen
email: nguyen@itd.nrl.navy.mil

Anthony Ephremides

MULTIPLE ACCESS FOR MULTIPLE DESTINATIONS IN AD HOC NETWORKS

In this paper, we address the medium-access control (MAC) problem in systems with multiple destinations, a problem that is crucial to the development of ad hoc networks. We focus on collision-resolution algorithms, specifically on First-Come, First-Serve (FCFS) splitting algorithms. We illustrate the difficulties that arise because of the interference between two populations of nodes that share a common channel, the most serious consequence of which is the possibility of deadlock that occurs when a well-known improvement for the single-destination case is used. Experimental results verify that the straightforward use of the FCFS algorithm in two-destination environments (without using this improvement) leads to deadlock-free performance. We show how performance can be improved by applying a time-division mechanism to separate interfering groups of nodes.

INTRODUCTION

One of the stumbling blocks to the implementation of ad hoc networks is the absence of medium-access control (MAC) schemes that are suitable for multihop networks. Throughout the long history of multiple-access research, virtually all studies of this problem have been limited to the case of a single destination. Although this framework was appropriate to address channel access problems that are characteristic of many cellular and satellite-based systems, the ad hoc networking environment presents challenges to the MAC problem that do not arise in the single-destination case. Several papers on routing in ad hoc networks include ideas and suggestions for handling the MAC issue, but do not examine the fundamental limitations of MAC in ad hoc networks (see e.g., [START_REF] Goldsmith | Design Challenges for Energy-Constrained Ad Hoc Wireless Networks[END_REF], which addresses many of the cross-layer interrelationships in ad hoc networks).

In this paper, we formulate the two-destination multiple-access problem, and present performance results relating to the throughput that can be achieved. We consider a scenario in which some nodes are within communication range of only one of these destinations, others are within range of the other, and some are within range of both (but each of their transmissions are intended for only one of the two destinations). We discuss the complications introduced by multiple destinations, most significantly the impact on the channel feedback process and the resulting impact on algorithm operation and performance. We focus on the fundamental problem of MAC in slotted systems; thus, we do not address carriersense multiple access (CSMA) -based approaches. 1 This work was supported by the Office of Naval Research. 1 Some issues associated with the use of IEEE 802.11 (which is

We consider two classes of channel access protocols, namely Slotted Aloha and splitting (collision-resolution) algorithms. In both cases, we illustrate the dependence of performance on the fraction of nodes that can communicate directly with both destinations. For Slotted Aloha, we evaluate the performance of uncontrolled systems. Our primary focus in this paper is on splitting algorithms for two or more destinations. We study the performance of the standard First-Come, First-Serve (FCFS) splitting algorithm [START_REF] Bertsekas | Data Networks[END_REF]. We illustrate the difficulties that are introduced by the need to resolve collisions among two populations that interfere with each other, and specifically show how an improvement that was developed for the single-destination case results in deadlock in the two-destination case. Finally, we propose improved versions that incorporate the information structures that are associated with the two-destination problem.

PROBLEM DEFINITION

Figure 1 shows a network with two destinations (D1 and D2) and a number of nodes that are within communication range of them (assuming the same fixed communication range for all nodes and the use of omnidirectional antennas). We define Group 1 to be the set of nodes within communication range of only D1; similarly Group 2 is the set of nodes within range of only D2. Nodes that are within range of both D1 and D2 (i.e., the intersection of the two circular regions) are said to be in Group 3. We use a simple interference model in which a node can cause interference only if it is within communication range, i.e., the combined effect of transmissions from multiple nodes outside the communication range is not considered. Nodes in Group 1 randomly generate packets for transmission to D1; similarly nodes in Group 2 randomly generate packets for transmission to D2. Nodes in Group based on CSMA with collision avoidance) in ad hoc networks are addressed in [START_REF] Xu | Does the IEEE 802.11 MAC Protocol Work Well in Multihop Wireless Ad Hoc Networks?[END_REF]. or D2 with equal probability; once a destination is chosen for a particular packet, the packet must be delivered to that destination (i.e., no reward is received for successful reception at the wrong destination).

Since packets from Group 3 reach both destinations, they are potentially a source of collisions at both the intended and non-intended destination. The challenge presented by the two-destination problem is to correctly model the impact of collisions caused by packets intended for another destination, and to develop channel-access protocols that are suitable for use in such environments.

We make the common assumption of an infinite population of unbuffered users in each group (e.g., as in [START_REF] Bertsekas | Data Networks[END_REF]). We define f = fraction of total traffic generated by nodes in Group 3.

Thus, when f = 0, Group 3 is empty and the collisionresolution process can proceed independently for the two destinations. At the other extreme, when f = 1, all packets reach both destinations.

The packets are classified into four types as follows, according to their intended destinations and group membership. Packets that are generated by Group 1's nodes are called type-1 packets; thus, type-1 packets are intended for D1 only, and never reach D2. Similarly, packets generated by Group 2's nodes are called type-2 packets. Recall that packets generated by Group 3's nodes are intended for either D1 or D2, and these packets will reach both D1 and D2. Packets generated by Group 3's nodes that are intended for D1 are called type-31 packets, while packets intended for D2 are called type-32 packets. These packets are generated according to a random (e.g., Poisson) process.

Each destination broadcasts ternary feedback information (i.e., success, collision, or empty slot) immediately following each slot. We assume that nodes in Group 1 and Group 3 receive D1's feedback, and that nodes in Group 2 and Group 3 receive D2's feedback. However, nodes in Group 1 (Group 2) do not receive feedback from D2 (D1).

SLOTTED ALOHA WITH TWO DESTINATIONS

We consider uncontrolled Slotted Aloha, where backlogged users retransmit with probability arbitrarily set to the fixed value of 0.01 in each slot. 2 When f = 0, the two regions are disjoint and the problem reduces to that of two independent Aloha systems. Our simulation results provide a maximum throughput of approximately 2/e = 0.736 packets per slot, or equivalently 1/e per destination, as expected. When f = 0.4, and thus 40% of the users are in Group 3, the maximum throughput per destination is approximately 0.26. Of course, such uncontrolled Aloha systems with infinite populations are unstable. They can be stabilized by choosing the retransmission probability as a function of estimates of the number of backlogged users.

However, better (and stable) performance can be achieved by using splitting algorithms, which have been studied in great detail for single-destination problems, but not previously for multiple destinations.

THE FIRST-COME FIRST-SERVE ALGORITHM WITH ONE DESTINATION

Splitting algorithms for contention resolution, which are inherently stable algorithms, provide a higher throughput than that which can be achieved using stabilized slotted Aloha. Splitting can be implemented in several ways, including coin toss, node ID, or time of arrival. Recently, the idea of splitting based on residual energy has been introduced [START_REF] Sagduyu | Energy-Efficient Collision Resolution in Wireless Ad-Hoc Networks[END_REF]. In this paper, we consider the use of First-Come First-Serve (FCFS) schemes. Splitting algorithms are based on the immediate availability of ternary feedback information (success, collision, or empty slot), transmitted by the destination, which describes the outcome of the slot.

We first consider the FCFS splitting algorithm for a single destination, which is discussed in detail in [START_REF] Bertsekas | Data Networks[END_REF]. Under this scheme, a collision-resolution epoch is initiated when a collision occurs. All packets that arrive within a specified "allocation interval" are transmitted in the first slot of this epoch. If there is a collision, this window is shortened, and the process is repeated.

The achievable throughput can be increased by means of two improvements [START_REF] Bertsekas | Data Networks[END_REF]. The first addresses the case of a collision that is followed immediately by an idle slot. Straightforward application of the FCFS algorithm guarantees that a collision will occur in the next slot. To avoid such wasting of a slot, the algorithm acts as if the collision had already occurred, and shortens the allocation interval accordingly. The second improvement (which has a somewhat smaller impact on achievable throughput) incorporates the right (i.e., later) half of the allocation interval involved in the original collision into subsequent collision resolution periods. The combined effect of using both improvements provides a maximum stable throughput of 0.4878.

It is noted in [START_REF] Bertsekas | Data Networks[END_REF] (as originally observed in [START_REF] Massey | Collision-Resolution Algorithms and Random-Access Communications[END_REF]), that use of the first improvement makes the algorithm vulnerable to errors in the feedback process. Specifically, consider the case in which an idle slot is erroneously interpreted as a collision. Since it is believed that a collision has occurred, the allocation interval is then split. However, since there are actually no packets in the allocation interval, the result is an empty slot. Use of the first improvement then results in further reduction of the empty allocation interval, a process that continues indefinitely, resulting in deadlock. (By contrast, it was observed in [START_REF] Massey | Collision-Resolution Algorithms and Random-Access Communications[END_REF] that, when the first improvement is not used, splitting algorithms are extremely insensitive to feedback errors.) In applications where such errors are rare, it may be appropriate to simply suspend the splitting of the interval after observing several such idle slots. However, it may be best not to use the first improvement at all when such errors are frequent. The case of multiple destinations creates a situation that is similar to that of frequent feedback errors, as we discuss in the next section.

THE FIRST-COME FIRST-SERVE ALGORITHM WITH TWO DESTINATIONS

We first consider the most straightforward implementation of the FCFS algorithm with two destinations, in which the first improvement is not used. For example, D1's neighbors are all the nodes in Group 1 and Group 3. D1 hears all of its neighbors' transmissions (even though some of the Group 3 transmissions are intended for D2). As in the single-destination case, D1's feedback for a slot in which none of its neighbors transmits is 0 (for an empty slot). If exactly one of its neighbors transmits, we must distinguish between two cases; the feedback is 1 (for a success) if the packet was, in fact, intended for D1; it is 0 if the packet is a Group 3 packet that was intended for D2. Whenever more than one of D1's neighbors transmits (regardless of their intended destinations), the feedback is e (for collision). In this simple implementation, nodes in Group 3 make use of the feedback transmitted by only their intended destination, even though they receive feedback from both destinations (which can potentially be used to improve performance). There is no coordination of the collision-resolution processes at the two destinations. Their allocation intervals are determined independently, as are the beginning and end points of the collision resolution intervals.

The difference between this problem and that of a single destination results from the interference (collisions) caused by packets in Group 3 that are intended for D2. We noted above that splitting algorithms (for a single destination) are extremely insensitive to feedback errors when the first improvement is not used. We have observed in our simulations of the two-destination problem that the straightforward implementation of the FCFS algorithm functions correctly (albeit at reduced throughput levels) despite the introduction of collisions with packets intended for the other destination. Before discussing the throughput levels that are achievable, we address the deadlock that can arise if the first improvement is used in the two-destination problem, which is similar to that which results from feedback errors, as discussed above. We then discuss a simple technique that eliminates the deadlock mechanism. Figure 2 illustrates an example that leads to deadlock if the first improvement is used. Suppose that we start with an empty system, and that three packets (A, B, and C) are generated during time slot 0. In this case, to simplify the figure, we implicitly assume that the initial value of the allocation interval for both D1 and D2 is equal to a slot length, although there is in general no direct relationship between allocation interval (a continuous variable) and slot length (which is fixed). What is important here is that A and B arrive during the first half of the allocation interval corresponding to D1 and that C arrives during the second half of the allocation interval corresponding to D2 (where in both cases, the initial allocation period for both D1 and D2 coincides with slot 0).

In our example, both A and B are from Group 3 and are intended for D1, while C is from Group 2 and is intended for D2. All three packets attempt transmission in slot 1, and both D1 and D2 experience collisions (note that all three transmissions reach D2). As a result of these collisions, both D1 and D2 will cut their allocation intervals in half; therefore, packets that arrived during the first half of slot 0 (namely, A and B, but not C) will transmit during slot 2. Again, a collision will be observed at both D1 and D2. D1's allocation interval will then be halved again, which may permit the successful transmission of A in slot 3 and B in slot 4 (depending on their actual arrival times within the first half of slot 0). However, D2's allocation interval will also be halved, with the first half (quarter of the original allocation interval) being assigned to slot 3. Upon finding this interval to be empty, and thus assuming that a collision is guaranteed in the next slot, the second half will be split (the first half of which will again be found to be empty), and so on. Therefore, the allocation interval for D2 will never extend beyond the third quarter of slot 0. Consequently, transmission of packet C, and thereafter all packets in Group 2, will be suspended indefinitely. This type of deadlock is not a rare event, and the high likelihood of its occurrence, even in networks with no feedback errors, suggests that it is not advisable to use the first improvement in multiple-destination scenarios, unless a means is developed to eliminate the possibility of such deadlocks.

D1 D2 Group 1 Group 2 Group 3 A B C A B C 0 1 2 3 slots ...
One possible approach to eliminate the complications created by the two destinations is to impose a timedivision structure that creates a time orthogonality among groups of nodes that interfere with each others' feedback process. For example, consider Fig. 1. The nodes in Group 1 and Group 2 can transmit simultaneously without any interference between the feedback of the two groups. The nodes in Group 3 would interfere with those in both Group 1 and Group 2; moreover, Group 3 nodes transmitting to D1 would interfere with those that transmit to D2. Let us assume that our goal is to determine an appropriate slot assignment that maximizes the throughput that can be achieved. 3 This situation suggests that three classes of slots be allocated proportionately to the packet rate in each class, namely: fraction of slots allocated to Group 1 and Group 2 nodes:

1-f 1+ f fraction of slots allocated to Group 3 nodes transmitting to D1: f 1+ f fraction of the slots allocated to Group 3 nodes transmitting to D2: f 1+ f Such a time-division structure would eliminate deadlock situations such as the one described above, and would therefore permit the use of the first improvement (assuming that feedback errors caused by noise are negligible).

Let us consider the problem of determining the maximum throughput per destination that can be supported in the network of Fig. 1. We assume a symmetrical system with an infinite number of nodes, in which the combined arrival rate of Group 1 nodes is equal to that of Group 2 nodes. We also specify the fraction f of the overall arrival rate of nodes in Group 3, half of which is destined for D1 and half for D2. It can be shown that the maximum achievable throughput per destination is [START_REF] Nguyen | Contention-Resolution Algorithms for Multiple Destinations in Wireless Networks[END_REF] S norm = S 1+ f where S is the maximum throughput achievable by the FCFS splitting algorithm (i.e., 0.4878 if both improvements are used). Another slot allocation can yield throughput that is bounded above by (1-f/2)S [START_REF] Nguyen | Contention-Resolution Algorithms for Multiple Destinations in Wireless Networks[END_REF].

The upper curve in Fig. 3 shows the maximum throughput per destination that can be achieved by allocating the time slots according to the fractions discussed above, as a function of f. The lower curve represents the use of the FCFS splitting algorithm independently at the two nodes, where only the second improvement is used (to avoid the deadlock situation discussed above). We observe that the orthogonal time division of transmissions can provide better performance, largely because it permits use of the first improvement. We also observe that reasonably good performance can be achieved by a very simple implementation, in which the FCFS algorithm is implemented without any effort to coordinate the collision-resolution process at the two destinations.

SUMMARY AND CONCLUSIONS

In this paper, we have addressed the MAC problem for multiple destinations, and have described the issues associated with using time-of-arrival-based splitting algorithms in such environments. We have demonstrated the susceptibility of the FCFS splitting algorithm to deadlock when a well-known throughput-increasing improvement is incorporated. However, we have also demonstrated, via simulation, that a straightforward implementation (without this improvement) is not subject to this deadlock. Furthermore, we have shown that the benefits of the improvement can be exploited if time division is used to partition the groups of users so that they do not interfere with each other.

Our future work will address techniques to improve the logic of the collision resolution algorithm, e.g., the nodes in Group 3 may be able to exploit the feedback information from both nodes, rather than only the node to which they are transmitting. We will also address the case of multiple destinations. The extension of our timedivision approach to multiple destinations is related to the link-activation problem (see, e.g., [START_REF] Baker | The Design and Simulation of a Mobile Radio Network with Distributed Control[END_REF], [START_REF] Ephremides | A Design Concept for Reliable Mobile Radio Networks with Frequency Hopping Signaling[END_REF], [START_REF] Ephremides | Scheduling Broadcasts in Multihop Radio Networks[END_REF]), in which non-interfering nodes (i.e., nodes that are sufficiently far from each other) are permitted to transmit simultaneously.

3 Fig. 1 -

 31 Fig. 1 -Multiple access in two overlapping regions.

Fig. 2 -

 2 Fig. 2 -Example of deadlock when the first improvement is used in a two-destination network.

Fig. 3 -

 3 Fig. 3 -Throughput per destination for two-destination collision resolution, as a function of the fraction of nodes in the intersecting region.

This assumption deviates slightly from the true classic Aloha scheme; usually the retransmission probability is chosen to either stabilize or control Aloha's behavior. However, here, the attention to Aloha is solely for establishing benchmark behavior and is not of interest per se. Thus the assumed constant retransmission probability represents a minor inconsequential deviation.

A variety of similar problems can also be considered. For example, if the fraction of total traffic associated with nodes in each Group is specified, it would be appropriate to determine the slot allocation that provides the maximum achievable throughput, as well as the throughput that can actually be achieved.


A predictive robust cascade position-torque control strategy for Pneumatic Artificial Muscles Lotfi Chikh † , Philippe Poignet ‡ , Franc ¸ois Pierrot ‡ and Micaël Michelin † Abstract-A robust cascade strategy is proposed and tested on an electropneumatic testbed for parallel robotic applications. It is applied on Pneumatic Artificial Muscles (PAMs). Nonlinear models are developed and presented. By specifying the pressure average between the two muscles, it is possible to control the torque by controlling the pressure in each muscle. A constrained LMI based H∞ controller is synthesized for the pressure inner loop. A computed torque is applied on the torque loop and a Generalized Predictive Controller (GPC) controller is then applied to control position of the muscles. Experimental results show the feasibility of the control strategy and good performances in terms of robustness and dynamic tracking.

Index Terms-GPC, cascade control, pneumatic actuators, modeling, LMI, feedback linearization, robust control, H∞, PAMs.

I. INTRODUCTION

The paper is motivated by pick-and-place parallel robotic applications: the study of 1-dof pneumatic actuation system is a preliminary stage which enables us to evaluate performances of Pneumatic Artificial Muscles (PAMs). Parallel robots have a lot of advantages which have made their success. Particularly, they are very fast as the actuators are transferred to the rigid frame reducing considerably the inertia of the moving links and increasing speeds and accelerations of the end effector. For instance, the fastest industrial robot in the world -the Quattro-has a parallel structure and can reach an acceleration of 15g [START_REF]Adept[END_REF]. Recently, the Par2 robot which is not industrialized yet has reached 43g while keeping a low tracking error [START_REF] Pierrot | Above 40g acceleration for pick-and-place with a new 2-dof pkm[END_REF]. However, a major obstacle for parallel robots expansion is their expensive price due to the cost of the embedded motors which they use. In this context, the inspection of other actuation type such as pneumatic actuators is a matter of fact as they are cheap actuators with low maintenance costs, and good force/weight ratio. As the big obstacle of industrial use of pneumatic actuation in robotics is difficulty in their control, this paper focuses on this particular point.

Numerous techniques have been studied in literature and among them a large part concerns robust control. Robust controllers are necessary to deal with uncertainties and to ensure high precision positioning. They have been applied as a nonlinear feedback using sliding modes by [START_REF] Van Damme | Sliding mode control of a 2-dof planar pneumatic manipulator[END_REF] for a planar 2 DOF manipulator. In [START_REF] Xiaocong | Integrated direct/indirect adaptive robust posture trajectory tracking control of a parallel manipulator driven by pneumatic muscles[END_REF], authors explored adaptive control techniques for parallel robotic applications. One of the problems of nonlinear controllers is difficulty in the synthesis of the control law and the high computational amount. That is why some authors have tested linear robust control techniques such as [START_REF] Osuka | H infinity control of a certain nonlinear actuator[END_REF] where PAMs are tested in an antagonist manner. In [START_REF] Osuka | H infinity control of a certain nonlinear actuator[END_REF], an H ∞ controller is synthesized basing on a linear identified nominal system which greatly limits the control performances. This explains why [START_REF] Kimura | H infinity control with minor feedback for a pneumatic actuator system[END_REF] has synthesized an H ∞ controller after application of feedback linearizing [START_REF] Isidori | Nonlinear Control Systems[END_REF] [8] of the nonlinear system. However, only one muscle is controlled and the antagonist case is not studied.

In this paper, a novel robust cascade strategy which combines an outer position predictive controller and an inner torque controller is proposed. The torque controller is based on the H ∞ pressure control in each muscle after specifying an average pressure into the muscles and a desired torque reference. This desired torque reference is given by the outer position loop. Cascade control has been studied in case of PAMs by [START_REF] Schröder | Improved control of a humanoid arm driven by pneumatic actuators[END_REF] and [START_REF] Hildebrandt | Cascaded control concept of a robot with two degrees of freedom driven by four artificial pneumatic muscle actuators[END_REF]. In this two studies, only simple controllers are applied. In our case, we extend the cascade control concept to robust controllers and predictive ones. The different experimental tests show the advantages of this approach, in terms of good time responses and ability of tracking of extended range references (up to ±25 • ). Another motivation of predictive control is that in a precedent study on pneumatic cylinders [START_REF] Chikh | A mixed GPC-H∞ robust cascade position-pressure control strategy for electropneumatic cylinders[END_REF], it was shown that predictive cascade concept can improve robustness performances comparing to state of art control techniques such as sliding modes. The control synthesis is simplified as system is linearized using feedback linearizing techniques [START_REF] Michel | Flatness and defect of non-linear systems: introductory theory and examples[END_REF] [START_REF] Isidori | Nonlinear Control Systems[END_REF]. A Generalized predictive controller (GPC) is synthesized for the outer loop and as the system is linear, an explicit solution exists for the predictive optimization problem. The obtained controller is an easy-to-implement one. As far as we know, only two experimental studies of GPC -on pneumatic cylinders and not on PAMs-have been carried out so far [START_REF] Song | Improved control of a pneumatic actuator pulsed with pwm[END_REF] [START_REF] Norgaard | Intelligent predictive control of nonlinear processes using neural networks[END_REF]. In [START_REF] Song | Improved control of a pneumatic actuator pulsed with pwm[END_REF], the model used is a linear one which limits greatly the performances of the controller. In [START_REF] Norgaard | Intelligent predictive control of nonlinear processes using neural networks[END_REF], the authors used a model estimation based on neural network theory [START_REF] Norgaard | Intelligent predictive control of nonlinear processes using neural networks[END_REF]. No application of GPC basing on an explicit nonlinear model has been found in literature. An intuitive motivation of using predictive theory in pick and place robotic applications is that the trajectories are determined a priori and therefore, future trajectory is known. By specifying the pressure average between the two muscles, it is possible to control the torque by controlling pressure in each muscle. A constrained LMI based H ∞ controller is synthesized for the pressure inner loop. It is based on LMI optimization [START_REF] Gahinet | The lmi control toolbox[END_REF] [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]. In addition to the classical advantages of H ∞ [START_REF] Duc | Commande Hinf et mu-analyse un outil pour la robustesse[END_REF] [17] control in terms of robustness, disturbance rejection, systematic synthesis of MIMO controllers and powerful combination of both frequency domain synthesis and state space synthesis, the LMI approach enables the addition of constraints in an intuitive manner. Therefore, pole placement constraints have been added to the H ∞ performances in order to have a better control of the transitory temporal behavior of the pressure controller.

This paper is organized as follows. Section II presents the versatile electropneumatic testbed for robotic applications and its nonlinear modeling. Section III deals with the control of the muscles. After introducing the feedback linearizing equations, both of GPC controller background and LMI based H ∞ multi-objective approach are introduced. The cascade strategy which combines these two control techniques is then introduced. Finally in section IV, the robust cascade strategy is implemented experimentally and different control tests, including robustness tests, are presented and deeply analyzed.

II. EXPERIMENTAL PNEUMATIC TESTBED AND NONLINEAR MODELING

The versatile electropneumatic setup is presented in Fig. 1. It includes two 5/3-way proportional valves1 and three kinds of actuators. The first one is a standard double acting cylinder which is very widespread in industry and one of the cheapest actuators. The second one is the rodless cylinder which has the advantage of being symmetric, that is to say, the maximum force provided in one direction equals the one in the inverse direction. Finally, Pneumatic Artificial Muscles2 (PAMs) which are deeply studied in this paper. Both of the two muscles are connected to a moving arm which may represent a robot arm. In the long run, our aim is to design efficient parallel robot which rely on pneumatic actuation with robust control: the moving arm of this test bed is very similar to arms used for various parallel robots such as Delta or Quattro robots. At its extremity, different masses can be attached for robustness tests of the controllers. Three types of sensors are used; a high resolution incremental encoder, two pressure sensors and two force sensors. The real time prototyping environment is xPC Target T M from Mathworks. In the sequel, we present the nonlinear model of the pneumatic muscles. Every electropneumatic positioning device includes an actuation element (the pneumatic cylinder), a command device (the valve), and a mechanical part and position, pressure and/or force sensors. A schematic representation of the electropneumatic system is given in Fig. 2. Supply pressure p s is supposed to be constant. p 0 denotes the atmospheric pressure. It is supposed that any variation of the muscle's volume or pressure can be described by the polytropic gas law [START_REF] Hildebrandt | A cascaded tracking control concept for pneumatic muscle actuators[END_REF]:

p 1 V γ 1 = p 2 V γ 2 (1)
Where p i is the pressure in one of the muscle chamber (indexes 1 and 2 are related to two pressure states) and V i is the volume in the muscle chamber, γ is the polytropic constant. The ideal gas equation describes the dependency of the gas mass:

m = pV rT ( 2
)
where m is the gas mass inside the cylinder chamber or the muscle, T is the air temperature which is considered to be equal to the atmospheric temperature and r is the specific gas constant. Therefore, combining ( 1) and ( 2) leads to the pressure dynamic expression:

dp dt = γ V (s) [rT q m (u, p) -p dV ds ṡ] (3)
where u represents the input voltage of the valve, s is the position of the piston in case of the cylinders. q m (u, p) represents the mass flow rate ( dm dt = q m (u, p)). Dynamic effects of the underlying position controller for the valveslide stroke are neglected. These hypothesis justify why the mass flow is a function of the input voltage and the pressure in the actuator chamber. The experimental curve representing the volume/stroke dependency is represented on figure 3. According to this characteristic, the volume/stroke model has been approximated by a third order polynomial using least squares:

V (s) = 3 i=0 b i s i (4)
The second term of ( 3) is now completely defined; the first term requires modeling of the valve. The valve model has been approximated -after identification-by the following expression (see [START_REF] Belgharbi | Analytical model of the flow stage of a pneumatic servodistributor for simulation and nonlinear control[END_REF] for details):

q m (u, p) = ϕ(p) + ψ(p)u (5)
ϕ and ψ defines 5 th -order polynomials with respect to p. For cascade control experiments, it is mandatory to determine the force/pressure dependency. For the muscles, force depends on pressure and also on contraction. The force is maximum when there is no contraction. As the contraction increases the generated force decreases in a non linear way. The experimental measured force characteristic is represented on Fig. 4. It is approximated by the following function:

F m (p, s) = (p -p 0 ) 5 i=0 c i s i (6)

III. CONTROL OF THE PNEUMATIC ARTIFICIAL MUSCLES

In this section, we first present the feedback linearization equations for the PAMs. The approach handled here can be justified by differential geometry concepts [START_REF] Isidori | Nonlinear Control Systems[END_REF] and differential flatness theory [START_REF] Michel | Flatness and defect of non-linear systems: introductory theory and examples[END_REF]. Once the linearized system obtained, two advanced control techniques are presented: GPC and LMI based constrained H ∞ .

A. Feedback linearization equations

It can be easily proved -in case of mass flow rate rate m g as an input and the pressure difference as an output-that differential flatness criteria is satisfied and that the system is completely linearizable (using differential geometry concepts [START_REF] Isidori | Nonlinear Control Systems[END_REF]). The pressure dynamic in each muscle is given by:

dp dt = γrT V (s) ṁg - γ V (s) dV ds ṡp (7)
If we take,

ṁg = V γrT (u aux + γ V dV ds ṡp) (8)
This leads to the following linearized system:

ṗ = u aux (9)
The linearized system is then a first integrator. The pressure controller which is synthesized basing on this first integrator is presented in the next section.

B. Multi-objective output feedback pressure control via LMI optimization

An LMI is any constraint of the form:

A(x) = A 0 + x 1 A 1 + . . . + x N A N < 0 (10)
where A 0 . . . A N are given symmetric matrices and x T = (x 1 . . . x N ) is the vector of unknown variables. In [START_REF] Hildebrandt | Cascaded control concept of a robot with two degrees of freedom driven by four artificial pneumatic muscle actuators[END_REF], the symbol < refers to negative definite 3 . One way of tuning simultaneously the H ∞ performance and transient behavior is to combine H ∞ and pole placement objectives using LMI optimization techniques. Poles are clustered in regions which can be expressed in terms of LMIs. The class of LMI region defined below has been introduced for the first time by [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]. It turns out to be suitable for LMI-based synthesis.

Def. LMI Regions. A subset D of the complex plane is called an LMI region if there exists a symmetric matrix

α = [α kl ] ∈ R m×m and a matrix β = [β kl ] ∈ R m×m such that: D = {z ∈ C : f D (z) < 0} (11)
with:

f D (z) := α + zβ + zβ T = [α kl + β kl z + β lk z] 1≤k,l≤m.
For instance, we use a disk LMI region centered at (-q, 0) with radius r. It is defined below [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]:

f D (z) = -r q + z q + z -r < 0 (12)
The constrained H ∞ problem under consideration can be stated as follows [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]. Given an LTI plant: T is composed by the controlled output y and the control signal u. In our case, and as defined latter in (37), y is the pressure difference between the chambers and u is the input voltage.

ẋ(t) = Ax(t) + B 1 ω(t) + B 2 u(t) e(t) = C 1 x(t) + D 11 ω(t) + D 12 u(t) y(t) = C 2 x(t) + D 21 ω(t) + D 22 u(t) (13
The controller transfer function is denoted K(s) and can be represented in the following state-space form by:

ẋK (t) = A K x K (t) + B K y(t) u(t) = C K x K (t) + D K y(t) (14)
Then, T wz (s

) = D cl + C cl (sI -A cl ) -1 B cl with 4 : A cl = A + B 2 D K C 2 B 2 C K B K C 2 A K , B cl = B 1 + B 2 D K D 21 B K D 21 , C cl = (C 1 + D 12 D K C 2 , D 12 C K ) and D cl = D 11 + D 12 D K D 21 .
We first examine each specification separately. It is shown in [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF] that the pole placement constraint is satisfied if and only if there exists X D > 0 such that:

[α kl X D + β kl A cl X D + β lk X D A T cl ] 1≤k,l≤m < 0 (15)
Meanwhile, the H ∞ constraint is expressed in terms of LMIs:

  A cl X ∞ + X ∞ A T cl B cl X ∞ C T cl B T cl -γI D T cl C cl X ∞ D cl -γI   < 0 (16)
Then the problem formulation of H ∞ synthesis with pole placement -assuming that the same Lyapunov matrix X > 0 is required -is:

Find X > 0 and a controller K(s) > 0 ≡ Ω K that satisfy (15) and ( 16)

with X = X D = X ∞ (17)
The controller matrix is denoted by:

Ω K = A K B K C K D K (18)
The difficulty in output feedback is that relations (15) and ( 16) involve nonlinear terms of the form BΩ K CX. This means that problem formulation is not convex and then, can not be handled by LMIs. Chilali and Gahinet [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF] solved this problem by taking the following change of variables of the controller:

B K = N B K + SB 2 D K C K = C K M T + D K C 2 R A K = N A K M T + N B K C 2 R + SB2C K M T + S(A + B 2 D K C 2 )R (19)
where R, S, N and M correspond to the following partition of X and it inverse as

X = R M M T U , X -1 = S N N T V R ∈ R n×n , S ∈ R n×n (20)
The proposed procedure is summarized in the following theorem [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]: Theorem1: Let D be an arbitrary LMI region contained in the open left-half plane and let [START_REF] Chikh | A mixed GPC-H∞ robust cascade position-pressure control strategy for electropneumatic cylinders[END_REF] be its characteristic function.

Then, the modified problem ( 17) is solvable if and only if the following system of LMIs is feasible. Find R = R T ∈ R n×n , S = S T ∈ R n×n , and matrices A K , B K ,C K and D K such that

R I I S > 0 (21)
α kl R I I S + β kl Φ + β lk φ T k,l < 0 (22)
Ψ 11 Ψ T 21 Ψ 21 Ψ 22 < 0 (23)
with the shorthand notation

Φ := AR + B 2 C K A + B 2 D K C 2 A K SA + B K C 2 (24)
Ψ 11 , Ψ 12 and Ψ 22 terms are detailed in [START_REF] Chilali | H-infinity design with pole placement constraints: an lmi approach[END_REF]. Given any solution to this LMI system:

• Compute via Singular Values Decomposition (SVD) a full-rank factorization M N T = I -RS of the matrix I -RS (M and N are then square invertible) • Solve the system of linear equations [START_REF] Belgharbi | Analytical model of the flow stage of a pneumatic servodistributor for simulation and nonlinear control[END_REF] for B K , C K and A K (in this order).

• Set K(s) = D K + C K (sI -A K ) -1 B K .
Then K(s) is an nth order controller that places the closedloop poles in D and such that T wz ∞ < γ.

C. GPC outer position Controller

The GPC algorithm is based on a CARIMA model which is given by:

A(z -1 )y(t) = z -d B(z -1 )u(t -1) + C(z -1 )e(t) (25)
where u(t) and y(t) are respectively the control and output sequences of the plant and e(t) is a zero-mean white noise.

A, B and C are polynomials of the backward shift operator z -1 . They are given by;

A(z -1 ) = 1 + a 1 z -1 + a 2 z -2 + ... + a na z -na , B(z -1 ) = b 0 + b 1 z -1 + b 2 z -2 + ... + b nb z -nb and C(z -1 ) = 1 + c 1 z -1 + c 2 z -2 + ... + c nc z -nc and = 1 -z -1 .
For simplicity, it is admitted that C -1 equals 1. It is important to mention that most SISO plants can be described by a CARIMA model after linearization. The GPC algorithm consists in applying a control sequence that minimizes a multistage cost function of the form:

J(N 1 , N 2 , N u ) = N2 j=N1 δ(j)[ŷ(t + j|t) -w(t + j)] 2 + N2 j=1 λ(j)[ u(t + j -1)] 2 (26)
where ŷ(t + j|t) is an optimum j step ahead prediction of the system output on data up to time t, N 1 and N 2 are the minimum and maximum costing horizons, N u is the control horizon, δ(j) and λ(j) are weighting sequences and w(t + j) is the future reference trajectory. The minimization of cost function leads to a future control sequence u(t), u(t + 1), ... where the output y(t + j) is close to w(t + j). Therefore, in order to optimize cost function, the best optimal prediction of y(t + j) (for N 1 ≤ j ≤ N 2 ) has to be determined. This needs the introduction of the following Diophantine equation:

1 = A j (z -1 ) Ã(z -1 ) + z -j F j (z -1 ) (27)
with à = A(z -1 ) and polynomials E j and F j are uniquely defined with degrees j -1 and na respectively. is defined as = 1 -z -1 By multiplying (25) by E j (z -1 )z j and considering (27), we obtain:

y(t + j) =F j (z -1 )y(t) + E j (z -1 )B(z -1 ) u(t + j -d -1) + E j (z -1 )e(t + j) (28)
Since the noise terms in (28) are all in the future (this is because degree of polynomial E j (z -1 ) = j -1), the best prediction of y(t + j) is:

ŷ(t + j|t) = G j (z -1 ) u(t + j -d -1) + F j (z -1 )y(t) (29)
where G j (z -1 ) = E j (z -1 )B(z -1 ) Polynomials E j and F j can merely be obtained recursively (demonstration can be found for instance in [START_REF] Camacho | Model Predicitve Control. Second Edition[END_REF]).

In the future, it will be referred only to N = N 2 = N u as the prediction horizon. N 1 is chosen equal to 0. Let's consider the following set of j ahead optimal predictions:

ŷ(t + d + 1|t) = G d+1 u(t) + F d+1 y(t) . . . ŷ(t + d + N |t) = G d+N u(t + N -1) + F d+N y(t) (30)
It can be written in the following compact form:

y = G u + F(z -1 )y(t) + G (z -1 ) u(t -1) (31)
where terms y, u, G and F can be defined in [START_REF] Camacho | Constrained generalized predictive control[END_REF]. Equation (31) can be rewritten in this form:

y = Gu + f (32)
Where f refers to the last two terms in (31) which only depend on the past. Now, we are able to rewrite (26) as:

J = (Gu + f -w) T (Gu + f -w) + λu T u (33)
where w = [w(t + d + 1), w(t + d + 1) . . . w(t + d + N )] T Equation (33) can be written as:

J = 1 2 u T Hu + b T u + f 0 (34) with H = 2(G T G + λI), b T = 2(f -w) T G and f 0 = (f -w) T (f -w)
Therefore, the minimum of J can simply be found by making the gradient of J equal to zero, which leads to:

u = -H -1 b = (G T G + λI) -1 G T (w -f ) (35)
Since the control signal that is actually sent to the process is the first element of vector u (receding strategy), it is given by:

u(t) = K(w -f ) ( 36
)
where K represents the first element of matrix (G T G + λI) -1 G T . Contrary to conventional controllers, predictive ones depend only on future errors and not past ones.

D. Position/torque control strategy for the pneumatic cylinders

The cascade concept for the muscles consists in an outer loop which controls the position and an inner one for the torque. As far as we know, it has been studied in literature in case of the muscles only by [START_REF] Hildebrandt | Cascaded control concept of a robot with two degrees of freedom driven by four artificial pneumatic muscle actuators[END_REF] and [START_REF] Schröder | Improved control of a humanoid arm driven by pneumatic actuators[END_REF]. In our case, we are investigating first the possibility of applying predictive control in such a scheme and the improvement that can be brought by such robust control techniques like GPC and constrained H ∞ . The presented approach is described on Fig. 6. As presented before, the linearized system in case of pressure as an output is an integrator. Therefore, equations of the system can be stated as follows:

ẋ = (0)x + 1 0 b v + (1)u y u = -1 0 x + 1 0 0 0 b v + 0 1 u y = (1)x + 0 1 b v + (0)u (37)
where ω = (b v) T and e = (y u) T have been defined before. Classical assumptions for solving H ∞ problem can easily be verified 5 . Using Newton second law yields to the following equation:

τ m = I mov α + M mov gdcos(α) ( 38
)
where α is the rotation angle, I mov represents the inertia of the moving link (in our case it can be the arm with the attached load), d is the distance between the center of rotation and center of mass of moving link, g is the gravity and τ m the torque generated by the two antagonist muscles. The computed torque [START_REF] Schilling | Fundamentals of robotics, analysis and control[END_REF] consists in the application of the following control signal:

τ m = I mov v + M mov gdcos(α) (39)
where v represents the new control input. Therefore, the obtained linearized system is a double integrator:

α = v (40)
The GPC controller will be then synthesized based on a double integrator.

For the inner control loop, the objective is to track the reference torque which is set by the outer position controller.

In order to set the pressure reference in each muscle, we introduce the pressure average which is defined by:

P m = (P m1 + P m2 )/2 (41)
Combining this equation with:

τ m = F m1 (s 1 , p 1 )s 1 (α) -F m2 (s 2 , p 2 )s 2 (α) (42)
Therefore, using 41, 42 and s i = αR, where R is the radius of the motor between the two muscles. The reference pressure in the second muscle can be calculated by the following relation:

P 2 ref = p m ref - T ref 2R 5 i=0 c i s i (43) 5 (A, B 2 ) controllable and (C 2 , A) observable

IV. EXPERIMENTAL RESULTS

The cascade GPC/H ∞ is implemented on the PAMs. Nominal tests and robustness are presented respectively. Some typical tests: On Fig. 7, a test of a square signal with an amplitude of +/ -10 • and a period of 8s is performed. In this test, the arm is not attached yet. The rising time equals 0.113s and settling time equals 0.148s. The overshoot can be accepted because the controller input will never a step of this size in industrial tasks. For the H ∞ controller, the LMI region for each muscle is a disk centered at (-100, 0) with a radius of 15. The choice of these values is motivated by the fact that big closed pole location values will induce an oscillatory control signal with chattering effect. Small pole values will not insure good tracking results. A tradeoff has to be found. For the synthesis of the position outer controller, the linearized system after computer torque is a double integrator. Therefore, and for a sampling time of 1ms, polynomials A and B of the CARIMA system introduced in Eq. 25 are given by A(z -1 ) = 1-2z -1 +z -2 and B(z -1 ) = 5•10 -7 (1+z -1 ). λ is chosen equal to 0.09. The average pressure has been chosen equal to 4bar. Indeed, there is a relation between the average pressure and the stiffness and consequently, this can have good or bad effects concerning the control performances. In order to evaluate this dependency, two tests with different average pressures and the same control tests are handled.

Results are represented in Fig. 8. As related on table I, the Average pressure influence on control results

Reference trajectory

Low average pressure (2bar) High average pressure (4bar) Fig. 8. Studying the effect of the average pressure influence performances obtained with a high pressure average are better. First, both negative and positive overshoots are reduced in case of p m = 4bar. The same remark can be done for the settling time which is smaller for a high pressure average. The difference is about 0.3s which is considerable. These results can be explained by the fact that a high pressure average increases the stiffness of the muscles which can bring an important improvement to the control performances. However, we are not able to claim that the higher stiffness is the better the performances are. Indeed some control tests with dynamic reference trajectories have given better results in case of low pressure averages. This can be explained by the fact that in the dynamic case, lower stiffness becomes an advantage for the system in terms of fast reactivity to reference changes. In Fig. 9, a dynamic test with the muscles for a fast sine reference is done. The period equals T = 2s and the amplitude is between +/-10 • . The maximum error equals 2.24 • during the switching time and less than 1 • during the dynamic tracking. On Fig. 10, the arm is attached and a high angle range test is done.. Indeed, one of the drawbacks of muscles is their small contraction capacity which does not exceed 20% of the length. In this +/ -25 • test, the tracking is insured but the problem for tests with arm is the appearance of vibrations that disturb the performances. Maximum error equals 3.25 • . This elasticity problem of muscles that limits the tracking performances especially for high angle range tests with loads is our big challenge up to now. An asymmetry is observed in the control results, possible explanations can be the friction influence, calibration problems with the valve, gravity effects or a backlash error. These are also some directions that will be investigated in the next months.

Robustness tests: In order to study the robustness of the developed control strategy, the following test is handled. The objective is to study the influence of the time delay for ap- plications where the valves are located far from the actuators. Indeed, one of the advantages of pneumatic actuation is the possible use in explosive environments. This is not true unless the electric auxiliary components -valves and sensors-are mounted far from the actuator. For the test shown in Fig. 11, a 10m long flexible tube is added between valve and cylinder. Actually, in parallel robotic applications, the valves will be mounted close to the actuators and this problem of time delay will not be crucial. This explains why time delay has not been included in the model. The objective of this test is just to study the robustness performances of the controller.

Results show effectively a deterioration in performances but stability is maintained and tracking error is always less than 1 • . Another robustness mass variation test is handled. The error signal is represented in figure 12. As a conclusion to this series of tests, following observations comparing to state-of-the-art former results can be made. First, the robust cascade approach brings good results in terms of robustness which is a very important requirement in industry. It is possible to track high amplitudes reference signals of (+/-25 • ) which is substantial regarding to the limited muscle contraction abilities. Finally, another inherent advantage of the proposed approach is that it enables a bigger flexibility in the adjusting of the controller gains as it is possible to act either on the position controller, pressure one and also on the pressure average reference.

V. CONCLUSION

A robust cascade predictive/H ∞ control strategy for PAMs is introduced. The goal is to test the muscles for parallel robotics applications. Experimental results have shown good tracking performances and good precision with a static error less of 0.04 • . Robustness tests have been also handled and the tracking is still insured when a 10m distance exists between valve and actuators. Extended angular displacement tests are experimented for the muscles and give satisfactory results. Future work will concern the elasticity problem which appears for fast tests with 5-10kg loads. A possible direction should be a deeper study of the average pressure influence on performances in order to design of optimal average pressure references which insure the best accuracy and robustness performances. Another possibility is to take this elasticity into account explicitly in the model.

Fig. 1 .

 1 Fig. 1. Versatile electropneumatic setup for robotic application. Fig. 2. Schematic representation of the experimental setup.

Fig. 2 .

 2 Fig. 1. Versatile electropneumatic setup for robotic application. Fig. 2. Schematic representation of the experimental setup.

  Length of the muscle (in m) Volume of the muscle (in m 3 ) Volume/stroke dependency

Fig. 3 .

 3 Fig. 3. Volume Stroke dependency

Fig. 4 .

 4 Fig. 4. MAS20 force characteristic

Fig. 5 .

 5 Fig. 5. System representation and notations

Fig. 6 .

 6 Fig.6. General block diagram of the cascade Position/torque strategy for the muscles. The inner loop consists in the control of the torque generated by the two muscles. This is done by controlling the pressure in each muscle. An LMI based H∞ controller is implemented. The outer loop consists in the computed torque, a position generalized predictive controller is then implemented.

Fig. 7 .

 7 Fig. 7. Pulse response of the muscles in case of a cascade Position/torque control: No load case

Fig. 9 .Fig. 10 .

 910 Fig. 9. Sine response of the muscles in case of a cascade Position/torque control

Fig. 11 .

 11 Fig. 11. Robustness position test; valve and pressure sensors have been transferred at 10m far from the actuator, vibrations appear but stability is preserved and tracking is still correct with a tracking error less than 1 •

Fig. 12 .

 12 Fig. 12. Robustness mass variation test: a 1kg and a 5kg loads are added at t = 10s.

MPYE-5-1/4-010-B from FESTO

PAMs: FESTO MAS -20-450N-AA-MC-K

we assume that the D 22 = 0, this assumption considerably simplifies the formulas. (Note that it is always possible to remove the D 22 term by a mere change of variables)

ACKNOWLEDGMENT

The authors would like to thank the Fatronik Tecnalia Foundation for its support with the experimental setup.


Lotfi Chikh
email: lchikh@fatronik.cometmmichelin@fatronik.com

Philippe Poignet
email: poignet@lirmm.fretpierrot@lirmm.fr

Micaël Michelin

François Pierrot

Commande predictive robuste des vérins

Keywords: actionneurs pneumatiques, GPC, LMI, H∞, linéarisation

teaching and research institutions in France or abroad, or from public or private research centers.

I. Introduction

Cette étude rentre dans le cadre de la commande des actionneurs pneumatiques dans des tâches de prise-dépose des robots parallèles. Le but étant de remplacer les actionneurs électriques coûteux utilisés dans ce type de robots mais aussi de permettre le déplacement de charges plus élevées, le rapport poids/puissance de ce type d'actionneurs étant largement supérieur à celui des moteurs électriques. Le prix réduit des actionneurs électropneumatiques se justifie par leur principe de fonctionnement simple qui ne nécessite pas l'emploi de procédés de fabrication complexe. Cependant, l'énergie pneumatique et du fait de la compressibilité de l'air est une technologie fortement non linéaire difficile à modéliser et à commander. Cet article traite de la problématique de commande robuste des actionneurs électropneumatiques. L'état de l'art met en évidence l'utilisation d'une multitude de techniques de commande de ces actionneurs. De nombreuses commandes non linéaires ont été synthétisées comme les modes glissants [START_REF] Girin | Position-pressure robust control of an electropneumatic actuator[END_REF] [2] [START_REF] Van Damme | Sliding mode control of a 2-dof planar pneumatic manipulator[END_REF] ou encore la com-mande adaptative [START_REF] Zhu | Adaptive robust posture control of a parallel manipulator driven by pneumatic muscles[END_REF], [START_REF] Xiaocong | Integrated direct/indirect adaptive robust posture trajectory tracking control of a parallel manipulator driven by pneumatic muscles[END_REF]. D'autres auteurs ont privilégié la synthèse de contrôleurs linéaires après avoir effectué une linéarisation entrée/sortie [START_REF] Isidori | Nonlinear Control Systems. Third edition. Communications and Control Engineering Series[END_REF] [START_REF] Michel | Flatness and defect of non-linear systems : introduc-Comparaison de résultats de commande Commande en cascade Modes glissants Signal d'erreur[END_REF]. C'est le cas notamment de [START_REF] Osuka | H infinity control of a certain nonlinear actuator[END_REF] et [9] pour la commande H ∞ . Cependant, la linéarisation exacte avec comme sortie de commande la position implique l'existence d'une dynamique des zéros pour laquelle aucune preuve de stabilité globale n'existe encore [START_REF] Brun | Commandes linéaires et non linéaires en électropneumatique[END_REF]. C'est l'une des raisons qui conduisent à privilégier une stratégie de commande en cascade avec une linéarisation exacte dans la boucle de pression. A notre connaissance, seulement deux études de commande en cascade des actionneurs pneumatiques ont été menées jusqu'à présent [START_REF] Xiang | Block-oriented approximate feedback linearization for control of pneumatic actuator system[END_REF] [START_REF] Lee | A study on tracking position control of pneumatic actuators[END_REF]. Les auteurs utilisent uniquement de simples régulateurs PID dans les deux boucles de commande. Dans notre cas, deux contrôleurs avancés ont été synthétisés. Il s'agit d'abord de la commande prédictive qui est particulièrement adaptée aux tâches de robotique industrielle (les références futures étant connues a priori ). Il s'agit aussi de prouver l'applicabilité de la commande prédictive pour les actionneurs pneumatiques en se basant sur un modèle non linéaire (sans pour autant nécessiter une optimisation en ligne pour calculer le gain du contrôleur prédictif). Dans les deux seules applications de la commande prédictive sur les vérins électropneumatiques trouvées dans la littérature, les auteurs utilisent soit un modèle linéarisé autour d'un point d'équilibre [START_REF] Song | Improved control of a pneumatic actuator pulsed with pwm[END_REF]. ou alors un modèle déduit par apprentissage en utilisant les réseaux de neurones artificiels [START_REF] Norgaard | Intelligent predictive control of nonlinear processes using neural networks[END_REF]. Pour la boucle interne de pression, les inégalités matricielles (ou LMI pour Linear Matrix Inequalities) sont utilisées car elles permettent de combiner de façon intuitive des contraintes de placement de placement et des performances robustes H ∞ . Cet article est organisé comme suit. La section 2 introduit brièvement le banc d'essai expérimental ainsi que les modèles utilisés. En section 3, les techniques de commande sont présentées ainsi que la stratégie de commande en cascade qui sera appliquée. Les résultats expérimentaux sont finalement exposés en section 4. La pression atmosphérique est notée p 0 . Toute variation de volume de la chambre ou de la pression est décrite par la loi polytropique suivante [START_REF] Hildebrandt | A cascaded tracking control concept for pneumatic muscle actuators[END_REF] :

p 1 V γ 1 = p 2 V γ 2 (1)
dp dt = γ V (s) [rT q m (u, p) -p dV ds ṡ] (3)
où u représente l'entrée de commande de la valve, s est la position du piston dans le cas des cylindres. q m (u, p) représente le débit massique ( dm dt = q m (u, p)) qui est fonction de la pression et de la tension d'entrée. La relation

V i (s) = V i (0) ± A i s où V i (0) = A i l 2 est
le volume initial de la chambre i (La position initiale est le milieu du vérin). l représente la course du vérin. A i est la section du piston des deux côtés de la chambre. Dans le cas du vérin sans tige, les sections sont égales (A 1 = A 2 ). Le modèle de la valve a été approximé -après identification au sens des moindres carrés-par l'expression suivante (plus de détails dans [START_REF] Belgharbi | Analytical model of the flow stage of a pneumatic servodistributor for simulation and nonlinear control[END_REF]) :

q m (u, p) = ϕ(p) + ψ(p)u (4)
ϕ et ψ sont des fonctions polynômiales d'ordre 5 en p.

Pour la partie mécanique, le principe fondamental de la dynamique est appliqué en négligeant les frottements 5 :

s = 1 M [A 1 p 1 -A 2 p 2 ] ( 5
)
où M représente la masse équivalente des parties en mouvement ramenées au cylindre. Lors de la synthèse de la commande, M est considérée constante. La dépendance force pression -relation nécessaire lors de la synthèse de la commande-est linéaire et est donnée par :

F standard = A 1 p 1 -A 2 p 2 + F g (6)
F g est le terme prenant en compte les effets de la gravité.

Pour le vérin sans tige, nous avons la même expression avec

A 1 = A 2 = A.

III. Stratégie de commande pour les vérins électropneumatiques

La linéarisation exacte est présentée dans un premier temps. Celle-ci est basée sur des concepts de géométrie différentielle [START_REF] Isidori | Nonlinear Control Systems. Third edition. Communications and Control Engineering Series[END_REF] ou encore de platitude [START_REF] Michel | Flatness and defect of non-linear systems : introduc-Comparaison de résultats de commande Commande en cascade Modes glissants Signal d'erreur[END_REF]. Dans un second temps et une fois le système linéarisé obtenu, les techniques de commande utilisées (GPC et H ∞ avec contraintes LMI) ainsi que le schéma en cascade sont développés et expliqués dans cette section.

A. Linéarisation exacte

Il est très facile de prouver -dans le cas de la tension de la valve comme entrée de commande et la différence de pression entre les deux chambres du vérin comme sortie à commander-qu'aussi bien les critères de platitude ainsi que ceux de géométrie différentielle (linéarisation complète du système) sont satisfaits. Cependant, quand la position est choisie comme variable de sortie à commander, le degré relatif du système est égal à trois. Il est donc inférieur à la dimension du système d'état (qui est égale à 4). Les équations d'états sont les deux équations de variation de pression [START_REF] Michel | Flatness and defect of non-linear systems : introduc-Comparaison de résultats de commande Commande en cascade Modes glissants Signal d'erreur[END_REF], l'équation ds dt = ṡ et l'équation [START_REF] Xiaocong | Integrated direct/indirect adaptive robust posture trajectory tracking control of a parallel manipulator driven by pneumatic muscles[END_REF]. Une dynamique des zéros existe donc et le système n'est pas complètement linéarisable. Comme il n'existe pas de preuve de stabilité globale de cette dynamique résiduelle (voir [START_REF] Brun | Commandes linéaires et non linéaires en électropneumatique[END_REF]), on a donc tout intérêt à prendre la différence de pression comme sortie pour la linéarisation.

dp 1 dt = γrT V 1 (s) φ(p 1 ) - γ V 1 (s) p 1 dV 1 ds ṡ + γrT V 1 (s) ψ(p 1 )u dp 2 dt = γrT V 2 (s) φ(p 2 ) - γ V 2 (s) p 2 dV 2 ds ṡ - γrT V 2 (s) ψ(p 2 )u (7)
5. Une compensation de frottement qui ne se base pas sur un modèle de frottement est privilégiée Dans ces équations, on a considéré que u = u 1 = -u 2 ce qui veut dire que les chambres utilisent la même source de pression et que la commande utilisée dans une chambre est l'opposée de celle utilisé dans la seconde. Cette hypothèse implique de fait deux avantages : le système devient monoentrée ce qui simplifie la synthèse du contrôleur. L'utilisation d'une seule servovalve devient alors suffisante (comme cela est représenté en Fig. 2).

Et en remplaçant la variation de volume par son expression (ie, dV1 ds = A = -dV2 ds ), les équations suivantes sont obtenues :

dp 1 dt = γrT V 1 (s) φ(p 1 ) - γA V 1 (s) p 1 ṡ + γrT V 1 (s) ψ(p 1 )u dp 2 dt = γrT V 2 (s) φ(p 2 ) + γA V 2 (s) p 2 ṡ - γrT V 2 (s) ψ(p 2 )u (8)
La différence de pression est donnée par l'expression suivante :

d dt (p 1 -p 2 ) = [ γrT V 1 φ(p 1 ) - γrT V 2 φ(p 2 ) - γA V 1 p 1 ṡ - γA V 2 p 2 ṡ] f1(s, ṡ,p1,p2) + [ γrT V 1 ψ(p 1 ) + γrT V 2 ψ(p 2 )] f2(s, ṡ,p1,p2) u (9)
Ainsi donc, la linéarisation exacte est obtenue en appliquant la commande suivante :

u = [ γrT V 1 ψ(p 1 ) + γrT V 2 ψ(p 2 )] -1 • [u aux -( γrT V 1 φ(p 1 ) - γrT V 2 φ(p 2 ) - γA V 1 p 1 ṡ - γA V 2 p 2 ṡ)] (10)
Où, u aux est la nouvelle entrée de commande. Ceci conduit au système linéarisé suivant :

ṗ1 -ṗ2 = u aux (11)
Le contrôleur commandant la différence de pression en se basant sur ce système linéarisé est synthétisé. Pour cela des performances H ∞ combinées à des placements de pôles sous forme de contraintes LMI ont été utilisés. Cette méthode est décrite dans la prochaine sous-section.

B. Synthèse de commande H ∞ avec contraintes de placement de pôles avec utilisation de LMI

Une LMI est une contrainte de la forme :

A(x) = A 0 + x 1 A 1 + . . . + x N A N < 0 ( 12
D = {z ∈ C : f D (z) < 0} (13) avec : f D (z) = α + zβ + zβ T = [α kl + β kl z + β lk z] 1≤k,l≤m.
Par exemple, les régions LMI de type disque centré en (-q, 0) et de rayon r qui seront utilisées dans cette étude, sont définies ci-dessous [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF] :

f D (z) = -r q + z q + z -r < 0 (14)
Le problème de synthèse H ∞ sous contraintes peut être formalisé comme suit [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF]. En prenant comme système

ẋ(t) = Ax(t) + B 1 ω(t) + B 2 u(t) e(t) = C 1 x(t) + D 11 ω(t) + D 12 u(t) y(t) = C 2 x(t) + D 21 ω(t) + D 22 u(t) (15
ẋK (t) = A K x K (t) + B K y(t) u(t) = C K x K (t) + D K y(t) (16)
Ainsi, T we (s

) = D cl + C cl (sI -A cl ) -1 B cl avec 7 : A cl = A + B 2 D K C 2 B 2 C K B K C 2 A K , B cl = B 1 + B 2 D K D 21 B K D 21 , C cl = (C 1 + D 12 D K C 2 , D 12 C K ) and D cl = D 11 + D 12 D K D 21 .
Dans [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF], il est montré que la contrainte de placement de pôles est satisfaite si et seulement si il existe X D > 0 tel que :

[α kl X D + β kl A cl X D + β lk X D A T cl ] 1≤k,l≤m < 0 (17)
Quant à elle, la contrainte assurant les performances H ∞ est donnée sous forme de la LMI suivante :

  A cl X ∞ + X ∞ A T cl B cl X ∞ C T cl B T cl -γI D T cl C cl X ∞ D cl -γI   < 0 (18)
Ainsi, le problème de formulation du problème de synthèse H ∞ combiné aux problèmes de placement de pôles -en supposant que la même matrice de Lyapunov X > 0 est requise-est :

Trouver X > 0 et un contrôleur K(s) > 0 ≡ Ω K qui satisfont (17) et (18) avec X = X D = X ∞ ( 19
)
La matrice du contrôleur est notée de la manière suivante :

Ω K = A K B K C K D K ( 20
)
La difficulté dans cette formulation est que la combinaison des relations ( 17) et ( 18) implique des termes non linéaires de la forme B Ω K CX conduisant à une formulation du problème qui n'est plus convexe. Elle ne peut donc pas être formulée sous forme de LMI. Chilali et Gahinet [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF] ont résolu ce problème en proposant le changement de variable suivant :

B K = N B K + SB ∈ D K C K = C K M T + D K C ∈ R A K = N A K M T + N B K C ∈ R + SB∈C K M T + S(A + B 2 D K C 2 )R (21)
où R, S, N et M correspondent à la partition de X et de son inverse sous la forme

X = R M M T U , X -1 = S N N T V R ∈ R n×n , S ∈ R n×n (22)
La résolution du problème est donnée sous la forme du théorème suivant [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF] : Théorème : Soit D une région LMI arbitraire qui se trouve dans le demi-plan gauche du plan complexe et soit ( 13) sa fonction caractéristique. Dans ce cas, le problème modifié [START_REF] Maciejowski | Predicitve Control with Constraints[END_REF] admet -après changement de variables-une solution si l'ensemble des LMIs suivant est faisable.

Trouver R = R T ∈ R n×n , S = S T ∈ R n×n , et les matrices A K , B K ,C K et D K telles que R I I S > 0 (23) α kl R I I S + β kl Φ + β lk φ T k,l < 0 (24) Ψ 11 Ψ T 21 Ψ 21 Ψ 22 < 0 (25) avec la notation Φ = AR + B 2 C K A + B 2 D K C 2 A K SA + B K C ∈ (26)
les termes Ψ 11 , Ψ 12 et Ψ 22 sont détaillés dans [START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF]. Pour toute solution à cette LMI : -Calculer via la décomposition en valeurs singulières (SVD) la factorisation de rang plein M N T = I -RS de la matrice I -RS (M et N sont des matrices carrées inversibles) -Résoudre le système d'équations linéaires [START_REF] Ning | High steady-state accuracy pneumatic servo positioning system withPVA/PV control and friction compensation[END_REF] pour

B K , C K etA K (dans cet ordre). -calculer K(s) := D K + C K (sI -A K ) -1 B K .
Ainsi, K(s) est un contrôleur d'ordre n qui place les pôles en boucle fermée dans D tout en assurant T we ∞ < γ.

C. Commande prédictive généralisée (GPC)

La motivation de la commande prédictive est assez triviale, le problème de connaissance des trajectoires futures ne se posant pas pour des tâches de robotique industrielle. A cela, s'ajoutent tous les avantages inhérents à la commande prédictive, à savoir : l'applicabilité à une large classe de procédés dont les systèmes instables (c'est le cas de la fonction de transfert à commander), réduction du temps de réponse de façon intuitive et existence d'une solution explicite dans le cas monovariable linéaire. Toute commande prédictive s'appuie sur trois principes : l'utilisation d'un modèle de prédiction sous forme explicite, l'utilisation d'un critère de prediction à optimiser, et enfin le recalage de l'horizon de prédiction à chaque instant t, la commande étant réactualisée à chaque temps d'échantillonnage. Dans le cas de la GPC, l'algorithme utilise le modèle de prédiction CARIMA donné par :

A(z -1 )y(t) = z -d B(z -1 )u(t -1) + C(z -1 )e(t) (27)
où u(t) et y(t) sont respectivement les séquences de commande et de sortie du procédé et e(t) est un bruit blanc.

A, B et C sont des polynômes en z -1 . Ils sont donnés par :

A(z -1 ) = 1 + a 1 z -1 + a 2 z -2 + ... + a na z -na , B(z -1 ) = b 0 + b 1 z -1 + b 2 z -2 + ... + b nb z -nb et C(z -1 ) = 1 + c 1 z -1 + c 2 z -2 + ... + c nc z -nc et = 1 -z -1
. Afin de simplifier les calculs, il est supposé que C -1 est égal à 1. L'algorithme de la GPC consiste à appliquer une séquence de commande qui minimise la fonction coût multidimensionnelle suivante :

J(N 1 , N 2 , N u ) = N2 j=N1 δ(j)[ŷ(t + j|t) -w(t + j)] 2 + N2 j=1 λ(j)[ u(t + j -1)] 2
(28) où ŷ(t + j|t) est l'optimum sur un horizon j calculé à l'instant t, N 1 et N 2 sont les horizons de prédiction minimum et maximum, N u représente l'horizon de commande, δ(j) et λ(j) sont les séquences de pondération et w(t + j) est la référence future. La minimisation de la fonction coût conduit à l'obtention de la future séquence de commande u(t), u(t + 1), ... où la sortie y(t + j) est proche de w(t + j).

L'expression de la sortie prédite peut être mise sous la forme compacte suivante :

y = Gu + f (29)
Le terme G et f sont définis dans [START_REF] Bordons | Model Predicitve Control. Second Edition[END_REF] ou [START_REF] Maciejowski | Predicitve Control with Constraints[END_REF]. Le critère peut être ainsi réécrit :

J = (Gu + f -w) T (Gu + f -w) + λu T u (30)
où w = [w(t + d + 1), w(t + d + 1) . . . w(t + d + N )] T . L'équation (30) prend la forme :

J = 1 2 u T Hu + b T u + f 0 (31) avec H = 2(G T G + λI), b T = 2(f -w) T G et f 0 = (f -w) T (f -w)
Ainsi, le minimum de J est déterminé en annulant le gradient de J, ce qui donne :

u = -H -1 b = (G T G + λI) -1 G T (w -f ) (32)
Puisque le signal de commande qui est appliqué au procédé est le premier élément du vecteur u (recalage), l'expression de la commande est donnée par :

u(t) = K(w -f ) (33)
où K représente le premier élément de la matrice (G T G + λI) -1 G T . Contrairement aux contrôleurs conventionnels, les régulateurs prédictifs dépendent uniquement des entrées futures.

D. Stratégie de commande en cascade position/pression

Le schéma global de la stratégie de commande en cascade est représenté en Fig. 4. La boucle interne est celle qui commande la différence de pression entre les chambres. Pour cela, le régulateur H ∞ sous contraintes présenté précédemment a été implementé. Pour la boucle externe, le régulateur GPC a été synthétisé en se basant sur la fonction transfert en boucle fermé de la boucle de pression en série avec la mécanique du bras (double intégration pour obtenir en sortie la position). Comme présenté précédemment, le système linéarisé de la boucle de pression est un intégrateur. Ainsi, les équations du système prennent la forme suivante :

ẋ = (0)x + 1 0 b v + (1)u y u = -1 0 x + 1 0 0 0 b v + 0 1 u y = (1)x + 0 1 b v + (0)u (34)
où ω = (b v) T et e = (y u) T ont été définis préalablement. Les hypothèses classiques pour la synthèse H ∞ sont facilement vérifiables 8 . Le résultat de l'optimisation LMI [START_REF] Gahinet | The lmi control toolbox[END_REF] donne le contrôleur de la forme suivante :

K(s) = K 1 s + K 2 s + K 3 (35) 8. (A, B 2 ) commandable et (C 2 , A) observable
Ainsi, la fonction de transfert équivalente en boucle fermée est de la forme

G clpressure (s) = K 1 s + K 2 s 2 + (K 1 + K 3 )s + K 2 (36)
En se basant sur cette fonction de transfert, dont la sortie est la différence de pression, il est facile de reconstruire la position en divisant par la masse et en intégrant deux fois la sortie [START_REF] Xiaocong | Integrated direct/indirect adaptive robust posture trajectory tracking control of a parallel manipulator driven by pneumatic muscles[END_REF]. Ainsi, le régulateur GPC est synthétisé en se basant sur la fonction de transfert suivante :

G position (s) = K 1 s + K 2 M s 2 (s 2 + (K 1 + K 3 )s + K 2 ) = G clpressure M s 2 (37)
IV. Résultats expérimentaux Différents tests ont été menés de façon à prouver expérimentalement l'intérêt d'une telle approche de commande en cascade. Durant les différents tests, le terme F g dans ( 6) a été négligé. Il a été observé qu'il n'a pas un grand impact sur les performances. Ainsi, le processus de synthèse est simplifié et ceci nous permet de dire que la commande synthétisée est robuste vis-à-vis des effets de la gravité et des variations de l'orientation du bras. En

FigFig. 2 .

 2 Fig. 1.Photo du banc expérimental.

  ) une région LMI D, et des performances H ∞ données par γ > 0, il s'agit de trouver la loi de commande u = K(s)y telle que : (1) les pôles en boucle fermée sont dans D et (2) T we ∞ < γ où T we (s) est la fonction de transfert en boucle fermée ayant pour entrée ω et comme sortie e. Ceci est schématisé en Fig. 3. avec le vecteur d'entrée

Fig. 3 .

 3 Fig. 3. Schéma de commande et notations

7 .

 7 D 22 = 0, cette hypothèse simplifie beaucoup les différents calculs. (Il est à noter qu'il est toujours possible d'avoir D 22 = 0 en faisant de simples changements de variables)

Fig. 5 ,

 5 un test de commande avec le bras sans masse est réalisé. La fréquence d'échantillonnage est de 0, 2khz. Pour le contrôleur H ∞ , une région LMI de type disque a été Suivi de la différence de pression [bar] Commande en position [degrés]

Fig. 5 .

 5 Fig. 5. Réponse à un signal carré (vérin sans tige)

  ) où A 0 . . . A N sont des matrices symétriques données et x T = (x 1 . . . x N ) est le vecteur des variables inconnues. Dans[START_REF] Lee | A study on tracking position control of pneumatic actuators[END_REF], le symbole < veut dire définie négative6 . Une sous région D du plan complexe est appelée région LMI s'il existe une matrice symétrique α = [α kl ] ∈ R m×m et une matrice β = [β kl ] ∈ R m×m telles que :

Une façon intuitive de synthèse d'un contrôleur H ∞ tout en agissant sur le régime transitoire est de combiner performances H ∞ et placement de pôles. Les pôles seront contraints à être dans des régions appelées régions LMI et introduites pour la première fois dans

[START_REF] Chilali | H-infinity design with pole placement constraints : an lmi approach[END_REF]

.

6. la plus grande valeur singulière est négative Définition : Régions LMI.

Le contrôleur prédictif a été synthétisé en se basant sur la fonction de transfert donnée dans (37). Le système a été discrétisé avec un temps d'échantillonnage de T s = 5ms. Les polynômes du modèle CARIMA ainsi obtenus sont :

La pondération λ a été choisie égale à 5 • 10 -6 et δ a été fixé à 1. L'horizon de prédiction est égal à 50. Sur la Fig. 5, un signal carré d'amplitude égale à -+ 15˚a été testé. Le temps de montée est très faible et est égal à 0, 045s. Le temps de réponse est quant à lui égal à 0, 24s. L'erreur statique est inférieure à 0, 01˚. Le saut observé à t = 3s est dû au phénomène du redécollage des actionneurs électropneumatiques (voir chapitre 3 de [START_REF] Brun | Commandes linéaires et non linéaires en électropneumatique[END_REF]). Ce problème a été résolu par deux approches différentes. La première consiste à agir sur les gains de la boucle de pression via l'emplacement des pôles en boucle fermée. Ceci est illustré en Fig. 6 où les pôles ont été déplacés vers un cercle centré en (-450, 0). L'apport de cette approche est montré en Fig. 7 où l'erreur statique est inférieure à 0, 01˚. Une autre approche a consisté à intégrer au schéma de commande de la Fig. 4 l'algorithme de compensation de frottement proposé par [START_REF] Ning | High steady-state accuracy pneumatic servo positioning system withPVA/PV control and friction compensation[END_REF]. Les résultats obtenus sont similaires à ceux de la méthode précédente. Un test de robustesse est effectué en rajoutant une masse de 5kg qui n'est pas prise en compte dans le modèle de commande. Les résultats sont montrés en Fig. 8. Les résultats sont satisfaisants et il est intéressant de voir la réaction du système par un saut dans la différence de pression (et donc dans la force) qui correspond à la charge additionnelle qui a été rajoutée. Le contrôleur en pression réagit automatiquement à ce changement de consigne et de ce point de vue, la structure en cascade assure déjà une robustesse naturelle qui est renforcée par le choix des techniques de commande. Un test de comparaison avec les modes glissants [START_REF] Smaoui | Commandes non linéaires robustes mono et multidimensionnelles de dispositifs électropneumatiques : synthèse et applications[END_REF] a été V. Conclusion Une stratégie de commande en cascade intégrant une boucle externe de position et une boucle interne de différence de pression est proposée pour les vérins électropneumatiques. Pour la position, un contrôleur de type GPC est synthétisé. Le choix de la commande prédictive permet de réduire de façon naturelle le temps de réponse du système qui anticipe tout changement dans la consigne. Pour la boucle interne, un régulateur H ∞ avec contraintes de placement de pôles est synthétisé par optimisation LMI. La synthèse des régulateurs est faite en se basant sur le modèle linéarisé du vérin et ce après linéarisation exacte de ce modèle.

Les lois de commande ont été implémentées expérimentalement sur banc de commande et les résultats montrent bien l'apport d'une telle approche en termes de performances et de robustesse.


Thierry Ranchin
email: ranchin@cenerg.cma.fr

Lucien Wald
email: wald@cenerg.cma.fr

Benefits of fusion of high spatial and spectral resolutions images for urban mapping

This communication presents a method allowing the combination of high spatial and spectral resolution images to improve the mapping of urban area. It makes use of two mathematical tools, the wavelet transform and the multiresolution analysis which are shortly introduced. This method, called ARSIS is then described. A qualitative comparison between a classical merging technique and the ARSIS method is presented through an example, the merging of SPOT XS images and the russian sensor KVR-1000. Then, the process used for urban mapping is explained and a map of a district in Riyadh, enhancing the roads and the buildings is presented.

I. INTRODUCTION

Urban mapping seems to be a promising market for the earth observation satellites in the next decade. For this application, on the one hand, the high spatial resolution images are necessary for an accurate geometrical description of the cities; on the other hand, the high spectral resolution images discriminate the different types of urban structures. Unfortunately, no civil satellite sensor is available to provide high spatial and spectral resolutions images at the same time. Airborne solutions are expensive and are often restricted by flight authorizations. However, mathematical solutions exist and sensor fusion techniques provide high spectral and spatial resolutions images. Furthermore, sensor fusion is the best solution to minimize the costs.

In section II, two mathematical tools are shorthly introduced and the ARSIS method (from its french name "amélioration de la résolution spatiale par injection de structures") is presented. An example of application of this method to the merging of SPOT XS images (high spectral resolution; spatial resolution: 20 m) and to russian images KVR-1000 (panchromatic band; spatial resolution: 2 m) and a qualitative comparison between the resulting images and the images obtained by a Intensity Hue Saturation (IHS) method are presented in section III. The process used for urban mapping is described in section IV and an extract of a urban map showing the roads and buildings in a district of the town of Ryadh (Saudi Arabia) presented.

II. ARSIS METHOD

For the application of a sensor fusion technique, images need to represent the same area, to be superposable, and it is required that no major change has occured in the observed landscape. The set of images to merge is supposed to have different spatial and spectral resolutions. Many methods have been proposed to enhance the spatial resolution of images owing the presence of one or more images of the same scene of better spatial resolution (see for example [START_REF] Carper | The use of Intensity-Hue-Saturation transformations for merging SPOT panchromatic and multispectral image data[END_REF][START_REF] Chavez | Comparison of three different methods to merge multiresolution and multispectral data: Landsat TM and SPOT panchromatic[END_REF]. But, not a lot of them take care of the multispectral content when increasing the spatial resolution. The ARSIS method was designed to synthetize in a set of images with different spatial and spectral resolutions, multispectral images with the best resolution available in the set of images and to preserve the radiometric content of original images. This method makes use of the wavelet transform and the multiresolution analysis. The multiresolution analysis was introduced by [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. It allows the computation of successive approximations of the same image with coarser and coarser spatial resolutions. Combined with the multiresolution analysis, the wavelet transform allows the description of the difference of information between two successive approximations of the same image. Figure 1 present the general scheme of the ARSIS method. The best resolution image A is decomposed by a multiresolution analysis using a wavelet transform. The differences of information provided by the wavelet transform are modelized by the wavelet coefficients. The second image B with a coarser spatial resolution than image A is also decomposed through the multiresolution analysis with the wavelet transform. A model describing the transformation of the wavelet coefficients of the image A to the wavelet coefficients of the image B is established. This model takes into account the physics of both images and the correlation or anti-correlation existing between both coefficients images. This model is used to computed the wavelet coefficients needed to synthetize the image B at the spatial resolution A.

Resolution n°1

Resolution n°2 Resolution n°3 This synthesis is made by a reconstruction, WT -1 in the general scheme which is the inverse operation of the multiresolution analysis. ARSIS allows the preservation of the spectral content of each image and an improvement of the spatial resolution of each images up to the best one available in the original set. A more complete description of the ARSIS method can be found in [START_REF] Ranchin | Applications de la transformée en ondelettes et de l'analyse multirésolution au traitement des images de télédétection[END_REF]. In the case of SPOT imagery, ARSIS was shown to give the best achievable results in terms of preservation of the original content [START_REF] Mangolini | Evaluation de la qualité des images multispectrales à haute résolution spatiale dérivées de SPOT[END_REF].

Spectral

III. EXAMPLE AND COMPARISON

An application of the ARSIS method is presented in this section. The set of images is composed by a SPOT XS scene of the town of Riyadh (Saudi Arabia) acquired the 16th of May 1993 and a russian image KVR-1000 of the same area acquired the 7th of September 1992. The three XS images have a spatial resolution of 20 m and a spectral range of 0.5 -0.59 µm for XS1 band, of 0.61 -0.68 µm for XS2 band and of 0.79 -0.89 µm for XS3 band. The KVR-1000 image has a spatial resolution of 2 m and a spectral range of 0.51 -0.71 µm equivalent to a panchromatic band. Two merging processes are applied to this set of images. The ARSIS method described in the previous section, and the Intensity-Hue-Saturation (IHS) method [START_REF] Carper | The use of Intensity-Hue-Saturation transformations for merging SPOT panchromatic and multispectral image data[END_REF]). To allow a comparison of the two methods and the benefits provided by the sensor fusion, three extracts of the same scene are presented. Figure 2 shows a composition of the XS images, Figure 3 the result of the IHS method and Figure 4 the result of the ARSIS method. The big crossroads at the upper left corner of image is of interest. One can see that the resulting image provided by the ARSIS method is close to the XS original composition than the resulting image of the IHS method. Due to the gap of resolution between the XS and the KVR images, it is difficult to propose a method to quantitatively estimated the resulting images. Even if the geometrical quality of both images seems to be very close, the images resulting from the ARSIS method allows to see all the roads on the crossroads even the lower left loop which exists and is also represented in the original images. Hence the images resulting provided by ARSIS where prefered for the urban mapping. This classification allows to discriminate two classes: the roads and the buildings. This result is used to help for the manual photo-interpretation of the area, and the urban map is derived.

V. CONCLUSION

In this communication, we have shown the benefits of sensor fusion for urban mapping. The possibility of having multispectral images at high spatial resolution with a preservation of their spectral content by the use of the ARSIS method, enables the application of classification algorithms. This allows to save time of manual photointerpretation and to improve the quality of the results. The ground-truth performed in this area as shown that all the roads where well classified and that close to all buildings also. The differences observed are principally due to the differnce of date between the ground-truth and the images.

Fig. 1 .

 1 Fig. 1. General scheme of the ARSIS method

Fig. 2 .

 2 Fig. 2. XS original composition


Thierry Ranchin
email: ranchin@cenerg.cma.fr

Lucien Wald
email: wald@cenerg.cma.fr

Merging SPOT-P and KVR-1000 images for updating urban maps

This communication presents a method allowing the combination of images acquired on the same area with different dates for updating urban maps. This method makes use of two mathematical tools, the wavelet transform and the multiresolution analysis which are shortly introduced. An example is provided in Jeddah (Saudi Aradia) for the updating of a KVR-1000 image by a SPOT panchromatic image acquired at a more recent date. Due to the difficulty to obtain KVR-1000 images at recent dates, this method allows the combination of the high spatial resolution of this kind of images and the regular acquisition of images by the different Earth observation satellites.

I. INTRODUCTION

High spatial resolution images are necessary to urban mapping, due to the great accuracy needed for this type of application. Till the end of 1992, it is possible to acquire images from the russian sensor KVR-1000. These images are provided in the panchromatic band with a spatial resolution of 2 m. Before the end of the century, new civil missions with high spatial resolutions close to 1 or 2 m are planned and will answer to the need of accuracy for urban mapping. Unfortunately, these missions have a limited life time (KVR-1000) or may be limited in access because of cloud cover, repetitivity, site access, data access, operational constraints, ... These limitations may prevent the update of information content of the area of interest. This communication proposes to benefit from the operational SPOT missions to provide an update of urban area by the means of sensor fusion techniques. The merging process is based on the multiresolution analysis and the wavelet transform presented in section II. The merging process itself is presented in section III. An example of the merging of a KVR-1000 image with a more recent panchromatic image acquired by the SPOT satellite is explained in section IV. A comparison of the two dates of acquisition is provided, enhancing the necessity of the updating process. The benefits of using such a fusion approach are discussed.

II. WAVELET TRANSFORM AND MULTIRESOLUTION ANALYSIS

The merging process proposed is based on the wavelet transform and the multiresolution analysis. The multiresolution analysis was introduced by [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. It allows the computation of successive approximations of the same image at coarser and coarser spatial resolutions. This can be represent by a pyramid, where the basis is the original image and the different approximations of the original image, the floors of the pyramid. The theoretical limit of this multiresolution analysis is a pixel that will represent the mean of the original image. In this representation, the difference of information existing between two successive approximations is represented by the wavelet coefficients computed from a wavelet transform. The wavelet basis is computed by dilatation and translation of an unique function called the mother wavelet. The image is decomposed on this basis, and a space-scale representation is obtained. The multiresolution analysis using the wavelet transform is an inversible operation and provided an exact reconstruction of the original image. Hence, from one approximation of the image and the different wavelet coefficients representing the differences of information between the original image and this approximation, it is possible to reconstruct exactly the original image. Ranchin and Wald (1993) have presented some possibilities of use of these tools for the field of remote sensing. For a more mathematical introduction, one can refer for example to [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF], [START_REF] Rioul | Wavelets and signal processing[END_REF], [START_REF] Daubechies | Ten lectures on wavelets[END_REF] and [START_REF] Meyer | Wavelets. Algorithms & Applications[END_REF].

III. MERGING METHOD

The data available for the urban mapping are:

-a high spatial resolution image A acquired at the first date, -a most recent image B acquired with a spatial resolution lower than the first image one.

The two images are supposed to represent the same area and to be superposable. If we considered that the image B contains the most update information, we need to give more importance to this image. But for urban mapping, the high spatial resolution is of great importance due to the scales needed to well represented the different kind of structures in the area. So it is needed to preserve the fine details provided by the image A and to introduce them in the most recent representation, the image B. As explained in the previous section, the wavelet transform and the multiresolution analysis are two mathematicals tools very well-adapted to represent an image at different spatial resolutions and to modelize the difference of information between these resolutions. The merging process will be the following. A multiresolution analysis using the wavelet transform is first applied to the image A, providing an approximation of the original image at a spatial resolution close to the image B resolution. The wavelet coefficients represent the difference of information between the spatial resolution of the image A and the one of the approximation computed. The differences represented by the wavelet coefficients and the image B, are used to reconstruct an image with the spatial resolution of the image A, the finest structures available in this image and the information content of the image B. In this approach, the geometrical information is privileged. The aim is not to preserve the spectral content of the original image as proposed in the ARSIS method (Ranchin et al. 1993), but to represent the most complete geometrical information.

IV. EXAMPLE

To illustrate this sensor fusion technique, an example is provided on a district of the town of Jeddah (Saudi Arabia). Figure 1 presents the KVR-1000 image acquired in March 16, 1991. This image has a spatial resolution of 2 m and a spectral range in the panchromatic band and represents an area of 1024 m by 1024 m. This district of the town is in a very modern part with perpendicular roads and square blocks of buildings. One can notice the very fine details in this image giving the possibility to photo-interpret this area to construct urban map. Figure 2 presents the SPOT panchromatic image acquired in September 5, 1993. To allow a comparison of the two images, this image was processed using a nearest-neighbour interpolation to obtain the same spatial resolution. The effective spatial resolution of the image is 10 m but the pixels have a resolution of 2 m. The big lower right square is of interest. In the middle of this square, major changes can be noticed between the image of 1991 and the image of 1993. It seems that new buildings were constructed even if it is sometime difficult to evaluate how many. On the upper right part, the same phenomenon is visible. One building has been constructed in a block where nothing exists in the KVR-1000 image. Figure 1 illustrates the necessity to have a good geometrical description of the area, allowing to distinguish the different buildings. Even if the SPOT panchromatic image has a lower spatial resolution, the regular acquisition of images allows to update the information on a given area. Hence it is necessary to perform an updating of the information, to work on the most complete image. To obtain the merged image, a multiresolution analysis using the wavelet transform was applied to the KVR-1000 image. The wavelet coefficients represent the fine structures missing to the SPOT panchromatic image. Hence, these differences are introduced in the panchromatic pyramid and the merged image is constructed. The multiresolution analysis algorithm used in this merging process is called the "à trous" algorithm and was proposed by [START_REF] Dutilleux | An implementation of the "algorithme à trous" to compute the wavelet transform[END_REF]. Figure 3 presents the result of the merging process. The different structures missing in the original KVR-1000 image are now available on this image. The different buildings constructed between March 1991 and September 1993 are visible in this image and the photo-interpretation can be realized only from this image which contains all the structures existing in both images. This image allows to build the urban map, enhancing the areas where a survey is needed or describing the differences between the two dates. Figure 4 shows the urban map realized on the lower part of the image. The photointerpreter has chosen to propose a survey for the area where changes have occured.

V. CONCLUSION

In this communication, a merging process was proposed to allow an updating of an image by a more recent one. It was shown that the merging process preserves the best of the two images, in one hand the high spatial resolution of the KVR-1000 image and in the other hand the most update information provided by the SPOT panchromatic image. This merging process can be used for a periodical updating of urban maps, allowing the reduction of the time necessary for the survey and the cost of this survey. It allows also the construction of a single representation of an area with all the informations available from different images.

Fig. 1 .

 1 Fig.1. original KVR-1000 image acquired in March 1991

Fig. 3 .

 3 Fig. 3. Image resulting of the merging process.

Fig. 4 .

 4 Fig. 4. Urban map derived from Fig. 3. The blocks where changes have occured are designed by areas to survey.