Matchable--Observable Linear Models for Multivariable Identification ...

2015 IEEE 54th Annual Conference on Decision and Control (CDC) December 15-18, 2015. Osaka, Japan

Matchable–Observable Linear Models for Multivariable Identification: Structure Selection and Experimental Results Rodrigo Alvite Romano, Member, IEEE and Felipe Pait, Senior Member, IEEE, and Rafael Corsi Ferrão Abstract— Identification of linear time–invariant multivariable systems can best be understood as comprising three separate problems: selection of system model structure, filter design, and parameter estimation itself. In previous contributions we approached the first using matchable–observable models originally developed in the adaptive control literature, and used direct or derivative–free optimization to design filters. In this paper we show a simple and effective structure– selection method and demonstrate its accuracy, robustness and moderate computational demands using data from an industrial evaporator and experimental results with a twin rotor.

I. M ULTIVARIABLE L INEAR I DENTIFICATION Efficient and reliable techniques for identifying multivariable plants, which often present significant interactions between the multiple input–output channels, have attracted considerable research effort. Challenges include the development of robust and accurate algorithms capable of dealing with large scale systems, especially when the amount of available data is limited. Methods that have been employed successfully include subspace methods, data driven local coordinates, and expectation– maximization algorithms. In the data driven local coordinates (DDLC) method [1], [2], at each iteration the optimization is regularized, so that the search is performed in a subspace which is of minimal dimension. The expectation– maximization (EM) algorithm [3] computes ML estimates of fully parameterized state–space models without using gradient–based optimization. Subspace system identification [4] methods are effective in finding state–space MIMO models for several applications, without requiring specific structural decisions. The subspace approach may also be used to provide an initial estimate for nonlinear optimization algorithms that compute ML estimates [5], [6]. The quasi–canonical MIMO state–space representations proposed in [7] in the context of adaptive control have the properties that allow parameter estimation to be accomplished without the need for iterative or nonlinear optimization. The present paper continues the development (see [8], [9]) of a multivariable system identification algorithm based on a discrete–time version of these matchable–observable design models, which, in spite of their desirable features, have to the best of our knowledge not yet been employed in system identification. R. A. Romano and R. C. Ferrão are with Escola de Engenharia Mauá, Instituto Mauá de Tecnologia, São Caetano do Sul, SP, Brazil; e-mail: [email protected], [email protected]. F. Pait is with Escola Politécnica da Universidade de São Paulo, SP, Brazil; e-mail: [email protected].

978-1-4799-7886-1/15/$31.00 ©2015 IEEE

A crucial point in using these representations is the choice of lists of observability indices. Structural identification for MIMO dynamic systems parameterized through quasi– canonical forms is the subject of a literature developed since the 1970s which focuses on methods for estimating the Kronecker indices [10], [11]. In the pioneering approach of [12], determination of structural indices of an controllable canonical form was treated using singularity tests on input– output data matrices. The structural parameters of a canonical state–space representation can also be determined by testing the determinant of a certain matrix whose entries are the correlation functions of the output sequence [13]. In the sequel we give an overview of the rationale behind the identification scheme we advocate. Section III gives an abridged description of the matchable–observable linear identification models presented in [8]. Section IV proposes a model structure selection algorithm, and Section V describes applications to real data. II. T HE MOLI–ZOFT A PPROACH Identification of linear time–invariant multivariable input– output systems can best be understood as comprising three separate problems. First is selection of a structure containing system models capable of approximating the input–output behavior of the process to be identified. Second is the design of filters that separate process signals from disturbances and measurement noise, in a manner conducive to parameter estimation. Third is parameter estimation itself. In single–input single–output (SISO) identification, model selection consists simply in picking the order of the transfer function to be estimated. In MIMO identification, linearly parameterized state–space realizations may contain process models of varying McMillan degree, leading to pole–zero cancellations and other problematic phenomena. Some of these issues appear when using equation error methods based on autoregressive models such as ARX and its generalizations. It is convenient to consider a more fundamental invariant: the list of observability indices of the models under consideration. Rather than trying to estimate the “true” indices of the plant, we start with low–order models and augment the elements of the list of observability indices as needed to improve identification. Our approach to the varied signals collectively referred to as “noise,” employs the barycenter method for direct or derivative–free optimization, as presented in [9]. Finally, parameter estimation itself reduces to linear regression and can be approached by standard least–squares methods [8]. The present paper focuses on the structure selection mechanism.

