Using ODE1 in Simulink yields different results for an equivalent model run in a for loop
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
I've run into an interesting situation with Simulink. I am using Simulink to evaluate alternative numerical integration routines for some time-varying ordinary differential equations. The Simulink model that I have built is provided in the attached screen shot. The block, quaterionKinematics, encodes the continuous-time, time-varying differential equation that governs the evolution of quaternion kinematics. As you can see, the output of this block is the time derivative of the quaternions, we integrate to get q(t), and then perform a normalization step.
function q_dot = quaternionKinematics(wx, wy, wz, q)
q_dot = 0.5*[0, -wx, -wy, -wz;...
wx, 0, wz, -wy;...
wy, -wz, 0, wx;...
wz, wy, -wx, 0]*q;
end
The gyroscope data, given by the wx, wy, wz is loaded from my workspace. In evaluating these results against alternative integraiton schemes that I will have to implement in another programming language, I have built the equivalent for loops that perform this same integration. Equivalent for loop:
%% Euler Integration
h = Ts; %sample period
Omega(:,:,1) = 0.5*[0 -wx(1) -wy(1) -wz(1);...
wx(1) 0 wz(1) -wy(1);...
wy(1) -wz(1) 0 wx(1);...
wz(1) wy(1) -wx(1) 0];
q_last = initQuat;
Omega_last = Omega(:,:,1);
for i = 2:length(time_i)
Omega(:,:,i) = 0.5*[0 -wx(i) -wy(i) -wz(i);...
wx(i) 0 wz(i) -wy(i);...
wy(i) -wz(i) 0 wx(i);...
wz(i) wy(i) -wx(i) 0];
q_Euler(:,i) = q_last + h*Omega_last*q_last;
q_Euler(:,i) = q_Euler(:,i)/norm(q_Euler(:,i),2);
q_last = q_Euler(:,i);
Omega_last = Omega(:,:,i);
end
When using ODE1, which is the forward Euler method that I have depicted above, my results disagree around certain points in my data...and by quite a bit. I have included a figure that shows the difference between one of the propagated quaternion elements.
When I move toward higher-order integration schemes, such as ODE2 (Heun), the output of Simulink tends to track the output of what should be an equivalent for loop a bit better, however, there is still inconsitency between the two outputs, which I would have thought should be equivalent. The most significant disagreement seems to be specific to the results that I am getting with ODE1.
Can you provide some guidance on why this may be? I am hesitant to use the results that I am obtaining from Simulink in making design decisions, now. Could this be related to how I am encoding my differential equations, where I have used an embedded Matlab function? This approach has not cuased problems for me in the past, and again, gross disagreement in results really only seem to be specific to use of ODE1 as my solver.
3 commentaires
Réponses (1)
Paul
le 10 Juin 2021
I think some of the equations are out of order. I got a match (obviously using my inputs) with this:
h = Ts; %sample period
q_Euler = zeros(4,numel(out.tout));
q_Eulerdot = zeros(4,numel(out.tout));
q_Euler(:,1) = out.qsim(1,:).'; % out.qsim is the normalized quaternion, output from simulink model
q_int = out.qsim1(1,:).'; % out.qsim1 is the output of the integrator (non-normalized quaternion) output from simulink model
for i = 1:length(out.tout)
q_Euler(:,i) = q_int/norm(q_int); % normalize the integrator output
% compute Omega with wx, wy wz outputs from Simulink model
Omega = 0.5*[0 -out.wx(i) -out.wy(i) -out.wz(i);...
out.wx(i) 0 out.wz(i) -out.wy(i);...
out.wy(i) -out.wz(i) 0 out.wx(i);...
out.wz(i) out.wy(i) -out.wx(i) 0];
q_Eulerdot(:,i) = Omega*q_Euler(:,i);
q_int = q_int + h*q_Eulerdot(:,i); % propagate integrator to next time step
end
It will be important that the quat_normalization function implements the /norm() as in the code above.
Historical note: a long time ago I ran into a similar problem where Simulink didn't exactly match the Matlab implementation, and was told by Tech Support that Simulink and Matlab could be linked with different libraries, so innoccuous operations like matrix multiply might not even give the exact same result. Don't know if that's still true.
0 commentaires
Voir également
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!