Skip to main content

Chapter 10: Physics Simulations for Robotics

10.1 Physics Engine Fundamentals

10.1.1 Introduction to Physics Simulation

Physics simulation forms the backbone of modern robotics testing and development, enabling virtual prototypes to behave according to real-world physical laws. In robotics, physics engines provide the computational framework for simulating object interactions, dynamics, and environmental constraints that govern robot behavior.

Physics simulation in robotics differs from game physics in its emphasis on accuracy and predictability. While games prioritize visual appeal and performance, robotics simulations require precise physical modeling for reliable algorithm development and testing.

10.1.2 Core Physics Components

A comprehensive physics simulation system consists of several interconnected components:

Rigid Body Dynamics

Rigid body dynamics govern the motion of solid objects that do not deform under stress. The fundamental equations of motion for rigid bodies are:

M * dv/dt = F_ext + F_contact$$ Where: - M = Mass matrix (including inertia tensor) - v = Velocity vector (linear and angular) - F_ext = External forces (gravity, user-applied forces) - F_contact = Contact and constraint forces #### Collision Detection Collision detection determines when and where objects intersect in the simulated environment. This occurs in two phases: **Broad Phase**: Quickly eliminates pairs that cannot possibly collide using spatial partitioning structures. ```python # Broad phase collision detection using AABB tree class AABBNode: def __init__(self, bounds, objects=None): self.bounds = bounds # (min_x, min_y, min_z, max_x, max_y, max_z) self.objects = objects or [] self.left = None self.right = None def potential_collisions(self): """Returns pairs of objects that might collide""" if self.is_leaf(): return [(obj1, obj2) for i, obj1 in enumerate(self.objects) for obj2 in self.objects[i+1:]] left_pairs = self.left.potential_collisions() right_pairs = self.right.potential_collisions() cross_pairs = [(obj1, obj2) for obj1 in self.left.get_all_objects() for obj2 in self.right.get_all_objects() if self.aabb_intersect(obj1.aabb, obj2.aabb)] return left_pairs + right_pairs + cross_pairs ``` **Narrow Phase**: Precise collision detection using geometric algorithms. <Equation> d(P1, P2) = ||P1 - P2|| < r1 + r2 </Equation> Where: - P1, P2 = Object positions - r1, r2 = Object radii (for sphere-sphere collision) #### Constraint Solving Constraints enforce relationships between objects, such as joints, motors, and contact constraints. The constraint solver ensures these relationships are maintained while respecting physical laws. ```cpp // Position constraint for a distance joint class DistanceJoint : public Constraint { private: RigidBody* body1; RigidBody* body2; Vector3 anchor1, anchor2; // Anchor points in local space float rest_length; public: void solve(float dt) override { Vector3 world_anchor1 = body1->transform_point(anchor1); Vector3 world_anchor2 = body2->transform_point(anchor2); Vector3 delta = world_anchor2 - world_anchor1; float current_length = delta.length(); if (current_length > 0.0f) { Vector3 normal = delta / current_length; float position_error = current_length - rest_length; // Calculate correction impulses Vector3 correction = normal * position_error * 0.5f; body1->position += correction; body2->position -= correction; } } }; ``` ### 10.1.3 Physics Engine Architecture Modern physics engines use a pipeline architecture to process simulation steps: ```cpp class PhysicsEngine { private: std::vector<RigidBody*> bodies; std::vector<Constraint*> constraints; std::unique_ptr<BroadPhase> broad_phase; std::unique_ptr<NarrowPhase> narrow_phase; std::unique_ptr<ConstraintSolver> solver; float time_accumulator = 0.0f; public: void step(float delta_time) { time_accumulator += delta_time; const fixed_dt = 1.0f / 60.0f; // 60 Hz physics update while (time_accumulator >= fixed_dt) { // Simulation pipeline apply_forces(fixed_dt); integrate_velocities(fixed_dt); // Collision detection auto contacts = detect_collisions(); // Constraint solving solve_constraints(contacts, fixed_dt); // Position integration integrate_positions(fixed_dt); time_accumulator -= fixed_dt; } } }; ``` <Diagram> Physics Engine Pipeline: 1. Force Application → 2. Velocity Integration → 3. Collision Detection → 4. Constraint Solving → 5. Position Integration </Diagram> ## 10.2 Collision Detection and Response ### 10.2.1 Collision Detection Algorithms #### Sphere-Sphere Collision The simplest collision detection case, perfect for wheeled robots and simple sensors: ```python def sphere_sphere_collision(sphere1, sphere2): """Detect collision between two spheres""" distance_vector = sphere1.position - sphere2.position distance = distance_vector.length() combined_radius = sphere1.radius + sphere2.radius if distance < combined_radius: # Collision detected penetration = combined_radius - distance normal = distance_vector / distance if distance > 0 else Vector3(1, 0, 0) return CollisionInfo( normal=normal, penetration=penetration, contact_point=sphere1.position - normal * sphere1.radius ) return None ``` #### Box-Box Collision (OBB) Oriented Bounding Box (OBB) collision is crucial for robotic arms and manipulators:

