Motion Equation Derivation for Constrained Variable-Mass Systems Using the Extended Gibbs-Appell Formulation

Modeling the dynamics of constrained mechanical systems is fundamental to advancements in robotics, aerospace, and automation. These systems, whether robotic manipulators or wheeled vehicles, often operate under constraints that restrict their motion, such as joint limits or no-slip conditions. Variable-mass systems, defined as mechanical systems where mass changes over time due to material expulsion (e.g., fuel in rockets) or depletion (e.g., liquid in agricultural sprayers), pose significant challenges, as seen in spacecraft expending fuel or agricultural robots depleting liquid tanks [1]. Such scenarios demand sophisticated methods to derive accurate motion equations, particularly when constraints and variable mass intertwine. Traditional approaches, like Lagrange’s equations, falter in these complex cases, often yielding unwieldy expressions or struggling with nonholonomic constraints. The G-A formulation offers a promising alternative, leveraging quasi-velocities to streamline derivations, yet its standard form falls short for systems with time-varying mass [2], [3], [4].

Variable-mass systems, characterized by time-varying mass, are ubiquitous in engineering dynamics. Canonical examples include rockets, where mass decreases due to fuel combustion, and robotic systems handling dynamic payloads, such as agricultural sprayers with depleting tanks [5]. Recent research has expanded the study of variable-mass systems to complex scenarios [6], such as spacecraft with fuel sloshing and manipulators with variable payloads. These works highlight the challenges of time-varying inertia, which introduces nonlinear terms into the motion equations [7]. However, the literature predominantly focuses on unconstrained or holonomically constrained systems, with limited attention to variable-mass systems under nonholonomic constraints. In precision agriculture, mobile robots like sprayers exemplify this complexity, as their tank mass decreases during operation, altering dynamics critical for navigation and spraying accuracy [8], [9]. The scarcity of methods addressing variable mass alongside diverse constraints motivates this current study, as such systems are increasingly prevalent in robotics and automation.

Constraints are the backbone of mechanical system dynamics, dictating permissible motions. Holonomic constraints, such as those imposed by robotic arm joints, are integrable, expressible as algebraic relations, and reduce the system’s degrees of freedom [10]. Nonholonomic constraints, prevalent in wheeled mobile robots, involve non-integrable velocity relations, preserving degrees of freedom but complicating dynamics due to their differential nature [11]. In robotics, constraints are ubiquitous. Mobile robots, like agricultural sprayers, navigate under nonholonomic constraints from wheeled motion, while holonomic constraints arise from mechanical components like articulated frames [12], [13]. Modeling these systems requires handling constraint forces, often through Lagrange multipliers in traditional methods. However, multipliers increase computational complexity, particularly for nonholonomic systems, where constraints resist integration. Recent advances, such as Maggi’s equations [14] or the Udwadia-Kalaba formulation [15], offer alternatives but assume constant mass, limiting their applicability to variable-mass systems. The interplay of constraints and variable mass—crucial for robots with dynamic payloads—remains a largely uncharted frontier, necessitating innovative approaches.

The G-A formulation offers a robust framework for deriving motion equations for constrained systems, leveraging quasi-velocities to simplify dynamics [16]. Unlike Lagrange’s method, which relies on generalized coordinates and multipliers, G-A employs quasi-velocities—linear combinations of velocities tailored to the system’s constraints. The G-A formulation has gained traction in multibody dynamics, with applications in manipulators, vehicles, and spacecraft [17], [18], [19]. Its efficiency stems from reducing the number of equations to the system’s degrees of freedom, making it ideal for real-time control. However, the standard G-A formulation assumes constant mass, rendering it inadequate for systems like agricultural sprayers, where mass variation alters the inertia matrix. Recent studies have applied G-A to nonholonomic systems but have rarely addressed variable mass, highlighting a critical gap this current study seeks to fill.

