djdgf4 发表于 2023-12-30 23:56

关于惯性导航和四元数的matlab代码,是否存在不正确的参数数据类型或缺少参数

这是一份实验,

要求:完成姿态四元数的更新(运载体为汽车,沿地面运动,高程h=0)。
已知:上一时刻的姿态四元数、上一时刻及当前时刻的角增量、当前纬度、当前速度。
1.构造一个函数,完成两个四元数的乘法。
按四元数乘法的原始定义来编写。
2.由地球转速、汽车速度及位置构造两个四元数
b系的更新四元数由等效旋转矢量构造。
n系的更新四元数构造

当我执行QuaternionScript.m文件的时候,报错
>> updateAttitudeScript
检查对函数
'updateQuaternionCarrier' 的调用中
是否存在不正确的参数数据类型或缺少参数。

出错
updateAttitudeScript>updateAttitude
(第 26 行)
    q_update_b =
    updateQuaternionCarrier(angularIncrement,
    angularIncrement_pri);

出错 updateAttitudeScript (第 10
行)
updatedQuaternion =
updateAttitude(initialQuaternion,
angularIncrement,
angularIncrement_pri, latitude,
velocity, timeIncrement);但是我确定没有缺少参数.不知道是为什么.

下面是我的代码:
C:\Users\name\Documents\MATLAB\QuanternionupdateAttitude.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Filename:      QuaternionRotation.m
%      
%Description:      要求:完成姿态四元数的更新(运载体为汽车,沿地面运动,高程h=0)。
%                     已知:上一时刻的姿态四元数、上一时刻及当前时刻的角增量、当前纬度、当前速度。
%                     1.构造一个函数,完成两个四元数的乘法。
%                     按四元数乘法的原始定义来编写。
%                     2.由地球转速、汽车速度及位置构造两个四元数
%      
%      Version:      1.0
%      Created:      13/29/23 19:36:24
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quaternion Class Definition
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% quaternion =
% also, quaternion =       
% _v repersent a vector, _v=
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
classdef Quaternion
      properties
      s
      x
      y
      z
      end
      
      methods
                function obj_quaternion = Quaternion(s,v_x,v_y,v_z)
                        obj_quaternion.s = s;
                        obj_quaternion.x = v_x;
                        obj_quaternion.y = v_y;
                        obj_quaternion.z = v_z;
                end
      
      function q_inv = inverse(obj)
            norm_squared = obj.s^2 + obj.x^2 + obj.y^2 + obj.z^2;
            q_inv = Quaternion(obj.s/norm_squared, -obj.x/norm_squared, -obj.y/norm_squared, -obj.z/norm_squared);
      end
      
      function q_conj = conjugate(obj)
            q_conj = Quaternion(obj.s, -obj.x, -obj.y, -obj.z);
      end
      
      function q_norm = norm(obj)
            q_norm = sqrt(obj.s^2 + obj.x^2 + obj.y^2 + obj.z^2);
      end
               
                % Define the rule of multiplity               
                function q_prod = times(obj,q)
                        q_prod = Quaternion(obj.s*q.s - obj.x*q.x - obj.y*q.y - obj.z*q.z,...
                                             obj.s*q.s + obj.x*q.x - obj.y*q.y + obj.z*q.z,...
                                             obj.s*q.s + obj.x*q.x + obj.y*q.y - obj.z*q.z,...
                                             obj.s*q.s - obj.x*q.x + obj.y*q.y + obj.z*q.z);
                end
               
                function q_div = mrdivide(obj, q)
            q_div = obj.times(q.inverse());
      end
      
      function q_sum = plus(obj, q)
            q_sum = Quaternion(obj.s + q.s, obj.x + q.x, obj.y + q.y, obj.z + q.z);
      end
      
      function q_diff = minus(obj, q)
            q_diff = Quaternion(obj.s - q.s, obj.x - q.x, obj.y - q.y, obj.z - q.z);
      end
      
      function disp(obj)
            fprintf('%.2f + %.2fi + %.2fj + %.2fk\n', obj.s, obj.x, obj.y, obj.z);
      end
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Quaternion rotation function $$Q^R_b$$
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % the input theta is the rotation angle (right-hand system)
      % u is the rotation axis vector
      % v is the original point position vector
      % the output V is the rotated point position vector
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function V = rotate(theta,u,v)
            u=u./sqrt(sum(u.*u));                                        % Normalized
            q=Quaternion(cos(theta/2),u(1)*sin(theta/2),u(2)*sin(theta/2),u(3)*sin(theta/2));
            qt=conjugate(q);                                                % Take the conjugate
            w=Quaternion(0,v(1),v(2),v(3));                        % Pure quaternion
            W = q.times(w).times(qt);               
            V=;                              
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % update (carrier coordinate system) function
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % Phi_k is equivalent rotation vector corresponding to the attitude change in       -- $$\Phi_k$$
      % q_new_k_B is the quaternion corresponding to the attitude change in       -- $$q^{b(k-1)}_{b{k}}$$
      % angularIncrement : 当前时刻的角增量
      % angularIncrement_pri : 前一时刻的角增量
      % Phi_k : 双子样计算的等效旋转矢量
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function q_update_B = updateQuaternionCarrier(angularIncrement,angularIncrement_pri)
            % Twin sample calculation method
            Phi_k = angularIncrement + cross(angularIncrement,angularIncrement_pri) / 12;
            length_norm = norm(0.5*Phi_k);
            % create the quaternion corresponding to the attitude change in
            q_update_B = Quaternion(cos(length_norm),...
                     sin(length_norm)/length_norm * 0.5 * Phi_k(1),...
                     sin(length_norm)/length_norm * 0.5 * Phi_k(2),...
                     sin(length_norm)/length_norm * 0.5 * Phi_k(3));
      end

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % update (Navigation coordinate system) function
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % O - NUE
      % latitude 纬度
      % velocity : 速度
      % V_theta : 速度偏向角
      % timeIncrement : 时间增量
                % h : 高程
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      function q_update_N = updateQuaternionNavigation(latitude, velocity, timeIncrement, V_theta, h)
                        if nargin < 4
                              V_theta = 0;
                        end
                        
                        if nargin < 5
                              h = 0;
                        end
                        
            % Earth's rotation angular speed
            w_earth = 2 * pi / 24 / 3600;               
            % 地球椭球长半轴和第一偏心率
            a = 6378136.49;
            e = 0.0066943;

            V_E = velocity*cos(V_theta);
            V_N = velocity*sin(V_theta);
            R_M = (a*(1-e*e)) / power( 1-e*e*power(sin(latitude),2) , 3/2 );
            R_N = a / power( 1-e*e*power(sin(latitude),2) , 1/2 );
            w_nie = ;
            w_nen = [V_E/(R_N+h);
                     -V_N/(R_M+h);
                     -V_E*tan(latitude)/(R_N+h)];
            zeta_k = (w_nen+w_nie)*timeIncrement;
            length_norm = norm(0.5*zeta_k);
            % create the quaternion corresponding to the attitude change in
            q_update_N = Quaternion(cos(length_norm),...
                     -sin(length_norm)/length_norm * 0.5 * zeta_k(1),...
                     -sin(length_norm)/length_norm * 0.5 * zeta_k(2),...
                     -sin(length_norm)/length_norm * 0.5 * zeta_k(3));
      end
      end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End Define
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


