DRS-LIP: Linear Inverted Pendulum Model for Legged Locomotion on Dynamic Rigid Surfaces

Legged robot locomotion on a dynamic rigid surface (i.e., a rigid surface moving in the inertial frame) involves complex full-order dynamics that is high-dimensional, nonlinear, and time-varying. Towards deriving an analytically tractable dynamic model, this study theoretically extends the reduced-order linear inverted pendulum (LIP) model from legged locomotion on a stationary surface to locomotion on a dynamic rigid surface (DRS). The resulting model is herein termed as DRS-LIP. Furthermore, this study introduces an approximate analytical solution of the proposed DRS-LIP that is computationally efficient with high accuracy. To illustrate the practical uses of the analytical results, they are used to develop a hierarchical planning framework that efficiently generates physically feasible trajectories for DRS locomotion. The effectiveness of the proposed theoretical results and motion planner is demonstrated both through simulations and experimentally on a Laikago quadrupedal robot that walks on a rocking treadmill.


I. INTRODUCTION
Legged robots have the potential to traverse various challenging surfaces, including stationary (uneven or discrete) surfaces [1]- [4] and dynamic rigid surfaces (i.e., rigid surfaces that move in the inertial frame) [5], [6]. Common real-world examples of the dynamic rigid surface (DRS) are ships, public transportation vehicles, trains, and elevators. This study aims to model and analyze the essential dynamic behaviors of a legged robot that walks on a DRS by building a reduced-order model and deriving its analytic approximate solution, so as to provide physical insights into the robot dynamics as well as to produce analytic results that can be used for efficient planning of legged locomotion. Yet, reduced-order modeling of DRS locomotion is fundamentally complex due to the nonlinear and high-dimensional nature of legged locomotion dynamics [7] and the timevarying movement of the surface-foot contact points [5].

A. Reduced-order models of stationary surface locomotion
A reduced-order dynamic model of legged locomotion captures the robot's essential dynamic behaviors. One of the most widely studied reduced-order model for stationary surface locomotion is the linear inverted pendulum (LIP) model [8], which models a legged robot as a point mass atop a massless leg. Although the LIP does not capture the complete robot dynamics, many of today's legged robots can be relatively accurately described by the LIP, including bipeds [9], [10] and quadrupeds [4], [11]. This is because they typically have a heavy upper body and lightweight legs.
Due to its simplicity, the LIP model is analytically tractable and can provide physical insights into the essential robot behaviors. It also explicitly reveals the relationship between the center of pressure (CoP), which can infer the feasibility of ground contact forces (i.e., whether the robot rolls about any edge of the region of contact), and the center of mass (CoM). Thus, the LIP can be exploited to represent the robot dynamics in motion planning for ensuring the computational efficiency and physical feasibility of planning. Beside motion planning, the LIP model has also been used as a basis to form balance [12] and stability [9] criteria, and its inherent connection with another widely studied reducedorder model, hybrid zero dynamics [13], has recently been analyzed [14].
The LIP model has been extended to more complex scenarios by considering a varying CoM height [15], CoM movement on a 3-D plane [10], nontrivial centroidal angular momentum [9], and hybrid robot dynamics [16]. Since these models are built upon the assumption that the foot-surface contact point is stationary, they may not be valid for legged locomotion on a DRS, especially when the surface motion is relatively significant.

B. Reduced-order models of dynamic surface locomotion
For legged locomotion on a rigid surface whose motions are affected by the robot (e.g., passive and relatively lightweight surfaces), several reduced-order models have been derived and analyzed, including the LIP [17], [18], centroidal dynamics [19], and rimless-wheel model [20]. Yet, these models may not be valid for robot walking over a rigid surface whose motion is not affected by the robot (e.g., trains, vessels, and elevators). For these substantially heavy or rigidly actuated surfaces, the effects of the surface motion on a spring-loaded inverted pendulum has been numerically studied [21]. Still, it is unclear under what conditions this model is valid for real-world DRS locomotion.
Beyond the scope of legged locomotion, the modeling and analysis of an inverted pendulum with a vertically oscillating support is a classical physics problem. A well-known example is the Kapitza pendulum [22], which has an intriguing property that under high-frequency support oscillation, the pendulum's upper equilibrium becomes stable whereas its lower equilibrium is unstable. However, it remains unknown whether and when the Kapitza pendulum is a reasonable simplification of DRS locomotion. Also, in practical realworld applications, such as legged locomotion on a vessel, the surface motion frequency is typically too low to meet the conditions underlying the Kapitza pendulum.

