Adrien Herubel

Venceslas Biri

S Deverly

Gaspard Monge

Duran Duboi

Autonomous Lighting Agents in Photon Mapping

come   L'archive ouverte pluridisciplinaire

Introduction

Physically-based global illumination algorithms such as photon mapping [START_REF] Jensen | Realistic Image Synthesis Using Photon Mapping[END_REF] have a linear progression between complexity and quality. To a given quality, rendering time scales linearly with computer performances. With Moore's law call in question and increasing demand in quality, those algorithms need more and more optimisations. Classical optimisations such as irradiance caching [START_REF] Ward | A ray tracing solution for diffuse interreflection[END_REF] or shadow photons [START_REF] Jensen | Efficiently Rendering Shadows using the Photon Map[END_REF] are themselves linear in performance gain. They usually consist of adding knowledge about the scene and using it to interpolate previously computed values. Data is gathered in large volumes although due to heavy redundancy we observe a low density of useful data (see Fig 1.B and C).

We propose a brand new optimisation method for the photon mapping global illumination algorithm which uses Artificial Intelligence structures and concepts to resolve the sparse data problem in global illumination presented above. More precisely, we create our own autonomous agent structure called Autonomous Lighting Agents (ALA) that performs an agent-based scene discovery. They can efficiently gather and store large amounts of useful data (see Fig 1 .D andE) thanks to various sensors. Moreover ALA exist within multiple graph representations that depict relations between agents such as neighbourhood or light paths. The resulting structures are networks of agents storing data non-uniformly with increased density in area of interest. Then any ALA and the associated networks can be queried during rendering to make decisions regarding, for example, irradiance caching or shadow casting. Therefore, fewer photons are cast reducing computing time during both the photon casting and the rendering phases of the algorithm.

In the next section, we present the photon mapping algorithm and overview the main optimisations on that method. Then we present our new criterion to evaluate usefulness of local data density. To match closely this particular data density, we detail our method showing how ALA discover the scene and how they are queried for final rendering. Finally, we present our first results showing an important decrease of memory occupation and slightly shorter rendering times, compared to optimized photon mapping.

Global illumination and photon mapping

We call global illumination the simulation of all light scattering phenomena in a virtual scene. Photon mapping presented in [START_REF] Jensen | Importance driven path tracing using the photon map[END_REF][START_REF] Jensen | Global Illumination using Photon Maps[END_REF] is currently one of the most efficient physically based algorithm capable of global illumination rendering. Compared to other global illumination methods like radiosity [START_REF] Goral | Modeling the interaction of light between diffuse surfaces[END_REF][START_REF] Keller | Instant radiosity[END_REF] or Metropolis Light Transport [START_REF] Veach | Metropolis light transport[END_REF][START_REF] Segovia | Metropolis Instant Radiosity[END_REF], photon mapping is a robust and consistent two-pass algorithm [START_REF] Hachisuka | Progressive photon mapping[END_REF], is able to handle many light effects including caustics, and is modular since it separates illumination in several layers [START_REF] Jensen | Importance driven path tracing using the photon map[END_REF]. The first pass traces photons from all lights through the scene and stores them, at hit points, in so called photon map. The second pass uses ray tracing to render the image. In the ray tracing pass the photon map is used to estimate the radiance at different locations within the scene. This is done by locating the nearest photons and performing a nearest neighbor density estimation (see Fig. 2). A detailed presentation of this method can be found in [START_REF] Jensen | Realistic Image Synthesis Using Photon Mapping[END_REF].

Despite its numerous qualities, photon mapping is still a slow and complex algorithm. Classical optimisations include shadow photons, direct computing and importance sampling of the indirect illumination evaluation. Shadow photons, Fig. 2. The two passes of photon mapping as shown in [START_REF] Jensen | Efficiently Rendering Shadows using the Photon Map[END_REF], are stored in a separate photon map which is used to speedup ray traced shadow computing. The photon map can also quickly compute indirect illumination since it has a low variance. To achieve this, diffusely reflected photons must be stored in a separate photon map that will be queried directly to compute indirect illumination at a given point. It is advised to use precomputed irradiance described in [START_REF] Christensen | Faster photon map global illumination[END_REF] to speed up the photons gathering. Finally, for the complete indirect illumination evaluation, the photon map can be used to accelerate computing. This is done thanks to importance sampling [START_REF] Jensen | Importance driven path tracing using the photon map[END_REF] by gathering the diffusely reflected photons stored near the hit point and retrieving their incoming direction. Ward et al [START_REF] Ward | A ray tracing solution for diffuse interreflection[END_REF] proposed also irradiance caching, used to accelerate the indirect diffuse illumination. The key idea is to interpolate, if possible, irradiance value between several points which are relatively close and share the same orientation. Unfortunately, this optimization is itself costly and not adapted to multithreaded implementations of path tracing.

Advanced optimization address the drawback of the density of the photon map. Indeed, the photon distribution is driven by lights, and for large models where only a small part is visible from the viewer, like the example scene presented in [START_REF] Veach | Metropolis light transport[END_REF], it may be better to cast photon only where they are needed. Techniques like [START_REF] Pietrek | Importance driven construction of photon maps[END_REF][START_REF] Fan | Metropolis photon sampling with optional user guidance[END_REF] focus on the camera to influence the direction of photons casted by lights. However this techniques, like Importons, may induce high variance in the photon map by generating high-power photons as explained in [START_REF] Christensen | Faster photon map global illumination[END_REF] (pp 146). Density control, introduced in [START_REF] Suykens | Density control for photon maps[END_REF], consists of limiting the presence of photons in bright area allowing more photons in darkest areas. Recent algorithms [START_REF] Havran | Fast final gathering via reverse photon mapping[END_REF][START_REF] Hachisuka | Progressive photon mapping[END_REF], inverting the two passes of photon mapping, allow to minimize the bias, but they are very dependent on the view position and still use a traditional photon mapping inducing the same drawbacks in photon density. Nevertheless, the locality of photon maps in [START_REF] Hachisuka | Progressive photon mapping[END_REF] involves a limited density for each iteration of the progressive photon mapping.

Local useful data density

In photon mapping, lights are responsible for the distribution of photons and an inadequate repartition will keep the photon map from playing its role in rendering. The first pass sends an important amount of photons with only slight variations in their own data. Unfortunately, they are required to achieve sufficient visual quality, or noise can appear as shown in [START_REF] Jensen | Realistic Image Synthesis Using Photon Mapping[END_REF] (pp 150).

Classical optimization, as density control or inversion techniques, presented in 2, address partially this issue by restricting the number of photons using local density and visual importance. But they are very view dependant, still send a lot of photons or assume strong variance in lighting in order to be efficient, especially for the density issue.

Therefore, our first goal is to found a mean of characterizing the usefulness of data in each area of the scene. Indeed, if we observe Fig. 1.B. we find a high density photons which are locally very similar. It leads us to introduce the concept of local useful data density (LUDD) in the photon map to represent the local density of useful data for any criterion. First, for a particular criterion k, we define the local density of useful k-data, denoted by δ k . It is defined as the standard deviation of k computed for N included samples {x kj } on a local sphere.

δ k = 1 N -1 N i=1 x ki -1 N N j=1 x kj 2 1 N N j=1 x kj (1)
The whole local density of useful data δ for K criteria is simply :

δ = K k=1 δ k N (2)
Criteria used to perform this measure include surface normal, light power or direction. This measure can be expanded for every additional data carried by the photon map such as visual importance. The δ value computed Fig. 1.C. is particularly low which implies that added value of each photon is very small. However in Fig. 1.E., δ is higher which means our agent map has locally a strong useful data density. Our observations show that local useful data density is not evenly distributed in the photon map and that this measure is strictly correlated to local conditions in the scene.

4 Our method using Autonomous Lighting Agents

An autonomous agent approach

To maximize δ measure, we need to create a new structure dedicated to scene discovery instead of carrying light flux with an evenly strong δ. We want it to be independent of human interaction, of any tessellation, and of visual importance. The ability to discover local phenomena as geometry gradient and light variations leads us to the autonomous agent software paradigm [START_REF] Maes | Designing autonomous agents: theory and practice from biology to engineering and back[END_REF]. In consequence our structure is not a monolithic entity but a swarm of agents able to perceive their environment through sensors [START_REF] Maes | Modeling adaptive autonomous agents[END_REF] and to make behavioural decisions based on the given criteria. We define our autonomous lighting agents (ALA) as entities with a temporal persistence coexisting in a scene with other ALA and interacting with them [START_REF] Treuil | Modélisation et simulation à base d'agents[END_REF].

We use formal description as described in litterature [START_REF] Lerman | A macroscopic analytical model of collaboration in distributed robotic systems[END_REF] to represent mathematically our model. We give each agent an immutable and unique identity i. An ALA has three abilities, it can be cast by an arbitrary light source, it can discover the area around its location and it can answer simple questions related to this area. The model also specifies attributes or sensors owned by agents and denoted by a. Therefore, our new ALA improved photon mapping is organized as follows :

1. Replace the classical photon casting pass with ALA casting pass (see 4.2) 2. Make ALA discover their environment (see 4.3) 3. Do the rendering pass using ALA decision making for different layers (see 4.4)

ALA casting

We choose to cast ALA from the lights, uniformly on a sphere. When hitting a surface an ALA acts like both a photon and a shadow photon in the photon mapping algorithm, as it is stored and duplicated. Like shadow photon, the obtained ALA clone passes through the surface. In consequence ALA are present in occluded zones. Agents are stored using a kd-tree in the same way than photon maps. A visualization of the ALA network is shown in Fig. 1.D. After this stage, agents are evenly distributed in the scene at location {x i }.

ALA deployment

Local discovery with sensors Each agent can deploy sensors to gather local data about the scene and other agents. The photon map only provides data at each discrete position of each photon whereas our ALA deploy sensors on a local parametric discovery area. A sensor is identified with an identity j and is attached to an agent i. Each sensor corresponds also to an attribute a, so we define a value discovered by a sensor as a ij . Each agent can deploy more than one sensor for each kind of attribute a. Sensor position are denoted by x a ij . Finally, for each sensor a, we define a coefficient of interest σ a , computed by the difference of the value perceived by the sensor with an arbitrary chosen median value σ m a . Sensor used are light sensor, penumbra sensor, geometry gradient sensor and proximity sensor. Table 1 shows the different measures and median values used for each sensor.

