Identification of a Flexible Beam Setup SC4040 - Filtering and Identification Roxana Bucsa - 4330110 Radu Florea - 4330358 Remy Kabel - 4132165
n o i t a c fi i t n e d I d n a 0 g n 4 i 0 r e 4 t l i C S F
Delft Center for Systems and Control
Table of Contents
1 Introdu Introducti ction on
1
1-1 Descri Descripti ption on of of the system system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1-2 Descriptio Description n of of the the assig assignment nment . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
2 Data Data PrePre-pr proces ocessin sing g
3
2-1 Data Data analys analysis is
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2-2 Spectru Spectrum m analys analysis is . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5 6
2-3 Detrending, Detrending, Data Data Compres Compression sion and and Order Estima Estimation tion . . . . . . . . . . . . . . . 3 Model Estima Estimation, tion, Valida Validation tion and and Optimizat Optimization ion
7
3-1 Hanke Hankell Matri Matrixx Test Test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3-2 Estimatio Estimation n of States-Spa States-Space ce Matrices Matrices . . . . . . . . . . . . . . . . . . . . . . . . 3-3 3-4 3-5 3-6 3-7 3-8
Model Model Valid Validati ation on . . . . Auto-Co Auto-Correlat rrelation ion Test . . Cross-Co Cross-Correlat rrelation ion Test . . Cross-Va Cross-Validat lidation ion Test . . Variance Variance Accounte Accounted d For For . Model Model Optimi Optimizat zation ion . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
7 8 8 8 9 11 11 12
4 Conc Conclu lusi sion on
15
A Identifica Identification tion Matla Matlab b Script Script
17
i
ii
Table of Contents
ii
Chapt hapteer 1
Introduction
As a finalizing project for the TU Delft course SC4040 Filtering and System Identification, the task has been assigned to us to identify, validate, and optimize several data sets for a certain experimental model. This puts the knowledge gathered throughout the course to the test in a practical setting.
1-1 1-1
Desc Descri ript ptio ion n of of the the sy syst stem em
Figure 1-1: The experimental setup of the flexible beam
The experimental setup considers a flexible beam structure as depicted in Figure 1-1 1-1.. Six piezos are bonded to the beam, three of which are visible in the photograph and three of which which are located at the rear of the beam. These piezos piezos can either either be used to develop develop a strain strain when a voltage is applied to them (causing the beam to bend) or they can detect a strain by 1
2
Introduction
emitti emitting ng a volta voltage ge when the beam is bent. In the setup setup used here, here, piezos piezos 1 and 3 are are used used as actuators, whereas piezos 2, 4, 5 and 6 are used as sensors. For this system the maximum voltage voltage that can be supplied supplied is 250 Volt Volt before the actuators actuators saturate. saturate. The constraint constraint on the sensor outputs is 10 Volt before the electronic amplifiers saturate.
1-2 1-2
Desc Descri ript ptio ion n of the the ass assig ignm nmen entt
As a starting point, seven different data sets are given with inputs and corresponding outputs for all piezos. The channels (rows) in the data sets are defined as: •
1 - time (seconds)
•
2-5 - measurements from piezos 2, 4, 5 and 6 (Volts)
•
6-7 - signals sent to piezos 1 and 3 (Volts)
•
8 - sample counter (not used)
The task at hand is to identif identify y the real life system. system. First, First, howev however, er, all data data sets sets must must be pre-processed to find any unsuitable sets. In the following chapters we will follow the method of identification as described in figure 4.1 in [ 1, p.65].
2
Chapt hapteer 2
Data Pre-processing
As noted previously, before starting the identification of the data gathered in the experiment, the dataset datasetss must must first be analyz analyzed ed to discar discard d any any unusa unusable ble sets. sets. As such, such, this this chapt chapter er describes how the available data sets are processed in order to decide upon the discarding of the sets.
2-1 2-1
Data Data anal analys ysis is
Whereas one might like to jump in to start doing any such processing operations, it is a good idea to first do a general analysis of the collected data and try to observe any particularities within it. To this end, the obtained obtained input-output input-output data has been plotted; to better study it. The obtained plots are presented in Figure 2-1 and Figure 2-2 Figure 2-2.. It can immediately be observed that, for dataset 3 (Figure 2-2), 2-2 ), the output electronic amplifiers have been saturated. This renders the whole data set unusable due to the fact that it is not possib p ossible le to extract extract any any knowledg knowledgee on the system’s system’s beha b ehavior vior beyond the saturation saturation point. This is information which might be crucial in obtaining an accurate estimate of the system. As such, dataset 3 will be dismissed for the rest of the identification process. Another way of analyzing the datasets is by looking at the provided .mat data file. From this we can obtain the actual number of samples that the different data sets contain, as well as the sampling period that was used to obtain the data. Using the notation: ts,i for the time sample of data set i, and N s,i s,i for the number of samples of data set i , with i = 1 ,..., 7, the following values are observed:
ts,1 = 0.001 001ss N s, s,1 = 50000
ts,2 = 0.001 001ss N s, s,2 = 50000
ts,5 = 0.001 001ss N s, s,5 = 750
ts,3 = 0.001 001ss N s, s,3 = 20000
ts,6 = 0.00005 00005ss N s, s,6 = 10000
3
ts,4 = 0.001 001ss N s, s,4 = 20000
ts,7 = 0.001 001ss N s, s,7 = 10000
4
Data Pre-processing
100
100
0
0
−100 0 200
5
10
15
20
25
30
35
−100 0 200
40
0
10
15
20
25
30
35
40
5
10
15
20
25
30
35
40
0
−200 0 500
] V [ e d u t i l p m A
5
5
10
15
20
25
30
35
−200 0 500
40
0
0
−500 0 10
5
10
−500 0 10
15
0
5
10
15
5
10
15
0
−10 0 100
5
10
−10 0 100
15
0
0
−100 0 100
0.1
0.2
0.3
0.4
0.5
−100 0 100
0.4
−100 0 100
8
−100 0
0
0.1
0.2
0.3
0.4
0.5
0
−100 0 100
0.05
0.1
0.15
0.2
0.25
0. 3
0.35
0
0. 05
0. 1
0.15
0.2
0.25
0.3
0.35
0.4
1
2
3
4
5
6
7
8
0
−100 0
1
2
3
4
5
6
7
Time [s]
Figure 2-1: Input signals of the data sets.
Piezo 1 (left) and piezo 3 (right). From top to bottom are datasets 1 to 7.
10
10
10
10
0
0
0
0
−10 0 10
20
40
0
−10 0 10
20
40
0
−10 0 10
20
40
0
] −10 V [ 0 e 1 d u 0 t i l p −1 0 m 5 A
20
−10 0 10
20
40
20
−10 0 1
−1
10
20
0
10
20
−10 0 10
0
0
0
−5
−5
0
0. 2
0. 4
20
40
10
20
0
0. 2
0. 4
10
−10 0 10
−10 0 1
0
10
20
−1 0 5
0
0.2
0.4
10
0
0
0
0
−10 0 5
−10 0 5
0. 4
−10 0 5
10
−5 0
0.4
0.2
0.4
0
0
0
−5
−5
−5
0
5
10
0
5
10
20
40
10
20
10
20
0 −5 0 10
−10 0 5
0.2
40
0
5
−5
20
0
−10 0 1
−1
−10 0 10 0
0
5
10
40
0
0
10
20
0
0
10
−10 0 10
0. 2
0.2
0.4
0. 2
0. 4
5
10
0
0
5
Time [s]
Figure 2-2: Output signals of the data sets.
From left to right: piezo 2, piezo 4, piezo 5, and piezo 6 From top to bottom are datasets 1 to 7.
4
2-2 Spectrum analysis
5
Through looking at the values for the running time of the experimental analyses and the numbe numberr of sample sampless used used for this timesp timespan, an, we can exclud excludee two two datase datasets: ts: datase datasett 5 and 6. Similar findings can be gathered from the data presented in Figure 2-2 2-2.. Despit Despitee datas dataset et 6 having a large number of samples, the sample frequency used is very large resulting in the shortest shortest running time of the experiment, experiment, namely 0.5 seconds. Similarly Similarly,, dataset dataset 5 also only encompasses a running time of 0.75 seconds. Basically, such running times of gathering experimental data are too short to analyze the slower dynamics of the system. As such, part of the information of the experiment is lost upon using either of these datasets.
2-2 2-2
Spect Spectru rum m anal analys ysis is
In order to reduce the influence of noise in certain areas of the spectrum, we pre-filter the remaining remaining four signals. signals. Before Before pre-filtering pre-filtering the data, it should be noted that the useful signal and the disturbance can be separated because the signal has signal power mainly in the low low freque frequenci ncies, es, while while the noise noise contai contains ns mostly mostly high high freque frequenci ncies. es. This This ma may y be observ observed ed by plotting the spectrum of the signal with the help of the MATLAB function pwelch. This returns the power spectral density estimate of the discrete-time signal vector using the Welch’s averaged, modified periodogram method. The results have been plotted in Figure 2-3. 2-3 .
30
30
30
30
20
20
20
20
10
10
10
0
0
0
−10
−10
−10
−20
−20
−20
−30
−30
−30
−40
−40
−40
−50
−50
−50
−60
−60
−60
) 10 e l p m 0 a s / d a r −10 / B d ( y−20 c n e u−30 q e r f / r −40 e w o P−50
−60
−70 0
0.5
1
−70 0
0.5
1
−70 0
Normalized Frequency (
0.5
1
−70 0
0.5
1
rad/sample)
× π
Figure 2-3: Power spectrum of the output signals. From left to right: dataset 1, dataset 2, dataset 4, and dataset 7.
Based on the spectrum in Figure 2-3 and on output signals from Figure 2-2, a rather low power and low output amplitude can be observed for data set 4. This gives an indication that this data set can be highly sensible to noise, thus it is not as reliant as the rest, therefore it will be eliminated from the data sets considered for the identification process. 5
6
Data Pre-processing
4
4
10
3
3
10
2
10
2
10
2
10
1
0
10
3
10
10
4
10
10
1
10
20
S−value
30
10
0
1
10
20
30
S−value
10
0
5
10
15
20
S−value
Figure 2-4: Singular values
From left to right: dataset 1, dataset 2, dataset 4, and dataset 7)
2-3
Detren Detrendin ding, g, Data Data Compres Compressio sion n and Order Order Estimat Estimation ion
Before starting the subspace identification of the model, first, some more preliminary actions have to be taken. To start, the signals will be detrended. This pre-filtering is done using the MATLAB function detrend and removes unknown offsets and/or linear trends from signals. Furthermore, in order to later be able to do the cross-validation procedure (presented in Chapter 3 Chapter 3), ), the measured data is split into two parts: an identification part (representing the 2 /3 of the data set) and a validation part (which is the remaining 1/ 1 /3 of the data). first 2/ As a method of subspace identification, we are looking to implement the PO-MOESP using the supplied MATLAB MATLAB toolbox. toolbox. For this, several several steps have to be taken. taken. Firstly Firstly,, the order of the system system should be estimated. estimated. Consecutiv Consecutively ely,, the A and C matrices matrices may be estimated estimated and finally, the B and D matrices of the system subspace may be estimated. In this chapter we will restrict restrict ourselves ourselves to the necessary necessary preliminar preliminary y actions. actions. The actual estimation estimation of the state-space of the system will be discussed in Chapter 3. 3 . In pursuance of the order of the system, an estimate may be drawn from the singular values. To obtain the singular values, we use the available "ord" function - dordpo in the case of PO-MOE PO-MOESP SP - in order to compre compress ss the data data and to genera generate te a model model order order estimate estimate.. For this, we only need to specify the block-size parameter s, which should be considered larger than the expected order of the system, s=30. The obtain obtained ed vector vector S contains the singular values values and will be presented shortly. Additionally, Additionally, a compressed matrix R is obtained which will further be used to estimate the system matrices A and C in the next chapter. The singular values are plotted by the use of the MATLAB function semilogy. Based on the results plotted in Figure 2-4, 2-4 , we can choose to estimate a model of order 9-16.
6
Chapt hapteer 3
Model Estimation, Validation and Optimization
In the previous previous chapter, chapter, the datasets datasets have have been analyzed and pre-processed. pre-processed. Remaining Remaining are the three datasets (data set 1, data set 2 and data set 7) which are usable and, supposedly, grant enough information to construct an accurate model of the entire system. As such, this chapter chapter will entail the estimation estimation and validati validation on of the model. Upon the completion completion of this, the gathered results may be optimized.
3-1 3-1
Hank Hankel el Matr Matrix ix Test est
In order to verify the richness property of the input, which is mathematically called persistency of excitation , the Hankel matrices specific to the inputs of every data set need to be verified. This actually means that for the available available data sets,the order of every Hankel matrix is checked checked to have full rank s, where s = 3n and n is the desired system order. The Hankel matrix has the form [2] form [2]::
U 0,s,N =
u (0) u (1) .. .
u (1) · · · u (N − 1) u (2) · · · u (N ) N ) .. .. .. . . . u (s − 1) u (s) · · · u (N + s + s − 2)
(3-1)
It is computed using the MATLAB function hankel. Running Running the Hankel Hankel matrix test for a desired model of order 9 ( n = 9, s = 27) and order 16 ( n = 16, s = 48), we obtain the ranks for the Hankel matrices, presented in Table 3-1. From the gathered results it is clear that the inputs used in data set 2 are not persistently [27,, 48] (i.e. n [9, [9, 16]). Therefore exciting of order s , when s [27 Therefore,, it is safe to conclude that data set 2 is of no use in estimating a system of the desired order. 7
8
Model Estimation, Validation and Optimization Table 3-1: Hankel matrix rank for different desired model orders
Data set 1 2 7
3-2
Hankel Hankel matrix matrix rank for input 1 n = 9 ( s = 27) n = 16 (s = 48) 27 48 14 14 27 48
Hankel Hankel matrix matrix rank for input 2 n = 9 ( s = 27) n = 16 ( s = 48) 27 48 14 14 27 48
Estima Estimatio tion n of States States-Sp -Space ace Matric Matrices es
As described before in paragraph 2-3 2-3,, after estimating the order of the system, the a model of the system may be created by constructing the state-space matrices. First, the system matrices A, C as well as the Kalman gain K are estimated using the toolbox function dmodpo. These estimates estimates are of the system in innov innovation form. form. In order to establish establish this, the inforinformation previously obtained, namely the data matrix R, is used. Consecutiv Consecutively ely,, the toolbox function dac2bd is used to estimate B and D .
3-3 3-3
Model Model Valid alidat atio ion n
( Ae,Be,Ce,De)) is available, the next step is to validate Considering that an estimated model (Ae,Be,Ce,De this this model model on the basis basis of a more more detail detailed ed inspectio inspection. n. In order order to do this, we use the tests tests available within the class of residual tests, namely the auto-correlation and cross-correlation tests. tests. These validat validation ion methods are applied to innov innovation models in which which an optimal model means that the variance variance of the one-step ahead ahead predictor predictor error is minimal. For a parameter parameter estimate ˆ of , the one-step-ahead prediction error ˆ k, θ) can be estimated via the Kalmanfilter equations [2, p.38 [2, p.387]. 7]. If the system system under consideratio consideration n can be represen represented ted by a model in the innov innovation form, the prediction prediction-err -error or sequence sequence for equal to the global optimum, calculated using data from an open-loop experiment, should be the following:
•
The sequence ˆ k, θ) is a zero-mean white-noise sequence.
•
u (k ). The sequence ˆ k, θ) is independent of the input sequence u(
In order to obtain the vector ˆ k, θ) we use the toolbox function dltisim. For this, we have to provide the innovation state-space model previously obtained and we obtain the residual vector, as the difference between the actually measured output y(k) and the innovation model output yˆ k, θ).
3-4
Auto-Co uto-Corre rrelat lation ion Test
After obtaining the residual vector ˆ k, θ), we have to check the whiteness of the residual vector. vector. We can do this by computing its auto-corr auto-correlati elation on function. function. To do this, we use the MATLAB Signal Proceesing Toolbox function xcorr. The output of this function function will give give us the cross-correlation between two signals, and the cross correlation of ˆ k, θ) and ˆ k, θ) equals 8
3-5 Cross-Correlation Test
9
4
2
x 10
0 −2 −3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0 Lag
1000
2000
3000
4
1
x 10
0 −1 −3000 4 x 10 2 0 −2 −3000 4 x 10 1 0 −1 −3000
Figure 3-1: The auto-correlation function for dataset 1.
From top to bottom: outputs 1 to 4.
an auto-correlation. The auto-correlation is defined as a signal of length 2N+1, corresponding to the auto-correlation from lag -MAXLAG up to and including lag +MAXLAG . The results have been plotted in Figure 3-1 for dataset 1 and in Figure 3-2 3-2,, for dataset 7. The decision to first estimate an order 9 model has been taken in order to be able to evaluate the auto-correlation and cross-correlation functions specific to the estimated models from dataset 1 and dataset 7. This will help in deciding if both the datasets are suited for system identification.
3-5
CrossCross-Co Corre rrelat lation ion Test
As previously mentioned, the cross-correlation test is a residual test that can be used on innov innovation state-spa state-space ce models. The role of this test is to check check the second second property of the residual vector, as specified in the previous section: the residual vector has to be independent of the input signal u(k). For this, we are going going to use the same residua residuall vecto vectorr comput computed ed before. The only difference is the fact that we now now compute compute the cross-co cross-correl rrelation ation between between the residual and the input signal using the same MATLAB function xcorr. We use the same value for the lag as before. The results results have been plotted in Figure 3-3 3-3 for for dataset 1 and in Figure ?? for datas dataset et 7. Based Based on these, these, it can clear clearly ly be seen seen that the residua residuall sequence sequence ˆ (k, θ), of the model estimated using dataset 7, is not a zero-mean white-noise sequence and the cross-correlation between it and the input is higher than the one of the residual vector obtained obtained from the estimated estimated model using dataset dataset 1. Thus, Thus, dataset 7 does not give out a good model estimate, therefore it will be discarded. 9
10
Model Estimation, Validation and Optimization
50 0 −50 −3000 100
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0 Lag
1000
2000
3000
0 −100 −3000 100 0 −100 −3000 100 0 −100 −3000
Figure 3-2: The auto-correlation function for dataset 7.
From top to bottom: outputs 1 to 4.
5
2
5
x 10
2
5
x 10
2
5
x 10
2
1
1
1
1
0
0
0
0
−1
−1
−1
−1
−2 −5000
0
5000
−2 −5000
5
2
0
5000
−2 −5000
5
x 10
2
0
5000
−2 −5000
5
x 10
2
2
1
1
1
0
0
0
0
−1
−1
−1
−1
0
5000
−2 −5000
0
5000
−2 −5000
Lag
5000
0
5000
0
5000
x 10
−2 −5000
Figure 3-3: The cross-correlation function for dataset 1.
From top to bottom: outputs 1 to 4.
10
0 5
x 10
1
−2 −5000
x 10
3-6 Cross-Validation Test
11
1000
1000
1000
1000
500
500
500
500
0
0
0
0
−500
−500
−500
−500
−1000 −5000
0
5000
−1000 −5000
0
5000
−1000 −5000
0
5000
−1000 −5000
1000
1000
1000
1000
500
500
500
500
0
0
0
0
−500
−500
−500
−500
−1000 −5000
0
5000
−1000 −5000
0
5000
−1000 −5000
Lag
0
5000
−1000 −5000
0
5000
0
5000
Figure 3-4: The cross-correlation function for dataset 7.
From top to bottom: outputs 1 to 4.
3-6
CrossCross-Va Valid lidati ation on Test
It may happen that the we cannot know for sure that the model adequetly describes the system system dynamics. dynamics. This may happen due to overfiittin overfiitting, g, which actually means that the model complexity or the number of model parameters has become so large with respect to the length of the data sequence that the predicted output matches very accurately the identification data. In order to prevent this from happening, the measured data is split in two parts as touched upon in section 2-3 section 2-3:: an identification part and a validation part. This means that the identification part (we consider the first 2/3 of all samples), after which we can apply the validation tests on the second part (the remaining 1/3 of all samples) of the measured data. This This procedu procedure re has been applie applied d to all datasets datasets before startin startingg any any other other procedu procedure. re. The presented auto and cross-correlation plots are obtained using the validation part of the given data.
3-7
Varian ariance ce Accounte ccounted d For
The variance accounted for (VAF) is defined as a scaled variant of the cost function J N N θ ) and it is defined as (formula page 4.4 page 78 companion). This can be easily computed in MATLAB by using the toolbox function vaf for all the components components of the output output signal. It requires the measured measured output y(k) (only the validation part) and the predicted output yˆ k ). The value valuess obtained obtained are between between 0 % and 100%; the higher the value of the VAF, the lower the prediction error and the better the model obtained through identification.
In order to decide the order of the most accurate model, the system has been estimated for ...16), after which the VAF, for each one of these models, has eight different orders ( n = 9...16 11
12
Model Estimation, Validation and Optimization Table 3-2: VAF for different estimated model orders
Model Model order order 9 10 11 12 13 14 15 16
VAF of Output Output 1 90.4186 90.3580 90 90.3500 90 90.2027 90 90.0076 90 89.8875 89 89.7378 89 89.5581 89
VAF of of Output Output 2 91.4660 91.5305 91.5012 91.6933 91.7611 91.7555 91.7518 91.5508
VAF of Output 3 93.0888 93.2522 93.2521 93.5027 93.5970 93.6187 93.5852 93.6436
VAF of Output 4 93.7086 93.8219 93.8259 94.0585 94.1724 94.1894 94.1758 94.2237
been computed. The results are presented in Table 3-2. 3-2 . As can be seen, the estimated model of order 9 gives a very good accuracy, its output VAFs being comparable comparable to the the ones of the higher order order estimates. estimates. This is a very very good reason reason to use the order 9 model estimate, as it reduces system complexity and, at the same time, computation time for the rest of the operations, such as estimate optimization, future controller computations (a system’s model is estimated so that a controller can be designed afterwards) etc.
3-8 3-8
Model Model Opti Optimi miza zati tion on
( Ai , Bi , C i , Di , K i ) matrices, the model is through an optimization Using the initial estimated (A algor algorithm ithm so that that better better results results can be obtain obtained. ed. This This is achiev achieved ed throug through h the use of the provided MATLAB toolbox command doptlti. The obtained model matrices are give by:
Ao =
0.7724 −0.7920 0.0337 0.0929 −0.0605 −0.1919 −0.0078 0.1477 −0.0508
0.7935 0.6872 0.0122 0.1910 −0.0489 0.2100 −0.0195 0.0034 0.0477
0.1698 −0.1471 0.9382 0.3808 0.1664 −0.1458 0.3554 −0.4658 −0.0049
Bo =
0.107 0.0899 0.0448 0.8952 −0.0404 −0.2375 −0.0923 0.0546 −0.0717
0.0007 0.0011 0.0026 −0.0062 0.0028 −0.0030 0.0174 0.0082 −0.0013
12
0.1296 0.0351 −0.1057 0.1502 0.9939 0.0054 −0.0404 −0.1428 0.0481 0.0000 0.0011 −0.0008 −0.0010 0.0028 0.0013 0.0138 0.0122 −0.0013
0.2650 −0.1153 −0.0803 0.2970 0.0251 0.8045 0.0707 −0.0705 0.0242 −
0.1336 0.0925 0.0207 0.0050 0.0151 0.2785 0.1141 0.9261 −0.0725 −
0.0374 0.0610 0.0516 0.0101 −0.0543 −0.0420 −0.7133 0.3279 0.0129
0.3688 0.3206 0.1421 0.1580 −0.0050 −0.2024 −0.0453 −0.0087 0.9184
3-8 Model Optimization
13
2000 1000 0 −1000 −3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
−2000
−1000
0
1000
2000
3000
4000 2000 0 −2000 −3000 2000 1000 0 −1000 −3000 2000 1000 0 −1000 −3000
Lag
Figure 3-5: The auto-correlation function the optimized system.
From top to bottom: outputs 1 to 4.
C o =
0.2020 0.1542 −0.4502 −0.4579
0.2522 −0.2048 0.2364 0.2260 −
0.2812 −0.2983 −0.0584 −0.0498
Do =
K o =
0.2944 −0.3286 0.8820 −0.2048 0.3168 0.2226 −0.2730 −0.4411 −1.0180
0.2393 0.0734 −0.3457 −0.3501 −
0.1023 −0.2876 −0.0774 −0.0713
0.0088 0.0019 0.0015 −0.0005 0.2694 0.3281 −0.2919 −0.3820 0.0674 −0.2087 0.5600 0.5226 −0.2229 −
0.0070 0.0038 0.0003 0.0011 0.3719 0.2746 0.5999 0.0147 0.1973 −0.4069 −0.2664 0.1155 0.7588 −
0.3575 0.1158 −0.1791 −0.1689
0.0723 −0.0334 0.1279 0.1247
0.1073 0.7015 0.5601 0.5389
1.1787 0.2465 0.0420 0.0594
−
0.0525 0.0259 −1.3204 −0.2319 −0.6010 0.3511 0.2181 1.3118 −0.6026 −
98..6287 6287,, 98 98..1546 1546,, 99 99..2903 2903,, 99 99..2927[], which The optimized model estimate gives better VAFs: 98 is a significat significat increase increase in model accuracy accuracy. This can be also observed observed in the auto-correla auto-correlation tion and cross-correlation plots presented in Figure 3-5 and Figure 3-6 Figure 3-6..
13
14
Model Estimation, Validation and Optimization
4
4
x 10
6
6
4
x 10
5
5
x 10
1
x 10
4
4
0.5 2
2
0
0 −2
0
0
−2 −0.5
−4
−4
−6 −5000
0
5000
−6 −5000
4
0
5000
−5 −5000
4
x 10
4
4
5000
5000
0
5000
x 10
4
4
2
2
0
0
−2
−2
1 0
0
4
x 10
2
0
−1 −5000
4
x 10
3
2
0
−1 −2
−2 −3
−4 −5000
0
5000
−5000
−4 0
5000
−5000
−4 0
5000
−5000
Lag
Figure 3-6: The cross-correlation function the optimized system.
From top to bottom: outputs 1 to 4.
14
Chapt hapteer 4
Conclusion
The task at hand was to identify the real life system as presented in the introduction. Through data pre-processing, model estimation and model validation, an accurate data set has been identified, modeled and optimized. Firstly, the data sets have been pre-processed after which four out of seven datasets have been discarded for various reasons. Through a simple analysis of the datasets in terms of running time and sampling frequency, datasets 5 and 6 were discarded due to short running times of only 0.75s and 0.5s respectively. Such short running times of gathering experimental data are too short to analyze the slower dynamics of the system. Consecutively, the remaining four datasets have been analyzed using spectrum analysis. The power spectra found showed that a rather low power and low output amplitude can be observe served d for data set 4. This This gives gives an indica indicatio tion n that that this this data data set can be highly highly sensibl sensiblee to noise, thus it is not as reliant as the rest. With only three datasets remaining the final pre-processing started by detrending and compressing the data after which the order is estimated through singular value analysis, showing an order of around around 9-16. Consideri Considering ng this fact, a Hankel matrix matrix test shows shows that the inputs used in data set 2 are non persistently exciting of sufficient order and as such it is discarded. Continu Continuing, ing, the model’s state-spac state-spacee matrices matrices have been b een estimated and validated. alidated. Furtherurthermore, auto- and cross-correlation tests could be processed for datasets 1 and 7. This showed that for dataset 7, the residual sequence of the estimated model is not a zero-mean white-noise sequence sequence and the cross-corr cross-correlati elation on is higher higher than the one obtained obtained from dataset 1. As such the 7th dataset is also discarded and we are left with solely dataset 1. Finally, the model is optimized to improve accuracy which is analyzed using the VAF-values. For the first output an increase of about 8% is found and for the second output and increase of about 6% is found.
15
16
Conclusion
16
Append ppendix ix A
Identification Matlab Script
As explained in Chapter 2 Chapter 2,, the datasets have been preprocessed using the following MATLAB script. Upon its completion, the results have been used as argument to reject several sets. 1 2 3
cl c ; clear all ; close all ;
4 5
load ( ’datasets.mat’ )
6 7 8 9 10 11 12 13
% % % % % % %
data3 - cli data3 clippe pped d (out (output puts s sat satura urated) ted) data6 dat a6 - sh short ort exp experi erimen ment t d at at a5 a5 - n ot ot e no no ug ug h i nf nf or or ma ma ti ti on on data2 dat a2 - no not t per persis sisten tently tly exc exciti iting ng d at at a4 a4 - l ow ow p ow ow er er ( s pe pe ct ct ru ru m ) , l ow ow a mp mp li li tu tu de de ( m or or e s en en si si bl bl e t o n oi oi se se distur dis turban bances), ces), not as den dense se d at at a7 a7 - c or or el el la la ti ti on on i s n ot ot g oo oo d
14 15 16 17
u = [ ] ; y = [ ] ; t = [ ] ;
18 19 20 21
u_val = [ ] ; y_val = [ ] ; t_val = [ ] ;
22 23 24 25
% Spl Split it da data ta fo for r cr cros oss-valid s-validati ation on fo r i = 1 : 7 data = eval ( sprintf ( ’data%d’ , i ) ) ;
26 27
N = length ( data ) ;
28 29 30 31
% Ge Get t ide identi ntific ficati ation on da data ta u { i } = data ( 6 : 7 , 1 : floor ( 2 ∗ N / 3 ) ) ’ ; y { i } = data ( 2 : 5 , 1 : floor ( 2 ∗ N / 3 ) ) ’ ;
17
18
Identification Matlab Script t { i } = data ( 1 , 1 : floor ( 2 ∗ N / 3 ) ) ’ ;
32 33 34 35 36 37 38
en d
% Ge Get t eva evalua luatio tion n da data ta u_val { i } = data ( 6 : 7 , floor ( 2 ∗ N /3)+1: N ) ’ ; y_val { i } = data ( 2 : 5 , floor ( 2 ∗ N /3)+1: N ) ’ ; t_val { i } = data ( 1 , floor ( 2 ∗ N /3)+1: N ) ’ ;
39 40 41 42 43 44
% % P l ot ot d at at a % Input figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ; fo r i = 1 : 7 N = length ( data ) ;
45 46 47 48
subplot ( 7 , 2 , 2 ∗ i −1) plot ( t { i } , u { i } ( : , 1 ) ) ; grid on ;
49 50 51 52 53
subplot ( 7 , 2 , 2 ∗ i ) plot ( t { i } , u { i } ( : , 2 ) ) ; grid on ; en d
54 55 56 57 58 59 60
% Out Output put figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ; fo r i = 1 : 7 subplot ( 7 , 4 , 4 ∗ i −3) plot ( t { i } , y { i } ( : , 1 ) ) ; grid on ;
61 62 63 64
subplot ( 7 , 4 , 4 ∗ i −2) plot ( t { i } , y { i } ( : , 2 ) ) ; grid on ;
subplot ( 7 , 4 , 4 ∗ i −1) plot ( t { i } , y { i } ( : , 3 ) ) ; grid on ;
65 66 67 68 69 70 71 72 73
subplot ( 7 , 4 , 4 ∗ i ) plot ( t { i } , y { i } ( : , 4 ) ) ; grid on ; en d
74 75 76 77
%% Eli Elimin minate ate da data ta3, 3, da data5 ta5 an and d da data6 ta6 u{3} = [ ] ; u{5} = [ ] ; u{6} = [ ] ; y{3} = [ ] ; y{5} = [ ] ; y{6} = [ ] ;
78 79 80
u_val { 3 } = [ ] ; u_val { 5} 5} = [ ] ; u_val { 6 } = [ ] ; y_val { 3 } = [ ] ; y_val { 5} 5} = [ ] ; y_val { 6 } = [ ] ;
81 82 83
u = u ( ~ cellfun ( ’isempty’ , u ) ) ; y = y ( ~ cellfun ( ’isempty’ , y ) ) ;
84
18
19
85 86
u_val = u_val ( ~ cellfun ( ’isempty’ , u_val ) ) ; y_val = y_val ( ~ cellfun ( ’isempty’ , y_val ) ) ;
87 88 89 90 91 92 93
%% Obt Obtain ain sig signal nal spe spectr ctrum um figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ; fo r i = 1 : 4 subplot ( 1 , 4 , i ) ; pwelch ( y { i } ) ; en d
94 95 96 97
%% Eli Elimin minate ate da data4 ta4 ud { 3 } = [ ] ; yd { 3 } = [ ] ;
98 99 100
ud_val { 3} 3} = [ ] ; yd_val { 3} 3} = [ ] ;
101 102 103
ud = ud (~ cellfun ( ’isempty’ , ud ) ) ; yd = yd (~ cellfun ( ’isempty’ , yd ) ) ;
104 105 106
ud_val = ud_val ( ~ cellfun ( ’isempty’ , ud_val ) ) ; yd_val = yd_val ( ~ cellfun ( ’isempty’ , yd_val ) ) ;
107 108 109 110
%% Det Detren rend d da data ta ud = [ ] ; yd = [ ] ;
111 112 113
ud_val = [ ] ; yd_val = [ ] ;
114 115 116 117
fo r i = 1 : 3 ud { i } = detrend ( u { i } ) ; yd { i } = detrend ( y { i } ) ;
118 119 120 121
ud_val { i } = detrend ( u_val { i } ) ; yd_val { i } = detrend ( y_val { i } ) ; en d
122 123 124 125
% % C he he ck ck s in in gu gu la la r v al al ue ue s t o o bt bt ai ai n a n e st st im im at at e o f t he he s ys ys te te m o rd rd er er s = 30; R = [ ] ;
126 127 128 129
figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ; fo r i = 1 : 3 [ S , R { i } ] = dordpo ( ud { i } , yd { i } , s ) ;
130 131 132 133 134
subplot ( 1 , 4 , i ) semilogy ( S , ’* ’ ) ; grid on ; en d
135 136 137
% % O bt bt ai ai n H an an ke ke l m at at ri ri ce ce s o f t he he i np np ut ut s ( c he he ck ck i f U i s r an an k 9 - > i n pu pu t i s % per persis sisten tently tly exci excitin ting) g)
19
20
138 139 140
Identification Matlab Script
n = 9; s = 3∗ n ; r = [ ] ;
141 142 143 144 145 146 147
fo r i = 1 : 3 N = length ( ud { i } ) ; U1 = hankel ( ud { i } ( 1 : s , 1 ) , ud { i } ( s : N , 1 ) ) ; U2 = hankel ( ud { i } ( 1 : s , 2 ) , ud { i } ( s : N , 2 ) ) ; r ( i , 1 ) = rank ( U1 ) ; r ( i , 2 ) = rank ( U2 ) ;
148
if r ( i , 1 ) ~= ~= s st r = sprintf ( ’ Da Da ta ta s et et % d i s n o t g oo oo d ’ , i ) ; disp ( st r ) ; en d
149 150 151 152 153
if r ( i , 2 ) ~= ~= s st r = sprintf ( ’ Da Da ta ta s et et % d i s n o t g oo oo d ’ , i ) ; disp ( st r ) ; en d
154 155 156 157 158
en d
159 160
r
161 162 163 164
%% Eli Elimin minate ate da data2 ta2 ud { 2 } = [ ] ; yd { 2 } = [ ] ;
165 166 167
ud_val { 2} 2} = [ ] ; yd_val { 2} 2} = [ ] ;
168 169 170
ud = ud (~ cellfun ( ’isempty’ , ud ) ) ; yd = yd (~ cellfun ( ’isempty’ , yd ) ) ;
171 172 173
ud_val = ud_val ( ~ cellfun ( ’isempty’ , ud_val ) ) ; yd_val = yd_val ( ~ cellfun ( ’isempty’ , yd_val ) ) ;
174 175 176
%% PO-MOESP Subs Subspace pace Ident Identificat ification ion yi = [ ] ;
177 178 179 180 181
fo r i = 1 : 2 [ S , R ] = dordpo ( ud { i } , yd { i } , s ) ; [ Ai , Ci , Ki ] = dmodpo ( R , n ) ; [ Bi , Di ] = dac2bd ( Ai , Ci , ud { i } , yd { i } ) ;
182
Ae Be Ce De
183 184 185 186
= = = =
Ai−Ki ∗ Ci ; [ Bi−Ki ∗ D i K i ] ; Ci ; [ Di −ey e ( 4 , 4 ) ] ;
187 188 189 190
en d
% Obt Obtain ain res residu idual al vec vector tor epsilon { i } = dltisim ( Ae , Be , Ce , De ,
20
[ ud { i } yd { i } ] ) ;
21
191 192 193
%% Obt Obtain ain th the e au auto-cor to-correl relati ation on plo plots ts maxlag = 3 0 0 0 ;
194 195 196
fo r i = 1 : 2 figure ( ’Position’ , [ 4 0 0 1 5 0 8 20 20 4 2 5 ] ) ;
197
fo r j = 1 : 4 ac i = xcorr ( epsilon { i } ( : , j ) , epsilon { i } ( : , j ) , maxlag ) ; subplot ( 4 , 1 , j ) plot (− ma xl ag : maxlag , ac i ) ; grid on ; en d
198 199 200 201 202 203 204
en d
205 206 207 208
%% Obt Obtain ain th the e cr cros oss-corre s-correlat lation ion plots plots fo r i = 1 : 2 figure ( ’Position’ , [ 4 0 0 1 5 0 8 20 20 4 2 5 ] ) ;
209
fo r j = 1 : 4 xc i = xcorr ( epsilon { i } ( : , j ) , ud { i } ( : , 1 ) , maxlag ) ; subplot ( 2 , 4 , 2 ∗ j −1) plot (− ma xl ag : maxlag , xc i ) ; grid on ;
210 211 212 213 214 215
xc i = xcorr ( epsilon { i } ( : , j ) , ud { i } ( : , 2 ) , maxlag ) ; subplot ( 2 , 4 , 2 ∗ j ) plot (− ma xl ag : maxlag , xc i ) ; grid on ;
216
217 218 219
en d
220 221
en d
222 223 224 225
%% Eli Elimin minate ate da data7 ta7 ud { 2 } = [ ] ; yd { 2 } = [ ] ;
226 227 228
ud_val { 2} 2} = [ ] ; yd_val { 2} 2} = [ ] ;
229 230 231
ud = ud (~ cellfun ( ’isempty’ , ud ) ) ; yd = yd (~ cellfun ( ’isempty’ , yd ) ) ;
232 233 234
ud_val = ud_val ( ~ cellfun ( ’isempty’ , ud_val ) ) ; yd_val = yd_val ( ~ cellfun ( ’isempty’ , yd_val ) ) ;
235 236 237
%% Obt Obtain ain VAFs n = 9:16;
238 239 240
fo r i = 1 : length ( n ) s = 3∗ n ( i ) ;
241 242 243
[ S , R ] = dordpo ( ud { 1 } , yd { 1 } , s ) ; [ Ai , Ci , Ki ] = dmodpo ( R , n ) ;
21
22
Identification Matlab Script
[ Bi , Di ] = dac2bd ( Ai , Ci , ud { 1 } , yd {1}) ;
244 245
Ae = Ai−Ki ∗ Ci ; Be = [ Bi−Ki ∗ D i K i ] ; Ce = Ci ; De = [ Di zeros (4 ,4 ) ] ; % Sim Simula ulate te th the e ide identi ntifie fied d sys system tem yi { i } = dltisim ( Ae , Be , Ce , De , [ ud_val { 1} yd_val { 1 } ] ) ;
246 247 248 249 250 251 252
% Cal Calcul culate ate th the e VA VAF F disp ( ’ VA VA F f or or o rd rd er er ’ ) ; n ( i ) vaf_po = va f ( yd_val { 1 } , yi { i } ) en d
253 254 255 256 257
% % R un un o pt pt im im iz iz at at io io n f or or m od od el el o f o rd rd er er 9 n = 9; s = 3∗ n ;
258 259 260 261
[ S , R ] = dordpo ( ud { 1 } , yd { 1 } , s ) ; [ Ai , Ci , Ki ] = dmodpo ( R , n ) ; [ Bi , Di ] = dac2bd ( Ai , Ci , ud { 1 } , yd {1}) ;
262 263 264 265
[ Ao , Bo , Co , Do , ~ , Ko ] = doptlti ( ud { 1 } , yd { 1 } , Ai , Bi , Ci , Di , [ ] , Ki ) ;
266 267
Ae = Ai−Ki ∗ Ci ; Be = [ Bi−Ki ∗ D i K i ] ; Ce = Ci ; De = [ Di zeros (4 ,4 ) ] ; % Sim Simula ulate te th the e ide identi ntifie fied d system system yi c = dltisim ( Ae , Be , Ce , De , [ ud_val {1 } yd_val { 1 } ] ) ;
268 269 270 271 272 273 274 275
Ae o = Ao−Ko ∗ Co ; Be o = [ Bo−Ko ∗ D o K o ] ; Ce o = Co ; De o = [ Do zeros (4 ,4 ) ] ; % Sim Simula ulate te th the e opt optimi imized zed sys system tem yo c = dltisim ( Aeo Ae o , Beo Be o , Ceo Ce o , Deo De o ,
276 277 278 279 280 281
[ ud_val {1 } yd_val { 1 } ] ) ;
282 283 284 285 286 287
% Cal Calcul culate ate the VA VAF F disp ( ’VAF non-opti non-optimized’ mized’ ) ; vaf_po = va f ( yd_val { 1 } , yi c ) disp ( ’VAF optim optimized’ ized’ ) ; vaf_po = va f ( yd_val { 1 } , yo c )
288 289 290 291 292 293
%% Ae Be Ce De
Obtain au Obtain auto to an and d cr cros oss-corre s-correlat lation ion fo for r opt optimi imized zed sys system tem Ao−Ko ∗ Co ; [ Bo−Ko ∗ D o K o ] ; Co ; [ Do −ey e ( 4 , 4 ) ] ;
= = = =
294 295 296
% Obt Obtain ain res residu idual al vec vector tor epsilon = dltisim ( Ae , Be , Ce , De ,
[ ud {1 } yd { 1 } ] ) ;
22
23
297 298
figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ;
299 300 301 302 303 304 305
fo r j = 1 : 4 ac i = xcorr ( epsilon ( : , j ) , epsilon ( : , j ) , maxlag ) ; subplot ( 4 , 1 , j ) plot (− ma xl ag : maxlag , ac i ) ; grid on ; en d
306 307 308 309 310 311 312
figure ( ’Position’ , [ 4 0 0 1 5 0 8 2 0 4 2 5 ] ) ; fo r j = 1 : 4 xc i = xcorr ( epsilon ( : , j ) , ud { 1 } ( : , 1 ) , maxlag ) ; subplot ( 2 , 4 , 2 ∗ j −1) plot (− ma xl ag : maxlag , xc i ) ; grid on ;
313 314 315 316 317 318
xc i = xcorr ( epsilon ( : , j ) , ud { 1 } ( : , 2 ) , maxlag ) ; subplot ( 2 , 4 , 2 ∗ j ) plot (− ma xl ag : maxlag , xc i ) ; grid on ; en d
23
24
Identification Matlab Script
24
Bibliography
Filteringg and System System Identi Identific ficatio ation: n: An [1] M. Verhaeg erhaegen, en, V. Verdult erdult,, and N. Bergboer Bergboer,, Filterin Introduction to using Matlab Software . Delft University of Technology, August 2007.
[2] M. Verhaegen Verhaegen and V. Verdult, Verdult, Filtering and System Identification . Cambridge University Press, 2007.
25
