Mircea Parpalea

Nicoleta Avesalon

Eleonor Ciurea

Minimum parametric flow over time

Keywords: dynamic networks, parametric flow, partitioning algorithm

The paper presents a dynamic solution method for dynamic minimum parametric networks flow. The solution method solves the problem for a special parametric dynamic network with linear lower bound functions of a single parameter. Instead directly work on the original network, the method implements a labelling algorithm in the parametric dynamic residual network and uses quickest paths from the source node to the sink node in the time-space network along which repeatedly decreases the dynamic flow for a sequence of parameter values, in their increasing order. In each iteration, the algorithm computes both the minimum flow for a certain subinterval of the parameter values, and the new breakpoint for the maximum parametric dynamic flow value function.

Introduction

Dynamic flow problems where networks structure changes depending on a scalar parameter λ are widely used (see for example [START_REF] Aronson | A survey of dynamic network flows[END_REF]) to model different network-structured, decision-making problems over time. These types of problems are arising in various real applications such as communication networks, air/road traffic control, and production systems. Moreover, in many applications of graph algorithms, including communication networks, graphics, assembly planning, and scheduling, graphs are subject to discrete changes, such as additions or deletions of arcs or nodes. In the last decade there has been a growing interest for such dynamically changing graphs, and a whole body of algorithms and data structures for dynamic graphs has been discovered: [START_REF] Fathabadi | A simple algorithm for reliability evaluation in dynamic networks with stochastic transit times[END_REF], [START_REF] He | Model and a solution algorithm for the dynamic resource allocation problem for large-scale transportation network evacuation[END_REF], [START_REF] Nassir | Efficient negative cycle-canceling algorithm for finding the optimal traffic routing for network evacuation with nonuniform threats[END_REF], [START_REF] Rashidi | Vehicle Scheduling in Port Automation: Advanced Algorithms for Minimum Cost Flow Problems[END_REF] or [START_REF] Zheng | On the system optimum dynamic traffic assignment and earliest arrival flow problems[END_REF]. Further on, the next section presents some basic discrete-time dynamic networks terminology and notations and Section 3 introduces the parametric minimum flow over time problem. In Section 4 the algorithm for solving the parametric minimum flow in discrete dynamic networks is presented and Section 5 ilustrates how the algorithmon work on a given dynamic network. An arc a ∈ A from node i to node j is denoted by (i, j). Parallel, as well as opposite arcs are not allowed in graph G. The following time-dependent functions are associated with each arc a = (i, j) ∈ A: the upper bound (capacity) function u(i, j; θ), u : A × H → + , representing the maximum amount of flow that can enter the arc (i, j) at time θ, the lower bound function (i, j; θ), : A × H → + , i.e. the minimum amount of flow that must enter the arc (i, j) at time θ, and the transit time function h(i, j; θ), h : A × H → ℵ. Time is measured in discrete steps, so that if one unit of flow leaves node i at time θ over the arc a = (i, j) that one unit of flow arrives at node j at time θ + h(i, j; θ), where h(i, j; θ) is the transit time of the arc (i, j). The time horizon T represents the time limit until which the flow can travel in the network. The network has two special nodes: a source node s and a sink node t.

Time-space network

For a given discrete-time dynamic network, the time-space network is a static network constructed by expanding the original network in the time dimension, by considering a separate copy of every node i ∈ N at every discrete time step θ ∈ H. A node-time pair (NTP) (i, θ) refers to a particular node i ∈ N at a particular time step θ ∈ H, i.e., (i, θ) ∈ N × H.

Definition 1

The time-space network G T of the original dynamic network G is defined as follows:

(a) N T := {(i, θ)|i ∈ N, θ ∈ H}; (b) A T := {a θ = ((i, θ), (j, θ + h(i, j; θ)))|0 ≤ θ ≤ T -h(i, j; θ), (i, j) ∈ A}; (c) u T (a θ ) := u(a; θ) f or a θ ∈ A T ; (d) T (a θ ) := (a; θ) f or a θ ∈ A T .
For every arc (i, j) ∈ A with traversal time h(i, j; θ), capacity u(i, j; θ) and lower bound (i, j; θ), the time-space network G T contains the arcs ((i, θ), (j, θ + h(i, j; θ))), θ = 0, 1, • • • , T -h(i, j; θ) with capacities u(i, j; θ) and lower bounds (i, j; θ). A (discrete-time) dynamic path P is defined as a sequence of distinct, consecutively linked NTPs.

Parametric dynamic network

Definition 2 A discrete-time dynamic network G = (N, A, T ) for which the lower bounds (i, j; θ) of some arcs (i, j) ∈ A are functions of a real parameter λ is referred to as a parametric dynamic network and is denoted by Ḡ = (N, A, u, ¯ , h, T ).

For a parametric dynamic network Ḡ, the parametric lower bound function ¯ : A × H × [0, Λ] → + associates to each arc (i, j) ∈ A and for each of the parameter values λ in an interval [0, Λ] the real number ¯ (i, j; θ; λ), referred to as the parametric lower bound of arc (i, j):

¯ (i, j; θ; λ) = 0 (i, j; θ) + λ • L(i, j; θ), λ ∈ [0, Λ], θ ∈ H, (1)
where L : A × H → is a real valued function, associating to each arc (i, j) ∈ A and for each time-value θ ∈ H the real number L(i, j; θ), referred to as the parametric part of the lower bound of the arc (i, j). The nonnegative value 0 (i, j; θ) is the lower bound of the arc (i, j) for λ = 0, i.e., ¯ (i, j; θ; 0) = 0 (i, j; θ). For the problem to be correctly formulated, the lower bound function of every arc (i, j) ∈ A must respect the condition u(i, j; θ) ≥ ¯ (i, j; θ; λ) for the entire interval of the parameter values, i.e., ∀(i, j) ∈ A, ∀θ ∈ H and ∀λ ∈ [0, Λ]. It follows that u(i, j; θ) ≥ 0 (i, j; θ) and the parametric part of the lower bounds L(i, j; θ) must satisfy the constraints: L(i, j; θ) ≤ [u(i, j; θ) -0 (i, j; θ)]/Λ, ∀(i, j) ∈ A.

Parametric flow over time

The parametric dynamic flow value function v : N × H × [0, Λ] → associates to each of the nodes i ∈ N , at each time moment θ ∈ H, a real number v(i; θ; λ) referred to as the value of node i at time θ, for each of the parameter λ values.

Definition 3 A feasible parametric flow over time f (i, j; θ; λ) in a parametric dynamic network Ḡ = (N, A, u, ¯ , h, T ) is a function f : A × H × [0, Λ] → + that satisfies the following constraints for all λ ∈ [0, Λ]:

