ELEC4632 Lab 1

An introduction to linear least square method and system identification

In this lab, you will use linear least square method to identify a model for an input/output data collected from a physical system.

Pre-lab Exercise: MATLAB Refresher

You should write a MATLAB code (m-file) to create a set of sinusoidal data and plot them exactly as shown inFig. 1, through the following steps.

Fig. 1. Comparison of a set of sinusoidal data, (a) Original data, (b) Cut-off data.

1.   Create a column vector representing time from 0 to 100 seconds with time spacing of 0.1 second and name it t. Use either linspace function or use direct vector definition like t = a:b:c. Make sure to transpose it as this form. of variables are row vectors in MATLAB by default. Type in help + function_name, such as linspace, in MATLAB command window to learn how to use the function, and you can also see full help on any function by searching them in the main help page of MATLAB. Do not know which function can achieve your objective? Google is a good searching tool.

2.   Use sin and cos functions to generate two sinusoids, y1 and y2, respectively, with a period of 100 seconds, and add a uniformly distributed noise using rand function with the bound of ±0.2.

3.   Cut off the first 200 samples from all data and assign them to new variables t_new, y1_new, and y2_new, respectively. Make sure t_new starts from zero. You can extract part of a data stored in a variable by using colon operator in indexing rows and columns, i.e., A(a:b,:) would extract data in matrix A from row a to row b for all columns (recall that A(c,d) in MATLAB returns the element stored in row c and column d of A)

4.   Store these three cut-off vectors under a new matrix variable named data_new. Use data_new = [t_new y1_new y2_new] to attach column vectors with the same size and create an n-by- 3 matrix, i.e., data_new n×3 . This method is known as matrix concatenation.

5.   Plot the original data you created in Step 2 against time, i.e., y1 versus t and y2 vs t, on the top side of the figure with the same colors as shown in Fig. 1(a) using the following functions: figure, subplot, plot, stairs, hold on, hold off, xlim, ylim, title, grid, xlabel, ylabel, legend. To display Greek letters, like τ , read the Text Properties in MATLAB help. Limit x-axis and y-axis between (0, 140) and (一1.5, 1.5).

6.   Use the same functions above (except figure) to plot the cut-off data at the bottom of the figure as shown in Fig. 1(b). The color code for cos(0.02*pi*t) is  [0.6 0.7 0.8] (see the help for plot and RGB color codes)

Introduction to System Identification

In this lab, you will learn how to determine or identify a suitable dynamic model for a process using linear least squares method [1], [2]. In control systems theory, to design a control system for a process using a so-called Model-Based Control method, a mathematical model of the process is needed. In general, there are two methods to model a process (also known as system identification). One is by using the physical relationships that describe the dynamics of the process. For instance, using Newton’s Iaws to obtain equations of motion for a mass-spring-damper process, or using Kirchhoff’s VoItage and Current Iaws in addition to Newton’s Iaws to derive a permanent magnet DC motor  equation [3].  This  method  can  sometimes  result  in  complex  dynamic  equations  for complicated  processes.  Another  way  of  finding  a  dynamic  model  for a  process  is  by  using experimental input and output data, and then trying to match them with a mathematical equation, either linear or nonlinear, depending on the nature of the system and its operating conditions. This method is useful for the cases where there is little information about the physics of the process, or due to the limitation in accessing different parts of the process for parameter measurements, only input and output signals are available to be used. We all know that a process would reveal its dynamic characteristics through its outputs and states variables if a proper input signal is applied to excite all internal modes of the process.

One of the most common methods of empirical system identification in industry using input/output data is the so-called Step Response modelling. This method is suitable for processes with slow dynamics which are mostly controlled for regulation purposes or set-point control. An example is provided in Appendix section showing how to find a first order model using a recorded data in continuous-time domain.

Parametric System Identification Using Linear Least Squares Method

Most of the processes in industry are controlled in a way to keep the systems near their operating points, and as a result, they can be represented by linear differential or difference equations (rational transfer function). Therefore, we can determine the finite number of parameters that characterize such a dynamic model, i.e., polynomial coefficients or zeroes and poles. This method is known as parametric system identification. However, it is not always possible to model these processes with a first-order transfer function using one simple step response as discussed before, since their dynamics can be more complicated. Therefore, by choosing a suitable input signal and recording the corresponding output using a digital computer, we can determine causal discrete-time models, as a dynamic model for a process to be controlled. The reason for identifying the models in discrete- time form. is that the control system will eventually be implemented on a digital computer, so it would be beneficial to have a discrete-time model of the process and design the control system directly in discrete-time domain, even though the real nature of the process is in continuous-time domain. In this lab, our focus is on second-order discrete-time models with strictly proper form. It should be mentioned that, when no information is available about the process, we need to vary the order of the model to get the best approximation of the actual process in terms of closeness of the responses.

Linear discrete-time dynamic models in second-order form. are given as follows,

Transfer function in Z domain (z-1  is called “delay operator”):

