ACIM speed estimation

In case if we would like to know the induction motor rotor speed, we can use Extended Kalman Filter(EKF) to estimate the unknown states. Here is a good tutorial where you can find more information about the EKF. Also you can read the article on the Wikipedia. I'm using the FOC block in Matlab to get the motor control current, voltage and the angle between the fluxes, just for testing purpose.
I assume that the rotor speed is constant so $\frac{d\omega_r}{dt}=0$. In this case we have the following model

$$A= \begin{bmatrix} -\frac{R_s+\frac{L_m^2R_r}{L_r^2}}{L_s\sigma} & 0 &  \frac{L_m}{L_s\sigma L_r\tau_r} & \frac{\omega_rL_m}{L_s\sigma L_r} & 0 \\ 0 & -\frac{R_s+\frac{L_m^2R_r}{L_r^2}}{L_s\sigma}  & -\frac{\omega_rL_m}{L_s\sigma L_r} & \frac{L_m}{L_s\sigma L_r\tau_r} & 0 \\ \frac{L_m}{\tau_r} & 0 & -\frac{1}{\tau_r} & -\omega_r & 0 \\ 0 & \frac{L_m}{\tau_r} &   \omega_r & -\frac{1}{\tau_r}  & 0 \\ 0 & 0 & 0 & 0 & 0\end{bmatrix}\tag{1}$$

$$B=\begin{bmatrix} \frac{1}{L_s\sigma} & 0 \\ 0 & \frac{1}{L_s\sigma} \\ 0 & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix}\tag{2}$$

We can discretize the model using the following approach:

$$A_d=e^{AT} \approx I+AT\tag{3}$$
$$B_d=\int_0^T e^{At}B dt \approx BT\tag{4}$$

The Simulink model is the following

In the block above, the second and third module are the Clarke's and Park's transformation, in the previous articles I already defined this transformations. The code of the speed and flux estimator is

function [xx_out,PP_out,error,W]=ekf(I,U,Ts)
% EKF   Extended Kalman Filter  
% induction motor's parameters   
persistent PP;
persistent xx;

if isempty(PP)
    xx = [0 0 0 0 0 ]';
    %[i_alfa , i_beta , Fi_r_alfa , Fi_r_beta , W_r ]
    PP = 1e-2*diag([1 1 1 1 1]);
end
Q =  diag([1e-4 1e-4 1e-9 1e-9 5e-4]);
R =  1e1*diag([1 1]);

L_r=0.479;
L_s=0.423;
L_m=0.421;
R_s=5.27;
R_r=5.07;

tau_r = L_r/R_r;
gamma_r = 1-(L_m^2)/(L_s*L_r);
T_r = R_s+(L_m^2*R_r)/L_r^2;
 
Ak =[ 1-Ts*(T_r/(L_s*gamma_r)),          0,        Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),        0;
    0,                  1-Ts*(T_r/(L_s*gamma_r)), -Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),     Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        0;
    Ts*(L_m/tau_r),                  0,                   1-Ts*(1/tau_r),                  -Ts*xx(5),                                    0;
    0,                      Ts*(L_m/tau_r),                  Ts*xx(5),                      1-Ts*(1/tau_r),                         0;
    0                           0                           0                    0                                           1];

Bk = [Ts/(L_s*gamma_r)              0;
            0                Ts/(L_s*gamma_r);
            0                       0;
            0                       0;
            0                       0];

%Predictor
xx = Ak*xx+Bk*U;


Gk =[ 1-Ts*(T_r/(L_s*gamma_r)),          0,        Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),        Ts*((xx(4)*L_m)/(L_s*gamma_r*L_r));
    0,                  1-Ts*(T_r/(L_s*gamma_r)), -Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),     Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        -Ts*((xx(3)*L_m)/(L_s*gamma_r*L_r));
    Ts*(L_m/tau_r),                  0,                   1-Ts*(1/tau_r),                  -Ts*xx(5),                                    -Ts*xx(4);
    0,                      Ts*(L_m/tau_r),                  Ts*xx(5),                      1-Ts*(1/tau_r),                         Ts*xx(3);
    0                           0                           0                           0                                           1];

H = [1 0 0 0 0;
     0 1 0 0 0;];

PP=Gk*PP*Gk'+Q;