C. Contributions
This study aims to analytically extend the LIP model [8] to substantially heavy or rigidly actuated DRSes (e.g., vessels and elevators) and to demonstrate the uses of the theoretical results in motion planning. The main contributions are: (a) Deriving a reduced-order model that describes the essential dynamic behaviors of a legged robot during DRS locomotion, by expanding the LIP model to DRS locomotion, which we term as "DRS-LIP". (b) Forming the analytical approximate solution of the DRS-LIP under a vertical, sinusoidal DRS motion, and giving physical insights into the model's stability. (c) Developing a hierarchical planner that utilizes the proposed DRS-LIP and its solution to efficiently generate feasible reference trajectories for DRS locomotion. (d) Assessing the accuracy and computational cost of the proposed solution with MATLAB simulations. (e) Validating the planner efficiency and feasibility through both physics-based simulations and hardware experiments, highlighting the practical usefulness of the proposed analytical results for motion planning. The derivation of the DRS-LIP model, i.e., item (a), has been partly reported in [23]. The new, substantial contributions are items (b)-(e), which are all previously missing.
The paper is organized as follows. Section II introduces the derivation of the proposed DRS-LIP model. Section III presents the analytical approximate solution of the DRS-LIP model under a vertical, sinusoidal surface motion. Section IV develops an efficient motion planner based on the analytical results. Section V reports the validation results. Section VI discusses the capabilities and limitations of the proposed methods. Section VII provides the conclusion remarks.
II. REDUCED-ORDER MODEL OF DRS LOCOMOTION This section introduces a reduced-order model that captures the essential dynamics of a legged robot that walks on a DRS. The model is derived by extending the classical LIP model from static surfaces [8] to a DRS.
Many of today's legged robots have a heavy upper body and lightweight legs. Their center of mass (CoM) dynamics can be approximately described by an LIP, i.e., a point mass atop a massless leg [8], under the assumption that: (A1) The robot's rate of whole-body angular momentum about the CoM is negligible. Assumption (A1) is reasonable for real-world robot locomotion because the trunk is typically controlled to maintain a steady posture for housing sensors (e.g., cameras).
In this study, we use a three-dimensional (3-D) LIP to capture the essential dynamics of a 3-D legged robot walking on a DRS (see Fig. 1). The point mass and support point S in Fig. 1 correspond to the robot's CoM and CoP. Let r wc = [x wc , y wc , z wc ] T and r ws = [x ws , y ws , z ws ] T respectively denote the position of the CoM and point S in the world frame. Then, the CoM position relative to point S, denoted as r sc , is defined as: r sc = r wc − r ws = [x sc , y sc , z sc ] T .
The CoM dynamics during DRS locomotion is given by: x wc = f a x sc mr sin θ ,ÿ wc = f a y sc mr sin θ ,z wc = f a m cos θ − g. (1) Here, m is the total mass, θ is the angle of the vector r sc relative to the vertical axis, g is the norm of the gravitational acceleration vector, and r is the projected length of r sc on the horizontal plane. The scalar variable f a is the norm of the ground contact force.

