Addition and intersection of linear time-invariant behaviors

We define and analyze the operations of addition and intersection of linear time-invariant systems in the behavioral setting, where systems are viewed as sets of trajectories rather than input-output maps. The classical definition of addition of input-output systems is addition of the outputs with the inputs being equal. In the behavioral setting, addition of systems is defined as addition of all variables. Intersection of linear time-invariant systems was considered before only for the autonomous case in the context of"common dynamics"estimation. We generalize the notion of common dynamics to open systems (systems with inputs) as intersection of behaviors. This is done by proposing trajectory-based definitions. The main results of the paper are 1) characterization of the link between the complexities (number of inputs and order) of the sum and intersection systems, 2) algorithms for computing their kernel and image representations and 3) a duality property of the two operations. Our approach combines polynomial and numerical linear algebra computations.


Introduction
The behavioral setting [1,2] is an approach to system theory where systems are defined as sets of trajectories rather than input-output maps. This fact has the following far-reaching consequences.
• The system and its representations are separated: a representation is an equation while the system is the solution set. Equivalence of representations is naturally defined then as equality of the corresponding solution sets.
• The variables of interest are not a priori separated into inputs and outputs. An input/output partitioning is in general not unique and may not be a priori given, however, the classical approach imposes a fixed one. This may lead to inconsistencies as shown in [2].
In this paper we analyze the basic operations of addition and intersection of linear time-invariant (LTI) systems in the behavioral setting. We provide trajectory-based definitions of these operations and algorithms that compute their representations starting from representations of the original systems.
The notion of addition in the behavioral setting differs from the one in the classical setting. While in the latter we only sum the outputs by leaving the inputs the same, in the former, we add all variables (inputs and outputs). The notion of intersection is not even well-defined in the input-output setting. Only the particular case of intersection of autonomous systems is considered in the context of "common dynamics" estimation [3,4,5,6].

Literature overview
The operations of addition and intersection are not new to the literature. For two-dimensional systems [7,8], conditions which allow representing of the system as a direct sum of two systems are given in [9]. For one-dimensional systems (the topic of the paper), the sum is used decompositions such as • controllable and autonomous [1]; • stable and unstable [10]; and • zero input and zero initial condition.
The decomposition into subsystems also appears in the modeling by tearing, zooming and linking [2].
The intersection operation is used in control in the behavioral setting to restrict the plant behavior by a controller. Control by behaviors intersection is equivalent to what is called full interconnection [11], that is the case when the sets of control variables and to be controlled variables coincide. This is a special case of interconnection [12] in the behavioral setting, where the constraints are imposed only on a subset of the variables.
Computation of a kernel representation of the sum of two systems is given in [13,Lemma 2.14]. Similar result for the intersection operation is given in Section 3. We remark that all the existing works involving the addition and intersection of linear time-invariant systems compute them via their kernel representations only.

Contribution and organization of the paper
We study the addition and intersection operations in the behavioral framework. In particular, we propose a trajectory-based definition for these operations, which allows performing these computations directly from the observed data. We characterize their representations in terms of the original systems representations and thus derive new algorithms for their computations. Also the approach that we use is new: we use numerical linear algebra operations based on structured matrices rather than polynomial algebra operations. This makes our results readily implementable in practice.
The rest of the paper is organized as follows: Section 2 reviews the main results and definitions of behavioral system theory used in the paper. We define the trajectory-based operations of addition and intersection of behaviors in Section 3, and we propose new algorithms for the computation of their representation in Section 4. Illustrative examples with simulations are given in Section 5.