The integration of variable mass with constraints introduces time-dependent inertia, complicating the mass matrix and adding non-linear terms [20], [21]. For example, in a mobile robot sprayer, the tank’s mass depletion shifts the center of mass, affecting acceleration and constraint forces [22], [23]. Nonholonomic constraints, such as no-slip conditions for wheels, impose differential equations that resist integration, while holonomic constraints, like joint limits, add algebraic complexity. Traditional methods like Lagrange’s equations yield lengthy expressions with multipliers, and Newton-Euler approaches require extensive numerical computations, hindering real-time applications. In robotics, this complexity is particularly pronounced. Mobile robots operate under stringent constraints—wheels impose nonholonomic limits, and mechanical components add holonomic ones [24]. A sprayer robot [25], tasked with uniform pesticide application, must maintain precise navigation despite a depleting tank. Existing methods often simplify by assuming constant mass or neglecting constraints, leading to inaccurate models that compromise performance. The need for a unified method to handle variable mass and constraints is evident, especially for precision agriculture, where efficiency and accuracy are paramount.

To address this gap, an extended G-A formulation was proposed for variable-mass systems in this study, accommodating both holonomic and nonholonomic constraints. Its efficacy was demonstrated through application to a mobile robot sprayer in precision agriculture, a domain where precise dynamics modeling is paramount for efficient navigation and resource application. This study focuses on a mobile robot sprayer in precision agriculture, a field revolutionizing farming through automation. These robots navigate fields, dispersing fertilizers or pesticides from tanks whose mass decreases over time. Nonholonomic constraints from wheels ensure stable motion, while holonomic constraints from mechanical frames maintain structural integrity. The variable-mass tank complicates dynamics, as its depletion alters inertia, impacting path planning and spraying precision. Standard G-A cannot address this, as it lacks provisions for mass variation. The proposed extended G-A formulation overcomes this limitation, adapting quasi-velocities to accommodate time-varying mass and diverse constraints, delivering a precise model for such robots. This contribution advances theoretical dynamics and offers practical tools for robotics and automation. This study is structured below. Section 2 details the extended G-A formulation. Section 3 models the reactive force from fluid expulsion for incorporation into the extended G-A formulation. Section 4 derives the motion equations and applies the method to the sprayer robot. Section 5 validates the approach through numerical simulations for different paths. Section 6 concludes with implications and future directions.

Consider a mechanical system of $N$ particles, each with potentially variable mass. These particles move within an inertial coordinate system. The configuration of each particle was described using $n$ independent generalized coordinates, denoted by $q \equiv\left\{q_1, q_2, \ldots, q_n\right\}$. The position vector of the $i$-th particle, denoted by $r_i$, can be expressed as a function of these generalized coordinates and time, i.e., $r_i=r_i\left(q_i t\right)$. The velocity of the $i$-th particle, $v_i=\dot{r}_j$, is then given by:

To obtain the acceleration $a_i=\dot{V}_i$, the velocity evaluated in Eq.(1) with respect to time was differentiated:

Consider again $N$ particles, which are subject to known forces ($G$), restraining counter-forces ($H$), and counter-thrust forces $(R)$. The differential motion equation for a particle is:

where, $m_i$ denotes the mass of the $i$-th particle. If the system's constraints are ideal, they perform no work under virtual displacements. Mathematically, this means that for any virtual displacement $\delta r_i$ :

As a result, when the dot product of Eq. (3) with the virtual displacement was performed and then summed over all particles $i$, the reaction forces $H$ were eliminated. Then, it concluded:

Eq. (5) represents a new form of the motion equation for mechanical variable-mass systems, applicable to both holonomic and nonholonomic systems. If mass is made constant, then Eq. (5) becomes as simple as the fundamental equation of dynamics.

The constraints imposed on the mechanical system's motion are represented by the general equation: $f_\beta(q, \dot{q}, \ddot{q}, t)=0$, where $f$ is a function of the generalized coordinates $q$, with their time derivatives $\dot{q}$ and $\ddot{q}$ and time $t$. The nature of the constraint is then determined by the properties of the function $f$. These constraints can be classified as follows:

If $f$ is independent of $\dot{q}$, the constraint is holonomic.

