Implement a tire model and compute the normal force

[tjstienstra]

This PR aims to implement a basic setup for creating more advanced tire models. The most difficult issue to solve here is to reliably compute the normal force by introducing auxiliary speeds and adding noncontributing loads to the system.

The main problem of these auxiliary speeds is to ensure correct propagation throughout the velocity graph and to take them into account in the velocity constraints. Related is the introduction of the velocity_constraints attribute in sympy#26374.

I have designed a new utility called AuxiliaryDataHandler to tackle this issue. The basic idea is to keep all auxiliary speeds separate and introduce them at the end of the define kinematics step. The position graph of points and the orientation graph of reference frames do not contain any cycles. Therefore, you can convert them into a tree representation with the inertial frame as the root of the orientation graph and a fixed point in the inertial frame as the root of the position graph. The advantage of this tree representation is that if you know that an auxiliary speed should be applied to a certain point, then you can easily determine to what points this speed should propagate (namely its children). The auxiliary speed of each point can is also remembered such that they can be requested when determining the velocity constraints. The load objects of the noncontributing loads are defined as disconnected points with the auxiliary velocity in the inertial frame. It is known that they do no work, so all other terms (if there are any) must fall out anyway.

The main disadvantage in the current setup is that this auxiliary data handler must be shared among all BRiM objects. This does break a bit with the idea that every model is its own system, but I do not really see another way as auxiliary speeds from other systems will have to be able to propagate throughout everything. Some consequences of this decision involve the introduction of the ModelBase.is_root attribute and that the noncontributing loads/auxiliary speeds are only added to the root system.

Another fun feature of the AuxiliaryDataHandler is its _compute_velocity method, which tries to compute the velocity of a point based on v2pt_theory or v1pt_theory before falling back to the time derivative of the position vector. I am not entirely sure if I'll keep this feature or will remove its usage. This should be decided based on some benchmarks.

TODO:

• [x] Correctly compute the normal force of both the rear and front tire.
• [x] Implement some basic testing.
• [x] Compare benchmarks of different methods to compute the velocity of points.
• [x] Fix failing tests.
• [x] Design a generic tire model interface.
• [x] Ensure full line coverage with the tests.
• [x] Add an explanation section to the AuxiliaryDataHandler.
• [x] Add a disclaimer that one should aim to compute velocities based on other points according to the tree structure of the position tree. It is kind of something you never encounter to violate anyways, but it is one of the requirements of the auxiliary data handler to function correctly.
• [x] Validate tire model.

Future features could be:

• Support for noncontributing torques.
• More advanced and specific tire models.
• Make a tutorial on how to work with noncontributing forces. Users will have to know how to solve and use the auxiliary equations after forming the equations of motion with Kane's method.

[moorepants]

Will this implement the functionality to be able to make a tire model as well as the tire model? For the later, there are some considerations to make, like what conventions to use, whether it is a linear or nonlinear model, what model it is, etc.

For a full nonlinear tire model the forces all depend on the normal force, but if that is a noncontributing force for the basic formulation we have. In very general tire force models, the lateral forces and various moments can be nonlinear in the normal force which means you can't simply solve for the noncontributing force algebraically, you'd need to have an iterative numeric solver to handle that. I think that some tire models implement the vertical tire compliance to get around this issue. Then the normal force is a function of the new state variables and there is no need to try to solve for this normal force.

[tjstienstra]

Will this implement the functionality to be able to make a tire model as well as the tire model?

The main aim is the functionality to be able to make a tire model. Implementing a tire model is more like a second.

I specifically had in mind to first implement a more generic tire model where the tire forces and torques are just dynamicsymbols. I would expect it to result in something like the following:

class InContactTireTireBase):
    """Generic tire model that is in contact with the ground."""

    def __init__(self, name: str) -> None:
        super().__init__(name)
        self.compute_normal_force: bool = False
        self.no_longitudinal_slip: bool = False
        self.no_lateral_slip: bool = False

    def _define_objects(self) -> None:
        """Define the objects of the tire model."""
        super()._define_objects()
        self.symbols.update({
            name: dynamicsymbols(self._add_prefix(name))
            for name in ("Fx", "Fy", "Fz", "Mx", "Mz")
        })
        if self.compute_normal_force:
            self.u = Matrix([dynamicsymbols(self._add_prefix("uaux_z"))])
        if self.no_longitudinal_slip:
            self.symbols["Fx"] = S.Zero
        if self.no_lateral_slip:
            self.symbols["Fy"] = S.Zero
            self.symbols["Mz"] = S.Zero
        if self.no_lateral_slip and self.no_longitudinal_slip:
            for name in ("Fx", "Fy", "Mx", "Mz"):
                self.symbols[name] = S.Zero

    def _define_kinematics(self) -> None:
        super()._define_kinematics()
        self._set_pos_contact_point()
        super()._define_kinematics()
        self._set_pos_contact_point()
        if self.compute_normal_force:
            direction = -self.ground.get_normal(self.contact_point)
            if self.on_ground:
                direction = -direction
            self.auxiliary_handler.add_noncontributing_force(
                self.contact_point, direction, self.u[0], self.symbols["Fz"])

    def _define_loads(self) -> None:
        super()._define_loads()
        # TODO add the loads based on the tire model.

    def _define_constraints(self) -> None:
        super()._define_constraints()
        if self.no_longitudinal_slip and self.no_lateral_slip:
            # TODO add both nonholonomic constraints for no slip.
            pass
        elif self.no_longitudinal_slip:
            # TODO compute nonholonomic constraint for no longitudinal slip.
            pass
        elif self.no_lateral_slip:
            # TODO compute nonholonomic constraint for no lateral slip.
            pass
        if not self.on_ground:
            # TODO add the holonomic constraints.
            pass