3391

The overall procedure we advocate combines the 3 phases hierarchically as follows. We start with a low–order model, with a given set of observability indices. For each value of the filter parameters, the least–squares optimal values of the model parameters can be computed. The estimation procedure results in an index of merit, the user defined criterion J, which depends on the filter. We search for filter parameters via direct optimization of J, and the model parameters are determined as a consequence. We then increase the McMillan degree of the model by raising each observability index in the list. For each list of indices, the procedure of direct filter optimization, which in turn calls for the solution of a sequence of least–squares problems, is repeated. The one list of indices which results in the largest improvement of J is selected. The process of augmenting the list of observability indices one by one is repeated until either the desired performance is achieved or the model dimension reaches an user defined upper bound. We call the method MOLI–ZOFT, for Matchable Observable Linear Identification—Zero–order Oracle Filter Tuning.1 The procedure for picking the model structure by augmenting the list of observability indices and thus increasing the model’s McMillan degree is straightforward and does not lead to a combinatorial explosion of the number of candidate model structures. Moreover, direct optimization of the filters is computationally light and robust, and parameter estimation via linear least–squares can be done in closed–form for batch identification, or in a recursive manner, so that finite–time convergence is guaranteed by construction. III. M ATCHABLE – OBSERVABLE LINEAR IDENTIFICATION MODEL STRUCTURE

Any linear system of McMillan degree n with m inputs (uk ∈ Rm ) and p outputs (yk ∈ Rp ) can be represented by the following state–space model structure −1 xk+1 = A + D(θ) (I − G(θ)) C xk + B(θ)uk (1) yk = (I − G(θ))

−1

Cxk ,

(2)

whose coefficients are in turn parameterized via a parameter vector η. The observer state ξk is composed of delayed versions of the input–output data filtered by 1/α(q), where q is the shift operator. Thus, given η, the estimation of θ (which is mapped to CI and GI ) can be carried out through (3)–(4) using a linear regression. The roots of α(q) may be interpreted as filter poles, and tuned via derivative–free optimization. For a more detailed explanation we refer the reader to [7], [8], and to research to be published [15]. IV. E STIMATION ALGORITHMS We have argued that system identification can be profitably decomposed into three separate tasks. A previous contribution discussed parameter estimation itself, approached via proven LS methods [8], and another one [9] showed how filter tuning, the adjustment of the vector η containing the parameters of α, can be effected via direct optimization. The present paper shows recently obtained results concerning structure selection. A. Filter tuning Rather than treating all coefficients of α(q) as free parameters, we simplify the design and search for the prefilter’s dominant poles only, that is, α is parameterized in terms of discretized versions of second–order polynomials with cutoff frequency ωc and damping ζ, the elements of the 2–dimensional vector η. Then, filter tuning may be treated as a 2–dimensional optimization problem in a convex set. This will be accomplished by means of derivative–free direct optimization. Consider a set Ω = {η1 , . . . , ηnη } of candidate values of η, and let J (η, Z) be a functional which quantifies the performance of the values η ∈ Ω given a data set Z = {y1 , u1 , . . . , yN , uN } of length N . We shall pick the values of the filter parameters according to the formula Pnη ηr e−µJ(ηr ,Z) ∗ , (5) η = Pr=1 nη −µJ(ηr ,Z) r=1 e

