System Reduction and System Identification of Rational Systems Jana Nˇ emcov´ a, ∗ Mih´ aly Petreczky, ∗∗ Jan H. van Schuppen ∗∗∗ ∗

Department of Mathematics, Institute of Chemical Technology Prague, Technick´ a 5, 166 28 Prague 6, Czech republic (email: [email protected]) ∗∗ Department of Computer Science and Automatic Control, Ecole des Mines de Douai, 941, rue Ch. Bourseul, 59508 Douai, France (email: [email protected]) ∗∗∗ CWI, P.O. Box 94079, 1090 GB Amsterdam, The Netherlands (email: [email protected]) Abstract: System identification of polynomial and of rational systems is inherently an approximation problem. In the area of systems biology and also in the area of engineering, there is a need for system identification of these classes of systems. Procedures and algorithms for system identification and model reduction of rational systems are proposed. Examples illustrate the algorithms. Keywords: rational systems, realization theory, identification, model reduction. 1. INTRODUCTION The reader finds in this paper algorithms for system reduction and for system identification of rational or polynomial systems, and an explanation why these algorithms may be useful. The motivation for this paper is the need in systems biology but also in engineering for realistic models of systems based on measurement data. The systems are usually nonlinear. However, for a number of important applications, such as biochemical reaction networks, genetic regulatory networks, etc. the systems of interest can be modelled by so called rational systems. Rational systems are non-linear systems in continuous-time such that the entries of the corresponding vector fields and readout maps are rational functions. By rational functions we mean quotients of two polynomials. Hence, rational systems have a very particular algebraic structure, which can be exploited for performing analysis, control and system identification. The goal of this paper is to present algorithms for system identification and model reduction which exploit the algebric structure of rational systems. More precisely, both algorithms can be obtained from the realization procedure Nˇemcov´ a and van Schuppen (2009, 2010), by making some of its steps computationally effective. The identification algorithm corresponds to choosing the state variables based on time series data. The model reduction procedure corresponds to eliminating unnecessary state variables. There is an extensive literature about system identification of linear systems, of Gaussian systems, and of various subclasses of nonlinear systems Ljung (1999). Most of existing non-linear system identification methods are parametric. The only black-box non-parametric methods we are aware

of is the subspace method for bilinear and LPV systems Verdult (2002). However, the latter method is not directly applicable to rational systems. The parametric identification techniques are based on finding the parameter value as a minimum point of a cost function. The cost function expresses how well the data predicted by the model fits the observations. This method requires the user to choose a parameterization in advance. In addition, computing the minimum points might be difficult, occasionally not possible. The interpretation of the results also requires attention. The algorithms are inspired by the subspace identification algorithm but depend on the algebraic structure of rational functions and of polynomials. A summary of the paper follows. Section 2 contains the formal definition of rational systems and the statement of the model reduction and system identification problems. A procedure for system identification is proposed in Section 4 and a procedure for system reduction is developed in Section 5. Illustrative examples are provided in Section 6.

2. PROBLEM FORMULATION Below we will formally define the class of rational systems and state the problems of system identification and model reduction addressed in this paper. By rational systems we refer to dynamical systems with inputs and outputs, which have a specific structure. The framework used within this paper is adopted from Bartosiewicz (1988, 1987); Nˇemcov´a and van Schuppen (2009). However, for the sake of simplicity we consider the systems only on Rn instead of on irreducible real affine varieties.