The nice thing about the above setup is that this class can be inherited and one can simply specify self.symbols["Fx] = some_function(self.symbols["Fz"]). Here one can choose for more details, the only constrain is that no compliance is modeled.

About implementing a linear or nonlinear model. I don't expect to implement full nonlinear models right now as that would result in a lot of complexity as you mention. A model I will probably start with is the same one used in the kickplate model.

[moorepants]

A model I will probably start with is the same one used in the kickplate model.

Ok.

[tjstienstra]

Benchmark results of the Whipple bicycle following Moore's convention:		Time	#operation before CSE	#operations after CSE
main 2ed1713f737d4ec1dc1203cb3d558a27f5d81552	3.04	352210	2291
this PR with custom velocity computation da346de773e073ec912afeab853b16db68c661cf	3.49	368849	2762
this PR with Point.vel	2.85	266654	2526
this PR with custom velocity computation and ICR assumption front wheel	3.34	388994	2484
this PR with Point.vel and ICR assumption front wheel	3.10	352210	2291

ICR assumption front wheel refers to whether the front wheel center's velocity is computed based on the assumption that the contact point is an instantaneous center of rotation. The problem with wanting to apply this assumption is that the velocity tree doesn't match the shape of the position tree resulting in the addition of the wrong auxiliary velocities. I will start thinking about how to resolve this problem, but it will require more active tracking of the graphs and that users/modelers notify the auxiliary data handler when they apply these assumptions.

At least for now, I can conclude that it is best to use Point.vel as taking time derivatives of position vectors seems to give really good results.

[tjstienstra]

The benchmark results when using v2pt_theory above just differentiating the vector have always surprised me, as you would normally expect that it is cheaper to use such a neat rule instead of expressing a vector into a frame and then differentiating it. Therefore, I had a bit of a deeper dive into this issue. The problem seems to be a bad choice of how we apply our cross-product. In v2pt_theory we apply our cross-product using the unit vectors of the frames used in the position vector, while it would be way more efficient to use the unit vectors of the angular velocity. I just experimented with this a little by changing self.set_vel(outframe, v + (omega ^ dist)) in v2pt_theory to self.set_vel(outframe, v - (dist ^ omega)). This change resulted in 2.59 seconds, 232020 operations before CSE and 1919 operations after CSE.

[moorepants]

I think that Peter may have come across this before, also showing that equations were longer when using a method that you would think would be shorter.

[tjstienstra]

I think that Peter may have come across this before, also showing that equations were longer when using a method that you would think would be shorter.

Yeah, wouldn't surprise me.

I may open an issue in SymPy on whether we shall make this change. I just have to think a bit more carefully about how to explain/prove that this is an improvement, while also hoping that it won't break a lot of tests.

[tjstienstra]

It just removed the yaw frame in the formulation of the Whipple bicycle following Moore's convention. This causes a reduction in the number of operations in the EoMs after CSE from almost ~2500 to ~1900 (the number of operations before CSE only slightly increased from 230224 to 232020). This change seems a bit odd to me, as orient_body_fixed just does two rotations and stores them as one. It would be nice to investigate what calculation exactly profits from this change. If I have to guess then I would start with investigating vector differentiation.

[moorepants]

I think there is an open issues somewhere with this exact finding from Peter.

[tjstienstra]

I think there is an open issues somewhere with this exact finding from Peter.

Found it sympy#23075. The issue is a bit long, so I quickly skimmed through it. The advantage gained by orient_body_fixed is that it safes the angular velocity expressed in the child frame. The following two give the same number of operations:

yaw_frame = ReferenceFrame("yaw_frame")
roll_frame = ReferenceFrame("roll_frame")
yaw_frame.orient_axis(self.ground.frame, self.ground.frame.z, self.q[2])
roll_frame.orient_axis(yaw_frame, yaw_frame.x, self.q[3])
roll_frame.set_ang_vel(self.ground.frame, roll_frame.ang_vel_in(self.ground.frame).express(roll_frame))

roll_frame = ReferenceFrame("roll_frame")
roll_frame.orient_body_fixed(self.ground.frame, (*self.q[2:4], 0), "zxy")

[tjstienstra]

An optional feature to be implemented in InContactTire is a property and accompanying attribute to define tire model-specific equations. The implementation would be something as follows:

class InContactTire(TireBase):
    ...
    def __init__(self, name: str):
        ...
        self.substitute_load_eqs: bool = True

    def _define_loads(self) -> None:
        ...
        def get_symbol(load_str: str) -> Expr:
            sym = self.symbols.get(load_str, 0)
            if self.substitute_load_eqs:
                return self.load_equations.get(sym, sym)
            return sym
        ...
        tire_force = (
            get_symbol("Fx") * self.longitudinal_axis +
            get_symbol("Fy") * self.lateral_axis
        )
        ...

    @property
    def load_equations(self) -> dict[Function, Expr]:
        """..."""
        return {}

Two reasons why I would propose the above structure:

• The equations will probably depend on the slip angle and/or chamber angle. These can only be computed after the define kinematics step.
• You don't always want to make this substitution. It can be favorable to leave tire forces specified as functions of time.