|n·(c2 - c1)| ≤ |n·a1| + |n·a2|$$

Where:

  • n = Separating axis
  • c1, c2 = Box centers
  • a1, a2 = Box extents along axis
def obb_collision(box1, box2):
    """Separating Axis Theorem for OBB collision"""
    # Get all potential separating axes
    axes = []
    axes.extend(box1.axes)  # Box1 local axes
    axes.extend(box2.axes)  # Box2 local axes

    # Cross products of axes
    for a1 in box1.axes:
        for a2 in box2.axes:
            cross = a1.cross(a2)
            if cross.length() > 0.0001:  # Avoid parallel axes
                axes.append(cross.normalized())

    # Test each axis
    min_overlap = float('inf')
    min_axis = None

    for axis in axes:
        proj1 = box1.project_onto_axis(axis)
        proj2 = box2.project_onto_axis(axis)

        overlap = min(proj1.max, proj2.max) - max(proj1.min, proj2.min)

        if overlap <= 0:
            return None  # No collision

        if overlap < min_overlap:
            min_overlap = overlap
            min_axis = axis

    return CollisionInfo(normal=min_axis, penetration=min_overlap)

Mesh-Mesh Collision

Complex geometries require triangle mesh collision detection:

class TriangleMesh {
public:
    struct Triangle {
        Vector3 vertices[3];
        Vector3 normal;

        Vector3 get_support(const Vector3& direction) const {
            Vector3 support = vertices[0];
            float max_dot = direction.dot(support);

            for (int i = 1; i < 3; ++i) {
                float dot = direction.dot(vertices[i]);
                if (dot > max_dot) {
                    max_dot = dot;
                    support = vertices[i];
                }
            }
            return support;
        }
    };

private:
    std::vector<Triangle> triangles;
    std::unique_ptr<BVH> bvh;  // Bounding Volume Hierarchy

public:
    CollisionResult collide_with(const TriangleMesh& other) {
        // Use GJK algorithm for convex meshes
        return gjk_collision(*this, other);
    }
};

10.2.2 Collision Response

Once collisions are detected, the physics engine must calculate appropriate responses:

Impulse-Based Response

Impulse-based collision response simulates the instantaneous change in velocity during collision:

J = -(1 + e) * v_rel · n / (1/m1 + 1/m2 + (r1×n)²/I1 + (r2×n)²/I2)$$ Where: - J = Impulse magnitude - e = Coefficient of restitution (bounciness) - v_rel = Relative velocity - n = Collision normal - m1, m2 = Object masses - r1, r2 = Contact points relative to centers of mass - I1, I2 = Moments of inertia ```python def apply_collision_impulse(body1, body2, contact): """Apply collision impulse between two bodies""" # Calculate relative velocity at contact point v1 = body1.velocity + body1.angular_velocity.cross(contact.r1) v2 = body2.velocity + body2.angular_velocity.cross(contact.r2) v_rel = v1 - v2 # Don't resolve if velocities are separating velocity_along_normal = v_rel.dot(contact.normal) if velocity_along_normal > 0: return # Calculate restitution (bounciness) e = min(body1.restitution, body2.restitution) # Calculate impulse scalar j = -(1 + e) * velocity_along_normal j /= 1/body1.mass + 1/body2.mass # Apply impulse impulse = j * contact.normal body1.velocity += impulse / body1.mass body2.velocity -= impulse / body2.mass # Apply angular impulse body1.angular_velocity += contact.r1.cross(impulse) / body1.inertia body2.angular_velocity -= contact.r2.cross(impulse) / body2.inertia ``` #### Position Correction To prevent objects from sinking into each other, position correction (also called "baumgarte stabilization") is applied: ```cpp void PositionConstraint::solve(float dt) { Vector3 world_anchor1 = body1->transform_point(local_anchor1); Vector3 world_anchor2 = body2->transform_point(local_anchor2); Vector3 delta = world_anchor2 - world_anchor1; float distance = delta.length(); if (distance > 0.0f) { Vector3 normal = delta / distance; float position_error = distance - rest_length; // Baumgarte stabilization float baumgarte = 0.2f; // Stability factor float correction = position_error * baumgarte; Vector3 correction_vector = normal * correction * 0.5f; body1->position += correction_vector; body2->position -= correction_vector; } } ``` ## 10.3 Robot-Specific Physics ### 10.3.1 Wheel and Track Physics Wheeled robots require specialized physics for accurate ground interaction: #### Rolling Resistance Rolling resistance opposes wheel motion and depends on surface deformation:

F_rr = C_rr * m * g * cos(θ)$$

Where:

  • F_rr = Rolling resistance force
  • C_rr = Coefficient of rolling resistance
  • m = Mass
  • g = Gravitational acceleration
  • θ = Incline angle
class Wheel {
    def __init__(self, radius, mass, friction_coeff=0.8):
        self.radius = radius
        self.mass = mass
        self.friction_coeff = friction_coeff
        self.angular_velocity = 0.0
        self.steering_angle = 0.0

    def update_physics(self, torque, dt, ground_normal):
        """Update wheel physics based on applied torque"""
        # Calculate linear force from torque
        linear_force = torque / self.radius

        # Apply rolling resistance
        normal_force = self.mass * 9.81 * ground_normal.y
        rolling_resistance = 0.01 * normal_force  # C_rr ≈ 0.01 for wheels

        # Net force considering friction
        max_friction = self.friction_coeff * normal_force
        net_force = max(-max_friction, min(max_friction, linear_force - rolling_resistance))

        # Update velocities
        angular_acceleration = net_force / (self.mass * self.radius)
        self.angular_velocity += angular_acceleration * dt

        # Return linear velocity at wheel center
        return self.angular_velocity * self.radius

Slip and Traction

Wheel slip occurs when the applied force exceeds available traction:

F_max = μ * N$$ Where: - F_max = Maximum available traction force - μ = Coefficient of friction - N = Normal force ```python def calculate_wheel_slip(wheel, applied_torque, normal_force): """Calculate wheel slip ratio""" # Theoretical linear velocity (no slip) theoretical_velocity = wheel.angular_velocity * wheel.radius # Maximum traction force max_traction = wheel.friction_coeff * normal_force # Required force for applied torque required_force = applied_torque / wheel.radius if abs(required_force) <= max_traction: # No slip slip_ratio = 0.0 else: # Slip occurs if wheel.actual_velocity > 0: slip_ratio = (theoretical_velocity - wheel.actual_velocity) / wheel.actual_velocity else: slip_ratio = 1.0 # Full slip return slip_ratio ``` ### 10.3.2 Joint and Actuator Physics Robotic joints require complex physics modeling including constraints, motors, and compliance: #### Revolute Joint Physics Revolute joints constrain rotation to a single axis: ```cpp class RevoluteJoint : public Constraint { private: RigidBody* body1; RigidBody* body2; Vector3 anchor; // Joint position in world space Vector3 axis; // Rotation axis float lower_limit; // Lower angle limit float upper_limit; // Upper angle limit float motor_torque; // Applied motor torque float motor_speed; # Target motor speed public: void solve_velocity(float dt) override { // Get anchor points in world space Vector3 r1 = anchor - body1->position; Vector3 r2 = anchor - body2->position; // Calculate angular velocities around joint axis float w1 = body1->angular_velocity.dot(axis); float w2 = body2->angular_velocity.dot(axis); float relative_w = w1 - w2; // Motor control if (motor_torque != 0.0f) { float motor_impulse = motor_torque * dt; // Apply motor torque body1->apply_torque(axis * motor_impulse); body2->apply_torque(axis * -motor_impulse); } // Apply angular constraint float angular_error = relative_w * 0.5f; # Soft constraint Vector3 correction = axis * angular_error; body1->angular_velocity -= correction; body2->angular_velocity += correction; } }; ``` #### PD Controller for Joint Control PD (Proportional-Derivative) controllers are commonly used for joint position control:

τ = K_p * (θ_target - θ_current) + K_d * (0 - ω_current)$$

class PDController:
    def __init__(self, kp, kd):
        self.kp = kp  # Proportional gain
        self.kd = kd  # Derivative gain
        self.prev_error = 0.0