Gold miner decision algorithm

As we aim to achieve an optimal LUDD repartition across the scene, ALA are entitled with the decision of choosing the quantity of data they carry. In the classical photon map, LUDD is mainly related to local photon density. ALA strongly differ on this point, as data is not related by the presence of the agent itself. A single agent can deploy dozens of sensors on large distances as well as many agents can choose not to use any. Each agent will choose which types of sensors are deployed, how many of them and how far they can go. We use a decision algorithm making the agent acting like a gold miner (see Alg. Table 1. Different sensors used with their measurement and median value used

A sensor reports its coeficient of interest σ, an interesting result will increase the stimulation factor γ of the agent. On the contrary, a non-pertinent result will decrease this factor. Stimulation factors include : little number of sensors, strong irradiance gradient, high geometric gradient, high or low number of neighbours. Discouragement factors include great number of sensors or planarity of surface. Factors are of different natures and some of them implies computing the gradient of two given criteria a ij . Between each deployment we test if an excitement or unexciting factor is triggered. The agent will continue to deploy sensor and gather data until the stimulation factor drops below a certain level.

Computing the proximity network When all ALA have been cast and all sensors are deployed, we compute the proximity network on the whole agent set. The network is a graph in which each agent is a node related to its closest neighbours in the scene space. The graph is computed quickly using the kd-tree to find the neighbours for each agent.

Algorithm 1 Gold miner algorithm γ = 0.0, σ = 0.0; for all sensor types do while γ > 0 do direction = random() distance += ǫ σ = deploy sensor(direction,distance) γ += σ -( σ/10) end while end for

Rendering pass using ALA

The ALA structure allows an efficient agent based local discovery of the scene, data is then used during the rendering pass. The ALA network aim to replace the photon map for all kind of decision making optimisations therefore avoiding the costly generation of direct photons, indirect photons and shadow photons.

We define our model as an agent based decision maker called the oracle algorithm. Decision making optimisations can be seen as "Should I?" questions such as "Should I trace shadow ray?" or "Should I recompute incoming irradiance?". We formalize thoose questions in an Artificial Intelligence procedure.

Avoiding shadow rays The oracle algorithm, answers the question : "Should I trace shadow rays ?". This is used each time we want to trace a shadow ray from a particular position. Practically, the oracle will locate the closest ALA using the ALA kd-tree and compute the difference between the sum of its penumbra sensors p ij (as defined in Table . 1) and the median penumbra sensor value.

p ij -σ m p > 0 → {Y ES, N O}
To achieve more robustness it is possible to use the median response of neighbours agents of i using the proximity network.

Caching Irradiance Our implementation allows to compute an approximation of diffuse indirect illumination at a position x ij . We implement a procedure which compute the mean distance to closests irradiance sensors g ij at positions x g ij and compare it to the median value of this sensors.

x g ij -x ij -σ m i > 0 → {Y ES, N O}
If the procedure answers negatively, indirect illumination is computed by interpolating irradiance values of closest sensors. In application due to homogeneous repartion of ALA few indirect diffuse illumination values need to be computed during the rendering phase.

Results and discussion

In this section we present results based on our implementation of the ALA structure and algorithms. For comparison we implemented a standard photon mapping algorithm with shadow photons, importance sampling and irradiance caching. Images have been rendered on a 2.67GHz Intel Core i7 920 using four cores with a resolution of 1600×1200. The test scene shown in Fig. 4 is a common Cornell box with a mesh and a glass ball. The scene is illuminated by a spherical light source on the ceiling. Fig. 4 shows a rendering of our test scene including only direct lighting, caustics and shadows. We can see that difference is subtle and due uniquely to some aliasing near the silhouette. In Fig. 5 1M photon are cast leading to a 4.9Mb shadow photon map. We observe that the ALA method, while achieving faster rendering, avoids the usage of shadow photons and produces nearly identical shadows with only 50000 ALA. The indirect diffuse illumination is computed using importance sampling on the indirect photon map and irradiance caching using the direct and indirect photon map. The ALA allow to avoid computation of the two photon maps while providing an accurate but perfectible estimation of the indirect illumination in slightly shorter rendering times as shown in Fig. 3.

We observe in Fig. 3 that the local useful data density considering the irradiance criteria is uniformly low in the photon maps. This is explained by the local redundancy of the irradiance data. Our structure show a considerably better density, as the number of sensors perfectly fits with the redundancy of data. This difference leads to an ALA casting phase shorter than the photon casting phase as well as a significant reduction of memory usage as show in 2.

Conclusions and future work

In this paper, we stated that photon density is not necessarily related to useful data for some light estimation. We therefore propose a new estimation of the local density of useful data we call LUDD for local useful data density. We then presented a new agent-based model to discover efficiently a scene, in the sense of matching correctly LUDD. To implement the method we introduced a new structure : Autonomous Lighting Agents (ALA), which are an efficient container to handle data about the scene. We show how the ALA network can be used as an oracle or estimation to lower rendering times and reduce memory overhead in the photon mapping algorithm.

ALA concept and networks could be extended to answer other critical question related to rendering, like halting or not path tracing recursion or handle importance sampling in each intersection related to the whole light distribution in the scene. ALA network could also force photon casting in importance direction, leading light path. We plan also to improve greatly irradiance estimation with ALA networks to correctly replace irradiance caching. Finally, this approach of concentrating efforts on local area, where useful data can be found, could be tested on other kind of vizualisation of 3D data or on other global illumination algorithm like Metropolis Light Transport.

Fig. 1 .

 1 Fig. 1. A. Our test scene with classical indirect photon map (B. middle top) and ALA map (D. middle bottom) and corresponding local useful data density (respectively C. top right for classical photon map and D. bottom right for ALA map); In false color, red means high useful data when blue means low

1 d

 1 Number of intersections with surfaces along the ray between p sensor location and light source divided by distance to source light k , d k being the distance to the closest intersection from sensor to any surface in a arbi-

Fig. 3 .

 3 Fig. 3. Local useful data density for the indirect on ALA (left) and on classical the photon map (right)

Fig. 4 .

 4 Fig. 4. Rendering test scene without indirect diffuse lighting with standard photon mapping (left), ALA photon mapping (center) and difference multiplyed by 5 between two rendered image (right)

Fig. 5 .

 5 Fig. 5. Rendering test scene with standard photon mapping (left), ALA photon mapping (center) and difference multiplyed by 5 between two rendered image (right)

Table 2 .

 2 Comparative results with ALA in memory occupation and casting times

	Map	Total	Memory
		casting times (s) occupation (mb)
	Direct PhotonMap		18
	Indirect PhotonMap Caustic PhotonMap	46	9.4 0.004
	Shadow PhotonMap		4.9
	ALA Network	27	5.2

Table 3 .

 3 Comparative results with ALA in rendering times

Acknowledgements

We would like to thank the DuranDuboi R&D team for their support and for providing the rendering infrastructure.


Benjamin Raynal
email: raynal@univ-mlv.fr

Michel Couprie
email: coupriem@esiee.fr

Venceslas Biri
email: biri@univ-mlv.fr

Homeomorphic Alignment of Weighted Trees

Keywords: Graphs, homeomorphism, alignment, matching algorithm

HAL is

Introduction

Motivation

Motion capture without markers is a highly active research area, as shown by Moeslund et al. [START_REF] Moeslund | A survey of advances in vision-based human motion capture and analysis[END_REF]: between 2000 and 2006, more than 350 papers on this topic were published. Motion capture is used in several applications which have different aims and constraints: 3D models animation, for movies FX or video games for example, requests an highly accurate model, but does not need real-time computation (offline video processing is acceptable).

Real-time interaction, for virtual reality applications, requests a fast computation, at the price of a lower accuracy.

This paper is placed in the context of real-time interaction. Moeslund et al. describe in a previous work [START_REF] Moeslund | A survey of computer vision-based human motion capture[END_REF] the different steps of motion capture.

The first step (called initialization step) consists of finding the initial pose of the subject, represented here by a 3d shape (visual hull) constructed using a multi-view system with an algorithm of Shape From Silhouette [START_REF] Laurentini | The visual hull concept for silhouette-based image understanding[END_REF]. From the right to the left: the original 3D shape, its skeleton, the data tree computed from the skeleton, and the pattern tree. The numbers represent the weights of the edges (i.e. the geodesic distances in the skeleton for the data tree, and the expected ones for the pattern tree).The grey dashed lines represent the expected alignment between the data tree and the model.

crotch, hands and feet), and each edge represents the link between this parts, associated to a weight, representing the distance between two parts. It can be seen as a very simplified kinematic structure, without definition of degrees of freedom.

Our approach to find the initial pose of the subject is therefore to find the best alignment between the pattern tree and the data tree, that is to say, the matching which involves as few modifications as possible, to transform the data tree in the pattern tree (see Fig. 1 for an example of our complete pipeline, and expected alignment).

As we know which part of the shape correspond to each vertex of the pattern tree, and which are the 3D coordinates of each vertex of the data tree, the alignment will give the position of each part of the 3D shape.

Problems

Several kinds of noise and deformities can appear in the data tree:

Ghosts limbs: Due to the reconstruction from silhouettes, parts of space cannot be carved, resulting in "limbs" of the object which do not exist on the real model. Our method must be accurate enough to distinguish these ghosts limbs from real ones.

Spurious branches: Due to the skeletonization algorithm and to the amount of noise of the shape surface, branches without important topological signification can appear on skeleton. For example, in Fig. 1, the edges {g, h}, {l, m}, {i, j}, {j, k} are spurious branches. Our method must be robust enough to work on data trees with consequent amount of spurious branches.

Useless vertex: Vertices with exactly two neighbors are not useful to describe the topology of a shape, and then uselessly split an edge (and its weight) into two parts, making difficult a good matching. This kind of vertices can appear when removing spurious branches or ghosts limbs. For example, in Fig. 1, vertices j, k, m are useless after spurious branches deletion. Our method must be able to match two edges joined by this kind of vertex, with a unique edge.