We will consider only piecewise-constant functions as inputs of rational systems. For a given input-space U ⊆ Rm we consider the set Upc of all piecewise-constant functions u : [0, tu ] → U , tu ≥ 0. In the sequel, for any u ∈ Upc , tu denotes the real number such that [0, tu ] is the domain of u. Any u ∈ Upc can be identified with a finite sequence u = (α1 , t1 ) · · · (αn , tn ) where αi ∈ U , ti ∈ [0, ∞) for i = Pi Pi+1 1, . . . , n. Then, for t ∈ [ j=0 tj , j=0 tj ), u(t) = αi+1 ∈ U Pn for i = 0, 1, . . . , n − 1, t0 = 0, and u( j=0 tj ) = αn . Note that the same u can have several different representations as a sequence of tuples (α, t) where α ∈ U , t ∈ [0, ∞). Every u = (α1 , t1P ) · · · (αn , tn ) ∈ Upc has a time-domain n [0, tu ] where tu = j=1 tj . If u = (α1 , t1 ) · · · (αn , tn ), v = (β1 , s1 ) · · · (βk , sk ) ∈ Upc then by (u)(v) ∈ Upc we denote the input (α1 , t1 ) · · · (αn , tn )(β1 , s1 ) · · · (βk , sk ). To express that the input u was applied only on the time-domain [0, t] ⊆ [0, tu ] we write a subindex [0, t] to u like u[0,t] . The empty input e is the input with te = 0. Definition 2.1. A rational system Σ with an input-space U and an output-space Rr is a dynamical system with (1) the state-space X = Rn , (2) the dynamics specified as x(t) ˙ = f (x(t), u(t)), where u ∈ Upc , and where f : X × U → Rn is such that for every α ∈ U the components fα,i : X → R, i = 1, . . . , n of the vector field fα : X 3 x → f (x, α) ∈ Rn are rational functions (defined almost everywhere) on X, (3) with the output defined by the map h : X → Rr the components of which are rational functions (defined almost everywhere) on X, (4) the initial state x0 = x(0) ∈ X such that the components of h are defined at x0 . We use the notation Σ = (X, f, h, x0 ). The state trajectory of Σ = (X, f, h, x0 ) corresponding to an input u ∈ Upc is a continuous piecewise-differentiable function xΣ (·; x0 , u) : [0, tu ] → X such that xΣ (0; x0 , u) = d x0 and dt xΣ (t; x0 , u) = f (xΣ (t; x0 , u), u(t)) for t ∈ [0, tu ]. For Σ one defines the set Upc (Σ) ⊆ Upc of admissible inputs u so that the trajectories xΣ (·; x0 , u) are well-defined. In particular, the vector fields fα are defined at x0 , where α refers to the initial constant parts of the inputs of Upc (Σ). In order to define the input-output behavior of a rational system, we have to define the notion of a response map. To this end, we have to recall the notion of an admissible g input set. A subset U pc of the set Upc is called a set of admissible inputs if: g g (i) ∀u ∈ U pc ∀t ∈ [0, tu ] : u[0,t] ∈ Upc , g g (ii) ∀u ∈ Upc ∀α ∈ U ∃t > 0 : (u)(α, t) ∈ U pc , g (iii) ∀u = (α1 , t1 ) · · · (αk , tk ) ∈ Upc ∃δ > 0 ∀ti ∈ [0, ti + δ], i = 1, . . . , k : g u = (α1 , t1 ) · · · (αk , tk ) ∈ U pc . The first two conditions of the definition above assure the well-definedness of the derivative at switching times of differentiable functions which are defined on a set of admissible inputs. In the sequel, unless stated otherwise, g U pc denotes a set of admissible inputs.

r g We say that a function ϕ : U pc → R is a response map if for every component ϕi , i = 1, . . . , r of ϕ and for every sequence of input values α1 , . . . , αk ∈ U , k > 0, the function ϕi;α1 ,...,αk , to be defined below, is analytic.