    def calculate_torque(self, target_position, current_position, current_velocity, dt):
        """Calculate control torque using PD control"""
        error = target_position - current_position

        # Proportional term
        p_term = self.kp * error

        # Derivative term
        if dt > 0:
            derivative = (error - self.prev_error) / dt
            self.prev_error = error
        else:
            derivative = 0.0

        d_term = self.kd * (-current_velocity)  # Damping

        # Total torque
        torque = p_term + d_term

        # Clamp torque to motor limits
        max_torque = 10.0  # N⋅m
        torque = max(-max_torque, min(max_torque, torque))

        return torque

10.3.3 Sensor Physics

Sensors interact with the physics environment to generate realistic readings:

Lidar Physics Simulation

Lidar sensors require ray casting and surface interaction modeling:

class LidarSensor:
    def __init__(self, position, rotation, num_rays=360, max_range=10.0):
        self.position = position
        self.rotation = rotation
        self.num_rays = num_rays
        self.max_range = max_range
        self.noise_std = 0.01  # 1cm standard deviation

    def cast_rays(self, physics_world):
        """Cast rays and return distance measurements"""
        measurements = []

        for i in range(self.num_rays):
            # Calculate ray direction
            angle = (i / self.num_rays) * 2 * math.pi
            direction = Vector3(
                math.cos(angle),
                0,
                math.sin(angle)
            ).rotate(self.rotation)

            # Cast ray through physics world
            hit_distance = self._cast_ray(physics_world, direction)

            # Add sensor noise
            noisy_distance = hit_distance + random.gauss(0, self.noise_std)

            measurements.append({
                'angle': angle,
                'distance': max(0, min(noisy_distance, self.max_range))
            })

        return measurements

    def _cast_ray(self, physics_world, direction):
        """Cast single ray and return hit distance"""
        # Use physics engine's raycast functionality
        hit_result = physics_world.raycast(
            origin=self.position,
            direction=direction,
            max_distance=self.max_range
        )

        if hit_result.hit:
            return hit_result.distance
        else:
            return self.max_range

Camera Physics Simulation

Camera sensors require optical physics modeling for realistic images:

class CameraSensor:
    def __init__(self, position, rotation, fov=60, resolution=(640, 480)):
        self.position = position
        self.rotation = rotation
        self.fov = fov  # Field of view in degrees
        self.resolution = resolution
        self.focal_length = resolution[0] / (2 * math.tan(math.radians(fov/2)))

    def render_frame(self, physics_world):
        """Render camera view using ray tracing"""
        image = np.zeros((*self.resolution, 3), dtype=np.uint8)

        # Generate rays for each pixel
        for y in range(self.resolution[1]):
            for x in range(self.resolution[0]):
                # Convert pixel to ray direction
                ray_direction = self._pixel_to_ray(x, y)

                # Cast ray through scene
                color = self._trace_ray(physics_world, ray_direction)

                image[y, x] = color

        return image

    def _pixel_to_ray(self, x, y):
        """Convert pixel coordinates to ray direction"""
        # Normalize pixel coordinates
        nx = (x - self.resolution[0]/2) / self.focal_length
        ny = -(y - self.resolution[1]/2) / self.focal_length  # Negative because image y is inverted

        # Create ray direction in camera space
        ray_direction = Vector3(nx, ny, 1.0).normalized()

        # Transform to world space
        return ray_direction.rotate(self.rotation)

    def _trace_ray(self, physics_world, direction, max_bounces=3):
        """Trace ray through scene and return color"""
        if max_bounces <= 0:
            return self._get_background_color(direction)

        # Cast ray
        hit = physics_world.raycast(self.position, direction)

