In sim-to-real transfer we train policies in simulation and then deploy them on physical robots. A key failure mode (see Lesson 1: Why Policies Fail in Reality) is the reality gap: differences in dynamics, sensing, and appearance between simulator and real world. Domain randomization addresses this by replacing a single, fixed simulator with a distribution of simulators.

Let the simulator be parameterized by a vector \( \boldsymbol{\phi} \in \Phi \subset \mathbb{R}^d \) collecting masses, frictions, contact coefficients, sensor delays, and rendering parameters. Domain randomization chooses, at the start of each episode, a random parameter vector \( \boldsymbol{\phi} \sim p(\boldsymbol{\phi}) \) and trains a policy that performs well on average across this distribution.

\[ \text{Train policy } \pi_{\boldsymbol{\theta}} \text{ to maximize } J_{\text{DR}}(\boldsymbol{\theta}) = \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[\, J(\boldsymbol{\theta}, \boldsymbol{\phi}) \,\right], \]

where \( J(\boldsymbol{\theta}, \boldsymbol{\phi}) \) is the expected return of policy \( \pi_{\boldsymbol{\theta}} \) on the simulator with parameters \( \boldsymbol{\phi} \). Intuitively, the policy must hedge against many plausible worlds, making it more likely to succeed in the true, unknown real-world parameter setting.

flowchart TD
  A["Sample parameters phi ~ p(phi)"] --> B["Reset simulator with phi"]
  B --> C["Run episode with policy pi_theta"]
  C --> D["Collect trajectories and rewards"]
  D --> E["Update theta using RL/optimizer"]
  E --> A["Repeat until convergence"]
        

We now connect domain randomization to the Markov Decision Process (MDP) formulation introduced in Chapter 12 (Reinforcement Learning for Robotics). A parametric MDP is a family \( \{ \mathcal{M}_{\boldsymbol{\phi}} : \boldsymbol{\phi} \in \Phi \} \) where each MDP is

\[ \mathcal{M}_{\boldsymbol{\phi}} = (\mathcal{S}, \mathcal{A}, P_{\boldsymbol{\phi}}, r_{\boldsymbol{\phi}}, \gamma, \rho_{\boldsymbol{\phi}}). \]

Here:

• \( \mathcal{S} \): state space (joint positions, velocities, etc.),
• \( \mathcal{A} \): action space (torques, velocities, etc.),
• \( P_{\boldsymbol{\phi}}(s' \mid s,a) \): transition kernel determined by dynamics parameters \( \boldsymbol{\phi} \),
• \( r_{\boldsymbol{\phi}}(s,a) \): reward function (may also depend on \( \boldsymbol{\phi} \) via, e.g., contact penalties),
• \( \gamma \in (0,1) \): discount factor,
• \( \rho_{\boldsymbol{\phi}} \): initial state distribution.

For a fixed policy \( \pi_{\boldsymbol{\theta}}(a\mid s) \), the return on MDP \( \mathcal{M}_{\boldsymbol{\phi}} \) is

\[ J(\boldsymbol{\theta}, \boldsymbol{\phi}) = \mathbb{E}\Bigg[ \sum_{t=0}^{\infty} \gamma^t r_{\boldsymbol{\phi}}(s_t, a_t) \,\Big|\, \mathcal{M}_{\boldsymbol{\phi}}, \pi_{\boldsymbol{\theta}} \Bigg]. \]

Standard simulation training (without domain randomization) corresponds to using a Dirac distribution \( p(\boldsymbol{\phi}) = \delta_{\boldsymbol{\phi}_0} \) at some nominal parameter \( \boldsymbol{\phi}_0 \). Domain randomization replaces this with a broader distribution that (ideally) covers the true real-world parameter \( \boldsymbol{\phi}_{\mathrm{real}} \).

Given a parameter distribution \( p(\boldsymbol{\phi}) \), the domain-randomized objective is

\[ J_{\mathrm{DR}}(\boldsymbol{\theta}) = \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[ \mathbb{E}\Bigg[ \sum_{t=0}^{\infty} \gamma^t r_{\boldsymbol{\phi}}(s_t,a_t) \,\Big|\, \mathcal{M}_{\boldsymbol{\phi}}, \pi_{\boldsymbol{\theta}} \Bigg] \right]. \]

Using the likelihood-ratio (REINFORCE) trick from Chapter 12, we can derive the policy-gradient of the domain-randomized objective. Assume the usual regularity conditions allowing interchange of gradient and integral. Then

\[ \nabla_{\boldsymbol{\theta}} J_{\mathrm{DR}}(\boldsymbol{\theta}) = \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[ \mathbb{E}\left[ \sum_{t=0}^{\infty} \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \, G_t \,\bigg|\, \mathcal{M}_{\boldsymbol{\phi}}, \pi_{\boldsymbol{\theta}} \right] \right], \]

where \( G_t \) is any return estimator (e.g., Monte Carlo return or critic-based advantage). Key point: the gradient has the same form as for a single MDP, but the expectation now also spans simulator parameters \( \boldsymbol{\phi} \). In practice, this means:

• Sample \( \boldsymbol{\phi} \) at episode reset,
• Run rollouts and compute policy-gradient estimates as usual,
• Average gradients across different sampled \( \boldsymbol{\phi} \).

Thus, any policy-gradient or actor–critic algorithm can be extended to domain randomization by adding one outer loop over randomized parameters.

We want the policy optimized under randomization to perform well on the real robot with parameters \( \boldsymbol{\phi}_{\mathrm{real}} \). Suppose the return is Lipschitz-continuous in the parameters:

\[ \forall \boldsymbol{\phi}_1, \boldsymbol{\phi}_2 \in \Phi:\quad \big| J(\boldsymbol{\theta}, \boldsymbol{\phi}_1) - J(\boldsymbol{\theta}, \boldsymbol{\phi}_2) \big| \leq L \, \|\boldsymbol{\phi}_1 - \boldsymbol{\phi}_2\|_2, \]

for some Lipschitz constant \( L \) depending on the dynamics, reward, and policy class. Then the gap between expected performance under domain randomization and performance on the real robot satisfies:

\[ \big| J(\boldsymbol{\theta}, \boldsymbol{\phi}_{\mathrm{real}}) - J_{\mathrm{DR}}(\boldsymbol{\theta}) \big| \leq L \, \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[\,\|\boldsymbol{\phi} - \boldsymbol{\phi}_{\mathrm{real}}\|_2\,\right]. \]

Proof (sketch). Using Lipschitz continuity with \( \boldsymbol{\phi}_1 = \boldsymbol{\phi}_{\mathrm{real}} \) and arbitrary \( \boldsymbol{\phi} \):

\[ \big| J(\boldsymbol{\theta}, \boldsymbol{\phi}_{\mathrm{real}}) - J(\boldsymbol{\theta}, \boldsymbol{\phi}) \big| \leq L \, \|\boldsymbol{\phi} - \boldsymbol{\phi}_{\mathrm{real}}\|_2. \]

Taking the expectation over \( \boldsymbol{\phi} \sim p(\boldsymbol{\phi}) \) on both sides and applying Jensen's inequality to the absolute value directly yields the claimed bound. The inequality suggests that:

• Covering the real parameters well (small expected distance) is crucial.
• Overly narrow \( p(\boldsymbol{\phi}) \) risks excluding \( \boldsymbol{\phi}_{\mathrm{real}} \),
• Overly broad \( p(\boldsymbol{\phi}) \) may slow learning but still provides a bound if the real parameters remain in a region with bounded Lipschitz constant.

flowchart TD
  P["Choose parameter distribution p(phi)"] --> C["Compute/estimate coverage around phi_real"]
  C --> L["Assume Lipschitz constant L for J(theta,phi)"]
  L --> B["Bound |J(theta,phi_real) - J_DR(theta)|"]
  B --> D["Adjust p(phi) (range, shape) if bound too loose"]
        

The practical strength of domain randomization depends critically on how we choose \( p(\boldsymbol{\phi}) \). Suppose \( \boldsymbol{\phi} = (m,\, \mu,\, k,\, d) \) collects:

• \( m \): link masses,
• \( \mu \): contact and friction coefficients,
• \( k \): stiffness constants (e.g., joint compliance),
• \( d \): actuator damping or delay parameters.

A simple factorized distribution is:

\[ p(\boldsymbol{\phi}) = p(m)\,p(\mu)\,p(k)\,p(d), \]

with (for example) uniform ranges identified via engineering tolerance:

\[ m \sim \mathrm{Unif}\big[(1-\alpha_m)\,m_0,\,(1+\alpha_m)\,m_0\big],\quad \mu \sim \mathrm{Unif}\big[(1-\alpha_{\mu})\,\mu_0,\,(1+\alpha_{\mu})\,\mu_0\big], \]

\[ k \sim \mathrm{LogUnif}\big[k_{\min},\,k_{\max}\big],\quad d \sim \mathcal{N}(d_0,\,\sigma_d^2) \text{ truncated to } d \ge 0. \]

For visual randomization (textures, lighting, camera pose), we can treat textures as discrete categories and continuous parameters for lighting and pose. For example, a camera elevation angle \( \theta_{\mathrm{cam}} \) might be drawn from:

\[ \theta_{\mathrm{cam}} \sim \mathrm{Unif}[\theta_{\min}, \theta_{\max}],\quad \theta_{\min} \le \theta_{\mathrm{real}} \le \theta_{\max}, \]

where \( \theta_{\mathrm{real}} \) is the physical camera configuration. Coverage of the real configuration is again essential.

We illustrate domain randomization in a simplified 2-DOF planar arm environment, assuming a Gym-style interface (as introduced in previous labs). We randomize link masses and friction at every environment reset.

import numpy as np

class PlanarArmEnvDomainRandomized:
    def __init__(self,
                 mass_nominal=(1.0, 1.0),
                 mass_rel_range=0.3,
                 friction_nominal=0.05,
                 friction_rel_range=0.5,
                 seed=None):
        self.rng = np.random.default_rng(seed)
        self.mass_nominal = np.array(mass_nominal, dtype=float)
        self.mass_rel_range = mass_rel_range
        self.friction_nominal = float(friction_nominal)
        self.friction_rel_range = friction_rel_range

        # Parameters that will be randomized each reset
        self.link_masses = self.mass_nominal.copy()
        self.joint_friction = np.array([self.friction_nominal,
                                        self.friction_nominal], dtype=float)

        # State: [q1, q2, dq1, dq2]
        self.state = np.zeros(4, dtype=float)
        self.dt = 0.02

    def _sample_parameters(self):
        # Uniform multiplicative perturbation
        low_m = (1.0 - self.mass_rel_range) * self.mass_nominal
        high_m = (1.0 + self.mass_rel_range) * self.mass_nominal
        self.link_masses = self.rng.uniform(low_m, high_m)

        low_f = (1.0 - self.friction_rel_range) * self.friction_nominal
        high_f = (1.0 + self.friction_rel_range) * self.friction_nominal
        self.joint_friction[:] = self.rng.uniform(low_f, high_f, size=2)

    def reset(self):
        self._sample_parameters()
        # Small random initial state
        self.state = self.rng.normal(0.0, 0.05, size=4)
        return self.state.copy()

    def step(self, u):
        # u: torques (2-dim)
        u = np.clip(u, -5.0, 5.0)

        # Very simplified dynamics: dqdot = (u - friction * dq) / mass
        q1, q2, dq1, dq2 = self.state
        dq = np.array([dq1, dq2])
        # element-wise "effective mass" per joint for illustration
        eff_mass = self.link_masses + 1e-3
        ddq = (u - self.joint_friction * dq) / eff_mass

        dq_new = dq + self.dt * ddq
        q_new = np.array([q1, q2]) + self.dt * dq_new

        self.state = np.concatenate([q_new, dq_new])

        # Reward: negative distance to target plus small control penalty
        target = np.array([0.0, 0.0])
        pos_error = np.linalg.norm(q_new - target)
        r = -pos_error - 0.01 * np.dot(u, u)

        done = False  # finite horizon handled outside
        info = {
            "link_masses": self.link_masses.copy(),
            "joint_friction": self.joint_friction.copy()
        }
        return self.state.copy(), float(r), done, info
      

To integrate with an RL algorithm (e.g., policy gradient), simply create this environment and call env.reset() at the beginning of each episode. Because the parameters are resampled there, trajectories will be drawn from the mixture of simulators required for \( J_{\mathrm{DR}}(\boldsymbol{\theta}) \).

In C++, domain randomization usually appears as parameter sampling code around a physics engine (e.g., Bullet, DART, MuJoCo C API). Below is a minimal skeleton (no specific engine) that randomizes link mass and contact friction at reset:

#include <random>
#include <array>

struct ArmParams {
    std::array<double, 2> link_masses;
    double contact_friction;
};

class PlanarArmEnvDR {
public:
    PlanarArmEnvDR()
    : gen_(std::random_device{}()),
      unif01_(0.0, 1.0)
    {
        mass_nominal_[0] = 1.0;
        mass_nominal_[1] = 1.0;
        friction_nominal_ = 0.5;
        mass_rel_range_ = 0.3;
        friction_rel_range_ = 0.5;
    }

    void reset() {
        sampleParameters();
        // reset simulator state, call engine API, etc.
        // engine_.setLinkMass(0, params_.link_masses[0]);
        // engine_.setLinkMass(1, params_.link_masses[1]);
        // engine_.setContactFriction(params_.contact_friction);
    }

    const ArmParams& currentParams() const { return params_; }

    // step(...) would call physics engine with params_

private:
    void sampleParameters() {
        for (int i = 0; i < 2; ++i) {
            double low = (1.0 - mass_rel_range_) * mass_nominal_[i];
            double high = (1.0 + mass_rel_range_) * mass_nominal_[i];
            params_.link_masses[i] = uniform(low, high);
        }
        double lowf = (1.0 - friction_rel_range_) * friction_nominal_;
        double highf = (1.0 + friction_rel_range_) * friction_nominal_;
        params_.contact_friction = uniform(lowf, highf);
    }

    double uniform(double a, double b) {
        return a + (b - a) * unif01_(gen_);
    }

    std::mt19937 gen_;
    std::uniform_real_distribution<double> unif01_;

    std::array<double, 2> mass_nominal_;
    double friction_nominal_;
    double mass_rel_range_;
    double friction_rel_range_;

    ArmParams params_;
};
      

The environment wrapper used by the RL algorithm calls reset() at the start of each episode, thereby sampling new dynamics parameters and pushing them into the underlying simulator.

Java-based robotics simulators (e.g., jBullet-based or custom engines) can be wrapped similarly. Below is a minimal Java skeleton illustrating the idea:

import java.util.Random;

public class PlanarArmEnvDR {
    private final Random rng;
    private final double[] massNominal = {1.0, 1.0};
    private final double massRelRange = 0.3;
    private final double frictionNominal = 0.5;
    private final double frictionRelRange = 0.5;

    private double[] linkMasses = new double[2];
    private double contactFriction;

    public PlanarArmEnvDR(long seed) {
        this.rng = new Random(seed);
    }

    private double uniform(double low, double high) {
        return low + (high - low) * rng.nextDouble();
    }

    private void sampleParameters() {
        for (int i = 0; i < 2; ++i) {
            double low = (1.0 - massRelRange) * massNominal[i];
            double high = (1.0 + massRelRange) * massNominal[i];
            linkMasses[i] = uniform(low, high);
        }
        double lowf = (1.0 - frictionRelRange) * frictionNominal;
        double highf = (1.0 + frictionRelRange) * frictionNominal;
        contactFriction = uniform(lowf, highf);
    }

    public double[] reset() {
        sampleParameters();
        // Here one would call engine APIs, e.g.:
        // engine.setLinkMass(0, linkMasses[0]);
        // engine.setLinkMass(1, linkMasses[1]);
        // engine.setContactFriction(contactFriction);

        // Return initial state (here, all zeros)
        return new double[]{0.0, 0.0, 0.0, 0.0};
    }

    public StepResult step(double[] action) {
        // Call simulator integration with current randomized params
        // ...
        double[] nextState = new double[]{0.0, 0.0, 0.0, 0.0};
        double reward = 0.0;
        boolean done = false;
        return new StepResult(nextState, reward, done);
    }

    public static class StepResult {
        public final double[] nextState;
        public final double reward;
        public final boolean done;

        public StepResult(double[] nextState, double reward, boolean done) {
            this.nextState = nextState;
            this.reward = reward;
            this.done = done;
        }
    }
}
      

As in Python and C++, the crucial point is that reset() samples new parameters from a distribution representing uncertainty about the real robot.

In MATLAB/Simulink, it is common to expose physical parameters as workspace variables used by Simulink blocks (e.g., joint mass or friction in a manipulator model). Domain randomization then becomes randomized assignment of these variables before each simulation run.

function [phi, simOut] = runEpisodeDR(seed)
    % seed: scalar for reproducibility
    rng(seed);

    % Nominal parameters
    m1_nom = 1.0; m2_nom = 1.0;
    mu_nom = 0.05;
    mass_rel_range = 0.3;
    mu_rel_range   = 0.5;

    % Sample randomized parameters
    m1 = m1_nom * (1 + mass_rel_range * (2*rand - 1));
    m2 = m2_nom * (1 + mass_rel_range * (2*rand - 1));
    mu = mu_nom * (1 + mu_rel_range   * (2*rand - 1));

    % Assign to base workspace (used by Simulink model "planar_arm_model.slx")
    assignin('base','m1',m1);
    assignin('base','m2',m2);
    assignin('base','mu_contact',mu);

    % Run simulation for one episode
    simOut = sim('planar_arm_model', 'StopTime', '5.0', ...
                 'SaveOutput', 'on', 'SaveState', 'off');

    % Collect return (assuming block outputs reward signal "r")
    r = simOut.r;  % timeseries
    dt = r.Time(2) - r.Time(1);
    gamma = 0.99;
    G = sum((gamma .^ (0:numel(r.Data)-1))' .* r.Data);

    % Output parameter vector and return for analysis
    phi = [m1; m2; mu];
    fprintf('Episode return: %f\n', G);
end
      

The Simulink model planar_arm_model.slx would use parameters m1, m2, and mu_contact in its dynamics blocks. An outer RL loop calls runEpisodeDR many times, each time with different seeds, and updates the policy using the observed returns and trajectories.

To demonstrate domain randomization in Wolfram Mathematica, consider a simple pendulum with uncertain length and mass. The equation of motion is

\[ \ddot{\theta}(t) + \frac{b}{m \ell^2} \dot{\theta}(t) + \frac{g}{\ell} \sin \theta(t) = \frac{u(t)}{m \ell^2}, \]

where \( m \), \( \ell \), and \( b \) are random parameters. In Mathematica we can integrate trajectories for many random draws and evaluate a control law \( u(t) = -k \theta(t) \).

ClearAll["Global`*"];

g = 9.81;
gamma = 0.99;
Tfinal = 5.0;

sampleEpisode[k_] := Module[
  {m, ell, b, u, sol, theta, dtheta, tgrid, r, G},
  (* Sample parameters *)
  m   = RandomReal[{0.8, 1.2}];
  ell = RandomReal[{0.8, 1.2}];
  b   = RandomReal[{0.05, 0.15}];

  (* Control law u(t) = -k * theta(t) *)
  u[theta_] := -k * theta;

  (* Dynamics: theta'' = f(theta, theta') *)
  sol = NDSolve[
    {
      theta''[t] + (b/(m ell^2)) theta'[t] + (g/ell) Sin[theta[t]] ==
        u[theta[t]]/(m ell^2),
      theta[0] == 0.3,
      theta'[0] == 0.0
    },
    {theta},
    {t, 0, Tfinal}
  ];

  tgrid = Range[0, Tfinal, 0.01];
  dtheta = theta' /. sol[[1]];
  r = - (theta[t]^2 + 0.01 dtheta[t]^2) /. sol[[1]] /. t -> # & /@ tgrid;

  G = Total[(gamma^Range[0, Length[r]-1]) r];

  <<| "m" -> m, "ell" -> ell, "b" -> b, "Return" -> G |>>
];

(* Example: evaluate expected return under domain randomization *)
With[{k = 3.0},
  Mean[Table[sampleEpisode[k]["Return"], {50}]]
]
      

An outer optimization loop can adjust the gain \( k \) to maximize the expected return across random parameter draws, thus implementing a simple form of domain randomization.

Problem 1 (Parametric MDP and Domain-Randomized Objective). Let \( \{ \mathcal{M}_{\boldsymbol{\phi}} \}_{\boldsymbol{\phi} \in \Phi} \) be a parametric MDP family, and fix a policy \( \pi_{\boldsymbol{\theta}} \). Show that training on a single nominal parameter \( \boldsymbol{\phi}_0 \) is equivalent to domain randomization with \( p(\boldsymbol{\phi}) = \delta_{\boldsymbol{\phi}_0} \).

Solution. By definition, the domain-randomized objective is

\[ J_{\mathrm{DR}}(\boldsymbol{\theta}) = \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[ J(\boldsymbol{\theta}, \boldsymbol{\phi}) \right]. \]

If \( p(\boldsymbol{\phi}) = \delta_{\boldsymbol{\phi}_0} \) (a Dirac measure at \( \boldsymbol{\phi}_0 \)), then for any integrable function \( f \) we have \( \mathbb{E}_{\boldsymbol{\phi} \sim \delta_{\boldsymbol{\phi}_0}} [f(\boldsymbol{\phi})] = f(\boldsymbol{\phi}_0) \). Applying this with \( f(\boldsymbol{\phi}) = J(\boldsymbol{\theta}, \boldsymbol{\phi}) \) yields

\[ J_{\mathrm{DR}}(\boldsymbol{\theta}) = J(\boldsymbol{\theta}, \boldsymbol{\phi}_0), \]

which is exactly the standard simulation objective with fixed parameters. Thus, classical simulation training is a degenerate special case of domain randomization.

Problem 2 (Lipschitz Bound Interpretation). Assume the Lipschitz condition from Section 4, and suppose \( \boldsymbol{\phi} \) is drawn uniformly from a ball of radius \( R \) around \( \boldsymbol{\phi}_{\mathrm{real}} \) in \( \mathbb{R}^d \). Show that \( \mathbb{E}\big[\|\boldsymbol{\phi}-\boldsymbol{\phi}_{\mathrm{real}}\|_2\big] \le R \) and interpret the bound on \( |J(\boldsymbol{\theta}, \boldsymbol{\phi}_{\mathrm{real}}) - J_{\mathrm{DR}}(\boldsymbol{\theta})| \).

Solution. If \( \boldsymbol{\phi} \) is uniformly distributed in the closed ball \( B(\boldsymbol{\phi}_{\mathrm{real}}, R) \), then by definition \( \|\boldsymbol{\phi} - \boldsymbol{\phi}_{\mathrm{real}}\|_2 \le R \) almost surely. Therefore

\[ \mathbb{E}\left[ \|\boldsymbol{\phi} - \boldsymbol{\phi}_{\mathrm{real}}\|_2 \right] \le \mathbb{E}[R] = R. \]

Plugging into the Lipschitz-based generalization bound from Section 4,

\[ \big| J(\boldsymbol{\theta}, \boldsymbol{\phi}_{\mathrm{real}}) - J_{\mathrm{DR}}(\boldsymbol{\theta}) \big| \le L \, \mathbb{E}\left[ \|\boldsymbol{\phi} - \boldsymbol{\phi}_{\mathrm{real}}\|_2 \right] \le L R. \]

Hence, reducing the radius \( R \) (i.e., making the randomization ball tighter around the real parameters) tightens the worst-case gap. However, if \( R \) is too small and the true parameters are mis-estimated, domain randomization might completely miss important regions of parameter space.

Problem 3 (One-Dimensional Linear System). Consider a scalar system \( x_{t+1} = a x_t + b u_t \) with state \( x_t \in \mathbb{R} \) and action \( u_t \in \mathbb{R} \). The parameter \( a \) is uncertain, and assumed \( a \sim \mathrm{Unif}[a_{\min}, a_{\max}] \), while \( b \) is known. We fix a linear policy \( u_t = -k x_t \). The stage cost is \( \ell(x_t, u_t) = x_t^2 + \lambda u_t^2 \), and the discounted cost is \( J(k,a) = \mathbb{E}\big[\sum_{t=0}^{\infty} \gamma^t \ell(x_t,u_t)\big] \). Assuming the closed-loop system is stable for all \( a \in [a_{\min}, a_{\max}] \), write the domain-randomized objective \( J_{\mathrm{DR}}(k) \) explicitly as an integral over \( a \).

Solution. With feedback \( u_t = -k x_t \), the closed-loop dynamics are \( x_{t+1} = (a - b k) x_t \). Let \( \alpha(a,k) = a - b k \). If the system is stable for all \( a \) in the interval, then the cost for fixed \( a \) has the standard quadratic form

\[ J(k,a) = \frac{q(a,k)}{1 - \gamma \alpha(a,k)^2}, \]

where \( q(a,k) \) is the expected stage cost at \( t=0 \) under the initial-state distribution (here we do not need its exact form to answer the question). Because \( a \sim \mathrm{Unif}[a_{\min}, a_{\max}] \), the domain-randomized objective is

\[ J_{\mathrm{DR}}(k) = \mathbb{E}_{a}[J(k,a)] = \frac{1}{a_{\max} - a_{\min}} \int_{a_{\min}}^{a_{\max}} \frac{q(a,k)}{1 - \gamma (a - b k)^2}\, \mathrm{d}a. \]

Analytic optimization of this integral may be difficult, but it illustrates the domain-randomized objective as an average over closed-loop costs for different dynamics parameters.

Problem 4 (Gradient Estimator with Domain Randomization). Suppose we implement a Monte Carlo policy-gradient algorithm with domain randomization. At each iteration, we sample \( \boldsymbol{\phi}^{(i)} \sim p(\boldsymbol{\phi}) \), roll out trajectory \( \tau^{(i)} \), and compute an episodic gradient estimate

\[ \hat{g}^{(i)} = \sum_{t=0}^{T^{(i)}-1} \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t^{(i)} \mid s_t^{(i)}) \, \hat{G}_t^{(i)}. \]

Show that \( \mathbb{E}[\hat{g}^{(i)}] = \nabla_{\boldsymbol{\theta}} J_{\mathrm{DR}}(\boldsymbol{\theta}) \), assuming the trajectories are sampled from the parametric MDPs under \( p(\boldsymbol{\phi}) \).

Solution. Conditioned on a sampled parameter \( \boldsymbol{\phi}^{(i)} \), the standard policy-gradient theorem (Chapter 12) implies that

\[ \mathbb{E}\left[ \hat{g}^{(i)} \,\big|\, \boldsymbol{\phi}^{(i)} \right] = \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}, \boldsymbol{\phi}^{(i)}). \]

Taking expectation over \( \boldsymbol{\phi}^{(i)} \sim p(\boldsymbol{\phi}) \), and applying the tower property of conditional expectation,

\[ \mathbb{E}[\hat{g}^{(i)}] = \mathbb{E}_{\boldsymbol{\phi}^{(i)}}\left[ \mathbb{E}\left[ \hat{g}^{(i)} \mid \boldsymbol{\phi}^{(i)} \right] \right] = \mathbb{E}_{\boldsymbol{\phi} \sim p(\boldsymbol{\phi})} \left[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}, \boldsymbol{\phi}) \right] = \nabla_{\boldsymbol{\theta}} J_{\mathrm{DR}}(\boldsymbol{\theta}), \]