Splitted vertex: Vertices with more than three neighbors in the pattern tree can correspond to a cluster of vertices linked by weakly weighted edges in the data tree, due to the skeletonization algorithm. For example, in Fig. 1, vertex T of pattern tree matches with vertices b and e in data tree.

Our method must be able to match them.

Our contribution

Approaches found in the literature (see Sect. 3) do not permit to achieve a robust matching, with respect to all these perturbations. However, some existing edit-based distance, the alignment distance [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF], is an interesting way to solve the problem of splitted vertices, and it preserves the topology during the matching. In addition, an operation described in [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF], the cut operation, is specially designed for considering only a subpart of the tree. This operation is exactly what we need to avoid the problems of ghost limbs and spurious branches, which can be considered as useless parts of the data tree.

The problem of useless vertices cannot be solved by methods found in the literature. We introduce in this paper a new kind of alignment, which solves this problem by considering the homeomorphism between trees instead of the isomorphism.

This paper is organized as follow: In Sect. 2, we give the basic definitions and notations related to edge-weighted graphs that will be used in the sequel. In Sect. 3, we determine which edit-based distance of the literature is the most appropriate to base our method. In Sect. 4, we introduce our main contribution, the homeomorphic alignment. We also introduce algorithms to compute it efficiently for rooted trees, as well as unrooted trees. In Sect. 5, we show how to use the cut operation [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF] with homeomorphic alignment. Finally, in Sect. 6, we show the results of different experimentations, and a comparison between our homeomorphic alignment distance and the classical alignment distance.

Basics notions and notations

Undirected graphs

An undirected graph is a pair (V, E), where V is a finite set, and E a subset of {{x, y}, x ∈ V, y ∈ V, x = y}. An element of E is called an edge, an element of V is called a vertex. If {x, y} ∈ E, then x and y are said to be adjacent or neighbors. The set of all neighbors of x is denoted by N (x). The number of vertices adjacent to a vertex v is called the degree of v, and is denoted by deg(v). Let G = (V, E) be an undirected graph, and let x, y be in V , a path

from x to y in G is a sequence of vertices v 0 , ..., v k such that x = v 0 , y = v k and {v i-1 , v i } ∈ E, 1 ≤ i ≤ k.
The number k is called the length of the path. If k = 0 the path is called a trivial path. The path is closed if x = y. The path is simple when no vertex occurs more than once in the sequence of vertices of the path (except possibly x = y). A non-trivial simple closed path in which all edges are distinct is called a cycle. A graph is connected if for all {x, y} ⊂ V , a path from x to y exists in G. A tree is a connected graph with no cycles. A simple path from x to y in a tree is unique and is denoted by π(x, y). A graph with no cycles is called a forest, each of its connected components being a tree.

Directed graphs

A directed graph is a pair (V, A), where V is a finite set, and A is a subset of V × V . An element of A is called an arc, an element of V is called a vertex. Let G = (V, A) be a directed graph, and let x, y be in V , a path from x to y in G is a sequence of vertices s 0 , ..., s k such that

x = v 0 , y = v k and (v i-1 , v i ) ∈ A, 1 ≤ i ≤ k. The undirected graph associated to G is the undirected graph G ′ = (V, E), such that {x, y} ∈ E if and only if (x, y) ∈ A or (y, x) ∈ A. A vertex r ∈ V is a root of G if for all x ∈ V \ {r}, a path from r to x in G exists.
The graph G is antisymmetric if for all (x, y) ∈ A such that x = y, (y, x) / ∈ A. The graph G is a rooted tree (with root r) if r is a root of G, G is antisymmetric and if the undirected graph associated to G is a tree. A graph, where each of its connected components is a tree, is called a rooted forest.

Let G = (V, A) be a rooted tree. If (y, x) ∈ A, we say that y is the parent of x (denoted by par(x)), and that x is a child of y. The set of all children of y is denoted by C(y). The maximum length of a path between the root and any node is called the height of the tree. The vertices on the path from the root to a vertex x are called the ancestors of x. We denote the set of the ancestors of x by anc(x).

Common definitions

Unless otherwise indicated, all the other definitions and notations in this paper are similar for the two kinds of graphs. We give them for directed graphs, the versions for undirected graphs can be obtained by replacing arcs by edges.

Two graphs G = (V G , A G ) and G ′ = (V G ′ , A G ′ ) are said to be isomorphic if there exists a bijection f :

V G → V G ′ , such as for any pair (x, y) ∈ V G × V G , (x, y) ∈ A G if and only if (f (x), f (y)) ∈ A G ′ .
A weighted graph is a triplet (V, A, ω), where V is a finite set, A a subset of V × V , and ω a mapping from A to R. In a weighted tree, the weight of the unique path from x to y, denoted by ω(x, y), is the sum of the weights of all arcs traversed in the path.

In this paper, we say that two weighted graphs (V, E, ω) and (V ′ , E ′ , ω ′ ) are isomorphic whenever the graphs (V, E) and (V ′ , E ′ ) are isomorphic (regardless of the weights).

The merging is an operation that can be applied only on arcs sharing a 2degree vertex. The merging of two arcs (u, v) and (v, w) in a weighted graph G = (V, A, ω) consists of removing v in V , replacing (u, v) and (v, w) by (u, w) in A, weighted by ω((u, w)) = ω((u, v)) + ω((v, w)). [START_REF] Whitney | A set of topological invariants for graphs[END_REF] if there exists an isomorphism between a graph obtained by mergings on G and a graph obtained by mergings on G ′ . See Fig. 2 for some examples of isomorphism and homeomorphism.

Two weighted graphs

G = (V G , A G , ω G ) and G ′ = (V G ′ , A G ′ , ω G ′ ) are home- omorphic
U 1 U 2 U 3 U 4 D 1 D 2 D 3 D 4 D 5

Edit-based distances

The problem of comparing graphs (in particular trees) occurs in diverse areas such as computational biology, image analysis and structured databases. However, the graphs considered in these domains are most often with labeled vertices. Here, each notion will be introduced in the case of graphs with weighted edges/arcs.

In this section, we review different edit-based distances proposed in the literature, and select the one that is the best adapted our aims.

Edit operations

A lot of methods have been developed for tree comparison, as the maximal common subtree [START_REF] Pelillo | Many-to-many matching of attributed trees using association graphs and game dynamics[END_REF][START_REF] Torsello | Polynomial-time metrics for attributed trees[END_REF], or the approximate subtree homeomorphism [START_REF] Pinter | Approximate labelled subtree homeomorphism[END_REF]. An approach widely used to compare trees is to search for a sequence of simple primitive operations (called edit operations) that transforms a tree into the other and that has a minimal cost.

For a graph G = (V, A, ω), classical edit operations are:

Resize: Change the weight of an arc a = (u, v) ∈ A.

Delete: Delete an arc a = (u, v) ∈ A and merge u and v into one vertex.

Insert: Split a vertex in two vertices, and link them by a new arc. See Fig. 3 for illustrations.

The cost of these edit operations is given by a cost function γ(w, w ′ ), where w (respectively w ′ ) is the total weight of the arcs involved in the operation before (respectively, after) its application. As a consequence, the cost of a resizement can be denoted by γ(w, w ′ ), where w is the former weight of the resized arc and w ′ is the new weight, the cost of a deletion can be denoted by γ(w, 0), where w is the weight of the deleted arc, and the cost of an insertion, γ(0, w), where w is the weight of the created arc. Furthermore, we assume that γ is a metric. Typically, γ(w, w

′ ) = |w -w ′ | or (w -w ′ ) 2 .

Choice of edit-based distance

Various edit-based distances have been defined, using different constraints on the edit operations order or on the resulting correspondence between the vertices sets, and different definitions of operations. These edit-based distances can be classified, as proposed by Wang et al. [START_REF] Wang | Finding similar consensus between trees: an algorithm and a distance hierarchy[END_REF]: edit distance [START_REF] Tai | The Tree-to-Tree Correction Problem[END_REF], alignment distance [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF][START_REF] Jansson | A fast algorithm for optimal alignment between similar ordered trees[END_REF], isolated-subtrees distance [START_REF] Tanaka | The tree-to-tree editing problem[END_REF][START_REF] Valiente | An efficient bottom-up distance between trees[END_REF][START_REF] Zhang | Algorithms for the constrained editing distance between ordered labeled trees and related problems[END_REF], and top-down distance [START_REF] Selkow | The Tree-to-Tree Editing Problem[END_REF][START_REF] Yang | Identifying Syntactic differences Between Two Programs[END_REF][START_REF] Wang | Identifying approximately common substructures in trees based on a restricted edit distance[END_REF].

Edit distance

Let G 2 be the graph that results from the application of an edit operation s to graph G 1 ; this is written G 1 ⇒ G 2 via s. Let S be a sequence s 1 , s 2 , ..., s k of edit operations. We say that S transforms graph G to graph

G ′ if there is a sequence of graphs G 0 , G 1 , ..., G k such that G = G 0 , G ′ = G k and G i-1 ⇒ G i via s i for 1 ≤ i ≤ k.
The cost of the sequence S, denoted by γ(S), is the sum of costs of the constituent edit operations. The distance from G to G ′ , denoted by δ(G, G ′ ), is the minimum cost of all sequences of edit operations taking G to G ′ .

For our purpose, this kind of edit-based distance cannot be used, because the associated matching does not preserve topological relations between trees. In addition, the algorithm based on this distance is the one with the highest time complexity.

Isolated-subtrees distance

Isolated-subtrees distance is an edit-based distance such that two disjoint subtrees in the pattern tree will always be matched with two disjoint subtrees in the data tree.

Top-down distance

Top-down distance is an edit-based distance such that an arc (p, p ′ ) in the pattern tree can match an arc (d, d ′ ) in the data tree, only if (par(p), p) matches (par(d), d).

Isolated-subtrees distance and top-down distance cannot always match all the model tree, but only subparts, most often unconnected. However, we will see in the next subsection that it is not the case for alignment distance.

Alignment distance

In [START_REF] Jiang | Alignment of trees -an alternative to tree edit[END_REF], Jiang et al. propose a similarity measure between vertex-labeled trees, that we transpose here for edge-weighted graphs.