If $f$ depends on $\dot{q}$, and cannot be integrated to remove the dependence of the velocity, the constraint is nonholonomic.

If $f$ is linear in q and $\dot{q}$, the constraint is linear; otherwise, it is nonlinear.

The order of the constraints is determined by the highest derivative of $q$ present in $f$.

Having established the constraints governing the mechanical system, an appropriate set of variables needs to be defined to facilitate the application of the G-A formulation. This process requires a careful selection of generalized coordinates, quasi-velocities, and, where applicable, quasi-coordinates, tailored to the specific nature of the constraints. The defined variables must effectively incorporate the constraints into the G-A formulation. For nonholonomic constraints, this often entails expressing the constraints in terms of the quasi-velocities, leading to simplified constraint equations. For holonomic constraints, the constraints may be used to reduce the number of independent generalized coordinates, thereby simplifying the system.

Furthermore, the representation of accelerations $\ddot{q}$ in terms of the chosen variables is necessary for the application of the G-A formulation. This may involve deriving relationships between the time derivatives of the quasi-velocities and the generalized accelerations. The objective is to obtain an expression for the G-A function in terms of the selected variables, enabling the derivation of the motion equations via Eq. (10).

where, $Q$ represents the generalized non-conservative forces, including any contributions from the constraint forces. In comparison to Eq. (5), it’s equal to the right-hand side expressions.

Finally, it is imperative to provide a clear justification for the selection of variables, emphasizing their physical relevance and their ability to streamline the analysis. The chosen variables should ideally reflect the inherent symmetries and physical characteristics of the constrained system, facilitating a more transparent and insightful understanding of its dynamics. By systematically defining these variables, a robust foundation was established for the derivation of the G-A equations, enabling a comprehensive analysis of the constrained mechanical system.

Now consider again Eq. (2); the next goal is to express $\ddot{q}_s$ in terms of the chosen variables (e.g., generalized coordinates, quasi-velocities, etc.), considering the constraints as follows:

As for the holonomic constraint as per Eq. (6), the following equation is given:

From this equation, $\ddot{q}_s$ can be solved in terms of the other accelerations and velocities. This allows us to eliminate a redundant acceleration from the G-A function.

As for the nonholonomic constraint as per Eq. (7), it is beneficial to express the constraint in terms of quasi-velocities as follows:

The resulting equation in terms of quasi-velocities and the time derivative of the quasi-velocities is as follows:

For a higher-order constraint, the constraint equation with respect to time was differentiated as many times as necessary to obtain an expression that relates the highest-order derivative to lower-order derivatives. Next, by substituting the expression into the acceleration equations for the particles in the system, the higher-order constraints were introduced into the acceleration expressions.

Now, it’s the turn to importing the details of constraints in G-A formulation extraction. In terms of generalized coordinates $q$, virtual displacements are expressed as:

\text { where, } $\delta q_s$ \text { are the virtual displacements of the generalized coordinates. }

Again, consider that the expression on the right-hand side of Eq. (5) is the virtual work performed by all external known forces ($G$) and counter thrusts $(R)$ under the virtual displacements $\delta r_i$, it may be expressed as:

where, it is considered that

where, $Q_s$ represents general forces corresponding to the generalized coordinate $q_s$, and $\psi_s$ qrepresents the term related to the mass variation. Eqs. (16)-(17) can be replaced with the following:

where, $\frac{\partial \mathrm{r}_i}{\partial q_s}=\frac{\partial \ddot{\mathrm{r}}_i}{\partial \ddot{q}_s}$ . Eq. (14) was recalled to reformulate the fundamental equation of dynamic as follows:

Given the presence of different constraints (holonomic, nonholonomic, and higher-order), the method for evaluating virtual displacements must be adjusted accordingly. Each constraint equation effectively reduces the number of independent virtual displacements, which in turn limits the number of independent quasi-velocities or generalized coordinates.

Following the preceding considerations, the quantity ($G$) is defined as the Gibbs function, also known as the energy of acceleration, for a mechanical system subject to constraints and potentially exhibiting mass variation.