        if hit.hit:
            # Calculate lighting
            light_direction = Vector3(1, 1, 1).normalized()

            # Check if in shadow
            shadow_ray = physics_world.raycast(
                hit.position + hit.normal * 0.001,
                light_direction
            )

            if shadow_ray.hit:
                light_intensity = 0.3  # Ambient light
            else:
                light_intensity = max(0, hit.normal.dot(light_direction)) * 0.7 + 0.3

            # Apply material properties
            object_color = hit.object.get_color()
            color = object_color * light_intensity

            # Add reflection for shiny objects
            if hit.object.reflectivity > 0:
                reflect_dir = direction - hit.normal * 2 * direction.dot(hit.normal)
                reflect_color = self._trace_ray(
                    physics_world,
                    reflect_dir,
                    max_bounces - 1
                )
                color = color * (1 - hit.object.reflectivity) + reflect_color * hit.object.reflectivity

            return (color * 255).astype(np.uint8)

        return self._get_background_color(direction)

10.4 Optimization Techniques

10.4.1 Spatial Partitioning

Efficient collision detection requires spatial data structures to reduce pairwise checks:

Uniform Grid

Simple grid-based partitioning for uniform object distributions:

class UniformGrid:
    def __init__(self, world_size, cell_size):
        self.world_size = world_size
        self.cell_size = cell_size
        self.grid_size = (
            int(world_size[0] / cell_size),
            int(world_size[1] / cell_size),
            int(world_size[2] / cell_size)
        )
        self.cells = {}

    def insert_object(self, obj):
        """Insert object into grid cells"""
        cells = self._get_object_cells(obj)
        for cell_key in cells:
            if cell_key not in self.cells:
                self.cells[cell_key] = []
            self.cells[cell_key].append(obj)

    def get_potential_collisions(self):
        """Get all potential collision pairs"""
        collision_pairs = set()

        for cell_objects in self.cells.values():
            # Check pairs within cell
            for i, obj1 in enumerate(cell_objects):
                for obj2 in cell_objects[i+1:]:
                    if id(obj1) < id(obj2):  # Avoid duplicate pairs
                        collision_pairs.add((obj1, obj2))

        return list(collision_pairs)

    def _get_object_cells(self, obj):
        """Get all grid cells occupied by object"""
        min_cell = self._world_to_grid(obj.bounds.min)
        max_cell = self._world_to_grid(obj.bounds.max)

        cells = []
        for x in range(min_cell[0], max_cell[0] + 1):
            for y in range(min_cell[1], max_cell[1] + 1):
                for z in range(min_cell[2], max_cell[2] + 1):
                    cells.append((x, y, z))

        return cells

Dynamic Bounding Volume Tree (DBVT)

Hierarchical structure for efficient broad-phase collision detection:

class DBVT {
    struct Node {
        AABB bounds;
        Node* left;
        Node* right;
        Object* object;  // Leaf nodes only

        bool is_leaf() const { return object != nullptr; }
    };

private:
    Node* root = nullptr;

public:
    void insert(Object* obj) {
        Node* new_node = create_leaf_node(obj);

        if (!root) {
            root = new_node;
        } else {
            Node* sibling = find_best_sibling(new_node);
            Node* parent = create_internal_node(sibling, new_node);
            replace_node(sibling, parent);
            refit_hierarchy(parent);
        }
    }

    void remove(Object* obj) {
        Node* leaf = find_leaf(obj);
        if (leaf == root) {
            root = nullptr;
        } else {
            Node* parent = leaf->parent;
            Node* grandparent = parent->parent;
            Node* sibling = (parent->left == leaf) ? parent->right : parent->left;

            if (grandparent) {
                replace_node(parent, sibling);
                refit_hierarchy(grandparent);
            } else {
                root = sibling;
                sibling->parent = nullptr;
            }

            destroy_node(parent);
            destroy_node(leaf);
        }
    }