Notation and preliminaries
In this section, we recall some useful notions and definitions of the behavioral approach needed in the rest of the paper.
A dynamical system is defined by a triple (T , W, B), where T ⊆ N is the time span, W is the signal space, the space of values for the system trajectories, and B ⊆ W T is the set of admissible trajectories for the system, the behavior. The choice of T highlights that, in the paper, we focus on discrete time systems.
The system (T , W, B) will be denoted simply by its behavior B. B ∈ L q means that B is a linear time-invariant system with q variables, where q is the sum of the number of inputs (the free variables), m = m(B), and the number of outputs, p = p(B).
Linear time-invariant behaviors are described by systems of difference equations, and they admit several representations (remember that the behavior and its representations are decoupled in this framework, but representations are helpful from the computational point of view). The kernel representation is the most important in the paper. Every B ∈ L q admits a kernel representation [14], i.e., B consists of the set of trajectories that belong to the kernel of a certain matrix polynomial operator R(σ), where σ is the shift operator (σw)(t) = w(t + 1).
where the operator R(σ) is represented by the matrix polynomial The variables can be partitioned element-wise into inputs and outputs w = Π [ u y ], where Π is a permutation matrix. This leads to the input-output representation where In Lemma 1, the matrix polynomial R(z) has exactly p(B) rows, in which case the kernel representation is called minimal. A minimal kernel representation corresponds to a full row rank R(z). Every kernel representation can be reduced to a minimal one by a suitable transformation. Remark 1. A kernel representation R(σ) is not unique. An equivalent one can be obtained by premultiplication of R(z) with a unimodular matrix polynomial U (z), which is a square matrix polynomial whose determinant is a nonzero constant. In a minimal kernel representation deg R = ℓ. In the follows, we assume that the kernel representations are minimal. Definition 1. The lag of a system, ℓ = ℓ(B), is the (minimum) degree of the matrix polynomial R(z).
The order of a system, n = n(B), is the number of zeros of det(P (z)), where P (σ) is the square matrix polynomial operator in (2).
The following concept of controllability (in the behavioral setting) is needed to define an alternative system representation.

Definition 2.
A behavior B is controllable if for all w 1 , w 2 ∈ B, there exists ā t > 0 and a w ∈ B such that A system is uncontrollable if it is not controllable. In terms of the kernel representation, a system is controllable if and only if the matrix polynomial R(z) is left prime [15].
An alternative behavior representation for B ∈ L q can be obtained by introducing some auxiliary variablesw. The corresponding representation is a hybrid representation, and it has the form R(σ)w = M (σ)w for a certain full column rank matrix polynomial operator M (σ). An image representation is a special hybrid representation where R(σ) is the identity.
Lemma 2. Given a controllable behavior B ∈ L q and a set of auxiliary variables w, there exists a matrix polynomial operator M (σ) such that We underline that image representations exist only for controllable behaviors. Moreover, by eliminating the auxiliary variablesw, we can always switch from an image to a kernel representation of the same system [1].
In the following sections, the matrix polynomial operators R(σ) and M (σ) and their rows and columns, respectively, will play an important role. Given B ∈ L q , its complexity is defined as the pair (m, n) = (number of inputs, order). B| L denotes the behavior B restricted to the interval [1, L], which means all the trajectories are truncated to the time L. Its dimension is given by the following formula [16] in terms of the system complexity Some important results link the behavior B| L with the image and the kernel of some special structured matrices.
Definition 4 (Hankel and Toeplitz matrices). Given a time series w(t) ∈ R q , the block-Hankel matrix with L block rows H L (w) is defined as Given a polynomial of degree ℓ, p = p 0 + p 1 z + · · · + p ℓ z ℓ , the Toeplitz matrix with L columns (for L ≥ ℓ + 1) is the following (L − ℓ) × L matrix: We can now state the connection between a (finite length) behavior and the block-Hankel matrix built from an observed trajectory w [17].
As a consequence of Lemma 3, we have that (7) allows the computation of the complexity of B directly from an observed trajectory w (that is, without writing down a system representation) by solving a linear system of equations. The key point is to equate the right-hand sides of (4) and (7): The system of equations is built by choosing two suitable values for L in (8).
We summarize how to perform such a computation in Algorithm 1.

