Laurent Bako

René Vidal

Algebraic Identification of MIMO SARX Models

We consider the problem of identifying the parameters of a multipleinput multiple-output switched ARX model with unknown number of submodels of unknown and possibly different orders. This is a very challenging problem because of the strong coupling between the unknown discrete state and the unknown model parameters. We address this challenge by algebraically eliminating the discrete state from the switched system equations. This algebraic procedure leads to a set of hybrid decoupling polynomials on the input-output data, whose coefficients can be identified using linear techniques. The parameters of each subsystem can then be identified from the derivatives of these polynomials. This exact analytical solution, however, comes with an important price in complexity:

The number of coefficients to be identified grows exponentially with the number of outputs and the number of submodels. We address this issue with an alternative scheme in which the input-output data is first projected onto a low-dimensional linear subspace. The projected data is then fit with a single hybrid decoupling polynomial, from which the classification of the data according to the generating submodels can be obtained. The parameters of each submodel are then identified from the input-output data associated with each submodel.

Introduction

Hybrid systems are mathematical models of physical processes which exhibit both continuous and discrete behaviors. Such systems can be thought of as a collection of dynamical submodels with interacting behavior resulting from switching among all the submodels. The switches can be exogenous, deterministic, state-driven, event-driven, time-driven or totally random. Given input-output data generated by such a system, the identification problem consists of determining the parameters of each dynamical submodel as well as those of the switching mechanism (if any). Prior work. Most of the existing hybrid system identification methods have been developed for the class of piecewise auto-regressive exogenous (PWARX) systems [START_REF] Roll | Identification of piecewise affine systems via mixedinteger programming[END_REF][START_REF] Ferrari-Trecate | A clustering technique for the identification of piecewise affine systems[END_REF][START_REF] Ferrari-Trecate | Single-linkage clustering for optimal classification in piecewise affine regression[END_REF][START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF][START_REF] Juloski | A Bayesian approach to identification of hybrid systems[END_REF][START_REF] Bemporad | A bounded-error approach to affine system identification[END_REF], for which the regressor space is partitioned into polyhedral regions with one ARX submodel associated with each polyhedron. For a comprehensive review of hybrid system identification techniques, we refer the readers to the survey paper [START_REF] Paoletti | Identification of hybrid systems: A tutorial[END_REF]. The optimization based method [START_REF] Roll | Identification of piecewise affine systems via mixedinteger programming[END_REF] solves the identification problem as a linear or quadratic mixed integer programming problem. The clustering based procedures [START_REF] Ferrari-Trecate | A clustering technique for the identification of piecewise affine systems[END_REF][START_REF] Ferrari-Trecate | Single-linkage clustering for optimal classification in piecewise affine regression[END_REF][START_REF] Nakada | Identification of piecewise affine systems based on statistical clustering technique[END_REF] use clustering to separate the data into different groups, linear regression to find the boundaries of the polyhedral regions, and linear identification to determine a submodel for each region. Other methods alternate between assigning the data to submodels and estimating simultaneously their parameters by performing a weights learning technique on a fuzzy parameterized model [START_REF] Ragot | Parameter estimation of switching piecewise linear systems[END_REF], solving a Minimum Partition into Feasible Subsystems (MinPFS) problem [START_REF] Bemporad | A bounded-error approach to affine system identification[END_REF] or resorting to Bayesian inference [START_REF] Juloski | A Bayesian approach to identification of hybrid systems[END_REF]. The algebraic approach [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF][START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF] is applicable to the class of Switched ARX (SARX) models, where the switching mechanism can be arbitrary. This approach uses a single decoupling polynomial that vanishes on all the data regardless of their generating submodel. Once this polynomial is computed, the problem reduces to that of recovering the system parameters from the derivatives of the polynomial evaluated at a subset of the regressors.

Unfortunately, most of the aforementioned identification methods can only deal with single-input single-output (SISO) systems. While a few identification methods for multiple-input multiple-output (MIMO) switched linear [START_REF] Bako | Online subspace identification of switching systems with possibly varying orders[END_REF][START_REF] Huang | Identification of hybrid linear time-invariant systems via subspace embedding and segmentation[END_REF][START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF] and piecewise affine [START_REF] Münz | Identification of hybrid systems using a priori knowledge[END_REF][START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF][START_REF] Münz | Continuous optimization approaches to the identification of piecewise affine systems[END_REF] systems in state-space form do exist, they generally require the restrictive assumption of a minimum dwell time in each discrete state. In addition, they often iterate between data clustering and model estimation, which is quite sensitive to initialization. Paper contributions. We present an algebraic solution to the problem of identifying MIMO SARX models. The orders of the submodels are unknown and possibly different and the number of submodels is not available. Our method is based on a technique called Generalized Principal Component Analysis (GPCA) [START_REF] Vidal | Generalized Principal Component Analysis (GPCA)[END_REF], which can cluster data into multiple subspaces by polynomial fitting and differentiation. In contrast to the identification of SISO SARX models [START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF], where only one vanishing polynomial is used to embed the data lying in a mixture of hyperplanes, the identification of MIMO SARX models involves a potentially unknown number n h ≥ 1 of independent homogeneous polynomials that vanish on subspaces of co-dimension higher than one. In order to conveniently construct the regressors to which the embedding is applied, we first estimate the orders of the submodels and the number of discrete states from a rank constraint on the input-output data. Then, given the number of submodels, we compute the number of vanishing polynomials n h and subsequently identify the ARX parameters from the derivatives of these polynomials. However, the number of coefficients to be estimated grows exponentially with the number of outputs and the number of submodels, thereby making the method computationally expensive. We thus propose an alternative method that first partitions the data according to each submodel using a single vanishing polynomial. Given the classification of the data according to each submodel, the parameters of each submodel are then identified using linear techniques.

Problem Statement

We consider a MIMO SARX model of the form

y(t) = n λ t i=1 A i λt y(t -i) + n λ t i=0 B i λt u(t -i) + e(t), (1)
where y(t) ∈ R ny is the output vector, u(t) ∈ R nu is the input vector, λ t ∈ {1, . . . , s} is the discrete state, n λt is the order of the j-th submodel for λ t = j, s is the number of submodels of the SARX system and

A i j i=1,••• ,nj j=1,••• ,s ∈ R ny×ny and B i j i=1,••• ,nj j=1,••• ,s ∈ R ny×nu
are the associated parameter matrices. The modeling error or process noise is represented by e(t) ∈ R ny . In this representation, there may exist for certain models j an integer δ j < n j such that B i j = 0 for i > δ j but we require that A nj j = 0 for all j.