A. DRS-LIP under a General Vertical Surface Motion
We consider the following assumption on the vertical distance z sc between the CoM and point S (see Fig. 1): (A2) The CoM maintains a constant height z 0 above the support point S (i.e., z sc = z 0 ). This assumption is analogous to the simplifying assumption of the classical LIPM model that the point-mass height over the stationary surface is constant [8].
Under assumption (A2), the relationshipsż wc =ż ws and z wc =z ws hold, and then the axial force f a becomes f a = m(z ws + g)/ cos θ . Thus, the horizontal LIP dynamics is: x wc = (z ws + g) x sc z 0 andÿ wc = (z ws + g) y sc z 0 .
Then, by substitutingẍ wc =ẍ ws +ẍ sc andÿ wc =ÿ ws +ÿ sc into, (2), the horizontal LIP dynamics becomes: When there is no slippage between the support point S and the surface, the acceleration of point S, i.e., (ẍ ws ,ÿ ws ,z ws ), equals the DRS's acceleration at that point. Given that realworld DRSes (e.g., vessels) are typically equipped with high-accuracy, real-time motion monitoring systems [24], we assume the time profile of (ẍ sw ,ÿ sw ,z sw ) is known. Accordingly, they are treated as explicit time functions. Thus, the dynamics in (3) is linear, nonhomogenous, and time-varying.
Since DRSes, such as cruising ships in regular sea waves, have relatively small horizontal acceleration compared with vertical acceleration [25], [26], we assume that the horizontal acceleration of point S is sufficiently small to be ignored: (A3) The horizontal acceleration of point S (i.e.,ẍ ws andÿ ws ) are negligible.
Then, the forcing terms in (3) (i.e., −ẍ ws and −ÿ ws ) can be approximated as zero, and the horizontal LIP dynamics in (3) becomes linear, time-varying, and homogeneous: x sc − (z ws + g) z 0 x sc = 0 andÿ sc − (z ws + g) z 0 y sc = 0. (4) Note that the vertical CoM trajectory is given by z sc = z 0 . Remark 1: The LIP in (4), along with z sc = z 0 , describes the simplified dynamics of a robot walking on a DRS under assumptions (A1)-(A3), which we call "DRS-LIP".

B. DRS-LIP under a Vertical Sinusoidal Surface Motion
A real-world DRS, such as a vessel in regular sea waves, typically exhibits a vertical, sinusoidal motion with a constant amplitude and frequency [26]. Thus, we focus on such motions for further analystis of the DRS-LIP.
Under a vertical, sinusoidal surface motion, the vertical accelerationz ws of point S is sinusoidal, and (4) becomes the well-known Mathieu's equation [27].
Without loss of generality, the vertical sinusoidal motion of the DRS at the surface-foot contact point is assumed as: where the real, scalar parameters A and ω are the amplitude and frequency of the vertical motion, respectively. Therefore, the surface accelerationz ws at the support point is expressed as:z ws := −Aω 2 sin ωt. Substituting the expression ofz ws in (4), we have: As the equations in the x-and y-directions are decoupled and share the same structure, their solutions share the same form. Thus, for brevity, we focus on deriving the solution along the x-direction, x sc , in the next section.
The DRS-LIP in (6) can be transformed into the standard Mathieu's equation [27] by introducing a new time variable τ := π+2ωt 4 and rewriting (6) in terms of τ as: where the real scalar coefficients c 0 and c 1 are defined as c 0 := − 4g ω 2 z 0 and c 1 := 2A z 0 .

