BUG: u_dot_generalized applies the center of mass lever arm with the wrong sign

[thc1006]

Describe the bug

Flight.u_dot_generalized reads r_CM and r_NOZ with the opposite orientation to the one the formulas around them assume, so the angular solve ends up with the wrong lever arm. A lateral force applied at the center of dry mass pitches the rocket the other way from what Newton-Euler gives.

Rocket.evaluate_com_to_cdm_function and Rocket.evaluate_nozzle_to_cdm both carry a leading minus, so they run CM to CDM and NOZ to CDM. The Equations of Motion v1 page defines r_CM the other way. Its linear equation is

m (v' + w' x r_CM + w x (w x r_CM) + r_CM'' + 2 w x r_CM') = ...

and that left hand side is m a_CM only when r_CM runs CDM to CM. Rocket.__evaluate_single_surface_cp_to_cdm has no leading minus, so surfaces_cp_to_cdm already matches the documented orientation and the aerodynamic moments come out about the CDM as the derivation expects.

The documented block system is the usual rigid body pair for a body fixed reference point, and the two coupling matrices are transposes of one another:

M v' + [m r]x^T w' = T20
I w' + [m r]x   v' = T21

In the code the same one shows up in both:

w_dot = I_CM.inverse @ (T21 + (T20 ^ r_CM))
v_dot = K @ (T20 / total_mass - (r_CM ^ w_dot)) - 2 * (w_earth ^ v)

Taking r_CM as the property actually returns it, the v_dot line is the consistent one and w_dot is not. I mention that because flipping the sign in v_dot is the change the algebra suggests first, and on its own it makes things worse.

After the cancellations the difference is

w_dot_code - w_dot_true = 2 * I_CM^-1 * (r_CM x (thrust + aerodynamic force))

Thrust is always axial here, so r_CM x thrust vanishes and the aerodynamic part is what survives. The lever from the center of mass to a center of pressure comes out as cp + r_CM rather than cp - r_CM. On calisto at ignition r_CM is 0.132 m against a 0.261 m arm, so the restoring moment acts through roughly twice the distance it should.

To Reproduce

Constant mass, no gravity, no drag, w = 0, and one lateral force applied at the center of dry mass with no couple. Everything except the term in question is exactly zero, so the expected answer is closed form.

import numpy as np
from rocketpy.mathutils.vector_matrix import Matrix, Vector
from rocketpy.simulation.flight import Flight

MASS, COM_TO_CDM = 40.0, 0.30
INERTIA = [[60.0, 0.0, 0.0], [0.0, 60.0, 0.0], [0.0, 0.0, 2.0]]
FORCE = np.array([100.0, 0.0, 0.0])
ZERO = Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])


class Const:  # every Function u_dot_generalized touches, held constant
    def __init__(self, v):
        self.v = v

    def get_value_opt(self, t):
        return self.v

    def differentiate_complex_step(self, t):
        return 0.0

    def differentiate(self, t, order=1):
        return 0.0


class Surface:  # returns a prescribed force and moment, so aero is exact
    reference_length = 1.0

    def compute_forces_and_moments(self, *a, **k):
        return (*FORCE, 0.0, 0.0, 0.0)


class Motor:
    burn_start_time = burn_out_time = 0.0
    nozzle_radius, thrust = 0.05, Const(0.0)

    def pressure_thrust(self, p):
        return 0.0


class Rocket:
    area, radius = 1.0, 0.5
    cp_eccentricity_x = cp_eccentricity_y = 0.0
    thrust_eccentricity_x = thrust_eccentricity_y = 0.0
    air_brakes, motor = [], Motor()
    total_mass, total_mass_flow_rate = Const(MASS), Const(0.0)
    com_to_cdm_function, nozzle_to_cdm = Const(COM_TO_CDM), 1.25
    nozzle_gyration_tensor = ZERO
    aerodynamic_surfaces = [(Surface(), None)]
    surfaces_cp_to_cdm = {aerodynamic_surfaces[0][0]: Vector([0, 0, 0])}

    def get_inertia_tensor_at_time(self, t):
        return Matrix(INERTIA)

    def get_inertia_tensor_derivative_at_time(self, t):
        return ZERO

    def power_off_drag_7d(self, *a):
        return 0.0

    power_on_drag_7d = power_off_drag_7d