Given input-output data {u(t), y(t)} N t=1 generated by an SARX system of the form (1), and upper bounds on the system orders n ≥ max(n j ) and on the number of submodels s ≥ s, the identification problem can be formulated as follows: identify the number of submodels s, their orders {n j } s j=1 and their parameters

A i j , B i j i=1,••• ,nj j=1,••• ,s .

Algebraic Identification of MIMO Switched ARX Systems

To begin with the identification procedure, let us define the parameter matrices

Γ j = B nj j A nj j • • • B 1 j A 1 j B 0 j A 0 j ∈ R ny×(nj +1)(nu+ny) , P j = 0 ny×qj Γ j ∈ R ny×K , j = 1, • • • , s, (2)
and the regressor vector

x n (t) = u(t -n) ⊤ y(t -n) ⊤ • • • u(t -1) ⊤ y(t -1) ⊤ u(t) ⊤ -y(t) ⊤ ⊤ ∈ R K , (3)
with n = max j (n j ), A 0 j = I ny , q j = (nn j )(n u + n y ) and K = (n + 1) (n u + n y ). For now, assume that the data is not corrupted by noise i.e. e(t) = 0 in (1). Then, the equations defining an SARX system of the form (1) may be re-written as

(P 1 x n (t) = 0) ∨ • • • ∨ (P s x n (t) = 0) , (4)
where ∨ refers to the logical or operator. To eliminate the discrete state from this set of sn y equations, similarly to the case of SISO SARX models [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF], we take the product of one equation per submodel. The advantage of doing so is that we obtain a set of polynomial constraints s j=1 θ ⊤ ij x n (t) = 0, with θ ⊤ ij = P j (i, :) for i = 1, . . . , n y and j = 1, . . . , s, that are satisfied by all the data regardless of their generating submodel. Consequently, the equations in (4) are equivalent to a set of up to n s y (not necessarily independent) homogeneous polynomials p i1,••• ,is on x n (t) of the form

p i1,••• ,is z = s j=1 θ ⊤ ij z = h n1,••• ,nK i1,••• ,is z n1 1 • • • z nK K = h ⊤ i1,••• ,is ν s z . (5)
Here, ν s : R K → R Ms(K) , with M s (K) = K+s-1 s , is the Veronese map which associates to z ∈ R K the vector of all monomials of degree s, z s1

1 • • • z sK K , s 1 + • • • + s K = s,
organized in a descending lexicographic order. Therefore, each p i1,••• ,is is a homogeneous polynomial of degree s with coefficient vector h i1,••• ,is ∈ R Ms(K) and all monomials of degree s in K variables stacked as a vector in ν s (z) ∈ R Ms(K) .

Known Number of Submodels of Known and Equal Orders

In this subsection, we assume that the number of submodels s is known, and that the orders of all the submodels are also known and equal to n. Note that the regressor vectors x n (t) generated by the hybrid model (1) lie in the union of the s subspaces {null(P j )} s j=1 . A basis for each one of these subspaces can be estimated using the GPCA algorithm [START_REF] Vidal | Generalized Principal Component Analysis (GPCA)[END_REF] as follows. From the entire set {u(t), y(t)} N t=1 of input-output data available, construct the matrix of embedded regressor vectors

L(n, s) = ν s x n (n + 1) • • • ν s x n (N ) ⊤ ∈ R (N -n)×Ms(K) . (6)
Then the coefficient vectors h i1,••• ,is of the vanishing polynomials must satisfy

L(n, s)h i1,••• ,is = 0. (7)
In order to solve for the parameters h i1,••• ,is from [START_REF] Paoletti | Identification of hybrid systems: A tutorial[END_REF], one needs to compute the null space of the embedded data matrix L(n, s). Note that h i1,••• ,is is the symmetric part of the tensor product of an indexed set of rows θ ij s j=1 taken from {P j } s j=1 , i.e. ) , where ⊗ denotes the Kronecker product. The linear span of all these coefficient vectors gives a subspace of R Ms(K) that we will refer to as the space of homogeneous polynomials of degree s vanishing on the data. By computing the null space of L(n, s), we obtain a basis for this subspace. In what follows, we will denote such a basis of dimension n h as H = h 1 • • • h n h . Notice that the elements of this basis need not have the structure of a symmetric tensor product.