The definition of ϕi;α1 ,...,αk is as follows. Denote by Tα1 ,...,αk the set of all k-tuples (t1 , . . . , tk ) ∈ [0, +∞)k such g that the input (α1 , t1 )(α2 , t2 ) · · · (αk , tk ) belongs to U pc . We define the map ϕi;α1 ,...,αk : Tα1 ,...,αk → R as follows: for all (t1 , . . . , tk ) ∈ Tα1 ,...,αk , ϕi;α1 ,...,αk (t1 , . . . , tk ) = ϕi ((α1 , t1 ) · · · (αk , tk )). g One can prove that the set A(U pc → R) of all real functions g defined on Upc which are analytic at switching times of the g inputs of U pc is an integral domain. This allows taking the g g quotient field of A(U pc → R), denoted by Q(Upc → R). In the sequel, we will often need the following notation. g If ϕ ∈ A(U pc → R) then by Dα , α ∈ U derivative d of ϕ we mean (Dα ϕ)(u) = dt ϕ((u)(α, t))|t=0+ , where g (u)(α, t) ∈ Upc . A rational system Σ = (X, f, h, x0 ) is said to be a r g realization of the response map p : U pc → R if g g U pc ⊆ Upc (Σ) and p(u) = h(xΣ (tu ; x0 , u)) for all u ∈ Upc . Next, we define the problems of realization, system identification and model reduction. In addition, we also explain the relationship between these problems. Problem 2.2. Realization problem For a specified response r g map p : U pc → R , find a rational system Σ which is a realization of p. A problem closely related to the realization problem is the following problem, which we will call the identification problem. Problem 2.3. Identification problem Consider a response r g map p : U pc → R , a finite collection v1 , . . . , vK ∈ g U pc of piecewise-constant input signals and a finite collection Uc ⊆ U of input values and an integer M . Find an algorithm for computing a rational system Σ from the data {(p(vi ), vi ), (Dα1 Dα2 · · · Dαk p(vi ), vi ) | i = 1, . . . , K, α1 , . . . , αk ∈ Uc , k ≤ M } such that Σ is a realization of p. In other words, we require that a realization of p can be constructed from finitely many values of p and finitely many of the derivatives of p. Note that in the formulation of the identification problem one often requires that the collection {p(v1 ), . . . , p(vK )} represents a time series, i.e. g there exists an input u ∈ U pc and time instances 0 ≤ t1 < . . . < tK < tu such that vi = u|[0,ti ] and one does not take into account the derivatives of p. The latter formulation is a special case of Problem 2.3 and it is more practical. However, we cannot solve this special case yet. What we will show in the sequel is that for each response map p there exists a choice of v1 , . . . , vK , Uc and M such that Problem 2.3 has a solution. Note that in the identification problem, we seek a system Σ which not only re-creates the response of p and its derivatives to vi , i = 1, . . . , K, but it re-creates all the

values of p. In other words, based on responses to finitely many inputs, we would like to find a rational system which correctly describes even the responses to infinite number of inputs. In the definitions above, we were interested in finding rational systems which could realize a certain response map. However, in practice, any model represents the system behavior only in an approximate sense. Hence, it makes sense to replace rational systems with a large number of states with rational system with a smaller number of states such that the response map of the latter system is an approximation of the response map of the original system. This leads us to the problem of model reduction. Problem 2.4. System reduction of rational systems. Consider a rational system Σ = (Rn , f, h, x0 ). Construct anˆ x ˆ = (Rn1 , fˆ, h, other rational system Σ ˆ0 ) such that

g (1) Find ϕ1 , . . . , ϕk ∈ A(U pc → R) such that Qobs (p) = R(ϕ1 , . . . , ϕk ). (2) Compute viα ∈ R(X1 , . . . , Xk ) for all α ∈ U , i = 1, . . . , k such that Dα ϕi = viα (ϕ1 , . . . , ϕk ). (1) The existence of such viα follows from the fact that Qobs (p) is closed with respect to Dα derivatives. (3) Compute wj ∈ R(X1 , . . . , Xk ) for j = 1, . . . , r such that pj = wj (ϕ1 , . . . , ϕk ). (2) The existence of such wj follows from the fact that pj ∈ Qobs (p). (4) Define the the rational system Σ where X = Rk , fα =

k X

viα

∂ , α ∈ U, ∂xi

