Recursive Least Squares Estimation for Dynamic Systems in Real time

11 views (last 30 days)
Hello. Actually I am looking to implement and simulate Recursive Least Squares Estimation with White noise of zero mean and unit variance. I have written the following codes:
rlsupdate.m
function theta_hat = rlsupdate(data)
% RLS
% MC 02/05/2020
%
persistent P theta phi
%
if isempty(P)
P = eye(2);
theta = 0;
phi = 0;
end
u = data(1);
y = data(2);
% recursive equations
K = P*phi*inv(eye(1)+phi'*P*phi);
theta = theta + K*(y-(phi'*theta));
P = P - P*phi*inv((eye(1) + (phi'*P*phi)))*phi'*P;
phi = [y u]';
return
h = 1.0; % sample period [s]
a = 0.9;
b = 0.1;
uA = 1.0; % amplitude of input (square)
uP = 75*2; % period of input (square) [s]
data = [a b];
rlsupdate(data)
After running the above code, I am getting the following warning:
Warning: Matrix is singular to working precision.
> In rlsupdate (line 15)
In rlsparams (line 8)
Warning: Matrix is singular to working precision.
> In rlsupdate (line 17)
In rlsparams (line 8)
Below is my Simulink model:
Can anyone please tell me why this warning is coming and how to get rid of it.

Answers (1)

Abolfazl Chaman Motlagh
Abolfazl Chaman Motlagh on 2 Mar 2022
The first time rlsupdate is called the value of P,phi are set to eye(2) and 0. in the line 15 of code you have inverse of eye(1)+phi'*P*phi. the eye(1) is just 1, so it becomes :
which is singular matrix (doesn't have inverse). ()
so the K become NaN. and after that P become NaN in line 17. and all goes wrong after that.
you only need to handle this problem in line 15 to solve this issue.
i think you mean:
K = P*phi*inv(eye(2)+phi'*P*phi);
theta = theta + K*(y-(phi'*theta));
P = P - P*phi*inv((eye(2) + (phi'*P*phi)))*phi'*P;
but there will be another problem after doing this. not for first time. first time it will run without problem. after first time because of line 20 :
phi = [y u]';
phi become 2x1 vector. but the calculation in line 15,16 and 17 are for scalar phi. you shoud handle them correctly.
maybe post the equation for K,theta and P so that i can help.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!