A. Computation of Approximate Analytical Solutions
The DRS-LIP model in (6) generally does not have an exact, closed-form analytical solution. One approach to derive an approximate analytical solution is to utilize the fundamental solution matrix based on the Floquet theory [27]. Alternatively, we choose to exploit the existing analytical results of the well-studied Mathieu's equation to obtain a more computationally efficient solution.
There are various existing analytical approximate solutions of Mathieu's equation, including periodic solutions [28] and those expressed through power series [27]. In this study, we adopt the general, exact analytical solution from [29] because of its generality and computationally efficiency: Theorem: The exact, general (periodic or non-periodic) analytical solution of the DRS-LIP in (7) takes the following form: Here, µ is the characteristic exponent of (7). The coefficients α 1 and α 2 are real scalars, the scalar n is an integer, and the coefficients C 2n 's are complex scalars. The proof of Theorem 1 can be readily obtained based on [29]. To use (8) to compute an approximate solution, we need to determine the number of terms to keep in the approximate solution as well as the values of the parameters µ, α 1 , α 2 , and C 2n 's [29], which is explained next.
1) Computing characteristic exponent µ: Substituting the solution (8) into (7) yields a recurrence relationship: where the scalar, complex coefficient function β n is given The derivation of (9) is given in [29] and the supplementary material.
2) Determining number of terms kept: The infinite series in the exact solution (8) is convergent because lim Specifically, C 2n is proportional to β n (µ) and lim n→∞ β n (µ) = 0 (see supplementary material for proof). Thus, we can approximate the exact solution with the sum of finite terms.
With N terms kept, the approximate solution is given by: To ensure a sufficient solution accuracy, we can determine N by numerically evaluating its effect on the accuracy.
3) Computing coefficients C 2n , α 1 , and α 2 : With µ computed, we can determine the value of the coefficients C 2n recursively from (9) by setting C 2N = 0 and C 0 = A [29]. The values of α 1 and α 2 can be obtained based on the initial The last step of determining the approximate analytical solution is to transform the solutionx sc (τ) into the usual time domain t using τ = (π/2 + ωt)/2.

B. Stability Analysis
By the Floquet theory, the DRS-LIP model in (6) is called "stable" if all its solutions are bounded for all t > 0, and is "unstable" if an unbounded solution exists for t > 0. The stability properties of the DRS-LIP model can be determined with the characteristic exponents µ. Since the DRS-LIP is a linear, second-order ordinary differential equation, it has two characteristic exponents, denoted as µ 1 and µ 2 . By the Floquet theory, the model is stable if and only if Re(µ 1 ), Re(µ 2 ) < 0. Suppose that Re(µ 1 ) < Re(µ 2 ). Validation of this physical insight is given in Sec. V.

IV. DRS-LIP BASED HIERARCHICAL PLANNING
To illustrate the practical uses of the DRS-LIP model and its approximate analytical solution, this section presents a hierarchical planner that exploits them to enable efficient and feasible planning of quadrupedal walking on a DRS.
The planner is designed for quadrupedal walking [5], [11] whose complete gait cycle comprises four continuous foot-swinging phases and four discrete foot-landing events, as illustrated in Fig. 3. This planner also assumes the DRS motion is known, which is realistic for real-world DRSes such as vessels because they typically possess motion monitoring systems [31].
The planner has two layers (see Fig. 2). The higher layer produces kinematically and dynamically feasible CoM trajectories for the DRS-LIP model of a legged robot by incorporating necessary feasibility constraints. The lower layer uses trajectory interpolation to translate the CoM trajectories into the desired motion for all degrees of freedom of the full-order robot model.

A. Higher-Layer CoM Trajectory Planner
The higher planner uses the DRS-LIP model to efficiently generate feasible reference CoM position trajectories for the complete gait cycle, by solving the following nonlinear optimization problem.
1) User-defined gait parameters: The input to the higherlayer planner is the user-defined gait parameters (which specify the desired gait features) and the known DRS motion (which is vertical and sinusoidal). The gait parameters are commonly chosen as: (i) average walking velocity (i.e., horizontal CoM velocity), (ii) foot-contact sequence, (iii) stance-foot positions, (iv) constant CoM height above the surface (for respecting assumption (A2)) and (v) gait period. The values of parameters (i)-(iv) are typically set to help ensure a kinematically feasible gait. The value of parameter (v) is selected such that the quotient of the DRS's motion period and the desired gait period is an integer (i.e., the desired CoM motion complies with the DRS motion).
2) Optimization variables: We choose the optimization variables α α α of the planner as the initial horizontal CoM position and velocity of each continuous phase. The rationale for this choice is that given the DRS-LIP model parameters these variables completely determine the horizontal CoM trajectories. The vertical CoM position is not included as an optimization variable because it can be readily obtained from the known DRS motion and the user-defined CoM height.
3) Constraints: We choose to design the constraints to help enforce gait feasibility and to respect the desired gait features as specified by the user-defined parameters. Note that these constraints are formed based on the proposed analytical approximate solution. The equality constraints include: (i) continuity of the CoM trajectories at the footlanding events and (ii) the user-specified walking velocity. The inequality constraints are: (i) friction cone constraint for avoiding foot slipping, (ii) confinement of CoM trajectories within the polygon of support for approximately respecting the CoP constraint, and (iii) upper and lower bounds on the optimization variables α α α.
To meet the constraints, α α α is solved through: where h(α α α) is a scalar cost function (e.g., energy cost of transport), and the vector functions f eq and g ineq are the sets of all equality and inequality constraints, respectively.

