Interval System Identification for MIMO ARX Models of Minimal Order
53rd IEEE Conference on Decision and Control December 15-17, 2014. Los Angeles, California, USA
Interval System Identification for MIMO ARX Models of Minimal Order Stefan Zaiser, Michael Buchholz, and Klaus Dietmayer Abstract— Focus of this paper is on system identification of models in AutoRegressive with eXogenous inputs form from data with unknown, but bounded measurement errors. Hereby, these errors in data as well as the resulting uncertainty in parameters are represented by intervals. The proposed method can be applied to linear, time invariant systems with multiple inputs and multiple outputs. The main contribution are algorithms to determine the minimal order of a discrete-time model description in ARX form with interval parameters. In addition, an approach for using multiple sequences of measurement data is introduced. Finally, the method is demonstrated and discussed on a simulation example.
I. I NTRODUCTION In many fields of engineering, mathematical models of the underlying system are required. These descriptions of the dynamics can be obtained via first-principles modeling or via black-box modeling, if measurements of the system are available. The latter, also called system identification, usually results in high quality models describing the input/output behavior. While these models are computational cheap, they do not allow any insight into internal system processes. In contrast to classical system identification, where uncertainties in measurement data are described by stochastic distributions, uncertainties are described as unknown, but bounded errors in this contribution. These errors are represented by intervals [1], i.e. an upper and lower bound enclosing the measured value. When using interval data for system identification, the resulting model parameters are intervals themselves, instead of fixed values in the classical approach. Such interval system models are e.g. used to generate guaranteed decisions in fault-detection and diagnosis, see e.g. [2], [3]. In literature, interval modeling approaches from measurement data concentrate on the identification of model parameters. Although sometimes called system identification, the order of the model is assumed to be known and the model parameters are enclosed using intervals (see e.g. [4] for an overview) or tighter set-membership descriptions like polytopes or ellipsoids. The survey paper [5] summarizes various aspects of the latter approaches. Other recent publications in this area are e.g. [6], [7] for SISO systems or [8] for MIMO systems. In [9], the order and an initial model are determined using classical point-real system identification, and then the model is extended to an interval transfer function model. This work is supported by the German “Federal Ministry of Education and Research” (BMBF) within the priority program “ICT2020” under project funding reference number “16N11869”. S. Zaiser, M. Buchholz, and K. Dietmayer are with the Institute of Measurement, Control and Microtechnology, Ulm University, 89081 Ulm, Germany. E-mail: [email protected]
978-1-4673-6088-3/14/$31.00 ©2014 IEEE
Recently, we presented the (to our knowledge) first method for interval black-box modeling, i.e. order and parameter identification, of LTI systems with multiple inputs and a single output (MISO) [10], [11]. It determines the actual order of the system and results in a discrete-time interval ARX model, which can be directly re-written as a state-space model in the MISO case. It is based on one single sequence of interval data of inputs and the output. In the MIMO case, different ARX models of possibly different orders exist due to internal couplings, which all exactly represent the correct input/output behaviour of the system. In this paper, an extension of the method from [10], [11] to MIMO systems is proposed, which is capable of identifying one of these possible ARX models with the minimal order from interval data. Once the order is determined by the proposed algorithm, the parameters of the model are identified by well-known interval methods from literature. Additionally, an extension of the interval system identification method for multiple (possibly shorter) data sequences is derived. This paper is organized as follows: Section II introduces the mathematical notation, fundamentals, and prerequisites. In Section III, the method from [10], [11] is recapitulated. Based on this method, Section IV presents the new MIMO interval identification method including an extension for multiple sequences of measurement data. In Section V, a simulation example is given. The paper closes with conclusions in Section VI. II. F UNDAMENTALS AND N OTATION This section introduces the fundamentals and the notation used in this paper. First, the basics of intervals are briefly introduced. Afterwards, the used notation of ARX models is given. Finally, the general prerequisites for interval system identification are presented, as well as some fundamentals of interval parameter identification from literature. A. Intervals A closed interval of real numbers [x] = [x ; x] := {x ∈ R | x ≤ x ≤ x } ∈ IR
(1)
describes a variable with unknown, but bounded uncertainties and is defined by its infimum x and its supremum x. Classical point-real values are a special case of intervals with x = x = x. Alternatively, an interval can be represented by midpoint and radius
1774
[x] = [xc − x∆ , xc + x∆ ]
(2a)
with x+x x−x mid [x] = xc = , rad [x] = x∆ = . (2b) 2 2 Interval extensions of vectors and matrices (denoted in boldface, e.g. [A], [b]) are obtained by replacing the pointreal entries by interval entries. Calculations with intervals are done independent of the distribution within the bounds using interval arithmetics, cf. e.g. [12]. Toolboxes like I NTLAB [13] allow validated interval calculations in M ATLAB. For the presented methods, the rank of an interval matrix [A] plays an important role. It is defined as the minimal rank of all matrices enclosed by the interval matrix: o n rank [A] := min r ∀A ∈ [A], r = rank(A) . (3) r
Row and column interval rank are defined analogously. As shown in [14], an interval matrix can be tested for full column rank by evaluating the inequality −1 T Ac A∆ < 1 , (4) ρ ATc Ac where ρ(·) denotes the spectral radius of a matrix.
The sensor distortions (i.e. measurement errors) ξ(k) are unknown, but bounded. No distribution, but only the lower and upper bounds ξ(k) and ξ(k) have to be known for each sample k. As a result, the true value of each sample of the output yi is guaranteed to be within the interval of measured data •
ytrue,i (k) ∈ [ymeas,i (k)−ξ∆,i (k) , ymeas,i (k)+ξ∆,i (k)], (6) and analogously for each input ui . If inputs are exactly known, this results in a radius of u∆ = 0. Furthermore, two requirements concerning the interval bounds of the measurements and the sampling time have to be fulfilled, which are introduced in Section III-B. D. Interval Parameter Identification If the order nARX is known, the remaining task of system identification reduces to identify the parameter matrices Ai and B i of the ARX model. As discussed in the introduction, various methods to compute an enclosure of these parameters can be found in literature, e.g. [4], [5], [7], [8]. One popular approach is to formulate the parameter estimation task in the form of an overdetermined interval linear system (ILS) [M ][P ] = [N ]
B. ARX Models The interval identification methods discussed in this paper are based on measurement data of inputs u ∈ Rnu and outputs y ∈ Rny of a linear time-invariant (LTI) system. The data is gathered at sampling points k Ts . The resulting discrete-time system model with sampling time Ts ∈ R+ is described in the AutoRegressive with eXogenous inputs (ARX) form and time (or sampling point) index k ∈ Z: y(k) = −
n ARX X i=1
Ai y(k − i) +
n ARX X i=1
B i u(k − i) + B 0 u(k). (5)
nARX is called the order of the ARX model in the following, assuming that AnARX 6= 0 holds, where 0 is the zero matrix of the appropriate dimensions. The direct feedthrough, if available, is represented by B 0 u(k). This definition easily extends to the interval case by replacing the real inputs, outputs, and parameters by intervals of appropriate dimensions. C. Prerequisites for Interval System Identification Persistent excitation of the inputs is one of the main common prerequisites for classical identification methods of non-autonomous systems [15]. From the collected data, only the controllable and observable part of the system can be identified. Furthermore, the sampling time Ts has to be chosen with respect to the system dynamics, and a sufficiently high number of sampling points is required to obtain accurate results. These basic requirements also hold for interval system identification. In addition, the following assumptions have to be fulfilled: • The real system generates the output sequences ytrue,i (k) according to the system dynamics and the sequences of input values utrue,i (k). Both sequences are possibly unknown due to measurement errors.
(7)
with the data matrix [M ] ∈ IRr×m , the output matrix [N ] ∈ IRr×n and the unknown parameter matrix [P ] ∈ IRm×n for r > m for the MIMO case, which reduces to the data vector [N ] ∈ IRr and the unknown parameter vector [P ] ∈ IRm The solution set of this ILS is given by n o [P ] := P M P = N for some M ∈ [M ], N ∈ [N ] . (8) All measured interval samples are compatible with this resulting set of parameters. In general, the exact solution set has a complicated shape, and its calculation is very costly. Therefore, a box enclosing this exact solution usually is calculated. [4] gives an overview on methods to calculate such solution set enclosures. In this paper, parameter identification is done using the following method, which is explained in more detail in [11]. The overdetermined parameter identification problem is divided into several square ILSs (SILS), and each square problem is solved individually. [12] gives an overview of methods to solve SILS, and most methods for overdetermined ILS can be applied to the square case as well. The intersection of all individual solutions yields an enclosure of the wanted solution set. Although this procedure requires a comparatively high computational effort, and methods yielding tighter enclosers than interval boxes exist, this method is chosen for the sake of simplicity in this paper, because it can be applied directly to the data matrices used for order determination. However, it could be easily exchanged by any other parameter identification method. III. BASIC I DEA OF I NTERVAL S YSTEM I DENTIFICATION This section recapitulates the basic idea of the underlying interval system identification method for MISO systems, which was originally published in [10], [11].
1775
uT (1) uT (2) .. .
y(1) y(2) .. .
uT (2) uT (3) .. .
uT (N − p) y(N − p) uT (N − p + 1)
|
y(2) y(3) .. .
· · · uT (p) · · · uT (p + 1) .. .
y(N − p + 1) · · · uT (N − 1) {z Z
A. Undisturbed Point-Real Data of a MISO system Assume a MISO LTI system without any measurement errors. The measured data (i.e. the true data in this case) from inputs and the single output is combined in data vectors T z(k) = u1 (k) . . . unu (k) y(k) (9) for each time step k. For a single sequence of N time steps, a block Hankel matrix can be built from the vectors z(k): T z (1) z T (2) ... z T (p) z T (2) z T (3) . . . z T (p+1) e = Z . (10) .. .. .. . . . zT (N −p) zT (N −p+1)
...
zT (N −1)
e are time series of inputs and the The block columns of Z output, shifted by one sample for each block column. The parameter p denotes the number of block columns, which has to be chosen greater than the unknown system order p > n. A sufficiently large data set (cf. Section II-C) is assumed, i.e. N − p > p. For systems with direct feedthrough, an additional block column only consisting of input signals has e resulting in an extended matrix to be added to Z, T u (p + 1) .. e (11) Z= . . Z T u (N ) In the rest of this paper, only the general case with direct feedthrough is regarded, which can be easily reduced to the case without direct feedthrough. Based on the matrix Z, the parameters Γ of an ARX model of order p with direct feedthrough can be estimated using Eq. (12), which is shown on the top of this page. Because an ARX model of order p > n is not of interest, a method to determine the unknown order n is required in a first step. In [10], [11], it is shown for MISO systems that this can be done by analyzing the column rank of Z: colrank(Z) = n + p · nu + nu .
(13)
Comparing this to the number of columns of Z cZ = p + p · nu + nu ,
(14)
it can be seen that only the first summand differs. If the parameter p is chosen greater than the system order n, the matrix Z is rank deficient in the absence of any errors. Thus, the system rank can be determined by building Z for p > n and then calculate its column rank.
bp T ap y(p) u (p + 1) y(p + 1) . T y(p + 1) u (p + 2) .. y(p + 2) = .. .. . .. . . b1 y(N − 1) uT (N ) y(N ) a1 | {z } } b Y 0 | {z } (12) Γ
B. Interval Data of a MISO System In the interval case, the matrix [Z ] is built from the interval time series [u] and [y ] analogously to Eq. (9), (10). The matrix Z of the (unknown) true values is always enclosed by [Z ] (cf. Section II-C), thus the interval column rank (cf. Section II-A) of [Z ] can not be greater than the column rank of Z. To avoid an unwanted reduction of the interval column rank compared to the case without uncertainty, two additional prerequisites have to be fulfilled by the interval data: • Firstly, if the interval bounds are chosen so large that some columns have an overlapping interval in each row, the interval column rank is reduced. In this case, a pointreal matrix with additional linear dependant columns is included in the interval matrix. • Secondly, the interval column rank can be reduced by choosing a too small sampling time in comparison to the system dynamics. In that case, the intervals of two consecutive sampling points of one output signal could overlap for all sample points due to the slow dynamics, thus an overlap would exist in the shifted columns of the same signal. To meet these requirements, which is assumed in the following, sensors with sufficiently precise measurements and a reasonable sampling time have to be chosen. Then, the column rank of [Z ] and Z are identical and the system order can be detected by rank inspection. For the interval case, the test for full column rank (4) is conducted for matrices [Z ] built with increasing p, until rank deficiency is detected, i.e. a computational expensive calculation of the interval rank can be avoided. Thus, the order is nARX = p − 1. IV. I NTERVAL MIMO S YSTEM I DENTIFICATION It is not possible to extend the MISO method directly to MIMO systems, as the different outputs can be coupled. These couplings allow the computation of some output values from a combination of other delayed output values and the inputs, which leads to a linear dependency. Therefore, when building the matrix [Z ] as in (11) from the vector [z(k)] with all outputs [y1 (k)] . . . [yny (k)], the matrix can be rank deficient. This prevents the determination of a valid MIMO ARX model order by this naive approach. Additionally, MIMO ARX models are not unique in general. Therefore, a more sophisticated way to obtain the minimal order of a MIMO ARX model is proposed and discussed in the following. Additionally, it is shown how multiple data sequences can be used for interval system identification.
1776
A. Basic Order Determination of MIMO Systems
global: pmax ; eq ]; global: ∀q = 0, . . . , pmax : [U e q,i ]; global: ∀i = 1, . . . , ny , ∀q = 1, . . . , pmax : [y
The basic idea of the proposed algorithm is to determine combinations of shifted input and output sequences which lead to a final data matrix [Z ] with full column rank. Due to the assumption that the inputs are persistently excited (cf. Section II-C), only the combinations of shifted output signals have to be investigated. The combinations are investigated up to a user-defined maximum delay pmax , which should be greater than nARX . To reduce the number of investigations, the algorithm works recursively by excluding impossible combinations. It requires shifted signal sequences of all inputs and each of the outputs1 . For the inputs, sequence matrices are built for shifts q = 0, . . . , pmax : [uT (pmax − q + 1)] [uT (p − q + 2)] max e . [Uq ] = (15a) .. . [uT (N − q)] For each of the outputs (i = 1, . . . , ny ), shifted vectors with q = 1, . . . , pmax are set up: [yi (pmax − q + 1)] [y (p − q + 2)] i max . [y e q,i ] = (15b) .. . [yi (N − q)] e p ] represents the Within the algorithm, the matrix [Z temporary data matrix under investigation up to a maximum delay of p, and Ωp ∈ {0, 1}p×ny represents the combination e p ]. Each row q of the delayed output signals included in [Z of this matrix corresponds to the according time delay of the outputs and each column i to the respective output. Ωq,i = 1 means that yi delayed for q steps is selected in the current e p ]. e q,i ] is included in [Z combination, i.e. [y The algorithm is initially started with p = 1 and Ω0 = [ ]. Each cycle is based on the previously selected signals Ωp−1 . Firstly, the matrix Ωp−1 is extended by each of the 2ny −1 possible combinations ω j ∈ {0, 1}ny , ω j 6= 0,
j = 1, . . . , 2ny −1
(16)
e p,i ], i = 1, . . . , ny , of the output columns with delay p [y yielding the matrices T Ωp,j = ΩTp−1 ω j . (17) e 0 ], . . . , [U e p ] and the e p,j ] is built from [U Each possible [Z e p,j ] selected output columns in Ωp,j . Then, each matrix [Z is tested for full column rank using criterion (4). For each j e p,j ], the next cycle of the with full column rank of [Z algorithm is started with Ωp,j as Ωp−1 and p = p + 1. e p,j ] has full column rank, all If none of the matrices [Z e p,i ] can be expressed as combinations of output columns [y a linear combination of the output columns already included in the temporary data matrix Ωp−1 and the inputs, i.e. no 1 Note that with analogous point-real definitions of Eq. (15), the matrix Z e0 . ep y e1 y in Eq. (12) can be written as Z = U e p,1 · · · U e 1,1 U
function [Ωpmin , pmin ] = ord det(p, Ωp−1 ) { if p == pmax return Ωpmin =NaN, pmin =NaN; %invalid branch for j = 1, . . . , 2ny −1 { Ωp,j = Ωp−1 ; ω Tj ; %add j-th combination of outputs e p,j ] = [U e 0 ], . . . , [U ep ], all [y [Z e q,i ] selected in Ωp,j ; e p,j ] == full %use criterion (4) if colrank [Z [Ωpmin ,j , pmin,j ] = ord det(p+1,Ωp,j ); }; e p,j ] 6= full if ∀j = 1, . . . , 2ny −1 : colrank [Z return Ωpmin = Ωp−1 , pmin = p−1; j = arg(minj (pmin,j )); if isunique(j) %test if minimal solution is unique return Ωpmin = Ωpmin ,j , pmin = pmin,j ; c =user criterion Ωpmin ,j ; %select one of the minimal solutions return Ωpmin = Ωpmin ,c , pmin = pmin,c ; }; Alg. 1.
Pseudo-code of the basic recursive function
additional information on the system is added by the output columns of this delay p. Any further delays of any output(s) in this branch can also be iteratively substituted by the respective linear functions. Therefore, the algorithm terminates the current cycle and returns p − 1 as the wanted order for this branch together with the corresponding matrix Ωp−1 to the superordinate cycle. The recursion also terminates in one branch if no minimal order p < pmax could be determined and this branch is marked as invalid. If different combinations of output columns yield a solution in one cycle, the solution with the lowest order is chosen for this branch and returned to the superordinate cycle. In general, multiple solutions with this minimum order occur, because MIMO ARX models are not unique. This allows a further optimization of the selected model structure. Here, a criterion to maximize the number of zeros in Ωp is used, because each zero entry leads to ny point-real zeros instead of interval parameters to be estimated in the resulting ARX model, thus reducing the uncertainty of the model. If still multiple solutions exist, these are equivalent for our application, and any of them can be chosen. After all recursive cycles are terminated, the overall minimal order pmin is found together with the corresponding matrix Ωpmin . If all branches are marked as invalid, the algorithm returns an error. This can happen if pmax was chosen too small or if assumptions are violated. The recursive function of the algorithm is summarized in Alg. 1. Thereby, the variables marked as global are read-only values available to each cycle of the recursion. B. Optimized Order Determination of MIMO Systems The number of recursions of the algorithm can be further reduced, which is especially helpful for larger numbers of outputs. To do so, the possible combinations ω j of the output e p,i ] have to be processed in a specific columns with delay p [y order within the algorithm. First, the combinations with only
1777
e p,i ] one selected signal are tested. If one single column [y leads to rank deficiency, the true data in this column can be expressed as a linear combination of the true data contained in the outputs selected in Ωp−1 and the true input data by some linear function. Due to this relationship, the rank of matrices in this branch, i.e. based on the same selection of signals Ωp−1 , containing this column will always be deficient. Thus, all other combinations including this column can be removed from the set of possible combinations. Combinations based on the same Ωp−1 without this column of delay p, but with higher delays of the same output occur in other subbranches and do not have to be addressed in this step. If a combination leads to a data matrix of full rank, the next recursion cycle is processed as in the basic algorithm. After all combinations with only one column have been tested, the same approach can be followed for all remaining combinations of two columns. Again, all further combinations including these two columns can be removed if a combination leads to a rank deficient data matrix. This is repeated for all remaining combinations with an increasing number of columns, until all combinations either have been tested or removed from the set of possible combinations for the current delay p. Additionally, the recursion depth (i.e. the maximum delay) e p,j ] can be reduced during runtime. If none of the matrices [Z has full column rank, and p − 1 is returned as minimum order for that branch, the maximum delay can be reduced to pmax = p. Other solutions with the same minimal order p − 1 would be ignored if pmax is reduced to p − 1, not allowing for further optimization of the model structure. The result of this optimized algorithm equals the result of the basic algorithm from the previous subsection, but the number of recursions will be reduced for many problems. C. Discussion of the Order Determination Algorithms Both recursive algorithms always identify the minimal order nARX of the ARX model if pmax > nARX and all other prerequisites are met (cf. Section II-C and III-B). In the basic algorithm, many combinations of possible model configurations are recursively investigated, resulting in a high number of column rank checks. However, because the criterion (4) is used instead of an interval rank calculation, this check is computationally cheap, and the number of checks can be limited by selecting an appropriate termination limit pmax . The maximal number of rank checks is limited to (2ny − 1)pmax , however the actual number is far below this theoretical limit for most applications. The optimized algorithm even further reduces this number, e.g. (2ny −1 − 1) combinations can be removed if one single column combination leads to a rank deficient data matrix. D. MIMO Interval Parameter Identification After the recursive algorithm is successfully completed, the minimal order of the ARX representation of the system as well as the corresponding required signals are known. With this information, the parameters of the model can be identified from the interval data with any method from literature
(cf. Section II-D). The choice of the parameter identification method has significant influence on the overestimation of the parameters, but is not the focus of this paper. Here, the interval parameter identification method from [11], which was summarized in Section II-D, is used for sake of simplicity. For this purpose, a new data matrix [Z ] is constructed based on sequence matrices of all inputs and sequence vectors of the selected outputs in Ωpmin , which are built as in Eq. (15), now using the determined order pmin instead of pmax to make use of all data points. The corresponding unknown parameters to these input and output signals, which are subsumed in the matrix [Γ], are estimated using Eq. (12), where [Y ] is now a matrix with one column for each output signal. All other parameters of the ARX model are set to point-real zero, because the corresponding signals are not required within the model. E. Multiple Data Sequences If multiple data sets should be used for identification instead of one single sequence, these data sets cannot be simply concatenated due to the dynamics of the system. Instead, the vectors (15) have to be built for each sequence individually, as well as the matrix [Y ] for parameter identification. Afterwards, these shifted vectors, and thus the matrices [Z ] for each sequence, can be stacked. This assures that for the solution of the interval linear systems only data points from the same sequence are handled in one equation, i.e. one row of the matrices [Z ] and [Y ]. This extension can be applied analogously to the existing method from [10], [11]. V. E XAMPLE To demonstrate the proposed interval identification method for MIMO systems, it was implemented using the Intlab package [13] in M ATLAB. In the following, a second order system with two outputs, three inputs and direct feedthrough is regarded. The parameters of this ARX system (5) are − 0.5 0 0.2 0 , (18a) , A2 = A1 = 1 0 0 0.3 − 1 −1 0 0 0 1 B1 = , B2 = , (18b) −2 1 0 1 0 0 0 1 0 B0 = . (18c) 0 0 0 To generate a set of identification data of the true (undisturbed) system, this system is simulated to obtain 4 sequences of 2 500 sampling points each. Then, uniformly distributed random errors with a maximum amplitude of 1‰ of the respective maximum signal amplitude is superimposed on the simulated signals (inputs and outputs) to reproduce measurement errors. The bounds of the interval data sets are set to the same values, i.e. all prerequisites are fulfilled. The proposed interval system identification method results in an ARX model of the correct minimal order nARX = 2 with the interval parameter matrices [ 0.17 ; 0.23 ] [ −0.02 ; 0.02 ] [A 1 ] = , (19a) [ −0.04 ; 0.03 ] [ 0.28 ; 0.32 ]
1778
1
This is an enclosure of the true parameters (18). The pointreal zeros in the second column of the matrix [A2 ] result from the finally used combination of output signals 1 1 Ω2 = . (20) 1 0 This indicates that the second output delayed for two steps has no unique influence on both outputs and thus is not used for identification. Nevertheless, it could be added without causing rank deficiency, because the alternative solution with the same ARX order 1 1 Ω2,rejected = (21) 1 1 was rejected by the selection criterion to maximize the number of point-real zeros (cf. Section IV). Initialized with pmax = 10, i.e. (2ny − 1)pmax = 59 049, both algorithms had a maximum recursion depth of four, and 90 column rank checks (4) were performed by the basic algorithm, 63 by the optimized algorithm. To analyze the behavior of the set of systems included in the interval model, a Monte-Carlo analysis to calculate the eigenvalues was carried out with 100 000 point-real models randomly drawn from the interval model. Fig. 1 shows the true eigenvalues marked by blue asterisks. The convex hulls of the eigenvalues from the randomly drawn systems are shown as filled red areas. One eigenvalue of the interval model is point-real and located in the origin due to the pointreal zeros in [A2 ]. Thus, all randomly drawn models also have an eigenvalue exactly in the origin. Another eigenvalue is always real, but varies on the real axis. Together with the zoomed plot, it can be seen that the true system dynamics is included in the identified interval model, and the latter describes a set of systems with similar dynamics. Thus, the identified model is a valid description for the true system. VI. C ONCLUSIONS In this paper, a method for the identification of MIMO LTI systems from interval measurement data was presented. The main contribution are the determination algorithms for the minimal order of an ARX model of the system. It was also shown how multiple data sequences can be used for this class of interval system identification methods. The recursive order determination algorithms also select a minimum number of required delayed output signals, resulting in a maximum number of point-real zero entries
−0.67
0.5 Im
[ −0.51 ; −0.49 ] 0 [A 2 ] = , (19b) [ 0.99 ; 1.01 ] 0 [ −0.03 ; 0.03 ] [ 0.97 ; 1.03 ] [ −0.03 ; 0.02 ] [B 0 ] = , [ −0.03 ; 0.03 ] [ −0.04 ; 0.04 ] [ −0.03 ; 0.03 ] (19c) [ −1.03 ; −0.97 ] [ −1.04 ; −0.96 ] [ −0.02 ; 0.02 ] [B 1 ] = , [ −2.03 ; −1.97 ] [ 0.96 ; 1.04 ] [ −0.03 ; 0.03 ] (19d) [ −0.07 ; 0.07 ] [ −0.04 ; 0.04 ] [ 0.98 ; 1.02 ] [B 2 ] = . [ 0.92 ; 1.08 ] [ −0.05 ; 0.05 ] [ −0.03 ; 0.03 ] (19e)
−0.73 0.07
0 −0.5 −1
0.13
true system rand. drawn systems −1
Fig. 1.
−0.5
0 Re
0.5
1
Eigenvalues of true and identified systems
in the parameter matrices of the model, thus minimizing the uncertainty in the model. Besides a basic algorithm, an optimized version was presented which further decreases the number of recursions. With the results of these algorithms, an interval parameter identification method from literature is used to determine the non-zero model parameters. The overall method was demonstrated on a simulation example, showing the low number of recursions in practical use. R EFERENCES [1] R. E. Moore, Interval analysis. Englewood Cliffs: Prentice-Hall, 1966. [2] P. Planchon, Guaranteed Diagnosis of Uncertain Linear Systems Using State Set Observation. Berlin: Logos-Verlag, 2007. [3] F. Wolff, Konsistenzbasierte Fehlerdiagnose nichtlinearer Systeme mittels Zustandsmengenbeobachtung, ser. Schriften des Instituts f¨ur Regelungs- und Steuerungssysteme, Karlsruher Institut f¨ur Technologie. Karlsruhe: KIT Scientific Publishing, 2010, vol. 09. [4] J. Hor´acˇ ek and M. Hlad´ık, “Computing Enclosures of Overdetermined Interval Linear Systems,” Reliable Computing, vol. 19, no. 1, pp. 142– 155, 2013. [5] M. Milanese and A. Vicino, “Optimal estimation theory for dynamic systems with set membership uncertainty: An overview,” Automatica, vol. 27, no. 6, pp. 997–1009, 1991. [6] V. Cerone, D. Piga, and D. Regruto, “Enforcing stability constraints in set-membership identification of linear dynamic systems,” Automatica, vol. 47, no. 11, pp. 2488–2494, 2011. [7] ——, “Set-Membership Error-in-Variables Identification Through Convex Relaxation Techniques,” IEEE Transactions on Automatic Control, vol. 57, no. 2, pp. 517–522, 2012. [8] M. Pouliquen, E. Pigeon, and O. Gehan, “Output error identification for multi-input multi-output systems with bounded disturbances,” in 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC), 2011, pp. 7200–7205. [9] L. H. Keel, J. S. Lew, and S. Bhattacharyya, “System identification using interval dynamic models,” in American Control Conference, vol. 2, 1994, pp. 1537–1542. [10] S. Zaiser, M. Buchholz, and K. Dietmayer, “Black-Box-Modellierung mit unsicheren Parametern aus Messdaten mit unbekannten, aber beschr¨ankten Fehlern,” in Tagungsband GMA-Fachausschuss 1.30 ”Modellbildung, Identifikation und Simulation in der Automatisierungstechnik”, O. Sawodny and J. Adamy, Eds., VDI/VDE-GMA. Darmstadt: Technische Universit¨at Darmstadt, Institut f¨ur Automatisierungstechnik und Mechatronik, 2013. [11] ——, “Black-box modeling with uncertain parameters from measurement data with unknown, but bounded errors,” at - Automatisierungstechnik, vol. 62, no. 9, pp. 607–618, 2014. [12] R. E. Moore, R. B. Kearfott, and M. J. Cloud, Introduction to Interval Analysis. Society for Industrial and Applied Mathematics, 2009. [13] S. Rump, “INTLAB - INTerval LABoratory,” in Developments in Reliable Computing, T. Csendes, Ed. Dordrecht: Kluwer Academic Publishers, 1999, pp. 77–104. [14] J. Rohn, “A Handbook of Results on Interval Linear Problems.” ˇ Prague, Czech Republic, Tech. Institute of Computer Science AS CR, Rep. V-1163, 2012. [Online]. Available: http://hdl.handle.net/11104/ 0212095 [15] L. Ljung, System Identification: Theory for the User, 2nd ed. Upper Saddle River: Prentice Hall, 1999.
1779
Interval System Identification for MIMO ARX Models of Minimal Order Stefan Zaiser, Michael Buchholz, and Klaus Dietmayer Abstract— Focus of this paper is on system identification of models in AutoRegressive with eXogenous inputs form from data with unknown, but bounded measurement errors. Hereby, these errors in data as well as the resulting uncertainty in parameters are represented by intervals. The proposed method can be applied to linear, time invariant systems with multiple inputs and multiple outputs. The main contribution are algorithms to determine the minimal order of a discrete-time model description in ARX form with interval parameters. In addition, an approach for using multiple sequences of measurement data is introduced. Finally, the method is demonstrated and discussed on a simulation example.
I. I NTRODUCTION In many fields of engineering, mathematical models of the underlying system are required. These descriptions of the dynamics can be obtained via first-principles modeling or via black-box modeling, if measurements of the system are available. The latter, also called system identification, usually results in high quality models describing the input/output behavior. While these models are computational cheap, they do not allow any insight into internal system processes. In contrast to classical system identification, where uncertainties in measurement data are described by stochastic distributions, uncertainties are described as unknown, but bounded errors in this contribution. These errors are represented by intervals [1], i.e. an upper and lower bound enclosing the measured value. When using interval data for system identification, the resulting model parameters are intervals themselves, instead of fixed values in the classical approach. Such interval system models are e.g. used to generate guaranteed decisions in fault-detection and diagnosis, see e.g. [2], [3]. In literature, interval modeling approaches from measurement data concentrate on the identification of model parameters. Although sometimes called system identification, the order of the model is assumed to be known and the model parameters are enclosed using intervals (see e.g. [4] for an overview) or tighter set-membership descriptions like polytopes or ellipsoids. The survey paper [5] summarizes various aspects of the latter approaches. Other recent publications in this area are e.g. [6], [7] for SISO systems or [8] for MIMO systems. In [9], the order and an initial model are determined using classical point-real system identification, and then the model is extended to an interval transfer function model. This work is supported by the German “Federal Ministry of Education and Research” (BMBF) within the priority program “ICT2020” under project funding reference number “16N11869”. S. Zaiser, M. Buchholz, and K. Dietmayer are with the Institute of Measurement, Control and Microtechnology, Ulm University, 89081 Ulm, Germany. E-mail: [email protected]
978-1-4673-6088-3/14/$31.00 ©2014 IEEE
Recently, we presented the (to our knowledge) first method for interval black-box modeling, i.e. order and parameter identification, of LTI systems with multiple inputs and a single output (MISO) [10], [11]. It determines the actual order of the system and results in a discrete-time interval ARX model, which can be directly re-written as a state-space model in the MISO case. It is based on one single sequence of interval data of inputs and the output. In the MIMO case, different ARX models of possibly different orders exist due to internal couplings, which all exactly represent the correct input/output behaviour of the system. In this paper, an extension of the method from [10], [11] to MIMO systems is proposed, which is capable of identifying one of these possible ARX models with the minimal order from interval data. Once the order is determined by the proposed algorithm, the parameters of the model are identified by well-known interval methods from literature. Additionally, an extension of the interval system identification method for multiple (possibly shorter) data sequences is derived. This paper is organized as follows: Section II introduces the mathematical notation, fundamentals, and prerequisites. In Section III, the method from [10], [11] is recapitulated. Based on this method, Section IV presents the new MIMO interval identification method including an extension for multiple sequences of measurement data. In Section V, a simulation example is given. The paper closes with conclusions in Section VI. II. F UNDAMENTALS AND N OTATION This section introduces the fundamentals and the notation used in this paper. First, the basics of intervals are briefly introduced. Afterwards, the used notation of ARX models is given. Finally, the general prerequisites for interval system identification are presented, as well as some fundamentals of interval parameter identification from literature. A. Intervals A closed interval of real numbers [x] = [x ; x] := {x ∈ R | x ≤ x ≤ x } ∈ IR
(1)
describes a variable with unknown, but bounded uncertainties and is defined by its infimum x and its supremum x. Classical point-real values are a special case of intervals with x = x = x. Alternatively, an interval can be represented by midpoint and radius
1774
[x] = [xc − x∆ , xc + x∆ ]
(2a)
with x+x x−x mid [x] = xc = , rad [x] = x∆ = . (2b) 2 2 Interval extensions of vectors and matrices (denoted in boldface, e.g. [A], [b]) are obtained by replacing the pointreal entries by interval entries. Calculations with intervals are done independent of the distribution within the bounds using interval arithmetics, cf. e.g. [12]. Toolboxes like I NTLAB [13] allow validated interval calculations in M ATLAB. For the presented methods, the rank of an interval matrix [A] plays an important role. It is defined as the minimal rank of all matrices enclosed by the interval matrix: o n rank [A] := min r ∀A ∈ [A], r = rank(A) . (3) r
Row and column interval rank are defined analogously. As shown in [14], an interval matrix can be tested for full column rank by evaluating the inequality −1 T Ac A∆ < 1 , (4) ρ ATc Ac where ρ(·) denotes the spectral radius of a matrix.
The sensor distortions (i.e. measurement errors) ξ(k) are unknown, but bounded. No distribution, but only the lower and upper bounds ξ(k) and ξ(k) have to be known for each sample k. As a result, the true value of each sample of the output yi is guaranteed to be within the interval of measured data •
ytrue,i (k) ∈ [ymeas,i (k)−ξ∆,i (k) , ymeas,i (k)+ξ∆,i (k)], (6) and analogously for each input ui . If inputs are exactly known, this results in a radius of u∆ = 0. Furthermore, two requirements concerning the interval bounds of the measurements and the sampling time have to be fulfilled, which are introduced in Section III-B. D. Interval Parameter Identification If the order nARX is known, the remaining task of system identification reduces to identify the parameter matrices Ai and B i of the ARX model. As discussed in the introduction, various methods to compute an enclosure of these parameters can be found in literature, e.g. [4], [5], [7], [8]. One popular approach is to formulate the parameter estimation task in the form of an overdetermined interval linear system (ILS) [M ][P ] = [N ]
B. ARX Models The interval identification methods discussed in this paper are based on measurement data of inputs u ∈ Rnu and outputs y ∈ Rny of a linear time-invariant (LTI) system. The data is gathered at sampling points k Ts . The resulting discrete-time system model with sampling time Ts ∈ R+ is described in the AutoRegressive with eXogenous inputs (ARX) form and time (or sampling point) index k ∈ Z: y(k) = −
n ARX X i=1
Ai y(k − i) +
n ARX X i=1
B i u(k − i) + B 0 u(k). (5)
nARX is called the order of the ARX model in the following, assuming that AnARX 6= 0 holds, where 0 is the zero matrix of the appropriate dimensions. The direct feedthrough, if available, is represented by B 0 u(k). This definition easily extends to the interval case by replacing the real inputs, outputs, and parameters by intervals of appropriate dimensions. C. Prerequisites for Interval System Identification Persistent excitation of the inputs is one of the main common prerequisites for classical identification methods of non-autonomous systems [15]. From the collected data, only the controllable and observable part of the system can be identified. Furthermore, the sampling time Ts has to be chosen with respect to the system dynamics, and a sufficiently high number of sampling points is required to obtain accurate results. These basic requirements also hold for interval system identification. In addition, the following assumptions have to be fulfilled: • The real system generates the output sequences ytrue,i (k) according to the system dynamics and the sequences of input values utrue,i (k). Both sequences are possibly unknown due to measurement errors.
(7)
with the data matrix [M ] ∈ IRr×m , the output matrix [N ] ∈ IRr×n and the unknown parameter matrix [P ] ∈ IRm×n for r > m for the MIMO case, which reduces to the data vector [N ] ∈ IRr and the unknown parameter vector [P ] ∈ IRm The solution set of this ILS is given by n o [P ] := P M P = N for some M ∈ [M ], N ∈ [N ] . (8) All measured interval samples are compatible with this resulting set of parameters. In general, the exact solution set has a complicated shape, and its calculation is very costly. Therefore, a box enclosing this exact solution usually is calculated. [4] gives an overview on methods to calculate such solution set enclosures. In this paper, parameter identification is done using the following method, which is explained in more detail in [11]. The overdetermined parameter identification problem is divided into several square ILSs (SILS), and each square problem is solved individually. [12] gives an overview of methods to solve SILS, and most methods for overdetermined ILS can be applied to the square case as well. The intersection of all individual solutions yields an enclosure of the wanted solution set. Although this procedure requires a comparatively high computational effort, and methods yielding tighter enclosers than interval boxes exist, this method is chosen for the sake of simplicity in this paper, because it can be applied directly to the data matrices used for order determination. However, it could be easily exchanged by any other parameter identification method. III. BASIC I DEA OF I NTERVAL S YSTEM I DENTIFICATION This section recapitulates the basic idea of the underlying interval system identification method for MISO systems, which was originally published in [10], [11].
1775
uT (1) uT (2) .. .
y(1) y(2) .. .
uT (2) uT (3) .. .
uT (N − p) y(N − p) uT (N − p + 1)
|
y(2) y(3) .. .
· · · uT (p) · · · uT (p + 1) .. .
y(N − p + 1) · · · uT (N − 1) {z Z
A. Undisturbed Point-Real Data of a MISO system Assume a MISO LTI system without any measurement errors. The measured data (i.e. the true data in this case) from inputs and the single output is combined in data vectors T z(k) = u1 (k) . . . unu (k) y(k) (9) for each time step k. For a single sequence of N time steps, a block Hankel matrix can be built from the vectors z(k): T z (1) z T (2) ... z T (p) z T (2) z T (3) . . . z T (p+1) e = Z . (10) .. .. .. . . . zT (N −p) zT (N −p+1)
...
zT (N −1)
e are time series of inputs and the The block columns of Z output, shifted by one sample for each block column. The parameter p denotes the number of block columns, which has to be chosen greater than the unknown system order p > n. A sufficiently large data set (cf. Section II-C) is assumed, i.e. N − p > p. For systems with direct feedthrough, an additional block column only consisting of input signals has e resulting in an extended matrix to be added to Z, T u (p + 1) .. e (11) Z= . . Z T u (N ) In the rest of this paper, only the general case with direct feedthrough is regarded, which can be easily reduced to the case without direct feedthrough. Based on the matrix Z, the parameters Γ of an ARX model of order p with direct feedthrough can be estimated using Eq. (12), which is shown on the top of this page. Because an ARX model of order p > n is not of interest, a method to determine the unknown order n is required in a first step. In [10], [11], it is shown for MISO systems that this can be done by analyzing the column rank of Z: colrank(Z) = n + p · nu + nu .
(13)
Comparing this to the number of columns of Z cZ = p + p · nu + nu ,
(14)
it can be seen that only the first summand differs. If the parameter p is chosen greater than the system order n, the matrix Z is rank deficient in the absence of any errors. Thus, the system rank can be determined by building Z for p > n and then calculate its column rank.
bp T ap y(p) u (p + 1) y(p + 1) . T y(p + 1) u (p + 2) .. y(p + 2) = .. .. . .. . . b1 y(N − 1) uT (N ) y(N ) a1 | {z } } b Y 0 | {z } (12) Γ
B. Interval Data of a MISO System In the interval case, the matrix [Z ] is built from the interval time series [u] and [y ] analogously to Eq. (9), (10). The matrix Z of the (unknown) true values is always enclosed by [Z ] (cf. Section II-C), thus the interval column rank (cf. Section II-A) of [Z ] can not be greater than the column rank of Z. To avoid an unwanted reduction of the interval column rank compared to the case without uncertainty, two additional prerequisites have to be fulfilled by the interval data: • Firstly, if the interval bounds are chosen so large that some columns have an overlapping interval in each row, the interval column rank is reduced. In this case, a pointreal matrix with additional linear dependant columns is included in the interval matrix. • Secondly, the interval column rank can be reduced by choosing a too small sampling time in comparison to the system dynamics. In that case, the intervals of two consecutive sampling points of one output signal could overlap for all sample points due to the slow dynamics, thus an overlap would exist in the shifted columns of the same signal. To meet these requirements, which is assumed in the following, sensors with sufficiently precise measurements and a reasonable sampling time have to be chosen. Then, the column rank of [Z ] and Z are identical and the system order can be detected by rank inspection. For the interval case, the test for full column rank (4) is conducted for matrices [Z ] built with increasing p, until rank deficiency is detected, i.e. a computational expensive calculation of the interval rank can be avoided. Thus, the order is nARX = p − 1. IV. I NTERVAL MIMO S YSTEM I DENTIFICATION It is not possible to extend the MISO method directly to MIMO systems, as the different outputs can be coupled. These couplings allow the computation of some output values from a combination of other delayed output values and the inputs, which leads to a linear dependency. Therefore, when building the matrix [Z ] as in (11) from the vector [z(k)] with all outputs [y1 (k)] . . . [yny (k)], the matrix can be rank deficient. This prevents the determination of a valid MIMO ARX model order by this naive approach. Additionally, MIMO ARX models are not unique in general. Therefore, a more sophisticated way to obtain the minimal order of a MIMO ARX model is proposed and discussed in the following. Additionally, it is shown how multiple data sequences can be used for interval system identification.
1776
A. Basic Order Determination of MIMO Systems
global: pmax ; eq ]; global: ∀q = 0, . . . , pmax : [U e q,i ]; global: ∀i = 1, . . . , ny , ∀q = 1, . . . , pmax : [y
The basic idea of the proposed algorithm is to determine combinations of shifted input and output sequences which lead to a final data matrix [Z ] with full column rank. Due to the assumption that the inputs are persistently excited (cf. Section II-C), only the combinations of shifted output signals have to be investigated. The combinations are investigated up to a user-defined maximum delay pmax , which should be greater than nARX . To reduce the number of investigations, the algorithm works recursively by excluding impossible combinations. It requires shifted signal sequences of all inputs and each of the outputs1 . For the inputs, sequence matrices are built for shifts q = 0, . . . , pmax : [uT (pmax − q + 1)] [uT (p − q + 2)] max e . [Uq ] = (15a) .. . [uT (N − q)] For each of the outputs (i = 1, . . . , ny ), shifted vectors with q = 1, . . . , pmax are set up: [yi (pmax − q + 1)] [y (p − q + 2)] i max . [y e q,i ] = (15b) .. . [yi (N − q)] e p ] represents the Within the algorithm, the matrix [Z temporary data matrix under investigation up to a maximum delay of p, and Ωp ∈ {0, 1}p×ny represents the combination e p ]. Each row q of the delayed output signals included in [Z of this matrix corresponds to the according time delay of the outputs and each column i to the respective output. Ωq,i = 1 means that yi delayed for q steps is selected in the current e p ]. e q,i ] is included in [Z combination, i.e. [y The algorithm is initially started with p = 1 and Ω0 = [ ]. Each cycle is based on the previously selected signals Ωp−1 . Firstly, the matrix Ωp−1 is extended by each of the 2ny −1 possible combinations ω j ∈ {0, 1}ny , ω j 6= 0,
j = 1, . . . , 2ny −1
(16)
e p,i ], i = 1, . . . , ny , of the output columns with delay p [y yielding the matrices T Ωp,j = ΩTp−1 ω j . (17) e 0 ], . . . , [U e p ] and the e p,j ] is built from [U Each possible [Z e p,j ] selected output columns in Ωp,j . Then, each matrix [Z is tested for full column rank using criterion (4). For each j e p,j ], the next cycle of the with full column rank of [Z algorithm is started with Ωp,j as Ωp−1 and p = p + 1. e p,j ] has full column rank, all If none of the matrices [Z e p,i ] can be expressed as combinations of output columns [y a linear combination of the output columns already included in the temporary data matrix Ωp−1 and the inputs, i.e. no 1 Note that with analogous point-real definitions of Eq. (15), the matrix Z e0 . ep y e1 y in Eq. (12) can be written as Z = U e p,1 · · · U e 1,1 U
function [Ωpmin , pmin ] = ord det(p, Ωp−1 ) { if p == pmax return Ωpmin =NaN, pmin =NaN; %invalid branch for j = 1, . . . , 2ny −1 { Ωp,j = Ωp−1 ; ω Tj ; %add j-th combination of outputs e p,j ] = [U e 0 ], . . . , [U ep ], all [y [Z e q,i ] selected in Ωp,j ; e p,j ] == full %use criterion (4) if colrank [Z [Ωpmin ,j , pmin,j ] = ord det(p+1,Ωp,j ); }; e p,j ] 6= full if ∀j = 1, . . . , 2ny −1 : colrank [Z return Ωpmin = Ωp−1 , pmin = p−1; j = arg(minj (pmin,j )); if isunique(j) %test if minimal solution is unique return Ωpmin = Ωpmin ,j , pmin = pmin,j ; c =user criterion Ωpmin ,j ; %select one of the minimal solutions return Ωpmin = Ωpmin ,c , pmin = pmin,c ; }; Alg. 1.
Pseudo-code of the basic recursive function
additional information on the system is added by the output columns of this delay p. Any further delays of any output(s) in this branch can also be iteratively substituted by the respective linear functions. Therefore, the algorithm terminates the current cycle and returns p − 1 as the wanted order for this branch together with the corresponding matrix Ωp−1 to the superordinate cycle. The recursion also terminates in one branch if no minimal order p < pmax could be determined and this branch is marked as invalid. If different combinations of output columns yield a solution in one cycle, the solution with the lowest order is chosen for this branch and returned to the superordinate cycle. In general, multiple solutions with this minimum order occur, because MIMO ARX models are not unique. This allows a further optimization of the selected model structure. Here, a criterion to maximize the number of zeros in Ωp is used, because each zero entry leads to ny point-real zeros instead of interval parameters to be estimated in the resulting ARX model, thus reducing the uncertainty of the model. If still multiple solutions exist, these are equivalent for our application, and any of them can be chosen. After all recursive cycles are terminated, the overall minimal order pmin is found together with the corresponding matrix Ωpmin . If all branches are marked as invalid, the algorithm returns an error. This can happen if pmax was chosen too small or if assumptions are violated. The recursive function of the algorithm is summarized in Alg. 1. Thereby, the variables marked as global are read-only values available to each cycle of the recursion. B. Optimized Order Determination of MIMO Systems The number of recursions of the algorithm can be further reduced, which is especially helpful for larger numbers of outputs. To do so, the possible combinations ω j of the output e p,i ] have to be processed in a specific columns with delay p [y order within the algorithm. First, the combinations with only
1777
e p,i ] one selected signal are tested. If one single column [y leads to rank deficiency, the true data in this column can be expressed as a linear combination of the true data contained in the outputs selected in Ωp−1 and the true input data by some linear function. Due to this relationship, the rank of matrices in this branch, i.e. based on the same selection of signals Ωp−1 , containing this column will always be deficient. Thus, all other combinations including this column can be removed from the set of possible combinations. Combinations based on the same Ωp−1 without this column of delay p, but with higher delays of the same output occur in other subbranches and do not have to be addressed in this step. If a combination leads to a data matrix of full rank, the next recursion cycle is processed as in the basic algorithm. After all combinations with only one column have been tested, the same approach can be followed for all remaining combinations of two columns. Again, all further combinations including these two columns can be removed if a combination leads to a rank deficient data matrix. This is repeated for all remaining combinations with an increasing number of columns, until all combinations either have been tested or removed from the set of possible combinations for the current delay p. Additionally, the recursion depth (i.e. the maximum delay) e p,j ] can be reduced during runtime. If none of the matrices [Z has full column rank, and p − 1 is returned as minimum order for that branch, the maximum delay can be reduced to pmax = p. Other solutions with the same minimal order p − 1 would be ignored if pmax is reduced to p − 1, not allowing for further optimization of the model structure. The result of this optimized algorithm equals the result of the basic algorithm from the previous subsection, but the number of recursions will be reduced for many problems. C. Discussion of the Order Determination Algorithms Both recursive algorithms always identify the minimal order nARX of the ARX model if pmax > nARX and all other prerequisites are met (cf. Section II-C and III-B). In the basic algorithm, many combinations of possible model configurations are recursively investigated, resulting in a high number of column rank checks. However, because the criterion (4) is used instead of an interval rank calculation, this check is computationally cheap, and the number of checks can be limited by selecting an appropriate termination limit pmax . The maximal number of rank checks is limited to (2ny − 1)pmax , however the actual number is far below this theoretical limit for most applications. The optimized algorithm even further reduces this number, e.g. (2ny −1 − 1) combinations can be removed if one single column combination leads to a rank deficient data matrix. D. MIMO Interval Parameter Identification After the recursive algorithm is successfully completed, the minimal order of the ARX representation of the system as well as the corresponding required signals are known. With this information, the parameters of the model can be identified from the interval data with any method from literature
(cf. Section II-D). The choice of the parameter identification method has significant influence on the overestimation of the parameters, but is not the focus of this paper. Here, the interval parameter identification method from [11], which was summarized in Section II-D, is used for sake of simplicity. For this purpose, a new data matrix [Z ] is constructed based on sequence matrices of all inputs and sequence vectors of the selected outputs in Ωpmin , which are built as in Eq. (15), now using the determined order pmin instead of pmax to make use of all data points. The corresponding unknown parameters to these input and output signals, which are subsumed in the matrix [Γ], are estimated using Eq. (12), where [Y ] is now a matrix with one column for each output signal. All other parameters of the ARX model are set to point-real zero, because the corresponding signals are not required within the model. E. Multiple Data Sequences If multiple data sets should be used for identification instead of one single sequence, these data sets cannot be simply concatenated due to the dynamics of the system. Instead, the vectors (15) have to be built for each sequence individually, as well as the matrix [Y ] for parameter identification. Afterwards, these shifted vectors, and thus the matrices [Z ] for each sequence, can be stacked. This assures that for the solution of the interval linear systems only data points from the same sequence are handled in one equation, i.e. one row of the matrices [Z ] and [Y ]. This extension can be applied analogously to the existing method from [10], [11]. V. E XAMPLE To demonstrate the proposed interval identification method for MIMO systems, it was implemented using the Intlab package [13] in M ATLAB. In the following, a second order system with two outputs, three inputs and direct feedthrough is regarded. The parameters of this ARX system (5) are − 0.5 0 0.2 0 , (18a) , A2 = A1 = 1 0 0 0.3 − 1 −1 0 0 0 1 B1 = , B2 = , (18b) −2 1 0 1 0 0 0 1 0 B0 = . (18c) 0 0 0 To generate a set of identification data of the true (undisturbed) system, this system is simulated to obtain 4 sequences of 2 500 sampling points each. Then, uniformly distributed random errors with a maximum amplitude of 1‰ of the respective maximum signal amplitude is superimposed on the simulated signals (inputs and outputs) to reproduce measurement errors. The bounds of the interval data sets are set to the same values, i.e. all prerequisites are fulfilled. The proposed interval system identification method results in an ARX model of the correct minimal order nARX = 2 with the interval parameter matrices [ 0.17 ; 0.23 ] [ −0.02 ; 0.02 ] [A 1 ] = , (19a) [ −0.04 ; 0.03 ] [ 0.28 ; 0.32 ]
1778
1
This is an enclosure of the true parameters (18). The pointreal zeros in the second column of the matrix [A2 ] result from the finally used combination of output signals 1 1 Ω2 = . (20) 1 0 This indicates that the second output delayed for two steps has no unique influence on both outputs and thus is not used for identification. Nevertheless, it could be added without causing rank deficiency, because the alternative solution with the same ARX order 1 1 Ω2,rejected = (21) 1 1 was rejected by the selection criterion to maximize the number of point-real zeros (cf. Section IV). Initialized with pmax = 10, i.e. (2ny − 1)pmax = 59 049, both algorithms had a maximum recursion depth of four, and 90 column rank checks (4) were performed by the basic algorithm, 63 by the optimized algorithm. To analyze the behavior of the set of systems included in the interval model, a Monte-Carlo analysis to calculate the eigenvalues was carried out with 100 000 point-real models randomly drawn from the interval model. Fig. 1 shows the true eigenvalues marked by blue asterisks. The convex hulls of the eigenvalues from the randomly drawn systems are shown as filled red areas. One eigenvalue of the interval model is point-real and located in the origin due to the pointreal zeros in [A2 ]. Thus, all randomly drawn models also have an eigenvalue exactly in the origin. Another eigenvalue is always real, but varies on the real axis. Together with the zoomed plot, it can be seen that the true system dynamics is included in the identified interval model, and the latter describes a set of systems with similar dynamics. Thus, the identified model is a valid description for the true system. VI. C ONCLUSIONS In this paper, a method for the identification of MIMO LTI systems from interval measurement data was presented. The main contribution are the determination algorithms for the minimal order of an ARX model of the system. It was also shown how multiple data sequences can be used for this class of interval system identification methods. The recursive order determination algorithms also select a minimum number of required delayed output signals, resulting in a maximum number of point-real zero entries
−0.67
0.5 Im
[ −0.51 ; −0.49 ] 0 [A 2 ] = , (19b) [ 0.99 ; 1.01 ] 0 [ −0.03 ; 0.03 ] [ 0.97 ; 1.03 ] [ −0.03 ; 0.02 ] [B 0 ] = , [ −0.03 ; 0.03 ] [ −0.04 ; 0.04 ] [ −0.03 ; 0.03 ] (19c) [ −1.03 ; −0.97 ] [ −1.04 ; −0.96 ] [ −0.02 ; 0.02 ] [B 1 ] = , [ −2.03 ; −1.97 ] [ 0.96 ; 1.04 ] [ −0.03 ; 0.03 ] (19d) [ −0.07 ; 0.07 ] [ −0.04 ; 0.04 ] [ 0.98 ; 1.02 ] [B 2 ] = . [ 0.92 ; 1.08 ] [ −0.05 ; 0.05 ] [ −0.03 ; 0.03 ] (19e)
−0.73 0.07
0 −0.5 −1
0.13
true system rand. drawn systems −1
Fig. 1.
−0.5
0 Re
0.5
1
Eigenvalues of true and identified systems
in the parameter matrices of the model, thus minimizing the uncertainty in the model. Besides a basic algorithm, an optimized version was presented which further decreases the number of recursions. With the results of these algorithms, an interval parameter identification method from literature is used to determine the non-zero model parameters. The overall method was demonstrated on a simulation example, showing the low number of recursions in practical use. R EFERENCES [1] R. E. Moore, Interval analysis. Englewood Cliffs: Prentice-Hall, 1966. [2] P. Planchon, Guaranteed Diagnosis of Uncertain Linear Systems Using State Set Observation. Berlin: Logos-Verlag, 2007. [3] F. Wolff, Konsistenzbasierte Fehlerdiagnose nichtlinearer Systeme mittels Zustandsmengenbeobachtung, ser. Schriften des Instituts f¨ur Regelungs- und Steuerungssysteme, Karlsruher Institut f¨ur Technologie. Karlsruhe: KIT Scientific Publishing, 2010, vol. 09. [4] J. Hor´acˇ ek and M. Hlad´ık, “Computing Enclosures of Overdetermined Interval Linear Systems,” Reliable Computing, vol. 19, no. 1, pp. 142– 155, 2013. [5] M. Milanese and A. Vicino, “Optimal estimation theory for dynamic systems with set membership uncertainty: An overview,” Automatica, vol. 27, no. 6, pp. 997–1009, 1991. [6] V. Cerone, D. Piga, and D. Regruto, “Enforcing stability constraints in set-membership identification of linear dynamic systems,” Automatica, vol. 47, no. 11, pp. 2488–2494, 2011. [7] ——, “Set-Membership Error-in-Variables Identification Through Convex Relaxation Techniques,” IEEE Transactions on Automatic Control, vol. 57, no. 2, pp. 517–522, 2012. [8] M. Pouliquen, E. Pigeon, and O. Gehan, “Output error identification for multi-input multi-output systems with bounded disturbances,” in 50th IEEE Conference on Decision and Control and European Control Conference (CDC-ECC), 2011, pp. 7200–7205. [9] L. H. Keel, J. S. Lew, and S. Bhattacharyya, “System identification using interval dynamic models,” in American Control Conference, vol. 2, 1994, pp. 1537–1542. [10] S. Zaiser, M. Buchholz, and K. Dietmayer, “Black-Box-Modellierung mit unsicheren Parametern aus Messdaten mit unbekannten, aber beschr¨ankten Fehlern,” in Tagungsband GMA-Fachausschuss 1.30 ”Modellbildung, Identifikation und Simulation in der Automatisierungstechnik”, O. Sawodny and J. Adamy, Eds., VDI/VDE-GMA. Darmstadt: Technische Universit¨at Darmstadt, Institut f¨ur Automatisierungstechnik und Mechatronik, 2013. [11] ——, “Black-box modeling with uncertain parameters from measurement data with unknown, but bounded errors,” at - Automatisierungstechnik, vol. 62, no. 9, pp. 607–618, 2014. [12] R. E. Moore, R. B. Kearfott, and M. J. Cloud, Introduction to Interval Analysis. Society for Industrial and Applied Mathematics, 2009. [13] S. Rump, “INTLAB - INTerval LABoratory,” in Developments in Reliable Computing, T. Csendes, Ed. Dordrecht: Kluwer Academic Publishers, 1999, pp. 77–104. [14] J. Rohn, “A Handbook of Results on Interval Linear Problems.” ˇ Prague, Czech Republic, Tech. Institute of Computer Science AS CR, Rep. V-1163, 2012. [Online]. Available: http://hdl.handle.net/11104/ 0212095 [15] L. Ljung, System Identification: Theory for the User, 2nd ed. Upper Saddle River: Prentice Hall, 1999.
1779