class Env:
    earth_rotation_vector = [0.0, 0.0, 0.0]
    density = gravity = pressure = Const(0.0)
    wind_velocity_x = wind_velocity_y = Const(0.0)
    speed_of_sound, dynamic_viscosity = Const(340.0), Const(1.8e-5)


flight = Flight.__new__(Flight)
flight.rocket, flight.env = Rocket(), Env()
u = [0, 0, 1000, 0, 0, 10, 1, 0, 0, 0, 0, 0, 0]
omega_dot = np.array(flight.u_dot_generalized(0.0, u)[10:13])

# constant mass, no gravity, no rotation, force applied at the CDM: the moment
# about the CM is (CM -> CDM) x F, so this is a closed-form Newton-Euler answer
lever = np.array([0.0, 0.0, COM_TO_CDM])
inertia_cm = np.array(INERTIA) - MASS * (
    lever @ lever * np.eye(3) - np.outer(lever, lever)
)
expected = np.linalg.solve(inertia_cm, np.cross(lever, FORCE))

print(f"u_dot_generalized : {omega_dot}")
print(f"Newton-Euler      : {expected}")

On develop at 7e785a6:

u_dot_generalized : [ 0.         -0.53191489  0.        ]
Newton-Euler      : [ 0.          0.53191489  0.        ]

Expected behavior

The two lines agree. More generally, the same rocket described about its center of dry mass and about its center of mass should give the same angular acceleration, which is the check I found most useful because it does not depend on which orientation you consider authoritative.

Additional context

Apogee barely moves, which is probably why this has gone unnoticed. The lateral state does. On flight_calisto_custom_wind, at apogee, velocity x goes from -14.826366 to -11.634095 m/s and the aerodynamic moment M1 from -0.509418 to -0.652606 N·m, while apogee above ground goes from 3423.2700 to 3427.7198 m.

What I have locally, if it is useful: negate both vectors where they are read, keep w_dot as it is, and flip the coupling sign in v_dot. Four changes that have to move together. It comes with eleven tests in two layers, four of which do not depend on the documented orientation at all, plus a control test that survives every one of the four reverts while at least one named test goes red for each. test_rocket_csys_equivalence still passes, so nose_to_tail is covered.

I would rather ask before sending a patch, since it moves six recorded values and changes everyone's lateral trajectories. Happy to open the PR, or to leave it with you if you would rather take it a different way. This code came in with #364 and the derivation there is sound as far as I can tell; what I think went wrong is only the orientation the two properties hand over.

One thing to watch if this does get changed: it interacts with #1185. For a rocket whose com_to_cdm_function ends up with a callable source, the - total_mass * r_CM_ddot term carries a very large spurious value at ignition, and this change flips its sign. If a flight is integrated through the ignition instant rather than sitting on the rail, that turns a large upward impulse into a downward one.

[thc1006]

/cc @zuorenchen

[thc1006]

A data point I picked up while working on #1188, since it bears on this one and I would not have gone looking for it deliberately.

tests/acceptance/test_defiance_rocket.py carries a measured number, MEASURED_APOGEE_AGL = 9308.32, from the 2024 Defiance flight. Running that example on three trees:

	apogee AGL	error vs measured
develop	9242.3701 m	0.709 %
nozzle gyration tensor fixed (#1188)	9243.8350 m	0.693 %
both that and the sign convention here	9247.9733 m	0.648 %

Each correction moves the simulation toward the measured apogee rather than away from it. It is one flight and the margin is small, so I would not lean on it on its own, but it is the only real-flight comparison I have found that either change touches, and it does not argue against them.

The impact point moves a good deal more, from (1625.57, 81.79) on develop to (1550.23, 104.80) with both changes in. There is no measured impact point in that file to check it against.