Let

G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , A 2 , ω 2 ) be two weighted graphs. Let G ′ 1 = (V ′ 1 , A ′ 1 , ω ′ 1 ) and G ′ 2 = (V ′ 2 , A ′ 2 , ω ′
2 ) be weighted graphs obtained by inserting arcs weighted by 0 (zero) in G 1 and G 2 , such that there exists an isomorphism

I between (V ′ 1 , A ′ 1 ) and (V ′ 2 , A ′ 2 ). The set of all couples of arcs A = {(a 1 , a 2 ); a 1 ∈ A ′ 1 , a 2 ∈ A ′ 2 , a 2 = I(a 1 )} is called an alignment of G 1 and G 2 .
The cost C A of A is the sum of the costs of all operations used to align G 1 and G 2 : the insertions of arcs weighted by 0 (zero) (which is free, due to non-variation of weights), and the resizement for each arc a 1 ∈ A ′ 1 to the weight of I(a 1 ). More formally :

C A = (a1,a2)∈A γ(ω ′ 1 (a 1 ), ω ′ 2 (a 2 )) .
The minimal cost of all alignments from G 1 and G 2 , called the alignment distance, is denoted by α(G 1 , G 2 ). The alignment distance is a special case of the edit distance, where any insertion occurs before any deletion. As a result,

α(G, G ′ ) ≥ δ(G, G ′ ).
Alignment distance is interesting in our case for three reasons: it takes into account topological relations between trees, it can be computed in polynomial time, and it enables to "remove edges", regardless of the rest of the graph, solving the problem of splitted vertices.

Illustration

We consider only the problem of splitted vertices in this example. In Fig. 4, we show the way to obtain an alignment with minimal cost: only one edge insertion in the pattern tree, on the vertex T , representing the torso.

The resulting alignment, represented by dotted lines, has a cost of 12. The corresponding vertices matching is the following: H matches with a, T with {b, e}, A 1 with d, A 2 with c, C with f , F 1 with n, and F 2 with o.

H T C A 2 A 1 F 1
G 1 G ′ 1 G ′ 2 G 2 Figure 4: Example of alignment: G ′ 1 (resp. G ′ 2 ) is obtained from G 1 (resp. G 2 ) by insertions of edges (here G ′ 2 = G 2 ).
The dotted lines represent one of the possible alignments of G ′ 1 and G ′ 2 .

Homeomorphic alignment distance

If we now take into account the problem of useless 2-degree vertices, we can see that the alignment strategy does not work anymore (see Fig. 5). In this example, the split of the "leg" edges of the data tree in several parts involves that the best alignment is obtained by matching the "arms" of the pattern (edges {T, A 1 } and {T, A 2 }) with the parts of the legs in the data tree which are closest of them, in regard of their weights (edges {m, o} and {j, n}). The cost of the alignment, shown by the dotted lines, is equal to 26 and is minimal.

For the purpose of solving the useless vertex problem, we propose a new alignment strategy.

f 1 f 2 0 A 2 A 1 T 1 T 2 T 3 T 4 H C 1 C 2 C 3 F 1

Merging kernel

Considering that a merging on a vertex v on the graph G = (V, A, ω) does not affect the degree of any vertex in V \ {v} (by definition of merging operation) and therefore the possibility of merging this vertex, the number of possible mergings decreases by one after each merging. In consequence, the maximal size of a sequence of merging operations, transforming G into another graph G ′ = (V ′ , A ′ , ω ′ ) is equal to the initial number of possible mergings in G. It can be remarked that any sequence of merging operations of maximal size yields the same result. The graph resulting of such a sequence on G is called the merging kernel of G, and is denoted by MK(G).

The following proposition is straightforward:

Proposition 1. Two graphs G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , G 2 , ω 2 ) are home- omorphic iff MK(G 1 )
and MK(G 2 ) are isomorphic.

Homeomorphic alignment distance

Let

G 1 = (V 1 , A 1 , ω 1 ) and G 2 = (V 2 , A 2 , ω 2 ) be two weighted graphs. Let G ′ 1 = (V ′ 1 , A ′ 1 , ω ′ 1 ) and G ′ 2 = (V ′ 2 , A ′ 2 , ω ′
2 ) be weighted graphs obtained by deleting arcs in G 1 and G 2 , such that there exists an homeomorphism between

G ′ 1 and G ′ 2 (not necessarily unique). Let G ′′ 1 = (V ′′ 1 , A ′′ 1 , ω ′′ 1 ) and G ′′ 2 = (V ′′ 2 , A ′′ 2 , ω ′′ 2 ) be the merging kernel of G ′
1 and G ′ 2 , respectively. From proposition 1, there exists an isomorphism I between G ′′ 1 and G ′′ 2 . The set of all couples of arcs 6). The graph G ′′ 1 is called the left graph of H. The graph G ′′ 2 is called the right graph of H. The cost C H of H is the sum of the costs of all operations used to homeomorphically align G 1 and G 2 : the deletion of arcs in G 1 and G 2 , to obtain G ′ 1 and G ′ 2 respectively, and the resizement for each arc a 1 ∈ A ′′ 1 to the weight of H(a 1 ). More formally :