By differentiating with respect to quasi accelerations, the following equation can be derived:

where,

where, $g$ represents the number of restricted quasi-velocities resulting from the known constraints. Let's consider that $G$ is a function of the form $G=G(q, \dot{q}, t)$ and the fundamental equation of dynamic, i.e., Eq. (5). The right-hand side can be derived by considering Eq. (15) as follows:

Hence, it can be concluded that

Since the variations $\delta q_s$ are independent, the resulting equation is equal to the following:

These are the G-A motion equations for the constrained variable-mass systems.

When the robot sprayer ejects water (or pesticide/fertilizer), it experiences a reactive force due to the momentum of the expelled fluid, governed by Newton’s third law and variable-mass dynamics. This force opposes the robot’s motion or alters its trajectory, depending on the spray direction and velocity. Additionally, the spraying process may introduce the following resistive forces:

• Drag or viscous resistance: The spray’s interaction with air or crops may create a drag force opposing the robot’s motion.

• Frictional effects: The mechanical spraying mechanism (e.g., nozzles and pumps) may introduce internal resistance.

• Mass depletion: The tank’s mass decreases, reducing the robot’s inertia and affecting acceleration, which must be modeled dynamically.

These effects are particularly significant for a mobile robot with nonholonomic constraints (e.g., wheeled motion) and a variable-mass tank, as they influence navigation and spraying precision in agricultural fields.

The robot’s tank loses mass as water is sprayed, making it a variable-mass system. The reactive force from spraying is derived from the Tsiolkovsky equation adapted for a multibody system. Assuming that the water is ejected at a relative velocity $u_i$ (relative to the robot) with a mass flow rate $\dot{m}_i=-\frac{d m_i}{d t}$ (negative since mass decreases), the following equations can be derived:

• Reactive force: The force is caused due to spraying as follows:

where, $u_i$ is the velocity vector of the ejected water (e.g., from nozzles, typically opposite to the spray direction).

• Mass variation: The robot’s mass, denoted as $m_i(t)$, decreases over time, e.g., linearly $\left(\dot{m}_i=\right.$ constant $)$ or exponentially, depending on the spraying rate.

In the G-A formulation, the variable mass affects the Gibbs function, as the mass matrix becomes time-dependent. The reactive force $F_{\text {react}_i}$ is included as an external generalized force.

• Drag force: If the spray interacts with air or crops, it may produce a drag force proportional to the spray velocity or flow rate. A simplified model is:

where, $k_d$ is the drag coefficient, and $\widehat{u}_i$ is the unit vector of the spray direction.

• Mechanical resistance: The spraying mechanism (e.g., pump friction) may introduce a resistive force, modeled as:

where, $k_f$ is the friction coefficient.

These forces were incorporated as generalized forces in the G-A formulation, acting alongside $F_{\text {react }}$. They contribute to Eq. (26) in term $Q$.

To derive the motion equations for a mobile robot sprayer equipped with multiple nozzles, the extended G-A formulation must account for variable-mass dynamics, reactive forces from multiple nozzles, nonholonomic constraints from wheeled motion, and resistive forces (e.g., drag and friction) introduced by spraying. The multiple nozzles introduce complexity, as each nozzle contributes a reactive force based on its spray direction, velocity, and mass flow rate, affecting the robot’s linear and angular motion.

The kinematics of the mobile platform with $N$ nozzles was investigated in this study. The robot's motion was analyzed with respect to a body-fixed coordinate frame ($X_0, Y_0, Z_0$) located at the center of mass, where the $x$-axis aligns with the forward direction, the $y$-axis is lateral, and the $z$-axis is vertical. In addition, the general coordinates were selected as $(X, Y, Z)$.

Figure 1 provides the view of the mobile spraying robot, annotated with the kinematic vectors to illustrate its motion characteristics. The primary kinematic vectors are defined as follows:

• Linear Linear velocity $\left(v_0\right)$ : The translational velocity of the robot's center of mass, directed along the $x$-axis, is denoted by $v_0$. This vector represents the robot's forward or backward motion.