(1)

State Space:

(2)

(Canonical Controllable Form)

Difference Equation:

y(k) + a1y(k -1) + a2y(k - 2) = b1u(k -1) + b2u(k - 2).                                                                 (3)

Thus, the unknown parameters of the process model to be determined are {a1 , a2, b1 , b2}. In order to estimate these parameters, the so-called linear least squares method is used in this lab, which is a well-known and commonly used method for this purpose. More details on least squares method can be found in[1]and[2]. From the three different discrete-time models above, difference equations in Eq. (3) is more suitable since it can be written as a set of linear equations containing discrete- time values of the measured input/output data. Later on, for the purpose of controller design, state space model will be used primarily. Thus, the difference equation can be rearranged as follows,

y(k ) = -a1y(k - 1) - a2 y(k - 2) + b1u(k - 1) + b2 u(k - 2),

(4)

The model y(k) = φT(k)θ is called regression model and the variables  inside φT(k) are called regression variables or regressors.  They  contain  the  previous  values  of  the   input  and corresponding output values at time steps k − 1 and k − 2 (k is the discrete time step related to continuous time as t = kTs for k = 1, 2, 3 … with Ts as the sampling time). Since we use computer to apply input signal with a fixed sampling rate and then measure the corresponding output, N linear equations are constructed, as shown in Fig. 2, with only four unknown parameters for a second- order model, where N is the number of samples.

Fig. 2. System identification with input/output data.

Therefore, the problem of solving a set of linear equations is shown as below,

(5)

(6)

(7)

As can be seen in Eq.(7), vector Y and matrix Φ contain all recorded input and output information, and θ is the vector of unknown model parameters to be found. The result is similar to the case where we are dealing with solving a set of linear equations having more unknown variable than the number of equations, which indicates that there is no unique solution to this problem. Least squares method can offer the best solution for model parameters {a1 , a2 , b1 , b2} in the sense of minimum error between the left hand side and the right hand side of those linear equations. Moreover, we know that the real process is not perfectly linear, and the actual measurements y(k) in each time step does not have a linear relationship with the previous measurements and the input values as the model implies, i.e., Y ≈ Φθ. Thus, the error between the actual output values and the model output values, i.e., E = Y 一 Φθ, includes both measurement noise and model uncertainty as well. Finally, the least squares solution that gives the best estimate for model parameters θ(ˆ) is obtained as follows,

(8)

The obtained solution θ(ˆ) in Eq.(8)is the best estimate of the model parameters in the sense

of least squares (note that θ = [一a1 一a2 b1 b2]T whereas we seek [a1 a2 b1 b2]T ). It is clear, however, that the main condition for solving the above matrix equation (also known as normal equations) is that the matrix Φ TΦ has full rank and it is positive-definite (well-conditioned), which is known as Persistent Excitation condition. More details on least square solution in Eq. (8) are provided in Appendix.

Selection of a Suitable Input

So far, we have learned how to find the best estimate of coefficients for a system model using linear least  squares.  In  the  next  step,  we  want  to  know  how  to  perform data collection  for  system identification. The most recommended types of input signals for system identification are Pseudo Random Binary Signal (PRBS), Square-Wave Signals, or the sum of many sinusoidal waves. The amplitude of the input signal should be bounded to keep the system in its linear operation region. It should be mentioned that any process has constraints on its input in terms of the maximum range of the input amplitude allowed to be applied before the process is damaged. In addition, the resulting output y(k) should be much larger than measurement noise. Otherwise, the output would not be usable for system identification and it might create a so-called ill-conditioned ΦTΦ; and therefore, the persistent excitation condition is no longer met.

DC Offset Compensation

The  output  of  a  linear  time-invariant  (LTI)  system  is  always  zero  when  zero  input  is applied (assuming zero initial conditions). Thus, if a system responds to zero input with the nonzero value, the system is called to have offset value in its output (DC offset), which probably exists in industrial processes. The other form. of offset in systems is when the input range is limited to only positive or negative values. In this case, the value in the middle of the input range is called input offset, and the corresponding output to this input is called output offset. In both above-mentioned cases, the offsets should be removed from both input and output data before filling the matrix Φ to be able to find an LTI  model  accepting  both  negative and  positive  input  values.  However,  in the actual control operation, the control input u(k) should be calculated from offset-free output (output offset yoffset should be subtracted from the measured output ym (k) to be used in the computation of u(k)). Then, input offset uoffset should be separately added to the calculated control input to find the actual process input uin(k) as shown in Fig. 3.

Fig. 3. Compensation of the input and output offsets.

Model Verification

After performing data collection and then identifying the unknown parameters of a process model using  linear  least  squares,  it  is  time  to  verify  the  obtained  dynamic  model  to  see  whether  it demonstrates similar behavior. to the actual process or not. Hence, the same input signal should be applied to both the model and the process, and then the output responses should be compared as illustrated inFig. 4. It is preferred to use a different set of input/output data for validation than those used for system identification. The smaller the error between the output responses, the better the quality of the identified model will be. One good criterion for verification of the identified model is the application of Mean Squared Error method (MSE) on the difference between output responses or the error signal.