Algorithm 1 Computation of system complexity from a trajectory
Proposition 1. Given a trajectory w ∈ B| T , Algorithm 1 computes the complexity of B.
Proof. The integer L is the maximum number such that the block-Hankel matrix H| L (w) has more columns than rows, hence both the Hankel matrices H| L (w), H| L−1 (w) are rectangular (i.e., (8) holds). The number of inputs and the order come from a double application of the formula (8) with the two Hankel matrices H| L (w) and H| L−1 (w): Assuming that m(B), n(B) are unknowns, we can compute them by solving a linear system of equations. The solution is possibly rounded since we are looking for integers numbers.
Remark 2. Lemma 3 is a classical result. It has been recently shown [18] that the assumptions of Lemma 3 (controllability and persistency of excitation of the input component) can be replaced by the condition on the rank (8).
One more result which links a finite-length behavior and a structured matrix, follows: Lemma 4. Given a behavior B expressed by its kernel representation B = ker R(σ), let R 1 (σ), . . . , R p (σ) be the rows of R(σ). Then we have . . .
where ℓ is the degree of R(σ) and T L (R 1 ), . . . , T L (R p ) are Toeplitz matrices with L block columns.
Proof. Given a time series w = w(1), . . . , w(T ) , assume w ∈ B| T . If R(σ) is the matrix polynomial operator of degree ℓ, which represents the behavior B, then for each row R i (σ) of R(σ) we have Equation (9), written in matrix form, leads to T L (R)w = 0, i.e., the elements of B| T are in the kernel of the generalized Toeplitz matrix. On the other side, if the coefficients generating the Toeplitz matrix are the ones of the matrix polynomials operator representing the behavior B, the elements in the kernel of such structured matrix belong to the behavior B| T . Observe that the degrees of the annihilators determine the size of the associated Toeplitz blocks.
Lemmas 3 and 4 are useful because they link behaviors, trajectories, representations, and structured matrices. Hence these results connect system theory using the behavioral approach, linear algebra, and matrix computation (this highlights the differences with the modules-based approach of some works mentioned in Section 1.1). Lemma 3, also known in the literature as fundamental lemma, provides the conditions to retrieve the originating behavior directly from an observed trajectory. The fact that we can construct finite-length trajectories using the image of a Hankel matrix is effectively used in data-driven control (see [19] and the references therein for a summary of the recent works on the topic). The interest in data-driven algorithms also motivated the extension of such a result to several classes of nonlinear systems [20,21,22].
On the other hand, the result in Lemma 4 allows the construction of finitelength trajectories starting from the system kernel representation. This approach is not as direct as the previous one, but computing the kernel of a numerical matrix is usually easier than computing the kernel of a matrix polynomial operator.
Remark 3. This work deals with deterministic systems and exact (noiseless) data. If the data are affected by noise, the Hankel matrix H L (w) is full rank independently on the value of L; a possible noise removal approach can be preprocessing the data via Hankel low-rank approximation [23,24] to satisfy the rank condition in (8).

Addition and intersection of behaviors
We define the operations of addition and intersection of two linear timeinvariant behaviors by looking at the relation between their dimensions and the ones of the original systems. Then, we state how to compute the representations of the addition and intersection systems starting from the ones of the original systems. These results can be naturally extended to more than two behaviors. Definition 6. Given two behaviors, A and B, with the same number of variables, their sum is naturally defined as the set of the sums of the elements of A and B, while their intersection is defined as the set of elements that belong to both A and B: If A and B are linear and time-invariant, their sum and intersection are linear and time-invariant too.
The following result states the link between the starting systems' dimensions and the ones of their sum and intersection.
Proof. Equation (11)  The systems A + B and A ∩ B admit kernel representations, which can be expressed in terms of the kernel representations of A and B. If the systems are controllable, they admit also image representations, so that in this case we can also characterize their image representations in terms of the image representations of A and B. Theorem 1. Let A, B be controllable behaviors having the same number of variables. Let P a , P b be the image representations of A, B, respectively, and R a , R b be their kernel representations. The following hold true: 1. An image representation of the sum is given by the union of generators:

A kernel representation of the intersection is given by the union of anni-
hilators: Proof. First of all, we underline that the controllability assumption is needed only to guarantee the existence of the image representations; hence it is not required for the second point. Consider the first point. If w a = P a (σ)ℓ a is an element of A, and Hence, the image representation of the sum is obtained by stacking next to each other the two image representations. For the second point, let 0 = R a (σ)z = R b (σ)z for a certain z ∈ A∩B. Then Therefore it follows the expression for the kernel representation of the intersection of the two behaviors (13).
Theorem 1 shows a duality between the addition and the intersection of behaviors and the corresponding representations as union of generators and intersection of annihilators. This means that the representations of the sum and intersection systems can be computed by opposite operations by switching annihilators with generators, union with intersection, row-wise with column-wise operations on matrices. We expect that these relations still hold true by reversing the computations of the sum and intersection systems as intersection of annihilators and generators, respectively, as shown in Table 1. But the computation of the intersection of annihilators or generators needs to be implemented in an algorithm.

Intersection of annihilators and generators
The union of generators and annihilators is easy to be computed. But we may need to find the sum or intersection systems representations by the intersection of annihilators or generators, respectively. We describe, in the following, the problem of intersection of annihilators. Then, we can apply the observed duality property to get the dual computational algorithm for the intersection of generators.

Problem 1. Given minimal kernel representations of the two behaviors
By solving this problem, we obtain a computational method for (R a , R b ) → R + , which is a direct way to get the kernel representation of the sum from the kernel representations of the starting systems. The computational algorithm for this problem is Algorithm 2, whose correctness is proved below.

Algorithm 2 Sum of two behaviors by intersection of annihilators
Require: R a , R b kernel representations of the starting systems Ensure: R + kernel representation of the sum system Compute the number of outputs of the sum p(B + ) = q − m(B + ) Build the Sylvester matrix by choosing L such that the left kernel of S L has dimension p(B + ) Compute the left kernel basis Z a −Z b of the Sylvester matrix S L Define R + as Proposition 2. Algorithm 2 computes the kernel representation R + of the sum of the two starting behaviors A + B.
Proof. The number of outputs of the sum system is needed to understand the number of rows of the sought representation, and it can be computed from the starting systems by (12). The Sylvester matrix has two Toeplitz blocks generated by the starting representations, and its dimension depends on a parameter L: if p(B + ) > 0, we can always choose L such that S L has a non-trivial left kernel of dimension p(B + ). Such a value of L can be computed from the dimensions of the starting representations; if the sum system is nontrivial (p(B + ) > 0), the size of S L (the difference between the number of rows and the number of columns) depends on the value L.
Consider then two trajectories a ∈ A ⊂ A + B and b ∈ B ⊂ A + B. If the Sylvester matrix S L has a nontrivial left kernel (of suitable dimension), the sought representation R + (σ) satisfies Observe that, by construction, Z a −Z b has (at least) p(B + ) rows. We can define the coefficients of R + (z) as Since R + is the kernel representation of both the behaviors A and B, it is also a representation of their sum A + B.
The computed kernel representation can, in general, be nonminimal depending on the choice of L in the Sylvester matrix and on the possible nonminimality of the starting representations.
Under the assumption p(B + ) > 0, possibly by changing the value of L, we can always find p(B + ) elements in the left kernel of the Sylvester matrix S. The algorithm always computes a nontrivial solution for the kernel representation. But, if p(B + ) = q − m(B + ) = 0, the kernel representation is trivial (the behavior has no annihilators different from the zero polynomial). In this case, the Sylvester matrix S always has more columns than rows for each value of L (or possibly it is square), so the left kernel is trivial (we assume the matrix polynomials R a and R b have no common factors, see [25]).
Next, we state the dual problem for the intersection of generators, where the controllability of the behaviors is assumed to guarantee the existence of the image representation.
The computational procedure to approach this problem is dual with respect to the one described in Algorithm 2, so we omit the corresponding algorithm and its proof of correctness.

Examples
We show here some examples that illustrate the results derived in the previous sections. We propose simple analytical computations dealing with the addition and intersection of some linear time-invariant systems, followed by numerical computations.