B. Lower-Layer Full-Body Trajectory Planner
The lower-layer planning is essentially trajectory interpolation that translates the reference CoM position trajectory supplied by the higher-layer planner into the full-order trajectories of a quadrupedal robot in Cartesian space. To impose a steady trunk/base pose and to avoid swing foot scuffing on the surface, we choose these full-order trajectories to be the absolute base pose and relative swing-foot position with respect to the base.
The input to the lower-layer planner (Fig. 2) is: the known DRS motion that is vertical and sinusoidal, the CoM trajectories provided by the higher-layer planner, and userdefined parameters (including robot's CoM height above the CoP, stance foot locations, and maximum swing foot height). Furthermore, as inspired by previous quadrupedal robot planning [4], a brief four-leg-in-support phase is inserted upon a foot-landing event when the two consecutive polygons of support only share a common edge (i.e., "Switching 1 → 2" and "Switching 3 → 4" in Fig 3, so as to ensure smooth and feasible transitions during these events. 1) Base trajectories: The CoM of the robot is approximated as the base (i.e., the geometric center of the trunk) because a quadruped's trunk typically has a symmetric mass distribution and is substantially heavier than the legs. Thus, the desired base position trajectories are provided by the higher-level planner. To avoid overly stretched leg joints for ensuring kinematic feasibility, the desired base orientation trajectories are designed to comply with the DRS orientation.
2) Swing foot trajectories: The swing foot trajectories are designed to agree with the user-defined stance-foot locations and to respect the kinematic limits of the robot's leg joints. Specifically, we obtain the desired swing foot trajectory during a continuous phase by using Bézier polynomials [5] to connect the adjacent desired stance-foot positions. Remark 3: The dynamic feasibility of the planned trajectories depends on the closeness of the DRS-LIP and the actual robot dynamics. The DRS-LIP model is a relatively faithful representation of an actual DRS-robot system when the system behavior reasonably respects assumptions underlying the proposed model and its solution. Indeed, assumption (A3) holds when the known surface motion is vertical and sinusoidal, and the planner explicitly imposes assumption (A2). Moreover, although the planner enforces the desired base orientation to comply with the surface orientation for kinematic feasibility, assumption (A1) is still reasonably respected by the planned motion. This is because the rate of the robot's centroidal angular momentum is negligible compared with that of its linear momentum under the typical angular movement range of real-world DRSes (e.g., vessels [31]).

