Model Building and Parameter Estimation

See the attached file.

This problem refers to the model building and parameter estimation in the presence of noise.

For full details, see the attachment.
The task is to implement a simulation model that generates data (simulating real-world measurements), and then fit two types of models to these data. One of the models is a simple polynomial model, and the other is the correct model, but with an unknown parameter.

The goal of this exercise is to get more familiar with some of the tools for model-building.
The exercise is to be performed in MATLAB as per instructions.

The system consists of an input signal x(t), known to be generated as an exponential function of time, t, but where the parameter a is unknown. This signal is then corrupted with experimental noise and then sampled using a sampling time Ts. The resulting sequence is recorded in the vector y.

The task will be to implement both the data generation part and the model estimation block, and then to compare the performance of two different models.

We are going to simulate the case where data originate from a system which can be assumed to have the following underlying physical model:
x(t) = e^{at}, (1)
where t is time (in seconds) and a is an unknown parameter. For this exercise we also assume that 0 < t < 10.
Any measurement will also contain noise, here assumed to be additive white Gaussian noise with zero mean and variance 2. This means that the measured data (or signal) y(t)
will be
y(t) = x(t) + n(t), (2)
where n(t)  N(0, 2) and x(t) is given by Eq. (1).
Throughout this work we also assume that the signals have been sampled with a sampling
time Ts, so that
t = 0, Ts, 2Ts, . . . , (N -1)Ts. (3)

The goal of this exercise is to implement and evaluate the different steps of the system described.
In MATLAB this means that all signals are represented by vectors, so that
x =  [x(0) x(Ts) x(2Ts) .... x((N - 1)Ts) T
(4)
is a column vectors of size N x 1.

Write a MATLAB function, that given the sampling time Ts and the noise standard deviation  returns x, y, and t. Set the parameter a in the model to a = 0.5.

Modify the function so that it returns K realizations of y, in a matrix Y, where each column corresponds to one realization of y.

Plot the noise-free signal x and one realization of the noisy signal y in the same plot, for a noise standard deviation  = 10 and a sampling time Ts = 0.05 s, for 0 < t < 10.
This should look something like Fig. 2.

We are now going to fit two different models to the data generated in Section 3. The first one is a simple polynomial model, which has the advantage that it is linear in the unknown parameters and hence the estimator has a closed-form solution. The other model is the true exponential model, as given by Eq. (1). This model has only one unknown parameter, but we need to use some iterative optimization technique to find it.

We are now assuming that x(t) can be approximated as
bx(t) = b0 + b_1 t + b_2 t² + . . . + b_p t_p, (5)
i.e. a polynomial of order p. As explained in [1], this can be expressed in matrix notation as
y = Ab, (6)
where y is an N ? 1 vector containing the measured data, A is an N x (p + 1) matrix, and b is the (p + 1) x 1 vector with the unknown polynomial coefficients.
The least-squares estimate of b is now given by
bLS = (A^{T}A)^{-1} *A^T y, (7)
And the modeled x, let?s call it bx is given by:
x = Ab_{LS}. (8)

Write a MATLAB function that, given the vector t and the polynomial model order p, creates the matrix A.
Write a MATLAB function that, given K realizations of the measured sequence y, returns the K estimates, bx, as columns of a matrix bY (similar to the assignment in
Section 3.2).

Plot some examples of the estimated polynomial together with the true signal x(t), for polynomial model orders p = 3, 4, and 5. This could look like in Fig. 3.

Now you will implement the Gauss-Newton method for solving the non-linear least-squares problem of finding a in Eq. (1), i.e. an iterative approach so that
baj+1 = baj + !HT (baj)H(baj)
&#8722;1
HT (baj) (y &#8722; x(baj)) , (9)
where j is the iteration number, baj is the estimate of the unknown model parameter a for iteration j. H(bak) is the model gradient, and x(baj) is the model, evaluated for the current value of baj .

Derive an expression for the gradient.
Write a MATLAB function that, given the time vector t and a value of a, returns a vector containing the gradient.
Implement the Gauss-Newton algorithm as a MATLAB function that returns the estimated value of a and the resulting model y(t).
Here is a pseudo-code example of the algorithm. All you need to do is to fill in the blanks and write it as proper MATLAB code. See [1] for more details.
j = 0; % Iteration number
a = 0.2; % Starting guess of a
update = inf; % Initial change of the parameter
thresh = 1e-10; % Stop criterion
while (j < max no. of iter.) and (update > thresh),
H = model_gradient(a)
da = inv(H?*H)*H?*(y-signal_model(t,a)
update = a;
a = a + da;
update = update - a; % Check how much a has changed
j = j + 1;
end
5

Now it is time to summarize and compare the two models. We know that the polynomial model is wrong, although the results seem quite good. The advantage of this model is primarily its simplicity in terms of parameter estimation. It is also fairly easy to evaluate its do some simulation-based evaluation. The other model, which in this case happens to be the correct one, is non-linear with respect to the unknown parameter, a. Because of this we used an iterative method to find the parameter. The questions we should consider now are:
1. Is the polynomial model just as good? If not, how?
2. How does the accuracy (or bias) of the two models compare?
3. How does the uncertainty (standard deviation) of the two models compare?
We will consider these questions in relation to:
1. Noise variance.
2. Sampling time (i.e. the number samples available for fitting the models).
3. Computational complexity.

Write a MATLAB script that, based on K = 300 realizations, returns the matrices
bYpoly and bYexp, containing the estimates using the polynomial and exponential models,
respectively.

Evaluate the repeatability (uncertainty) of the two models. Do it by generating and
discussing the following plots:
1. The average of the K estimated polynomial models and the ?2 interval, for a
sampling time Ts = 0.01 s.
2. Same as above, but for the exponential models.
3. Same as the two above, but now using a sampling time of Ts = 0.05 s. Describe
what happens?

Evaluate the model mismatch of the two models, by generating the following plots:
1. The average error as function of time, i.e.
Note: Plot the average error of both models in the same plot.

2. All model errors as function of time. Make one plot for the polynomial model and
one for the exponential model. Does the model mismatch appear to be approximately
the same over the whole time interval? If not, discuss what you see and
what you think causes this.

3. Repeat the two above, but now look at what happens if the polynomial order
p = 3, 4, or 5.