Pick and place a cube with a scripted grasp in PyBullet
Define a five-joint arm with a parallel gripper as an inline URDF, solve its inverse kinematics in closed form and script a grasp in PyBullet that moves a 40 mm cube to within 1 mm of its target.
| Runs on | Laptop CPU (macOS, Linux or Windows), PyBullet in DIRECT mode, under 5 seconds |
|---|---|
| Time | about 3 hours |
| Stack | python, numpy, pybullet, urdf |
What you will build
This project writes a complete pick-and-place program for a simulated arm, with nothing borrowed from a model library. You define the robot yourself: a five-joint arm with a two-finger parallel gripper, as a URDF string of boxes and cylinders. You solve its inverse kinematics in closed form so that the gripper always points straight down. Then you script the classic sequence: approach, close, lift, move, place, release and retreat. The cube starts rotated 17 degrees from the robot’s axes, so the wrist has to turn to line the fingers up with its faces.
PyBullet runs headless (p.DIRECT), with its plane.urdf as the floor. The program prints a line per phase and finishes with the cube’s final position and a pass or fail against a 10 mm tolerance. In our run the cube lands 0.6 mm from the target, upright, after 10.8 s of simulated time.
It helps to have done the two-link arm project first, because this arm’s IK reuses the same law-of-cosines solution in a vertical plane. The URDF and MJCF cheat sheet covers the model format, and the manipulation and grasping hub puts scripted grasps in context.
What you need
- Python 3.10 or later with
pybulletandnumpy(pip install pybullet numpy). We ran it with PyBullet 3.2.5 and NumPy 2.5.4 on an Apple Silicon Mac mini. pybullet_data, which ships inside thepybulletpackage. Onlyplane.urdfis used from it.- No GPU and no display. The run takes a few seconds of CPU.
Design
The robot
The arm is a common layout for small manipulators. A base yaw joint turns the whole arm about the vertical. Three pitch joints (shoulder, elbow, wrist) about the turret’s y axis keep the arm in one vertical plane. A wrist roll about the hand’s own axis lets the fingers line up with an object at any yaw. The gripper has two prismatic fingers that slide along the hand’s y axis. A fixed joint places a massless-looking tcp link (1 g, so PyBullet does not invent an inertia for it) 55 mm below the hand plate, between the finger pads.
Every movable joint in the URDF has a <limit> with lower and upper bounds, effort and velocity. The reference URDF parser, urdfdom, rejects a revolute or prismatic joint without one. Every link with collision geometry has an <inertial> block, because PyBullet gives a link with no inertial data a mass of 1 kg and unit inertia and prints a warning. The geometry is collision shapes only. With p.DIRECT nothing is drawn, so visual shapes would only cost loading time.
Closed-form IK with the gripper pointing down
Pick and place from above needs the tcp at a position (x, y, z), the hand axis pointing straight down, and the fingers at a chosen yaw. That is five constraints for five joints, so the solution is closed form:
- Base yaw. q1 = atan2(y, x) turns the arm’s plane through the target.
- Wrist position. The hand points down, so the wrist pitch axis is L_hand = 0.095 m straight above the tcp. In the arm’s plane the wrist is at horizontal distance r = √(x² + y²) and height h = z + 0.095 - 0.20 relative to the shoulder.
- Shoulder and elbow. This is the two-link problem again, with links of 0.30 m and 0.25 m. The angles are measured from the vertical, because the arm points straight up at zero. The elbow angle comes from the law of cosines and the shoulder from two
atan2calls. The code always takes the branch with the elbow above the wrist. - Wrist pitch. The three pitch angles add up to the hand’s angle from the vertical. For the hand to point down that sum must be π, so q4 = π - q2 - q3.
- Wrist roll. Pointing the hand down flips it over. With the hand flipped, its yaw in the world becomes q1 - q5 rather than q1 + q5, a sign that is easy to get wrong. We got it wrong in the first draft. The fingers closed at the wrong angle, met the cube near its corners and twisted it as they gripped. A cube looks the same every 90 degrees, so the code wraps q5 into the range ±45 degrees and never rolls the wrist further than needed.
move_tcp interpolates the tcp along a straight line in Cartesian space with a smooth start and stop (3s² - 2s³), solves IK at every 1/240 s step and sends the joint angles to PyBullet’s position controllers. Straight-line moves matter on the way down and up: a joint-space move between the hover and grasp poses would swing the fingers sideways through the cube.
The grasp sequence
The cube is 40 mm on a side and 0.1 kg. It starts at (0.42, -0.12) m with a yaw of 0.3 rad and goes to (0.25, 0.35) m. The tcp hovers at 0.15 m and grasps at 0.03 m, which leaves the finger tips 5 mm above the floor. To close, the program commands each finger to the centre line with a force limit of 20 N. The fingers stop on the cube faces and keep squeezing at that limit. With a friction coefficient of 1.0 on the fingers and the cube, the friction available is far more than the 0.98 N weight of the cube. On release the tcp is 5 mm higher than at the grasp, so the cube drops about 5 mm onto the floor rather than being pushed into it.
Coupled fingers
Two independently force-limited fingers will squeeze the cube, but nothing makes them squeeze it in the middle. One finger arrives first and pushes the cube sideways until the other meets it. Then, during the move, the pair can drift sideways together, because each finger is only holding a force. Our first version did exactly this. The cube was 7 mm off-centre after the lift and 15 mm off at the place position, and the run failed the 10 mm tolerance.
Real parallel grippers avoid this with a mechanism, a rack and pinion or a linkage, that ties the fingers together. The code does the same with a PyBullet gear constraint between the two finger joints, with a ratio of -1. Both joint axes point outwards, so the constraint forces the two openings to be equal. After that change the cube sits 0.4 mm off-centre after the lift and lands within 1 mm of the target.
Code
The whole program is one file. Save it as pick_place.py.
# pick_place.py: a 5-joint arm with a parallel gripper picks a cube and places it, in PyBullet (DIRECT mode).
import math
import os
import tempfile
import numpy as np
import pybullet as p
import pybullet_data
H_SHOULDER = 0.20 # m, shoulder axis height above the floor
L_UPPER, L_FORE = 0.30, 0.25
L_HAND = 0.095 # m, wrist pitch axis to the tool centre point (between the finger pads)
ARM_URDF = f"""<?xml version="1.0"?>
<robot name="pick_arm">
<link name="base_link">
<inertial><origin xyz="0 0 0.05"/><mass value="3.0"/><inertia ixx="0.01" iyy="0.01" izz="0.01" ixy="0" ixz="0" iyz="0"/></inertial>
<visual><origin xyz="0 0 0.05"/><geometry><cylinder radius="0.08" length="0.1"/></geometry></visual>
<collision><origin xyz="0 0 0.05"/><geometry><cylinder radius="0.08" length="0.1"/></geometry></collision>
</link>
<link name="turret">
<inertial><origin xyz="0 0 0.05"/><mass value="0.5"/><inertia ixx="0.001" iyy="0.001" izz="0.001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0 0.05"/><geometry><cylinder radius="0.05" length="0.1"/></geometry></collision>
</link>
<link name="upper_arm">
<inertial><origin xyz="0 0 0.15"/><mass value="0.5"/><inertia ixx="0.004" iyy="0.004" izz="0.0002" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0 0.15"/><geometry><box size="0.04 0.04 {L_UPPER}"/></geometry></collision>
</link>
<link name="forearm">
<inertial><origin xyz="0 0 0.125"/><mass value="0.4"/><inertia ixx="0.002" iyy="0.002" izz="0.0001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0 0.125"/><geometry><box size="0.035 0.035 {L_FORE}"/></geometry></collision>
</link>
<link name="wrist">
<inertial><origin xyz="0 0 0.02"/><mass value="0.1"/><inertia ixx="0.0001" iyy="0.0001" izz="0.0001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0 0.02"/><geometry><cylinder radius="0.025" length="0.04"/></geometry></collision>
</link>
<link name="hand">
<inertial><origin xyz="0 0 0.01"/><mass value="0.1"/><inertia ixx="0.0001" iyy="0.0001" izz="0.0001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0 0.01"/><geometry><box size="0.03 0.12 0.02"/></geometry></collision>
</link>
<link name="finger_left">
<inertial><origin xyz="0 0.005 0.03"/><mass value="0.02"/><inertia ixx="0.00001" iyy="0.00001" izz="0.00001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 0.005 0.03"/><geometry><box size="0.02 0.01 0.06"/></geometry></collision>
</link>
<link name="finger_right">
<inertial><origin xyz="0 -0.005 0.03"/><mass value="0.02"/><inertia ixx="0.00001" iyy="0.00001" izz="0.00001" ixy="0" ixz="0" iyz="0"/></inertial>
<collision><origin xyz="0 -0.005 0.03"/><geometry><box size="0.02 0.01 0.06"/></geometry></collision>
</link>
<link name="tcp">
<inertial><mass value="0.001"/><inertia ixx="1e-7" iyy="1e-7" izz="1e-7" ixy="0" ixz="0" iyz="0"/></inertial>
</link>
<joint name="base_yaw" type="revolute">
<parent link="base_link"/><child link="turret"/><origin xyz="0 0 0.1"/><axis xyz="0 0 1"/>
<limit lower="-3.14" upper="3.14" effort="50" velocity="2.0"/>
</joint>
<joint name="shoulder" type="revolute">
<parent link="turret"/><child link="upper_arm"/><origin xyz="0 0 0.1"/><axis xyz="0 1 0"/>
<limit lower="-1.6" upper="1.6" effort="50" velocity="2.0"/>
</joint>
<joint name="elbow" type="revolute">
<parent link="upper_arm"/><child link="forearm"/><origin xyz="0 0 {L_UPPER}"/><axis xyz="0 1 0"/>
<limit lower="-2.6" upper="2.6" effort="30" velocity="2.0"/>
</joint>
<joint name="wrist_pitch" type="revolute">
<parent link="forearm"/><child link="wrist"/><origin xyz="0 0 {L_FORE}"/><axis xyz="0 1 0"/>
<limit lower="-2.6" upper="2.6" effort="10" velocity="3.0"/>
</joint>
<joint name="wrist_roll" type="revolute">
<parent link="wrist"/><child link="hand"/><origin xyz="0 0 0.04"/><axis xyz="0 0 1"/>
<limit lower="-3.14" upper="3.14" effort="5" velocity="3.0"/>
</joint>
<joint name="finger_left_joint" type="prismatic">
<parent link="hand"/><child link="finger_left"/><origin xyz="0 0 0.02"/><axis xyz="0 1 0"/>
<limit lower="0.0" upper="0.05" effort="20" velocity="0.1"/>
</joint>
<joint name="finger_right_joint" type="prismatic">
<parent link="hand"/><child link="finger_right"/><origin xyz="0 0 0.02"/><axis xyz="0 -1 0"/>
<limit lower="0.0" upper="0.05" effort="20" velocity="0.1"/>
</joint>
<joint name="tcp_joint" type="fixed">
<parent link="hand"/><child link="tcp"/><origin xyz="0 0 0.055"/>
</joint>
</robot>
"""
ARM_JOINTS = [0, 1, 2, 3, 4] # base_yaw, shoulder, elbow, wrist_pitch, wrist_roll
FINGERS = [5, 6]
TCP = 7
DT = 1.0 / 240.0
def arm_ik(x: float, y: float, z: float, grasp_yaw: float) -> list[float]:
"""Analytic IK with the gripper pointing straight down and its fingers closing across grasp_yaw + 90 deg."""
q1 = math.atan2(y, x)
r = math.hypot(x, y) # horizontal distance in the arm's vertical plane
h = z + L_HAND - H_SHOULDER # wrist pitch axis height relative to the shoulder
c3 = (r * r + h * h - L_UPPER ** 2 - L_FORE ** 2) / (2 * L_UPPER * L_FORE)
if abs(c3) > 1:
raise ValueError("target out of reach")
q3 = math.acos(c3) # elbow bent forward, elbow above the wrist
q2 = math.atan2(r, h) - math.atan2(L_FORE * math.sin(q3), L_UPPER + L_FORE * math.cos(q3))
q4 = math.pi - q2 - q3 # hand axis points down
# With the hand flipped to point down, its yaw in the world is q1 - q5.
# A cube looks the same every 90 deg, so pick the smallest roll that lines the fingers up with its faces.
q5 = q1 - grasp_yaw
q5 = (q5 + math.pi / 4) % (math.pi / 2) - math.pi / 4
return [q1, q2, q3, q4, q5]
sim_steps = 0
def step(n: int = 1) -> None:
global sim_steps
for _ in range(n):
p.stepSimulation()
sim_steps += 1
def set_gripper(robot: int, opening: float, force: float) -> None:
"""opening is the distance of each finger from the centre line, in metres."""
p.setJointMotorControlArray(robot, FINGERS, p.POSITION_CONTROL, targetPositions=[opening] * 2, forces=[force] * 2)
def move_tcp(robot: int, start: np.ndarray, goal: np.ndarray, yaw: float, seconds: float) -> None:
"""Straight-line Cartesian move of the tool centre point, solving IK at every control step."""
n = int(seconds / DT)
for i in range(1, n + 1):
s = i / n
s = 3 * s * s - 2 * s ** 3 # smooth start and stop
q = arm_ik(*((1 - s) * start + s * goal), yaw)
p.setJointMotorControlArray(robot, ARM_JOINTS, p.POSITION_CONTROL, targetPositions=q,
forces=[50, 50, 30, 10, 5], positionGains=[0.3] * 5)
step()
step(60) # 0.25 s to settle
def tcp_pos(robot: int) -> np.ndarray:
return np.array(p.getLinkState(robot, TCP, computeForwardKinematics=True)[4])
p.connect(p.DIRECT)
p.setAdditionalSearchPath(pybullet_data.getDataPath())
p.setGravity(0, 0, -9.81)
p.setTimeStep(DT)
p.loadURDF("plane.urdf")
urdf_path = os.path.join(tempfile.mkdtemp(), "pick_arm.urdf")
with open(urdf_path, "w") as f:
f.write(ARM_URDF)
robot = p.loadURDF(urdf_path, useFixedBase=True)
# Couple the fingers like the linkage in a real parallel gripper, so the cube is always centred between them.
gear = p.createConstraint(robot, FINGERS[0], robot, FINGERS[1], p.JOINT_GEAR, [0, 1, 0], [0, 0, 0], [0, 0, 0])
p.changeConstraint(gear, gearRatio=-1, erp=0.1, maxForce=50)
CUBE = 0.04 # m, edge length
pick = np.array([0.42, -0.12, CUBE / 2])
place = np.array([0.25, 0.35, CUBE / 2])
cube_yaw = 0.3 # rad, the cube is not square to the robot
col = p.createCollisionShape(p.GEOM_BOX, halfExtents=[CUBE / 2] * 3)
cube = p.createMultiBody(baseMass=0.1, baseCollisionShapeIndex=col, basePosition=pick.tolist(),
baseOrientation=p.getQuaternionFromEuler([0, 0, cube_yaw]))
for body, link in [(cube, -1), (robot, FINGERS[0]), (robot, FINGERS[1])]:
p.changeDynamics(body, link, lateralFriction=1.0)
# Start above the cube with the gripper open.
GRASP_Z = 0.03 # TCP height when grasping: finger tips 5 mm above the floor
HOVER_Z = 0.15
home = np.array([pick[0], pick[1], HOVER_Z])
for j, q in zip(ARM_JOINTS, arm_ik(*home, cube_yaw)):
p.resetJointState(robot, j, q)
set_gripper(robot, 0.035, 20)
step(120)
print(f"start: tcp at {np.round(tcp_pos(robot), 4)}, cube at {np.round(p.getBasePositionAndOrientation(cube)[0], 4)}")
# Approach, close, lift, move, place, release, retreat.
above_pick, at_pick = home, np.array([pick[0], pick[1], GRASP_Z])
above_place, at_place = np.array([place[0], place[1], HOVER_Z]), np.array([place[0], place[1], GRASP_Z + 0.005])
move_tcp(robot, above_pick, at_pick, cube_yaw, 1.0)
print(f"approach: tcp error {np.linalg.norm(tcp_pos(robot) - at_pick) * 1000:.1f} mm")
set_gripper(robot, 0.0, 20) # aim past the cube faces so the fingers squeeze with 20 N
step(240)
grip = p.getJointState(robot, FINGERS[0])[0] + p.getJointState(robot, FINGERS[1])[0]
contacts = len(p.getContactPoints(robot, cube))
print(f"close: finger gap {grip * 1000:.1f} mm (cube is {CUBE * 1000:.0f} mm), contact points {contacts}")
move_tcp(robot, at_pick, above_pick, cube_yaw, 1.0)
offset = np.array(p.getBasePositionAndOrientation(cube)[0]) - tcp_pos(robot)
print(f"lift: cube centre {-offset[2] * 1000:.1f} mm below the tcp, {np.linalg.norm(offset[:2]) * 1000:.1f} mm off-centre sideways")
move_tcp(robot, above_pick, above_place, cube_yaw, 2.0)
move_tcp(robot, above_place, at_place, cube_yaw, 1.0)
set_gripper(robot, 0.035, 20)
step(240)
print("release: gripper open")
move_tcp(robot, at_place, above_place, cube_yaw, 1.0)
step(240)
pos, orn = p.getBasePositionAndOrientation(cube)
pos = np.array(pos)
roll, pitch, yaw = p.getEulerFromQuaternion(orn)
xy_err = np.linalg.norm(pos[:2] - place[:2])
TOL = 0.01
print(f"cube final position: x={pos[0]:.4f} y={pos[1]:.4f} z={pos[2]:.4f} m")
print(f"cube final tilt: roll {math.degrees(roll):.2f} deg, pitch {math.degrees(pitch):.2f} deg, yaw {math.degrees(yaw):.1f} deg")
print(f"xy error {xy_err * 1000:.1f} mm, tolerance {TOL * 1000:.0f} mm: {'PASS' if xy_err < TOL and abs(pos[2] - CUBE / 2) < 0.005 else 'FAIL'}")
print(f"simulated {sim_steps * DT:.1f} s in {sim_steps} steps")
p.disconnect()
The URDF is written to a temporary file because p.loadURDF reads from disk. useFixedBase=True bolts base_link to the world. In PyBullet the link index of a child link is the same as the index of the joint that creates it. For this model p.getJointInfo lists base_yaw as joint 0 through to tcp_joint as joint 7, which is why TCP = 7. Check the indices with p.getJointInfo whenever you change a model.
Run
python -m pip install pybullet numpy
python pick_place.py
On our Mac the run takes about 0.13 s. PyBullet prints its build time first, on standard error. That line is not an error, and its date depends on your PyBullet build.
Test
This is the output from our run. PyBullet in DIRECT mode is deterministic on one machine, so you should see the same lines. On another platform the last digit of some values may differ.
pybullet build time: Oct 21 2025 16:30:45
start: tcp at [ 0.42 -0.12 0.15], cube at [ 0.42 -0.12 0.02]
approach: tcp error 0.0 mm
close: finger gap 39.9 mm (cube is 40 mm), contact points 6
lift: cube centre 10.5 mm below the tcp, 0.4 mm off-centre sideways
release: gripper open
cube final position: x=0.2506 y=0.3503 z=0.0200 m
cube final tilt: roll 0.00 deg, pitch -0.00 deg, yaw 17.2 deg
xy error 0.6 mm, tolerance 10 mm: PASS
simulated 10.8 s in 2580 steps
Line by line:
- start. The arm is reset onto the IK solution above the cube, and the tcp is where the IK said it would be.
- approach. After the 1 s descent and 0.25 s of settling, the tcp error rounds to 0.0 mm.
- close. The finger gap of 39.9 mm against a 40 mm cube shows the fingers are pressing on the faces. PyBullet allows a small overlap at a contact. The contact count confirms that the robot touches the cube.
- lift. The cube came up with the hand. Its centre is 10.5 mm below the tcp (grasp height 30 mm, cube centre 20 mm), and only 0.4 mm off-centre sideways, so it did not slip.
- final position and tilt. The cube rests on the floor (z = 0.0200 m is half its size), with zero roll and pitch. Its yaw is unchanged at 17.2 degrees, the 0.3 rad it started with.
- PASS. The horizontal error is 0.6 mm, well inside the 10 mm tolerance.
If you get FAIL, look at the lift line first. A large sideways offset means the fingers are not coupled. A cube that is not below the tcp at all means the grasp failed, usually because of a sign error in the wrist roll or low friction.
On a real robot
The sequence is the same on hardware, and the sequence is the easy part. Here is what changes.
- You do not know where the cube is. In simulation the cube’s pose is a variable. On a robot it comes from a camera and a perception step, with errors of millimetres and a degree or two. Make the gripper opening wide enough to tolerate that, and check the grasp with the finger position or a force reading before you lift.
- The arm does not track perfectly. Real joints have backlash, compliance and calibration errors. The repeatability that datasheets quote measures how closely the arm returns to a taught point. It says nothing about how closely the arm reaches a point computed from a model. Teach the pick and place points on the robot if the task allows.
- Gripper force is a choice. 20 N holds a 0.1 kg cube with a large margin, but it would crush a paper cup. Many electric grippers let you set the force. Pick the lowest force that holds the part under the accelerations of the move.
- Collisions are real. Add a collision check against the table and fixtures before every move. Lower speeds near the workspace boundaries, and keep the emergency stop in reach the first time the program runs.
The block below shows the shape of the same sequence on hardware. It needs a physical robot. RobotArm and Gripper stand in for your robot’s own driver or SDK, and the method names are placeholders. Every vendor’s interface is different, so it is not runnable as written.
# needs-hardware
from my_robot_driver import Gripper, RobotArm # placeholder for your vendor's SDK
HOVER_Z, GRASP_Z = 0.15, 0.03 # m
arm, gripper = RobotArm("192.168.1.10"), Gripper()
def pick_and_place(pick_xyz, place_xyz, yaw_rad):
gripper.open()
arm.move_linear(pick_xyz[0], pick_xyz[1], HOVER_Z, yaw_rad, speed=0.1)
arm.move_linear(pick_xyz[0], pick_xyz[1], GRASP_Z, yaw_rad, speed=0.02)
gripper.close(force=20.0)
if not gripper.holding_object():
raise RuntimeError("grasp failed, stopping before the lift")
arm.move_linear(pick_xyz[0], pick_xyz[1], HOVER_Z, yaw_rad, speed=0.05)
arm.move_linear(place_xyz[0], place_xyz[1], HOVER_Z, yaw_rad, speed=0.1)
arm.move_linear(place_xyz[0], place_xyz[1], GRASP_Z + 0.005, yaw_rad, speed=0.02)
gripper.open()
arm.move_linear(place_xyz[0], place_xyz[1], HOVER_Z, yaw_rad, speed=0.05)
The vertical moves are slower than the horizontal ones. That is the usual practice near contact, and a grasp check sits between the close and the lift.
Things to try
- Uncouple the fingers. Delete the two gear-constraint lines and run again. In our run the cube landed 15.1 mm from the target and the program printed
FAIL. - Less friction. Set the lateral friction to 0.1 and then 0.05, and find the value at which the cube slips out during the lift. Compare it with the force you would expect from 20 N per finger.
- Use PyBullet’s IK. Replace
arm_ikwithp.calculateInverseKinematics(robot, TCP, target, targetOrientation)and a downward-pointing quaternion. PyBullet quaternions are (x, y, z, w). Compare the tcp error after the approach. - Move faster. Halve every move time. Watch the lift line for slip and the approach line for tracking error, then raise
positionGainsto compensate. - Stack two cubes. Add a second cube and place it on the first, at a place height 40 mm higher.
- Random starts. Draw the cube’s start position and yaw from a seeded random generator inside the reachable area, and count how many of 50 runs pass.
Summary
- A URDF of boxes and cylinders, with
<limit>on every movable joint and<inertial>on every link, is enough for a working arm and gripper in PyBullet. - With the gripper pointing straight down, five joints give a closed-form IK: base yaw, a two-link solution for shoulder and elbow, wrist pitch = π - q2 - q3, and a wrist roll that must account for the flipped hand.
- Move the tcp in straight Cartesian lines near the object, and solve IK at every control step.
- Couple the fingers of a parallel gripper, in simulation with a gear constraint and on hardware with a mechanism. Two independently force-limited fingers let the object drift sideways.
- Print what you can check: tcp error, finger gap, slip after the lift and the final pose against a tolerance.
The PyBullet Quickstart Guide documents loadURDF, setJointMotorControlArray, createConstraint and the contact queries used here. The URDF glossary entry and the simulation and digital twins hub cover the model format and simulator choice.