H = {(a, a ′ ); a ∈ A ′′ 1 , a ′ ∈ A ′′ 2 , a ′ = I(a)} is called an homeomorphic alignment of G 1 with G 2 (see figure
C H = (a,a ′ )∈H γ(ω ′′ 1 (a), ω ′′ 2 (a ′ )) + a d ∈A1\A ′ 1 γ(ω 1 (a d ), 0) + a ′ d ∈A2\A ′ 2 γ(0, ω 2 (a ′ d )) .
This minimal cost of all homeomorphic alignments between G 1 and G 2 , called the homeomorphic alignment distance, is denoted by η(G 1 , G 2 ).

Illustration

Let us consider again the pattern and data trees used in Fig. 5. In Fig. 6 is shown the way to obtain a homeomorphic alignment with minimal cost: only one edge deletion in the data tree (edge {b, e}), is necessary to obtain a homeomorphism. The merging kernel of

G ′ 1 is equal to G ′ 1 ,

Algorithm for rooted trees

In order to compute homeomorphic alignment distance with good complexity, we will use a bottom-up approach. This approach, widely used in the literature on edit-based distances, is based on a reformulation of the distance between trees in terms of both the distance between their subtrees and the matching of their roots. Then, using this reformulation, the computation starts by the distances between the simplest subtrees (i.e the leafs), and goes up to the root. The distances between subtrees are kept in memory, avoiding redundancy of computation.

However, there are two main differences to take into account in the case of homeomorphic alignment distance: the fact that information is on arcs, rather than of vertices, and the fact that an arc can be matched with a set of arcs merged together (due to merging kernel application).

The first problem involves to take into consideration several arc weights in the reformulation (one for each arc between the root and its children), instead of the unique value linked to the root, as in the literature.

The second problem involves to take into account not only the subtrees of the tree, but also those which can be generated by some merging.

Our approach to solve these two problems is to use a special kind of tree, the root of which has only one child. Then, we make a reformulation of the homeomorphic alignment distance between this kind of trees, in function of subtrees with the same property. By this way, there is only one new arc to take into consideration at each step, the one between the root and its child. Moreover, we consider that this arc can be the result of a merging. In combination with the recursive nature of the bottom-up approach, it will lead to consider subtrees with all possible mergings. Finally, we also need a reformulation of the homeomorphic alignment distance between "classic" trees in function of the one between the subtrees described above.

Definitions and notations

First, we need to define more formally the particular kind of tree described above.

Let T = (V, A, ω) be a weighted tree rooted in r T . We denote by T (v), v ∈ V , the subtree of T rooted in v (see Fig. 7). Let v a be an ancestor of v, we denote by T cut (v, v a ) the subgraph of T obtained from T (v a ) by removing all complete subtrees which do not contain vertices of T (v).

We denote by T (v, v a ) the tree obtained from T cut (v, v a ) by merging on each vertex n ∈ anc(v) \ {v a , v}. We say that T (v, v a ) is the subtree of T rooted in v a pruned in v, and we call this kind of trees a pruned tree (see Fig. 7 for an example).

We denote by F(T, v) the pruned forest, the connected components of which are the trees T (p, v), for all p ∈ C(v) (see Fig. 7 for an example). By abuse of notation we also denote by F(T, v) the set of all connected components of this forest.

Reformulations

In order to use a bottom-up approach, we have to express the homeomorphic alignment distance:

• between trees in function of the distances between their pruned forests (Prop. 2).

• between pruned trees in function of the distances between their pruned subtrees (Prop. 4).

• between pruned forests in function of the distances between their pruned subtrees (Prop. 5). In addition, we have to reformulate these distances in the special cases where at least one of the trees is empty (Prop. 3).

In the sequel we consider two weighted trees P = (V P , A P , ω P ) and D = (V D , A D , ω D ), rooted respectively in r P and r D .

Proposition 2. We have :

η(P, D) = η(F(P, r P ), F(D, r D )).
Proof. Since the root of a tree cannot be eliminated by any operation (merging, deletion) involved in the definition of η, it may be seen that any homeomorphic alignment H of P with D has a corresponding alignment H ′ of F(P, r P ) with F(D, r D ) of same cost, and the converse also holds. Proposition 3. Let i ∈ V P \ {r P }, j ∈ V D \ {r D }, i a ∈ anc(i), j a ∈ anc(j), η(∅, ∅) = 0 η(P (i, i a ), ∅) = η(F(P, i), ∅) + γ(ω P (i a , i), 0) η(F(P, i a ), ∅)

= i ′ ∈C(ia) η(P (i ′ , i a ), ∅) η(∅, D(j, j a )) = η(∅, F(D, j)) + γ(0, ω D (j a , j)) η(∅, F(D, j a )) = j ′ ∈C(ja) η(∅, D(j ′ , j a )). Proof. Straightforward. Proposition 4. Let i ∈ V P \ {p}, j ∈ V D \ {d}, i a ∈ anc(i), j a ∈ anc(j). η(P (i, i a ), D(j, j a )) = min       
η(P (i, i a ), ∅) + η(∅, D(j, j a )), γ(ω P (i a , i), ω D (j a , j)) + η(F(P, i), F(D, j)), min jc∈C(j) {η(P (i, i a ), D(j c , j a )) + η(∅, F(D, j) \ D(j c , j))}, min ic∈C(i) {η(P (i c , i a ), D(j, j a )) + η(F(P, i) \ P (i c , i), ∅)}.

Proof. Let H be an homeomorphic alignment of P (i, i a ) with D(j, j a ), and let

H L = (V L , A L , ω L ) and H R = (V R , A R , ω R )

the left and right graphs of H, respectively.

There are seven possible cases:

1. H is an empty set (it is cheaper to remove both P (i, i a ) and D(j, j a ) than to align them). 2. {(i a , i), (j a , j)} ∈ H. 3. ∃f ∈ A R , f being obtained by merging (j a , j) with other arcs. In this case, there is one and only one child j c of j, such (j, j c ) is merged with (j a , j) in f , and then all D(j ′ c , j), j ′ c ∈ C(j) \ {j c } are deleted. 4. ∃f ∈ A L , f being obtained by merging (i a , i) with other arcs. In this case, there is one and only one child i c of i, such (i, i c ) is merged with (i a , i) in f , and then all P (i ′ c , i), i ′ c ∈ C(i) \ {i c } are deleted. Cases 1,2,3,4 justify, respectively, the lines 1,2,3,4 of the expression of η(P (i, i a ), D(j, j a )) in the proposition.

The three last cases cannot lead to a better homeomorphic alignment:

5. The deletion of (i a , i) and (j a , j) cannot be preferred to the resizement (possible case 2), because γ(ω P (i a , i), 0) + γ(0, ω D (j a , j)) ≥ γ(ω P (i a , i), ω D (j a , j)). 6. If (i a , i) was deleted, and not (j a , j), then only one

P (i c , i), i c ∈ C(i) is
aligned with D(j, j a ), the other being removed. It is less expensive to merge (i a , i) with (i, i c ) (possible case 3), because γ(ω P (i a , i), 0) + γ(ω P (i, i c ), ω D (j a , j)) ≥ γ(ω P (i a , i c ), ω D (j a , j)). 7. The deletion of (j a , j) is more expensive than the merging of (j a , j) (possible case 4), for the same reasons as above.

Proposition 5. ∀A ⊆ F(P, i), B ⊆ F(D, j),

η(A, B) = min                        min D(j ′ ,j)∈B {η(A, B \ {D(j ′ , j)}) + η(∅, D(j ′ , j))}, min P (i ′ ,i)∈A {η(A \ {P (i ′ , i)}, B) + η(P (i ′ , i), ∅)}, min P (i ′ ,i)∈A,D(j ′ ,j)∈B {η(A \ {P (i ′ , i)}, B \ {D(j ′ , j)}) +η(P (i ′ , i), D(j ′ , j))}, min P (i ′ ,i)∈A,B ′ ⊆B {η(A \ {P (i ′ , i)}, B \ B ′ ), +η(F(P, i ′ ), B ′ ) + γ(ω P (i, i ′ ), 0)}, min A ′ ⊆A,D(j ′ ,j)∈B {η(A \ A ′ , B \ {D(j ′ , j)})+ η(A ′ , F(D, j ′ )j) + γ(0, ω D (j, j ′ ))}.
Proof. Let H be an homeomorphic alignment of A ⊆ F(P, i) with B ⊆ F(D, j), and let P (i ′ , i) ∈ A and D(j ′ , j) ∈ B. There are five possible cases:

1. D(j ′ , j) is not aligned with element of A, 2. P (i ′ , i) is not aligned with element of B, 3. P (i ′ , i) is aligned with D(j ′ , j), 4. disjoint subparts of P (i ′ , i) are aligned with elements of B, 5. disjoint subparts of D(j ′ , j) are aligned with elements of A.

Cases 1, 2, 3, 4, 5 justify, respectively, the lines 1, 2, 3, 4, 5 of the expression of η(A, B) in the proposition.

Algorithm

Using the reformulations defined above, we can now design our algorithm following the bottom-up approach. We denote by N the maximum degree of a vertex in the pattern tree or in the data tree, H the maximal height of both trees and S the maximal size (number of vertices) of both trees. We recall that a set of size n has 2 n subsets.

The time complexity of computing each reformulation used in the algorithm is:

• η(P (i, i a ), ∅) and η(∅, D(j, j a )) → O(1).

• η(A, ∅) and η(∅, B) → O(N ).

• η(P (i, i a ), D(j, j a )) → O(N ).

• η(A, B) → O(N * 2 N ).
Combining these complexities with the different loops of the algorithm, we obtain that the total computation is in O(S 2 * (N * 2 3 * N +H 2 * N )) time complexity.

If the maximal degree is bounded, the total computation is in O(S 2 * H 2 ) time complexity.

Algorithm for unrooted trees

First, let us give an expression of the homeomorphic alignment distance between unrooted trees, in function of the distances between all their possible rooted versions.

Let G = (V, E, ω) be a weighted tree, let r ∈ V , we denote by G r , the directed weighted tree rooted in r, such that G is the undirected graph associated to G r (see Fig. 8). Proposition 6. Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. We have: η(P, D) = min i∈VP ,j∈VD {η(P i , D j )}.

Proof. Let G = (V, E, ω) be a graph, and r ∈ V a vertex of G. Notice that:

• a merging occurring on a vertex v ∈ V \ {r} in G can occur in G r ,
• deletion, insertion, resizement, and division occurring in G can occur in G r .

On the other hand, for each optimal homeomorphic alignment H of P in D, it is easy to see that there exists p ∈ V P and d ∈ V D , such that p and d are not affected by a merging. For example, if

D ′ = (V ′ D , E ′ D , ω ′ D )
is a subgraph such as η(P, D) = η(P, D ′ ), p and d can be chosen as 1-degree vertices of V P and V ′ D , respectively. As a result, η(P, D) = η(P p , D d ). Since the homeomorphic alignment is more constrained in the case of rooted trees than in the case of unrooted trees, we can assure than η(P, D) ≤ η(P a , D b ), a ∈ V P , b ∈ V D . To sum up, knowing that η(P, D) ≤ η(P a , D b ), a ∈ V P , b ∈ V D , and that there exists p ∈ V P and d ∈ V D , such that η(P, D) = η(P p , D d ), we conclude that η(P, D) = min i∈VP ,j∈VD η(P i , D j ).

Naive algorithm

Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. From proposition 5, we propose a first, naive algorithm to compute η(P, D), consisting of computing the homeomorphic alignment distance for all couples of weighted rooted trees we can obtain from P and D, and keeping the minimum reached.

Complexity. P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. As the number of rooted trees we can obtain from an undirected tree is equal to the number of vertices of this graph, the number of couples of weighted rooted trees we can obtain from P and D is equal to

|V P | * |V D |.
The total computation is then in

O(|V P | 2 * |V D | 2 * (2 dP * 2 dD * (d D * 2 dP + d P * 2 dD ) + h P * h D * (d P + d D )
)) time complexity, where h P (respectively, h D ) is the maximal height of a rooted tree obtained from P (respectively, D).

If the maximal degree is bounded, the total computation is in

O(|V P | 2 * |V D | 2 * h P * h D ) time complexity.

Optimized algorithm

It is easy to see that the above algorithm computes the alignment of subparts of P and D more than one time. Using dynamic programming and an adapted order of navigation in the tree, we can avoid useless computation.

Let P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two undirected weighted trees.

We denote by F(P, a, b), a, b ∈ V P the set of rooted trees P r , r ∈ V P , such that b is an ancestor of a in P r . We denote by anc(P, a, b), a, b ∈ V P the set of vertices x ∈ V P such that x is an ancestor of a in at least one rooted tree in F(P, a, b). We denote by C(P, a, b), a, b ∈ V P the set of of vertices x ∈ V P such that x is a child of a in at least one rooted tree in F(P, a, b).

It is easy to see that for computing η(P p (i, i a ), D d (j, j a )) and η(F(P p , i), F(D d , j)), we need to know η(P p (i c , i), D d (j c , j)) for all i c ∈ C(P, i, p), j c ∈ C(D, j, d). We can start by computing η(P p (i, i a ), D d (j, j a )) and η(F(P p , i), F(D d , j)), for all i (respectively, j) being a leaf of P p (respectively, D d ), which have no child, by definition, and continue iteratively with all vertices which have all their children already computed.

An adapted order of navigation for a tree T = (V, E, ω) can be obtained by the following algorithm.

Algorithm 2: Order of navigation computation algorithm (compute-Order)

Data: unrooted tree T Result: list of couples of vertices L begin list of couples of vertices

L ← ∅ LIFO queue of vertices Q ← ∅ foreach v ∈ V do if deg(v) = 1 then push v in Q while Q = ∅ do pop v of Q foreach w ∈ N (v) do if (v, w) / ∈ L then if ∀x ∈ N (v) \ {w}, (x, v) ∈ L then push w in Q add (v, w) at the end of L end
Complexity. Let T = (V, E, ω) be an unrooted tree. Each vertex of degree 1 (one) will be put in the queue, then, for each edge, each of its vertices will be put twice in the queue. As |V | = |E| + 1, the complexity of computeOrder is in O(|E|).

Final algorithm. We can compute Homeomorphic Alignment for unrooted trees with improved complexity, using this order of navigation. Combining the time complexities of computing the reformulations (given in the Sec.4.3.4) with the loops of the algorithm, and taking into account the above observation, we obtain that the total computation time of this algorithm is in O(S 2 * (N * 2 3 * N + S 2 * N )) complexity, with N and S defined as in Sec. 4.3.4. If the maximal degree is bounded, the total computation is in O(S 4 ) time complexity.

Cut operation

For the purpose of removing spurious branches without any cost, we propose to integrate the cut operation in our alignment.

In [START_REF] Wang | Finding approximate patterns in undirected acyclic graphs[END_REF], Wang et al. propose a new operation allowing to consider only a part of a tree. Let G = (V, A, ω) be a weighted tree. Cutting G at an arc a ∈ A, means removing a, thus dividing G into two subtrees G 1 and G 2 . The cut operation consists of cutting G at an arc a ∈ A, then considering only one of the two subtrees. Let K a subset of A. We use Cut(G, K, v) to denote the subtree of G containing v and resulting from cutting G at all arcs in K. In the case of a rooted tree, we consider that the root r G of G cannot be removed by the cut operation, and then we use the notation Cut(G, K) = Cut(G, K, r G ). In the case of a rooted forest, we consider that the root of each rooted tree composing the rooted forest cannot be removed by the cut operation, and then we use the same notation than above: Cut(G, K).

See Fig. 9 for some examples of cut operation. Our main problem can be stated as follows: Given a weighted tree P = (V P , A P , ω P ) (the pattern tree) and a weighted tree G D = (V D , A D , ω D ) (the data tree), find η cut (P, D) = min K⊆AD,v∈VD {η(P, Cut(D, K, v)} and the associated homeomorphic alignment. In the case of rooted trees and rooted forests, η cut (P, D) = min K⊆AD {η(P, Cut(D, K))}.

Integration of cut operation in our algorithm

It can be seen that the cut operation can be integrated in the homeomorphic alignment by replacing the deletion of complete subtrees by their cut.

In order to obtain the algorithms for the computation of η cut , we just need to replace:

• η by η cut in the reformulations and algorithms of the Subsec. 4.3 and 4.4

• Prop. 3 by the following:

Proposition 7. Let i ∈ V P \ {p}, j ∈ V D \ {d}, i a ∈ anc(i), j a ∈ anc(j), η cut (∅, ∅) = 0 η cut (P (i, i a ), ∅) = η cut (F(P, i), ∅) + γ(ω P (i a , i), 0) η cut (F(P, i a ), ∅) = i ′ ∈C(ia) η cut (P (i ′ , i a ), ∅) η cut (∅, D(j, j a )) = 0 η cut (∅, F(D, j a )) = 0.
It is interesting to remark that the complexity is not modified.

Experimentation

Usage of homeomorphic alignment

There are several ways to use the homeomorphic alignment:

• If we have no a priori knowledge both on pattern tree P and on data tree D, we need to use the homeomorphic alignment on the two unrooted trees. In this case, the complexity is in O(|V P | 2 * |V D | 2 ).

• If we want to be sure that a specific vertex v in the data tree is aligned with a specific vertex w in the data tree (for example if we are sure than the vertex of the head in the data is the one with the highest z-coordinate), we can use the homeomorphic alignment between P v and D w . In this case, the complexity is in O(|V P | * |V D | * h P * h D ), where h P (respectively, h D ) represents the height of P (respectively, D).

• If we want to be sure that a specific vertex v in the model tree is aligned (i.e. there is no merging on v), we can use the homeomorphic alignment between P v and D, by successively computing the homeomorphic alignment between P v and D w , for all w ∈ V D , and using optimizations as in 4.4.2. In this case, the complexity is in O(|V P | * h P * |V D | 2 ).

• If we want to be sure that a specific vertex v in the data tree is aligned (i.e. there is no merging or cut on v), we can use the homeomorphic alignment between P and D v , using the same method as above. In this case, the complexity is in

O(|V P | 2 * |V D | * h D ).
In our application, consisting to find the initial pose of a subject, we can at least assume that the torso of the subject will be aligned.

Results

Our model tree contains seven vertices, representing head, torso, crotch, the two hands and the two feet. Experimentally, the data tree obtained from the skeleton of the visual hull has a degree bounded by 4, and its number of vertices is between seven and twenty, with a Gaussian probability repartition centered on ten.

All the results have been obtained on a computer with a processor Intel(R) Core(TM) 2 Quad Q8200 at 2.33 GHz and 3 Go of RAM, with Linux, Ubuntu 9.04. The algorithms have been implemented in C++.

Speed

Protocol. To find the average computation time of involved algorithms, we have randomly generated 32 pattern trees, and for each one, 32 data trees, yielding 1024 pairs of trees. Each pattern tree has seven vertices, one of which has a degree equal to 4. Each data tree has at least one 4-degree vertex. The results for both alignment and homeomorphic alignment algorithms, for the different kinds of trees, are shown in figure 10.

Discussion. For our purpose of initializing the pose of a subject, we can see that in the average case (|V d | = 12), the homeomorphic alignment can be computed very quickly (frequency ≥ 100), even in the case of unrooted trees. However, we prefer to use the homeomorphic alignment with a rooted pattern tree, because of its better speed, without loss of precision.

The main observation we can do, in regard of the results, is the weak difference of performances between alignment and homeomorphic alignment in practice, even if there is a significant difference of complexity between alignment and homeomorphic alignment. It is due to small size of the trees involved here, and to the fact that, even if the degree of the trees is bounded, and by this way, can be considered as a constant, the computation of the alignment (or homeomorphic alignment) of subforests take a large part of the computation time. This computation part being similar for both algorithms, the rest (where the complexity differs) is less significant in regard of the speed.

Accuracy

Protocol. Each experiment consists of the following:

1. We randomly generate a pattern tree P . 2. We generate a data tree D from P in two steps: weight variation: we randomly alter the edge weights by a given amount of variation.

structural noise: we randomly add new vertices by three ways, corresponding to the different types of noise: splitting an existing vertex, and linking the two parts by a 0-weighted edge (splitted vertices), adding a new 2-degree vertex by the split of an edge (useless 2-degree vertices), and adding a new 1-degree vertex, linked by an randomly weighted edge (spurious branches and ghost limbs).

3. We compute the alignment and the homeomorphic alignment between P and D. 4. The precision is given by the percentage of good matchings, that is the percentage of pattern tree vertices that match with their equivalent in the The spurious branches do not disturb the good alignment of the model on the data.

The only case of bad alignment has been obtained on a very low quality visual hull, where the length of legs was shorter than the length of arms. This case can be solved by rooting the model tree on the head, and the data tree on the vertex with highest z-coordinate.

H H H H T C A 1 A 2 F 1 F 2 T T T C C C A 1 A 1 A 1 A 2 A 2 A 2 F 1 F 1 F 1 F 2 F 2 F 2

Conclusion

In this paper, we have introduced a new type of alignment between weighted trees, the homeomorphic alignment, taking into account the topology and avoiding the noise induced by spurious branches, splitted and useless 2-degree vertices. We have also developed several robust algorithms to compute it with a good complexity, which enable its application in real time for motion capture purpose.

In future works, we will take into account more useful information on the model, such as spatial coordinates of data vertices, and include them in our algorithm, for a better robustness. Finally, using this alignment, we will propose a new fast method of pose initialization for motion capture applications.

Figure 1 :

 1 Figure 1: Example of data tree acquisition and expected alignment with pattern (model) tree.From the right to the left: the original 3D shape, its skeleton, the data tree computed from the skeleton, and the pattern tree. The numbers represent the weights of the edges (i.e. the geodesic distances in the skeleton for the data tree, and the expected ones for the pattern tree).The grey dashed lines represent the expected alignment between the data tree and the model.

Figure 2 :

 2 Figure 2: First row: U 1 is obtained from U 2 by merging on x. U 4 is obtained from U 3 by merging on y. U 2 and U 3 are isomorphic, U 1 and U 4 are homeomorphic. Second row: D 1 is obtained from D 2 by merging on x. D 5 is obtained from D 4 by merging on y. D 2 and D 4 are isomorphic, D 1 and D 5 are homeomorphic, D 2 and D 3 are not isomorphic (red arcs have different orientations).

4 Figure 3 :

 43 Figure 3: Examples of edit operations: on the top, on undirected graphs. On the bottom, on directed graphs.

2 G 1 Figure 5 :

 215 Figure 5: Example of bad alignment, due to useless 2-degree vertices.

2 Figure 6 :

 26 Figure 6: Example of homeomorphic alignment: G ′ 1 (resp. G ′ 2 ) is obtained from G 1 (resp. G 2 ) by deletions of edges. G ′′ 1 = MK(G ′ 1 ) and G ′′ 2 = MK(G ′ 2 ). The dotted lines represent one of the possible homeomorphic alignments of G 1 and G 2 .

Figure 7 :

 7 Figure 7: From the left to the right: a tree T , a subtree of T rooted in b, a subtree of T rooted in a and pruned in e, and a pruned forest of T with origin b.

Algorithm 1 : 4 end 4

 144 Homeomorphic Alignment Distance for Rooted Trees Data: pattern rooted tree P , data rooted tree D Result: η(P, D) = η(F(P, r P ), F(D, r D )); // Prop.2 begin foreach p ∈ V P , in suffix order do foreach A ⊆ F(P, p) do Compute η(A, ∅); // Prop.3 foreach p a ∈ anc(p) \ {p} do Compute η(P (p, p a ), ∅); // Prop.3 foreach d ∈ V D , in suffix order do foreach B ⊆ F(D, d) do Compute η(∅, B); // Prop.3 foreach d a ∈ anc(d) \ {d} do Compute η(∅, D(d, d a )); // Prop.3 foreach p ∈ V P , in suffix order do foreach d ∈ V D , in suffix order do foreach A ⊆ F(P, p) do foreach B ⊆ F(D, d) do Compute η(A, B); // Prop.5 foreach p a ∈ anc(p) \ {p} do foreach d a ∈ anc(d) \ {d} do Compute η(P (p, p a ), D(d, d a )); // Prop.

Figure 8 :

 8 Figure 8: A tree G and the rooted trees G a and G e .

Algorithm 3 : 1 .

 31 Homeomorphic Alignment Distance for Unrooted Trees Data: pattern rooted tree P , data rooted tree D Result: η(P, D) begin list of couples of vertices L P ← computeOrder(P ) list of couples of vertices L D ← computeOrder(D) foreach (p, p ′ ) ∈ L P do foreach A ⊆ F(P p ′ , p) do Compute η(A, ∅); // Prop.3 foreach p a ∈ anc(P, p, p ′ ) \ {p} do Compute η(P pa (p, p a ), ∅); // Prop.3foreach (d, d ′ ) ∈ L D do foreach B ⊆ F(D d ′ , d) do Compute η(∅, B); // Prop. 2 foreach d a ∈ anc(D, d, d ′ ) \ {d} do Compute η(∅, D da (d, d a )); // Prop.3 foreach (p, p ′ ) ∈ L P do foreach (d, d ′ ) ∈ L D do foreach A ⊆ F(P p ′ , p) do foreach B ⊆ F(D d ′ , d) do Compute η(A,B); // Prop.5 foreach p a ∈ anc(P, p, p ′ ) \ {p} do foreach d a ∈ anc(D, d, d ′ ) \ {d} do Compute η(P pa (p, p a ), D da (d, d a )); // Prop.4 Compute η(P, D); // Prop.2 and Prop.6 end Complexity. P = (V P , E P , ω P ) and D = (V D , E D , ω D ) be two weighted trees. The initialization of L P and L D is in O(|V P | * d P + |V D | * d T ) time complexity. Observe that, for p ∈ V P , p ′ ∈N (p) |anc(P, p, p ′ ) \ {p}| = |V P | -As a result, (p,p ′ )∈LP |anc(P, p, p ′ ) \ {p}| = (|V P | -1) × |V P | .

  {{c, e}}, e) Cut(G, {{c, e}, {h, j}}, e) , {(l, m), (o, p)})

