Daniel Leivant
email: leivant@indiana.edu

Jean-Yves Marion
email: jean-yves.marion@loria.fr

Evolving graph-structures and their implicit computational complexity

Dynamic data-structures are ubiquitous in programming, and they use extensively underlying directed multi-graph structures, such as labeled trees, DAGs, and objects. This paper adapts well-established static analysis methods, namely data ramification and language-based information flow security, to programs over such graph structures. Our programs support the creation, deletion, and updates of both vertices and edges, and are related to pointer machines. The main result states that a function over graph-structures is computable in polynomial time iff it is computed by a terminating program whose graph manipulation is ramified, provided all edges that are both created and read in a loop have the same label.

Introduction

The interplay of algorithms and data-structures has been central to both theoretical and practical facets of programming. A core method of this relation is the organization of data-structures by underlying directed multi-graphs, such as trees, DAGs, and objects, where each vertex points to a record. Such data structures are often thought of as "dynamic", because they are manipulated by algorithms that modify the underlying graph, namely by creating, updating and removing vertices and edges. Our imperative language is inspired by pointer machines [START_REF] Kolmogorov | On the definition of an algorithm[END_REF][START_REF] Schönhage | Storage modification machines[END_REF] and by abstract state machines [START_REF] Gurevich | Sequential abstract state machines capture sequential algorithms[END_REF].

In this work we propose a simple and effective static analysis method for guaranteeing the feasible time-complexity of programs over many dynamic datastructures. Most static analysis efforts have focused in recent years on program termination and on safety and security. Our work is thus a contribution to another strand of static analysis, namely computational complexity.

Static analysis of computational complexity is based on several methods, classified broadly into descriptive ones (i.e. related to Finite Model Theory), and applicative (i.e. identifying restrictions of programs and proof methods that guarantee upper bounds on the complexity of computation). One of the most fruitful applicative methods has been ramification, also referred to as tiering. Initially this method was used for inductive data, such as words and natural numbers, but lately the method has been applied to more general forms of data. The intuition behind ramification is that programs' execution time depends on the nature of information flow during execution, and that such flow can be constrained naturally and effectively by imposing a precedence relation regulating the information flow, e.g. from higher tier to lower tier. Here we refer to a ramification that pertains to commands of our imperative language as well as to expressions denoting graph elements.

Our main result states that a function over graph-structures is computable in polynomial time iff it is computed by a terminating program whose graph manipulation is ramified, provided all edges that are both created and read in a loop have the same label. This result considerably extends previous uses of ramification in implicit computational complexity, and this extension touches on some of the most important aspects of programming, which have been disregarded in previous research. A simple modification of the proof gives an analogous characterization of Logspace computation over graph structures. 3 We thus believe that this work is of both theoretical and practical significance. Our results raise interesting questions about relations between data ramification and typed systems for program security [START_REF] Volpano | A sound type system for secure flow analysis[END_REF][START_REF] Sabelfeld | Declassification: dimensions and principles[END_REF], where the concept of information flow is explicit.

2 Evolving graph-structures

Sorted partial structures

The framework of sorted structures is natural for the graph-structures we wish to consider, with vertices and data treated as distinct sorts. Data might itself be sorted, but that does not concern us here. Recall that in a sorted structure, if V and D are sorts, then a function f is of type V → D if its domain consists of the structure elements of sort V, and its range of elements of sort D.

The graphs we consider are essentially deterministic transition graphs: edges are labeled (to which we refer as "actions"), and every vertex has at most one out-edge with a given label. Such graphs are conveniently represented by partial functions, corresponding to actions. An edge labeled f from vertex u to vertex v is represented by the equality f (u) = v. When u has no f-labeled out-edges, we leave f (u) undefined. To represent function partiality in the context of structures, we post a special constant nil, assumed to lie outside the sorts (or, equivalently, in a special singleton sort) 4 . We write f : V D, and say that f is a partial function from sort V to sort D, if f : V → (D ∪ {nil}). We write V D for the type of such partial functions.

Graph structures

We consider sorted partial structures with a distinguished sort V of vertices. To account for the creation of new vertices, we also include a sort R of reserved-vertices. For simplicity, and without loss of generality, we assume only one additional sort D, which we dub data. A graph vocabulary is a sorted vocabulary for these sorts, where we have just five sets of identifiers: V (the vertex constants), D (the data constants), F (function-identifiers for labeled edges, of type V V), G (functionidentifiers for data, of type V → D), and R (the relation-identifiers), each of some type

τ × • • • × τ , where each τ is V or D. As syntactic parameters we use v ∈ V, d ∈ D, f ∈ F, g ∈ G, and R ∈ R.
Given a sorted vocabulary Σ as above, a Σ-structure S consists of a vertexuniverse V S which is finite, a potentially infinite reserve-universe R S , a datauniverse D S , a distinct object ⊥ to interpret nil, and a sort-correct interpretation

A S for each Σ-identifier A: v S ∈ V S ; d S ∈ D S ; g S ∈ [V S → D S ] (a data- function); f S ∈ [V S
V S ] (a partial function), and for a relation-identifier R, a relation R S ⊆ τ S × • • • τ S . Note that we do not refer to functions over data, nor to functions of arity > 1. Also, the fact that our graphs are edge-deterministic is reflected in our reprensetation of edges by functions.

Our graph structures bear similarity to the struct construct in the programming language C, and to objects (without behaviors or methods): a vertex identifies an object, and the state of that object is given by fields that are specified by the unary function identifiers. This is why Tarjan, in defining similar structures [START_REF] Tarjan | Reference machines require non-linear time to maintain disjoint sets[END_REF], talks about records and items rather than vertices and edges. The restriction of a graph-structure S to the sort V of vertices can be construed as a labeled directed multi-graph, in which there is an edge labeled by f from vertex u to vertex v exactly when f S (u) = v. Thus the fan-out of each graph is bounded by the number of edge-labels in Σ. Examples of graph structures abound, see examples in Section 5. Linked-lists of data is an obvious one, of which words (represented as linked lists of alphabet-symbols) form a special case.

Expressions

Expressions are generated from a set X of vertex-variables, a set Y of datavariables, and the vocabulary identifiers, as follows. Equality here does not conform strictly to the sort discipline, in that we allow equations of the form V = nil5 .

V ∈ VExpr ::= X | nil | v | f (V ) where X ∈ X D ∈ DExpr ::= Y | d | g(V ) where Y ∈ Y B ∈ BExpr ::= V = V | D = D | ¬(B) | R(E 1 . . . E n ) where R : τ n , E i : τ
3 Imperative programs over graph-structures

Programs

We refer to a skeletal imperative language, which supports pointers:

P ∈ Prg ::= X:=V | Y :=D | f (X):=V | g(X):=D | New(X) | skip | P ; P | if (B) {P } {P } | while (B) {P }
The boolean expression B of conditional and iterative commands is said to be their guard. We posit that each program is given with a finite set X 0 ⊂ X of input variables.

Example: Tarjan's Union Algorithm

The following graph-algorithm is due to Tarjan [11, p.21]. It refers to a representation of sets by linked-lists, whose initial vertex also serves as a name for the set.

The linked-list is represented by a partial-function next, and the function parent maps each node to the head of its linked-list. The algorithm generates, given as input two lists r, q representing disjoint sets, a list representing their union. It successively inserts right after r's head the entries of q; thus r is maintained as the name of the union.

while ( q = nil) { s a v e := next ( q ) ; parent ( q ) := r ; % parent and next a r e m o d i f i e d next ( q ) := next ( r ) ; next ( r ) := q ; q := s a v e }

We shall see that our tiering method admits the program above.

Evolving structures

In defining the semantics of an "uninterpreted" imperative program one refers to structures for that program's vocabulary, augmented with a store (i.e. environment, valuation). For programs over a Σ-structure S, a store consists then of a function µ = µ X ∪ µ Y , where µ X : X → V S ∪ {⊥}, and µ Y : Y → D S . A Σ-configuration is a pair consisting of a Σ-structure S and a store µ. We chose the phrase configuration to stress their dynamic nature, as computation stages of an evolving structure.

Commands of imperative languages are interpreted semantically as partialfunctions that map an initial configuration to a final one. Here we have two commands whose intended semantic interpretation is to modify the structure itself. A command New(X) modifies the sorts, by moving an element of R S to V S , and updating the store to have X point to the new element. We write (S, µ)[νX] for the resulting configuration. Thus X points to a fresh vertex. More importantly, a command of the form f (X) := V modifies the semantics of the partial-function f S , in that for a = µ(X) and w the value of V in (S, µ) (defined formally below), an f-edge v u is replaced by v w. We write (S, µ)[f (X) ← w] for the resulting interpretation.

Our structure updates are obviously related to Gurevich's abstract sequential machines (ASM) [START_REF] Gurevich | Sequential abstract state machines capture sequential algorithms[END_REF]. Gurevich divides identifiers into two classes: static identifiers, whose interpretation remains constant, and dynamic identifiers, whose interpretation may evolve during computation. Here the only dynamic identifiers are the edge functions. An ASM computation progresses through "states", where every state is a structure. In contrast, we refer to configurations, because program variables play a central role in our imperative programs. Thus, the execution of our programs progresses through configuration.

Our programming language is related, more broadly, to pointer machines. Tarjan [START_REF] Tarjan | Reference machines require non-linear time to maintain disjoint sets[END_REF] defined a pure reference machine, consisting of an expandable collection of records and a finite number of registers. Pure reference machines are easily simulated by our programs, and vice versa.

Our programs are also related to Schönhage machines [START_REF] Schönhage | Storage modification machines[END_REF]. Each such machine consists of a finite control program combined with a dynamic structure (which is essentially the same as our graph-structures). Schönhage machines are an extension of Kolmogorov-Uspensky machines [START_REF] Kolmogorov | On the definition of an algorithm[END_REF]. Another source of inspiration is the work of Jones & als. on blob model of computations [START_REF] Hartmann | Programming in biomolecular computation: Programs, self-interpretation and visualisation[END_REF].

