Final Project 1

Course	Hybrid Dynamics in Mechanical Systems
Course Number	00360087

Ido Fang Bentov
CLASSIFIED
CLASSIFIED

Declaration of independent work

I confirm that this submission reflects my own work and understanding. ChatGPT/Claude/Gemini were used solely to validate algebraic manipulations, and all results were independently reviewed.

MATLAB Code: All related code files can be found on GitHub as well as OneDrive.

---

---

Problem Description

The figure shows the “compass biped” walking robot model. This model consists of two rigid legs connected by a passive revolute joint. The two legs are identical slender rods of mass , length and moment of inertia with respect to their center-of-mass located at their midpoints. Another point mass is located at the passive “hip joint” connecting the two legs. One leg touches the ground, and is called “stance leg”, while the other leg is called “swing leg”. The endpoint of the stance leg, called “stance foot”, makes a unilateral point contact with an inclined plane having slope angle and friction coefficient (both static and dynamic). The stance foot is considered as a passive revolute joint, and may slip on the ground or stick, depending on friction.

The generalized coordinates of the system are , where denote the position of the stance foot, expressed in a fixed reference frame where and axes are aligned with the tangent and normal directions to the inclined plane, respectively. The angles are orientations of the stance and swing legs, respectively, with respect to the axis. Normal contact displacement of the stance foot is constrained as , whereas the tangential displacement satisfies for no-slip contact or during slippage.

The swing leg’s endpoint, called “swing foot”, is located at and its velocity can be written as . The swing foot moves until it passes ahead of the stance leg and hits the ground when . Then it undergoes inelastic impact with friction, and attaches to the ground. As a result of this collision, the rear foot lifts from the ground and the swing and stance feet switch their roles.

Note that during the phase where the swing leg passes across the stance leg, it actually penetrates into the ground for a short time, that is, . This effect, called “foot scuffing”, should be ignored in the analysis and not considered as a collision.

The compass biped passive dynamic walker on an inclined plane.

Sign Convention:

• (stance leg): positive clockwise — leg tips toward (downhill)
• (swing leg): positive counterclockwise — leg tips toward (downhill)

Parameter values for numerical simulations:

---

---

Task 1

Write the equations of motion in matrix form and the matrices , , , , .

Solution:

The positions of the key points in the slope-aligned coordinate system are determined by the generalized coordinates. The stance foot is located at . With measured clockwise from the vertical, the stance leg COM lies at distance along the leg from the foot:

The hip is at distance from the stance foot:

The swing leg COM is at distance from the hip, with measured counterclockwise from the downward vertical:

The swing foot position is:

The velocities are obtained by differentiating the positions with respect to time:

The kinetic energy consists of translational and rotational components for each body. The stance leg contributes , the hip point mass contributes , and the swing leg contributes . The total kinetic energy is .

The gravity vector in the slope-aligned coordinate system is . The potential energy for each body is , giving:

Using the Lagrangian , the equations of motion are derived in the standard form , where are the normal and tangential contact forces. The mass matrix is obtained from the Hessian of the kinetic energy with respect to the generalized velocities, the Coriolis/centrifugal vector is computed using Christoffel symbols, and the gravity vector is the gradient of the potential energy.

From symbolic computation in MATLAB, the resulting matrices are (using shorthand , , , , ):

For stance foot contact, the normal constraint is and the tangential constraint (no-slip) is . The constraint gradients give:

The swing foot velocity is , where the swing foot Jacobian is:

---

---

Task 2

Write two MATLAB functions for numerical simulations of the dynamic equations with no-slip contact:

A function defining the state-space differential equation:

function dXdt = sys_stick(t, X)

A function defining condition for stopping event due to onset of slippage:

[value, isterminal, direction] = events_stick(t, X)

Include also a possible event of “falling” when the hip joint hits the ground.

Solution:

Under no-slip contact, the constraints and reduce the system to two degrees of freedom. The reduced state vector is .

The function sys_stick computes the state derivative for ODE integration. It extracts the angles and angular velocities from the state vector and calls the symbolically-derived stick_accelerations function to compute and , returning .

The function events_stick detects terminating events during simulation. It monitors four conditions: forward slip onset when with increasing , backward slip onset when with decreasing , swing foot collision when the swing foot height with the swing foot ahead of the stance foot (to exclude foot scuffing), and falling when the hip reaches the ground (). The contact forces are computed using the symbolically-derived stick_forces function.

---

---

Task 3

Write MATLAB function of the form Xnew = impact_law(Xold) for calculating the change in the state vector due to collision of the swing foot with the ground. Use Chatterjee’s impact law and assume a frictional fully-plastic collision (). Note that the function should also include the relabeling of coordinates due to interchange between stance and swing legs, so that Xnew represents initial conditions for a new step. The function should also check that the rear foot separates from the ground at the post-impact state, otherwise it should return “failure” code due to double-foot impact.