Figure 9 :

 9 Figure 9: Examples of cut operations. Top row: on an unrooted tree. Bottom row: on a rooted forest.

Figure 10 :

 10 Figure 10: Frequencies of alignment and homeomorphic alignment for variable sizes of data tree, for the different rooting cases.

Figure 12 :

 12 Figure 12: Special case where alignment is better than homeomorphic alignment. First column: P is the pattern tree, D is obtained from P by a weight variation of 20% on {a, b}, and by adding a new 1-degree vertex c, linked to b by an edge weighted by 1. Second column: optimal alignment. P ′ = P and D ′ = Cut(D, {{b, c}}, b). Accuracy is equal to 100%. Last column: homeomorphic alignment. P and D are already homeomorphic. P ′′ = P and D ′′ = MK(D). Accuracy is equal to 50%.

Figure 13 :

 13 Figure 13: Examples of results obtained on 3D shapes with different resolutions and noises: green voxels represent the points matching with pattern tree.

  and the merging kernel of G ′2 is obtained by merging on vertices j, h and m.The resulting homeomorphic alignment, represented by dotted lines, has a cost of 12. The corresponding vertices matching is the following: H with a, T with {b, e}, A 1 with d, A 2 with c, C with f , F 1 with n, and F 2 with o.

  Then, we can use the matching with O(|V P | * h P * |V D | 2 ) for the initialization. For the tracking, if a part of the subject is static, we can use the last alignment of this part for obtain a matching in O(|V P | * |V D | * h P * h D ).