Semantics of expressions

We give next the evaluation rules for Σ-expressions E in a Σ-configuration (S, µ), writing S, µ |= E e ⇒ a to indicate that E evaluates to element a of S.

6 b ∈ V ∪ D S, µ b e ⇒ b S Z ∈ X ∪ Y S, µ Z e ⇒ µ(Z) S, µ E e ⇒ a h ∈ F ∪ G S, µ h(a) e ⇒ h S (a) S, µ E i e ⇒ a i a 1 ..a n ∈ R S S, µ R(E 1 , . . . , E n ) S, µ E i e ⇒ a i a 1 ..a n ∈ R S S, µ ¬R(E 1 , . . . , E n )

Semantics of programs

The semantics of programs is defined below:

S, µ E e ⇒ a S, µ Z:=E s =⇒ S, µ[Z ← a] skip S, µ New(X) s =⇒ (S, µ)[νX] skip S, µ X e ⇒ a S, µ V e ⇒ b S, µ f (X):=V s =⇒ S, µ[f(a) := b] skip S, µ P1 s =⇒ S , µ P 1 S, µ P1; P2 s =⇒ S , µ P 1 ; P2 S, µ P1 s =⇒ S , µ skip S, µ P1; P2 s =⇒ S , µ P2 S, µ B S, µ if (B){P0}{P1} s =⇒ S, µ P0 S, µ ¬B S, µ if (B){P0}{P1} s =⇒ S, µ P1 S, µ B S, µ while(B){P } s =⇒ S, µ P ; while(B){P } S, µ ¬B S, µ while(B){P } s =⇒ S, µ skip
The phrase S, µ P s =⇒ S , µ P conveys that evaluating a program P starting with configuration (S, µ) is reduced to evaluating P in configuration (S , µ ); i.e., P reduces to P while updating (S, µ) to (S , µ ).

An initial configuration is a configuration (S, µ) where µ(X) = nil for every non-input variable X. A program P computes the partial function [[P ]] with initial configurations as input, defined by: [[P ]](S, µ) = (S , ξ) iff S, µ P ( s =⇒) * T , ξ skip.

Run-time

We say that a program P runs in time t on input (S, µ), and write Time P (S, µ) = t, when S, µ P ( s =⇒) t T , ξ skip for some (T , ξ). We take the size |S, µ| of a configuration (S, µ) to be the number n of elements in the vertex-universe V . Since the number of edges is bounded by n 2 , we disregard them here. We also disregard the size of the data-universe, because our programs do not modify the data present in records. A program P is running in polynomial time if there is a k > 0 such that Time P (S, µ) k • |S, µ| k for all configurations (S, µ), 4 Ramifiable programs

Tiering

Program tiering, also referred to as ramification, has been introduced in [START_REF] Leivant | Predicative recurrence and computational complexity I: Word recurrence and poly-time[END_REF] and used in restricted form already in [START_REF] Bellantoni | A new recursion-theoretic characterization of the poly-time functions[END_REF]. It serves to syntactically control the runtime of programs. Here we adapt tiering to graph-structures. The main challenge here is the evolution of structures in course of computation. To address it, we consider a finite lattice T = (T, , 0, ∨, ∧), and refer to the elements of T as tiers. However, in order to simplify soundness proofs, and without loss of generality, we will focus on the boolean lattice T = ({0, 1}, , 0, ∨, ∧). We use lower case Greek letters α, β as discourse parameters for tiers.

Given T, we consider T-environments (Γ , ∆). Here Γ assigns a tier to each variable in V, whereas ∆ assigns to each function identifier f : V V one or several expressions of the form α → β, so that either 1. all types in ∆(f ) are of the form α → α, in which case we say that f is stable in the environment; or 2. all types in ∆(f ) are of the form α → β, with β ≺ α, and we say that f is reducing in the environment.

A tiering assertion is a phrase of the form Γ , ∆ V : α, where V is a vertexexpression and (Γ , ∆) a T-environment. The correct tiering assertions are generated by the tiering system in Figure 1.

Ramifiable programs

Given a lattice T, a program P is T-ramifiable if there is a T-environment (Γ , ∆) for which Γ, ∆ P : α for some α, and such that Γ (X ) = 1 for every input variable X ∈ X 0 . 7 Thus, ramifiable programs can be construed as programs decorated with tiering information.

Lemma 1 (Subject reduction).

If S, µ P s =⇒ S , µ P and Γ, ∆ P : α then Γ, ∆ P : α.

Lemma 2 (Type inference). The problem of deciding, given a program P and a lattice T, whether P is T-ramifiable, is decidable in polynomial time. 7 Recall that each program is assumed given with a set X0 of input variables.

Γ, ∆ c : α Γ (X ) = α Γ, ∆ X : α α → β ∈ ∆(f) Γ, ∆ V : α Γ, ∆ f(V ) : β Γ, ∆ Vi : α Γ, ∆ R(V1, . . . , Vn) : α Γ , ∆ Vi : α Γ , ∆ V0 = V1 : α
Fig. 1. Tiering rules for vertex and boolean expressions Proof. We associate with each vertex-variable X a "tier-variable" α X , and with each function f ∈ F two variables α f and β f , with the intent that α f → β f is a possible tiering of f. The typing rules for tiers give rise to a set of linear constraints on these tier-variables, a problem which is poly-time decidable.

