# Online Parameter Estimation using steepest descent

I have a first order system which is described by the following differential equation:

dx/dt = -a*x + b*u where u is the input u = 5*sin(3*t). In order to get the states x(t) of the system we are given that a = 2 & b = 1. I have written the following matlab code in order to estimate the parameters but results are not as expected. Identification error is not converging to zero as it is supposed to do and I can't figure out why. Any help would be really appreciated.

close all
clear all
clc

global am
global f
global estimator_gain
f = @take_value_of_input;
am = 2; % first order filter (s + am)
estimator_gain = 1;

tspan = 0:0.1:100;

[T,XY] = ode45(@sys, tspan, [0 0 0 0 0]);

state = XY(:,1);

filter = [XY(:,2) XY(:,3)];
theta_estimated = [XY(:,4) XY(:,5)];

x_estimated = zeros(length(tspan),1);

for i = 1:1:length(tspan)
x_estimated(i,1) = theta_estimated(i,:) * filter(i,:)';
end

error = state - x_estimated;

a_estimated = am - theta_estimated(end,1);
b_estimated = theta_estimated(end,2);

a_precision = (a_estimated / 2 ) * 100;
b_precision = (b_estimated / 1) * 100;

error_a = (am - 2) - (am - theta_estimated(end,1));
error_b = 1 - theta_estimated(end,2);

fprintf('\nEstimated values of parameters are:\n\n a = %.2f --> %.1f%% precision \n\n b = %.2f --> %.1f%% precision\n\n', a_estimated, a_precision, b_estimated, b_precision)

fprintf('Errors are:\n a : %.2f\n b : %.2f\n\n', error_a, error_b);

figure
plot(T,error)

figure
plot(T,state)
hold on
plot(T,x_estimated)
legend('state','x-estimated')
hold off

function f = take_value_of_input(t)

f = 5*sin(3*t);

end


And the function which solves the differential equations is:

function dxy = sys(t,xy)

global f
global am
global estimator_gain

u = f(t);

dxy = [-2*xy(1) + u ; xy(1) - am*xy(2) ; u - am*xy(3) ; estimator_gain*(xy(1) - xy(4)*xy(2))*xy(2) ; estimator_gain*(xy(1) - xy(5)*xy(3))*xy(3)];

% estimator_gain*(xy(1) - xy(4)*xy(2))*xy(2) --> This is the equation for θ1
% stimator_gain*(xy(1) - xy(5)*xy(3))*xy(3) --> This is the equation for θ2
% xy(1) - am*xy(2) --> Equation for φ1
% u - am*xy(3) --> Equation for φ2

end


$$θ = [θ_1~ θ_2]$$ is the vector which contains the estimates and $$φ = [φ_1 ; φ_2]$$ is the vector which contains the filter coefficients. The differential equation which describes the parameter estimation procedure is: $$\frac{dθ}{dt} = γ\,(x - θ^\top\,φ)\,φ$$ where $$x$$ is the computed state of the system.

• You need to explain what the other equations are that you're simulating. What's the filter variable and $\theta$ variable and how do they help your estimator converge to the true value? What's the basic theory behind it? – spektr May 21 '19 at 1:59
• Could you highlight where the gradient descent happens? What is the functional that you are minimizing? As here there is only one integration happening, are you testing if the gradient at the exact solution is zero? – Lutz Lehmann May 21 '19 at 7:22
• Note that in your equation you have a scalar product $θ^Tφ$. This is not reflected in your equations, where in the last two terms the second factor should be both times (xy(1) - xy(4)*xy(2) - xy(5)*xy(3)). If you repair this, does that solve your problem? – Lutz Lehmann May 21 '19 at 8:06
• But the equation for the estimates is: dθ1/dt = γ*(x - θ1*φ1)*φ1 and dθ2/dt = γ*(x - θ2*φ2)*φ2 . – Teo Protoulis May 21 '19 at 8:11
• But that is not what the product θ′∗φ means. Only one set of your formulas can be correct. Can you give a link to a derivation of this method? – Lutz Lehmann May 21 '19 at 8:19

The situation seems to be: You have some input function $$x$$ which more or less follows a model $$\dot x=-ax+bu$$. There may be noise involved, so the values of $$x$$ are not exact, and simply computing difference quotients will in general not be close to the right side of the differential equation. To filter out the noise some averaging is required. This you do by applying an approximate integrating factor with some estimate $$a_m$$ for $$a$$, so that theoretically $$e^{a_mt}x(t)=\int_0^te^{a_ms}(\dot x(s)+a_mx(s))ds =\int_0^te^{a_ms}((a_m-a)x(s)+bu(s))ds\\ x(t)=(a_m-a)φ_1(t)+bφ_2(t)$$ where $$φ_1(t)=e^{-a_mt}\int_0^te^{a_ms}x(s)ds$$ and $$φ_2(t)=e^{-a_mt}\int_0^te^{a_ms}u(s)ds$$ are solutions of the differential equations $$\dot φ_1=-aφ_1+x\\ \dot φ_2=-aφ_2+u.$$ This accounts for the first three equations.

Now one wishes to find estimates $$(θ_1, θ_2)$$ for $$(a_m-a,b)$$ in a way that $$\frac12(θ_1φ_1+θ_2φ_2-x)^2$$ decreases for increasing $$t$$. This can be achieved by moving the current estimate along the negative gradient, giving $$\pmatrix{\dot θ_1\\\dot θ_2}=-γ(θ_1φ_1+θ_2φ_2-x)\pmatrix{φ_1\\φ_2}$$ for some dampening factor $$γ$$ (too small $$γ$$ - the convergence is slow, too large $$γ$$ - not enough averaging, might become numerically unstable, stiff).

To summarize, in the last two terms the second factor should be both times (xy(1) - xy(4)*xy(2) - xy(5)*xy(3)). A better human readable ODE function could be

function dy = sys(t,y)

x = y(1); phi=y(2:3); theta=y(4:5);
diff = theta'*phi-x
u = f(t);

dy = [-2*x + u ; x - am*phi(1) ; u - am*phi(2) ; -estimator_gain*diff*phi(1) ; -estimator_gain*diff*phi(2)];

end

• Works like a charm!! Thank you very very much. Your help was incredible!! Except from the working solution, I also really learnt a lot !! – Teo Protoulis May 21 '19 at 12:21