• Angular velocity $\left(\omega_0\right)$ : The rotational velocity about the $Y_0$-axis, denoted by $\omega_0$, governs the robot's yaw motion (turning).

• Wheel angular velocities $\left(\dot{\theta}_{F R}, \dot{\theta}_{B R}, \dot{\theta}_{F L}\right.$, and $\left.\dot{\theta}_{B L}\right)$ : The four wheels, labeled $\theta_{F R}, \theta_{B R}, \theta_{F L}$, and $\theta_{B L}$, have individual angular velocities $\dot{\theta}_{F R}, \dot{\theta}_{B R}, \dot{\theta}_{F L}$, and $\dot{\theta}_{B L}$, respectively. These velocities constrain the direction of the wheel's plane due to the nonholonomic nature of the system, prohibiting lateral motion along the $y$-axis.

• Spray speed vectors $\left(u_{i R}\right.$ and $\left.u_{i L}\right)$ : The robot's base is equipped with $n$ nozzles on left- and righthand sides for spraying, each emitting a spray in a direction represented by vectors $u_{i R}$ and $u_{i L}$ (where $i=1,2, \ldots, n)$. These vectors depend on the nozzle orientation and speed of spraying, influencing the spraying coverage.

The nonholonomic constraints arise from the rolling-without-slipping condition of the wheels, which can be mathematically expressed as:

where, $\dot{z}_0$ and $\dot{x}_0$ are the velocity components of the robot's center of mass in the global frame, and $\theta_0$ is the robot's orientation angle. This constraint ensures that the robot's motion is restricted to forward/backward translation and rotation about the $Y_0$-axis.

Unlike a general mobile platform, the sprayer's operation often involves considerations such as maintaining a specific orientation of the spray boom relative to the ground or target, and potentially controlling the speed for consistent application. The generalized coordinate approach can be adapted, focusing on the essential variables for describing the sprayer's motion and the position of its spraying elements. But in this study, the sprayer speed was considered to be constant and the nozzle’s angle was defined in the beginning and remained fixed. Assuming the mobile robot sprayer has two parallel axles with a total of four wheels, the system's configuration can be defined using a set of generalized coordinates. Let the vector of the system’s generalized coordinates be the following:

where, $\boldsymbol{R}_b^T=\left(\begin{array}{ll}x_0 & z_0\end{array}\right)^T, \quad \boldsymbol{\Theta}_b^T=\left(\theta_0, \theta_{F R}, \theta_{B R}, \theta_{F L}, \theta_{B L}\right)^T$, and $\boldsymbol{\theta}_{\text {spray }}^T=\left(\psi_{i R} \quad \psi_{i L}\right)^T$ for $i=1, \ldots, n$. In addition, $\theta_{\text {spray }}^T$ represents the configuration variables specific to the spraying mechanism. This could include the following:

• The height of the spray booms relative to the base or ground.

• The angular orientation of the spray boom (e.g., pitch and yaw relative to the base if actively controlled).

• The configuration of individual spray nozzles (e.g., their opening angles and orientations if actively controlled).

Three nonholonomic constraints from rolling without slipping reduce the system to two independent generalized velocities, chosen as $\dot{\theta}_R=\frac{\theta_{F R}+\theta_{B R}}{2}$ and $\dot{\theta}_L=\frac{\theta_{F L}+\theta_{B L}}{2}$, representing the average angular velocities of the right and left wheels. These quasi-velocities simplify the G-A formulation by directly relating to motor inputs. These quasi-velocities can also be considered as $v_0$ and $\theta_0$, which represent the forward speed and the angular velocity of the mobile platform about its center, respectively. In this case, the configuration variables specific to the spraying mechanism are assumed to be fixed. Hence, the linear velocity of center points is:

In addition, the angular velocity is:

The velocity of point $A$ is:

Multiple nozzles at positions $r_{i L}$ and $r_{i R}$ emit spray with velocity vectors $u_{i R}$ and $u_{i L}$, contributing to reactive forces. The general emit sprays were evaluated with respect to general reference coordinates by considering the mobile platform speeds. In addition, the rates of spraying for both sides of the mobile platform are the same. The spray's exit speed from the $i$-th nozzle was evaluated as:

where, $\vartheta_{i R}$ represents the output speed from the $i$-th nozzle. Differentiating from Eqs. (36) and (37) can yield the accelerations:

\text { The vectors } $u_{i L} \text { and } \dot{u}_{i L} \text { were evaluated in the same manner as } u_{i R} \text { and } \dot{u}_{i R}$ \text {. }

The extended G-A formulation was applied to the mobile robot sprayer, modeled as a differential drive system with two degrees of freedom, described by quasi-velocities $v_0$ and $\theta_0$. The robot's total mass $m(t)=m_0-\sum_i \dot{m}_i t$ decreases due to liquid expulsion from $n$ nozzles, where $\dot{m}_i$ is the mass flow rate per nozzle. The Gibbs function represents the energy of acceleration. For simplicity, this study focuses on the mobile robot's point $A$ by considering the whole mass $m_0$ and using the acceleration vectors as per Eqs. (37) and (38). The Gibbs function for the mobile robot sprayer included the mobile platform $\left(G_b\right)$, robot wheels $\left(G_w\right)$ and the sprayer system $\left(G_s\right)$. The Gibbs function, which is the same as that in Eq. (21), for the mentioned system was evaluated as:

where,

To derive the motion equations, the partial derivatives of the Gibbs function with respect to the generalized quasi accelerations $\left(\dot{v}_0\right.$ and $\left.\ddot{\theta}_0\right)$ were computed.

The generalized forces $Q$ include motor torques $\tau_{F R}$ and $\tau_{B R}$ the right-side wheels and $\tau_{F L}$ and $\tau_{B L}$ for the left-side wheels, drag forces ( $F_{\text {drag }}$ ), and friction forces ( $F_{\text {fric }}$ ) as well as the reactive forces due to mass expulsion ( $F_{\text {react }}$ ) evaluated in Eqs. (27)-(29). These forces are projected onto the quasivelocities as follows:

where, $\dot{\Theta}=\left[\begin{array}{ll}v_0 & \dot{\theta}_0\end{array}\right]^{\mathrm{T}}$.

The remaining dynamics are terms in the G-A formulation not dependent on quasi-accelerations. They arise from velocity-squared terms and time-varying mass effects. These terms do not contribute directly. However, they appear in the expansion of acceleration vectors, affecting the left-hand side through cross-terms.

For $n$ nozzles, the total reactive force affects both linear and angular motions.

The motion equation of the platform is obtained as:

The nonholonomic constraints are implicitly satisfied through the quasi-velocities, eliminating the need for constraint forces in the equations.

To validate the motion equations derived using the extended G-A formulation for a mobile robot sprayer with variable mass ( Figure 2), numerical simulations were conducted in MATLAB. The simulations aimed to: (a) verify the dynamic consistency of the proposed model under variable-mass conditions and nonholonomic constraints, (b) import the PID control strategy for tracking predefined paths, and (c) assess the impact of reactive forces from spray ejection on navigation accuracy. Three representative paths—straight, circular, and zigzag—were simulated to emulate diverse operational scenarios in precision agriculture, such as linear spraying along crop rows, circular navigation around obstacles, and oscillatory motion for broad coverage. This section details the simulation setup, methodology, control implementation, and results, focusing on path-tracking performance, mass variation, and momentum balance.

Numerical simulations were conducted in MATLAB using the ODE45 solver, selected for its adaptive step-size Runge-Kutta method, ensuring high numerical accuracy with relative and absolute tolerances set to $10^{-6}$. The simulation spanned a time interval of 0 to 50 seconds, discretized into 200 equally spaced points ( $\Delta t \approx 0.1 s$ ), balancing computational efficiency with trajectory resolution. The robot's initial state was defined as $[s, \theta, \vartheta, \dot{\theta}]=[ 0,0,1,0]$, corresponding to a starting position at the origin, zero initial orientation, a forward speed of $1 \mathrm{~m} / \mathrm{s}$, and zero angular velocity, aligning with typical agricultural robot initialization.

