Clear Filters
Clear Filters

Recursive Least Squares Estimation for Dynamic Systems in Real time

8 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.

Categories

Find more on Simulink in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!