Solution:

The swing foot collides with the ground when its height reaches zero: . This gives the impact condition , which corresponds to for typical walking configurations where both angles have the same sign. The pre-impact state is .

For a fully-plastic collision (), Chatterjee’s impact law requires the post-impact contact velocity to be zero. The collision matrix relates the contact velocity to the impulse: , where and are evaluated at the impact configuration. The candidate impulse for fully-plastic impact is , where is the pre-impact swing foot velocity.

The impulse must satisfy Coulomb’s friction bound . If violated, the tangential impulse is projected onto the friction cone boundary with where , and is adjusted to ensure :

The post-impact generalized velocities are then .

After impact, the swing foot becomes the new stance foot. Since the positive-rotation conventions for (CW) and (CCW) are opposite, the relabeling includes a sign change:

For a valid step, the old stance foot (now rear foot) must lift off the ground. This requires , where is the common angle at impact. If this condition is violated, the function returns a failure status indicating double-foot impact.

---

---

Task 4

Write MATLAB function of the form Znew = Poincare_map(Zold) for calculating the Poincaré map of the system , where the Poincaré section is chosen as the post-impact state. Assume no-slip motion and impact. The function should return “failure” codes in cases of falling or double-foot impact.

Solution:

The Poincaré section is defined at the post-impact state, where the impact constraint reduces the full 4D state to a 3D reduced state . The function poincare_map takes this reduced state, converts it to the full state , and integrates the no-slip dynamics using ode45 with event detection until the next impact occurs (swing foot height reaches zero with the swing foot ahead of the stance foot). The impact law is then applied to obtain the post-impact velocities, followed by leg relabeling. The new reduced state is extracted and returned along with a status code: 0 for success, 1 for falling, 2 for double-foot impact, and 3 for slip onset.

---

---

Task 5

Use the MATLAB command fsolve for finding a fixed point of the Poincaré map, which satisfies . You can use initial guess from the MSc thesis of Benny Gamus (page 35). If needed, start with parameter values from the thesis and then change them in small increments and update the initial guess for fsolve iteratively, until reaching the current values.

For the periodic solution, present the following graphs, for the time span of a single period:

(a) Plot of the two angles versus time, , overlaid on the same plot. Mark the “scuffing” event on the plot.

(b) Phase plane trajectories of vs. , one curve for each leg in different colors on the same plot. Mark the points of collision by ‘x’ and the points of scuffing by ‘o’ on the curves. Mark the jumps due to impact + foot relabeling in dotted straight-line segments.

(c) Plot of the normal contact force at the stance foot as a function of time , overlaid with a dashed horizontal line at height 0.

(d) Plot of the force ratio at the stance foot as a function of time . At the final time of impact, add a marker ‘x’ for the value of the impulse ratio at the collision.

(e) Plot of the normal height of the swing foot as a function of time . Mark the duration of “scuffing” effect when by two dashed vertical lines at the endpoints of this time interval.

All graphs should have labels and units, and all lines should have width .

Solution:

To find the fixed point, we define the residual function and use MATLAB’s fsolve to find such that . If the Poincaré map fails (due to falling, slip, or double-foot impact), the residual returns a large penalty value to guide the solver away from infeasible regions.

Using the initial guess from Gamus’s thesis, fsolve converges in 4 iterations with residual norm . The fixed point is:

(a) Angles vs time for the periodic solution. Scuffing occurs in the mid-stride region where .

(b) Phase portrait of the periodic solution. Collision points marked with , scuffing points with . Impact jumps (including foot relabeling) shown as dotted lines.

(c) Normal contact force at the stance foot. The force remains positive throughout, confirming sustained contact.

(d) Force ratio with friction bounds . The impact impulse ratio is marked at the end of the period.

(e) Swing foot height during the periodic solution. The scuffing interval where is marked by dashed vertical lines.

---

---

Task 6

Find the minimal value of friction coefficient such that no-slip contact is maintained along the entire periodic solution for all , including the impact. Explain your result.

Solution:

The no-slip periodic solution’s dynamics are independent of — neither the constrained equations of motion nor the fully-plastic no-slip impact law involve the friction coefficient. Therefore, can be determined directly from the periodic solution found in Task 5 by computing the maximum friction demand across both the continuous phase and the impact event.

During the continuous phase, no-slip contact requires at all times. At the impact, the no-slip impulse must satisfy . The minimum friction coefficient ensuring no-slip throughout the entire cycle is therefore:

This corresponds to the peak of graph (d) in Task 5, which shows the force ratio over the period along with the impact impulse ratio. Running the computation yields:

The continuous-phase demand exceeds the impact demand, so the critical instant is at (immediately after impact). Therefore:

This result means that the stance foot is most prone to slipping right at the start of each step, when the post-impact contact forces have the highest tangential-to-normal ratio. For any , the entire periodic walking cycle proceeds without slippage.

---

---