0

我正在尝试模拟系统的旋转动力学。我正在测试我的代码以验证它是否正在使用模拟工作，但我从未恢复过传递给模型的参数。换句话说，我无法重新估计我为模型选择的参数。

我为此使用 MATLAB，特别是ode45. 这是我的代码：

% Load the input-output data
[torque outputs] = DataLogs2(); 
u = torque;

% using the simulation data
Ixx = 1.00;
Iyy = 2.00;
Izz = 3.00;

x0          = [0; 0; 0];
Ts          = .02;
t           = 0:Ts:Ts * ( length(u) - 1 );
[ T, x ]    = ode45( @(t,x) rotationDyn( t, x, u(1+floor(t/Ts),:), Ixx, Iyy, Izz), t, x0 );

w           = x';
N           = length(w);
q           = 1;           % a counter for the A and B matrices

% The Algorithm
for k=1:1:N
    w_telda    = [   0      -w(3, k)    w(2,k); ...
                   w(3,k)       0      -w(1,k); ...
                  -w(2,k)    w(1,k)       0  ];

    if k == N        % to handle the problem of the last iteration
        w_dash(:,k)  = (-w(:,k))/Ts;
    else
        w_dash(:,k)  = (w(:,k+1)-w(:,k))/Ts;
    end

    a            = kron( w_dash(:,k)', eye(3) ) + kron( w(:,k)', w_telda );
    A(q:q+2,:)   = a;           % a 3N*9 matrix
    B(q:q+2,:)   = u(k,:)';     % a 3N*1 matrix   % u(:,k)

    q = q + 3;
end

% Forcing J to be diagonal. This is the case when we consider our quadcopter as two thin uniform 
% rods crossed at the origin with a point mass (motor) at the end of each.
A_new                 = [A(:, 1) A(:, 5) A(:, 9)];
vec_J_diag            = A_new\B;
J_diag                = diag([vec_J_diag(1), vec_J_diag(2), vec_J_diag(3)])
eigenvalues_J_diag    = eig(J_diag) 
error                 = norm(A_new*vec_J_diag - B)

我的动态模型定义为：

function [dw, y] = rotationDyn(t, w, tau, Ixx, Iyy, Izz, varargin)

    % The output equation
    y = [w(1); w(2); w(3)];

    % State equation
    % dw = (I^-1)*( tau - cross(w, I*w) );
    dw = [Ixx^-1 * tau(1) - ((Izz-Iyy)/Ixx)*w(2)*w(3);
          Iyy^-1 * tau(2) - ((Ixx-Izz)/Iyy)*w(1)*w(3);
          Izz^-1 * tau(3) - ((Iyy-Ixx)/Izz)*w(1)*w(2)];
    end

实际上，这段代码应该做的是计算惯性矩阵 的特征值J，即恢复Ixx, Iyy，并且Izz我在一开始就传递给模型（1、2 和 3），但我得到的都是错误的结果。

使用有问题ode45吗？

4

1 回答 1

0

那么问题不在ode45指令中，问题在于在系统识别中可以n-1从样本信号创建样本n信号，因此循环必须在上述代码中的 N-1 处结束。

于 2015-10-28T12:28:27.533 回答