(3) i=1 ˆ (1) admissible inputs to Σ and to Σ are the same, hj = wj , j = 1 . . . r, (2) n1 < n (ideally n1  n), x0 = (ϕ1 (e), . . . , ϕk (e)), e is the empty input. ˆ are close with respect (3) The response maps of Σ and Σ If fα and h are defined at x0 then Σ is a rational system to some criterion C : Upc (Σ) → [0, +∞) realizing p. ˆ ˆ (tu ; x ∀u ∈ Upc (Σ) : |h(xΣ (tu ; x0 , u))−h(x ˆ0 , u)| < C(u) Σ

In fact, the problems of model reduction and system identification are often combined to the following problem: based on finitely many measurements, find a system model of certain size such that the model approximates the true input-output behavior with a certain accuracy. This problem can also be formulated in a formal manner for rational systems. We will not present the formal definition in this paper, since we will not deal with this problem explicitly. Note, however, that the combination of the system identification algorithm and the model reduction, to be presented in the sequel, always yields a solution to the latter problem. The planned approach is to use the realization procedure from Nˇemcov´ a and van Schuppen (2009). 3. REALIZATION PROCEDURE FOR RATIONAL SYSTEMS Below we recall from (Nˇemcov´ a and van Schuppen, 2009, Proposition 5.14) the procedure for computing a rational realization of a response map. This procedure is not computationally effective, but it will play an important role in the subsequent system identification and model reduction algorithms. In order to present this procedure, we need to introduce the following definition. r g Definition 3.1. Let p : U pc → R be a response map. The observation algebra Aobs (p) of p is the smallest subalgebra g of A(U pc → R) which contains the components pi , i = 1, . . . , r of p, and which is closed with respect to Dα derivatives for all α ∈ U . The observation field Qobs (p) of p is the field of quotients of Aobs (p), it is a subfield of g Q(U pc → R). r g Procedure 3.2. Consider a response map p : U pc → R which is realizable by a rational system. Thus, Qobs (p) is g a finite field extension of R by the elements of A(U pc → R), see (Nˇemcov´ a and van Schuppen, 2009, Theorem 5.16).

4. SYSTEM IDENTIFICATION ALGORITHM In this section we present an algorithm for solving the identification problem. More precisely, we will be able to show that the algorithm solves Problem 2.3 for some choice of inputs and higher-order derivatives of the response map. To this end, we would like to use the realization procedure from Procedure 3.2 to compute a realization from input-output data. In order to turn Procedure 3.2 into an algorithm, the step of finding a generator set {ϕ1 , . . . , ϕk } of the observation field Qobs (p), and finding the rational functions viα , wj , i = 1, . . . , k, j = 1, . . . , r has to be made computationally effective. We start with the following simple statement, which has already appeared in Nˇemcov´a and Pomet (2010); Wang and Sontag (1992). Theorem 4.1. (Nˇemcov´a and Pomet (2010),Wang and Sontag (1992)). Assume the set of inputs values U is finite. Then the following holds. (1) The field Qobs (p) is finitely generated if and only if the transcendence degree of Aobs (p) is finite. (2) Assume x1 , . . . , xd is a transcendence basis of Aobs (p). Then the field Qobs (p) is generated by the finite set G, where G = {p1 , . . . , pr , xi , Dα xi | i = 1, . . . d, α ∈ U }. (4) The proof of Theorem 4.1 yields in fact the following proposition. Proposition 4.2. Assume that U is finite and that G from (4) is of the form G = {ϕ1 , . . . , ϕk }, i.e. Qobs (p) = R(ϕ1 , . . . , ϕk ). If we apply Procedure 3.2 with G as the set of generators, then we obtain a rational system Σ of the form (3), where the rational functions wj , viα ∈ R(X1 , . . . , Xk ), j = 1, 2, . . . , r and i = 1, . . . , k from (1) and (2) can be chosen as follows. For every g ∈ G,denote by Xg the indeterminate Xl for which ϕl = g. Then • For each j = 1, . . . , r, wj = Xpj . • If ϕi = xl for some l = 1, . . . , d, then let viα = XDα xl .

• If ϕi = Dβ xl for some β ∈ U , l = 1, . . . , d or ϕi = pj , j = 1, . . . , r, then define viα as follows. Consider the polynomial Q(Z0 , Z1 , . . . , Zk ) in variables Z0 , . . . , Zk such that Q(ϕi , x1 , . . . , xk ) = 0. Denote by ∂i Q the partial derivative of Q with respect to Zi , i = 0, . . . , k. ∂0 Q(Xi , Xx1 , . . . , Xxk ) viα = − Pk . j=1 XDα xj ∂j Q(Xi , Xx1 , . . . , Xxd ) In the sequel, we will focus on finding a transcendence basis {x1 , . . . , xd } of Aobs (p) and finding the polynomials Qj , Qα,i in d + 1 variables such that Qj (pj , x1 , . . . , xd ) = 0, j = 1, . . . , p, Qα,i (Dα xi , x1 , . . . , xd ) = 0, i = 1, . . . , d, α ∈ U.

(5)

If these polynomials are known then Proposition 4.2 and Procedure 3.2 yield an identification algorithm. To this end, assume that U is finite and for every integer M > 0 define SM = {p1 , . . . , pr , Dα1 · · · Dαl pj | α1 , . . . , αl ∈ U, j = 1, . . . , r, l ≤ M }. It is clear that the algebra generated by SM is a subalgebra of Aobs (p) and that for some M , a transcendence basis of Aobs (p) can be chosen from the elements of SM . Procedure 4.3. (1) Assume that |SM | = D, |SM +1 | = H, and SM +1 = {ψ1 , . . . , ψH } and ψ1 , . . . , ψH are chosen in such a way that SM = {ψ1 , . . . , ψD }. Consider the ideals I = {P ∈ R[X1 , . . . , XH ] | P (ψ1 , . . . , ψH ) = 0} J = I ∩ R[X1 , . . . , XD ] Notice that the algebra generated by SM is isomorphic to R[X1 , . . . , XD ]/J and the algebra generated by SM +1 is isomorphic R[X1 , . . . , XH ]/I. (2) Using the Noether’s normalization theorem and the corresponding algorithm (Eisenbud (1995)), find algebraically independent polynomials b1 , . . . , bd , bd+1 , . . . , bD such that R[X1 , . . . , XD ] is integral over S = R[b1 , . . . , bD ] and the ideal S ∩ J is generated by bd+1 , . . . , bD . (3) Define xi = bi (ψ1 , . . . , ψD ), i = 1, . . . , d. Then {x1 , . . . , xd } is a transcendence basis of the algebra generated by SM . (4) Find the polynomials Pi,α , Pj ∈ R(X1 , . . . , XH ) such that Dα xi = Pi,α (ψ1 , . . . , ψH ) pj = Pj (ψ1 , . . . , ψH ) for all j = 1, . . . , r, α ∈ U , i = 1, . . . , d. This can be done as follows. For every g ∈ SM +1 , let Xg denote the indeterminate Xl if g = ψl . Then Pj = Xpj and Pi,α =

D X l=1

XDα xi

∂bi ∂Xl

(5) For every P ∈ {Pj , Pi,α | i = 1, . . . , d, α ∈ U, j = 1, . . . , r} define the ideal KP as the ideal of the ring R[Z0 , Z1 , . . . Zd , X1 , . . . , XH ] generated by P − Z0 , bi − Zi , i = 1, . . . , d and the elements of I. Here Z0 , . . . , Zd are new indeterminates. Define Qj as any element of KQj ∩ R[Z0 , Z1 , . . . , Zd ] and define Qi,α as any element of KPi,α ∩ R[Z0 , Z1 , . . . , Zd ]. Lemma 4.4. If Aobs (p) has finite transcendence basis, then there exists M ∈ N such that Procedure 4.3 returns a transcendence basis {x1 , . . . , xd } of Aobs (p) and polynomials Qj , Qi,α , i = 1, . . . , d, α ∈ U satisfying (5).

The only step of Procedure 4.3 which is not computationally effective is the computation of a finite generator set of I. All the other steps can be executed using standard algorithms from computer algebra Eisenbud (1995). However, I can be approximated as follows. Procedure 4.5. (Approximation of I). (1) Suppose that we can evaluate ψ1 , . . . , ψD for inputs v1 , . . . , vK and define the points di = (ψ1 (vi ), . . . , ψD (vi )), i = 1, . . . , K. Then every element P of I satisfies P (di ) = 0, i = 1, . . . , K, i.e. I is contained in the affine variety spanned by the points d1 , . . . , dK . (2) Consider then the ideal IK = (X1 − ψ1 (vi ), . . . , XD − ψD (vi ) | i = 1, . . . , K). (3) Fix an integer R and consider all the polynomials F R of IK of degree less or equal than R. (4) Return IR,K as the ideal generated by the elements of F R , i.e. take F R as the finite generator of I. From ((Ma et al., 2007, theorem 2.9)) we obtain the following result. Lemma 4.6. There exists a choice of R, K and values v1 , . . . , vK such that the ideal IR,K returned by Procedure 4.5 equals I. Hence, in principle Procedure 4.5 is capable of computing I. Note that the ideal IK computed by Procedure 4.5 will always contain I and hence for large enough R the result of Procedure 4.5 will always contain I. Hence, we can sum up the results as follows. Procedure 4.7. (System identification). (1) Use Procedure 4.5 for some R and K to obtain an approximation IR,K of the ideal I. (2) Use Procedure 4.3 with IR,K instead of I to compute a transcendence basis {x1 , . . . , xd } and polynomials Qj , Qi,α , j = 1, . . . , r, i = 1, . . . , d, α ∈ U . (3) Use the polynomials Qj , Qi,α j = 1, . . . , r, i = 1, . . . , d, α ∈ U to compute the rational maps wj , viα as described in Proposition 4.2 and return the rational system Σ = (X, f, h, x0 ) described in Proposition 4.2. Procedure 4.7 is completely effective and can be implemented using standard algorithms of computer algebra. Combining Proposition 4.2, Lemma 4.4 and Lemma 4.6 we obtain the main result of this section. Theorem 4.8. (Main result). Assume U is finite. There g exists a choice of inputs v1 , . . . , vK ∈ U pc and integers K, R, M such that Procedure 4.7 returns a rational system Σ which is a realization of p. Theorem 4.8 applies only to the case when the set of input values U is finite. This is quite a restrictive assumption. Below we show that Theorem 4.8 and Procedure 4.7 can be extended to the more general case of input-affine systems. Definition 4.9. (Input-affine systems). A rational system m is input-affine, Pm if for any α = (α1 , . . . , αm ) ∈ R = U, f (x, α) = j=0 fj (x)αj , where α0 = 1 and fj (x), j = 0, 1, . . . , m are rational functions. It then follows from Wang and Sontag (1992) that if p has a realization by an input-affine rational system, then p admits a Fliess-series expansion, i.e. with the notation of Wang and Sontag (1992) there exists a generating series

g such that Fc [u](tu ) = p(u) for any u ∈ U pc . This remark prompts us to introduce the following definition. Definition 4.10. We call a response map p input-affine, if the zero vector and the standard basis vectors ej ∈ Rm , j = 1, . . . , m all belong to U and there exists a convergent generating series c : {0, 1, . . . , m}∗ → Rr such that p(u) = g Fc [u](tu ) for all u ∈ U pc . Here we used the notation and terminology of Wang and Sontag (1992). If a response map admits a realization by an input-affine rational system, then this response map is input-affine. Note that input-affine systems often occur in biological models; if one thinks of α as an enzyme concentration, then it is clear that enzyme concentrations enter the equations of mass-action kinetics in a linear fashion. If p is an inputaffine response map, then denote by pˆ the response map which is obtained by restricting p to piecewise-constant g inputs u ∈ U pc such that u takes values in the set Ud = {0, e1 , . . . , em }. We can easily derive the following result. ˆ = (X, fˆ, h, x0 ) is Proposition 4.11. A rational system Σ a realization of the response map pˆ, if and only if Σ = (X, f, h, x0 ) is a realization of p, where m X ∀α = (α1 , . . . , αm ) ∈ U : fα = fˆ0 + αj fˆj . j=1

Notice that pˆ is defined only for input functions which take values from a finite set. Hence, Proposition 4.11 and Theorem 4.8 yield the following. Corollary 4.12. There exists a choice of inputs v1 , . . . , vK ∈ g U pc and integers K, R, M such that Procedure 4.7 returns ˆ which is a realization of pˆ. Hence, the a rational system Σ corresponding system Σ described in Proposition 4.11 will be a realization of p. In fact, by repeating the proof of Theorem 4.8, it can be shown that the choice of the input signals v1 , . . . , vK need not be restricted to those input signals which take their values from the set {0, e1 , . . . , em }. 5. MODEL REDUCTION Below we present a model reduction procedure based on Procedure 3.2. The main idea is as follows. For every rational system Σ = (X, f, h, x0 ), we can define the response map p : Upc (Σ) → Rr , defined by ∀u ∈ Upc (Σ) : p(u) = h(xΣ (tu ; x0 , u)). In model reduction, a candidate generator set for Qobs (p) is computed from the state-space realization Σ = (X, f, h, x0 ). This candidate generator set and the realization procedure are then used to compute a rational system Σr , which should approximate p. In order to compute a generator set of Qobs (p) from Σ, one can use the results from (Nˇemcov´ a and van Schuppen (2009)), stating that Qobs (p) is the image of Qobs (Σ) under the map τΣ,x0 . One can then compute a generator set of the observation field Qobs (Σ) and take its image by the map τΣ,x0 . To derive the generators of Qobs (Σ), notice that Qobs (Σ) is a finite field extension of R. Hence, there exists an integer I such that Qobs (Σ) is generated by the map of the form Lfα1 Lfα2 · · · Lfαl hj , j = 1, . . . , r, α1 , . . . , αl ∈ U , l ≤ I. Here, Lf g denotes the Lie-derivative of the map g with respect to the vector field f .

Finally, below we present a procedure which computes a rational system whose input-output map is an approximation of p. The procedure consists of the following steps. (1) It takes a finite set of generators ϕ1 , . . . , ϕk of Qobs (p) as inputs. (2) It eliminates one of these generators. More precisely, it eliminates the generator which is the closest one to the subfield generated by the other generators. (3) Using the remaining generators, the procedure computes a rational state-space realization Σr . Intuitively, Σr should be an approximate realization of the response map at hand. Unfortunately, the precise relationship between Σ and Σr is not known yet. In order to formalize the above procedure a norm or distance is needed between response maps. We will not elaborate on the choice of the norm here. One option is that to assume that the input-output maps of interest take bounded values and to use a version of supremum norm. In any case, assume that we managed to fix a norm k.k. Now we are ready to present the proposed procedure. Procedure 5.1. Reduction of a subfield by one generator. Assume that a generator set ϕ1 , . . . , ϕk ∈ Qobs (p) is specified. (1) Calculate for i = 1, . . . , k the norms kϕi − Fk\i k, where, Fk\i = R({ϕj , j = 1, . . . , k, j 6= i}). (2) Calculate, j ∗ = argmini=1,...,k kϕi − Fk\i k. (3) Set the new generator set as {ϕi | i = 1, . . . , k, i 6= j ∗ } and by possibly re-labelling the generator set, assume that j ∗ = k (4) Calculate for i = 1, 2, . . . , k − 1 and every α ∈ U , the r rational maps fα,i (X1 , . . . , Xk−1 ) which minimize the distance r kDα φi (.) − fα,i (ϕ1 , . . . , ϕk−1 )k

(5) For every i = 1, . . . , r, find the rational function hri (X1 , . . . , Xk−1 ) which minimizes the distance kpi − hri (ϕ1 , . . . , ϕk−1 )k (6) Calculate the reduced order system Σr = (Rk−1 , f r , hr , xr0 ). Pk−1 r ∂ r r r where fαr = i=1 fα,i ∂xi , α ∈ U , h = (h1 , . . . , hr ), r x0 = (ϕ1 (e), . . . , ϕk−1 (e)). The procedure can also be used with several iterations resulting in a set of generators which is two or more less than that of the original set. This procedure is not an algorithm, since it contains steps which are not computationally effective. Turning this procedure into an algorithm remains topic of future research. 6. EXAMPLES There are two examples presented in this section. The first example shows how the choice of the generators of the observation field of a response map influences the structure of its realization. The proposed system reduction

algorithm is demonstrated in the second example. We show that if one chooses the generators of the observation field to be algebraically independent rather than algebraically dependent, then one obtains a realization with lower dimension. Example 6.1. Consider the response map p given as j−1

p(u) =

X αi αj 2 2 2 (Tj − Tj−1 ) + tj + ( (Ti2 − Ti−1 ) + ti ) + 1 2 2 i=1

g for u = (α1 , t1 )(α2 , t2 ) · · · (αj , tj ) ∈ U pc and T0 = 0, Ti = Pi t , i = 1, . . . , j. Let us derive a rational realization Σ l=1 l of p by following the steps of Procedure 3.2. First, note that for arbitrary α, β, γ ∈ R it holds that (Dα p)(u) = αtu + 1, (Dβ Dα p)(u) = α, (Dγ Dβ Dα p)(u) = 0. If one considers ϕ1 = p, ϕ2 = tu and applies Procedure 3.2 to the field Qobs (p) = R(ϕ1 , ϕ2 ), one obtains the realization Σ1 = (R2 , f, h, x(0)) of p, where f (x1 , x2 , u) = (ux2 + 1, 1)T , h(x1 , x2 ) = x1 , x(0) = (1, 0)T . However, if one chooses ϕ1 = αtu + 1, ϕ2 = p it still holds that Qobs (p) = R(ϕ1 , ϕ2 ) and one gets ˆ x Σ2 = (R2 , fˆ, h, ˆ(0)),

Because Qobs (p) = R({p, Dα1 · · · Dαj p|α1 , . . . , αj ∈ R, j ∈ N}) and because 1 (Dα p)(u) = αp(u) , (1 + tu )2 −2 1 + αp(u) , (Dβ Dα p)(u) = αβp(u) (1 + tu )4 (1 + tu )3 one derives that Qobs (p) = R(p, 1 + tu ). Note that ϕ1 = p and ϕ2 = 1 + tu are algebraically independent over R. By following the steps of Procedure 3.2 with this choice of generators of Qobs (p) we derive the realization Σr = (R2 , f , h, x(0)), where f (x1 , x2 , u) = (u xx12 , 1)T , 2

h(x1 , x2 ) = x1 , x(0) = (1, 1)T . Since Σ and Σr have the same input-output behavior and nr = 2 < 3 = n, Σr is a reduced model of a phenomenon modeled by Σ. 7. CONCLUDING REMARKS The paper provides procedures for system reduction of rational systems and a procedure for system identification of rational systems. Examples of the procedures are provided. More research is required to clarify in particular the properties of the approximations. REFERENCES

C. A. Anthoulas. Approximation of Large-Scale Dynamical Systems. SIAM, 2005. ˆ 1 , x2 ) = x1 , x where fˆ(x1 , x2 , u) = (u, x1 )T , h(x ˆ(0) = Z. Bartosiewicz. Rational systems and observation fields. Systems & Control Lett., 9:379–386, 1987. (1, 1)T . Z. Bartosiewicz. Minimal polynomial realizations. Math. Control One can try to choose the generators of Qobs (p) so that the Signals Systems, 1:227 – 237, 1988. dynamics of the computed realization is relatively simple. S. Basu, R. Pollack, and M.-F. Roy. Algorithms in real algebraic geometry. Number 10 in Algorithms and Computation in MatheIn our example the choice ϕ1 = u, ϕ2 = tu simplifies the matics. Springer, Berlin, 2003. ISBN 3-540-00973-6. dynamics of the computed realization as follows: N. Bourbaki. Algebra II: Chapters 4-7. Springer, Berlin, 1990. 2 ¯ ¯ Σ3 = (R , f , h, x ¯(0)), David Eisenbud. Commutative algebra – With a view toward algebraic geometry. Number 150 in Graduate Texts in Mathematics. x1 x22 T ¯ ¯ where f (x1 , x2 , u) = (0, 1) , h(x1 , x2 ) = 2 + x2 + 1, Springer, Berlin, 1995. x ¯(0) = (0, 0)T . Note that the dynamics of the system got N. Jacobson. Basic algebra, Volumes 2. W.H. Freeman and Company, New York, 1980. simpler, but the output function got more complicated. Example 6.2. Consider the system Σ = (R3 , f, h, x(0)), L. Ljung. System Identification: Theory for the user (2nd Ed.). PTR Prentice Hall., Upper Saddle River, USA, 1999. where f (x1 , x2 , x3 , u) = (ux1 x2 , −2x3 , −3x22 )T , h(x1 , x2 , x3 ) =Y. Ma, A. Yang, H. Derksen, and R. Fossum. Estimation of subspace T x1 and x(0) = (1, 1, 1) . One can check that Σ is a arrangements with applictions in modeling and segmenting mixed realization of the response map data. SIAM Rev., vol.50, no. 3, August 2008. Z tu J. Nˇ emcov´ a and J.H. van Schuppen. Realization theory for rational u(s) ds). p(u) = exp( systems: Minimal rational realizations. Acta Applicandae Mathe2 (1 + s) 0 maticae, 110:605–626, 2010. To construct Σ from p by following Procedure 3.2 we need Jana Nˇemcov´a and Jean Baptiste Pomet. On the existence of various realizations. In Proceedings of the International Symto choose the following generators of Qobs (p): posium on the Mathematical Theory of Networks and Systems 1 1 (MTNS.2010), Budapest, 2010. University of Budapest. ϕ1 = p, ϕ2 = , ϕ3 = . 2 (1 + tu ) (1 + tu )3 Jana Nˇ emcov´ a and Jan H. van Schuppen. Realization theory for rational systems: The existence of rational realizations. SIAM J. Then Qobs (p) = R(ϕ1 , ϕ2 , ϕ3 ) and Control & Opt., 48:2840 – 2856, 2009. αp V. Verdult. Nonlinear System Identification: A State-Space ApDα ϕ1 = Dα p = = αϕ ϕ , 1 2 (1 + tu )2 proach. PhD thesis, University of Twente, 2002. Yuan Wang and Eduardo D. Sontag. Algebraic differential equations d 1 −2 Dα ϕ2 = |t=0+ = = −2ϕ3 , and rational control systems. SIAM J. Control & Opt., 30:1126– dt (1 + tu + t)2 (1 + tu )3 1149, 1992.

d 1 −3 |t=0+ = = −3ϕ22 , 3 dt (1 + tu + t) (1 + tu )4 w(ϕ1 , ϕ2 , ϕ3 ) = p = ϕ1 , which yields Σ. Note that ϕ2 and ϕ3 are algebraically dependent over R. Dα ϕ3 =