h i1,••• ,is = Sym (θ i1 ⊗ • • • ⊗ θ is ) ∈ R Ms(K
When the data are perfect and rich enough so that the dimension of the null space of L(n, s) is exactly equal to n h , the matrix of polynomial coefficients H can be computed as a basis for null(L(n, s)) using the Singular Value Decomposition (SVD) of L(n, s). A basis for span(P ⊤ λt ) can then be computed by differentiating the polynomials defined by H at x n (t). The parameter matrix P λt of the submodel generating x n (t) can then be computed as the basis of span(P ⊤ λt ) with an identity matrix at the end, as defined in (2). As we do not need to compute the parameter matrices at each time instant, we can alternatively choose s regressors z j ∈ null(P j ) (see §4.1) and obtain the s parameter matrices {P j } s j=1 from the derivatives of the vanishing polynomials at {z j } s j=1 . Algorithm 1 gives a basic version of the GPCA algorithm [START_REF] Vidal | Generalized Principal Component Analysis (GPCA)[END_REF] for computing the system parameter matrices {P j } s j=1 from input-output data in a deterministic framework. In practice the input-output data may be affected by noise. In this case, even with the assumption that the orders and the number of submodels are known, the matrix L(n, s) is likely to be full rank and so, one may not be able to get the right basis H of polynomials. Therefore, it is desirable to know in advance the dimension n h of this basis. In this way, H could be approximated by the right singular vectors of L(n, s) that correspond to its n h smallest singular values. But since the matrices P j are not known, it is not easy to compute n h in a general framework. However, under certain assumptions on the intersection between the null spaces of the matrices P j , we can derive a closed form formula for n h as outlined in Proposition 1.

Algorithm 1 (Identification of MIMO SARX systems using the GPCA algorithm)

Step 1: Compute a basis H for the null space of L(n, s) by SVD and let the corresponding basis of vanishing polynomials of degree s be Q

(z) = p1(z) • • • pn h (z) = νs(z) ⊤ H. Step 2: Let ∇Q(z) = ∂p 1 (z) ∂z • • • ∂pn h (z) ∂z = ∂νs(z) ∂z ⊤ H.
Step 3: Obtain by SVD a basis Tj ∈ R K×ny for span(P ⊤ j ) as the range space of ∇Q(zj),

j = 1, • • • , s, where zj ∈ null(Pj) but is not in null(Pi), for all i = j. Step 4: Let T ⊤ j = T 1 j T 2 j
be a partition of T ⊤ j such that T 2 j ∈ R ny ×ny . T 2 j is necessarily invertible and we can get Pj = T 2 j Proposition 1 Let H be the symmetric tensor product of a set of matrices B 1 , . . . , B s in R K×m . That is, H is the matrix whose columns are all vectors in R Ms(K) of the form

Sym (b i1 ⊗ • • • ⊗ b is ), where b i1 , . . . , b is are, respectively, columns of B 1 , • • • , B s . If s i=1 rank(B i ) -s < K and for all {i 1 , • • • , i q } ⊂ {1, • • • , s}, q ≤ s, rank( B i1 , • • • , B iq ) = min K, q j=1 rank(B ij ) , (8)
then rank(H) = s j=1 rank(B j ). Assumption ( 8) of Proposition 1 corresponds to an important property of the subspace arrangement ∪ s j=1 null(B ⊤ j ) that is known as transversality. This property states that the dimension of the intersection of any subset of subspaces in the arrangement ∪ s j=1 null(B ⊤ j ) is as small as possible [START_REF] Ma | Estimation of subspace arrangements with applications in modeling and segmenting mixed data[END_REF]. Under this assumption, the number of independent homogeneous polynomials that vanish on ∪ s j=1 null(B ⊤ j ) is equal to rank(H). If the same property holds for ∪ s j=1 null(P j ) and if (n + 1) (n u + n y ) > (s -1) n y , then it follows from Proposition 1 that n h is given by n h = s j=1 rank(P j ) = n s y since rank(P j ) = n y for all j. Although our formula is less general than the one derived in [START_REF] Derksen | Hilbert series of subspaces arrangements[END_REF], it is much easier to compute. In the rest of the section, we will assume that the conditions of Proposition 1 hold, unless stated otherwise.

To summarize, given n and s, the parameter matrices P j follow directly from Algorithm 1. If noise is present in the data, the same algorithm still applies but with the difference that the basis H is approximated by the singular vectors of L(n, s) that are associated with its n h = n s y smallest singular values.

Unknown Number of Submodels of Unknown and Possibly Different Orders

Consider now the more challenging case where neither the orders nor the number of submodels are known and the orders are possibly different. Consequently, n h is also unknown. This means that we need to derive all the parameters of the SARX model (1) directly from the data. In order to properly estimate these parameters, we shall first identify the orders and the number of submodels. Once this task is accomplished, Algorithm 1 can be applied to a certain submatrix of L(n, s) that will be defined later.

Before proceeding further, we need to introduce some notations. For r and l, positive integers, we use the same definitions for x r (t) and L(r, l) as before. Without loss of generality, we denote by n = n 1 ≥ n 2 ≥ • • • ≥ n s the orders of the different submodels that constitute the SARX system and let ρ = n 1 • • • n s ∈ N s be a vector consisting of all the orders enumerated in a non-increasing order. It follows from ( 2) and (3) that the equations defining the SARX model (1) may be re-written as

(Γ 1 x n1 (t) = 0) ∨ • • • ∨ (Γ s x ns (t) = 0) , (9)
where x nj (t) ∈ R Kj , K j = (n j + 1)(n u + n y ) and Γ j ∈ R ny×Kj for j = 1, . . . , s. As before, we may eliminate ∨ in ( 9) by taking the product of one equation per submodel. This leads to a set of polynomial equations on the input-output data of the form

θ ⊤ 1 x n1 (t) • • • θ ⊤ s x ns (t) = h ⊤ η ρ x n (t) , (10)
where θ ⊤ j ∈ R 1×Kj is a row of Γ j , for j = 1, . . . , s, and η ρ x n (t) is a vector obtained from ν s (x n (t)) after removing some of the monomials. η ρ x n (t) does not contain all the monomials, because n j ≤ n for all j = 1, . . . , s, hence x nj (t) is a sub-vector of x n (t), and so the product in [START_REF] Ma | Identification of deterministic switched ARX systems via identification of algebraic varieties[END_REF] does not give rise to all the monomials in ν s (x n (t)).

In order to define the set of monomials that are to be removed, let z = x n (t) and consider a monomial

z α1 1 • • • z αK K , α 1 + • • • + α K = s.
From the definition of x n (t) in (3), it can be seen that the element z αj j is contained in a monomial of η ρ x n (t) if the number of regressors x ni (t) with length K i ≥ K 1j + 1 (that is the number of regressors where z j shows up) is greater or equal to α j . Therefore, in order for the whole monomial z α1 1 • • • z αK K to be included in η ρ x n (t) , we must have that k j ≥ α j for all j = 1, • • • , K, where k j = card ({i :

K i ≥ K 1 -j + 1}).
In view of this analysis, it can shown that the set of monomials to be removed can be indexed by the set

I ρ of exponents (α 1 , • • • , α K ) satisfying α 1 + • • • + α j > k j for j ≤ K 1 -K s .
With this notation, we define a new embedded data matrix in R (N -n)×(Ms(K1)-|Iρ|)

V ρ := η ρ x n (n + 1) , • • • , η ρ x n (N ) ⊤ ( 11
)
that is simply the matrix L(n, s) with |I ρ | missing columns (n = ρ(1)). As before, the null space of V ρ contains the coefficients of the set of vanishing polynomials. However, we may not compute such coefficients directly, because we neither know the system orders ρ nor the number of models s. As it turns out, both ρ and s can be computed from the data under the assumption that the data are rich enough. More specifically:

Definition 1. We say that the data {u(t), y(t)} N t=1 are sufficiently exciting for the SARX system [START_REF] Roll | Identification of piecewise affine systems via mixedinteger programming[END_REF] if the null space of V ρ in [START_REF] Bako | Online subspace identification of switching systems with possibly varying orders[END_REF] is of dimension exactly equal to n h , i.e.

rank(V ρ ) = M s (K 1 ) -n h -|I ρ | . (12)
Notice that Definition 1 assumes implicitly that all the discrete states have been sufficiently visited. If we denote the matrix of data vectors related to the discrete state j by Xj = x n (t j 1 ) • • • x n (t j Nj ) , where the t j k , k = 1, . . . , N j , are the time instants t such that λ t = j, then Xj must span completely null(P j ). Otherwise, null(P j ) may not be identifiable from ∪ s j=1 null(P j ). We have the following result. Theorem 1 Let s ≥ s be an upper bound on the number of submodels and let r be an integer. Assume that the data are sufficiently exciting in the sense of Definition 1. Assume further that N j ≫ M s(K 1 ) for all j = 1, . . . , s. Then dim null(L(r, s)) = 0 if and only if r < max (n j ).

Proof. Assume r < n 1 and let q be the number of submodels whose orders are less than or equal to r.

Let X = x r (t o 1 ), • • • , x r (t o No ) ∈ R f ×No , with f = (r +1
)(n u +n y ), be a matrix whose columns are regressor vectors formed by data generated by the (sq) submodels of orders n j > r. Since the data are sufficiently exciting, X must be full row rank. It follows from Lemma 5 in [START_REF] Vidal | Recursive identification of switched ARX systems[END_REF]

that rank(ν s(X )) = min(N o , M s(f )) = M s(f ), where ν s (X ) = ν s(x r (t o 1 )), • • • , ν s x r (t o No . Consequently, L(r, s) is full column rank, because it is equal to a row permutation of ν s(X ), ν s(X s-q+1 ), • • • , ν s(X s ) ⊤ .
Assume now that r ≥ max(n j ). Then the row nullity of each data matrix X j is at least one. This means that, for all j = 1, . . . , s, there exists a nonzero b j ∈ R f satisfying b ⊤ j X j = 0. One can then verify that Sym(b

1 ⊗• • •⊗b s ⊗a s+1 ⊗• • •⊗a s) ∈ null(L(r, s)) for some a i ∈ R f . Hence, dim(null(L(r, s))) ≥ 1. ⊓ ⊔
Let s ≥ s and n ≥ max(n j ) be upper bounds on the number of submodels and their orders respectively. Thanks to Theorem 1, we can estimate both the number of submodels s and the orders {n j } from the rank of the embedded data matrix L(r, s). The basic idea is that, whenever r is less than one of the orders, there is no polynomial of degree s ≥ s vanishing on the entire data set, provided that N ≫ s and that the data is sufficiently exciting. Therefore, as shown in Algorithm 2, we can obtain the first order n 1 by setting ρ = r • • • r ∈ N s, so that V ρ = L(r, s), and then start decreasing r from r = n to r = 0 until null(V ρ) = {0} for some r * . We then have n 1 = r * + 1.

Given n 1 , we can set ρ = n 1 r • • • r ∈ N s and repeat the procedure starting from r = n 1 and so on, until all the orders of all the s submodels are identified. Notice that, once all the orders of the s submodels have been correctly estimated, r will go to zero for the ss remaining presumed submodels. Therefore, if one assumes that n j > 0 for all j = 1, . . . , s, then the number of submodels can be estimated as the number of orders n j strictly greater than zero.

One advantage of Algorithm 2 is that it does not require knowledge of the dimension n h of the space of vanishing polynomials. If all the orders are correctly identified, then the sufficiency of excitation condition in Definition 1 guarantees that the dimension of the null space of V ρ is exactly equal to n h . Given n h , we can use Algorithm 1 to compute a basis H ρ of null(V ρ ). We can then complete that basis with zeros to form a matrix H ∈ R Ms(K1)×n h such that the rows indexed by I ρ are null. The remaining steps of Algorithm 1 are then performed without additional change.

Algorithm 2 (Identification of the orders and the number of submodels)

Set jo ← 1, nj ← n for j = 1, . . . , s, K ← (n + 1) (nu + ny), V ← L(n, s), 1. Determine the maximum order n1 using Theorem 1 -While rank(V ) < Ms(K), do • nj ← n1 -1 for j = 1, . . . , s • K ← (n1 + 1) (nu + ny) • V ← last Ms(K) columns of V -EndWhile -the maximum order as n1 ← n1 + 1 then set nj ← n1 for j = 1, . . . , s -Set V ← L(n1, s) and K ← (n1 + 1) (nu + ny) 2. Find the remaining orders nj, j = 2, . . . , s using Theorem 1

-jo ← jo + 1 -While rank(V ) < Ms(K) -|Iρ| • nj ← nj o -1 for j = jo, . . . , s • ρ ← n1 • • • ns • Compute Iρ

Implementation of Algorithm 2 with Noisy Data

The algorithm proposed in the previous subsection will operate correctly in the absence of noise. When dealing with noisy data, however, the multiple rank tests required may cause Algorithm 2 to fail, because the involved matrices may always be full rank. In this subsection, we discuss some possible improvements of the algorithm in order to enhance its ability to deal with noisy data.

Recall first that the purpose of the rank test is to check whether or not the dimension of the null space of V ρ is zero for a given vector of orders Therefore, we do not need to know the rank of V ρ exactly. We just need a measure of how likely it is that there exists a nonzero vector h ρ satisfying V ρh = 0.

One possible way of approaching this problem is to inspect the smallest singular value of V ρ for different vectors ρ. For example, to compute n 1 , let ρ1,l = l • • • l ∈ N 1×s , l = 0, . . . , n, and define W ρ1,l . = 1 N -n V ρ1,l V ⊤ ρ1,l as the matrix obtained from 1 N -n L(n, s) ⊤ L(n, s) by removing its columns and rows indexed by I ρ1,l . Denote by σ 1,l , the smallest eigenvalue of the matrix W ρ1,l for l = 0, • • • , n. According to Theorem 1, W ρ1,l has at least one nonzero vector in its null space for all l ≥ n 1 and hence,

σ 1,n1 ≈ • • • ≈ σ 1,n ≈ ε 1,n1 . = 1 n-n1 (σ 1,n1+1 + • • • + σ 1,n
) and are small compared to σ 1,0 , • • • , σ 1,n1-1 . Therefore, to determine n 1 , one needs to look for the smallest integer l ∈ {0, • • • , n} for which σ 1,l ≈ ε 1,l in a certain sense.

Following this procedure, Algorithm 2 can be implemented in a more efficient way for determining the orders. With n0 = n, and given a user-defined decision threshold ε 0 , the following algorithm directly computes the orders starting from j = 1 through j = s, by avoiding the rank tests required in Algorithm 2.

ρj,l = n1 • • • nj-1 l • • • l , l = 0, • • • , nj-1 , σ j,l = min λ W ρj,l , l = 0, • • • , nj-1 , ε j,l = 1 nj-1 -l σ j,l+1 + • • • + σ j,nj-1 , l = 0, • • • , nj-1 , S j = {l = 0, • • • , nj-1 : |σ j,l -ε j,l | < ε o } , nj = min {l : l ∈ S j } , if S j = ∅ nj-1 otherwise, j ← j + 1,
where λ(W ρj,l ) is the set of all eigenvalues of the matrix W ρj,l . In the notation such as ρj,l , the index j indicates which submodel's order is being estimated while l is a possible value of the order sought.

Complexity reduction using a projection approach

The algebraic algorithm proposed in the previous section becomes computationally prohibitive when the dimensions of the SARX system are large. This is because the regressor x n (t) ∈ R K1 constructed from all n y outputs is large, and so it induces an exponential increase in M s (K 1 ), the dimension of the space of homogeneous polynomials space of degree s in K 1 variables. Moreover, the number n h of polynomials to be estimated is unknown, even when the orders and the number of submodels are given, unless one makes certain assumptions. In this section, instead of attempting to compute a potentially large and unknown number of polynomials, we propose a computationally simpler method to identify the model parameters. The idea is to transform the MIMO system into a multiple-input single-output (MISO) system, and hence use only one decoupling polynomial to partition the data according to different ARX submodels. Once all the data are correctly partitioned, the SARX system identification problem reduces to a standard regression problem for each discrete state.

To that end, notice that, without loss of generality, system (1) can transformed into the MISO system3

y(t) = n λ t i=1 a i λt y(t -i) + n λ t i=0 F i λt u(t -i) + e(t), (13)
where the

a i j j=1,••• ,s i=1,••• ,nj are the coefficients of the polynomial z nj -a 1 j z nj -1 -• • •-a nj j
that encodes the poles of the jth submodel as its roots. Let γ = γ 1 • • • γ ny ⊤ be a vector of real nonzero numbers and let y o (t) = γ ⊤ y(t) ∈ R be a weighted combination of all the system outputs. Then, ( 13) can be transformed into the following single output system

y o (t) = n λ t i=1 a i λt y o (t -i) + n λ t i=0 γ ⊤ F i λt u(t -i) + γ ⊤ e(t). (14)
Remark 1 To the purpose of separating the data according to their generating submodels, one may be tempted to consider a single output y j (t) from ( 13) instead of a combination of all the n y outputs. The problem with proceeding in this way is that, after pole-zero cancellation, the MISO system with output y j (t) may be common to many different modes and so, we may not be able to differentiate between those modes. By choosing a random linear combination of the outputs, such degenerate situations can be avoided almost surely.

By introducing the blended output y o (t), we obtain only one hybrid decoupling polynomial g(z) that is easier to deal with. However, at the same time the parameters of different submodels are combined. This raises the question of whether this combination of outputs preserves the distinguishability of the different submodels that constitute the SARX system. In fact, depending on the weights vector γ, two submodels which were initially distinct may reduce to the same submodel in [START_REF] Münz | Identification of hybrid systems using a priori knowledge[END_REF]. To analyze this risk, let

F j = F nj j • • • F 1 j F 0 j ∈ R ny×(nj +1)nu and a j = a nj j • • • a 1 j ⊤ ∈ R nj . ( 15
)
It follows from ( 14) that two different modes i and j become indistinguishable after the previous transformation by γ, if they have the same order (n i = n j ), the same dynamics (a i = a j ) and F ⊤ i -F ⊤ j γ = 0, i.e. when γ lies in null(F ⊤ i -F ⊤ j ). If the F j were known one could readily select a γ which does not satisfy this condition. But these matrices are precisely what we are looking for. The question is, without knowing the F j , how can we choose γ in such a that for any i = j, γ / ∈ null(F ⊤ i -F ⊤ j ). In fact, it is not hard to show that when γ is drawn randomly, this condition is satisfied with probability one. Therefore, two submodels that are distinct in the original system (13) remain so after the transformation. However, the separability of the modes, which is a measure of how close the different submodels are, may be affected.

From [START_REF] Münz | Identification of hybrid systems using a priori knowledge[END_REF], let us redefine the parameter vector θj and the regressor xn (t) as

θj = 0 ⊤ qj γ ⊤ F nj j a nj j • • • γ ⊤ F 1 j a 1 j γ ⊤ F 0 j 1 ⊤ ∈ R K , j = 1, • • • , s (16)
xn (t) = u(t -n) ⊤ y o (t -n) • • • u(t) ⊤ -y o (t) ⊤ ∈ R K , (17)
where K = (n + 1)(n u + 1). One can view the smallest singular value σ 0 (X(γ)) of

X(γ) = xn (n + 1) • • • xn (N )
, as a certain measure of how likely the data can be fitted to one subspace of R K . It is in fact intuitive that the more distinguishable the subspaces are, the larger σ 0 (X(γ)) should be. Therefore, to preserve the separability of the modes, we suggest to choose γ for example as γ * = arg max γ: γ ≤1 σ0(X(γ))

σmax(X(γ)) , where σ max (X(γ)) is the largest singular value of X(γ). Since this could be a hard optimization problem, an alternative is to choose several candidate γs in such a way that σ 0 (X(γ)) is in a certain proportion of σ max (X(γ)).

Once γ has been chosen, we can proceed with the identification procedure. As before, we eliminate the dependency of the system equation on the switches by considering the following decoupling polynomial which vanishes on the data independently of their generating submodel:

g xn (t) = s j=1 θ⊤ j xn (t) = h ⊤ ν s xn (t) = 0. (18)
Solving ( 18) is a particular and simpler case (n y = 1) of the case studied in section 3.

The procedure for the determination of θj is roughly the same:

1. Solve for the orders and number of submodels using Algorithm 2.

2. Obtain h ρ as any nonzero element in null(V ρ ) (which is expected to be one dimensional when the data are sufficiently exciting), and 3. Complete h ρ with zeros to form a h ∈ R Ms(K) so that the entries of h defined by I ρ are zero.

Given h, the parameters may be obtained from the derivative of g as shown in [START_REF] Vidal | An algebraic geometric approach to the identification of a class of linear hybrid systems[END_REF]:

θj = ∇g(z j ) e ⊤ K ∇g(z j ) , j = 1, . . . , s, (19)
where z j is a point in S j \ ∪ s i =j S i , S j = x ∈ R K : θ⊤ j x = 0 , e K is a vector of length K with 1 in its last entry and 0 everywhere else.

Classification of the Data

The computation of θj for each submodel, involves finding a point lying in S j but not in any other S i , i = j = 1, . . . , s. We find a point in S j as z j = xn (τ j ), where

τ j = arg min t∈Dj ∇g (x n (t)) ⊤ xn (t) e ⊤ K ∇g (x n (t)) (20)
D 1 = {t : ∇g(x n (t)) = 0} and D j = {t : ∇g(x n (t)) = 0, θ⊤ i xn (t) = 0, i = 1, ..., j-1}, for j > 1. Then one can compute the parameters by ( 19) using z j = xn (τ j ).

Recall that recovering the vectors θj s j=1 associated with the blended output y o (t) is only an intermediate step in achieving the goal of computing the parameters a j and F j that define each subsystem of the original system [START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF]. Now, from the parameters θj obtained, we can determine the discrete state of ( 14) which is the same as that of [START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF] and then, compute finally the system sought. In order to discard possible outliers in the data we set up a performance bound ε < 1 to the following decision rules:

If min j ∆( θj , xn (t)) > ε xn (t) , then λ t is undecidable. If min j ∆( θj , xn (t)) ≤ ε xn (t) , then λ t = arg min j ∆ θj , xn (t) .
Here ∆( θj , xn (t)) = θ⊤ j xn (t) θj is the distance from the point xn (t) to the linear hyperplane S j defined by its normal vector θj . We define X j = {t > n : λ t = j} = t j 1 , • • • , t j Nj , j = 1, . . . , s as the set of time instances in which the regressors are generated by the submodel j.

Estimation of the Submodel Parameters

Based on the results of the previous classification, we know which data correspond to each generating mode. Therefore, we are left with determining the parameters of each mode j from the data indexed by X j . To begin with, consider a single linear submodel j of order n j from [START_REF] Verdult | Subspace identification of piecewise linear systems[END_REF]. For any t ∈ X j , let us define

Φ y j (t) := y(t -1) • • • y(t -n j ) ∈ R ny×nj , (21)
φ u j (t) := u(t) ⊤ • • • u(t -n j ) ⊤ ⊤ ∈ R (nj +1)nu . (22)
The parameters of the submodels of system (13) can be computed as the solution to the following linear regression problem

y(t) = Φ y j (t) φ u j (t) ⊤ ⊗ I ny a j vec(F j ) + e(t), t ∈ X j . ( 23
)
This equation is obtained by making use of the identity vec(AXB) = (B ⊤ ⊗A)vec(X), where the symbol ⊗ refers to the Kronecker product and vec(•) is the vectorization operator. Notice that in the whole procedure, the vectors a j , j = 1, . . . , s, are estimated twice. The first estimate (obtained from θj ) is considered as a raw estimate that is required here just to be able to discriminate among the different modes. The second estimate from ( 23) is expected to be more accurate.

Numerical results

We test the performance of the proposed approach on an SARX system composed of two submodels of orders 2 and 1, with n u = 1 input and n y = 2 outputs. The system equations are given by

y(t) = a 1 j I ny y(t -1) + a 2 j I ny y(t -2) + b 0 j u(t) + b 1 j u(t -1) + b j u(t -2), (24)
where a 1 j and a 2 j , j = 1, 2, are scalar coefficients and b 0 j , b 1 j , b 2 j are vectors of dimension n y = 2. The coefficients a 2 j and b 2 j are zero for the second submodel. The system is driven by a zero-mean white Gaussian noise input with unit standard deviation and switches periodically from one discrete state to another every 10 samples. The output is corrupted with additive noise with a signal-to-noise ratio (SNR) of 30 dB.

The parameters of the two ARX models are given by the matrices

P 1 = 1.
which are defined with respect to the regressor vector

y(t -1) ⊤ y(t -2) ⊤ u(t) u(t -1) u(t -2) ⊤ .
Given input-output data generated by this system on a time window of size 1500, we are interested in extracting the number of constituent submodels, the orders of these submodels and the parameters that describe them. To demonstrate the performance of our algorithm we carried out a Monte-Carlo simulation of size 1000 with the following user-defined set of parameters: n = 3 and s = 3. For a threshold of ε 0 = 10 -3 in the algorithm of §3.3, the estimation of the orders of both submodels is realized with 100% of successes. Since we provided s = 3, the vector of orders is obtained as ρ = 2 1 0 . The means of the estimates P1 and P2 obtained across all the simulations are given by:

P1 = 1.
Figure 1 shows a histogram with the maximum angle between the column space of the hybrid parameter matrix H and that of its estimate Ĥ. Notice that for all simulations the cosine of this angle is larger than 0.99, implying a strong correlation between H and its estimate. For the second identification method, the result is much better since H consists of only one vector.

Figure 2 shows the relative errors between the true parameter matrices P j and the estimates Pj obtained by our algorithm. Observe that the percentage of simulations that give errors less than 0.05 is about 66% for the first submodel and about 85% for the second submodel. These percentages improve significantly (86% and 93%) when we use the classification approach described in Section 4.

Conclusions

We have presented an algebraic approach to the identification of MIMO SARX models with unknown number of submodels of unknown and possibly different orders. The number of submodels and their orders are estimated from a rank constraint on the inputoutput data, and the model parameters using a subspace clustering technique called GPCA. As the complexity of the method is exponential on the number of outputs and submodels, we proposed a simpler approach that applies GPCA to a MISO system built by projecting the original data. Future work includes developing recursive identification algorithms for MIMO SARX systems, such as the one in [START_REF] Vidal | Recursive identification of switched ARX systems[END_REF] for SISO systems.

  and remove the corresponding columns of V -EndWhile -Obtain the order nj o : nj ← nj o + 1 for j = jo, . . . , s -Set V ← L(n1, s) 3. Go to step 2 until jo = s or until one gets nj o = 0 4. Determine the number of submodels s = card ({j : nj > 0})

Fig. 1 :

 1 Fig.1: Histograms of the maximum subspace angle between span(H) and span( Ĥ).

Fig. 2 : 2 / P 2 2

 222 Fig. 2: Histograms of the errors P 1 -P1 2 / P 1 2 and P 2 -P2 2 / P 2 2 .

j ∈ R ny ×K , j = 1, • • • , s.

Note that the orders nj in (13) may be larger than the ones in[START_REF] Roll | Identification of piecewise affine systems via mixedinteger programming[END_REF]. By an abuse of notation, we will keep using the same notation for the orders.

Acknowledgements. The authors thank Mr. Dheeraj Singaraju for his help in proofreading this paper. This work has been funded by BOURSE-MOBILITE from the Regional Council of Nord-Pas-de-Calais (France), by Johns Hopkins startup funds, and by grants NSF EHS-05-09101, NSF CAREER IIS-04-47739 and ONR N00014-05-1083.


Forests are basically heterogonous, in tree size, tree species, ecological growth conditions etc…

So which sampling strategy should we adopt ?

Turn heterogeneity into Homogeneity

Global rules: stratify the variability ! (to get homogeneous entities on which measurements will be done)

For example, at the scale of a country With the objective of quantifying C stocks, Then the stratification should be made according to :

-ecological regions, -site index within each region, -age within each site index class, -and stand density within each age class…..

For example, at the scale of an eddy-correlation site

With the objective of comparing NEP estimates (from the flux tower and from destructive measurements) Then the stratification should be done according to :

-the age of the stand, -the season within each studied year, -the foot print of the flux tower….

Turn heterogeneity into Homogeneity

Global rules : Drawing of the sampling chain highlighting the objectives of the study

Sampling strategy Estimation of the overall stand volume or biomass

Size of the plot to be inventoried ?

In general, for a given sampling effort :

Number of plots

Plot area

Zone area

If the stand is heterogeneous, it is better to install numerous small plots; while in homogeneous stands, it is better to install few large plots

Frequency

  This step is then objective-dependent and the availability of time, staff and money (!) , case study relating tree volume to D 2 H This model makes the hypothesis that the relationship between tree volume and D 2 H is linear and that the error distribution ε follows a normal law of mean 0 and variance σ 2 It can be fitted using n trees on which volume and D 2 H were measured -This sample also gives The volume prediction for one single tree (knowing D* and H*) is then given by : If now, we want to assess how many trees (n) we should sample to have a precision of E on this tree (D* and H* are known) Making the choice of sampling randomly the trees into the stand, then With µ = and τ the standard deviation of D2H on the sample Taking E=5%, a tree volume of 2m 3 , µ=5m 3 and τ=1m 3 …… n=98 trees Sampling strategy Some theoretical considerations, case study relating tree volume to D 2 H Taking E=5%, a tree volume of 2m 3 , µ=5m 3 and τ=1m 3 …… n=98 trees Is there a way to lower this number of trees, ie, to optimize the sampling ? selection : 4 to 12 trees per stand according to their basal area Reference: Saint-André et al. 2005  Field Constraint Optimum from Laclau, 1996 Sampling strategy Practical applications : volume tables for the main species in France Reference: Vallet et al, 2006 It is recommended to use the tree-basal-area bar charts (equal width of basal area classes) because the largest trees contribute more to the total stand volume than the smallest ones (if volumes are normally distributed) and because the variance of volume, biomass or nutrient content increase with the tree size It is then better to sample the largest trees to catch this variability Note: for multiple-stem trees, use the sum of stem basal area to get the average diameter of the trees; in some particular cases (ex: dry zones) it is better to use the base of the trunk than Same number of trees by basal area classes ? Or proportionality to the number of trees within each class ?  This is the usual case when the objective is to get a robust model (ex: one equation for several ages): Extreme trees have to be sampled to act upon the building and the fitting of the model Same number of trees  This is the usual case when the objective is to get an equation devoted to one particular purpose (ex: tree volume at age 7years) Proportionality to the number of trees within each class Practically, we often use a mixture of these two options  Base with the same number of trees by classes + complementary sampling proportional to the number of trees within each classes Sampling strategy Selecting trees on the field Case of non-homogeneous standsHere again, the sampling to be applied is a function of the wanted accuracy and the available means It is also preferable to avoid particular trees (top broken, sinuous shapes, etc…)EXCEPT if these trees represent a big proportion of the stand, or if the objective is to quantify these defectsLastly it is necessary to avoid trees located near the roads or near the forest limitsSometimes, field constraints change drastically the sampling designed on the computerIt is preferable to select trees within plot that satisfy two main criteria: no high rates of mortality to ensure an homogeneous and representative tree growth clean plots to avoid unwanted roots (herbs, shrubs)Sampling strategy Estimation of the overall stand volume or biomassOnce the biomass equation is built, the question of the sampling strategy for forest inventory remains with and N the total number of trees measured in the plot It can be shown that the previous equation can be rewritten as follows with e denoting the values calculated on the sample used to calibrate the equation And the confidence interval of total volume is given by Sampling strategy Estimation of the overall stand volume or biomass Without stratification, the variance of total volume is given by With stratification, the variance of total volume is given by Where h denotes the strata, N h the total number of trees measured in the strata, n h the number of trees felled in the strata for building the allometric equation Sampling strategy Estimation of the overall stand volume or biomass Size of the plot to be inventoried ? Smaller is the plot= smaller is the representativity in natural Chave, J., R. Condit, et al. (2003). "Spatial and temporal variation of biomass in a tropical forest: results from a large census plot in Panama." Journal of Ecology 91: 240-252. Henry, M., S. Adu-Bredu, et al. (submitted). "Structure and functioning of rain forests ecosystems in Ghana " International Journal of forestry Research.


Ambly-Fleury Vatan

Saint-Cyr-En-Val Moyenvic

Champenoux Pleyber

Blois Nouzilly

Lamballe Kerlavic

Gennes

Michael Nahm

Frank Brodbeck

Laura Van Den Kerchove

Axel Weinreich

  Guémené-Penfao Estrées-Mons

• Wherever possible, exclude work performed by inhouse-efforts from cost-calculations.

• Difficult to give general recommendations for harvesting operations.

Each field site requires an individual approach.

• Develop larger tractor-mounted cutter chipper to reduce costs.

Transport logistics:

Difficult to quantify. Usually, farmers use their own and often very old equipment for short distances.

At present, largest potential to reduce costs: Plan better! Often, harvests were not organized well  raised costs. No apparent effect on stocks.

Decision guideline scheme for practitioners

All: n.s. The low level of implementation of SRC production is caused by:

Michael Nahm Frank Brodbeck

ERA

• lack of knowledge amongst farmers,

• undeveloped markets and unclear quality criteria for the final products,

• absence of regional business co-operations between producers and consumers.

Why farmers do not invest in SRC?

 Technical and organizational uncertainties leading to low efficiency of SRC value chains and high production costs

 Unfavorable profitability level compared to competing agricultural products

Laura Van den Kerchove Axel Weinreich

Leading hypothesis:

-•

  Sites considered together to evaluate the impact of the pedoclimate  And site-specific experiments to evaluate the impact of management practices harvest, spring vs. fall? Do we have to coppice at the end of the first year? Is there an optimal planting density? Is it relevant to chemically fertilize SRC? Is there an optimal planting density? Is sludge spreading a relevant practice? What is the best cultural antecedent? What are/is the "best" species for SRC? What are/is the "best" species for SRC? Do we have to use monoclonal blocks or mixtures? Is wastewater spreading a relevant practice? Without machine transport costs -Without biomass transport costs -Cutter-chipper = most cost-efficient -Still, can only be recommended on sites with good biomass productivity (take establishment, management & recultivation costs into account) In general, the documented harvesting systems work well, including innovative systems like the Stemster, the combination of a feller-buncher with a forwarder, and the small-scale tractor-mounted cutter chipper JF Z 20.

  cut of forage harvester (~ 9.9 cm above ground) 2. Clean cut with chainsaw (~ 2.7 cm above ground) 3. Stocks destroyed with axe 4. Max 4: cut of forage harvester, simulating snow layer (~ 28.5 cm above ground)



  Water content: High (average 58.5% poplar, 52.4% willow)  Calorific value: Gross calorific value high, net calorific value low  Ash content: Averages of poplar (2.4%) and willow (1.6%)  Ash melting behaviour: Deformation temperatures moderately high  Bulk density: Comparably low, values with high variation  Particle size distribution: High variability of shares of fines (<3.15 mm) between 1.1 to 18.8% (average 8.6%)  Bulk density and particle size with high dependency on chipper type Lab scale storage simulation device planned and constructed  Realized prototype is available for lab scale storage trails  For different chipped, shredded or pelletised wood materials  Simulation of different controllable microclimatic scenarios  Access points for diverse physical sensor modules and analysis hard-and soft-ware for assessment of dynamic material parameter development and interactions Findings of online Survey (semi-quantitative) -Survey amongst active and potential energetic or substantial SRC material consumers in Germany: " well known, SRC material not used in a high extent, high willingness of a future use  Rather older material with big diameters (> 3 years), poplar in front of willow  In most cases, chips delivered directly after harvest from short haul distance  High knowledge on SRC water content, but rather low on SRC ash content  No frequent technical problems cited  Storage done under roof or open, with high effects on water content and dry matter development  Price mostly affected by wood chip market and energy prices  Companies rather use own quality standards Conclusions  Compared to residual forest and wooden landscape management material, SRC material has good material properties  Promising fuel characteristics  High interests and rather positive estimations of energetic consumers  Existence of effective storage techniques in terms of improved material quality  Planning of foreseen end use and produced material assortments in front of plantation establishment  Harvesting technique should be adjusted to planned product / material assortment and potentially conducted storage type largest share of costs caused by the rent (24%) for agricultural land and by the establishment of the cultivation (21% of total cost). Also fixed and indirect costs determine (20%) the overall costs to a large extent 2. Small scale plantations more expensive. On a five hectare field the cultivation costs about 5% lower than on a one hectare fields 3. The specific cultivation costs per ha lower for willow (-15%) compared to poplar. However, yield also lower (assumption: 8 vs. 10 t/ha x a) (-20%) 4. Clear effects of the choice of harvesting technique and mode and distance of transport on the costs. Total costs for SRC cultivation may vary by up to 25% 5. Additional information and consultation offers should be provided on SRC, concentrating on giving a deeper understanding of the whole process chain and the long term character of the crop cultivations

Michael Nahm Frank Brodbeck Harvest on marginal field sites: ERA-Net Bioenergy project CREFF Helsinki, February 7-8, 2012 Excel-based tool for practitioners Input Output for planning and calculating harvest + transport on individual field sites, based on project data Michael Nahm Frank Brodbeck ERA-Net Bioenergy project CREFF Helsinki

	Recommendation	Short rotation periods,		Short or long rotation
	limited. Not suitable!	high plant density,		periods, high or low
	Short rotation periods,	large areas, harvest	Long rotation periods,	plant density, harvest
	short transport	with:	area < 1 ha, low plant	motor-manually by
	distances, harvest with	A) Forager	density,	in-house effort or
	forager, immediate	B) Stemster + chip in	high degree of in-	with the Stemster,
	use or interim storage.	summer or autumn.	house effort, motor-	chip in summer or
			manual harvest.	autumn

