System Identification assignment 2018


Due before 12:00, 9th March 2018. Please submit a pdf version of your report via blackboard.

Here you will also find the files

 >> loadass.m (will load the data into a suitable vector)
 >> crib.m (available as a skeleton for your program)
 >> dord2.m (a way of generating a test system)

loadass2.m crib.m dord2.m

A Matlab `diary' from the lesson on Wed 31-01-2018 is available here

The assignment has three parts, all relating to recursive least squares system identification.

Part A develops software to do recursive least squares identification.
Part B provides some real data and allows you to identify a mechanical "master" robot of a master/slave tele-manipulator.
Part C (worth 20% of the marks) extends the RLS algorithm to either include instrument variables or a PLS noise model based on the Moving Average Model, or both.

It is recommended that you use Matlab for the assignment, although you will not be penalised for using other programming languages.

Part A Create a recursive least squares identification program.

You may want to use the skeleton code crib.m to set up variables for your recursive least squares estimator (see crib.m above)

  1. Write a recursive least squares program (ARX) to identify a model given input and output data u and y, and the number of numerator and denominator coefficients.The model will also need initial values for P(0) and theta(0).
  2. Choose a stable second order sampled model. Models generated by drmodel do not compare well with the data for part B so I suggest you either either copy across the function dord2 and use >> [B,A]=dord2, or use the following model. B=[0 0.0416 0.0395] A=[1 -1.5363 0.8607]
  3. Generate any input (eg step u=ones(1,50), sampled sin wave u=5*sin(1:.2:70), random u=.5*randn(1,200), etc.) and get the model response to this. You can use the filter program >> y=filter(B,A,u).
  4. Test your recursive least squares program to see if it is predicting the correct coefficients.
  5. Test the response if Gaussian noise is added to the output. eg >> y=filter(B,A,u)+ 0.1*randn(size(u))

You should initialise your model with theta(0)=[0..0] and P(0)=LN*eye(p), where LN is a large number. You should plot how the coefficients in theta track towards the "true" values [A B] as each new iteration occurs.

Tabulate the true coefficients, the level of noise and the final value of the estimate. Comments on different types of inputs, and different lengths of data are encouraged.

Extra marks if you include a forgetting factor, and observe its effects.

Note that if you use the matlab filter function to compare estimateded and actual response, you will need to add a leading 1 to the a parameters and a leading 0 to the b parameters

Part B real data from a mechanical "master" robot

In this part you will build a model of the system parameters for the master of a head controlled telerobot. This was the master component of a head operated master/slave telerobot. To improve the performance of a master/slave telerobot it was necessary to determine an accurate model for the master dynamics. To do this the step response was studied with different conditions of Kp and Kd in the PID controller.


>> loadass2 

to load two data sets as specified by your name. Or alternatively download the data set from

Either way you should have a variable $y$ in the Matlab workspace containing the corresponding response data to a step input.

Files containing the response to a 60mm step input.

You will have to generate the corresponding step input, the command

>> u=60*ones(size(y)); 

will suffice.

The results

You should estimate parameters for the following four models

theta=[a1 a2 b1 b2]
theta=[a1 a2 a3 b1 b2]
theta=[a1 a2 a3 b1 b2 b3]
theta=[a1 a2 a3 a4 b1 b2 b3]

For each model plot the development of the coefficients with each sample.

Determine the best model by simulating the response to a step input of 60mm and plot the predicted and actual response on the same graph. (use either the filter or the dstep command ). This plot should be discussed in your report.

For all the models do a z transform and find the location of the final poles and plot these on the unit circle. Is the predicted model stable? Which model does the best job of predicting the mechanical system?

hint: poles depend only on the a coefficients. a0 is taken as 1. Thus you can form a vector

>>  a1=theta(1);
>>  a2=theta(2);
>>  a3=theta(3);
>>  den=[1 a1 a2 a3];

Then plot these on the complex plane with the following commands

>>  pzmap(1,den);
>>  hold on;
>>  axis([-1.5 1.5 -1.5 1.5]); 
>>  i=1:101;plot(cos(i*pi*2/100),sin(i*pi*2/100));

Part C extends the RLS algorithm (worth 20% of the marks)

to either include instrument variables or a PLS noise model based on the Moving Average Model, or both.

Selecting data files for Part B

Data files for Part B are available at :

Selecting a data file for Part B Use the table below and the associated algorithm to select two data files.

Using the first and last letters of your initials get two numbers by assigning A=1, B=2 C=3 etc. Multiply these two numbers together and take the remainder after division by 31 Locate the number in the table below eg my initials are W.H. -> W=23, H=8 -> 23*8=184 -> 184div 31 = 5 remainder 29 -> files are xp50d15.dat and xp30d03.dat

NumberFiles NumberFiles
0,30xp10d00.dat, xp10d03.dat 1xp10d00.dat, xp10d10.dat
2xp10d01.dat, xp10d15.dat 3xp10d01.dat, xp20d07.dat
4xp20d00.dat, xp20d15.dat 5xp20d00.dat, xp30d07.dat
6xp20d01.dat, xp30d10.dat 7xp20d01.dat, xp30d15.dat
8xp30d00.dat, xp40d10.dat 9xp30d00.dat, xp40d15.dat
10xp30d01.dat, xp50d10.dat 11xp30d01.dat, xp50d15.dat
12xp30d03.dat, xp50d20.dat 13xp30d03.dat, xp10d03.dat
14xp40d00.dat, xp10d10.dat 15xp40d00.dat, xp10d15.dat
16xp40d01.dat, xp20d07.dat 17xp40d01.dat, xp20d15.dat
18xp40d03.dat, xp30d07.dat 19xp40d03.dat, xp30d10.dat
20xp50d00.dat, xp30d15.dat 21xp50d00.dat, xp40d10.dat
22xp50d01.dat, xp40d15.dat 23xp50d01.dat, xp50d10.dat
24xp50d03.dat, xp50d15.dat 25xp50d03.dat, xp50d20.dat
26xp40d03.dat, xp30d07.dat 27xp30d00.dat, xp40d15.dat
28xp20d01.dat, xp30d15.dat 29xp50d15.dat, xp30d03.dat

File naming convention

The file name indicates the test axis, and the proportion and derivative values of the PID controller (The integral term was set to zero). For example file xp10d07.dat contains ascii data from an experimental run on the x axis where the proportional gain was set to 10 and the derivative gain was set to 7.


Parts A and B of your report should include graphs showing how the parameter values inevolve as the iteration becomes more accurate as well as graphs that comparing the actual output and the estimated output. Part A should compare two different inputs eg a step and random input, discuss how noise effects the estimates. As well as giving the comparison between the different models, Part B should show at least on plot of the parameters settling over time, (plot(Theta')). You should also show a plot of the unit circle that give the positions of poles and zeros of the estimated system. Part C should give some evidence that the algorithm is working, for example a test system and the estimated parameters, see examples given in class based on Soderstrom and Stoica. In any report, make sure you give a meaningful caption to any graphs as well as labeling axes. Your report should include listings of all your programs. Don't go overboard with pages of plots showing essentially similar data (put similar data onto one plot)