V. SIMULATION AND EXPERIMENT RESULTS
This section presents the validation results for the DRS-LIP model, analytical solution, and hierarchical planner.
A. Solution Validation 1) Validation of solution accuracy and efficiency: The accuracy and computational efficiency of the proposed analytical approximate solution in (12) is assessed through comparison with the highly accurate numerical solution. For fairness of comparison, both solutions are computed in MATLAB for the interval t ∈ [0, 0.5] sec. The approximate solution has ten terms kept (i.e., N = 10) for a reasonable trade-off between accuracy and computational efficiency. The comparative numerical solution is computed using MAT-LAB's ODE45 solver with an error tolerance of 10 −9 and at 1000 evenly distributed instants within the given interval.
To validate the proposed solution under different initial conditions, 1000 sets of initial conditions are randomly chosen within |x sc (0)| < 0.2 m and |ẋ sc (0)| < 0.2 m/s. The DRS-LIP model parameters are chosen to be within realistic ranges of DRS motions [26], [31] and quadrupedal robot dimensions: A = 7 cm, ω = π rad/s, and z 0 = 42 cm. Figure 4 shows the accuracy of the approximate analytical solution (with ten terms kept) compared with the numerical  solution for 100 out of the 1000 trials. Within the shown 100 trials, the maximum value of the absolute percentage error is lower than 0.02% in magnitude, indicating the reasonable accuracy of the proposed approximate solution. The absolute percentage error, measured by mean ± one standard deviation, is (0.0012 ± 0.005)% for all 1000 trials. Table I displays the comparison of the average computational time cost (measured by mean±SD) for 1000 trials. The approximate analytical solution is around 15 times faster to compute than the numerical one. 2) Validation of stability property: For typical ship motions in regular sea waves [31], the DRS-LIP model parameters take values as: A ≤ 100 cm and ω ≤ 2π rad/s. Also, the kinematically feasible CoM height z 0 of a typical quadrupedal robot (e.g., Unitree's Laikago) is within [0.3, 0.55] m. Under these parameter ranges, we use (11) to numerically compute the characteristic exponents and obtain that Re(µ 2 ) > 0 and Re(µ 1 ) < 0. Thus, by the Floquet theory, the DRS-LIP is unstable (i.e., an unbounded solution exists) under the considered operating condition. Figure 5 presents the validation results of this physical insight. Subplot (a) displays the analytical approximate solutions under different initial conditions (|x sc (0)| < 0.4 m and |ẋ sc (0)| < 0.4 m/s) and DRS-LIP parameters (ω = π rad/s, A = 7 cm, and z 0 = 42 cm). Subplot (b) shows the solutions under the same initial condition (x sc (0) = 0.02 m anḋ x sc (0) = 0.10 m/s) but different model parameters (0 < ω ≤ 2π rad/s, 0 < A ≤ 100 cm, and 30 ≤ z 0 ≤ 55 cm). In all cases except for the trivial initial condition x sc (0),ẋ sc (0) = 0, the solutions grow towards infinity as time increases, confirming the physical insight that the DRS-LIP is unstable under the considered operating condition.

B. Planner Validation
The efficiency and feasibility of the proposed planner are validated through both simulations and experiments.
1) Simulation and experimental setup: The validation of the planner utilizes a Laikago quadruped (Fig. 6) developed by Unitree Robotics. The dimension of the robot is 55 cm × 35 cm × 60 cm. The robot weighs 25 kg in total. Each leg weighs 2.9 kg, including its motors that weigh 1.5 kg and are located near the trunk. It has twelve independently actuated joints, and the torque limits of hip-roll, hip-pitch, and kneepitch motors are 20 Nm, 55 Nm, and 55 Nm, respectively. Gait parameters. Recall that the proposed planner takes user-defined gait parameters as its input. To assess the planner under different gait parameters, two sets of parameters (G1) and (G2) are used (see Table II). DRS motion. Three DRS motions are tested to assess the effectiveness of the planner under different surface motions: (DRS1) The DRS motion is vertical and sinusoidal with A = 10 cm and ω = π rad/s, which emulates vessel motions in regular sea waves [26]. Motek M-gait treadmill (see Fig. 6). During the testing, the robot is placed approximately 100 cm away from the treadmill's axis of rotation. The treadmill weighs 750 kg with a dimension of 2.3 m × 1.82 m × 0.5 m. A 4.5 kW servo motor powers each of of the treadmill's two belts. The treadmill can be pre-programmed to perform user-defined pitching (but not vertical) motions and belt translation. The belt speed is set to be the same as the desired walking speed. Step length (cm) 10 15 Max. step height (cm) 5 4 2) Validation of planner efficiency: To validate that using the proposed analytical solution improves planner efficiency compared with using the numerical solution, the higher-layer CoM trajectory planning problem is solved based on both solutions under the user-defined gait parameters (G1) and surface motion (DRS2). For simplicity, the cost function h in (13) is chosen as trivial. Also, a 6 th -order Bézier curve is used to design the desired swing-foot trajectory for allowing more freedom in trajectory design.
To demonstrate the improved efficiency under different common solvers, both MATLAB and C++ are used to solve the optimization-based planning problem with the same initial guess for 1000 runs. For fairness of comparison, the optimality and constraint tolerances are set as 10 −6 in all runs. In MATLAB, fmincon is used with an interior-point solver. For the C++ optimization, the nonlinear optimization solver of the Ipopt package [32] is utilized with the same tolerance setting as MATLAB. Table III shows that the time cost of the analytical solution based planning is approximately 5 times shorter than the numerical solution based one. Since the lower-layer planner is essentially trajectory interpolation, solving it is typically fast (e.g., MATLAB can solve it within < 2 ms).
Table III also indicates that the higher-layer planner takes 51.4 ± 3.5 ms to generate the desired CoM trajectory when it is solved by C++ using the approximate anaytical solution. Therefore, the corresponding time cost for solving both higher and lower layers will be less than 60 ms, which is much shorter than the gait period of 2 s. With such a planning speed, the planner will be capable of quickly re-planning the desired full-order trajectories in case of any significant changes in the DRS motion. This planning scheme is reasonable for real-world DRSes, such as vessels [31], because the changes of their motions are relatively slow due to low motion frequency (0.01∼1 Hz). 3) Validation of planner feasibility: Beside efficiency, the proposed DRS-LIP and its solution also help guarantee planning feasibility. To test the feasibility of the planned motion, our previous controller [5], which is derived based on input-output linearization and proportional derivative (PD) control, is utilized to track the planned full-order trajectories in Pybullet simulations and experiments. As this controller does not explicitly ensure the feasibility of ground contact forces, the planned trajectory needs to be physically feasible in order for the controller to be effective. Thus, if the controller is able to reliably track the planned motion and sustain walking on a DRS, then the physical feasibility of the proposed planner is confirmed. To help ensure a reasonable tracking performance, PD gains are tuned as 0.7 and 1.0 in simulations, and 5.5 and 0.15 on hardware.
To validate the planner feasibility under different gait parameters and surface motions, the gait parameters (G1) and the surface motion (DRS1), as well as (G2) and (DRS3), are tested in Pybullet simulations. Figure 7 shows the Pybullet simulation results under (G1) and (DRS1). The robot sustains walking for the entire testing duration, which is over 50 gait cycles. The base and joint trajectories closely track their reference values, as shown in subplots (a) and (b). Also, the actual robot motion respects the torque limits, as indicated in subplot (c). The simulation results under (G2) and (DRS3) show position tracking performance and torque profiles similar to the case under (G1) and (DRS1), and is given in supplementary material for space consideration.
The combination of (G1) and (DRS2) is tested in both simulations and experiments, with the results presented in Fig. 8. Experiment videos are included in the paper submission, and is also available at https://youtu.be/Amrb3WxLpPs. In both simulations and experiments, the robot walking is stable, as indicated by the trajectory tracking accuracy in subplots (a) and (b) as well as the experiment video. Moreover, as displayed in subplot (c), the joint torque limits are met in both the simulation and experiment. However, the torque profiles of the right leg's three joints show notable discrepancies between Pybullet and experiment results, which could be caused by the differences between the simulated and actual robot dynamics as well as the different inherent meanings of their effective PD gains. Also, from the experiment video, we can see that the robot experiences relatively notable rebounding and slipping at contact switching events when a rear leg lands on the surface. This violation of the planned contact sequence is essentially due to the robot's temporary loss of contact force feasibility, which could be mitigated through controller design as discussed in Sec. VI.
VI. DISCUSSION This paper has introduced a reduced-order dynamic model of a legged robot that walks on a DRS, by analytically extending the classical LIP model from stationary surfaces [8] to a DRS (e.g., a vessel). The resulting DRS-LIP model in (4) is a linear, second-order differential equation, similar to the classical LIP. However, the DRS-LIP is explicitly time-varying whereas the classical LIP is time-invariant. This fundamental difference is due to the time-varying movement of the surface at the surface-foot contact points. This study also provides the stability condition of the DRS-LIP based on the Floquet theory. Analogous to the classical  LIP [8], the DRS-LIP is unstable under the usual movement range of real-world DRSes such as vessels [31] (Sec. V-A).
The DRS-LIP is valid under the assumption that the actual robot's rate of whole-body angular momentum about the CoM is negligible (assumption (A1)). To relax this assumption, the point mass of the proposed DRS-LIP could be augmented with a fly wheel [9], [10] to account for nonzero rate of angular momentum. Also, the DRS-LIP can be generalized from a constant CoM height (as enforced by assumption (A2)) to a varying height by integrating with the variable-height LIP for stationary surfaces [15].
This study also derives the approximate analytical solution of the DRS-LIP for vertical, sinusoidal surface motions. Its improved accuracy and computational efficiency compared with numerical solutions are confirmed through MATLAB simulations, as shown in Fig. 4 and Table I. Yet, the solution is not valid for periodic surface motions that are not strictly sinusoidal or vertical. Under general periodic surface motions, the DRS-LIP in (4) becomes Hill's equation [27], and thus the solution could be derived based on existing results of Hill's equations. Also, as indicated by (3), the DRS-LIP could become nonhomogeneous under nontrivial horizontal DRS motions, with forcing terms induced by the horizontal DRS accelerations. The effects of horizontal DRS motions on robot's movement could be studied based on (3).
To highlight the usefulness of the analytical results, they have been used as a basis to synthesize a hierarchical planner that efficiently produces desire, physically feasible motions for quadrupedal robot walking on a DRS. The feasibility of the planned motion is validated by using our previous trajectory tracking controller [5] to command a quadrupedal robot to follow the planned motion during DRS walking. As discussed in Sec.V-B, both simulation and experimental results indicate the reasonable feasibility of the proposed planner under different gait parameters and surface motions (see Figs 7 and 8). To mitigate the temporary violation of the planned gait sequence observed in experiments, which is partly induced by the discrepancies between the DRS-LIP and the actual robot dynamics, the planned motion could be tracked by an optimization-based controller that explicitly ensures physical feasibility [23].
VII. CONCLUSION This paper has introduced a reduced-order model that captures the essential dynamic behavior of a legged robot that walks on a nonstationary rigid surface (i.e., a DRS). The model was derived based on the theoretical extension of the classical LIP model from a stationary surface to a DRS. The model was found to be unstable under the common operating conditions of a real-world DRS (e.g., a vessel in regular sea waves). Its approximate analytical solution was also provided and used as a basis to develop a motion planner that efficienctly generates feasible robot motions for quadrupedal walking on a DRS. The stability property of the model and the efficiency and accuracy of its solution were validated through MATLAB and C++ simulations.