where we used the fact that gradient and expectation can be interchanged under standard regularity assumptions. Thus, the Monte Carlo gradient estimator is unbiased for the domain-randomized objective.

Problem 5 (Qualitative Design Flow for Domain Randomization). Sketch a decision flow for choosing which parameters to randomize and the ranges for randomization in a manipulator simulation. Consider uncertainty in masses, friction, sensor noise, and camera pose.

Solution (flow sketch).

flowchart TD
  ST["Start"] --> U["List uncertain quantities (mass, friction, noise, camera pose)"]
  U --> P1["Estimate nominal values from CAD, datasheets"]
  P1 --> R1["Decide tolerances from manufacturing + experiments"]
  R1 --> C1["Is parameter critical to task performance?"]
  C1 -->|yes| DR1["Include in domain randomization with wider range"]
  C1 -->|no| DR2["Include with narrow range or fix"]
  DR1 --> F["Check that real robot values lie inside ranges"]
  DR2 --> F
  F --> E["Evaluate coverage vs training stability; adjust if needed"]
        

This process ensures that randomization covers realistic variability without unnecessarily enlarging the parameter space, which would slow optimization and risk under-training in relevant regions.

In this lesson we formalized domain randomization as training over a distribution of simulators modeled as a parametric MDP family. We defined the domain-randomized objective \( J_{\mathrm{DR}}(\boldsymbol{\theta}) \), showed that standard policy-gradient methods extend directly when we add an outer sampling loop over parameters, and used a Lipschitz continuity argument to bound the performance gap between randomized training and the real robot. We then discussed the design of parameter distributions and provided implementation skeletons in Python, C++, Java, MATLAB/Simulink, and Mathematica. These ideas underpin more systematic sim calibration and residual-learning methods in subsequent lessons.