data tree (which is known, according to the above protocol to generate the data tree).

The results, obtained by averaging precisions over 1000 experiments, are shown on Fig. 11. Discussion. The first important point which can observed is that homeomorphic alignment give always perfect matchings when only structural noise occurs.

In the other cases, the homeomorphic alignment globally gives a more accurate matching than the alignment, specially in regard of the structural noise. However, in case of high weight variation and low structural noise, the results of alignment are better than those of homeomorphic alignment. It is due to a special case of subtree and weight variation, which affects the accuracy of homeomorphic alignment but not the alignment (see Fig. 12 for an example).

In the case of our application, we empirically observe weight variation between 0 and 33 percent, and structural noise between 0 and 200 percent. For the majority of these cases, homeomorphic alignment give sensibly better results.

Pose Initialization

We have checked the alignment on several types of visual hulls, with different resolutions (64 3 and 128 3 ), with or without spurious branches and ghosts limbs. Some of these results are shown in Fig. 13.

The ghost limb can be matched by the algorithm only if its position on the skeleton is the same as another limb, and if it has approximately the same length. In the other cases, it is successfully removed.


Rogelio Lozano

Mazen Alamir

Jonathan Dumon

Ahmad Hably

Ahmad Hably Control

Control of a wind power system based on a tethered wing

Keywords: Kite, flight dynamics, aerodynamics, embedded control, renewable energies

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

Control of a wind power system based on a tethered wing R. Lozano Jr, M. Alamir, J. Dumon, A.Hably

INTRODUCTION

In the last few years new alternative systems for wind power generation [START_REF] Fletcher | Electricity Generation from Jet-StreamWinds[END_REF][START_REF] Manalis | Airborne Windmills: Energy Source for Communication Aerostats[END_REF][START_REF] Manalis | Airborne Windmills and Communication Aerostats[END_REF][START_REF] Riegler | Tethered Wind Systems for the Generation of Electricity[END_REF][START_REF] Riegler | Transformation of Wind Energy by a High-Altitude Power Plant[END_REF][START_REF] Roberts | Unmanned Rotorcraft to Generate Electricity Using Upper Atmospheric Winds[END_REF][START_REF] Fletcher | OntheRotaryWingConcept for Jet Stream Electricity Generation[END_REF][START_REF] Rye | Longitudinal Stability of a Hovering, Tethered Rotorcraft[END_REF] have appeared which overcome the drawbacks of standard power generation systems like nuclear power, gas turbine generator, hydraulic systems, etc. The wind turbine are certainly playing a major role in creating renewable energies. More recently some alternative wind power systems have been proposed such as those based on tethered kites attached to an electricity generator, see [START_REF] Loyd | Crosswind Kite Power[END_REF][START_REF] Williams | Optimal Cross-Wind Towing and Power Generation with Tethered Kites[END_REF][START_REF] Houska | Optimal Control of Towing Kites. Conference on Control and Decision[END_REF][START_REF] Williams | Modeling and Control of a Kite on a Variable Length Flexible Inelastic Tether[END_REF][START_REF] Williams | Nonlinear Control and Estimation of a Tethered Kite in Changing Wind Conditions[END_REF][START_REF] Lansdorp | Design and Testing of a Remotely Controlled Surfkite for the Laddermill[END_REF][START_REF] Lansdorp | Development and testing of a steering mechanism for Laddermill kites[END_REF][START_REF] Williams | Flexible Tethered Kite with Moveable Attachment Points, Part I: Dynamics and Control[END_REF][START_REF] Williams | Flexible Tethered Kite with Moveable Attachment Points, Part II: State and Wind Estimation[END_REF]. Such systems appear to be competitive because they are much less expensive, require less infrastructure to be installed, and have a much smaller ecological footprint. Nevertheless several aspects of kite systems are still to be studied like the control strategies which will enable the system to operate autonomously. Indeed a kite or a wing subject to wind and attached to the ground with a rope may produce oscillations which might lead to crash.

This paper presents a wind power system based on a tethered kite driving a dynamo on the ground, see figure 1 . We present the aerodynamical equations which model the system and propose a control strategy for stabilizing the angle of flight θ and the length of the rope r. The control objective is to keep the angle of the rope θ fixed to a desired value during an energy producing cycle in which the rope length r follows a periodic trajectory. As the kite moves upwards pushed by the wind, the dynamo generates power. When the kite is moving upwards the wing's angle attack is kept at an angle which produces a high lift force. During the period when the kite is moving downwards the angle of attack is smaller so that almost no energy is used bringing the kite to its lowest altitude. As a consequence, there will be a positive amount of energy generated at every cycle. The control strategy is based on the state feedback linearization technique [START_REF] Isidori | Nonlinear Control Systems[END_REF]. An experiment has been set up including a wind tunnel, a wing, and an embedded control system.

Fig. 1. Flight test for the wind power system in the GIPSA wind tunnel (wind speed : 8m/s).

This paper is organised as follows: Section II presents the model of the wind power generation system. The energy production is explained in the third section. The control laws are presented in section IV and V. The experimental set up is described in section VI with the stabilisation results. The energy production results are shown in section Fig. 2. Oscilations observed in the x-axis of the wing i.e. the angle θ, when no spoilers are used.

VII . Conclusions and future work are finally given in section VIII.

AERODYNAMICAL MODEL

This section presents the aerodynamical model of a kite wind power generator system as depicted in figure 3.

Fig. 3. Scheme of the wing in the airflow From Newton's second law we obtain the following nonlinear dynamical system:

r = r θ2 + 1 M [F r -T ] (1) θ = 1 r [- 1 M F t -2 ṙ θ] ( 2
)
where r is the rope length from the kite to the dynamo on the ground, θ is the angle the rope makes with respect to the horizontal line, M is the mass of the kite, T is the tension on the rope, F r and F t are respectively the radial and tangential forces acting on the kite due to aerodynamical forces and the weight w = M g. The aerodynamical forces can be expressed in terms of the lift L and the drag D as follows:

F r = L sin(θ -α w ) + D cos(θ -α w ) -w sin θ (3)
F t = L cos(θ -α w ) -D sin(θ -α w ) -w cos θ (4)
where α w is the relative wind angle [START_REF] Houska | Optimal Control of Towing Kites. Conference on Control and Decision[END_REF] due to the movement of the kite and

L = 1 2 ρSv 2 r C L (5) D = 1 2 ρSv 2 r C D ( 6
)
where ρ is the air density, S is the wing surface, v r is the relative wind speed and the lift coefficient C L and the drag coefficient C D are

C L = ∂C L ∂α (α w + α u ) + C L0 (7)
C D = C 2 L πeλ + C D0 + β θ u θ ( 8
)
where C L0 , C D0 , β θ are positive coefficients, e is the Oswald's efficiency factor, λ is the aspect ratio, α u is the control input of the angle of attack of the wing α. Note that α = α w + α u .The input u θ controls the spoilers which modify the drag of the wing, see figure 7 . Therefore the systems has three control inputs, α u , u θ and T which is the tension on the rope which is related to the torque on the dynamo/motor. Notice that the dynamo/motor will act as a motor or as a dynamo depending on whether it is generating power or pulling the kite down to restart a duty cycle. For simplicity we are neglecting the actuator dynamics.

The angle of attack α is equal to the sum of the pitch angle α u and the angle of apparent wind α w due to the movement of the kite.

The horizontal relative air speed (9) and the vertical airspeed (10) depend of the movements of the kite. When the kite is static, θ and ṙ are equal to 0 and the air speed relative to the kite is equal to the airspeed speed relative to the ground defined as V.