    std::vector<std::pair<Object*, Object*>> query_overlaps() {
        std::vector<std::pair<Object*, Object*>> pairs;
        query_overlaps_recursive(root, pairs);
        return pairs;
    }

private:
    Node* find_best_sibling(Node* leaf) {
        Node* best_node = root;
        float best_cost = std::numeric_limits<float>::max();

        // Traverse tree to find best insertion location
        std::stack<Node*> stack;
        stack.push(root);

        while (!stack.empty()) {
            Node* current = stack.top();
            stack.pop();

            if (current->is_leaf()) {
                continue;
            }

            // Calculate insertion cost
            AABB combined = leaf->bounds;
            combined.combine(current->bounds);
            float cost = combined.area() - current->bounds.area();

            if (cost < best_cost) {
                best_cost = cost;
                best_node = current;
            }

            // Continue traversal
            stack.push(current->left);
            stack.push(current->right);
        }

        return best_node;
    }
};

10.4.2 Parallel Processing

Modern physics engines leverage multi-core processors for parallel computation:

Task-Based Parallelism

Break simulation into independent tasks for parallel execution:

import multiprocessing as mp
from concurrent.futures import ThreadPoolExecutor

class ParallelPhysicsEngine:
    def __init__(self, num_workers=None):
        self.num_workers = num_workers or mp.cpu_count()
        self.bodies = []
        self.constraints = []

    def step(self, dt):
        """Parallel physics step"""
        with ThreadPoolExecutor(max_workers=self.num_workers) as executor:
            # Parallel force application
            futures = []
            for body in self.bodies:
                future = executor.submit(self._apply_forces, body, dt)
                futures.append(future)

            # Wait for force application
            for future in futures:
                future.result()

            # Parallel velocity integration
            futures = []
            for body in self.bodies:
                future = executor.submit(self._integrate_velocity, body, dt)
                futures.append(future)

            # Wait for velocity integration
            for future in futures:
                future.result()

            # Collision detection (parallel broad phase)
            contacts = self._parallel_collision_detection()

            # Sequential constraint solving (requires global state)
            self._solve_constraints(contacts, dt)

            # Parallel position integration
            futures = []
            for body in self.bodies:
                future = executor.submit(self._integrate_position, body, dt)
                futures.append(future)

            # Wait for position integration
            for future in futures:
                future.result()

    def _apply_forces(self, body, dt):
        """Apply external forces to body"""
        # Gravity
        body.force += Vector3(0, -9.81 * body.mass, 0)

        # Other forces (drag, user forces, etc.)
        body.force += body.drag_force

        # Clear torque if needed
        body.torque = Vector3(0, 0, 0)

    def _integrate_velocity(self, body, dt):
        """Integrate velocity using semi-implicit Euler"""
        # Linear velocity
        acceleration = body.force / body.mass
        body.velocity += acceleration * dt

        # Angular velocity
        angular_acceleration = body.inertia_tensor.inverse() * body.torque
        body.angular_velocity += angular_acceleration * dt

        # Apply damping
        body.velocity *= (1 - body.linear_damping * dt)
        body.angular_velocity *= (1 - body.angular_damping * dt)

    def _parallel_collision_detection(self):
        """Parallel collision detection"""
        # Broad phase can be parallelized
        with ThreadPoolExecutor(max_workers=self.num_workers) as executor:
            chunk_size = len(self.bodies) // self.num_workers
            futures = []

            for i in range(0, len(self.bodies), chunk_size):
                end = min(i + chunk_size, len(self.bodies))
                chunk = self.bodies[i:end]

                future = executor.submit(self._broad_phase_chunk, chunk)
                futures.append(future)

            # Collect potential collisions
            potential_pairs = set()
            for future in futures:
                potential_pairs.update(future.result())

        # Narrow phase (can also be parallelized)
        with ThreadPoolExecutor(max_workers=self.num_workers) as executor:
            futures = []

            for pair in potential_pairs:
                future = executor.submit(self._narrow_phase, pair)
                futures.append(future)

            contacts = []
            for future in futures:
                contact = future.result()
                if contact:
                    contacts.append(contact)

        return contacts

10.4.3 Level of Detail (LOD)

Reduce computational complexity by using simpler physics for distant objects:

class LODPhysicsObject:
    def __init__(self, high_detail_obj, medium_detail_obj, low_detail_obj):
        self.high_detail = high_detail_obj
        self.medium_detail = medium_detail_obj
        self.low_detail = low_detail_obj
        self.current_detail = "high"
        self.switch_threshold_high = 10.0  # Distance units
        self.switch_threshold_low = 50.0

    def update_lod(self, camera_position):
        """Update level of detail based on distance"""
        distance = (self.position - camera_position).length()

