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
$$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 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.
No comments:
Post a Comment