v h = V + r sin(θ) θ -ṙ cos(θ) (9)
v v = r cos(θ) θ + ṙ sin(θ) 2 (10)
Using ( 9) and ( 10)we obtain the wind angle α w which is the angle of the wind velocity vector mesured with respect to the horizon, and the kite's relative wind velocity v r which is to be used in ( 5) and [START_REF] Roberts | Unmanned Rotorcraft to Generate Electricity Using Upper Atmospheric Winds[END_REF].

α w = -arctan r cos(θ) θ + ṙ sin(θ) V + r sin(θ) θ -ṙ cos(θ) (11) v r = (r cos(θ) θ + ṙ sin(θ)) 2 + (V + r sin(θ) θ -ṙ cos(θ)) 2 (12)

FLIGHT PLAN AND ENERGY PRODUCTION

The amount of en energy generated during a cycle depends of the flight plan. In order to maximise it, we want to maximise F r during the production phase , and to reduce it to the lowest possible value during the recovery phase.

In order to optimise the power production, we work on the two mains inputs that are the angle of attack and the ṙ speed. We can define a basic flight plan with 6 parameters that are : ṙ1 the speed and F r1 the traction of the kite during the production phase, ṙ2 the speed and F r2 the traction of the kite during the recovery phase. r min is the length of the rope at the initial time t 0 and at the end of the cycle time t 2 . Finally r max is the maximum rope length attained at the time t 1 . The flight plan is showed in the figure 4.

The amount of energy produced during a cycle of this flight plan is given by:

E(t 0 , t 1 , t 2 ) = t1 t0 F r1 ṙ1 dt + t2 t1 F r2 ṙ2 dt (13)
the first term is the produced energy when the kite goes up, the second one is the consumed energy for used for pulling back the kite.

Therefore, we can write the mean power as:

P (t 0 , t 1 , t 2 ) = E(t 0 , t 1 , t 2 ) t 2 -t 0 (14)
We can see that the optimisation of the mean power depends of the ṙ terms,and the values of t 0 , t 1 and t 1 ,t 2 . We define r min as the length of the rope at the initial time t 1 and the end of the cycle time, t 2 . r max is the maximum rope length attained at the time t 1 . We can define a cyclic flight plan with the 6 following values: ṙ1 , ṙ2 , r max , r min , α 1 , α 2 .

variable production phase recovery phase ṙ 0.3 m/s 0.5 m/s α 19 degrees 5 degrees r 1 m (max value) 0.2m (min value)

This simple flight plan and the energy formula don't take into account the loss of energy due to a lot of other things like the the drag of the spoilers, the efficiency of the motors and actuators,the weight of the rope,the safety of the flight plan, the variation of the wind speed, friction loss , etc.. All these parameters are quite important and need to be studied in future work.

CONTROL OF THE ANGULAR POSITION θ

We first study the control strategy for the angular oscillations on θ which will be stabilized using the spoilers u θ .

We have noticed that without spoilers it is very difficult to stabilize the angular position θ because there exist sever oscillations, see figure 2. Such oscillations have been damped by introducing spoilers on the wing. To analyze the system we have to rewrite (2) in such a way that the spoiler control input u θ appears in the equation. Injecting (4) ( 5), ( 6) and ( 8) into (2) leads to θ = χ + Qu θ [START_REF] Lansdorp | Development and testing of a steering mechanism for Laddermill kites[END_REF] where

Q = ρSv 2 r β θ sin(α w -θ) 2M r χ = - 2 ṙ θ r - L cos(α w -θ) -w cos θ rM + ρSv 2 r ( C 2 L πeλ + C D0 ) sin(α w -θ) 2rM
The objective is to propose a control law for u θ such that the closed loop system behaves as

θ + 2 θ + θ = θ d (16)
where θ d is the desired value for θ which is typically 60 • . Using Laplace transform the above equation can be expressed as

(s + 1) 2 θ(s) = θ d
This means that the desired closed loop system has two real poles at -1. Introducing (15) into ( 16) it follows that the required control input is

u θ = (-2 θ -θ + θ d -χ)/Q (17)
The above control input linearizes the system (15) using a state feedback technique and guarantees that θ → θ d . The desired angular position θ d can be also computed to maximize the generated energy.

CONTROL OF THE ROPE LENGTH R

In this section we propose a control strategy such that the rope length r tracks a sinusoidal trajectory. We have chosen a sinusoidal trajectory in order to smooth the displacement of the kite. Notice that the kite velocity is zero, i.e. ṙ = 0, when the kite reaches its highest and lowest altitudes. Suppose for simplicity that we wish r to follow the desired trajectory: From ( 1) and the above it follows that the control input T should be chosen as

T = M -r d + 2( ṙ -ṙd ) + (r -r d ) -r θ2 + F r
The above control law is a state feedback control law which linearizes system (1) and forces it to track a sinusoidal trajectory. The amplitude of such a control input normally increases as A or f grow.

EXPERIMENTAL RESULTS

In order to validate our control laws, an experimental setup has been built in the Gipsa-lab laboratory for testing the proposed wind power generation system presented in the previous sections. A wind tunnel has been built so that an air flow is produced at speeds up to 12 m/s, see figure 1. The wind tunnel has been built to be able to test our prototype indoors so as to be independent of the weather conditions. The wind tunnel is composed of 9 brushless motors of 800W each distributed on a surface of 1 m 2 . A wing has been built using foam and carbon fiber tubes to increase the rigidity. We have chosen a wing profil NACA 0018. The yaw angular position ψ of the kite is controlled using a rudder (figure 6), which is controlled by PID gyroscopic control system FUTABA GY 5020 that stabilises the orientation ψ and the angular velocity ψ .

The wing has ailerons driven by servomechanisms which can be deployed to form spoilers. The spoilers are used to vary the drag D to introduce damping in the wing displacement.

Indeed, as it can be seen in figure 2, without spoilers there exist oscillations that might destroy the prototype. For simplicity purposes, the angle of attack α , is controlled using three strings attached to servomechanisms on the ground. Two strings are attached to the back tips of the wing and one is attached to the front of the wing. In future work we will consider controlling the angle of attack with ailerons which is a more realistic way to stabilize the angle of attack of the wing.

The coefficients of the experimental set-up are given in Table 1 Fig. We have observed experimentally that the control of the yaw angular displacement ψ works properly. Furthermore, for a fixed rope length r, we have observed that for a given angle of attack α,we can see on figure 8 that the rope angle θ reaches an almost constant value . The small oscillations are due to the fact that the wind in these experimental conditions is very turbulent.

The control of the rope length r works quite well, as il can be seen in figure 9. We have noticed that in the case of very low speed displacements, the behavior of the system has some problems due to the system's static friction.

RESULTS OF OUR SYSTEM AND FUTURE IMPROVEMENTS

We have smoothen the flight plan to avoid brutal variations by using a first order filter. We can see on figure 10 the mesurement of tracking of the desired trajectory during 4 cycles.

We can see that the amount of energy consumed during the recovery phase is about 50% of the energy generated during the production phase. This is due to many factors, the two most importants are that the angle of attack can't be lower than 5 degrees, the other factor is due to the system's friction. during the generation phase we can mesure a mean power of approximatively 1 W and the consumed power during the recovery phase is around 0.2 using the energy equation ( 13) and power formula [START_REF] Lansdorp | Design and Testing of a Remotely Controlled Surfkite for the Laddermill[END_REF], we can estimate that the mean power production is around 0.1 W. During our experimental work, we have tested various speed combinations between 0.1 and 1 m/s. We have noticed that for very hight ṙ speeds around 1 m/s, the spoilers are no longer able to stabilise the θ angle. We have also some structural problems, strong accelerations tend to bend the wing and the model is no longer valid. for these reasons , we do our tests with maximum speeds around 0.5 m/sec.

The angle of attack is controlled by the ropes lengths, it's accuracy depends of the tension of the rope. Therefore, we need to keep the tension greater than a minimum value which corresponds to the lift generated by 5 degrees of angle of attack.

This kite has proved that a producing energy during a controlled flight is possible. Building models with a higher aspect ratio, lower weight and a more adapted airfoil should enhance a lot the energy production, these are our next steps...

Conclusions

In this paper we have studied a wind power generation system mainly composed of a wing attached to a dynamo with a rope. The wing subject to airflow generates a lift force which drives the dynamo on the ground to generate power. We have focused our interest in such a configuration because it is less expensive than current wind power systems and could be used in rural areas without infrastructure. We have presented an aero dynamical non linear model of the wind power generation system. We have introduced spoilers in the wing so that damping introduced into the system cancels oscillations in the displacement of the wing. We have also proposed a control strategy to stabilize the wing displacement and the rope length. Such control laws are based on state feedback linearization. An experimental set-up has been built including a wind tunnel, a wing with spoilers and a system for controlling the wing angle of attack. The experimental prototype has been used to the control of the yaw angle, the angle of attack and the angle of the rope.

Current research is underway to study the control of the rope length r by using the angle of attack as an input control.

Fig. 4 .

 4 Fig. 4. flight plan .

Fig. 5 .

 5 Fig. 5. The Gipsa wind tunnel system, including the 7.2 kW wind tunnel, the wing, the ground station and the control computer . r d = A sin ωt ṙd = Aω cos ωt rd = -Aω 2 sin ωt where A and f = ω/2π are respectively the amplitude and frequency of the desired oscillation. Let us introduce the error e = r -r dand let us assume that we wish the closed loop system to behave as

Fig. 6 .Fig. 7 .

 67 Fig. 6. The wing with the spoilers in the minimum drag position. Notice the rudder used for the ψ angle stabilization.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Desired position tracking with Torque control .

  8. Response to a manual perturbation on the closed loop system.

	symbol	name	value
	M	mass	.1 Kg
	ρ	air density	1.225 Kg/m 3
	S	wing area	0.15 m 2
	e	Oswald's factor	0.7
	λ	aspect ratio	2
	∂Cl ∂α β θ	lift derivative w.r.t. α drag derivative w.r.t. u θ	0.05 deg -1 0.4
	C D0	zero lift drag	0.1
	V	mean air speed	8 m/s