Using the nomenclature presented in your Navigation Systems course (mathematics section), the vehicle body frame () is the object frame while the reference frame () is defined at some fixed point on the earth's surface as shown. The position of the vehicle is then denoted , where γ can be either or b. In this scenario, the vehicle velocity (), defined with respect to the body axes, yields an expression for the velocity observed in the object frame as,
By computing the appropriate velocity and acceleration vectors, write a MATLAB script that will track the position of the vehicle, resolved in earth axes, with the following control inputs:
assuming constant acceleration over each timestep, plot the resulting trajectory over 120s (2min). Use the integration schemes defined in the course notes to solve the ODE's for and , i.e. trapezoidal integration and the fourth-order scheme for .
以下是Matlab代码：
%% Define simulation parameters
dt = 0.01; % Simulation timestep (s)
t = 0:dt:120; % Simulation time points
om = 0.2; % Heading modulation frequency (rad/s)
V = 1.5; % Vehicle velocity
k = 0.05; % Heading modulation amplitude (rad)
om_ie_z = 7.292115e-5; % Earth spin rate (rad/s)
om_ie = [0;0;om_ie_z];
Om_ie = [0 -om_ie_z 0;om_ie_z 0 0;0 0 0];
%% Initialise variables
r_eb_e = zeros(3,length(t)); % Start at the origin
v_eb_e = zeros(3,length(t)); % Pre-allocate memory
vd_eb_e = zeros(3,length(t)); % Pre-allocate memory
om_ib_b_measured = zeros(3,length(t)); % Pre-allocate memory
f_ib_b_measured = zeros(3,length(t)); % Pre-allocate memory
psi0 = -0.24; % Initial vehicle Heading (rad)
Cbe = [cos(psi0) sin(psi0) 0;-sin(psi0) cos(psi0) 0;0 0 1]';
% Body-to-Earth co-ord transformation matrix
Cbe_est = Cbe; % Body-to-Earth co-ord transformation matrix (estimate)
u = [V;0;0]; % Absolute velocity resolved in body axes
v_eb_e(:,1) = Cbe*u; % Initialise estimated velocity
%% Calculate the true trajectory
% In this section, the true vehicle trajectory will be calculated using a
% fourth-order numerical integration scheme. Therefore, the true vehicle
% trajectory will be contained in the 'y' matrix.
%
y0 = [0;0;V*cos(psi0);V*sin(psi0);V;psi0;k*sin(psi0)];
[tx,y] = ode45(@(t,y)calcDerivatives(t,y,k,om,psi0),t,y0);
%% Define the error terms
errSwitch = 0; % Set this value to '1' to add constant errors and '0' otherwise
noiseSwitch = 0; % Set this value to '1' to add noise and '0' otherwise
noisea = 0.0001*noiseSwitch*rand(3,length(t)); % Accelerometer white noise
noiseg = 0.00001*noiseSwitch*rand(3,length(t)); % Gyro white noise
ba = 1e-4*[1;1;1]*errSwitch; % Accelerometer bias
sa = 1e-4*diag([1;1;1]); % Accelerometer scale factor error
ma = [0 -0.0001 0.0003;0.0001 0 0.0002;-0.0001 -0.0002 0]; % Acc. cross coupling error
Ma = (sa + ma)*errSwitch; % Combined accelerometer scale factor & cross coupling error matrix
bg = 1e-5*[1;1;1]*errSwitch; % Gyro bias
sg = 1e-5*diag([1;1;1]); % Gyro scale factor error
mg = [0 -0.0001 0.0003;0.0001 0 0.0002;-0.0001 -0.0002 0]; % Gyro cross coupling error
Mg = (sg + mg)*errSwitch; % Combined gyro scale factor & cross coupling error matrix
%% Loop over time
for i=2:length(t)
%
% ***********************************************
% * Calculate the estimated INU trajectory *
% ***********************************************
%
% First step, we need to calculate the true values of specific force
% and angular rate. These will be used to generate our simulated
% measurements.
% First, find the true angular rate
psi = psi0 + om*t(i);
om_ib_b_True = [0;0;k*sin(psi)] + Cbe_est*om_ie;
Om_ib_b_True = [0 -om_ib_b_True(3,1) 0;...
om_ib_b_True(3,1) 0 0;0 0 0];
% Next, calculate the true specific force
f_ib_b_True = (Om_ib_b_True + Cbe_est'*Om_ie*Cbe_est)*u;
% Add the errors onto the true angular rate and specific force vectors
% to create the measurements. (Hint: add the gyro random noise term
% using the following notation - noiseg(:,i). Sim for accelerometer
% noise). First, the gyro measurement
I=eye(3);
om_ib_b_measured(:,i) = bg+(I+Mg)*om_ib_b_True+noiseg(:,i);
Om_ib_b_measured = [0 -om_ib_b_measured(3,1) 0;...
om_ib_b_measured(3,1) 0 0;0 0 0]; % Skew-symmetric matrix form
% Now the accelerometer measurement
f_ib_b_measured(:,i) = ba+(I+Ma)*f_ib_b_True+noisea(:,i);
% Although not explicitly required for this formulation, calculate the
% strapdown equation in here i.e. find Cbd_dot_est.
Cbe_dot_est = Cbe_est*Om_ib_b_True-(Om_ie*Cbe);
% Navigation Equation goes in here
vd_eb_e(:,i) = Cbe*f_ib_b_measured(:,i)-2*Om_ie*v_eb_e(:,(i-1));
% Now integrate all equations using the techniques present in the
% notes i.e. fourth-order scheme for Cbe_est and trapezoidal
% integration for acceleration and velocity
% First, the DCM using the 4th order scheme. (Hint: when calculating
% alpha, remember to use om_eb_b_measured).
alpha = (om_ib_b_measured(:,i)-om_ie)*dt;
alpha_norm = norm(alpha);
if(alpha_norm > 0) % only update if there has been a rotation
skew_alpha = [0 alpha(3,1) 0;...
alpha(3,1) 0 0;0 0 0]; % Skew-symmetric matrix form
Cbe_est = Cbe;
end
% Finally, integrate the velocity and position ODE's using the
% trapezoidal rule. First, calculate the new velocity estimate,
v_eb_e(:,i) = v_eb_e(:,(i-1))+((vd_eb_e(:,i)+vd_eb_e(:,(i-1)))/2)*dt ;
% and the position estimate,
r_eb_e(:,i) = r_eb_e(:,(i-1))+((v_eb_e(:,i)+v_eb_e(:,(i-1)))/2)*dt;
% The following lines of code are used to store data for use in the assessments below.
if(i==2)
om_ib_b_check = om_ib_b_measured(:,i);
f_ib_b_check = f_ib_b_measured(:,i);
Cbe_dot_check = Cbe_dot_est;
vd_eb_e_check = vd_eb_e(:,i);
alpha_check = alpha;
Cbe_check = Cbe_est;
v_eb_e_check = v_eb_e(:,i);
r_eb_e_check = r_eb_e(:,i);
end
end
%
% Plot the results
figure(1);
% Plot the track
plot3(y(:,1),y(:,2),zeros(1,length(tx)),'b');
hold on;
plot3(r_eb_e(1,:),r_eb_e(2,:),-r_eb_e(3,:),'g');
xlabel('x (m)');
ylabel('y (m)');
zlabel('z (m)');
title('Vehicle trajectory');
legend('true','estimated','Location','NorthWest');
hold off;
% axis equal;
view(2);
grid on;
% Plot the error magnitude
figure(2);
plot(t,y(:,1)-r_eb_e(1,:)');
hold on;
plot(t,y(:,2)-r_eb_e(2,:)');
hold off;
xlabel('t (s)');
ylabel('error (m)');
title('Navigation errors vs. time');
legend('x-error','y-error','Location','SouthWest');
grid on;
%
figure(3);
plot(t,om_ib_b_measured(1:3,:));
xlabel('t(s)');
ylabel('omega_i_b^b (rad/s)');
title(' Measured Angular rates omega_i_b^b');
legend('p','q','r','Location','NorthWest');
grid on;
figure(4);
plot(t,f_ib_b_measured(1:2,:));
xlabel('t(s)');
ylabel('f_i_b^b (m/s^2)');
title('measured Specific force f_i_b^b');
legend('x-axis','y-axis','Location','SouthWest');
grid on;
function yd = calcDerivatives(t,y,k,beta,psi0)
%
yd = y;
yd(1,1) = y(5,1)*cos(y(6,1)); % xdot
yd(2,1) = y(5,1)*sin(y(6,1)); % ydot
yd(3,1) = -y(5,1)*sin(y(6,1))*y(7,1); % xddot
yd(4,1) = y(5,1)*cos(y(6,1))*y(7,1); % yddot
yd(5,1) = 0; % Vdot
yd(6,1) = k*sin(beta*t+psi0); % psidot
yd(7,1) = k*beta*cos(beta*t+psi0); % psiddot
end
其中的errSwitch和noiseSwitch可以设置为1，请尝试当errSwitch=1以及errSwitch和noiseSwitch都为1时图像的变化。