%Coorection
KalG=PP*H'*inv(H*PP*H'+R);
error = I-H*xx;
xx=xx+KalG*(error);
PP=(eye(5)-KalG*H)*PP;

xx_out = xx;
PP_out=PP;

W = xx(5);

I spent a lot of time to tune the filter matrices Q,R, but this is the best result what I got.
I wasn't satisfied with the result above so I extended this model with the torque, in this case the rotor speed can be expressed as follows:
$$\frac{d\omega_r}{dt} =\frac{3p^2L_m}{8JL_r}\left(i_{sq}\psi_{rd}-i_{sd}\psi_{rq}\right)-\frac{p}{2J}T_l\tag{5}$$

we assume that the torque is constant or the changes between two sampling time is very small

$$\frac{dT_l }{dt}=0\tag{6}$$
In this case the EKF's code is

function [xx_out,PP_out,error,W]=ekf(I,U,Ts)
% EKF   Extended Kalman Filter  
% induction motor's parameters   
persistent PP;
persistent xx;

if isempty(PP)
    xx = [0 0 0 0 0 0]';
    %[i_alfa , i_beta , Fi_r_alfa , Fi_r_beta , W_r ]
    PP = 1e-8*diag([1 1 1 1 1 1]);
end
Q =  diag([1e-3 1e-3 1e-9 1e-9 1e-5 1e-5]);
R =  1e-3*diag([1 1]);

L_r=0.479;
L_s=0.423;
L_m=0.421;
R_s=5.27;
R_r=5.07;
p = 2;
J =0.02;
fv = 0.0038;

m = (3*p^2*L_m)/(8*J*L_r);

tau_r = L_r/R_r;
gamma_r = 1-(L_m^2)/(L_s*L_r);
T_r = R_s+(L_m^2*R_r)/L_r^2;
 
Ak =[ 1-Ts*(T_r/(L_s*gamma_r)),          0,        Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),             0       0;
    0,                  1-Ts*(T_r/(L_s*gamma_r)), -Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),     Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),                0       0;
    Ts*(L_m/tau_r),                  0,                   1-Ts*(1/tau_r),                  -Ts*xx(5),                                       0       0;
    0,                      Ts*(L_m/tau_r),                  Ts*xx(5),                      1-Ts*(1/tau_r),                                 0       0;
    -Ts*m*xx(4)                  Ts*m*xx(3)                       0                               0                                   1           (-p*Ts)/(J*2);
           0                    0                               0                               0                                   0               1];

Bk = [Ts/(L_s*gamma_r)              0;
            0                Ts/(L_s*gamma_r);
            0                       0;
            0                       0;
            0                       0;
            0                       0];

%Predictor
xx = Ak*xx+Bk*U;


Gk =[ 1-Ts*(T_r/(L_s*gamma_r)),          0,        Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),        Ts*((xx(4)*L_m)/(L_s*gamma_r*L_r))   0;
    0,                  1-Ts*(T_r/(L_s*gamma_r)), -Ts*((xx(5)*L_m)/(L_s*gamma_r*L_r)),     Ts*(L_m/(L_s*gamma_r*L_r*tau_r)),        -Ts*((xx(3)*L_m)/(L_s*gamma_r*L_r))     0;
    Ts*(L_m/tau_r),                  0,                   1-Ts*(1/tau_r),                  -Ts*xx(5),                                    -Ts*xx(4)                          0;
    0,                      Ts*(L_m/tau_r),                  Ts*xx(5),                      1-Ts*(1/tau_r),                                 Ts*xx(3)                        0;
    -Ts*m*xx(4)              Ts*m*xx(3)                       Ts*m*xx(2)                       -Ts*m*xx(1)                                      1                           (-p*Ts)/(J*2);
    0                               0                           0                           0                           0                                           1];                         
             
H = [1 0 0 0 0 0;
     0 1 0 0 0 0;];

PP = Gk*PP*Gk'+Q;


KalG = PP*H'*inv(H*PP*H'+R);
error = I-H*xx;
xx = xx + KalG*(error);
PP = (eye(6)-KalG*H)*PP;

xx_out = xx;
PP_out = PP;

W = xx(5);

The result is

The error between the measured and the estimated is smaller than the previous model, also it was easier to tune the estimator. The disadvantage of the model, that is has six state which means that the calculations takes longer than the previous model.