        if distance < self.switch_threshold_high:
            self.current_detail = "high"
        elif distance < self.switch_threshold_low:
            self.current_detail = "medium"
        else:
            self.current_detail = "low"

    def get_physics_object(self):
        """Get current physics object based on LOD"""
        if self.current_detail == "high":
            return self.high_detail
        elif self.current_detail == "medium":
            return self.medium_detail
        else:
            return self.low_detail

    def should_collide_with(self, other):
        """High-detail objects only collide with nearby high-detail objects"""
        if self.current_detail == "high" and other.current_detail == "high":
            distance = (self.position - other.position).length()
            return distance < self.switch_threshold_high * 2
        return False

10.5 Performance Analysis and Tuning

10.5.1 Profiling Physics Performance

Measure and analyze physics simulation performance:

class PhysicsProfiler:
    def __init__(self):
        self.metrics = {
            'force_application': [],
            'collision_detection': [],
            'constraint_solving': [],
            'integration': [],
            'total_frame_time': []
        }

    def profile_step(self, physics_engine, dt):
        """Profile a single physics step"""
        import time

        total_start = time.time()

        # Profile force application
        start = time.time()
        physics_engine.apply_forces(dt)
        self.metrics['force_application'].append(time.time() - start)

        # Profile collision detection
        start = time.time()
        contacts = physics_engine.detect_collisions()
        self.metrics['collision_detection'].append(time.time() - start)

        # Profile constraint solving
        start = time.time()
        physics_engine.solve_constraints(contacts, dt)
        self.metrics['constraint_solving'].append(time.time() - start)

        # Profile integration
        start = time.time()
        physics_engine.integrate(dt)
        self.metrics['integration'].append(time.time() - start)

        self.metrics['total_frame_time'].append(time.time() - total_start)

    def get_performance_report(self):
        """Generate performance report"""
        report = {}
        for metric, values in self.metrics.items():
            if values:
                report[metric] = {
                    'average': sum(values) / len(values),
                    'max': max(values),
                    'min': min(values),
                    'std_dev': self._calculate_std_dev(values)
                }
        return report

    def identify_bottlenecks(self):
        """Identify performance bottlenecks"""
        report = self.get_performance_report()
        bottlenecks = []

        # Find the slowest stage
        slowest_stage = None
        slowest_time = 0

        for metric, stats in report.items():
            if metric != 'total_frame_time' and stats['average'] > slowest_time:
                slowest_time = stats['average']
                slowest_stage = metric

        if slowest_stage:
            bottlenecks.append(f"Bottleneck: {slowest_stage} averaging {slowest_time:.4f}s")

        # Check if exceeding real-time requirements
        if report.get('total_frame_time', {}).get('average', 0) > 0.016:  # 60 FPS = 16ms
            bottlenecks.append("Cannot maintain 60 FPS - performance optimization needed")