The robot’s physical parameters, as listed in Table 1, were chosen to reflect a realistic differential drive sprayer with a variable-mass tank.

Three paths were designed to test the robot’s dynamic response and control performance under varying conditions:

• Straight path: Represents linear motion along crop rows, with a reference trajectory defined as $x=t$ and $y=0$, and control inputs $u_1=1 \mathrm{~m} / \mathrm{s}$ and $u_2=0 \mathrm{rad} / \mathrm{s}$.

• Circular path: Emulates navigation around obstacles or circular fields, with a nominal radius of $R=u_1 / u_2=5 \mathrm{~m}$. The reference trajectory was computed via numerical integration of the kinematic equations $\dot{x}=u_1 \cos (\theta), \dot{y}=u_1 \sin (\theta)$, and $\dot{\theta}=u_2$, with $u_1=1 \mathrm{~m} / \mathrm{s}$ and $u_2=0.2 \mathrm{rad} / \mathrm{s}$.

• Zigzag path: Mimics oscillatory motion for wide-area spraying, with $u_1=1 \mathrm{~m} / \mathrm{s}$ and $u_2(t)=0.2 \sin (0.5 t) \mathrm{rad} / \mathrm{s}$. The reference orientation was derived analytically.

The reference trajectories were precomputed over the simulation duration to ensure smoothness and consistency with the kinematic model, addressing previous numerical issues with analytical approximations that led to high tracking errors.

A PID controller was implemented to regulate the quasi-velocities to their desired values for each path, compensating for disturbances from reactive forces, drag, friction, and mass variation. The control law for each quasi-velocity is:

where, $e_i=u_{i, desired}-u_i$ is the tracking error, and $K_p^i, K_d^i$ and $K_i^i$ are the proportional, derivative, and integral gains, respectively. Path-specific gain tuning was employed to optimize performance. An antiwindup mechanism limited integral terms $\pm 1$ to prevent overshoot, addressing earlier high tracking errors.

For each path, the motion equations were solved using the ODE45 solver, with the dynamics implemented. The simulation procedure, as presented in Figure 3, details the iterative process for each path, including initialization, dynamic integration, and validation.

The simulation results demonstrated the efficacy of the extended G-A formulation and the PID control strategy, with significant improvements over initial trials that exhibited high tracking errors and momentum mismatches. Key findings are summarized below, with visual representations proposed in Figures~\ref{fig4}-\ref{fig7}.

Figure 4 represents the mass decrease rate across all paths, with a computed mean of -2.0Kg/s, matching the expected value. The robot's mass decreased linearly from 150 Kg to 110 Kg over 20 seconds, respecting the minimum mass constraint. This validates the variable-mass model's accuracy in capturing tank depletion dynamics.

Figure 5 presents the robot’s axial speed across all paths. As it can be checked, the predefined path restricts the axial velocity to the specified value followed (zigzag and circular paths). However, variation exists for the straight path due to the initial condition which is more with respect to the required value.

Figure 6 discusses the path tracking in comparison to reference paths for three cases. The spray lines, nozzle positions, and start and final points are marked in this figure. The outcomes present that the model works well and the PID control algorithm prepares the inputs as the robot follows the path with the possible least error. It should be considered that these paths follow the changes in robot mass by considering the system constraints simultaneously.

The mean tracking errors and orientation errors for each path are presented in Table 2, with trajectories visualized in Figure 6. The robot achieved a mean tracking error of 0.08 m , reflecting precise adherence to the reference trajectory. The low error confirms the PID controller's ability to maintain constant forward speed under mass depletion and reactive forces. Initial simulations reported a high tracking error of 6.3888 m , indicating significant deviation from the 5 m radius circle. After implementing numerical integration for the reference path and importing PID gains, the tracking error reduced to 0.2346 m . The mean orientation error was 0.0023 rad , consistent with prior results ( 0.0460 rad ), confirming accurate angular velocity tracking. Earlier trials for zigzag exhibited an unacceptable tracking error of 0.35 m . Using an analytical orientation and tuned PID gains, the tracking error decreased to 0.201 m , with a mean orientation error of 0.0523 rad.