ξk+1 = AI ξk + DI yk + BI uk

(3)

where ηr denotes the rth element of Ω and µ ∈ R is an appropriately chosen positive constant. In formula (5), η ∗ represents the barycenter of the points ηr with weights e−µJ(ηr ,Z) . The intuition behind taking the barycenter as a substitute for the hard–to–compute value of η which minimizes J (η, Z) is that the points ηr which lead to low values of the cost criterion are given more weight than those that lead to high values. The complete procedure was presented in [9].

yk = CI (θ)ξk + GI (θ)yk ,

(4)

B. Model structure selection

where θ denotes the model parameters, I is an identity matrix of suitable dimensions, and (C, A) is a stable, observable and parameter–independent pair. The parameter matrix G(θ) ∈ Rp×p is strictly lower triangular. The other parameter matrices D(θ) and B(θ) take values in Rn×p and Rn×m , respectively. This structure is fully spelled out in [7]. It can be shown that the state–space realization

is input–output equivalent to the model structure (1)–(2).The construction of (3)–(4), which we shall call a Morse observer, is omitted from this paper for reasons of space. For our purposes here it is enough to consider that the eigenvalues of AI are all contained in a polynomial α(q), 1 The term zero–order oracle refers to the fact that the direct optimization algorithm employed for filter tuning is derivative–free: it makes use of the value of the function only [14].

Previous work described algorithms for a fixed list of the observability indices l. Here the determination of a suitable l is addressed. We employ the following strategy to select the list of observability indices. Rather than trying to determine the model internal structure for a certain dynamic order n, the structure selection procedure starts from the list l = {1, . . . , 1}, which corresponds to the simplest p-output model

3392

from a structural point of view (cases with outputs with observability index zero are of limited interest). Next, p new models are estimated using a set of linear regressions driven by the direct optimization method of [9], each one structured from a new list obtained by raising the index relating to one of the outputs. Then, one picks the list l which results in the best performance, according to a criterion J (·) established beforehand. The strategy is depicted in Figure 1. l = {1, 1, . . . , 1}

...

l = {n1 + 1, n2 , . . . , np }

l = {n1 , n2 , . . . , np + 1}

ηnη

η1 parameter estimation

...

J(·)

parameter estimation ... Filter tuning

direct optimization η∗ parameter estimation Filter tuning J(·)

pick the best model according to J(·) and update l with the corresponding list

Stop?

ˆ i is its where yi is the observed ith output sequence and y estimate, which is computed using only the deterministic part of each identified models. The operator var (·) denotes the variance of the argument. The closer VAFi is to 100%, the lower the energy of the modeling error in the ith output.

no

yes

Fig. 1.

