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

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,:)';

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);


hold on
hold off

function f = take_value_of_input(t)

    f = 5*sin(3*t);


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


$θ = [θ_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.

  • $\begingroup$ 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? $\endgroup$
    – spektr
    May 21, 2019 at 1:59
  • $\begingroup$ 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? $\endgroup$ May 21, 2019 at 7:22
  • $\begingroup$ 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? $\endgroup$ May 21, 2019 at 8:06
  • $\begingroup$ But the equation for the estimates is: dθ1/dt = γ*(x - θ1*φ1)*φ1 and dθ2/dt = γ*(x - θ2*φ2)*φ2 . $\endgroup$ May 21, 2019 at 8:11
  • $\begingroup$ 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? $\endgroup$ May 21, 2019 at 8:19

1 Answer 1


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)]; 

  • $\begingroup$ Works like a charm!! Thank you very very much. Your help was incredible!! Except from the working solution, I also really learnt a lot !! $\endgroup$ May 21, 2019 at 12:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.