Scalar autonomous systems with simple poles
Consider two scalar autonomous linear time-invariant systems A, B defined by their minimal kernel representations where R a , R b are scalar polynomials of degree n a = n(A), n b = n(B), respectively. Assuming that all the poles are simple, the trajectories of A and B are the sum of damped exponentials: for some coefficients a i , b i , where exp z (t) := z t . By (15) and the definition of addition of behaviors (10), the trajectories of the sum of two linear timeinvariant systems with simple poles are still sums of damped exponentials, i.e., where the poles λ + are the union of the poles of A and B: λ + = λ(A + B) = λ(A) ∪ λ(B). The order n + is the number of distinct elements in λ(A) ∪ λ(B), that is n a + n b − n c (where n c is the number of common poles). By (15) and the definition of intersection of behaviors (10), also the trajectories of the intersection of two linear time-invariant behaviors contains sums of damped exponentials whose poles are the common poles of the two behaviors: The order of A ∩ B is the number of common poles n c between A and B.
The previous results are summarized in the following lemma.
Lemma 6. Let A and B be two scalar autonomous linear time-invariant behaviors with minimal kernel representations R a (σ), R b (σ), defined by the scalar polynomials R a (z), R b (z). The polynomials for the minimal kernel representations R + (z) of A + B and R ∩ (z) of A ∩ B are given by, respectively, the least common multiple and the greatest common divisor of R a (z) and R b (z).
Proof. The expression of the trajectories of the sum and intersection systems are given in (16) and (17), respectively. The systems poles of the two starting systems are the exponents λ ai , λ bi , which are also roots of the (scalar) polynomials R a (z) and R b (z), respectively. From (16), we can see that the poles of the sum system are the union of λ ai and λ bi , hence the roots of R + (z) are the union (without repetitions) of the roots of R a (z) and R b (z). Similarly, from (17) follows that the roots of R ∩ (z) are the intersection of the roots of R a (z) and R b (z). Since all the poles are simple, the thesis follows.
Remark 5. The result of Lemma 6 can be naturally extended to the case of Multi-Input Multi-Output systems by replacing scalar with matrix polynomials.
The fact that λ + is the union of the poles of the two systems can also be checked in Matlab. Two (random) systems can be generated by the function drss, once we fix the number of inputs, outputs and the orders (in the inputoutput setting, the sum is well-defined only for systems with the same number of inputs and outputs). pole(sys1) % poles first system pole(sys2) % poles second system pole(sys1 + sys2) % poles of the sum While the addition of systems can be easily obtained by the sum, to compute the intersection of two systems we should use some suitable algorithm (only algorithms for the common dynamic estimation of scalar autonomous systems exist at the moment). Anyway, the Matlab function intersect can be called to check the presence of common poles among the poles of the two systems.
Simulation. We plan to show a simulation example which reproduces the results of Lemma 6. The literate programming style [26,27] is used in the following. This makes the experiment reproducible (the reader can copy and paste the code in the Matlab command line, but the numerical values will be different due to the randomness of the data). In the following code chunks we use the function blkhank(p, l, T) to build block-Hankel matrices with l rows and T columns from a time series p. It can be downloaded from the slra toolbox [24].
The first step is to generate two scalar autonomous systems with a given set of poles and the corresponding trajectories. For convenience, the systems are represented in state-space form. y1 = initial(sys1, randn(3, 1), 20); %trajectory first system y2 = initial(sys2, randn(3, 1), 20); %trajectory second system A trajectory of the sum is given simply by adding the two trajectories of the original systems. We compute, then, the kernel representations of sys1, sys2 and sys1+sys2 by computing the left null spaces of the Hankel matrices H 4 (y1), H 4 (y_2) and H 6 (y1+y2), respectively, using the function blkhank. The lag ℓ for each system is its number of poles (or possibly it can be computed as ℓ = n/p, where n, p are the outputs of Algorithm 1. We can finally check that the roots of the three polynomials are exactly the expected poles. Observe that the presence of the common pole in the starting systems was reflected in the polynomials degrees. The roots of the computed polynomials also reveal the non-trivial common factor between R1 and R2. Remark 6. If the coefficients of the given representations are inexact, the presence of common poles between the polynomials in the kernel representations could not be detected; they can be estimated by computing approximate common divisors, e.g., via the algorithms developed in [28] for scalar polynomials (in the SISO case) or in [25] for matrix polynomials (in the MIMO case).

A single-input single-output system and an autonomous system
Consider a single-input single-output system and an autonomous system with simple poles. The kernel representation of the first is given by a 1 × 2 matrix polynomial: The kernel representation of autonomous systems involves a square matrix polynomial R(σ) whose determinant is nonzero. The poles are then the roots of the determinant of R(σ). We consider the following kernel representation: .
By choosing an input / output partition of the set of variables 2 w = (u, y), the trajectories of the first system satisfy the equation where the star denotes the convolution product. We need to add to these trajectories the free response y a,f ∈ ker p a (σ), which are the trajectories corresponding to zero input. Hence, the trajectories of the system A have the general form Observe that both the free response y a,f as well as the impulse response h a are sums of damped exponential signals of the form na i=1 α i exp za,i (t) where n a is the degree of the polynomial p a , while z a,1 , . . . , z a,na are the roots of p a .
The trajectories of the system B satisfy the equation u p b (σ)y = 0. We see that the input can only be zero, so that the output is constrained to the free response, i.e., y b,f ∈ ker p b (σ). The trajectories of the system B have the general form w b = 0 y b,f . The free response y b,f is still a sum of damped exponentials signal whose exponents are the poles of the system B.
The sum A + B is a single-input single-output system whose trajectories have the form w + = w a + w b , and the poles λ(A + B) are the union of the poles λ(A) ∪ λ(B). But the poles λ(B) appear only in the free response and not in the convolution with the input. A kernel representation of the sum A + B is given by The kernel representation (18) shows that the sum of a single-input singleoutput system with an autonomous system is always uncontrollable because of the presence of the common factor p b (z) [1] (i.e., the matrix R + (z) is not left prime). The trajectories of the intersection A ∩ B should be of the form w a and w b at the same time. Hence, the input is constrained to be zero and the output contains only the free response y ∩,f , which should be in the kernels of both p a (z) and p b (z), i.e., in the kernel of their greatest common divisor. Hence, the intersection is an autonomous system whose kernel representation has the following expression: Simulation. As in the previous example, we illustrate the results on a numerical simulation. The first system is a Single-Input Single Output controllable system.
The second system is autonomous (that is, uncontrollable), as in the previous example. Observe that now we have two system variables, so to match the sizes of the system trajectories we have to add a zero input to the trajectory of the autonomous system. sys1 = drss(1, 1, 1); %first system sys2 = drss(1, 1, 0); %second system u1 = randn(20, 1); %random input y1 = lsim(sys1, u1); %output first system w1 = [u1 y1]; %trajectory first system y2 = initial(sys2, randn(1, 1), 19); %output second system w2 = [zeros(20, 1) y2]; %trajectory second system Then, we can compute the kernel representations starting from the system trajectories. We use the function blkhank to build the Hankel matrices H 2 (w1), H 2 (w2) and H 3 (w1+w2). About the autonomous system, we want a square matrix polynomial (as in (19)) but we will find more annihilators because of the presence of the zero rows corresponding to the input. Hence only a subset of the computed annihilators is selected. It holds true that p = conv(p2, p1) and q = conv(p2, q1) (up to a normalization), so we checked the expression (18) numerically.
Remark 7. We remark that the computation of a trajectory of the intersection between two behaviors is possible (up to the authors' knowledge) only for autonomous behaviors by using ad hoc algorithms; hence, it is intentionally skipped in the numerical simulations.

Conclusion
We studied the two basic operations of addition and intersection of linear time-invariant systems in the behavioral setting, showing that they are different from the classical definitions in the input-output setting. The proposed trajectory-based definition of sum and intersection allow to perform such computations directly from the data. Moreover, we saw how the resulting system representations depend on the representations of the starting systems, and we proposed algorithms for their computations based on structured matrices and polynomial computations. We summarize some of the main advantages due to the proposed definitions and results: 1. the intersection has been extended to open systems (systems with inputs); 2. the two operations can be performed directly using the system trajectories (observed data); the system representations, if needed for the problem, can be computed at a later stage; 3. we can sum and intersect systems with different numbers of inputs and outputs (but the same number of variables!).
The development and implementation of an algorithm to estimate the common dynamics among open systems can be object of future research.