The PO–MOESP (Past Output–Multivariable Output– Error State–sPace) subspace method [4] (implemented via the software described in [16]). • The expectation maximization (EM) algorithm developed in [3] and implemented in the University of Newcastle Identification Toolbox (UNIT) [6]. • Gradient–based search using the Data-Driven Local Coordinates (DDLC) parameterization [1]. The gradient–based approach is available in the widely used MATLABr System Identification ToolboxTM [17] through the pem function, which was used with the property Focus set to three different options: “Predicton”, “Stability” and “Simulation” (henceforth pred, stab and sim, respectively). The first incorporates np extra parameters to model the noise; the second forces the model to be stable; and the third employs an output–error method. The property nk was set to [1,1,1], in order to not consider the feedthrough term. The maximum number of iterations was limited to 40 for EM and gradient–based search methods. These iterative methods were initialized using a CVA–weighted subspace estimate achieved using the n4sid function as discussed in [17]. The PO–MOESP block-size parameter was set to 2n. To measure the identified models accuracy, consider the variance–accounted–for criterion (VAFi ) [16] î) var (yi − y VAFi = max 1 − , 0 × 100% , var (yi ) •

Model structure selection, filter tuning, and parameter estimation.

In order to keep the number of steps finite, an upper bound nmax to the McMillan degree is defined, and the model order is chosen over the interval between p ≤ n ≤ nmax as a compromise between complexity and performance. As the variance–accounted–for (VAF) is used in the present work as a performance measure, the lowest model order whose average VAF of the outputs is within 1% of the best attained one is selected. This simple heuristic circumvents the problem posed by the number of the candidate structures, which increases combinatorially with model degree. The price is that not all possible lists of observability indices are considered and one might end up with a more complex model than necessary. V. N UMERICAL AND EXPERIMENTAL RESULTS Data from a multivariable industrial process was used in V-A to test MOLI–ZOFT, including the structure selection algorithm, and compare it with well–established estimation methods. Some of the test results in Section V-A are contained in a paper under review [15], and are presented here for completeness and ease of reference. We compare our proposed algorithm with the following methods:

A. A real–data case study The MOLI–ZOFT algorithm was applied to the identification of a real four–stage industrial evaporator for reducing the water content of a product. The data set is available from the Daisy database [18]. The system inputs are feed flow (u1 ), vapor flow to the first evaporator stage (u2 ) and cooling water flow (u3 ). The three outputs are the dry matter content (y1 ), the flow (y2 ) and the temperature of the outcoming product (y3 ). Only 500 samples of the available data were used for estimation and another 500 for validation. The functional J (η) used for structure selection is an average VAFi value: p ˆ i (η)) 1X var (yi − y J (η) = min ,1 . p i=1 var (yi ) Model structure selection is performed using the procedure described in Section IV-B. A modified randomized version of the direct optimization method of [9] is used. The starting point of the procedure is the definition of a set Ω of candidate filter parameter vectors. At each step the barycenter is computed using seven values of the candidate filter parameters η: the barycenter computed in the previous step, three elements chosen at random from the original Ω, and three other random candidates centered on

3393

89

P O–MOE S P (+) ↓

DDL C si m (♦) ց

88 87 Average VAF

the barycenter, with normal distribution and 20% of variance of the sets ζ and ωc . In this case the direct optimization filter tuning is carried-out considering a set Ω of 18 candidates to η which are generated by ζ = {0.6, 0.7, 1.1} and ωc = {0.02π, 0.0355π, 0.0631π, 0.1122π, 0.1995π, 0.3548π}. Considering candidate structures up to n = 15 and starting from l = {1, 1, 1}, the following list sequence was obtained: {1, 2, 1}, {1, 3, 1}, {1, 4, 1}, {1, 5, 1}, {2, 5, 1}, {3, 5, 1},{3, 6, 1}, {4, 6, 1}, {4, 7, 1}, {4, 8, 1}, {5, 8, 1} and {6, 8, 1}. The average VAF computed using the validation data set for each of these models is presented in Figure 2. To evaluate the effectiveness of the search, the performance achieved by all possible lists of observability Ppindices (i.e., all combinations of p indices ni , such that i=1 ni = n) are shown using the symbols “◦”. For each McMillan degree, the result achieved with the best list l of all possibilities is highlighted by 4. The structure selection algorithm was able

86 85

↑ MOLI –Z OFT(∗)

84 83

տ E M(◦)

82

տ ⊓ D DLC p red (⊔)

81 6

9 Model order n

12

15

Fig. 3. Average VAF of the estimated outputs by models having orders 4 ≤ n ≤ 15. R u nti me to get th e mo d el wi th th e ch osen ord er

88

1

10

D D LC pr e d D D LC s im

Cumulative runtime (seconds)

Average VAF

86

84

82

80

0

10

3

6

9 Model order n

12

տ

տ

PO–MOESP

MOLI–ZOFT

MOL I –Z OFT: 1. 59s (n = 11) P O–MOE S P : 1. 14s (n = 8) E M: 4. 75s (n = 11) DDL C p red: 5. 05s (n = 10) DDL C si m: 3. 78s (n = 8)

10

Structure selection algorithm Best model when all possible lists are tested Arbitrary lists

ր

EM

−1

78

ց

ց

15 2

Fig. 2. Model accuracy for each candidate list of indices, for McMillan degree from 4 to 15.

to find the best list for most, although not all, McMillan degrees in this test. We compared MOLI–ZOFT to other multivariable state– space identification methods using the evaporator data. The average VAF provided by each method of order 4 to 15 is reported in Figure 3. The results with method DDLCstab are also not shown because they are almost identical to those with DDLCpred . Performance of PO–MOESP for n = 6 is not depicted because an unstable model was obtained. The average VAF attained by the proposed method is higher than 86% for virtually all models with n ≥ 5. This result is comparable to the accuracy achieved by PO– MOESP and DDLCsim , the algorithms which achieved the best performance for this dataset. We analyzed the computational effort needed to identify a real process in which the model order is not a priori known. Figure 4 presents the cumulative runtime necessary to estimate models from n = p to a certain order. The time spent to estimate the models with suitable orders is indicated by crosses. The test shows that MOLI–ZOFT is undoubtedly

4

6

8 10 Model order n

12

14

16

Fig. 4. Cumulative runtime (in seconds) spent to estimate models up to a specific order.

very competitive with the the other methods in terms of computing time. B. Twin rotor system The twin rotor system (TWR) is an academic workbench which is useful to evaluate control strategies, due to its highly coupled behavior. It is composed of two aerodynamic rotors (main and tail propellers) arranged at the ends of a free bar. These actuators enable the device to move around the horizontal (pitch angle) and vertical (yaw angle) axes.Therefore, the system inputs are the signals applied to the the main (u1 ) and tail rotor u2 . The pitch and yaw angle measurement are the system outputs, denoted as y1 and y2 , respectively. System identification has been applied in the “black-box” framework to model the TWR system (e.g., [19], [20]). The more usual approach consists in handling each input–output pair individually, that is, only one rotor is activated at a time. Furthermore, to the best of authors knowledge, there is

3394

no work concerning the application of black–box identified models to the multivariable feedback control. The plant is excited using binary signals (Figure 5) with an average switching time of 8s. The data is recorded using a sampling period of 0.1s because the system dominant modes are below 1Hz. In this case, we intend to evaluate if the

80

Average VAF

70

0.15

60

50

Main rotor

0.1 0.05

40 MOLI−ZOFT PO−MOESP DDLCPRED

0 −0.05

30

DDLCSIM

−0.1

EM 4

0 0.3

10

20

30

40

50

5

6 7 Model order n

60

8

9

Tail rotor

0.2

Fig. 6.

0.1

Average VAF for increasing McMillan degree.

0 −0.1 0.2

−0.2 20

30 Time (seconds)

40

50

60

1

10

pitch (y )

0.1 0

0 −0.1

TWR input excitation.

−0.2 0

10

20

30

40

50

60

real data simulation (MOLI−ZOFT)

0.3 0.2 2

MOLI–ZOFT algorithm is able to provide a system model which is suitable for multivariable control design, using a restricted amount of samples. Thus, both actuators are excited simultaneously in a 60s experiment, but only the first 15s of the test (150 samples) are used for parameter estimation, which is considerably lower than other works reported in the literature [19]. The parameter vector η considered in filter tuning is determined by combinations of ζ = {0.5; 0.625; 0.75} and ωc = {0.8π; 1.04π; 1.28π; 1.52π; 1.76π; 2π}. Unlike the evaporator identification, in this case lower values of ζ are assumed, due to the perception of behavior typical of low damping modes in the output data. Figure 6 presents the average VAF achieved using the proposed structure selection heuristic (Section IV-B). For comparison purposes, the VAF achieved by other methods are also shown. The list of observability indices obtained is: {1, 2}, {2, 2}, {3, 2}, {4, 2}, {4, 3}, {4, 4} and {4, 5}. For the MOLI–ZOFT algorithm, the obvious choice is l = {3, 2} for which VAF1 = 81.5% and VAF2 = 79%. On the other hand, the other methods suggest 6th order models, except for DDLCsim which is the worst alternative for this data set. Figure 7 illustrates the TWR response and outputs simulated by the model identified by the MOLI–ZOFT algorithm for the 5th order selected structure. The estimated model reproduces the plant behavior with reasonable fidelity, even after the first 150 samples of the experiment, which were used in the model estimation. The responses of the models obtained by the other methods are omitted, as they are quite similar to the reported one. The model identified using the MOLI–ZOFT algorithm (associated to the structure selection heuristic) was used to

yaw (y )

Fig. 5.

0.1 0 −0.1 −0.2 0

10

Fig. 7.

20

30 Time (seconds)

40

50

60

MOLI–ZOFT output validation.

design a state feedback controller with integral action to address three requirements: 1) improve the transient response damping, 2) minimize the coupling effect and 3) attain zero steady–state tracking error. Thus, the state feedback gains are determined so as to minimize X JC = zkT Qzk + uTk Ruk , k