DERIVATION DETAILS OF PROPOSED ANALYTICAL APPROXIMATE SOLUTION
This section gives the derivation details for the proposed analytical approximate solution that are omitted in Sec. III of the main paper.

Recurrence Relationship between µ and β n
This subsection presents the derivation of the recurrence relationship between µ and β n , which is given in Eq. (9) in the main paper. The original derivation is given in the reference [29] cited in the main paper.

SIMULATION VALIDATION UNDER GAIT PARAMETERS (G2) AND SURFACE MOTION (DRS2)
As presented in Sec.V-B-3 of the main paper, three combinations of gait parameters and surface motions are tested to validate the feasibility of the desired full-order position trajectories produced by the proposed planner. These combinations are: parameter (G1) and motion (DRS1), (G1) and (DRS2), and (G3) and (DRS2). The results from the first two cases are shown in Figs. 7 and 8 in the main paper. The Pybullet simulation results from the last case is shown in Fig. 1 of this supplementary material. Similar to the first case, the position and orientation tracking of base and swing foot trajectories is accurate (subplots (a) and (b)), and the torque limit is respected (subplot (c)). Thus, the robot is able to satisfactorily track the planned motion under a controller that does not explicitly ensure gait feasibility, indicating the feasibility of the planned motion. Recordings of the Pybullet simulation are included in the video submission.