, February 7-8, 2012 Small burners: Use of dry material Use of material? Own use Sale (plan well) Large burners: Use of fresh or dry material High water availability Low water availability, slope Low water availability High water availability Recommendation limited. Not suitable!

Moyenvic, July 2010 Rottenburg, October 2011 Freiburg, September 2009 Blois, October 2008 Selestat, June 2009 Rottenburg, February 2010 Vosges, February 2009

Pilot Storage Studies SRC material storage efficiency:

Whole shoot drying most efficient, followed by breathable covered chip heap storage  Uncovered chip heap storage is unfavourable Water content reduction (9 to 28%), Net calorific value increase (27 to 94%) Ash content rises, higher shares of fines Ash melting behaviour slightly negative  Under a overall energy balance point of view, instant material use after harvest can be more efficient  But storage can be necessary for many conversion techniques (e.g., small heating plants)  And storage can have a higher added value (higher prices for dry material)

Pilot co-operations

Test one potential approach to overcome the constraints:

The goal of the pilot co-operation is to:

-Bring producers and consumers together -Provide a market with long term contracts and guaranteed prices to the SRC producers -Facilitate information and exchange of experience -Reduce production costs by adapted and adjusted value chains (sharing of transport costs f.ex.) -Reduce individual business risks The low level of implementation of SRC production is caused by:

• Lack of knowledge amongst farmers  Survey: Lack of knowledge is seen by 30% of SRC owners and 53% of farmers without SRC. 52% wish to receive more information.  Pilot co-operations: Lack of know-how addressed as major obstacle for SRC investment for nearly all farmers.

• Undeveloped markets and unclear quality criteria for the final products  Survey: 95% of SRC owners and 67% of farmers without SRC stated "uncertain business conditions for SRC" as major obstacle.

• Lack of knowledge amongst industrial consumers  So far little was known about the quality and suitability of SRC biomass (WP3).

  

Laura Van den Kerchove Axel Weinreich

Confirmation of hypothesis

• Absence of regional business-and logistics co-operations between producers and consumers  Survey: 58% of "farmers without SRC" and 75% of the "SRC owners" strove for business concepts with a high degree of cooperation with the SRC consumer side.  Pilot co-operations: Farmers repeatedly emphasized the importance of the pilot co-operation for providing prices, market security and transfer of information and know-how.

 Survey: Lack of knowledge seen by 30% of SRC owners and 53% of farmers without SRC. 52% wish to receive more information.  Pilot co-operations: Lack of know-how addressed as major obstacle for SRC investment for nearly all farmers.

 Survey: 95% of SRC owners and 67% of farmers without SRC stated "uncertain business conditions for SRC" as major obstacle.

 So far little was known about the quality and suitability of SRC biomass (WP3).

 Survey: 58% of "farmers without SRC" and 75% of the "SRC owners" strove for business concepts with a high degree of co-operation with the SRC consumer side.  Pilot co-operations: Farmers emphasized the importance of the pilot co-operation for providing prices, market security and transfer of information and know-how. Conditioning: HFR Rottenburg -Jan Focke (focke@hs-rottenburg.de) Thorsten Beimgraben (beimgraben@hs-rottenburg.de)

Harvest: FVA Freiburg -Michael Nahm (Michael.Nahm@Forst.bwl.de) Frank Brodbeck (Frank.Brodbeck@forst.bwl.de)