where zk is the augmented state vector including the integral error states. The weighting matrices Q, R and S are adjusted through simulations, from initial values obtained using the Bryson method [21]. In a black–box approach, the identified model states do not have physical meaning, and consequently cannot be measured. Therefore, the state vector is reconstructed using a Luenberger observer, whose gain is obtained through the pole placement technique. Table I presents the plant model, closed–loop and observer poles mapped to continuous–time domain. The closed–loop response to set–point changes is presented in Figure 8. It can be seen that the controller is able to

3395

TABLE I I DENTIFIED MODEL , CLOSED – LOOP AND OBSERVER POLES . TWR model poles −0.11 ± 2.15 −0.236 ± 0.678 −1.14

Closed–loop poles −0.282 −0.317 −2.26 −1.166 ± 2.81 −2.52 ± 2.34

and should be reported in upcoming work. We are grateful to the anonymous reviewer for asking about this issue.

Observer poles −8.62 ± 9.4 −11.45 ± 7.19 −16.2

ACKNOWLEDGMENT R Romano was supported by Instituto Mauá de Tecnologia (IMT), and by Fundaça˜ o de Amparo a` Pesquisa do Estado de São Paulo (FAPESP), grant 2012/03719-3. R EFERENCES

track pitch and yaw set–point changes in almost 10s. In fact, small oscillations around the reference are due to sensor quantization noise, or even to changes in the plant behavior when it moves away from the nominal operating point. In addition, the oscillatory behavior verified in open–loop was significantly reduced.

set−point y1

pitch angle

0.35 0.3 0.25 0.2 0.15 20

30

40

50

60

70

80

90

0.8 set−point y2 yaw angle

0.7

0.6

0.5 20

30

40

Fig. 8.

50 60 Time (seconds)

70

80

90

Closed–loop response.

The system coupling may still be noticed in the oscillation of y2 , caused by the change in the pitch set–point at t = 25s. However, it should be stressed that the coupling presented by a decentralized control strategy [22] (i.e., with two PID controllers) is much more significant. VI. C ONCLUSION We continue the study of a novel multivariable identification method and report on procedures to find a suitable set of invariants, without the need to test all possible list of observability indices. The MOLI–ZOFT methodology splits the overall identification problem into three stages, namely, parameter estimation, filter tuning and model structure selection. It proved to be at least as effective and more robust than other multivariable identification algorithms described in the literature. The results reported on the TWR system indicate the applicability of the structure selection heuristic to the MOLI–ZOFT approach, because they were successfully applied to the design of a state feedback controller and a full–order observer using an identified model from a restricted length data set. Case studies with high–dimensional systems are being carried out

[1] T. McKelvey, A. Helmersson, and T. Ribarits, “Data driven local coordinates for multivariable linear systems and their application to system identification,” Automatica, vol. 40, no. 9, pp. 1629–1635, 2004. [2] A. Wills and B. Ninness, “On gradient-based search for multivariable system estimates,” IEEE Transactions on Automatic Control, vol. 53, no. 1, pp. 298–306, 2008. [3] S. Gibson and B. Ninness, “Robust maximum-likelihood estimation of multivariable dynamic systems,” Automatica, vol. 41, no. 10, pp. 1667–1682, 2005. [4] M. Verhaegen and V. Verdult, Filtering and system identification: a least squares approach. Cambridge University Press, 2007. [5] T. McKelvey, “Identification of state-space models from time and frequency data,” Ph.D. dissertation, Linköping University, Linköping, Sweden, 1995. [6] B. Ninness, A. Wills, and A. Mills, “UNIT: A freely available system identification toolbox,” Control Engineering Practice, vol. 21, no. 5, pp. 631–644, 2013. [7] A. S. Morse and F. M. Pait, “MIMO design models and internal regulators for cyclicly switched parameter-adaptive control systems,” IEEE Transactions on Automatic Control, vol. 39, no. 9, pp. 1809–1818, 1994. [8] R. A. Romano and F. Pait, “Linear multivariable identification using observable state space parametrizations,” in Proceedings of the 52nd IEEE Conference on Decision and Control, Firenze, 2013. [9] R. Romano and F. Pait, “Direct filter tuning and optimization in multivariable identification,” in Decision and Control (CDC), 2014 IEEE 53rd Annual Conference on, Dec 2014, pp. 1798–1803. [10] M. Gevers, “A personal view of the development of system identification,” IEEE Control Systems Magazine, vol. 26, no. 6, pp. 93–105, Dec 2006. [11] A. Królikowski, “Model structure selection in linear system identification — survey of methods with emphasis on the information theory approach,” Eindhoven University of Technology, Eindhoven, Tech. Rep. 82-E-126, 1982. [12] R. Guidorzi, “Canonical structures in the identification of multivariable systems,” Automatica, vol. 11, pp. 361–374, 1975. [13] H. El-Sherief, “Multivariable system structure and parameter identification using the correlation method,” Automatica, vol. 17, no. 3, pp. 541–544, 1981. [14] Y. Nesterov, Introductory Lectures on Convex Optimization: A Basic Course. Springer, 2004. [15] R. A. Romano and F. Pait, “Matchable–observable linear models and direct filter tuning: An approach to multivariable identification,” 2014, submitted to the IEEE Transactions on Automatic Control. [16] M. Verhaegen, V. Verdult, and N. Bergboer, Filtering and system identification: an introduction to using Matlab software, Delft University of Technology, 2007. TM [17] L. Ljung, System identification Toolbox User’s Guide, Version 7.4, The MathWorks, Inc., Natick, MA, 2011. [18] B. L. R. De Moor, “DaISy: Database for the identification of systems,” Department of Electrical Engineering, ESAT/STADIUS, KU Leuven. Leuven, Belgium. http://homes.esat.kuleuven.be/ smc/daisy/, Feb 2008. [19] S. M. Ahmad, A. J. Chipperfield, and M. O. Tokhi, “Parametric modelling and dynamic characterization of a two–degree–of–freedom twin–rotor multi–input multi–output system,” Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, vol. 215, pp. 63–78, 2001. [20] S. F. Toha and M. O. Tokhi, “PID and inverse–model–based control of a twin rotor system,” Robotica, vol. 29, pp. 929–938, 2011. [21] J. P. Hespanha, Linear Systems Theory. Princeton, New Jersey: Princeton University Press, 2009. [22] Twin rotor MIMO system 33-949S User Manual, Feedback Instruments Ltd., East Sussex, UK.

3396