ECEF 6DOF block - integrator initial conditions

The body rate integrator in the ECEF 6DOF block appears to generate nonzero outputs at the first timestep (t = 0) with initial conditions set to zero for all rates. That is, the first element of the logged output for the rates signal is nonzero when all three rates have started out as zero.

Is this behavior intended for this instance of the integrator?

Answer 1

[Paul]

Hi Thomas,

Are you referring to 6DOF ECEF (Quaternion)?

What exactly do you mean by the "body rate integrator"? Are you logging a signal at the output of an integrator inside the block (if so, which one)? Or are you referring to the block output port labeled omega_b?

Are the initial conditions you're referrring to the block parameter "Initial body rotation rates [p,q,r]"?

A. I'm going to assume that the initial condition in question is that block parameter. Keep in mind that p,q, and r are the components of what the doc page calls omega_rel.

B. I'm going to assume that the logged output in question is the output port omega_b. Keep in mind that omega_b is the angular velocity relative to inertial space.

We have the equation:

omega_b = omega_e + omega_ned + omega_rel

where omega_e is the angular velocity of ECEF relative to inertial and omega_ned is angular velocity of the NED frame relative to ECEF (and I've left out some DCM multiplications for simplicity).

So even if omega_rel = 0 at t = 0, omega_b would be non-zero at t = 0 unless the initial velocity and position of the vehicle is specified such that omega_ned is the exact opposite of omega_e.

If I'm looking at the wrong block or if either assumption is incorrect, please add more specificity as to the exact block (with a link to its doc page), block parameter, and signal in question.

[Thomas]

Hi Paul,

I have a custom planet model with rotation rate set to zero, so omega_e is zero. The output signal being referred to is the [p,q,r] vector (w_b in the 6DOF block).

[Paul]

According to the doc page linked in the Answer, as I read it w_rel = [p,q,r]. The components of w_b are not named. So "[p,q,r] vector (w_b in the 6DOF block)" is still unclear.

If you have w_e = 0, then

w_b = w_ned + w_rel (1)

The block parameters allow one to set the initial conditons of w_rel directly via the "Initial body rotation rates [p,q,r]".

The initial conditions on w_ned are derived from other initial conditions set via block parameters.

Are you seeing equation (1) as not being satisfied?

[Thomas]

In my case w_ned is zero as well, so w_rel and w_b are equivalent. I am seeing a nonzero initial w_rel when I specify zero as the initial condition for each of P/Q/R.

[Thomas]

I think I've isolated the issue to the execution order of the components within the ECEF 6DOF block. The w_rel signal initializes to zero with ICs set to zero, but the w_b signal has slightly nonzero values since the 'pdot/qdot/rdot' calculation occurs before the signal is output for each timestep. I am unsure if this is a function of the integrator or of some other aspect of the block.

[Paul]

I'm looking at the implementation in the block

If block parameter "Initial velocity in body axes [U,v,w] " is set to [0 0 0]

and

If block parameter "Initial body rotation rates [p,q,r]" is set to [0 0 0]

and

If the block parameter Rotational rate of the Custom planet is 0

Then

The initial output of w_rel (rad/sec) and w_b (rad/sec) should be zero.

I see exactly that result (using R2024a), so I can't recreate what you're reporting based on the information at hand.

The calculation of 'pdot,qdot,rdot' inside the block, which is is output dw_b/dt, should not influence the values of w_rel (rad/sec) and w_b (rad/sec) at t = 0.

P.S. I think it's just great that the notation used inside the block implementation is different than the notation used on the doc page.

[Thomas]

Hi Paul,

How were you able to run the block with a zero initial velocity vector? I am getting a divide-by-zero error at the integrator within 'Calculate Velocity in Body Axes' when I attempt this.

[Paul]

Here is the model:

Here are the first few frames of output from the this model:

>> [out.simout.Time(1:3), out.simout.Data(1:3,:)]

ans =

                         0                         0                         0                         0
     2.009509145207664e-06     2.009508830145579e-04     2.009509462393110e-04     2.009509145207664e-04
     1.205705487124598e-05     1.205705298087348e-03     1.205705677435867e-03     1.205705487124596e-03

>> [out.simout1.Time(1:3), out.simout1.Data(1:3,:)]

ans =

                         0                         0                         0                         0
     2.009509145207664e-06     2.009509145207664e-04     2.009509145207664e-04     2.009509145207664e-04
     1.205705487124598e-05     1.205705487124599e-03     1.205705487124599e-03     1.205705487124599e-03

I used the default settings for everything in the 6DOF ECEF (Quaternion) block, except I set the

Rotational rate [rad/sec]

to 0 for the Custom planet.

Everything looks as it should, at least at t = 0 (I didn't try to independently verify the results for t > 0.

"I am getting a divide-by-zero error at the integrator within 'Calculate Velocity in Body Axes' "

I don't see how that can happen because that subsystem doesn't have a divide operation, at least as far as I can tell (in R2024a).

Please verify that we are talking about the same block, which I've now linked to twice in this thread.

[Thomas]

Hi Paul,

I am referring to the block you linked above in all cases. I correct my previous statement regarding the divide-by-zero error; the error is that a derivative that is not finite is fed into the velocity integrator when I specify zero initial velocity.

[Paul]

The original issue (I think) in this question was that w_b was nonzero at t = 0, even though the block parameters were set such that w_e and w_ned were both zero at t = 0. When you posed the question originally, how were you getting w_ned to be zero at t = 0 without getting the error that you're seeing now?

I guess one way you'd see the error that you're seeing now would be if the "Initial mass" block parameter was set to 0.

In any event, I've shown that the block (R2024a) can be used with no errors and with the desired initialization at t = 0 using all of the default parameters except the Custom planet "Rotation rate (rad/sec)", so either you're using a different version of the block, or your block parameters aren't set to give the expected result (or both, I suppose).

[Paul]

Hi Thomas,

Here is a a simple model that demonsrates that w_b and w_rel are zero for all time if the input forces and moments are zero with a custom planet with W_e = 0.

bdclose('all');

hmdl = new_system('ECEF');

hf = add_block('simulink/Sources/Constant','ECEF/Forces');
set_param(hf,'Value','[0 0 0]');

hm = add_block('simulink/Sources/Constant','ECEF/Moments');
set_param(hm,'Value','[0 0 0]');

h6dof = add_block('aerolib6dof2/6DOF ECEF (Quaternion)','ECEF/6DOF');
set_param(h6dof,'ptype','Custom');
set_param(h6dof,'W_e','0');

houtwrel = add_block('simulink/Sinks/To Workspace','ECEF/w_rel');
set_param(houtwrel,'VariableName','w_rel');

houtwb = add_block('simulink/Sinks/To Workspace','ECEF/w_b');
set_param(houtwb,'VariableName','w_b');

add_line('ECEF','Forces/1','6DOF/1');
add_line('ECEF','Moments/1','6DOF/2');
add_line('ECEF','6DOF/9','w_rel/1');
add_line('ECEF','6DOF/10','w_b/1');

% open on canvas
% open_system('ECEF');

out = sim('ECEF');
figure
plot(out.w_b);

figure
plot(out.w_rel);

[Thomas]

Hi Paul,

A final point regarding this which I think resolves the issue.

I took another look through the documentation and realized that the initial body rates specified are in the NED frame (NOT in the inertial frame). If there is no initial velocity, then the rotation rates of the local horizontal frame will be zero (hence w_b = w_rel); however, if an initial velocity is present, the local horizontal frame will have a nonzero rotation rate as the body will be changing its latitude and/or longitude. Therefore, for a nonzero initial velocity, if zero INERTIAL body rates are desired, the value input for the initial body rates in the ECEF block must be [0 0 0] - DCM_be * w_ned (if w_e is zero). When I make this correction, the body rates with respect to the inertial frame are consistently zero as expected.

I think this is isn't the best way to input initial conditions (a better way would be to have the initial rates w.r.t. the inertial frame), but this is how the block is for now so we have to live with it.

[Paul]

Correct, as stated previously, the fundamental equaion is:

omega_b = omega_e + omega_ned + omega_rel

The block parameter

Initial body rotation rates [p,q,r]

specify the initial conditions on omega_rel, which is the angular velocity of the body relative to the instantaneous NED frame.

The initial condition on omega_ ned, the angular velocity of the NED frame relative to the earth-fixed frame, is determined from "the initial velocity and position of the vehicle."

So if we have a non-zero initial velocity s.t. omega_ned is non-zero, then the intial value of omega_rel would need to cancel that out in order for omega_b (angular velocity relative to inertial space) to be zero (assuming omega_e is zero as previously stated).

BTW, there is yet another error on the doc page. The equation in the Algorithms section for omega_ned is incorrect. It shows omega_ned as a function of N and M, which are defined directly underneath that equation as components of the applied moments. Of course, that equation for omega_ned is incorrect and the M and N should be replaced with the appropriate ellipsoidal earth parameters. Between this thread and your other thread, I'm finding that doc page to be very impressive (heavy sarcasm).