C:\Users\name\Documents\MATLAB\Quanternion\QuaternionScript.m
% 定义输入参数
initialQuaternion = Quaternion(1, 0, 0, 0); % 初始姿态四元数
angularIncrement = ;      % 角增量
angularIncrement_pri = ;       %前一时刻的角增量
latitude = 30;                           % 当前纬度
velocity = 10;                           % 当前速度
timeIncrement = 0.01;                     % 时间增量

% 调用 updateAttitude 函数
updatedQuaternion = updateAttitude(initialQuaternion, angularIncrement, angularIncrement_pri, latitude, velocity, timeIncrement);

% 查看结果
disp('Updated Quaternion:');
disp(updatedQuaternion);

function updatedQuaternion = updateAttitude(initialQuaternion, angularIncrement, angularIncrement_pri, latitude, velocity, timeIncrement)
    % initialQuaternion: 上一时刻的姿态四元数
    % angularIncrement: 当前时刻的角增量
    % angularIncrement_pri: 前一时刻的角增量
    % latitude: 当前纬度
    % velocity: 当前速度
    % timeIncrement: 时间增量
   
    % 步骤二:构造更新四元数
    % a. 汽车运动四元数的构造
    q_update_b = updateQuaternionCarrier(angularIncrement, angularIncrement_pri);
   
    % b. 地球自转四元数的构造
    q_update_n = updateQuaternionNavigation(latitude, velocity, timeIncrement,0,0);
   
    % 步骤三:更新姿态四元数
    % 计算更新四元数
    q_update = q_update_n.times(q_update_b.conjugate());

    % 计算最终姿态四元数
    updatedQuaternion = q_update.times(initialQuaternion);

end
页: [1]
查看完整版本: 关于惯性导航和四元数的matlab代码,是否存在不正确的参数数据类型或缺少参数