Γ, ∆ X : α Γ, ∆ V : α Γ, ∆ X :=V : α Γ, ∆ f(X) : α Γ, ∆ V : α Γ, ∆ f(X ):=V : α Γ, ∆ X : 0 Γ, ∆ New(X) : 0 Γ, ∆ B : α Γ, ∆ P : α 0 ≺ α Γ, ∆ while(B){P } : α Γ, ∆ skip : 0 Γ, ∆ P : α Γ, ∆ P : β Γ, ∆ P ; P : α ∨ β Γ, ∆ B : α Γ, ∆ Pi : α Γ, ∆ if (

Stationary and tightly-modifying loops

We say that a function identifier f is probed in P if it occurs in P either in some assignment X := V or in the guard of a loop or a branching command. For example, f is probed in X := f(V ), as well as in if (f(X ) = nil){P }{P }. The identifier f is modified in P if it occurs in an assignment f (X) := V in P . Fix a lattice T, and a T-environment (Γ , ∆). By the tiering rules, if a loop while(B){P } is of tier α then Γ , ∆ B : α. We say that the loop is stationary if no f ∈ F of type α → α is modified therein. The loop above is tightly-modifying if it has modified function-identifiers of type α → α, but at most one of those is also probed. In other words, all edges that are both created and read in a loop have the same label. For instance, in Example 3.2 above, next is both modified and probed, but parent is modified without being probed. Thus the loop, with its obvious tiering environment, is tightly-modifying. 8

Fig. 2 .

 2 Fig. 2. Tiering rules for programs

We also present extensions of our language to incorporate recursion.

This is a minor, but deliberate, departure from the usual ontological (i.e. Church-style) typing of sorted structure.

Negation is useful, however, in defining commands' semantics and dispensing with truth constants.

Here we consider equality as just another relation.

Note that next and parent are of type 1 → 1, and all variables are of tier 1. Set union can be iterated because the result r is of tier 1, unlike in most other works.

The construct forall X in R(u), which is "blind," in the sense that it does not depend on node ordering. As a result, no function identifier is probed except visited.

Main Characterization Theorem

Given a lattice T and Γ , ∆ P : α, we say that (Γ , ∆) is a tight ramification of P if Γ is an initial tiering, and each loop of P is stationary or tightly-modifying. We say that P is tightly-ramifiable if it has a tight T-ramification (Γ , ∆), with Γ initial, for some non-trivial T.

Theorem 1. A function over graph-structures is computable in polynomial time iff it is computed by a terminating and tightly-ramifiable program.

The Theorem follows from the Soundness Lemma 6 and the Completeness Proposition 1 below.

Examples of ramified programs

Tree insertion. The program below inserts the tree T into the binary search tree whose root is pointed-to by x. The input variables are x and T . i f ( x 1 = n i l ) {x 1 :=T : 1 ; } % then c l a u s e { % e l s e c l a u s e while ( ( x 1 = n i l ) and ( key (T 1 ) = key ( x 1 ) ) ) { i f ( key (T 1 ) < key ( x 1 ) ) {p 1 :=x

Note that neither left nor right is modified in the loop, so the loop is stationary.

Copying lists. Here we use New to copy a list, where the copy is in reverse order. Note that the source list is of tier 1 while the copy is of tier 0.

y 0 = n i l : 0 ; while ( x 1 = n i l ) { z 0 :=y 0 : 0 ; New( y 0 ) ; suc ( y 0 ):= z 0 : 0 ; x 1 := suc ( x 1 ) : 1 } : 1

The loop is stationary, because the updated occurrence of suc is of type 0 → 0.

Soundness of programs for feasibility: run-time analysis

We show next that every tightly-ramified program computes a PTime function over configurations. The proof is based on the following observations, which we articulate more precisely below. First, if we start with a configuration where no vertex is assigned to variables of different tiers, then all configurations obtained in the course of computation have that property. We are thus assured that vertices can be ramified unambiguously.

The tiering rules imply that a program P of tier 0 cannot have loops, and is therefore evaluated in |P | steps. At the same time, the value of a variable of tier 1 depends only on vertices of tier 1. This implies, as we shall see, that the number of iterations of a given loop must be bounded by the number of possible configurations that may be generated by its body. Our restriction to tightly-modifying ramification guarantees that the number is polynomial.

Non-interference

Lemma 3 (Confinement). Let (Γ, ∆) be an environment. If Γ, ∆ P : 0, then Γ (X) = 0 for every variable X assigned-to in P .

The proof is a straightforward structural induction. Note also that a program P of tier 0 cannot have a loop, and is thus evaluated within |P | steps.

We say that a vertex-tiering Γ is compatible with a store µ if Γ (X) = Γ (X ) implies µ(X) = µ(X ) for all X, X ∈ X. We say that Γ is an initial tiering if Γ (X) is 1 for X initial (i.e. X ∈ X 0 ), and 0 otherwise. Thus an initial tiering is always compatible with an initial store.

Lemma 4. Suppose that Γ , ∆ P : α and S,

The proof is straightforward by structural induction on P .

We show next that tiering, when compatible with the initial configuration, guarantees the non-interference of lower-tiered values in the run-time of highertiered programs. A similar effect of tiering, albeit simpler, was observed already in [START_REF] Leivant | Predicative recurrence and computational complexity I: Word recurrence and poly-time[END_REF]. This is also similar to the security-related properties considered in [START_REF] Volpano | A sound type system for secure flow analysis[END_REF]. Non-interference can also be rendered algebraically, as in [START_REF] Marion | A type system for complexity flow analysis[END_REF].

The (Γ , ∆)-collapse of a configuration (S, µ) is the configuration (S ∆ , µ Γ ), where µ Γ (X) = µ(X) if Γ (X) = 1, and µ Γ (X) is undefined otherwise; whereas S ∆ is the structure identical to S except that each f for which (1 → 1) ∈ ∆(f ) is interpreted as ∅. Thus (S ∆ , µ Γ ) disregards vertices that are not not reachable from some variable of tier 1 using edges of type (1 → 1).

The next lemma states that a program's output vertices in tier 1 do not depend on vertices in tier 0, nor on edges that do not have tier 1 → 1.

The proof is straightforward by structural induction on programs.

Polynomial bounds

(Soundness) Lemma 6. Assume that Γ , ∆ P : α, where P is tightly-modifying. There is a k > 0 such that for every graph-structure S and every store µ compatible with

In proving Lemma 6 we will use the following combinatorial observation. Let G be a digraph of out-degree 1. We say that a set of vertices C generates G if every vertex in G is reachable by a path starting at C. The following lemma provides a polynomial, albeit crude, upper bound on the number of digraphs with k generators.

Lemma 7. The number (up to isomorphism) of digraphs with n vertices, and a generator of size k, is n 2k 2 .

Proof. A connected digraph of out-degree 1 must consist of a loop of vertices, with incoming linear spikes. If there are just k generators, then there are at most k such spikes. There are at most k entry points on the loop to choose for these spikes, and each spike is of size n. So there are at most n k × n k non-isomorphic connected graphs with a generator of size k.

Also, with only k generating vertices we can have at most k connected components. In sum there are at most (n 2k ) k = n 2k 2 non-isomorphic graphs of size n with k generators.

Proof of Lemma 6. We proceed by structural induction on P . The only nontrivial observations are as follows. For program composition, we use the Compatibility Lemma 4.

The crucial case is, of course, where P is of the form while(B)Q. Say X 1 , . . . , X m are the vertex-variables in B. The tiering rules require that B, and therefore X 1 , . . . , X m , are all of tier 1.

If Q updates only edges that are not probed in P (including the guard B), then neither the execution of Q nor the evaluation of B is affected, that is all configurations in the computation have the same vertices of tier 1, with no change in edges that affect the execution of P , by Lemma 5. Thus the truth of B in each invocation of Q is determined by the combinations of values assigned to X 1 , . . . , X m , while the structural changes caused by Q do not affect the execution of Q in subsequent invocations. Since we assume that P terminates, it follows that the combinations of values for X 1 , . . . , X m must be all different. If n is the number of vertices of tier 1, then n |S|, and there are n m such combinations. By IH Q terminates in PTime, and therefore so does P .

Suppose Q does update edges that are probed in P . Since P is assumed tightlymodifying, all such updates are for the same f ∈ F. Let C be the set of initial values of variables occurring in P (including X 1 . . . X m and possibly others). Let U be the set of vertices reachable from C by some path of tier 1 (i.e. using edges labels by h ∈ F assigned 1 → 1 by ∆from C vertices). By Lemma 5, the execution of P , including all iterated invocations of B and Q, has only vertices in U as value of tier 1. Moreover, U is generated by C, whose size is fixed by the syntax of P . It follows, by Lemma 7, that the number of such configurations is polynomial in the size of U , which in turns in bounded by the size |S| of the vertex-universe of the structure S.

Completeness of tightly-ramifiable programs for feasibility

Proposition 1. Every polynomial-time function on graph-structures is computable by a terminating and tightly-ramifiable program.

Proof. Suppose f is a unary function on graph-structures, which is computed by a Turing machine M over alphabet Σ, modulo some canonical coding of graph structures by strings. Assume that M operates in time k • n k . We posit that M uses a read-only input tape, and a work-tape. We simulate M by a tightlyramifiable program Γ , ∆ P : 1 over graph-structures, as follows. The structures considered simulate each of the two tapes by a double-linked list of records, with the three fields left, right, and val, returning two pointers and an alphabet letter, respectively. We take our data constants to include each σ ∈ Σ. The input-tape is assigned tier 1, and the work-tape tier 0. Initially the work-tape is empty, and new cells are progressively created by using the command New at tier 0. The machine's yield-relation between configurations is simulated at tier 0 by nested conditionals. Finally, we include in our simulation a clock, consisting of k nested loops, as in the following nesting of two loops: Hofmann and Schöpp [START_REF] Hofmann | Pure pointer programs with iteration[END_REF] introduced pure pointer programs based on a uniform iterator forall, and related them to computation in logarithmic space. Our programs differ from these in supporting modification of the structure, in the guises of vertex creation and edge displacement and deletion. Moreover, our programs are based on a looping construct that treats individually each vertex of the structure. These differences notwithstanding, our characterization of PTime can be modified to a characterization of log-space, by restricting the syntax of ramifiableprograms. Say that a program P over graph-structures is a jumping-program if it uses no edge update. Jones [START_REF] Jones | Logspace and ptime characterized by programming languages[END_REF] showed that a simple cons-free imperative programming language While \Ro accepts precisely the languages decidable by Turing machines in logarithmic space. Since our jumping-programs can be rephrased in While \Ro they too accept only log-space languages (where input strings are represented as linked lists). Conversely, the store used by a jumping-program is essentially of size k • log(n) + log(Q), where n is the size of the vertex-universe, Q the size of the data-universe, and k the number of variables. Consequently, we have: Theorem 2. A language is accepted by a jumping-program iff it is decidable in Logspace.

Adding recursion

It is not hard to augment our programing language with recursion. Here is a procedure that recursively searches for a path from vertex v to w. 9 Proc s e a r c h ( v 1 ,w 1 ) )

{ i f ( v=w) 1 return t r u e : 1 ; v i s i t e d ( v ) = t r u e : 1 ; f o r a l l t 1 i n A d j L i s t ( v ) % L i s t o f a d j a c e n c y nodes o f t { i f ( v i s i t e d ( t ) 1 =f a l s e ) i f ( s e a r c h ( t , w) 1 =t r u e ) return t r u e : 1 ; } return f a l s e : 1 ; } A restricted form of recursion is linear recursion, where at most one recursive call is allowed in the definition of a recursive procedure. Moreover, we suppose that each function body is stationary or tightly-modifying. Theorem 3. On its domain of computation, a tightly-ramifiable program with linear recursive calls is computable in polynomial time.


Max M Bruciamacchie

Serge S Garcia

Anne Stenger

How timber harvesting and biodiversity are managed in uneven-aged forests

Introduction

Nonindustrial private forest (NIPF) landowners have been shown to be more multi-objective by nature than industrial landowners: they give more importance to standing timber and forestland for the amenity values they provide [START_REF] Newman | Production economics of private forestry: A comparison of industrial and non-industrial forest owners[END_REF]. Among analyses of forest landowner behaviour, the household production framework recognises the benefits associated with forest amenities, as first applied by [START_REF] Binkley | Timber Supply from Non-Industrial Forests: A Microeconometric Analysis of Landowner Behavior[END_REF]. These non-market services are jointly produced with timber and are a determinant in the landowner's utility function. NIPF landowners comprise close to 70% of land ownership in many U.S. states and significant land holdings throughout Europe [START_REF] Amacher | Econometric analysis of forest landowners: is there anything left to do?[END_REF]. In France, almost 75% of the total forestland is privately owned, and 96% of private landowners are nonindustrial. In this article, we investigate the joint production of timber and biodiversity for NIPF landowners using a micro-econometric household production model. Even though our model is situated within a standard framework where a non-marketed good is jointly produced with timber products, we consider here that biodiversity is not totally disconnected from market strategies. Biodiversity is measured by the diversity of tree species. This assumption is based on the theory of coevolution introduced by [START_REF] Ehrlich | Butterflies and Plants: A Study in Coevolution[END_REF]. Coevolution acts as an evolutionary engine and a vehicle for biological diversification. Thus, the diversity of trees or plants may not only tend to increase the diversity of insects and animals, but the converse may also be true. In our model, tree diversity is a determinant of consumer satisfaction and a joint product in the profit-maximisation problem. Tree diversity has an additional impact: it is closely related to some market aspects since the different species have different monetary values. The forest landowner can decide to favour one tree species over another, depending on its value on the market. Conversely, he can make the choice of species diversification to cope with the volatility of timber prices. We focus on a complete set of forest landowners' decisions in uneven-aged forests where landowners are assumed to value the tree diversity of their forests, as well as timber harvesting. Our economic model is based on the maximisation of their utility that depends on the revenues from harvesting and tree diversity with respect to technological and budgetary constraints. The global objective of the paper is to explain the links between some of the harvest strategies of forest owners, unit price variability and the observed diversity of trees.

More precisely, we analyse: (1) their demand for species diversity and their timber supply; and (2) the joint production of timber and species diversity. Timber supply and amenity demand functions are derived using first-order conditions of the maximisation problem for the landowner. The behaviour of the forest owner is also strongly dependent on the characteristics of the forest blocks in question. Moreover, his/her harvesting strategy should differ according to the tree species and its value (depending itself on the quality and the diameter of the trees). The issue of heterogeneity in this case is crucial and its omission may result in consequent biases in the estimation stage. The estimation of timber supply and diversity demand is made using a database on uneven-aged forests in France for which several economic and ecological variables are regularly collected. This database typically concerns several forest blocks within which different tree species cohabit. This makes it possible to consider the forest owner within a multi-product framework where each product corresponds to a particular tree species.

Methods

Biodiversity and the economic model

In the literature on NIPF landowners, recent models of timber supply have included non-monetary returns or amenities [START_REF] Binkley | Timber Supply from Non-Industrial Forests: A Microeconometric Analysis of Landowner Behavior[END_REF][START_REF] Hyberg | The behavior of nonindustrial private forest landowners[END_REF][START_REF] Lähde | Stand structure as the basis of diversity index[END_REF][START_REF] Pattanayak | Timber and amenities on nonindustrial private forest land[END_REF][START_REF] Pepper | Robust inferences from random clustered samples: An application using data from the panel study of income dynamics[END_REF]. The idea is to better understand the trade-off between timber harvesting and amenity benefits. In this study, we attempt to understand forest owners' decisions concerning timber harvesting and biodiversity. Indeed, different tree species have different monetary values, and the forest landowner has several alternatives: to favour one tree species over another, depending on its market value, or to diversify the tree species in order to cope with the volatility of timber prices. Our definition of biological diversity may appear to be restrictive due to the sole inclusion of trees (instead of global biodiversity). Nevertheless, tree diversity accounts for a large part of biodiversity: it is generally accepted that the mixture of species is the guarantee of a certain degree of diversity of other living communities (for invertebrates, see [START_REF] Greatorex-Davies | The influence of shade on butterflies in rides of coniferised lowland woods in southern england and implications for conservation management[END_REF], and for bats, see [START_REF] Mayle | A biological basis for bat conservation in british woodlands -a review[END_REF]). This is the principle of coevolution [START_REF] Ehrlich | Butterflies and Plants: A Study in Coevolution[END_REF]. The diversity of trees or plants may not only tend to increase the diversity of insects and animals, but the converse may also be true. 1 Even if the extrapolation of tree diversity to global biological diversity is still in debate, this makes it possible to take both biodiversity and strategies on the timber market into account with only one indicator. Furthermore, there is no consensus about the choice of the diversity indicator. This is why several measures were tested in our model. We, in fact, used two notions, richness and diversity, the latter being the Shannon diversity index computed as H = -∑ h p h ln p h , where h represent a species. Three diversity indices were calculated:

1. Tree richness, designated by RICH, is computed as the number of species in the forest compartment. This is the simplest and the most intuitive index used to measure biodiversity. However, this measure strongly depends on the area surveyed.

2. The Shannon diversity index on the basis of number, designated by SHANN, is computed from the number of stems (n h ) with p h = n h ∑ h n h . 3. The Shannon diversity index on the basis of volume, referred to as SHANV, is expressed in volume v h : p h = v h ∑ h v h . The Shannon diversity index based on number is often used by ecologists, but the Shannon diversity index based on volume is more effect for characterising the crown size of different species.

In our model, tree diversity is a determinant of consumer satisfaction and a joint product in the profit-maximisation problem. The landowner i is represented in the framework of the household production function by a utility function that depends on the total income and non-pecuniary attributes:

U i = U(I i , z i ), (1)
where I i represents the total income of the landowner i and z i is the forest biodiversity.

The forest landowner faces a budget constraint where the total income is the sum of timber production profit π and exogenous income E:

I i = π ij + E i . ( 2
)
The timber profit π ij depends on timber production y ij sold at the price p ij , where the subscript j designates the tree species. The profit function is the difference between the timber revenue and the multi-product cost function related to the production of the (marketable) timber output y ij and the tree diversity z i conditional on some exogenous variables x ij (including forest capital and ecological variables). It can be written as:

π ij = p ij × y ij -C(y ij , z i , x ij ). (3)
Timber production y ij and tree diversity z i are linked by the following transformation function:

T(y ij , z i , x ij ) = 0. ( 4
)
The forest landowner has to choose the level of decision variables (i.e., y, z and I) that maximizes the utility function (1) subject to constraints (2) and (3). This utility maximisation problem can be solved by substituting these constraints into the utility function. The resolution is done in two steps: the household first selects the optimal level of I and z and then chooses the level of production y. In order to obtain explicit solutions to this problem, we have imposed some simple functional forms on our model. We chose a Cobb-Douglas form for the utility and cost functions. With these particular functional forms and by deriving with respect to y, we obtain the timber supply function that depends on timber price p, non-timber product z and other variables x. Expressing the first-order condition in log-linear form, we find the following timber supply function:

ln y ij = α 0 + α 1 ln p ij + α 2 ln z i + α 3 ln x ij , (5)
where the unknown parameters α are to be estimated. Note that α 1 represents the price elasticity of supply. If α 1 is respectively <, = or > 1 then the supply is price inelastic, unit-elastic or price-elastic. α 2 measures the trade-off between tree harvesting and diversity in terms of elasticity. If α 2 is negative, there is a substitution effect, whereas a positive sign is synonymous with complementarity.

Entering the equation ( 5) in the utility function and deriving it with respect to z give us the diversity demand. Transforming it into log-linear form, we have:

ln z i = β 0 + β 1 ln p ij + β 2 ln x ij , (6)
where β are the unknown parameters of the demand function to be estimated. β 1 represents the elasticity of diversity demand with respect to timber price. If β 1 is respectively <, = or > 1 then the diversity is inelastic relative to the timber price, unit-elastic or price-elastic.

The econometric approach

A two-step estimation procedure is implemented by first estimating the diversity demand equation (at the forest level), followed by the timber supply equation in which the predicted value of diversity is entered as a regressor.

Harvest observations collected for different tree species in different forests lead to the use of methods specific to cluster sampling [START_REF] Wooldridge | Cluster-sample methods in applied econometrics[END_REF]. However, the diversity of tree species is observed at the forest compartment level and is therefore cluster-invariant. Supposing that all variables are exogenous, the tree diversity demand equation ( 6) is estimated by the Ordinary Least Squares (OLS) method.

Cluster specificity is taken into account in the estimation of the timber supply equation ( 5).

The units within each cluster (or forest) may be correlated, whereas independence across clusters is assumed. Specific methods applied to Fixed Effect (FE) and Random Effect (RE) models make it possible to control for unobserved forest heterogeneity while studying the effects of factors that vary across species and forests (e.g., price), and others specific to forests (e.g., tree species diversity). [START_REF] Moulton | Random group effects and the precision of regression estimates[END_REF] shows the consequences of inappropriately using OLS estimation in the presence of random group effects. In particular, he demonstrates that the OLS standard errors that are not adjusted in this case are biased.

Consider the following timber supply cluster-sample equation:

y ij = α + X ij β + Z i γ + u ij , i = 1, . . . , N, j = 1, . . . , J i , (7)
where i indexes the "cluster" (or forest), j indexes individual observations within the cluster (or tree species). There is a total number of N clusters. The number of species is not the same throughout the different forests i, so that J (i.e., the number of species in the case of balanced data) is indexed by i. The total number of observations is n = ∑ J i . Harvest in the forest i of the tree species j is designated by y ij . X ij is a (1 × K) vector of explanatory variables that vary with respect to i and j. Z i contains L explanatory variables that only depend on the cluster i. u ij is the error term. α is the constant, and β and γ are the parameter vectors associated with the X and Z to be estimated, respectively. We consider the following unbalanced one-way error component:2

u ij = μ i + ij , i = 1, . . . , N, j = 1, . . . , J i , ( 8
)
where μ i is the cluster specific effect, and ij represents the remaining unobservables. μ i and ij are assumed to be independent and respectively i.i.d. (0, σ 2 μ ) and (0, σ 2 ). In matrix form, the one-way cluster model can be written as:

y = αι n + Xβ + Zγ + u = Rδ + u, ( 9
)
where u = R μ μ + , with R = (ι n , X, Z) and ι n a vector of n ones. y and R are of dimensions n × 1 and n × (1 + K + L). δ = (α , β , γ ) is the vector of parameters to be estimated. Finally,

R μ = diag(ι J i ) with ι J i is a vector of ones of dimension J i .
Supposing that all variables are exogenous, the equation can first be estimated by pooled OLS from the unbalanced data. The OLS estimator is trivially given by δOLS = (R R) -1 R y. It is unbiased and consistent. However, according to the method proposed by [START_REF] Pepper | Robust inferences from random clustered samples: An application using data from the panel study of income dynamics[END_REF], we use an estimate of the asymptotic variance matrix that is robust to heteroscedasticity and within-cluster correlation of arbitrary forms:

Var( δOLS ) = (R R) -1 ∑ N i=1 R i ûi û i R i (R R) -1 , where ûi is the N × 1 vector of OLS residuals (Y i -δOLS R i ).
Other consistent methods exist (some of which are more efficient), which make it possible to take the presence of unobserved effects in the error term into account. Cluster samples and panel data sets (where i represents individuals and j time periods) can be treated with similar methods (FE and RE models). In our case, the database has the same structure as an unbalanced panel data set. This is why we based our estimation method on the work of [START_REF] Baltagi | Incomplete panels : A comparative study of alternative estimators for the unbalanced one-way error component regression model[END_REF]. We can first consider that μ i represents the unobserved heterogeneity related to the forest, and treat it as a constant parameter to be estimated for each cluster i. If the fixed effects are correlated with the explanatory variables, there is an endogeneity problem that implies a biased estimator of parameters α, β and γ. We can obtain a consistent estimator of β by removing these effects with a suitable transformation (within-group transformation). However, an important drawback is that the parameters (γ) associated with cluster-invariant variables cannot be identified. The within-group transformation matrix for the (unbalanced)

cluster-sample case is Q = diag(E J i ). E J i = I J i - ι J i ι J i J i , where I J i is an identity matrix of dimension J i . The Within (or FE) estimator of β is: βFE = (X QX) -1 X Qy, ( 10
)
under the assumption of non correlation between and X. A drawback of this method is that γ cannot be identified because the variables Z disappear after within transformation. In order to take any possible autocorrelation or heteroscedasticity into account, [START_REF] Arellano | Computing robust standard errors for within-groups estimators[END_REF] proposes the following variance-matrix estimator:

Var( βFE ) = (X QX) -1 N ∑ i=1 QX i e i e i QX i (X QX) -1 ,
with e i = Qy -QX βFE , which is fully robust.

If the specific effects are assumed to be non-correlated with the explanatory variables, then a random effects (Generalised Least Squares, GLS) estimation can be used. Even if OLS estimators provide consistent parameters, a heteroscedasticity-consistent variance matrix is necessary. The effect μ i is now treated as a (cluster-specific) error term and assumed to be i.i.d. (0, σ 2 μ ). In this model, we can identify all coefficients related to all variables (including those that are cluster-invariant). Hence, the matrix of explanatory variables is now R = (ι n , X, Z). The vector of parameters δ = (α , β , γ ) and the variance components (σ 2 μ , σ 2 ) are estimated. The variance-covariance matrix of error terms u is Ω ≡ E(uu ) = σ 2 Σ, where Σ = I n + ρZ μ Z μ , with I n an identity matrix of dimension n and ρ = σ 2 μ σ 2 . The GLS (or RE) estimator is:

δRE = (R Ω -1 R) -1 (R Ω -1 y). ( 11
)
The variance of the RE estimator is:

Var( δRE ) = σ 2 (R Ω -1 R) -1 .
Several methods of estimation of variance components (σ 2 μ , σ 2 ) exist. However, the solution the most often chosen is the method of [START_REF] Swamy | The exact finite sample properties of the estimators of coefficients in the error components regression models[END_REF] by using the Within and Between residuals. The RE estimator is asymptotically more efficient than pooled OLS under the usual RE assumptions. However, if the cluster effects are correlated with μ i are correlated with X or Z, this estimator is not consistent. This possible endogeneity can be tested for by performing a Hausman test. The Hausman test statistic is:

( βFE -βRE ) [Var( βFE ) -Var( βRE )] -1 ( βFE - βRE ).
Under the null hypothesis, this statistic has an asymptotic chi-square distribution with a number of degrees of freedom equal to the number of cluster-variant variables (K).

Results and discussion

Data sources and characteristics

The database of the AFI network (Association Futaie Irrégulière -Uneven-aged forest network) was used. Uneven-aged forest management is characterised by two fundamental principles: the use of natural dynamics of the ecosystem and the individual treatment of each tree. The first principle implies the use of all tree species on the site: forests are always mixed-species (with variations depending on the acidity of the soils). The second principle means that each tree is examined in order to assess its different functions (e.g., value-added wood, aesthetic aspect). Hence, the decision of tree harvesting or conservation does not result from the stand age but rather from its functionality: Does this tree "pay" for its place? [START_REF] Bruciamacchie | La Futaie Irrégulière: Théorie et Pratique de la Sylviculture Irrégulière[END_REF]. Uneven-aged forest management is practised in numerous forests worldwide with a multitude of variations in terms of species composition and stand structures under local ecological, social and economic constraints. The AFI network consists of 68 compartments in the northern part of France. The compartment is the management unit for uneven-aged forests (whereas the whole forest is the unit considered for even-aged forests) and corresponds to a block that varies from 5 to 15 ha. One compartment is made up of ten permanent plots that make it possible to monitor the individual growth of approximately 200 trees per compartment. These compartments also make it possible to monitor poles, coppice and regeneration. Some of them are good examples of successful transitions between even-aged and uneven-aged stands. Our sample is made up of forests whose stands are well-balanced in terms of forestry (consistent harvesting), which makes it possible to handle economic data that are uniform on the long term. As mentioned above, we consider a forest owner who maximises his utility that is a function of total income and diversity. The forest owner decides on the main orientations of his/her forest management (e.g., level of revenues, distribution of species, risks concerning species management). However, we wanted to introduce an important characteristic of forest management into the empirical model: in practice, forests are managed by the "owner/forest manager" pair. Indeed, the owner often delegates the management to a forest manager who implements the owner's choices and can thus have an influence on the harvesting decision and the distribution of species. This is why we include dummies that proxy the identity of the manager (see below). Among the 68 compartments, 39 were selected because all of the information in all of the categories of variables was available. We classified tree species into five classes: oak, beech, precious broad-leaved trees, other broad-leaved trees and conifers. five classes of species are not observed in all of the compartments, so that the total number of observations in our sample is 102. 3However, the number of species is greater and we compute the diversity for each compartment from the total number of species (varying from 2 to 14 in our sample). As presented above, we calculate three diversity indices. The first index used is tree richness, designated by RICH, simply computed as the number of species in the forest compartment. The last two are Shannon diversity indices computed as H = -∑ h p h ln p h , where h represent a species. 4 We compute a Shannon diversity index on the basis of number (SHANN) and a Shannon diversity index on the basis of volume (SHANV), already defined above. The variables used in the model are the following:

• Variables observed per compartment and broken down by species: harvested volume (y), unit price (p),5 stock inventory (INV), volume increment (VOLINCR).

• Percentage of quality (QU AL%) and average diameter (DI AM) are measured for standing timber.

• At the compartment level, seven dummy variables (from ST1 to ST7)) are built for seven different ecological conditions ranging from the more basic to the more acid soils. In fact, this set of dummies represents an ecological indicator built from the variables, pH and moisture.6

• The type of owners is represented by four dummy variables: institution (DUMO1), individual (DUMO2), group of owners (DUMO3) or joint ownership (DUMO4).

• The owner often delegates the management to a forest manager. He/she implements the owner's decisions but can have an influence on the distribution of species. Dummies DUME1 to DUME10 are used for the manager. DUME10 is the remaining sum of managers that are in charge of only one forest compartment.

Descriptive statistics are reported in Table 1.

Estimation results

We first estimate the tree diversity demand equation ( 6). The diversity of tree species is observed for each forest compartment and is cluster-invariant. Since some explanatory variables vary according to the forest compartment as well as to the species (such as the price), the diversity equation is estimated by a (between-type) OLS method. All variables can be considered as exogenous in this estimation (at least, on the short term). In particular, the price is determined by the market. Hence, there can be no doubt about the direction of the cause-effect relationship. For example, it is the timber price that explains the tree diversity in a compartment and not vice versa. As mentioned above, there are three different indices to proxy diversity. Three regressions were successfully run with the three different indices as dependent variables. The estimated coefficients are similar. However, the goodness of fit as well as the significance of parameters are better with the logarithm of the number of species (i.e., the richness index). The richness varies as soon as an individual of a new species is added or removed. Shannon indices are preferred by ecologists because they take the richness as well as the distribution of species into account at the same time. However, according to the managers, taking biodiversity into account tends to favour minority species. Estimation results are presented in The overall performance of the demand equation is good since the adjusted R 2 is equal to 0.602. The estimated parameters are all significantly different from zero, except for the stock inventory (i.e., the standing timber per ha) and a dummy variable that proxies a forest expert. However, other variables related to the state or trend of forest capital such as the average diameter of trees, the share of qualitative stand wood and the volume increment of forest are significant in our model. In particular, the negative sign for the coefficient associated with the percentage of quality (QU AL%) has an interesting interpretation. Forests with the highest percentage of quality correspond to ones with the lowest diversity. Since a high percentage of quality increases the revenues over time, this result would mean that in this case, species diversity is less favoured, showing a trade-off between quality and diversity. Moreover, the coefficient associated with the variable DI AM is significantly positive. In order to favour diversity, some trees were harvested early to diminish natural competition between species.

As expected, the site context has a significant impact on the diversity. Coefficients associated with dummies from ST1 to ST6 are all significant with positive signs (decreasing from 2.43 to 1.88, respectively) with respect to acid soils (ST7), confirming a decrease in richness when the context is acid. Furthermore, the estimated coefficients allow a classification of the site conditions that is in agreement with the observed ecological link between the chemical characteristics of the soil and tree (and flora) diversity. Some forest managers have a significant positive impact on tree diversity, while other ones have a negative impact that supports a short-term view. The variables for the type of forest owner have been removed because their coefficients were not significantly different from zero.

The unit price has a significant and positive influence on the diversity. Its value (0.3110) means that a 10% decrease in timber price implies a 3.11% decrease in tree diversity. This result highlights the effect of timber price on the abandonment of species. For example, in the ecological context where diversity is the highest (14 species in our sample), a 23% decrease in price could lead to the loss of one species. Unit prices for timber are exogenous. However, average unit price (for one species in one compartment) can vary according to the distribution in the stand with respect to its quality and its size. The forest owner can therefore adjust his/her revenues by acting on these variables. One of the principles in uneven-aged forest management is to concentrate the volume increment on the high-quality trees. Hence, low-quality trees are progressively cut and, at the same time, the unit price of standing timber as well as that of harvested timber increase. Once this unit price has increased, forest managers and owners are more inclined to maintain the minority species. The objective is to reduce economic risks by finding an optimal distribution among the different species.

Using the estimates of the demand equation, the fitted value of diversity was computed and used as an explanatory variable in the timber supply equation ( 5). The use of generated regressors may produce non-consistent estimated standard errors. This is why a vector of regressors was used that includes some or all exogenous variables already in the first regression [START_REF] Pagan | Econometric issues in the analysis of regressions with generated regressors[END_REF]. This second-step OLS leads, in fact, to a two-stage least squares procedure since the regressors are variables used in the first-step estimation (of the demand equation), and gives correct standard errors. Because the predicted diversity ln RICH can be approximated by a linear function of the explanatory variables in the demand equation and leads to a problem of collinearity, several exclusion restrictions were used in the supply equation. Some variables that do not appear to be significant to explain harvesting have thus been excluded, including the volume increment of stock (VOLINCR) and some dummies that proxy the forest manager. Finally, this estimation procedure is implemented with a robust variance-covariance matrix. Within and GLS methods (for FE and RE models, respectively) are implemented as described in the econometric method section. They are also conducted in two steps like the OLS method. A Hausman test was then computed to check for the exogeneity of explanatory variables. The value of the statistic is 3.217 (with a P-value of 0.5222) and is below the χ 2 (4) critical value at the 1% level. This result confirms the exogeneity of variables. Hence, the GLS method is the best adapted here for dealing with the cluster feature of our sample. R 2 is equal to 0.606 and indicates a good fitting of our model. Estimation results of (second-step) OLS, Within and GLS methods are reported in [START_REF] Pepper | Robust inferences from random clustered samples: An application using data from the panel study of income dynamics[END_REF] and [START_REF] Arellano | Computing robust standard errors for within-groups estimators[END_REF].

Table 3. Supply estimation -Cluster-sample econometric methods

As explained above, OLS is less efficient than GLS since it does not fully take the cluster feature of our sample into account, even if a robust variance-covariance matrix makes it possible to alleviate this problem. OLS coefficients are rather similar to those estimated by specific cluster methods. However, some interest coefficients such as those associated with price and diversity are slightly overestimated. For example, the coefficient associated with the price is 0.57 with OLS, compared to 0.43-0.47 with GLS. For the diversity, it is equal to 0.87, compared to 0.75-0.78 (in absolute value).

The coefficients associated with the variables QU AL% and DI AM are significantly negative (with estimates of -1.38 and -0.04, respectively). This means that high-quality trees with big diameters are harvested to be sold. Hence, the actual standing timber is characterised by a lower percentage of quality and a lower average diameter. Moreover, it is not surprising to see that whereas forest managers have a positive impact on tree diversity, this is not the case for timber harvest.

Results also show a positive and significant impact of both timber inventory and unit price. As expected, timber harvest increases with the standing volume of trees. The coefficient associated with the price (or price elasticity of timber supply) is estimated at 0.47, meaning that a 10% increase in price implies a 4.7% increase in harvesting.

The diversity is negatively and significantly correlated to the timber harvest, all things being equal. The estimated coefficient can be directly interpreted as a measure of substitution between tree diversity and the volume of timber harvested. The point estimate is equal to -0.78. This value is rather high. However, based on the standard error estimate, we can reject the hypothesis of a unitary elasticity substitution. An explanation for this negative sign is that when the site context is acid, the forest manager cannot influence the unit timber price interval per species. In this case, under acid soil conditions, the forest manager can only act on timber volume. On the contrary, the basic context allows for a greater variety of species. However, in order to favour all species, the forest manager cannot increase the standing volume and in some cases, may be forced to reduce it. Hence, the forest stock is low on the long term and this trend leads to a lower timber harvest.

Conclusion

In this study, a household production approach was used to model the behaviour of the NIPF owner in order to derive the structural econometric equations of timber supply and diversity demand and to estimate substitution and price elasticities. In the empirical application, a definition of diversity was chosen solely on the basis of the number of tree species. This diversity is simple to calculate and positively correlated with the diversity in flora and fauna. Moreover, the richness of data related to harvested species and the cluster-sample methods used in this context make it possible to deal with heterogeneity and variability within clusters.

In addition, Within and GLS estimation methods make it possible to test for the possible endogeneity problem of some variables. This study revealed that diversity demand and timber supply are negatively linked, meaning that an increase in tree diversity will lead to a decrease in timber harvesting. This result confirms that these two forest outputs are substitutes. Estimation also shows that timber price and tree diversity evolve in the same direction: the positive and significant coefficient associated with the timber price in the demand equation indicates that a price decrease has a negative effect on diversity. This result is certainly the consequence of the characteristics of uneven-aged forests and the strategies used to manage them. This could be explained by the fact that a part of the diversity not only procures some satisfaction for the forest owner, but that the price paid for this diversity is a decrease in timber production. Management strategies should therefore be aimed at finding a trade-off between timber production and tree diversity in a given ecological context.

Table 2 .

 2

	Variable	Coefficient s.e	Variable Coefficient	s.e.
	Dependent variable: ln z = ln RICH
	Constant -4.1585** 1.7642	DUME6 -0.3070***	0.0837
	ln p	0.3110*	0.1781	DUME8 1.2389***	0.2417
	ln DI AM 0.7357** 0.2985	ST1	2.4252***	0.4362
	QU AL% -2.4432*** 0.6085	ST2	2.3281***	0.3816
	ln INVD	0.0645	0.2351	ST3	2.2615***	0.3882
	VOLINCR 0.4216*** 0.1398	ST4	1.8120***	0.4080
	DUME2	0.3974*** 0.0970	ST5	1.7468***	0.4016
	DUME3	0.1129	0.1259	ST6	1.8774***	0.4155

Notes: n = 102, N = 39. Adjusted R 2 = 0.602. Heteroscedasticity-consistent s.e. ***: significant at 1%, **: significant at 5%, *: significant at 10%.

Table 2 .

 2 Demand estimation -OLS method

Table 3

 3

			.
		OLS		Within		GLS
		(Pooled)	(FE -Fixed Effects)	(RE -Random Effects)
	Variable	Coef.	Robust s.e. Coef. Robust s.e. Coef.	s.e.
	Dependent variable: ln y
	Constant -1.2893*	0.7816 -3.3204*** 0.2433 -1.3713**	0.5866
	ln p	0.5735***	0.0955	0.3467*** 0.1089	0.4659***	0.0699
	ln DI AM -0.0368***	0.0089 -0.0424*** 0.0072 -0.0387***	0.0066
	QU AL% -1.4822**	0.6484 -1.1813	0.7142 -1.3775***	0.4515
	ln INV	0.7178***	0.0771	0.9054*** 0.0756	0.8085***	0.0630
	ln RICH -0.8651**	0.3456	-		-0.7797***	0.2610
	DUME2 -0.8210*	0.4226	-		-1.0129***	0.2387
	DUME3 -0.3628	0.2651	-		-0.4820**	0.2119
	ST4	-0.5238**	0.2433	-		-0.5008**	0.1980
	ST5	-0.8125**	0.3761	-		-0.7456**	0.3018
	σ2					0.5065
	σ2 μ					0.2861
	Hausman test (P-value)			3.217 (0.5222)
	R 2	0.613		0.460		0.606

Notes : n=102, N=39. ***: significant at 1%, **: significant at 5%, *: significant at 10%. Robust s.e. for OLS and FE estimation are respectively computed following

Many references exist on this topic, see[START_REF] Lähde | Stand structure as the basis of diversity index[END_REF],[START_REF] Barbier | Influence of tree species on understory vegetation diversity and mechanisms involved -a critical review for temperate and boreal forests[END_REF],[START_REF] Schuldt | Communities of ground-living spiders in deciduous forests: Does tree species diversity matter?[END_REF],[START_REF] Mcdermott | Short-and long-term implications of clearcut and two-age silviculture for conservation of breeding forest birds in the central Appalachians, USA[END_REF], among others.

Sustainable Forest Management -Current ResearchHow Timber Harvesting and Biodiversity are Managed in Uneven-Aged Forests: A Cluster-Sample Econometric Approach

Only five species are observed and not within all forests. We can therefore not implement a two-way error component regression model. Moreover, each forest is observed only once since we only have cross-section data.

In a complete data cluster, the number of observations would be 195.

We use two different subscripts in our article. Subscript j refers to the (five) classes of species, whereas h refers to the species alone (the total number of species varies from 2 to 14 in our sample).

Unit price refers to the market price depending on species, diameter and quality. In the empirical model, we use the average unit price, i.e., the unit price for one species in one compartment.

In reality, a more in-depth ecological study would take pH, moisture and altitude into account. There is actually no significant variation in altitude since all forests observed in our sample are located at altitudes below 500 meters.

Sustainable Forest Management -Current Research


Balázs Németh
email: bnemeth@sztaki.hu

Péter Gáspár
email: gaspar@sztaki.hu

József Bokor
email: bokor@sztaki.hu

Olivier Sename
email: olivier.sename@gipsa-lab.grenoble-inp.fr

Luc Dugard
email: luc.dugard@gipsa-lab.grenoble-inp.fr

Fault-tolerant control design for trajectory tracking in driver assistance systems

Keywords: fault-tolerant control, reconfiguration, fault detection, linear parameter varying control, robust control, autonomous systems

The paper proposes a control system with the brake and the steering for developing a driver assistance system. The purpose is to design a cruise control method to track the road geometry with a predefined velocity and guarantee the road stability of the vehicle simultaneously. An actuator selection method is developed in the control design, in which the actuator limits, energy requirements and vehicle operations are taken into consideration. The method is extended with a fault-tolerant feature based on a robust LPV method, into which the actuator selection procedure and the detected fault information are incorporated. The operation of the reconfigurable control system is illustrated through various vehicle manoeuvres.

INTRODUCTION AND MOTIVATION

The purpose of trajectory tracking is to follow a road geometry with a velocity defined by the driver and guarantee the road stability of the vehicle simultaneously. Since the actuators affect the same dynamics of the vehicle, in the operation of control systems interference or conflicts may occur between the control components. In the control design the interaction between the actuators must be taken into consideration and a coordination between them must be achieved. An integrated control system is designed in such a way that the effects of a control system on other vehicle functions are taken into consideration in the design process.

The demand for vehicle control methodologies including several control components arises at several research centers and automotive suppliers. Recently, important survey papers have also been presented in this topic, see e.g. [START_REF] Yu | Integrated vehicle dynamics control: State-of-the art review[END_REF]. A vehicle control with four-wheel-distributed steering and four-wheel-distributed traction/braking systems was proposed by [START_REF] Ono | Vehicle dynamics integrated control for four-wheel-distributed steering and four-wheel-distributed traction/braking systems[END_REF]. A yaw stability control system in which an active torque distribution and differential braking systems were used was proposed by [START_REF] Zhang | Vehicle stability control strategy based on active torque distribution and differential braking[END_REF]. Differential braking and front steering to enhance the vehicle yaw stability and the lateral vehicle dynamics was proposed by [START_REF] Doumiati | Gain-scheduled LPV/Hinf controller based on direct yaw moment and active steering for vehicle handling improvements[END_REF]. An integrated control that involves both four-wheel steering and yaw moment control was proposed by [START_REF] Jianyong | Improvement of vehicle handling and stability by integrated control of four wheel steering and direct yaw moment[END_REF]. Active steering and suspension controllers were also integrated to improve yaw and roll stability [START_REF] Poussot-Vassal | Attitude and handling improvements through gain-scheduled suspensions and brakes control[END_REF]. A global chassis control involving an active suspension and ABS was proposed by [START_REF] Gáspár | Design of a reconfigurable global chassis control[END_REF]; [START_REF] Zin | An LPV/Hinf active suspension control for global chassis technology: Design and performance analysis[END_REF]. The driveline system and the brake were integrated in [START_REF] Rajamani | Demonstration of integrated longitudinal and lateral control for the operation of automated vehicles in platoons[END_REF]. A possible integration of the brake, steering and suspension system was presented by [START_REF] Trachtler | Integrated vehicle dynamics control using active brake, steering and suspension systems[END_REF].

The paper proposes a control system with two active components for developing a driver assistance system. The purpose of the control is to generate control inputs, such as the steering angle and the difference in brake forces. Since both the actuators affect the lateral dynamics of the vehicle, in the control design a balance and priority between them must be achieved. An actuator selection method is applied in the control design. Moreover, detected fault information is also considered in order to guarantee the reconfigurable and fault-tolerant operation of the vehicle.

The paper is organized as follows: in Section 2 the vehicle model and the longitudinal-lateral trajectory tracking are formalized. In Section 3 the closed-loop interconnection structure is formalized and an actuator selection method is applied. In Section 4 the architecture of the control system and the fault-tolerant control are presented. In Section 5 simulation results are presented.

VEHICLE MODEL FOR THE TRAJECTORY TRACKING

In the design of trajectory-tracking assistance systems it is necessary to guarantee that the vehicle must perform the desired motion of the driver. The control system of the lateral vehicle dynamics assists the driver in tracking road geometry. It has advantages in critical situations, in which the driver is not able to ensure vehicle stability. In trajectory tracking the vehicle is moving in the entire plane of the road, thus both the longitudinal and the lateral dynamics must be taken into consideration as Figure 1 shows.

α 1 + δ α 2 β ξ l 1 l 2 X g l Y g l X v Y v ψ y v y g l M b r F l

Fig. 1. Lateral dynamic model of vehicle

Two actuators are used in the system, i.e., the front-wheel steering angle δ and the differential brake torque M br . In most of the lateral control problems, the lateral dynamics of the vehicle can be approximated by the linear bicycle model of the vehicle:

J ψ = C 1 l 1 α f -C 2 l 2 α r + M br (1a) mv( ψ + β) = C 1 α f + C 2 α r (1b)
where m is the mass, J is the yaw-inertia of the vehicle, l 1 and l 2 are geometric parameters, C 1 , C 2 are cornering stiffnesses, ψ is the yaw rate of the vehicle, β is the sideslip angle. Moreover, α f = -β +δ -l 1 • ψ/v and α r = -β + l 2 • ψ/v are the tyre side slip angles at the front and rear, respectively.

Two control systems will be designed based on the state space representation of the vehicle: ẋ = Ax + Bu (2) where the state vector consists of the yaw-rate and the side-slip angle of the vehicle x = ψ β T . In the brake control case the input of the system is u = M br , while in the steering control case the input is u = δ. The measured output of both systems is the yaw-rate, y = ψ. This approach is suitable in the decentralized control concept, where the components are designed independently. The advantage of this solution is that the components with their sensors and actuators can be designed by the suppliers independently. Since the controllers guarantee performances only locally, the stability and performances of the entire closed-loop system must also be guaranteed. It is required to perform an analysis step in the robust control framework on a global level.

CONTROL DESIGN BASED ON WEIGHTING FUNCTIONS

Performance specifications

In the driver assistance system the performance is the minimization of the tracking error of the yaw-rate

z 1 = [ ψref -ψ] T → min! (3)
where ψref is the reference yaw rate defined by the driver. The reference yaw-rate of the controller can be calculated from the steering wheel angle, see [START_REF] Pacejka | Tyre and vehicle dynamics[END_REF].

Simultaneously, actuator saturations must be avoided. The maximum control input of the steering is determined by its physical construction limits, while in the case of the braking system the constraints are determined by the tyreroad adhesion. These constraints will be built into the weighting strategy applied in the control design. The other performance of the system in terms of the control input is formalized as

z 2 = |u| → min! (4)
The control design is based on a weighting strategy, which is formalized through a closed-loop interconnection structure, see Figure 2. In the trajectory tracking problem the yaw-rate reference signal is introduced in order to guarantee the tracking of the road geometry: R = ψref . The role of the weighting functions is to define performance specifications, reflect disturbances and uncertainty.

Since the coordination between the actuators and creating priority between them are in the focus of the paper, in the following the design of the weighting functions for actuators is presented.

G K (ρ ) W p W n w n F d W w Δ z 1 Δ P K u y R W a c t z 2 ρ W u Fig. 2. Closed-loop interconnection structure

Weighting function for the actuators

In this section an actuator selection method is developed in the control design, in which the actuator limits, energy requirements and vehicle operations are taken into consideration. Since the steering angle and the brake moment actuators affect the same dynamics of the vehicle, a balance between them must be achieved.

First the steering operation is analyzed. Steering has a construction limit, i.e., the value of front-wheel steering can not exceed an upper bound. In order to avoid a steering limit differential braking must be increased. During driving the steering angle is used to handle the lateral dynamics. Moreover, during close to the limit of skidding steering is also preferred to the brake.

Second the brake intervention is considered. The brake moment is limited by the adhesion value between the road and the tire. It is necessary to prevent the skidding of tires, thus in case of skidding the differential braking must be reduced, while the yaw-motion of vehicle must be controlled by front-wheel steering. By using differential braking the velocity of the vehicle is decreased. Thus, during driving the use of differential braking must be avoided and frontwheel steering is preferred. During deceleration, however, the brake is already being used, thus the lateral dynamics is handled by the braking for practical reasons.

Two weighting factors ρ st , ρ br are introduced in order to take into consideration the influence of the steering and the brake moment. These are built into the weighting functions applied to the control design. The weighting for the front wheel steering and that for the brake yaw-moment are W act,st = ρ st /δ max (5a) W act,M br = ρ br /M brmax (5b) respectively, where δ max is determined by the constructional maximum of the steering system, while M brmax is the maximum of the brake yaw-moment. Figure 3 shows the characteristics of the weighting factors. In the case of driving the front wheel steering is actuated, which is determined by factor ρ st , see Figure 3(a). The value is reduced between δ 1 and δ 2 , which represents the constructional criterion of the steering system. In the case of braking the tyre longitudinal slip angle affects factor ρ br , see Figure 3(b). In this interval differential braking is preferred for practical reasons. Reducing tyre skidding requires an interval. Therefore two parameters are designed: α 1 and α 2 are applied to prevent the skidding of tyres. An interval to prevent chattering between steering and differential braking: F l,1 and F l,2 is also required.

|δ| δ 2 δ 1 ρ s t 1 0 F l F l , 1 F l , 2 (a) Steering α α2 α1 ρb r 1 0 Fl Fl , 1 Fl , 2 (b) Brake yaw-moment
In the following it is assumed that the longitudinal slip and the longitudinal force are measured or estimated, these weighting factors are available during the journey. The model, which is the basis of the control design, is an LPV form and the control design is based on the LPV method. The purpose of the quadratic LPV design method is to choose the parameter-varying controller K(ρ) in such a way that the resulting closed-loop system is quadratically stable and the induced L 2 norm from the disturbance and the performances is less than the value γ. Stability and performance are guaranteed by the design procedure, see [START_REF] Bokor | Linear parameter varying systems: A geometric theory and applications[END_REF]; [START_REF] Packard | Theory and application of linear parameter varying control techniques[END_REF]; [START_REF] Wu | Induced l 2norm control for LPV systems with bounded parameter variation rates[END_REF].

DESIGN OF THE FAULT-TOLERANT CONTROL SYSTEM

4.1 Architecture of the control system

The purpose of control design is to calculate the necessary front steering angle and brake yaw moment. The design of this upper level controller is based on the LPV method.

The designed longitudinal force and brake yaw moment are distributed between the four wheels of the vehicle.

Moreover, a third layer is also necessary since the required control forces must be tracked by using a low-level controller. This controller transforms the wheel forces and the values of the steering angle into a real physical parameter of the actuator. These components are implemented by Electronic Control Units (ECUs).

The design of a low-level steering controller might use more specific techniques that fit the specific nonlinear properties of the actuator. The steer-by-wire front steering system transforms the steering angle into a real physical parameter of the actuator. The real physical input of the system is the Pulse Width Modulated (PWM) signal of the electric servo motor, which moves the rack. The physical construction of electric steering has several variations, see e.g. [START_REF] Claeys | Modeling and control of steering actuator for heavy duty vehicle[END_REF]. Figure 4 shows the architecture of the low-level steering controller. The architecture of the controlled supervisory system is shown in Figure 5. In the fault-tolerant scheme fault detection and isolation (FDI) filters for actuators are assumed to be used. In this paper two kinds of actuator faults are considered: the fault of the steering control system and the fault of the braking circuits. There may be various fault scenarios, e.g the leakage of the hydraulic systems in the braking or steering servo, or the steering mechanism becomes jammed. The different changes in the operation of an actuator make it possible to realize the detection of a fault.

The filters are able to detect different types of faults in the operation of the actuators. An H ∞ method to design a fault detection and isolation (FDI) LPV filter was presented by [START_REF] Edelmayer | Robust detection filter design in the presence of time-varying system perturbations[END_REF]. The geometric approaches often lead to successful detection filter design, for details see [START_REF] Bokor | Detection filter design for LPV systems -a geometric approach[END_REF]. The selection of the performance weights in the design of FDI filters has been applied to vehicle systems, see in [START_REF] Gáspár | Lpv design of faulttolerant control for road vehicles[END_REF]. This paper focuses only on the design of fault-tolerant control and it is assumed that the FDI filters have been designed and they are available.

Modification of the weighting functions

Two actuators are operated in cooperation in order to provide a reconfigurable fault-tolerant control system. In case of a detected fault either the brake yaw moment M br or the front wheel steering δ can be changed with similar dynamic effects.

When a fault occurs in the operation of the steering system, all the lateral control tasks must be realized by using the braking system with the generation of the brake yaw moment M br . If fatal error occurs in the operation of the steering system the weight of steering is masked:

ρ st = 0.
When a fault occurs in the operation of a brake circuit the actuated brake yaw-moment is reduced. Moreover, the reduction of the brake yaw-moment is asymmetric. For example, in the case of a fault of a brake circuit on the left-hand side of the vehicle, the generated positive brake yaw-moment is reduced, or it is zero. In this case steering is activated to substitute for the actuation of braking and provide trajectory tracking. However, the negative M br can be realized by the healthy right-hand-side brake circuits. Consequently, the weight of braking ρ br depends on the sign of the desired M br . In the case of a lefthand-side brake circuit fault, positive M br is not allowed, therefore ρ br = 0. However, if M br < 0 then ρ br > 0. The actual modification of ρ br is based on a design parameter: ρ br,new = κ i • ρ br m where κ i is a weighting factor.

Quadratic stability of the entire control system

The stability of the individual controllers is guaranteed by the design method. The global control system contains two controllers, the brake and the steering. The global system uses two scheduling variables ρ = [ρ br , ρ st ]. According to Figure 3, these factors have limits. When these controllers are used simultaneously it is necessary to guarantee the stability of the global closed-loop system.

A common Lyapunov function for the closed-loop systems must exist. The following affine parameter-dependent closed-loop system is given, see [START_REF] Scherer | Lecture Notes DISC Course on Linear Matrix Inequalities in Control[END_REF]; [START_REF] Boyd | Linear Matrix Inequalities in System and Control Theory[END_REF]: ẋ(t) = A(ρ) x(t).

(6) where A(ρ) = A 0 + ρ 1 A 1 + ...ρ 4 A 4 . For the stability of the system (6) it is necessary to guarantee that all trajectories of system A converge to zero as t → ∞. A sufficient condition for this is the existence of a quadratic function V (ξ) = ξ T P ξ, P > 0, which decreases along every nonzero trajectory of (6). If there exists such a P , then ( 6) is said to be quadratically stable and V is called a quadratic Lyapunov function. The necessary and sufficient condition for quadratic stability of system (6) for all of A i is A T cl,i P + P A cl,i < 0; P > 0; i = 1, . . . n (7) Therefore it is necessary to find a V common Lyapunov function for all of the closed-loop systems which can guarantee the global stability of the systems in every scheduling variable.

The matrices of the closed-loop system A cl,i are computed using the next formula:

A cl,i = A + B 2 D ci C 2i B 2 C ci B ci C 2 A ci (8)
where A, B 2 , C 2 are the state space representation of the plant, A ci , B ci , C ci , D ci are the state-space representations of the controllers. The aim is to find a solution to P > 0. To analyze the global stability of the LTI systems,

ρ br ρ st 1 1 0 vertex of convex hull

Fig. 6. Convex hull of LTI systems

Co{A 1 , . . . A 4 } is covered by the convex hull of finitely many matrices A cl,i . According to the system, the convex hull contains 4 LTI systems, see Figure 6. For the analysis of global stability this convex hull can be used.

SIMULATION RESULTS

In this section the fault tolerance of the control system is illustrated through simulation examples. Several software packages are used for the design and analysis of the controlled system. The control design is performed by using the Matlab/Simulink software. The verification of the designed controller is performed by using the CarSim software. In this package the model of the actual road vehicle dynamics is represented with high accuracy.

The vehicle is traveling along a predefined road, while the integrated control system supports the driver to guarantee trajectory tracking. During the simulations different faults occur and these faulty cases are compared with a healthy simulation. A typical E-Class automobile is applied in the simulation. The mass of the 6-gear car is 2023 kg, its engine power is 300 kW (402 hp). The width of the track is 1605 mm and the wheel-base is 3165 mm. In the simulation examples the vehicle is traveling along a section of Waterford Michigan Race Track, which is shown in Figure 7(a). The velocity of the vehicle changes along its route as Figure 7(b) shows.

In the first simulation a steering fault occurs in the controlled system. Note that the driver assistance system is not able to modify front wheel steering angle, but the driver can steer the front wheels. The control system actuates only brake yaw-moment M br . Figure 8 shows the faulty simulation case compared with a healthy one. The lateral error of the system and the yaw-rate tracking are illustrated in Figure 8(a) and Figure 8(b). The integrated control system can tolerate a steering fault, the lateral error and the yaw-rate of faulty simulation results are close to the healthy cases. The largest difference is at the last bend. In this bend the longitudinal slip of the wheels reach -1, while in the fault-free case it can be reduced using the actuation of the front wheel steering. The reason for the skidding is the increased brake pressures, compared to the fault-free case, see Figure 8 In the case of a steering fault the weight of steering ρ st is equal to zero, while the weight of braking is influenced by skidding. In the second simulation example one of the rear brake circuits fails. In Figure 9(a) the effect of brake faults is shown. In the first bend the vehicle turns right, which means that the rear right-hand-side wheel brake is actuated to perform the maneuver. Therefore rear righthand-side brake circuit fault increases the lateral error. In the case of bends to the left the fault of the rear lefthand-side wheel circuit increases the lateral error. Figure 9(b) illustrates the steering wheel angle, which is actuated by the driver. It can be seen that the fault in the brake system necessitates faster and more intensive intervention by the driver. A deterioration of the braking effect induces an increase in the front wheel steering to perform the maneuver, see Figure 9(c). If a fault occurs in the brake the actuated M br moment has a limitation, as shown in Section 4. In the case of a left-hand-side brake circuit fault the vehicle is not turned anti clockwise, therefore positive M br is not allowed and vice versa. The actuated brakeyaw moments can be seen in Figure 9(d). Figures 9(e) and 9(f) show the actuated brake pressures, which prove the limitation of the brake-yaw moment. In the third simulation example all of the rear brake circuits have leakage. This situation is compared to a faultfree case and an uncontrolled situation. Figure 10(a) shows the lateral errors of the vehicle in the three cases. The lateral error of the vehicle increases because of faults and the faulty controlled system tracks the trajectory more accurately than the uncontrolled vehicle. The steering wheel angle and the front wheel steering angle are illustrated in Figures 10(c

CONCLUSION

The paper has proposed the design of a supervisory integrated reconfigurable driver assistance system which is able to track road geometry. The actuators of the control system are the front-wheel steering and the brake yawmoment. The paper extends the control design with an actuator selection procedure, which is built in the design of the supervisor of the system. The control design of actuators is based on the robust optimal LPV method, in which both performance specifications and model uncertainties are taken into consideration. The quadratic stability of the closed-loop LPV system, which contains the individually designed controllers, is guaranteed by a common Lyapunov function. A possible realization of the required control system has also been presented. The integrated system

Fig. 3 .

 3 Fig. 3. Selection of weights ρ st , ρ br

Fig. 4 .

 4 Fig. 4. The low-level driveline control structure

  Fig. 5. Architecture of control system

  Fig. 7. Trajectories of vehicles

  (e) and Figure 8(f). In Figure 8(c) and Figure 8(d) the steering and braking actuations of the controller are shown. If a fault occurs in the steering the actuation of M br and the brake pressures are increased. Figures 8(g

  ) and Figures8(h) illustrate the change in the weighting ρ of controllers.

  Fig. 8. Steering fault compared to the fault-free integrated control yaw moment affects the increased actuation of the front wheel steering.

Fig. 9 .

 9 Fig. 9. Comparison of a fault in the left-hand-side brake circuit with a fault in the right-hand-side brake circuit

Fig. 10 .

 10 Fig. 10. Comparison of the faults in the rear brake circuits with the fault-free integrated control makes it possible to achieve a reconfigurable and faulttolerant system. The fault-tolerance of the controlled system is demonstrated through simulation examples. It can be established that the designed integrated supervisory control system tolerates steering and braking faults by

The research was supported by the Hungarian National Office for Research and Technology through the project "Innovation of distributed driver assistance systems for a commercial vehicles platform" (TECH 08 2/2-2008-0088).