        return bottlenecks

10.5.2 Tuning Guidelines

Optimize physics performance through parameter tuning:

Time Step Optimization

dt_optimal = min(dt_max, 1/f_target)$$ ```python class AdaptiveTimeStep: def __init__(self, base_dt=1.0/60.0): self.base_dt = base_dt self.current_dt = base_dt self.max_dt = 1.0/30.0 # Minimum 30 FPS self.min_dt = 1.0/120.0 # Maximum 120 FPS self.target_frame_time = base_dt def adapt_time_step(self, last_frame_time, num_objects): """Adapt time step based on performance""" # Calculate frame time ratio ratio = last_frame_time / self.target_frame_time # Adjust time step if ratio > 1.1: # Running too slow self.current_dt = max( self.min_dt, self.current_dt * 0.9 ) elif ratio < 0.9: # Running fast self.current_dt = min( self.max_dt, self.current_dt * 1.1 ) # Consider scene complexity complexity_factor = min(1.0, num_objects / 1000.0) self.current_dt *= (1.0 - complexity_factor * 0.3) return self.current_dt ``` #### Solver Iteration Tuning Balance accuracy vs performance in constraint solving: ```python class TunableConstraintSolver: def __init__(self): self.base_iterations = 10 self.max_iterations = 50 self.min_iterations = 4 self.error_threshold = 0.01 self.adaptive_enabled = True def solve_constraints(self, constraints, dt): """Solve with adaptive iteration count""" if not self.adaptive_enabled: return self._solve_fixed_iterations(constraints, dt, self.base_iterations) # Start with minimum iterations iterations = self.min_iterations previous_error = float('inf') while iterations <= self.max_iterations: # Solve with current iteration count error = self._solve_and_measure_error(constraints, dt, iterations) # Check convergence if error < self.error_threshold: break # Check if iterations are helping if abs(error - previous_error) < 0.001: break # Converged previous_error = error iterations += 2 # Increase by 2 for efficiency return iterations def _solve_and_measure_error(self, constraints, dt, iterations): """Solve and return final error""" for _ in range(iterations): for constraint in constraints: constraint.solve(dt) # Calculate constraint violation total_error = 0 for constraint in constraints: total_error += constraint.get_violation() return total_error / len(constraints) if constraints else 0 ``` ## Chapter Summary This chapter covered the essential aspects of physics simulation for robotics applications: ### Key Concepts Covered 1. **Physics Engine Architecture**: Pipeline-based processing with force application, collision detection, constraint solving, and integration 2. **Collision Detection**: Broad phase spatial partitioning and narrow phase geometric algorithms 3. **Collision Response**: Impulse-based resolution with position correction for stability 4. **Robot-Specific Physics**: Specialized models for wheels, joints, actuators, and sensors 5. **Optimization Techniques**: Spatial data structures, parallel processing, and LOD systems 6. **Performance Tuning**: Adaptive time steps and solver iteration management ### Practical Implementations - Complete physics engine pipeline implementation - Wheel slip and traction modeling - PD controller for joint actuation - Lidar and camera sensor simulation - Parallel collision detection using threading - Performance profiling and bottleneck identification ### Next Steps With a solid understanding of physics simulation, you're ready to explore: - Chapter 11: NVIDIA Isaac Sim & Synthetic Data Generation - Chapter 12: Digital Twin Development - Advanced simulation techniques for specific robot types --- ## Glossary Terms **Term**: **Broad Phase Collision Detection** **Definition**: Initial collision filtering stage that quickly eliminates non-colliding object pairs using spatial partitioning or simple bounding volume tests **Related**: **Narrow Phase**, **Bounding Volume Hierarchy** **Term**: **Constraint Solver** **Definition**: Algorithm that enforces physical constraints (joints, contacts, motors) while respecting conservation laws and physical principles **Related**: **Impulse-Based Response**, **Position-Based Dynamics** **Term**: **Level of Detail (LOD)** **Definition**: Technique that reduces computational complexity by using simpler representations for distant or less important objects **Related**: **Spatial Partitioning**, **Performance Optimization** **Term**: **Semi-Implicit Euler Integration** **Definition**: Numerical integration method that updates velocity before position, providing better stability than explicit Euler while maintaining efficiency **Related**: **Runge-Kutta**, **Verlet Integration** **Term**: **Baumgarte Stabilization** **Definition**: Position correction technique that prevents constraint drift by applying corrective forces proportional to position errors **Related**: **Constraint Solving**, **Position-Based Dynamics** --- ## Exercises ### Exercise 10.1: Basic Physics Engine Implement a simple 2D physics engine with: - Circle-circle collision detection and response - Basic gravity and damping - Position integration using semi-implicit Euler - Visual verification of physics behavior ### Exercise 10.2: Wheel Slip Model Create a simulation demonstrating wheel slip: - Implement wheel-ground interaction - Show slip ratio vs traction force relationship - Visualize slip conditions under different torques - Compare with real-world slip curves ### Exercise 10.3: Sensor Simulation Build a sensor physics simulation: - Implement ray-based lidar with noise - Create camera simulation with basic lighting - Test sensor performance in different environments - Validate sensor readings against ground truth ### Exercise 10.4: Performance Optimization Optimize a physics simulation: - Implement uniform grid spatial partitioning - Profile before and after optimization - Measure performance gains with increasing object count - Document optimization trade-offs ### Exercise 10.5: Adaptive Time Stepping Create adaptive time step simulation: - Implement frame rate-based time step adaptation - Test stability under varying load conditions - Compare fixed vs adaptive time stepping - Analyze accuracy vs performance trade-offs