Fig. 4. Model verification after performing system identification.

Lab Exercise (2 marks)

In this lab, you are going to identify a second-order discrete-time linear model as in Eq. (1)-(3) using a pre-collected data from one of the W-T setups through the following steps,

1. Data extraction and analysis (0.6 marks, checked at 45 minutes)

a.   Download the pre-collected data from Moodle and save it in the current directory of MATLAB. The data file is named SysIdenData_StudentVersion.mat. Load the data into MATLAB Workspace using load function. The loaded data should appear in Workspace with the name LogData in Structure format. Use the following code to extract individual data in vector array form,

t = LogData.time;

y_act = LogData.signals(1).values(:,2);

y_actm = LogData.signals(1).values(:,1);

u_act = LogData.signals(2).values;

As you can see, the data contain recorded time in t, actual noise-reduced output in y_act, original measured output in y_actm and actual input data in u_act. Then, find the sampling time for which the data were recorded and name it Ts or h.

b.   Plot these data using the functions you practiced in Pre-lab exercise as shown in Fig. 5.

Fig. 5. Original data, (a) Comparison of noise-reduced and measured output signals, (b) Actual input signal.

c.   Find the best estimate of the offset value from the noise-reduced output y_act and name it y_offset. Then, remove the output offset by subtracting y_offset from y_act and name it  y. Find input offset as well and name it u_offset and remove it from input data u_act and name it  u. Then, plot them in one figure as shown in Fig. 6. As you can see, both offset- free input and output signals in Fig. 6(a) and (b), respectively, begin from zero.

Fig. 6. Offset-free data, (a) Offset-free output signal (noise-reduced), (b) Offset-free input signal.

Hint: As explained in DC Offset Compensationsection, for a set of input-output data which only contain positive values, like u_act and y_act here, offset values have to be detected and removed from the data before filling the matrix Φ . If no information is given about the range of input signal, we would assume the first value of the input signal is input offset as you can clearly see inFig. 5(b). However, we cannot simply consider the same fact for output offset as there is always some noise in the output signal, even in the noise-reduced one y_act here. The best approach to estimate the closest value to actual output offset is by taking average from the first period of output data which corresponds to the first period of input data before the first change. Thus, you have to write a loop in MATLAB to automatically detect the number of samples in the first period of output to be used in average. Use while function for loop and mean function for average.

2. Identifying a second-order discrete-time linear model (0.6 marks, checked at 1 hour 30 minutes)

a.  Create matrix Φ as shown in Eq.(6)using the first half of the offset-free input-output data. If starting from k = 1, you need to choose the values for y(−1), y(0), u(−1), and u(0) as initial conditions. They should be chosen rationally. You do not necessarily need to start from k = 1. You can choose to start from some samples ahead, i.e., starting from k = 10, for example. You can use matrix concatenation to create matrix Φ as explained in pre-lab exercise.

b.   Find the solution of least squares method θ(^) as given in Eq.(8), to obtain the estimate of the second-order discrete-time model parameters {a1 , a2, b1 , b2}.

c.   Create a second-order discrete-time transfer function and state space representation of the identified model in MATLAB, as given in Eq. (1) and (2), and display them on Command Window using tf and ss functions.

3. Model verification and simulation (0.8 marks, checked at 2 hours 30 minutes)

a.   Using the transfer function or state space model you defined in previous part, simulate the identified model, which means finding the response of the identified model to the offset-free input u as explained inModel Verificationsection.

b.  You should simulate the model and plot the results, first by using the second half of the offset-free input u and compare the simulated output with the second half of the offset-free output y, and then, by using the entire offset-free input u and compare the simulated output with the entire offset-free output y as shown inFig. 7(a) and (b), respectively.  Use lsim or filter functions to find the simulated response. You could try Simulink to simulate the model too to practice on Simulink too.

c.   In Fig. 7(a), can you explain why the simulated output does not start from the same point as the offset-free output y?

Fig. 7. Model verification, (a) comparison of 2nd  half of the simulated output with the 2nd  half of actual offset-free output, (b) comparison of the entire simulated output with actual offset-free output.

Post-lab Exercise

In this exercise, we want to compare the effect of different model orders in the accuracy of system identification.

1.   Follow the same procedure you did in Exercise 2 to identify a first-order model using offset-free data u and y. The structure of a first-order difference equation is given as below,

y(k) = −a1y(k −1) + b1u(k −1).                                                                                             (9)

2.   Simulate the first-order model using the entire offset-free input and compare its simulated output with the second order one you obtained before and the offset-free output y as shown in Fig. 8.

Fig. 8. Comparison of different order models with the actual offset-free output for accuracy purposes.

3.   Use mean squared error method to numerically compare between these two models as can be seen in the top-left corner of Fig. 8. Which model would you choose for controller design and why?