j|(i,j)∈A f (i, j; θ; λ) - j|(j,i)∈A ϑ|ϑ+h(j,i;ϑ)=θ f (j, i; ϑ; λ) = { v(i; θ; λ), i = s, t 0, i = s, t , ∀i ∈ N, ∀θ ∈ H; (2)
¯ (i, j; θ; λ) ≤ f (i, j; θ; λ) ≤ u(i, j; θ), ∀(i, j) ∈ A, ∀θ ∈ H; (3)
f (i, j; θ; λ) = 0, ∀(i, j) ∈ A, ∀θ ∈ [T -h(i, j; θ) + 1, T ]; (4) θ∈H v(s; θ; λ) = - θ∈H v(t; θ; λ) = v(λ); (5)
where f (i, j; θ; λ) determines the rate of flow (per time unit) entering arc (i, j) at time θ, for the parameter value λ, ∀θ ∈ {0, 1,

• • • , T } and ∀λ ∈ [0, Λ].
Definition 4 The parametric minimum flow over time (PmFT) problem is to compute all minimum flows over time for every possible value of λ:

minimise v(λ) = θ∈{0,1,•••,T } v(s; θ; λ), f or all λ ∈ [0, Λ], (6)
under flow constraints (2)-( 5).

Definition 5 For the minimum flow over time problem, the parametric time-dependent residual network with respect to a given feasible parametric flow over time f is defined as

Ḡ( f ) := (N, A( f ), T ), with A( f ) := A + ( f ) ∪ A -( f ), where A + ( f ) := {(i, j)|(i, j) ∈ A, ∃θ ≤ T -h(i, j; θ) with f (i, j; θ; λ) -¯ (i, j; θ; λ) > 0}; (7) A -( f ) := {(i, j)|(j, i) ∈ A, ∃θ ≤ T -h(j, i; θ) with u(j, i; θ) -f (j, i; θ; λ) > 0}. (8)
The direct arcs (i, j) ∈ A + ( f ) in Ḡ( f ) have same transit times h(i, j; θ) with those in Ḡ while the reverse arcs (i, j) ∈ A -( f ) have negative transit times h(i, j; θ + h(j, i; θ)) = -h(j, i; θ), with (j, i) ∈ A and 0 ≤ θ + h(j, i; θ) ≤ T .

Definition 6 For the parametric minimum flow over time (PmFT) problem, the parametric residual capacities of arcs (i, j) in the parametric time-dependent residual network Ḡ( f ) are defined as follows:

r(i, j; θ; λ) = f (i, j; θ; λ) -¯ (i, j; θ; λ), (i, j) ∈ A, 0 ≤ θ + h(i, j; θ) ≤ T ; (9) r(i, j; θ + h(j, i; θ); λ) = u(j, i; θ) -f (j, i; θ; λ), (j, i) ∈ A, 0 ≤ θ + h(j, i; θ) ≤ T. ( 10
)
Definition 7 Given a parametric flow over time f (i, j; θ; λ), the parametric residual capacity r( P ; θ; λ) of a dynamic path P is the inner envelope of the parametric residual capacity functions r(i, j; θ; λ) for all arcs (i, j) composing the path and for all parameter λ values:

r( P ; θ; λ) = min{r(i, j; θ; λ)|(i, j) ∈ P }. ( 11
)
Definition 8 The transit time τ ( P ; θ) of a dynamic path P is defined by:

τ ( P ) = (i,j)∈ P h(i, j; θ). (12)
A dynamic path P is referred to as the quickest dynamic path if τ ( P ) ≤ τ ( P ) for all dynamic paths P in the parametric time-dependent residual network Ḡ( f ).

Parametric minimum dynamic flow

The approaches for solving the maximum parametric flow over time problem via applying classical algorithms can be grouped in two main categories: by applying a classical parametric flow algorithm (see [START_REF] Hamacher | Algorithms for flows with parametric capacities[END_REF] and [START_REF] Ruhe | Algorithmic Aspects of Flows in Networks[END_REF]) in the static time-space network (for this approach see the algorithm presented in [START_REF] Parpalea | Partitioning algorithm for the parametric maximum flow[END_REF]) or by applying a non-parametric maximum dynamic flow algorithm (see the algorithms presented in [START_REF] Parpalea | The quickest maximum dynamic flow of minimum cost[END_REF]) in dynamic residual networks generated by partitioning the interval of the parameter values (partitioning method presented by Avesalon et al.). Given the powerful versatility of dynamic algorithms, it is not surprising that these algorithms and dynamic data structures are often more difficult to design and analyse than their static counterparts (see [START_REF] Orlin | Max flows in o(nm) time, or better[END_REF]). [START_REF] Cai | Time-varying Network Optimization[END_REF] proved that the complexity of finding a shortest dynamic flow augmenting path, by exploring the forward and reverse arcs successively, is O(nmT 2 ). For algorithms which explores the two sub-networks simultaneously, Miller-Hooks and Patterson ( 2004) also reported a complexity of O(n 2 T 2 ). By using special node addition and selection procedures, [START_REF] Nasrabadi | Minimum cost time-varying network flow problems[END_REF] succeeded to reduce significantly the number of node time pair that needs to be visited. The worst-case complexity of their algorithm is O(nT (n + T )). As far as we know, the problem of minimum flow over time in parametric networks has not been treated yet.

Parametric minimum dynamic flow (PmDF) algorithm

The idea of the algorithm is that if the parametric residual capacities for all arcs in Ḡ( f ) are maintained linear functions of λ, with no break points, the problem can be solved via a slightly modified non-parametric algorithm. Firstly, if one exists, a feasible flow must be established. The most convenient is this to be done in the nonparametric network G = (N, A, u, , h, T ) obtained from the initial network by replacing the parametric lower bound functions with the non-parametric ones: (i, j; θ) = max{ (i, j; ; θ; λ)|λ ∈ [0, Λ]}, i.e. (i, j; θ) = 0 (i, j; θ) for L(i, j; θ) ≤ 0 and (i, j; θ) = 0 (i, j; θ) + Λ • L(i, j; θ) for L(i, j; θ) > 0. For finding a feasible flow f (i, j; θ) in Ḡ see the algorithms presented in [START_REF] Ahuja | Network Flows. Theory, algorithms and applications[END_REF], applied in the static time-space network. During the running of the algorithm, the subinterval of the parameter values is continuously narrowed so that the above restriction to remain valid. In each subinterval [λ k , λ k+1 ] of the parameter values, the linear parametric residual capacity of every arc (i, j) can be written as rk (i; j; θ

; λ) = α k (i, j; θ) + β k (i, j; θ) • (λ -λ k ) while the parametric flow is fk (i; j; θ; λ) = f k (i, j; θ) + F k (i, j; θ) • (λ -λ k ).
As soon as the parametric time-dependent residual network Ḡ( f ) contains no dynamic paths, the algorithm computes the minimum flow for the considered subinterval and then reiterates on the next subinterval of the parameter values, until Λ value is reached.

Quickest Dynamic Paths (QDP) procedure

The successive shortest (quickest) dynamic paths procedure repeatedly performs the following operations: -Finds a shortest (quickest) dynamic path P in the time-dependent residual network Ḡ( f ) based on the successor vector σ : N × H → (i, θ) ∈ N × H; -Computes the parametric residual capacity r( P ; θ; λ) of the dynamic path P ; -Computes the first (in increasing order) value of the parameter λ up to which the parametric residual capacity of the dynamic path remains linear without break points; -Updates the subinterval of the parameter values for which the minimum flow is determined; -Decreases the flow along the dynamic path and updates the time-dependent residual network. The algorithm ends when none of the source nodes, at any time moments θ ∈ {0, 1, • • • , T }, is reachable from any of the sink nodes, i.e. there is no dynamic path from s to t.

(01) PROCEDURE QDP(k, λ k , BP ); (02) BEGIN (03) FOR all θ ∈ H DO (04) BEGIN (05) FOR all i ∈ N DO σ(i, θ) := (0, 0); (06) FOR all (i, j) ∈ A DO α k (i, j; θ) := f (i, j; θ) -0 (i, j; θ) -λ k • L(i, j; θ); β k (i, j; θ) := -L(i, j; θ); α k (j, i; θ + h(i, j; θ)) := u(i, j; θ) -f (i, j; θ); β k (j, i; θ + h(i, j; θ)) := 0; f k (i, j; θ) := f (i, j; θ); F k (i, j; θ) := 0; (07) END; (08) C := 1; λ k+1 = Λ; (09) LS(σ, C);
(10) WHILE(C = 1) DO ( 11) BEGIN ( 12) build P based on σ starting from (s, θ); (13) α := min{α k (i, j; θ)|(i, j) ∈ P }; β := min{β k (i, j; θ)|(i, j) ∈ P , α k (i, j; θ) = α}; (14) (i, θ) := (s, θ); (15) WHILE(i = t) DO ( 16) BEGIN (17) (j, ϑ) := σ(i, θ); Algorithm 2: Quickest Dynamic Paths (QDP) procedure.

(18) IF(β k (i, j; θ) < β) THEN λ k+1 := min{λ k+1 , λ k + (α k (i, j; θ) -α)/(β -β k (i, j; θ))}; (19) α k (i, j; θ) := α k (i, j; θ) -α; β k (i, j; θ) := β k (i, j; θ) -β; α k (j, i; ϑ) := α k (j, i; ϑ) + α; β k (j, i; ϑ) := β k (j, i; ϑ) + β; (20) IF((i, j) ∈ A + ( f )) THEN f k (i, j; θ) := f k (i, j; θ) -α; F k (i, j; θ) := F k (i, j; θ) -β; (21) ELSE f k (j, i; ϑ) := f k (j, i; ϑ) + α; F k (j, i; ϑ) := F k (j, i; ϑ) + β; (22) (i, θ) := (j,

Labels Setting (LS) procedure

The label setting procedure uses transit time labels τ (i, θ) associated to all nodes at each discrete time values. At any step of the algorithm, a label is permanent once it denotes the length of shortest augmenting path to a node-time pair, otherwise it is temporary. The LS procedure maintains a set L of candidate nodes in increasing order of their temporary labels, which initially includes only the sink nodes (t, θ), θ ∈ {0, 1, • • • , T }. For every node-time pair (j, ϑ) selected from the list, the arcs with positive residual capacity connecting (i, θ) to (j, ϑ) are explored, where ϑ = θ + h(i, j; θ) if the arc connecting (i, θ) to (j, ϑ) is a forward arc and 0 ≤ θ = ϑ + h(j, i; ϑ) if (i, j) is a reverse arc. At any iteration, the algorithm selects the node-time pair (i, θ) with the minimum temporary label, makes its transit time label permanent, checks optimality conditions and updates the labels accordingly. A transit time label τ (i, θ) represents the length (transit time) of a shortest (quickest) dynamic path from (i, θ) if τ (i, θ) ≤ τ (j, ϑ) + h(j, i; ϑ), ∀(i, j) ∈ A( f ). The process is repeated until there are no more candidate nodes in L. The transit time of the shortest (quickest) path P computed based on successor vector σ is given by τ .

(01) PROCEDURE LS(σ, C);

(02) BEGIN (03) FOR all θ ∈ {0, 1, • • • , T } DO (04) BEGIN (05) FOR all i ∈ N -t DO τ (i, θ) := ∞; (06) τ (t, θ) := 0; L := L ∪ {(t, θ)}; (07) END; (08) τ := ∞; θ := ∞ (09) WHILE(L = ∅) DO (10) BEGIN (11) select the first (j, ϑ) from L; L := L -{(j, ϑ)}; (12)
FOR all i ∈ N with (i, j) ∈ A + ( f ) DO ( 13) FOR all θ with θ + h(i, j; θ) = ϑ DO ( 14) IF(τ (j, ϑ) + h(i, j; θ) < τ (i, θ)) THEN ( 15) BEGIN ( 16)

τ (i, θ) := τ (j, ϑ) + h(i, j; θ) (17) σ(i, θ) := (j, ϑ) (18) IF((i, θ) / ∈ L) THEN L := L ∪ {(i, θ)}; (19) END; (20) FOR all i ∈ N with (i, j) ∈ A -( f ) DO (21)
IF(ϑ + h(j, i; ϑ) ≤ T and τ (j, ϑ) -h(j, i; ϑ) < τ (i, ϑ + h(j, i; ϑ))) THEN ( 22) BEGIN ( 23)

τ (i, ϑ + h(j, i; ϑ)) := τ (j, ϑ) -h(j, i; ϑ); (24) σ(i, ϑ + h(j, i; ϑ)) := (j, ϑ); (25) IF((i, ϑ + h(j, i; ϑ)) / ∈ L) THEN L := L ∪ {(i, ϑ + h(j, i; ϑ))}; (26) END; (27) END; (28) τ := min θ∈{0,1,•••,T } {τ (s, θ)}; θ := min θ∈{0,1,•••,T } {θ|τ (s, θ) = τ }; (29) IF(τ = ∞) THEN C := 0; (30) END.
Algorithm 3: Labels Setting (LS) procedure.

Theorem 1 [Correctness of PmDF algorithm] PmDF algorithm computes correctly a parametric minimum flow over time for a given time horizon T and for λ ∈ [0, Λ].

Proof:

The partitioning type algorithm iterates on successive subintervals [λ k , λ k+1 ], starting with λ 0 = 0 and ending with λ kmax+1 = Λ and consequently, the correctness of the algorithm obviously follows from the correctness of the Quickest Dynamic Paths (QDP) procedure which ends when none of the source node-time pairs is reachable from any of the sink node-time pairs, i.e. when there is no dynamic path from the source node to the sink node in the time-depending residual network. According to the classical flow decreasing path theorem presented by [START_REF] Ahuja | Network Flows. Theory, algorithms and applications[END_REF], this means that the obtained flow is a minimum dynamic flow for the given time horizon. In fact, the algorithm ends with a set of linear parametric minimum flows and with the partition BP of the interval of the parameter values in their corresponding subintervals.

2

Theorem 2 [Time complexity of PmDF algorithm] The parametric minimum dynamic flow (PmDF) algorithm runs in O(Kn 2 mT 3 ) time, where K + 1 is the number of λ values in the set BP at the end of the algorithm.

Proof: The building, as well as the updating of the time-dependent residual network requires an O(mT ) running time, since for each of the time values 0 ≤ θ ≤ T all the m arcs must be examined. Labels Setting (LS) procedure investigates at most nT adjacent node-time pairs for each of the node-time pairs which are removed from the list L (i.e. O(nT ) times). Thus, the complexity of labels setting (LS) procedure is O(n 2 T 2 ). Considering that in each of the iterations of the QDP procedure, for one of the time values, one arc is eliminated from the dynamic residual network, the algorithm end in at most O(mT ) iterations.

On each of the iterations the procedure finds a quickest dynamic path with the complexity O(n 2 T 2 ) and updates the time-dependent residual network in O(mT ) time. Thus, the total complexity of the quickest dynamic paths (QDP) procedure is O(n 2 mT 3 ).

For each of the K subintervals [λ k , λ k+1 ], k = 0, 1, • • • , K -1 in which the interval [0, Λ] of the parameter λ values is partitioned in, the algorithm makes a call to procedure QDP. Consequently, the complexity of the parametric minimum dynamic flow (PmDF) algorithm is O(Kn 2 mT 3 ). 2

Example

In the discrete-time dynamic network presented in Figure 1(a), node 1 is the source node s and node 4 is the sink node t; the time horizon being set to T = 3, i.e. H = {0, 1, 2, 3}. For the interval of the parameter λ values, set to [0, 1], i.e. Λ = 1, the transit times h(i, j; θ), parametric lower bound functions ¯ (i, j; θ; λ) = 0 (i, j; θ) + λ • L(i, j; θ), upper bounds (capacities) u(i, j; θ) and feasible flow f (i, j; θ) for all arcs in Ḡ are indicated in Table 1.

Tab. 1: Characteristics of dynamic network Ḡ presented in Figure 1(a) (i, j) h(i, j; θ) 0 (i, j; θ) L(i, j; θ) u(i, j; θ) f (i, j; θ)

(1, 2) 1, θ = 0 2, θ ≥ 1 3, θ = 0 0, θ ≥ 1 -2, θ = 0 0, θ ≥ 1 5 5, θ = 0 0, θ ≥ 1 (1, 3) 1, 0 ≤ θ < 2 2, θ ≥ 2 1, 0 ≤ θ < 2 0, θ ≥ 2 4, θ = 0 1, θ = 1 0, θ > 1 5 5, θ = 0 2, θ = 1 0, θ ≥ 1 (2, 3) 1, θ ≥ 0 0, θ ≥ 0 3, θ = 1 0, θ = 1 5 3, θ = 1 0, θ = 1 (2, 4) 1, 0 ≤ θ < 2 2, θ ≥ 2 0, θ ≥ 0 0, θ ≥ 0 5 2, θ = 1 0, θ = 1 (3, 4) 2, 0 ≤ θ < 2 1, θ ≥ 2 0, θ = 0 2, 0 < θ ≤ 2 0, θ > 2 -2, θ = 2 0, θ = 2 5 0, θ = 0 5, 0 < θ ≤ 2 0, θ > 2
After the initialisation step, for k = 0 and for the corresponding initial value of the parameter λ0 = 0, procedure QDP is called for the first time. The successor vector is initialised for all nodes at all time values to σ(i, θ) := (0, 0) and the Labels setting (LS) procedure is called for finding a quickest dynamic path in the timedependent residual network Ḡ( f ). The distance labels are initialised to τ (4, θ) := 0, ∀θ ∈ {0, 1, 2, 3} and the set of candidate nodes is set to L := {(4, 0), (4, 1), (4, 2), (4, 3)}. After setting transit time labels and finding successor values for all node-time pairs, the procedure computes the minimum label of the source node, τ := min{τ (1, 0), τ (1, 1), τ (1, 2), τ (1, 3)} = min{2, 1, ∞, ∞} = 1 with θ := 1. Since τ (4) = ∞ the procedure LS ends with the variable C keeping its initial value C = 1. Based on successor vector, the quickest dynamic path P := ((1, 1), (3, 2), (2, 1), (4, 2)) is built and its residual capacity r( P ; θ; λ) = α + (λ -λ k ) • β is computed with α = min{α0(1, 3; 1), α0(3, 2; 2), α0(2, 4; 1)} = min{1, 2, 2} = 1 and β = -1. After the time-dependent residual network is updated and the parametric dynamic flow is decreased along the shortest dynamic path, the successor vector is reinitialised and procedure LS is called again. The next shortest dynamic path found by QDP procedure is P := ((1, 0), (2, 1), (4, 2)) with r( P ; θ; λ) = 1+λ and both the time-dependent residual network and the parametric dynamic flow are accordingly updated. Then QDP procedure finds the new dynamic path P := ((1, 0), (3, 1), (4, 3)) with α = 3 and β = 0. Since β0(1, 3; 0) = -4 < β = 0, the upper limit λ k+1 of the subinterval of the parameter is updated to λ1 = min{1, (4 -3)/(0 + 4)} = 1/4. Finaly, after updating the time-dependent residual network and the parametric dynamic flow, based on successor vector found by LS procedure, QDP finds the last path P := ((1, 0), (2, 1), (3, 2), (4, 3)) with α = 1 and β = 1. The validity of the upper limit of the subinterval of the parameter values is tested and it is maintained unchanged since for the arc (2, 3) at θ = 1, β0(2, 3; 1) = -4 < β = 1 but λ1 = min{1/4, (4 -1)/(1 + 4)} = min{1/4, 3/5} = 1/4. At this point, after the updating step, no other dynamic path can be found in the time-dependent residual network so that C := 0 is set, the breakpoints list is updated to BP = {0, 1/4} and the first iteration ends with the minimum parametric flow over time f 0 (i, j; θ; λ) computed for the subinterval [0, 1/4] of the parameter λ values. The PmDF algorithm increments variable k to the value k = 1 and a next iteration is performed. The evolution of the algorithm is presented in Table 2. The piecewise linear minimum flow over time value function for the discrete dynamic network Ḡ, computed by Parametric minimum Dynamic Flow (PmDF) algorithm, is presented in Figure 1(b).

Tab. 2: The evolution of PmDF algorithm for the dynamic network Ḡ k λ k P rk ( P ; θ; λ) λ k+1 (1, 1), (3, 2), (2, 1), (4, 2)

1 -λ 1 Iteration 1 : 0 0 (1, 0), (2, 1), (4, 2) 1 + λ 1 (1, 0), (3, 1), (4, 3) 3 1/4 (1, 0), (2, 1), (3, 2), (4, 3) 1 + λ 1/4 (1, 1), (3, 2), (2, 1), (4, 2) 1 -λ 1 Iteration 2 : 1 1/4 (1, 0), (2, 1), (4, 2) 1 + λ 1 (1, 0), (3, 1), (4, 3) 4 -4λ 1 (1, 0), (2, 1), (3, 2), (4, 3) 1 + λ 3/5 (1, 1), (3, 2), (2, 1), (4, 2)

1 -λ 1 Iteration 3 : 2 3/5 (1, 0), (2, 1), (4, 2) 1 + λ 1 (1, 0), (3, 1), (4, 3) 4 -4λ 1 (1, 0), (2, 1), (3, 2), (4, 3) 4 -4λ 1

  ISSN subm. to DMTCS c 2015 by the author(s) Distributed under a Creative Commons Attribution 4.0 International License 2 Discrete-time dynamic network A discrete dynamic network G = (N, A, T ) is a directed graph with N being a set of |N | = n nodes i, A being a set of |A| = m arcs a and T being the finite time horizon discretized into the set H = {0, 1, • • • , T }.

  (01) PmDF ALGORITHM; (02) BEGIN (03) find a feasible flow f (i, j; θ) in network Ḡ; (04) BP := {0}; k := 0; λ k := 0; (05) REPEAT (06) QDP(k, λ k , BP ); (07) k := k + 1; (08) UNTIL(λ k = Λ); (09) END. Algorithm 1: Parametric minimum dynamic flow (PmDF) algorithm.

  (a) The discrete-time dynamic network Ḡ. (b) The piecewise linear minimum flow over time value function for Ḡ.

Fig. 1 :

 1 Fig. 1: Illustrating PmDF algorithm.


Ali Fakra
email: fakra@univ-reunion.fr

Hamada

Harry Boyer

Eddy Lafosse

Philippe Lauret

THE NEW METHOD TO DETERMINATE INDOOR DAYLIGTHING THE BUILDING: IMPROVING "LUMEN MODELS"

. There are three method types to calculate the daylight luminance through windows and large openings inside the building. Among others, 'simplified methods' (''Lumen method'', ''Split Flux System method''…), 'numerical simulation method' (''Radiosity'', ''Ray Tracing''…) and finally 'experimental method' (''scale model'', ''building experimental ''…) are described in literature. We put some specific interest in the 'simplified methods' because they are easy to use and don't need a lot of data base for simulating. The big disadvantage for most of these models is these models depend on restricted conditions. For example, if we considered ''Split Flux System method'', we see that the models used only to overcast condition. Another example, the ''lumen method'' applied for any skies but the models take only five values in the space. Therefore they can only work under precise and limited circumstances. In this present paper, we purpose to see how we can circumvent these difficulties by improving the ''Lumen Model''. Indeed, we are going for introduce the notion of ''Equivalent Deep'' (E D ) (''profondeur fictive'' in French) into the ''Lumen method''. This ''Equivalent Deep'' will calculate any values and any types of skies contrary to the original ''Lumen method''. The original ''lumen method'' can only calculate indoors dayligthing in five values.

INTRODUCTION

We will present the methodology set up in order to improve the "Lumen model". The dayligthing within room is given by this model (*). In "Lumen method", the only parameters used for determinate the dayligthing within the room are the distance between a point enlightened and the window (the Depth from the room (D)), the coefficient of utilization from sky luminance (CU) which we will see more in detail, the net transmittance of the window factor (TL), the exterior vertical illuminance from the half sky, and finally the exterior horizontal illuminance from the ground. This method calculates dayligthing from five points on the workplane (a ratio of 10, 30, 50, 70 and 90 to 100 of the Depth of the room) for a depth D given (which represents either the length or the width of the room depending on the position of the window). The "lumen method" applied in every sky (cloudy, cloudless or intermediate sky), but we have not calculate for any points. In this works, we introduce a new concept who allows us to calculate the dayligthing on any points of workplane from "Lumen method". Indeed, the depth (D) of a point in workplane is absolutely essential to determine the coefficient CU ( [START_REF]Recommended Practice for the Lumen Method of Daylight Calculation Prepared by the IES Calculation Procedures Committee[END_REF] and [START_REF] Jone | Very simple Hand Calculations for Daylighting[END_REF]). This fixed depth limits (D) us the numbers of points on which one can calculate dayligthing in only five positions. Into our study we propose to introduce the "Equivalent Deep" to solve this problem. By definition, the "Equivalent Deep" will be the distance between window and an imaginary wall that would move along the real depth of the room. This imaginary wall will reduce the real depth of the room and will make it possible to recalculate dayligthing in five other points.

We will compare the results between the "Lumen method" and the new "improving Lumen method" developed.

CALCULATION PROCEDURE AND LIMIT FOR "LUMEN METHOD"

The "Lumen method" as adopted by the IES Committee on Calculation Procedure, (1989) is a procedure used to predict interior daylight illuminance under predetermined conditions [START_REF] Ies | Recommended Practice for the Lumen Method of Daylight Calculation[END_REF]. Most of the time, this method is used in the United States [START_REF] Jone | Very simple Hand Calculations for Daylighting[END_REF].

Principal of calculation

The prediction of interior daylight illumination from side lighting is simplified by using standard configuration shown in (Fig. 1). The floor cavity extends from the window sill to the floor and is assumed to have a reflectance of 30 %. The ceiling cavity extends from the top of the windows to the ceiling and is assumed to have a reflectance of 70 %. The room cavity extends in height (H) from the top of the floor cavity to the bottom of the ceiling cavity, in width (W) along the window wall, in depth (D) from the window wall to the rear wall, and its reflectance is assumed 50 %. The lumen method assumes 5 sky luminance distribution models, each with succeeding higher ratio of luminance of zenith to that of the horizon. The total interior daylight component is calculated for a few points 10%, 30%, 50%, 70% and 90% distance from the window. The sky models used and the inability to distinguish illuminance values of points at equal distance from window but unequal distance from wall are main deficiencies of the "Lumen method". The parameters used in this method are primarily the geometrical properties (width, height, positioning…) and physics (coefficient of luminous transmission, constituting pane…) of the glazed opening, (the date, the hour considered, the locality (longitude and latitude of the place) and) the weather conditions of illumination (light intensity provided by a weather file) and especially the depth of the room with respect to the position of the window. This method applies for any types of sky. The geometrical model used in the Lumen method for lighting in frontage ( [START_REF] Ies | Recommended Practice for the Lumen Method of Daylight Calculation[END_REF] and [START_REF]Estimating daylight in building : Part 1-2[END_REF] is given by figure (Fig. 1): 2.2 Limit for the method

    i xv E E CU
The "Lumen method" applied for any type of sky (clear, cloudy and intermediate sky), but he give only illumination values of five points of the workspace. Thus, we could not calculate the values of illumination everywhere on the workspace. In this paper, we has improved the "Lumen method" in order to calculate the dayligthing values for several points in workspace.

IMPROVING LUMEN METHOD

In this paragraph, we will present the step to be followed for "Improving Lumen method" and the calculation algorithm used. The programming language used is the C.

Approach of the proceeding

The new idea is to consider an "Imaginary moving wall" which would enable us to calculate illumination in any point of the room by using the "Lumen method ''. Indeed, starting from our ''Imaginary moving wall", we can recalculate the illumination on 5 other points of the room by keeping the same principle of calculation as the original ''Lumen method'. The only notorious difference will reside on the new points of calculations which different from those given by the real depth of the room. Let us take as example, the case where we want to know the value of illumination on a point which is on the perpendicular axis of the window at a distance from 3 meter. The room has 6 meters deep. The direct calculation of the "Lumen method" starting from the real depth of the room will not give us the value of illumination on this point. In order to be able to give the approximate value on this point (distant of 3m); we are thus to take a imaginary depth D virtuel = 3/0,9 for example (one could also choose D virtuel = 3/0,7 or divided by 0,5 or 0,3 or 0,1 this depth following the need) and to make calculations * by the "Lumen method". That way, we would have the value of desired illumination at the distance (i.e. 3m). Note that in this example, we consider that the imaginary wall has a imaginary depth of 3,33 m instead of the real depth (of 6 m). we are not interested in the values of illumination beyond 3 metres deep. Calculates them remain true As long as D virtuel ≤ D réel (what will be always the case as long as we seek to know illumination inside the room).

Fig. 2: proceeding for calculating sidelighting with improving "Lumen method"

Calculation algorithm of new method

We will see here how make algorithm code which allowed to Improving "Lumen method" (Fig. 3).

APPLICATIONS

Case study and suppositions of the simulation

In order to compare the results "Lumen method" and "improving lumen method", we will study the case from one simple room that has only one window. In this study, we consider only the dimension of the window (length and width of the window) and the depth from the room. The dimensions of the room (length and width of the room) are not necessary because the ''Lumen method'' don't use these characteristics. The dimensions of the window are 3 meters length 2 meters width. 20 meters depths zone.

Results and comparisons

We varied our depth from 0 to 20 meters. We then classified the values of illuminations obtained following these variations, according to five categories (Table 1). The model used to calculate these illuminations is well on the traditional ''Lumen method''.

We could note that the values remain the same ones for each category of depth: one thus has values of illuminations which do not change when one passes from 1 to 2 meters depth (Category A) and so on. This illustration (Fig. 4 and Fig. 5) shows us the values of illuminations according to our depths categories. The (Table 2) recapitulates the values of illuminations obtained according to the five points well defined by the Lumen method. These five points are calculated according to the various depths of the room. F and G The curve (Fig. 6) gives us the results of the values of illumination for given local points and various depths. We observed that with this method we can calculate values of illuminations in any point. This approach took into account each enlightened position points in comparison to windows of the room. Thus, we did not limit ourselves to the calculation of illumination into five point of the room only . The curve (Fig. 7) illuststrate the correlation of various values.

Fig. 7: Correlation between first values and second values

We note that, certains values of points near the window are not taken correctly because in the reality, the physical phenomenon is more complex and like we use a simplified model, it is obvious that these values will not be taken with accuracy. The linear line shows the correlation between illuminations values given by two depths rooms different for the same point.

CONCLUSIONS

The method lumen is a method simplified and efficacious to calculate illumination at the interior of a room according to any type of sky. The big defect of this method is that the calculations are limited in five positions of the room only.

We saw here how one could improve the model so that it can calculate any point of the room. The ''Improving method Lumen'' results obtained are promising.

Fig. 1 :

 1 Fig.1: Standard configuration in a room for calculating sidelighting (Lumen method)[START_REF] Ies | Recommended Practice for the Lumen Method of Daylight Calculation[END_REF]

Fig. 4 :Fig. 3 :

 43 Fig. 4: Case study (use DIALUX software for the picture): an simple room of various depths and one window

Fig. 4 :

 4 Fig. 4: Results of ''Lumen Method'' for only one window are and various depth room

Fig. 6 :

 6 Fig. 6: Comparison from illuminations values for various depths rooms

TABLE 1 : CATEGORY OF ILLUMINATIONS IN FIVE POINTS STARTING FROM THE LUMEN METHOD AND ACCORDING TO THE VARIOUS DEPTHS OF THE ROOM

 1

	A	B	C
	1 at 2 meters 3 at 4 meters 5 at 6 meters
	D	E	F
	7 at 8 meters 9 at 12 meters 13 at 16 meters
		G
		17 at 20 meters

TABLE 2 : ILLUMINATIONS VALUES IN VARIOUS POINTS OF THE ROOM TEST Local points (m) First values (lux)

 2

	Second
	values
	(lux)	Beetwen


Jean Philippe Praene
email: praene@univ-reunion.fr

Alain Bastide

Franck Lucas

François Garde

Harry Boyer

Simulation And Optimization of a Solar Absorption Cooling System Using Evacuated Tube Collectors

Summer air conditioning represents a growing market in commercial and residential buildings. Solar energy is a very interesting energy source because of its advantages. Instead of a compressor system, which uses electricity, an absorption cooling system needs only heat produced by a solar collector plant. Furthermore, Reunion Island has a high solar energy potential. The yearly average solar-radiation and is 5.4 kWh/m²/day. This paper presents a solar-powered, single stage, absorption cooling system, using a water-lithium bromide solution. The first part of this work deals with the dynamic modeling of an evacuated tube collector used for the simulation of heat production. In a second part, simulation and optimization of the system has been investigated in order to determine the optimum of solar collector plant surface, storage tank volume and nominal capacity of the absorption chiller. A new building code was added under TRNSYS16, to have dynamic coupling between building loads and cooling production.

INTRODUCTION

Reunion Island is a French overseas department located in the southern hemisphere characterized by a tropical humid climate. Conventional energy will not be enough to meet the continuously increasing need for energy in the future. In this case, renewable energy sources will become important. Solar energy is a very important energy source because of its advantages (insulation > 5 kWh/m²/day in Reunion Island). Of share our insularity an orientation towards renewable energies is today a perennial solution to obtain certain energy independence, cf. Figure 1. The building remains at the present time the most consuming sector in electricity. Under our latitudes, one of the most important demands for electricity relates to the air-conditioning of the buildings during the summer (period active of November 1 to April 30). Thus solar cooling system are particularly interesting as it is an application in which the demand for cooling energy closely matches the availability of solar energy, both in seasonal and the daily variations. One of the many categories of solar cooling systems is the solar absorption cooling. As no CFCs are used, absorption systems are friendlier to the environment. The possibility of cooling system using solar energy was initiated by the technological developments in the solar field, [START_REF] Tabor | Use of solar energy for cooling purposes[END_REF]. Absorption air-conditioning systems are similar to vapor compression air-conditioning systems, but differ in the pressurization stages. In general an absorbent in the low pressure side absorbs an evaporating refrigerant (H 2 O). The most usual combinations of chemical fluids used include lithium bromide-water (LiBr-H 2 O), where water vapor is the refrigerant, and ammonia-water (NH 3 -H 2 O) system where ammonia is the refrigerant, [START_REF]Handbook of fundamentals[END_REF]. The electric quantity of power consumed by the pump is almost negligible. One needs nevertheless a contribution of heat, cf. Figure 2. The single effect absorption chillers require for their operation a hot water temperature on the level of the generator of 80°C -150 °C. Moreover, the temperature of cooling water must lie between 7°C and 43°C, [START_REF]Handbook of fundamentals[END_REF]. The higher limit is founded in order to limit the differences in pressure between the generator and the absorber and the condenser and the evaporator. The lower limit makes it possible to avoid the crystallization of lithium bromide which is carried out at low temperature. An European research program, IEA Task 25, was set up in 1998, in order to study these installations. Two approaches are classically used within the framework of this working group. The first consists of the establishment of real pilots, one counts at the present time an about sixty systems installed, [START_REF] Henning | Solar-assisted Air-Conditioning in Buildings[END_REF]. The second consists to simulate and optimize existing or future systems via various codes. It is within the framework of second axis that this article is registered. The objective of this present work is to model the behaviour of the future installation at the Civil Engineering Department. Solar collectors represent the heart of the performances and the investments ( ~ 57%) of the installation. Thus, a good prediction of the performances of the solar loop is particularly judicious. The first part of this study was first interested in the dynamic modeling of the solar collectors. Then a coupling between the various components under TRNSYS has been investigate.

METHODS

There are three types of solar collectors which can be used within the framework of solar cooling: flat-plate, parabolic and vacuum collectors. Necessary heat for absorption chiller will be produced by evacuated tube collector. The literature contain numerous works on the modeling of solar collectors. These models developed have different levels of complexity. Usually, solar collectors are described by stationary models, considering the collector working under steady-state conditions. These approach are generally based on the work of Klein, [START_REF] Klein | Transient considerations of a flat-plate solar collectors[END_REF]. The details of the model developed and elements of validation associated have been previously described by Praëne, [START_REF] Praëne | Dynamic modelling and elements of validation of a solar evacuated tube collector[END_REF]. A comparison between model forecast and measurements is presented in Figure 5. The simulation is carried out at the minute time step. A new type has been created under TRNSYS, in order to coupling with the rest of the solar cooling system. The thermal behavior of the building was initially simulate over a whole year in order to observe cooling period. We used for that CODYRUN, [START_REF] Boyer | A multimodel approach to building thermal simulation for design and research purposes[END_REF], a thermal building code translated under TRNSYS to replace the TYPE56, this code has good accuracy for short time step. We evaluated over one year the level of comfort inside the classrooms and the building load to calibrate the absorption chiller refrigerating power.

We first simulate the evolution the inside classroom air temperature under natural ventilation during one year. The first condition of comfort is fixed for a temperature lower than 25°C. The graphic at Figure 6 shows that period spreading out May to October offers conditions of comfort. This period approximately corresponds to the southern winter. The period from November to April corresponds to the period with maximum needs of solar cooling system. In order to evaluate the building, a sequence of January has been used. An air temperature of 25°C and humidity of 55% are fixed as comfort criteria to dimensioning the air-conditioning system, cf. Figure 8. The maximum cold power called y the building is 9.5 kWf. This peak occurs after midday that is due to the accumulation of the walls inertia of the building and consequently restitution of this heat. Thus, for the four classrooms of the Civil Engineering Department, an absorption chiller of 35 kWf is necessary.

RESULTS AND DISCUSSION

This section presents the principal results obtain under TRNSYS16 for simulation of the global solar absorption cooling system. This program consists of the use of several "subroutines" which represents the components of the system described by ordinary or algebraic equations. Two types were added to the models existing under TRNSYS. It is about the model of simulation of the building and the model of vacuum collector. The useful surface area of collector used during simulation is 60 m ². The storage of hot water uses Type 38. That corresponds to a vertical roll of 800 L made up of copper insulated thermically with polyurethane. We consider a use of the machine with absorption without additional contribution of a contribution of auxiliary heat. We chose a sequence of January, to carry out our simulation, see Figure 9. The average temperature in classroom is 27°C during the day. There is difference about 4°C between inside and outside air temperature.

Figure 10. Evolution of cooling rate of the absorption chiller and outlet temperature of hot water storage tank.

The absorption chiller works without additional contribution of a contribution of auxiliary heat. As we can see on Figure 10, the 60 m² of ETC is enough to meet the minimum temperature of 80°C for the generator. One of the most important points from an economic point of view of a solar cooling plant is the solar loop. The field of solar collectors accounts for approximately 60% of the total investment (in particular if it is vacuum collector). So it is important to dimension in a first place the needs for the building then the total surface of the field of solar collectors. Three points go in general to the decision of the field of solar collectors: Ec onomic constraint: Budget dedicated to the project Spa ce constraint: surface available for the solar field Weather constraint: averag e numbers of good days during the hot period. Also, it is important to quantify the influence of the elements of the solar loop, of as much we are in the case of an operation of solar cooling, without auxiliary contribution. The effect of solar collector is evaluated against the energy from absorption chiller. The cooling rate from 60 m² to 40 m² is near. This production reaches to low value for a solar area less than 20 m². Thus, to have continuous condition for cold production a collector area of 60 m² is required.

CONCLUSION

The aim of the current work was to present a method of using solar energy for the airconditioning. The system is modeled with the TRNSYS program. For the requirements in cooling for our classrooms, an optimum of the components consists of the use of 60m ² of vacuum solar collectors associated a storage tank of 0.8 m 3 . We also fixed the nominal capacity of the absorption chiller with 33 kWf. The final objective of this work was to set up a simulation tool representing a solar cooling plant. This environment will be used as basic support with the simulation of various configurations. Theses simulations have been used to dimensioning an installation in our university. The cooling plant will be setup on may June 2007, in order to be in function at the beginning of summer in November 2007.

Figure 1 .

 1 Figure 1. Evolution of the repartition of energy production in Reunion.

Figure 2 .

 2 Figure 2. The basic principle of the absorption cooling system.

Figure 3 .

 3 Figure 3. Comparison of Evacuated tube / compound parabolic / flat-plate collector for Desiccant / Adsorption / Absorption cooling systems.

Figure 4 .

 4 Figure 4. Thermal networks of the ETC model

Figure 5 .

 5 Figure 5. Simulation of the outlet temperature of the solar collector model.

Figure 6 .

 6 Figure 6. Yearly evolution of air temperature in classroom.

Figure 7 .

 7 Figure 7. Yearly evolution of PMV* inside the classroom in natural ventilation condition.

Figure 8 .

 8 Figure 8. Sensible and latent power called by the classroom TD2 during January.

Figure 9 .

 9 Figure 9. Comparison between resultant temperature (Tres) in classroom and outside temperature (Tair).

Figure 11 .

 11 Figure 11. Effect of volume tank on the outlet of hot water temperature.

Proceedings of Clima 2007 WellBeing Indoors


Vinicius Mariano Gonc ¸alves

Philippe Fraisse

André Crosnier

Bruno Vilhena Adorno

Parsimonious Kinematic Control of Highly Redundant Robots

When a robot is highly redundant in comparison to the task to be executed, current control techniques are not "economic" in the sense that they demand, most of the time unnecessarily, all the joints to move. Such behavior can be undesirable for some applications. In this direction, this work proposes a new control paradigm based on linear programming that intrinsically provides a parsimonious control strategy, that is, one in which few joints move. In addition to a formal stability proof, the paper presents simulation and experimental results on the HOAP-3 humanoid robot. Finally, a comparison is made with a least-square method based on the pseudoinverse of the task Jacobian, showing that the proposed method indeed uses fewer joints than the classic one.

I. INTRODUCTION

The use of redundant robots is on the rise not only in industrial environments [START_REF] Fitzgerald | Developing baxter[END_REF] but also in human-robot interaction applications, where frequently humanoid robots or mobile manipulators are used. The robot kinematic control strategies are widely based on the pseudoinverse of the task Jacobian ( [START_REF] Liégeois | Automatic Supervisory Control of the Configuration and Behavior of Multibody Mechanisms[END_REF], [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], [START_REF] Chiaverini | Singularity-robust task-priority redundancy resolution for real-time kinematic control of robot manipulators[END_REF]) or, when there are constraints, usually convex quadratic programming solvers are used (see [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF], [START_REF] De Lasa | Feature-based locomotion controllers[END_REF], [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]).

Under redundancy, however, these strategies may generate movements that are not very economic in the sense that a huge amount of actuators must be needlessly activated to achieve the desired task. For instance, the traditional pseudoinverse approach generates a trajectory in which all joints move, but intuitively, in some cases a task represented by a certain amount of degrees of freedom (DOF) in the task-space can be solved with at most the same number of DOF in the joint space.

This paper is concerned with the development and analysis of a control paradigm that induces "economic" movements in face of highly redundant (overactuated) systems. The idea is to exploit the redundancy of the system to generate control laws that demand as few actuators as possible, which could provide several benefits. For instance, recent research suggests that parsimonious control can potentially provide more natural human-like movements [START_REF] Gauthier | A biomechanical inactivation principle[END_REF], [START_REF] Mehrdad | Fast Human Movements and Sparse Optimal Control Policies[END_REF], which could improve human-robot interactions. Moreover, parsimonious control also provides easier networked control by reducing the amount of control data exchanged during wireless communication (see [START_REF] Nagahara | Maximum Hands-Off Control: A Paradigm of Control Effort Minimization[END_REF] and the references therein) and, in addition, it could provide reduced feedback-induced negative effects, such as excitation of non-modeled dynamics. The This work has been supported by CAPES, CNPq, and FAPEMIG. latter is justified by the fact that in kinematic control the inertial parameters, which depend on the robot configuration, are not taken into consideration. Therefore, the less these configurations change the less the imprecise parameters will play a relevant role in the system and, consequently, the less will be the induced disturbance in the control. In this sense, some recent works, specially in model predictive control, exploit techniques borrowed from compressed sensing to generate control actions with the sparse property (see [START_REF] Gallieri | Lasso mpc: Smart regulation of overactuated systems[END_REF], [START_REF] Aguilera | Quadratic mpc with l0-input constraint[END_REF]).

A. Compressed sensing

The proposed approach is inspired on the recent developments in compressed sensing in the setting of signal processing (see [START_REF] Donoho | Compressed sensing[END_REF]). In this setting, it is desired to represent a signal as a linear combination of previously-defined basis functions, and the challenge is to represent the signal in the most memory-efficient way possible without losing too much (if any) information. The signal is then sampled in a way that the number of samples is much smaller than the number of basis signals. The values of the reconstructed signal and the original one must match in these samples and therefore the problem of finding the coefficients of the linear combination consists of solving a highly undetermined linear system Ax = b. There are many possible solutions to this equation, but it is of special interest to find one that is sparse, or parsimonious, which means a low number of non-zero coefficients. The reasoning for that in the signal processing setting is that memory can be saved by storing only the non-zero coefficients. Such sparse decomposition is often possible, as experience shows that usually practical signals have significant energy only in some sparse frequencies along the spectrum. There is a similar intuition in robotics: for highly redundant robots, tasks can often be completed without using all available joints, differently from what happens when conventional approaches are used, where usually all joints are used simultaneously.

The problem of compressed sensing and finding sparse control inputs for trajectory generation/control under the differential inverse kinematics approach have exactly the same mathematical nature, an undetermined system. Indeed, in the latter equations of the form J q = ėd appear, in which J is the task Jacobian, q is the target joint velocity (the control input) and ėd is a desired velocity in the task/operational space. When the system is redundant and the matrix J has a number of columns considerably greater than the number of rows, there is the possibility of finding sparse solutions q. In the specific context of robotics, the widely popular choice for solving this problem, which is based on the pseudoinverse of the Jacobian, most often generates nonparsimonious solutions. Although the mathematical structure used in the context of compressed sensing and in the control of redundant robots are essentially the same, the theoretical challenges are different. In the control context, for instance, there is the necessity to prove stability for the closed loop induced by this approach, an issue nonexistent in the signal processing setting.

Finding the smallest number of non-zero coefficients for the solution of a linear system Ax = b is NP hard (see [START_REF] Muthukrishnan | Data streams: Algorithms and applications[END_REF]). Mathematically, the formulation is given by

argmin x x 0 subject to Ax = b (1)
in which • 0 is the "norm" that counts the number of nonzero entries.

However, adequate solutions can be found very efficiently if relaxations are done. For instance, using a specific convex positive definite error metric between Ax and b, namely the 1-norm, the optimization problem can be reduced to a linear program. Linear programs can be solved by using the Simplex algorithm, which has the property of intrinsically generating reasonably sparse solutions, even though it may be not the minimum amount of nonzero entries as in Problem 1. This approach has the benefit that linear programs can be solved, in average, very efficiently with the Simplex method [START_REF] Shamir | The efficiency of the simplex method: A survey[END_REF]. The formulation will be valid even when it is not possible to have Ax = b (the system is not redundant), giving in this case the best possible solution in the 1-norm sense. Indeed, even when parsimony is not possible at all, the benefit of reduced computational time should remain.

B. Linear programming for inverse differential kinematics problems

It is not the first time that linear programming is proposed to solve inverse differential kinematic problems although, interestingly, none of the previous works neither made the connection with compressed sensing nor observed the potential of having sparse solutions. Formal guarantee of stability was also absent. Ho et al. [START_REF] Ho | Computing inverse kinematics with linear programming[END_REF] proposed a particular case of the general formulation proposed in the present paper, which aims at generating trajectories for human-like 3D models. They note that linear programming can be efficient, in terms of computation effort, to generate trajectories when the system has inequality constraints. Indeed, according to their simulations, in some settings even in the unconstrained case it outperforms the classic pseudoinverse approach by taking much less computational effort, and this is as truer as the number of constraints and degrees of freedom grow. The authors of [START_REF] Decre | Extending itasc to support inequality constraints and non-instantaneous task specification[END_REF] also propose that a general convex function can be used as a metric, and explicitly mention the • 1 norm, but also do not show simulation results with norms other than the Euclidean one.

C. Contributions

The main contribution of this paper is the formulation of the inverse differential kinematics as a linear program and the associated formal guarantees of Lyapunov stability. It is also shown that parsimonious behavior is obtained when the optimization problem is solved by using the Simplex method, which has low computational cost. Both simulations and experimental results in the HOAP-3 humanoid robot are shown to illustrate the proposed methodology.

II. LINEAR PROGRAMMING FOR ROBOT MOTION

GENERATION

For task-solving approaches using differential inverse kinematics and task functions ( [START_REF] Samson | Robot control. the taskfunction approach[END_REF]), control actions in the joint space often follow as a result of optimization problems. These optimizations problems seem to be, mostly, quadratic programming ones (see [START_REF] Liégeois | Automatic Supervisory Control of the Configuration and Behavior of Multibody Mechanisms[END_REF], [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], [START_REF] Chiaverini | Singularity-robust task-priority redundancy resolution for real-time kinematic control of robot manipulators[END_REF], [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF], [START_REF] De Lasa | Feature-based locomotion controllers[END_REF], [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF]), for which solutions can be computed explicitly when there are no constraints. To exemplify, let q be the variable in the joint space and e(q) = 0 a description of a desired configuration. Hence, the goal is to drive e(q(t)) to 0. Since ė(q)=J(q) q, setting the target ėd to -ηe(q), for η > 0, aiming for an exponential convergence of the error, one can obtain the target q(t) by solving the following optimization problem

min q J q -ėd 2 =⇒ min q J q + ηe 2 (2)
in which J J(q) and e e(q). The minimal 2-norm solution to the previous optimization problem can be given by using the pseudoinverse as q = -ηJ + e. However, as mentioned in Section I, this approach frequently implies a solution in which all joints move, since q in this case frequently has few, if any, non-zero entries.

In case that some constraints for q (e.g., speed limits or obstacle avoidance) are imposed in the form W(q) q ≤ w(q), dense solutions may arise in the solution of the problem min q J q + ηe 2 subject to W q ≤ w.

(

) 3
in which W W(q) and w w(q). Problem 3 has no closed-form solution and the problem, which is a convex quadratic program, needs to be solved with a numerical solver.

If parsimony is sought, an approach that can also potentially generate parsimonious behavior, when possible, will be discussed in the next section.

A. Parsimony

This section shows that the problem of inverse differential kinematics can be posed as a linear program, which in addition of being computationally efficient to solve, also provides parsimony as long as an appropriate solver is used. The problem can be seen as follows

min q J q + ηe 1 subject to W q ≤ w, (4)
where

q ∈ R m , J ∈ R k×m , e ∈ R k , W ∈ R s×m and w ∈ R s .
The formulation in Problem 4 by itself does not guarantee parsimony. However, since it can be transformed into a linear program, the well-known Simplex algorithm can be used to solve it and then, thanks to the intrinsic properties of that algorithm, parsimony appears whenever possible. To explain this behavior, let g, c ∈ R v , B ∈ R u×v , b ∈ R u , and c ∈ R. Consider the canonical form of linear programs

min g c T g + c subject to Bg = b; g ≥ 0.
(

) 5
It is important to note that all linear programs can be written in this form. If one uses the Simplex algorithm to solve linear programs, at most u entries of the v ones in g will be nonzero. This is due to how the Simplex algorithm works (see [START_REF] Murty | Linear Programming[END_REF]): at each step, the vector g is split in u basic variables and v -u nonbasic variables. All the nonbasic variables are set to 0 while the basic ones may or may not be 0. Thus, when possible, linear programming implements parsimony intrinsically when Simplex algorithm is used.

In order to use the Simplex algorithm, Problem 4 must first be transformed into a linear program. In statistics, this problem is known as constrained least absolute deviation regression and many different formulations as linear programs have been proposed. According to experiments done in [START_REF] Li | Solve least absolute value regression problems using modified goal programming techniques[END_REF], the most efficient formulation, for computational purposes, is the following one:

min ( q,y) 1 T (2y -(J q + ηe)) subject to J q -y ≤ -ηe; W q ≤ w; y ≥ 0, (6)
in which 1 is a vector of ones of appropriate dimension.

To write Problem 6 in the canonical form (5), q must be somehow rewritten, since all variables are nonnegative in [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF]. For handling this, let qP and qN be nonnegative variables. One can write without loss of generality that q = qP -qN . This decomposition is not only useful to transform Problem 6 into the canonical one, but also to avoid a "flaw" in the formulation of Problem 4: it does not guarantee that q = 0 when e = 0 (i.e., the task is achieved). This implies that when the robot achieves the task it may still continue to move its joints, which may be undesirable. Such a concern is nonexistent with the pseudoinverse approach because it implicitly minimizes q 2 and then, clearly, q = 0 when e = 0. There are many ways that this problem can be handled, but a simple one, and convenient in the linear programming case, is to choose a positive definite function β(e) and impose that q 1 ≤ β(e). This way, if e = 0 one has necessarily q = 0. This constraint can be implemented very simply as a linear constraint using the previouslyintroduced decomposition q = qP -qN . Indeed, clearly q 1 = qP -qN 1 ≤ qP 1 + qN 1 = 1 T qP + 1 T qN , because qP and qN are nonnegative. Thus, the constraint 1 T qP + 1 T qN ≤ β(e) implies q 1 ≤ β(e), as desired.

By introducing slack variables z

A , z B , z C , letting g = [ qT P qT N y T z T A z T B z C ] T
, and imposing the additional constraint 1 T qP + 1 T qN ≤ β(e), one can recast Problem 6 as

min g [-1 T J 1 T J 2•1 T 0 T 0 T 0]g -η1 T e subject to J -J -I I 0 0 W -W 0 0 I 0 1 T 1 T 0 0 0 1 g= -ηe w β(e)
g ≥ 0, (7) Problem 7 has the canonical form of the linear program [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF]. In this case, the number of rows of matrix B ∈ R u×v is u=k+s+1, which means that if the Simplex algorithm is used, at most k+s+1 entries of g will be non-zero and hence at most k+s+1 entries of ( qP -qN ) = q ∈ R m will be non-zero. This implies the sparsity property when m>k+s+1.

Note that in the previous analysis, from u=k+s+1 DOF that can possibly be used to solve the task, s+1 are used due to inequality constraints, in which s DOF are required by the s inequalities in W q ≤ w and one DOF is required by the constraint 1 T qP + 1 T qN ≤ β(e). However, it turns out that some of these constraints may be not overly stringent and be fulfilled with strict inequality when the problem is solved. In this case, no DOF for that particular inequality is taken and then the number of nonzero entries may be smaller. Such behavior is observed in simulation and experimental results, as shown in Section III.

It is important to stress that the previous result guarantees that q(t) always has at most k+s+1 non-zero entries, but it does not mean that during the entire trajectory the same k+s+1 entries of q(t) will be non-zero. In fact, along the trajectory different sets of k+s+1 variables may be non-zero, which means that the solution of (7) guarantees sparsity only locally, not globally. Nevertheless, as the experimental and simulation results of Section III show, this sparsity can be global and, if it is not, it may still be the case that there are some joints that are never used during the whole trajectory, achieving some partial global sparsity.

B. Stability analysis

A natural concern is about the closed-loop stability when the input q is generated by solving the optimization Problem 4 with the additional constraint q 1 ≤ β(e), where β is a positive definite function, which is equivalent to min q∈A(q) J q + ηe 1 [START_REF] Gauthier | A biomechanical inactivation principle[END_REF] with

A(q) = {p : Wp ≤ w} ∩ {p : p 1 ≤ β (e (q))} , (9)
in which W W(q) and w w(q). Considering the system as a pure integrator, q(t)= qd (t), in which the target velocity qd (t) is a minimizer of Problem 8, stability analysis can be performed. It is important to note that the target specification is written as e(q(t))=0. Also, q(t)= qd (t) is indeed a dynamical system that depends on q because when the optimization problem is solved in a given step by finding qd (t), this value is integrated (or, in the real robot, this velocity is applied to the joints and they move) and therefore q changes. When q changes, so does e(q), J(q) and the set A(q) (the latter changes because W(q), w(q) and β(e(q)) change). Hence, the parameters are changed in Problem 8 and it has to be solved again in the next loop.

However, the set of all minimizers qd (t) of Problem 8 is usually not a singleton (a set with a single member) so the dynamical system generated is not an ordinary differential equation of the form q(t)=f (q(t)), but a differential inclusion q(t) ∈ F(q(t)) since the function F : R m → 2 R m is set-valued (it is a correspondence between R m and R m ). More specifically, in the case of Problem 8,

F(q) = argmin p∈A(q) J(q)p + ηe(q) 1 . (10)
In order to study stability for this differential inclusion, it is necessary to establish three lemmas.

Lemma 1: Let e(q) and J(q) be continuous functions for all q. Let A(q) be continuous, non-empty, compact-valued and convex-valued for all q. Then, the correspondence F(q) in ( 10) is non-empty, convex-valued, compact-valued and upper semicontinuous.

Proof: The function (q, p) → J(q)p + ηe(q) 1 is clearly convex in p. Furthermore, due to the assumptions on J(q) and e(q) it is also continuous for all q and p.

Using this fact and the assumptions in A(q), the proof can be concluded directly from an extended form of Berge's maximum theorem found in [START_REF] Sundaram | A first course in optimization theory[END_REF]. 1Lemma 2: Let φ be a convex function, η a scalar, v a vector and e(q) a differentiable function with Jacobian J(q). Let V (q) φ(ηe(q)), and (∂φ)(u) be the subgradient of φ evaluated at u. Then, there exists a vector z L (q, v) ∈ (∂φ)(ηe(q)) such that lim inf

→0+ V (q + v) -V (q) ≤ η(J(q)v) T z L (q, v). (11)
Proof: It is known that for any convex function φ and vectors r, s (see [START_REF] Asplund | Gradients of convex functions[END_REF])

∀z ∈ (∂φ)(r) , φ(r) -φ(s) ≤ (r -s) T z. ( 12
)
Consider r=ηe(q + v) and s=ηe(q) in the previous expression, in order to conclude that ∀z ∈ (∂φ)(ηe(q + v)), V (q + v) -V (q) ≤ η(e(q + v) -e(q)) T z, [START_REF] Donoho | Compressed sensing[END_REF] in which the definition V (q) = φ(ηe(q)) was used. Take a positive sequence n , n ∈ Z + of 's that converges to 0, and a convergent sequence z n ∈ (∂φ)(ηe(q+ n v)). 2 Furthermore, divide both sides of (13) by n (the inequality will keep the direction because n > 0) and take the lim inf n→∞ , which also preserves the direction of the inequality. Hence

lim inf n→∞ V (q + n v) -V (q) n ≤ lim inf n→∞ η e(q + n v) -e(q) n T A z n B . (14)
The limit inferior in the right-hand side of ( 14) can be replaced by the ordinary limit and distributed between the two terms A and B because both limits of A and B exist (e(q) is differentiable and the sequence z n is convergent).

Using the chain rule for differentiation, the limit of A reduces to J(q)v since ∂e/∂q=J. In the calculation of the limit of B, one notes that subgradients of convex functions are upper semicontinuous correspondences (see [START_REF] Asplund | Gradients of convex functions[END_REF]). Since e is continuous, the correspondence (∂φ)(ηe(q)) is also upper semicontinuous on q. Therefore, since the sequence z n is convergent, lim n→∞ z n exists and, by the definition of upper semicontinuity of correspondences, this limit belongs to (∂φ)(ηe(q)) (remember that lim n→∞ q+ n v = q). Defining z L (q, v) lim n→∞ z n , the proof is concluded.

Lemma 3: If φ(u) is a convex function and (∂φ)(u) is its subgradient evaluated at u, then for any x and y ∀z ∈ (∂φ)(y) , x T z ≤ φ(x + y) -φ(y).

Proof: Let r = y and s = x + y in (12) and multiply both sides of the inequality by -1, and the proof is concluded.

The following result, established in [START_REF] Nakamura | Discontinuous control of nonholonomic systems using nondifferentiable Lyapunov functions[END_REF], is important to establish Lyapunov stability by using non-differentiable functions.

Theorem 1: (see [START_REF] Nakamura | Discontinuous control of nonholonomic systems using nondifferentiable Lyapunov functions[END_REF]) Let F(q) be a non-empty, compact-valued, convex-valued upper semicontinuous correspondence. If there exists a function V (q) which is positive definite to a set S (zero in the set and positive outside it) and locally Lipschitz such that ∀ q ∈ F(q) , lim inf

→0+ V (q + q) -V (q) ≤ 0 (16)
then the differential inclusion q ∈ F(q) is Lyapunov stable to S. The closed-loop stability is now proved in the next proposition.

Proposition 1: Consider the following differential inclusion q(t) ∈ argmin p∈B(q(t))

J(q(t))p + ηe(q(t)) 1 ,

q(0) = q 0 , (17)
2 Such convergent sequence always exists. Suppose there is a bounded sequence zn (clearly the sequence can be chosen to be bounded since all members of the subgradient are bounded) and it is not convergent. Then Bolzano-Weierstrass theorem guarantees that there exists a convergent subsequence of zn. This convergent subsequence can be taken along with the respective subsequence of n.

in which e(q) and J(q)=∂e(q)/∂q are continuous functions, e(q) is locally Lipschitz, η > 0, and B(q) is a continuous, compact-valued, convex-valued and non-empty correspondence with 0 ∈ B(q) for all q. Then the set S {q : e(q) = 0} is Lyapunov stable.

Proof: Let q be any particular member of the correspondence in the right side of [START_REF] Decre | Extending itasc to support inequality constraints and non-instantaneous task specification[END_REF], henceforth denoted by F(q).

From the definition of argmin, the optimal p = q from the feasible set B attains the minimum value of the objective function H(p)

Jp + ηe 1 (for the sake of notation, all dependencies on q will be dropped). In special, the value of H( q) is not greater than H(0), since p = q globally minimizes H and p = 0 is in the feasible set B by assumption. Consequently,

J q + ηe 1 H( q) ≤ ηe 1 H(0) =⇒ J q + ηe 1 -ηe 1 ≤ 0. (18)
Consider φ(u) = u 1 , which is a convex function. Apply Lemma 3 with x = J q and y = ηe, together with [START_REF] Samson | Robot control. the taskfunction approach[END_REF], to conclude that ∀z ∈ (∂φ)(ηe), (J q) T z ≤ J q + ηe 1 -ηe

1 ≤ 0. (19) Since η > 0, then ∀z ∈ (∂φ)(ηe) , η(J q) T z ≤ 0. ( 20
)
Considering Lemma 2 with v = q, and the fact that z L (q, q) ∈ (∂φ)(ηe), one can conclude together with (20) that lim inf →0+ V (q + q) -V (q) ≤ 0 [START_REF] Sundaram | A first course in optimization theory[END_REF] in which, according to the definition in Lemma 2, V (q)=φ(ηe(q)). It is important to note that (21) holds true for any q ∈ F(q) because in the beginning of the proof it was assumed that q is arbitrary. Now, consider V (q) = φ(ηe(q)) ηe(q) 1 as the Lyapunov function to the set S = {q : e(q) = 0}. This function is locally Lipschitz because both • 1 and e(q) are locally Lipschitz functions.

According to Lemma 1, F(q) is a non-empty, compactvalued, convex-valued and upper semicontinuous correspondence. Therefore, using Theorem 1 together with (21) (remember that it holds true for any q ∈ F(q)), one concludes that S is Lyapunov stable.

Note that the assumption 0 ∈ B(q)-which implies that the solution q = 0 is always admissible-is reasonable: it means that the robot is able to stop moving. The compactness is also reasonable, since the joint velocities always have natural limits and in practice all inequality constraints are non-strict.

Remark 1: Let B(q) = A(q), where A(q) is given in ( 9) with β being a positive definite continuous function. If w(q) ≥ 0 for all q, and W(q), w(q), J(q) and e(q) are continuous functions, the latter also being locally Lipschitz, then the assumptions in Proposition 1 hold. Indeed, the set A(q) is closed, since it is formed by non-strict linear inequalities on p. Furthermore, since the function β is continuous, and hence always finite, the constraint q 1 ≤ β(e(q)) by itself guarantees that the set is bounded. Thus A(q) is naturally compact. The property of being locally Lipschitz is also verified for the common functions used in kinematics models (sine, cosine, etc...). Consequently, the closed loop induced by solving the linear programming Problem 7, which is equivalent to Problem 8, is Lyapunov stable.

C. Computational complexity

The time needed for computing the control action can be crucial in real time applications on robotics, and thus the computational complexity of the proposed approach deserves special attention. Since the main burden of the method lies in solving Problem 7 with Simplex, this algorithm must be better analyzed.

According to [START_REF] Vanderbei | Linear Programming: Foundations and Extensions[END_REF], variants of the Simplex method take, on average, approximately a number of O(u + v) steps to solve linear programs as Problem 5, in which u and v are the number of constraints and variables, respectively. Taking into consideration that each step takes O(uv) operations [START_REF] Murty | Linear Programming[END_REF], this leads to the conclusion that the Simplex algorithm has in practice a complexity of O(uv

2 + u 2 v). Since in Problem 7 u = k + s + 1 and v = 2m + 2k + s + 1, 3 this implies in practice a time complexity of O((k + s)(2m + 2k + s) 2 + (k + s) 2 (2m + 2k + s)), which is polynomial.
Since the pseudoinverse approach cannot handle inequality constraints, in order to compare it with LP, let s = 0. Furthermore, in this paper only the redundant case is relevant; therefore, considering that the number m of DOFs is considerably greater than the number k of tasks, the practical complexity of LP reduces to O(m 2 k).

In comparison, the pseudoinverse of a matrix J ∈ R k×m can be also computed on O(m 2 k) (under the consideration that m is much greater than k) using, for instance, a singular value decomposition (see [START_REF] Golub | Matrix Computations, ser. Johns Hopkins studies in the mathematical sciences[END_REF]). Therefore, the asymptotic complexity of both approaches is comparable. However, it has been reported in the literature that LP can outperform the pseudoinverse (see [START_REF] Ho | Computing inverse kinematics with linear programming[END_REF]), specially when the Jacobian is sparse.

When there are inequality constraints, i.e. s = 0, algorithms based on pseudoinverse cannot be used and convex quadratic programs as Problem 3 must be solved with numerical solvers. Since convex quadratic programming is a strict superset of linear programming, the former is expected to be, in average, at least as hard to solve as the latter.

III. EXPERIMENTAL AND SIMULATION RESULTS

A. Simulation results on highly redundant robot

To illustrate the method, consider a serial 8-link planar robot in which each link has unit length and unit mass. The task consists of moving the end-effector to the point [1 7] T while keeping the x component of the robot's center of mass (COM) in 0. Clearly, the robot is highly redundant in comparison with the desired task and there are many possible movements. More specifically, the task Jacobian has 3 rows and 8 columns.

Both pseudoinverse and linear programming (Problem 7) approaches are compared in the solution to this task, with common parameter η = 0.5, which implies the same convergence rate for both methods. The only constraint in Problem 7 is the limitation of the speed with β(e) = 5 e 1 .

Fig. 1 shows the robot's configurations for both approaches. The initial configuration (cyan) corresponds to the manipulator being entirely in the vertical with the exception of the last joint that is tilted 30 degrees clockwise. Intermediate configurations are shown in intermediate colors (lighter is closer to the initial configuration and darker is closer to the final one). The final configuration is purple.

In the linear programming approach, at any specific time only three joints moved, whereas in the total movement four joints moved (1,4,5 and 8, see Fig. 2). More specifically, there is a switch between joints 4 and 5 at around 70 samples (see Fig 2). In the pseudoinverse approach, all of them moved during the whole trajectory. Table I shows a comparison of both approaches under four different metrics. For each task, the error dynamics in the first and second experiment are exactly the same, so the comparison is fair.

Since m=8, k=3 and s=0 (no constraints), the developments in Subsection II-A predict that for the linear programming approach, at a given time at most k+s+1=3+0+1=4 Fig. 2: Comparison of the control inputs q from the pseudoinverse (bottom) and linear programming (top). The joints are numbered so the first joint is nearest to the base and the last one is nearest to the end-effector (see Fig. 1).

joints should move. In this case, only three joints moved at a given time (see Fig. 2). This is expected if the function β(e) is not overly stringent, so the bound is always achieved with strict inequality and no degree of freedom is used to fulfill the constraint in the last row in Problem 7. The function defined for this simulation, β(e) = 5 e 1 , was loose enough for this to happen.

B. Experimental results on humanoid robot HOAP-3

To further illustrate the proposed methodology, an experiment was performed in the HOAP-3 humanoid robot (see Fig. 3). In this experiment, initially the robot is standing with the right foot on the ground and the other one raised (see Fig. 4). The position is statically stable, since the (x, y) coordinates of the COM, which is located at m 0 = [-0.8cm -6cm] T with respect to the reference frame of the right foot, is in the support polygon.

The task consists of moving the (x, y) coordinates of the COM to m d = [-3.5cm -3.6cm] T , which is also in the support polygon and defined in a reference frame located at the right foot. The desired task is e(q) = m(q) -m d = 0. The experiment was performed with two approaches: the classic pseudoinverse of the Jacobian, q = -ηJ + e, and the formulation in Problem 7. In the latter case, the only constraint was defined as β(e) = 40 e 1 . In both approaches, η = 0.1 to guarantee the same convergence rate.

Furthermore, the trajectories of the robot joints were generated offline, both for the pseudoinverse approach and LP, and then executed in the real robot. Therefore, only the trajectory tracking of the joints was executed in closed loop. The implementation of the LP algorithms in the robot and consequently closing the entire loop will remain for a future development.

Fig. 5a and Fig. 5b show the required input q in both scenarios: pseudoinverse (PINV) and linear programming (LP). Since the Jacobian matrix has two rows (k = 2) and there are no constraints (s = 0) other than the natural one of velocity, it was claimed in Subsection II-A that at most k+s+1=2+0+1=3 joints will move in a given moment in the LP approach. Indeed, only two joints moved at all, because the natural velocity constraint was not overly stringent and therefore was never saturated (always achieved with strict inequality) to require a DOF. Furthermore, global sparsity was achieved since at all time the same two joints, namely 5 and 6 (see Fig. 3 for the joint convention), move during the entire trajectory.

On the other hand, in the pseudoinverse approach all joints move, as expected. Some of them moved very little, so if one considers as "effectively moving" only those joints that have a velocity above 5% of the maximum velocity in that given instant, then during the whole time only seven joints effectively moved.

Table II shows the comparison of the demanded input q according to four different metric for both approaches. As before, the error dynamics are exactly the same for both approaches, and therefore the comparison is fair.

It is interesting to note that in the case of LP, the two joints that move (5 and 6, see Fig. 3) are the roll and pitch of the ankle of the right feet, the one on the ground. (b) LP approach: (joints 8 to 25 were omitted).

Fig. 5: PINV vs. LP: demanded velocities for each joint; omitted joints do not contribute to the motion (see Fig. 3 for the joints convention). This agrees with the human movement observed in postural control. Indeed, it has been shown by [START_REF] Nashner | The organization of human postural movements: A formal basis and experimental synthesis[END_REF] that in the case of a small disturbance, the postural system will recover balance in the sagittal plane by using mainly the ankle joints. Finally, it is important to highlight that, although the joints trajectories were generated offline, it is feasible to generate them in real time. In order to support this claim, numerical simulations were done in an Intel i5 GHz with 4 GB of RAM. Among 500 samples, the pseudoinverse took in average 0.2 milliseconds to be computed, whereas each linear program with Simplex took in average 1.2 milliseconds to be solved. If warm start is used in the Simplex (the solution found in the previous step is used as a feasible basis in the current one), the time in Simplex also reduces to 0.2 milliseconds in average, because in the huge majority of cases the initial basis is primal and dual feasible (thus both feasible and optimal) and only part of the first iteration is necessary in the Simplex method. This happens because the problem has a continuous nature and thus at each step the parameters of the optimization problem change little.

Given that the processor of HOAP-3 is capable of using the pseudoinverse approach in real time with a wide margin (see [START_REF] Adorno | Towards a cooperative framework for interactive manipulation involving a human and a humanoid[END_REF]), the timings presented here imply that it could also be capable of using the LP approach, specially if warm start is implemented.

IV. CONCLUSION AND FUTURE DEVELOPMENTS

This paper proposed a new paradigm to control highly redundant robots. The technique is based on a linear programming formulation, which when solved using the Simplex algorithm intrinsically reduces the number of non-zero components in the control vector. Formal results of stability were presented, and two examples-one in simulation and the other in a real humanoid robot-were implemented to illustrate the methodology.

The natural step for next works is the generalization of the proposed approach to hierarchy of tasks, as in [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], [START_REF] Kanoun | Kinematic control of redundant manipulators: generalizing the task priority framework to inequality tasks[END_REF], [START_REF] De Lasa | Feature-based locomotion controllers[END_REF], [START_REF] Escande | Hierarchical quadratic programming: Fast online humanoid-robot motion generation[END_REF], among others. In that case, the 1-norm is attractive because the associated optimization problemsi.e., lexicographical linear programs-can be solved by a very simple variant of the Simplex algorithm, as pointed out in [START_REF] Isermann | Linear lexicographic optimization[END_REF]. Furthermore, since the Simplex algorithm is very well understood and efficient, the lexicographical Simplex is also expected to be efficient because it can borrow many advances from its traditional counterpart. However, a theoretical challenge would be to prove stability for the hierarchical case, which would consist of a generalization of Proposition 1.

Finally, it turns out that Proposition 1 can also be generalized to encompass tasks in which it is not only desirable to converge to a set of configurations (e.g., a circle), but also to keep forever circulating in it [START_REF] Goncalves | Vector fields for robot navigation along time-varying curves in n -dimensions[END_REF]. This self sustained motion-a limit cycle in the configuration space-could be used to accomplish, for instance, a dance-like motion or even walking [START_REF] Hobbelen | Limit Cycle Walking. Humanoid Robots: Human-like Machines[END_REF].

Fig. 1 :

 1 Fig. 1: Comparison of the sequence of configurations with pseudoinverse (left) and linear programming (right).

Fig. 3 :

 3 Fig. 3: HOAP-3 humanoid robot (left) and convention for numbering the DOF used in this paper (right).

Fig. 4 :

 4 Fig. 4: Initial configuration of the robot for the experiment.

  PINV approach (joints 22 to 25 were omitted). Time (s) Angular speed (rad/s)

TABLE I :

 I Simulation of a serial 8-link planar robot: comparison of the input with respect to four different metrics.

	Metric	PINV	LP	PINV/LP
	∞ 0	q(t) 1 dt	1.644	1.264	1.301
	∞ 0	q(t) 2 2 dt	0.503	0.638	0.789
	∞ 0	q(t) 1 dt	0.021	0.016	1.329
	∞ 0	q(t) 2 2 dt	0.011	0.011	1.000

TABLE II :

 II Experiment with HOAP-3: comparison of the input with respect to four different metrics.

	Metric	PINV	LP	PINV/LP
	∞ 0	q(t) 1 dt	0.185	0.131	1.412
	∞ 0	q(t) 2 2 dt	0.040	0.067	0.597
	∞ 0	q(t) 1 dt	0.101	0.072	1.400
	∞ 0	q(t) 2 2 dt	0.289	0.339	0.852

Vinicius Mariano Gonc ¸alves ( mariano@cpdee.ufmg.br) and Bruno Vilhena Adorno (adorno@ufmg.br) are with Universidade Federal de Minas Gerais (UFMG), Brazil.

André Crosnier (crosnier@lirmm.fr) and Philippe Fraisse (fraisse@lirmm.fr) are with Laboratoire d'Informatique, de Robotique et de Microélectronique de Montpellier (LIRMM), France.

The theorem deals with maximum of functions, but an analogous to minima easily follows.

Recall that qP , qN ∈ R m , y, z A ∈ R k , z B ∈ R s and z C ∈ R.