As shown in Figure 7, the robot states include the position and orientation as well as the axial and angular velocities presented for three paths. The smooth outcomes guarantee that the derived model for the system works well while the PID control algorithm can reach the good results.

Figure 8 presents the tracking errors for each path along with the corresponding control inputs for axial and angular motion. The main challenges for the similar initial condition exist in the straight path due to the overvalued axial speed which made an error with respect to the set point.

The simulation results validate the extended G-A formulation’s ability to model the complex dynamics of a variable-mass mobile robot sprayer under nonholonomic constraints. The significant reduction in tracking errors highlights the effectiveness of the numerical reference paths and tuned PID control. The accurate momentum balance confirms the correct incorporation of reactive forces, a critical aspect for precision agriculture where spray-induced disturbances affect navigation. The straight path’s low error underscores the model’s reliability for standard spraying tasks, while the circular and zigzag paths demonstrate versatility for complex field maneuvers.

Compared to traditional methods (e.g., Lagrange with multipliers), the extended G-A approach reduced computational complexity by eliminating constraint forces via quasi-velocities. The simulations also addressed practical challenges in precision agriculture, where accurate path tracking ensures uniform pesticide application, enhancing crop yield and resource efficiency. The proposed visualizations provide intuitive insights into the robot’s performance, supporting the theoretical contributions with empirical evidence.

Looking ahead, the proposed extended G-A formulation opens several exciting paths. First, adapting it to handle higher-order nonholonomic constraints could address complex systems like multi-trailer vehicles or soft robots, where constraints involve higher derivatives. Second, incorporating stochastic mass variations—such as irregular tank depletion in sprayers due to uneven terrain—would enhance robustness for real-world applications. Third, experimental validation on a physical robot sprayer is a natural next step, translating the simulations into field tests to confirm accuracy in navigation and spraying. Fourth, integrating the formulation with real-time control algorithms, like model predictive control, could enable autonomous sprayers to adapt dynamically to mass changes, boosting precision in agriculture. Finally, exploring applications beyond robotics, such as spacecraft with variable fuel mass or underwater vehicles with changing payloads, could broaden the method’s impact. These directions not only extend the proposed theoretical framework but also promise practical advancements in automation, sustainability, and dynamics, inviting researchers to push the boundaries of constrained systems.

This study ventured into the intricate realm of constrained variable-mass systems, a domain where traditional dynamic modeling often stumbles. By introducing an extended G-A formulation, a robust solution was offered for deriving motion equations that seamlessly integrate time-varying mass with both holonomic and nonholonomic constraints. The proposed approach redefines the standard G-A framework, leveraging quasi-velocities to eliminate Lagrange multipliers and accommodate the complexities of mass variation. This results in simpler, more efficient equations, a significant step forward for multibody dynamics. The application to a mobile robot sprayer in precision agriculture underscores the practical value of this study. By accurately modeling the robot’s dynamics—accounting for a depleting tank and wheeled constraints—this study enables precise navigation and spraying, enhancing efficiency in crop management. Numerical simulations further validate the proposed method, demonstrating superior accuracy and computational efficiency compared to conventional approaches. The significance of this study lies in its dual contribution: advancing theoretical dynamics and addressing real-world challenges. The extended G-A formulation fills a critical gap in the literature, where variable-mass systems with nonholonomic constraints have been rarely tackled. Its versatility makes it a powerful tool for diverse applications, from robotic manipulators to aerospace systems. In precision agriculture, where automation is revolutionizing farming, the proposed model supports the development of smarter, more efficient robots, contributing to sustainable practices. This bridges the gap between abstract mechanics and tangible societal benefits, a hallmark of impactful research. The proposed extended G-A formulation marks a meaningful advance in modeling complex systems. By tackling the interplay of variable mass and constraints and demonstrating its utility in precision agriculture, this study is expected to inspire further innovation in robotics, dynamics, and beyond.