The lab centers on a 2-DOF planar manipulator performing a reaching or pick-and-place task in simulation. Students will:

• Specify the manipulation task as a continuous-state, continuous-action MDP with state \( s_t \), action \( a_t \), and reward \( r_t \).
• Implement or configure a simulator providing the stochastic transition \( s_{t+1} \sim p(\cdot \mid s_t,a_t) \).
• Train an off-policy actor–critic policy (e.g., SAC / TD3 as introduced in earlier lessons).
• Perform transfer tests by modifying dynamics parameters and evaluating the learned policy.

Conceptually, this implements the following pipeline:

flowchart TD
  A["Define manipulation task as MDP"] --> B["Implement planar manipulator simulator"]
  B --> C["Configure actor-critic RL algorithm (e.g. SAC)"]
  C --> D["Collect rollouts and train in simulation"]
  D --> E["Evaluate in nominal simulator"]
  E --> F["Perturb dynamics (mass, friction, delay)"]
  F --> G["Transfer test and measure performance drop"]
        

Throughout, we relate each programming exercise (Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica) to the same underlying mathematical formulation.

We model the manipulation problem as a continuous-state MDP \( \mathcal{M} = (\mathcal{S}, \mathcal{A}, p, r, \gamma) \). For an \( n \)-DOF planar arm, let joint coordinates be \( q_t \in \mathbb{R}^n \) and joint velocities \( \dot{q}_t \in \mathbb{R}^n \). We consider torque control \( \tau_t \in \mathbb{R}^n \).

A practical state design for a reaching task is

\[ s_t = \begin{bmatrix} q_t \\ \dot{q}_t \\ x_{\mathrm{ee}}(q_t) - x_{\mathrm{goal}} \end{bmatrix} \in \mathbb{R}^{2n + d_x}, \]

where \( x_{\mathrm{ee}}(q_t) \) is the end-effector position from forward kinematics and \( x_{\mathrm{goal}} \) is the goal pose (fixed during one episode).

We choose actions as joint torques: \( a_t = \tau_t \in \mathbb{R}^n \). The (discretized) manipulator dynamics can be abstracted as

\[ s_{t+1} = f(s_t, a_t) + w_t, \quad w_t \sim \mathcal{N}(0, \Sigma_w), \]

where the function \( f \) encodes the rigid-body dynamics (derived from the Euler–Lagrange equations, which you know from kinematics/dynamics) and the state update via numerical integration.

A common quadratic reward shaping for reaching is

\[ r(s_t, a_t) = - \big\| x_{\mathrm{ee}}(q_t) - x_{\mathrm{goal}} \big\|_2^2 - \lambda_{\tau} \, \| \tau_t \|_2^2, \]

where \( \lambda_{\tau} > 0 \) penalizes large torques. The RL objective, for a stochastic policy \( \pi_{\theta}(a \mid s) \) with parameters \( \theta \), is

\[ J(\theta) = \mathbb{E}_{\tau \sim \pi_{\theta}} \left[ \sum_{t=0}^{T-1} \gamma^t r(s_t, a_t) \right], \]

where the expectation is over trajectories \( \tau = (s_0, a_0, s_1, a_1, \dots) \) induced by the simulator and policy. The lab focuses on estimating and maximizing \( J(\theta) \) in simulation and then measuring the same quantity under perturbed dynamics.

We now recall the soft actor–critic (SAC) style objective instantiated for our manipulation MDP. Let \( Q_{\phi}(s,a) \) denote a critic network with parameters \( \phi \) approximating the soft state–action value \( Q^{\pi} \). The replay buffer \( \mathcal{D} \) stores transition tuples \( (s,a,r,s') \).

The soft Bellman target for continuous control is

\[ y(r,s') = r + \gamma \, \mathbb{E}_{a' \sim \pi_{\theta}(\cdot \mid s')} \big[ Q_{\bar{\phi}}(s', a') - \alpha \log \pi_{\theta}(a' \mid s') \big], \]

where \( Q_{\bar{\phi}} \) uses slowly updated target parameters and \( \alpha > 0 \) is the entropy temperature.

The critic loss is the mean-squared Bellman error:

\[ J_Q(\phi) = \mathbb{E}_{(s,a,r,s') \sim \mathcal{D}} \Big[ \big( Q_{\phi}(s,a) - y(r,s') \big)^2 \Big]. \]

The policy (actor) minimizes the soft value:

\[ J_{\pi}(\theta) = \mathbb{E}_{s \sim \mathcal{D}} \left[ \mathbb{E}_{a \sim \pi_{\theta}(\cdot \mid s)} \big[ \alpha \log \pi_{\theta}(a \mid s) - Q_{\phi}(s,a) \big] \right]. \]

In practice, we estimate these expectations via minibatches from \( \mathcal{D} \), using the simulator to fill the buffer. The lab implementations in all languages are structured so that:

• The simulator produces samples \( (s,a,r,s') \).
• The critic update minimizes \( J_Q(\phi) \).
• The actor update minimizes \( J_{\pi}(\theta) \).

For comparison, the deterministic policy gradient (DPG) for a deterministic policy \( \mu_{\theta}(s) \) is

\[ \nabla_{\theta} J(\theta) = \mathbb{E}_{s \sim \rho^{\mu}} \left[ \nabla_{\theta} \mu_{\theta}(s) \, \nabla_{a} Q^{\mu}(s,a) \big|_{a = \mu_{\theta}(s)} \right], \]

which you derived in previous lessons; some of the lab code (e.g., in C++ and Mathematica) uses this simpler DPG objective for clarity.

Let \( d \) denote a vector of dynamics parameters (link masses, inertias, friction coefficients, motor gains). The simulator implements

\[ s_{t+1} = f(s_t, a_t; d) + w_t. \]

During training we fix or randomize \( d \) according to a training distribution \( p_{ \mathrm{train} }(d) \) and optimize

\[ J_{ \mathrm{sim} }(\theta) = \mathbb{E}_{ d \sim p_{ \mathrm{train} } } \left[ J(\theta; d) \right]. \]

For transfer, we evaluate the same policy under a different parameter distribution \( p_{\mathrm{test}}(d) \) (e.g., shifted masses or friction), yielding

\[ J_{\mathrm{test}}(\theta) = \mathbb{E}_{d \sim p_{\mathrm{test}}} \left[ J(\theta; d) \right]. \]

A key quantitative transfer metric is the transfer gap

\[ \Delta J(\theta) = J_{\mathrm{test}}(\theta) - J_{\mathrm{sim}}(\theta), \]

estimated empirically by averaging returns over multiple episodes in each setting. This lab focuses on measuring and interpreting \( \Delta J(\theta) \) rather than minimizing it (which is studied in the next chapter on sim-to-real).

flowchart TD
  S1["Train with dynamics distribution p_train(d)"] --> P1["Learn policy parameters theta*"]
  P1 --> E1["Evaluate in nominal simulator \n(sample d ~ p_train)"]
  P1 --> E2["Evaluate in perturbed simulator \n(sample d ~ p_test)"]
  E1 --> M1["Estimate J_sim(theta*)"]
  E2 --> M2["Estimate J_test(theta*)"]
  M1 --> G["Compute transfer gap Delta J"]
  M2 --> G
        

In Python we use a planar 2-DOF manipulator simulated in PyBullet and a SAC implementation (e.g., stable-baselines3). The environment follows the Gymnasium API.

5.1 Environment Definition

import gymnasium as gym
import numpy as np
import pybullet as p
import pybullet_data

class PlanarReachEnv(gym.Env):
    metadata = {"render_modes": ["human"], "render_fps": 60}

    def __init__(self, render=False, dynamics_params=None):
        super().__init__()
        self.render = render
        self.dt = 1.0 / 60.0
        self.max_steps = 200
        self.step_count = 0

        # Dynamics parameters (mass, friction, etc.)
        default_params = {
            "link_masses": [1.0, 1.0],
            "joint_damping": [0.05, 0.05],
            "joint_max_torque": 5.0,
        }
        if dynamics_params is None:
            dynamics_params = default_params
        self.params = dynamics_params

        # Observation: [q1, q2, dq1, dq2, ee_x - goal_x, ee_y - goal_y]
        high = np.array([np.pi, np.pi, 5.0, 5.0, 2.0, 2.0], dtype=np.float32)
        self.observation_space = gym.spaces.Box(-high, high, dtype=np.float32)

        # Actions: joint torques
        max_torque = self.params["joint_max_torque"]
        self.action_space = gym.spaces.Box(
            low=-max_torque, high=max_torque, shape=(2,), dtype=np.float32
        )

        self._client = None
        self.robot_id = None
        self.goal = np.array([0.4, 0.4], dtype=np.float32)

    def _connect(self):
        if self._client is None:
            self._client = p.connect(p.GUI if self.render else p.DIRECT)
            p.setAdditionalSearchPath(pybullet_data.getDataPath())
            p.setGravity(0, 0, -9.81, physicsClientId=self._client)

    def _reset_sim(self):
        p.resetSimulation(physicsClientId=self._client)
        p.setGravity(0, 0, -9.81, physicsClientId=self._client)
        plane_id = p.loadURDF("plane.urdf")
        # Simple 2-link arm URDF with revolute joints
        self.robot_id = p.loadURDF("two_link_manip.urdf", basePosition=[0, 0, 0])

        # Set dynamics parameters
        for link_idx, m in enumerate(self.params["link_masses"]):
            p.changeDynamics(
                self.robot_id, link_idx,
                mass=float(m),
                linearDamping=self.params["joint_damping"][link_idx],
                angularDamping=self.params["joint_damping"][link_idx],
            )

    def _get_state(self):
        joint_states = [
            p.getJointState(self.robot_id, j, physicsClientId=self._client)
            for j in range(2)
        ]
        q = np.array([js[0] for js in joint_states], dtype=np.float32)
        dq = np.array([js[1] for js in joint_states], dtype=np.float32)

        # EE position via forward kinematics from PyBullet
        ee_state = p.getLinkState(self.robot_id, 2, computeForwardKinematics=True)
        ee_pos = np.array(ee_state[0][0:2], dtype=np.float32)  # (x, y)
        diff = ee_pos - self.goal
        return np.concatenate([q, dq, diff], axis=0)

    def compute_reward(self, state, action):
        pos_error = state[-2:]
        torque_penalty = 0.001 * np.sum(action * action)
        return -float(pos_error.dot(pos_error) + torque_penalty)

    def step(self, action):
        self.step_count += 1
        # Clip action to torque limits
        action = np.clip(action, self.action_space.low, self.action_space.high)

        p.setJointMotorControlArray(
            self.robot_id,
            jointIndices=[0, 1],
            controlMode=p.TORQUE_CONTROL,
            forces=action,
            physicsClientId=self._client,
        )
        p.stepSimulation(physicsClientId=self._client)

        obs = self._get_state()
        reward = self.compute_reward(obs, action)
        terminated = False
        if np.linalg.norm(obs[-2:]) < 0.02:
            terminated = True
        truncated = self.step_count >= self.max_steps
        info = {}
        return obs, reward, terminated, truncated, info

    def reset(self, seed=None, options=None):
        super().reset(seed=seed)
        self._connect()
        self._reset_sim()
        self.step_count = 0

        # Randomize initial joint angles
        for j in range(2):
            init_q = self.np_random.uniform(low=-np.pi / 2.0, high=np.pi / 2.0)
            p.resetJointState(self.robot_id, j, init_q, 0.0, physicsClientId=self._client)

        obs = self._get_state()
        info = {}
        return obs, info

    def close(self):
        if self._client is not None:
            p.disconnect(physicsClientId=self._client)
            self._client = None
      

5.2 Training and Transfer Test with SAC

from stable_baselines3 import SAC

# Nominal dynamics (train)
train_env = PlanarReachEnv(render=False)

# Train SAC in simulation
model = SAC(
    policy="MlpPolicy",
    env=train_env,
    learning_rate=3e-4,
    buffer_size=200000,
    batch_size=256,
    gamma=0.99,
    tau=0.005,  # target smoothing
    train_freq=1,
    gradient_steps=1,
    policy_kwargs=dict(net_arch=[256, 256]),
    verbose=1,
)
model.learn(total_timesteps=200000)

# Transfer: perturbed dynamics (heavier links, more friction)
test_params = {
    "link_masses": [1.3, 1.3],
    "joint_damping": [0.08, 0.08],
    "joint_max_torque": 5.0,
}
test_env = PlanarReachEnv(render=False, dynamics_params=test_params)

def evaluate(env, model, n_episodes=20):
    returns = []
    for ep in range(n_episodes):
        obs, info = env.reset()
        done = False
        truncated = False
        ep_ret = 0.0
        while not (done or truncated):
            action, _ = model.predict(obs, deterministic=True)
            obs, r, done, truncated, info = env.step(action)
            ep_ret += r
        returns.append(ep_ret)
    return np.mean(returns), np.std(returns)

J_sim_mean, J_sim_std = evaluate(train_env, model)
J_test_mean, J_test_std = evaluate(test_env, model)
print("J_sim:", J_sim_mean, "±", J_sim_std)
print("J_test:", J_test_mean, "±", J_test_std)
print("Delta J:", J_test_mean - J_sim_mean)
      

This Python lab should be run first; the subsequent C++, Java, MATLAB, and Mathematica snippets mirror the same conceptual structure but with lighter-weight implementations.

In C++ we implement a simplified deterministic policy gradient loop using a linear policy and a quadratic critic. Dynamics can be computed via your preferred rigid-body library (e.g., Pinocchio or RBDL), but here we abstract them into a function stepDynamics.

#include <Eigen/Dense>
#include <vector>
#include <random>

using Eigen::VectorXd;
using Eigen::MatrixXd;

struct Transition {
    VectorXd s;
    VectorXd a;
    double r;
    VectorXd s_next;
    bool done;
};

class PlanarManipDPG {
public:
    PlanarManipDPG(int state_dim, int action_dim)
        : state_dim_(state_dim),
          action_dim_(action_dim),
          W_(MatrixXd::Zero(action_dim, state_dim)),
          critic_P_(MatrixXd::Identity(state_dim, state_dim)),
          gamma_(0.99),
          alpha_actor_(1e-3),
          alpha_critic_(1e-3) {}

    VectorXd policy(const VectorXd& s) const {
        // Deterministic linear policy: a = W s
        return W_ * s;
    }

    double critic(const VectorXd& s) const {
        // Quadratic critic V(s) = s' P s
        return 0.5 * s.transpose() * critic_P_ * s;
    }

    void update(const std::vector<Transition>& batch) {
        // Critic update: TD(0) on V(s) ~= s' P s
        for (const auto& tr : batch) {
            double v_s = critic(tr.s);
            double v_next = critic(tr.s_next);
            double td_target = tr.r + (tr.done ? 0.0 : gamma_ * v_next);
            double td_error = td_target - v_s;

            // Gradient wrt P: dV/dP = 0.5 * (s s' + s s')
            MatrixXd grad_P = -td_error * (tr.s * tr.s.transpose());
            critic_P_ -= alpha_critic_ * grad_P;
        }

        // Actor update: DPG-style using critic gradient approximated by P
        for (const auto& tr : batch) {
            VectorXd a = policy(tr.s);
            // Here we fake Q gradient as proportional to -a and -s for illustration.
            // In practice, use a separate Q(s,a) critic.
            VectorXd grad_mu = a;  // d(mu)/dW * dQ/da; here simplified
            MatrixXd grad_W = grad_mu * tr.s.transpose();
            W_ += alpha_actor_ * grad_W;
        }
    }

private:
    int state_dim_;
    int action_dim_;
    MatrixXd W_;
    MatrixXd critic_P_;
    double gamma_;
    double alpha_actor_;
    double alpha_critic_;
};

// Pseudo-code for rollout and training loop
int main() {
    const int state_dim = 6;
    const int action_dim = 2;
    PlanarManipDPG agent(state_dim, action_dim);

    std::mt19937 rng(0);
    for (int episode = 0; episode < 500; ++episode) {
        VectorXd s = VectorXd::Zero(state_dim); // reset from simulator
        bool done = false;
        std::vector<Transition> traj;
        while (!done) {
            VectorXd a = agent.policy(s);
            // simulate one step: stepDynamics(s, a, s_next, r, done)
            VectorXd s_next(state_dim);
            double r = 0.0;
            stepDynamics(s, a, s_next, r, done); // user-implemented dynamics

            traj.push_back({s, a, r, s_next, done});
            s = s_next;
        }
        agent.update(traj);
    }
    return 0;
}
      

Although the update rules here are simplified, the structure matches the theoretical DPG formula and can be extended by plugging a genuine critic \( Q_{\phi}(s,a) \) using, for example, a small neural network library.

Java has fewer mainstream robotics libraries, but we can use jBullet for dynamics and ND4J (or DL4J) for linear algebra. Below is a skeleton of a policy-gradient loop.

import org.nd4j.linalg.api.ndarray.INDArray;
import org.nd4j.linalg.factory.Nd4j;
import java.util.ArrayList;
import java.util.List;
import java.util.Random;

class Transition {
    INDArray s;
    INDArray a;
    double r;
    INDArray sNext;
    boolean done;
}

public class PlanarManipPG {
    private final int stateDim;
    private final int actionDim;
    private INDArray W;      // actionDim x stateDim
    private double alphaActor = 1e-3;
    private double gamma = 0.99;
    private Random rng = new Random(0);

    public PlanarManipPG(int stateDim, int actionDim) {
        this.stateDim = stateDim;
        this.actionDim = actionDim;
        this.W = Nd4j.zeros(actionDim, stateDim);
    }

    public INDArray policy(INDArray s) {
        // a = W s
        return W.mmul(s);
    }

    public void update(List<Transition> batch) {
        // Simple REINFORCE-like update with baseline (omitted here)
        for (Transition tr : batch) {
            // Compute Monte Carlo return G_t
            double G = 0.0;
            double pow = 1.0;
            for (Transition tr2 : batch) {
                G += pow * tr2.r;
                pow *= gamma;
            }
            INDArray a = policy(tr.s);
            // Here we approximate grad log pi by a - mu(s)
            INDArray gradLogPi = a.neg(); // placeholder
            INDArray gradW = gradLogPi.mmul(tr.s.transpose()).mul(G);
            W.addi(gradW.mul(alphaActor));
        }
    }

    public static void main(String[] args) {
        int stateDim = 6;
        int actionDim = 2;
        PlanarManipPG agent = new PlanarManipPG(stateDim, actionDim);

        for (int episode = 0; episode < 500; episode++) {
            INDArray s = Nd4j.zeros(stateDim, 1); // get initial state from simulator
            boolean done = false;
            List<Transition> traj = new ArrayList<>();
            while (!done) {
                INDArray a = agent.policy(s);
                // Step jBullet-based manipulator (user-implemented)
                StepResult res = stepDynamics(s, a); // returns (sNext, r, done)
                Transition tr = new Transition();
                tr.s = s;
                tr.a = a;
                tr.r = res.reward;
                tr.sNext = res.nextState;
                tr.done = res.done;
                traj.add(tr);
                s = res.nextState;
            }
            agent.update(traj);
        }
    }
}
      

This code mirrors the Python SAC pipeline conceptually (interaction with a simulator, collection of trajectories, gradient-based updates) but keeps the policy and critic extremely simple to highlight the data flow rather than the network details.

MATLAB's Reinforcement Learning Toolbox and Robotics System Toolbox provide a convenient way to couple Simulink manipulator models with RL agents. Assume you have built a Simulink model TwoLinkRLModel with:

• Block inputs: joint torques tau (dimension 2).
• Block outputs: state vector obs and scalar reward.
• A stop condition when the episode terminates.

% Observation and action specifications
obsDim = 6; % [q1, q2, dq1, dq2, ee_x - goal_x, ee_y - goal_y]
actDim = 2;

obsInfo = rlNumericSpec([obsDim 1], ...
    'LowerLimit', -ones(obsDim,1)*inf, ...
    'UpperLimit',  ones(obsDim,1)*inf);
obsInfo.Name = 'observations';

actInfo = rlNumericSpec([actDim 1], ...
    'LowerLimit', -5*ones(actDim,1), ...
    'UpperLimit',  5*ones(actDim,1));
actInfo.Name = 'torques';

% Create Simulink environment
mdl = 'TwoLinkRLModel';
agentBlk = [mdl '/RL Agent'];
env = rlSimulinkEnv(mdl, agentBlk, obsInfo, actInfo);

% Actor and critic networks (SAC-style)
statePath = [
    featureInputLayer(obsDim, 'Normalization', 'none', 'Name', 'obs')
    fullyConnectedLayer(256, 'Name', 'fc1')
    reluLayer('Name', 'relu1')
    fullyConnectedLayer(256, 'Name', 'fc2')
    reluLayer('Name', 'relu2')];

criticHead = [
    fullyConnectedLayer(1, 'Name', 'value')];

criticNet = layerGraph(statePath);
criticNet = addLayers(criticNet, criticHead);
criticNet = connectLayers(criticNet, 'relu2', 'value');

actorNet = layerGraph([
    featureInputLayer(obsDim, 'Normalization', 'none', 'Name', 'obs')
    fullyConnectedLayer(256, 'Name', 'afc1')
    reluLayer('Name', 'arelu1')
    fullyConnectedLayer(256, 'Name', 'afc2')
    reluLayer('Name', 'arelu2')
    fullyConnectedLayer(actDim, 'Name', 'mean')]);

criticOpts = rlOptimizerOptions('LearnRate', 1e-3, 'GradientThreshold', 1.0);
actorOpts  = rlOptimizerOptions('LearnRate', 1e-3, 'GradientThreshold', 1.0);

critic = rlValueFunction(criticNet, obsInfo, ...
    'ObservationInputNames', 'obs', ...
    'OptimizerOptions', criticOpts);

actor = rlContinuousGaussianActor(actorNet, obsInfo, actInfo, ...
    'ObservationInputNames', 'obs', ...
    'ActionMeanOutputNames', 'mean', ...
    'OptimizerOptions', actorOpts);

agentOpts = rlSACAgentOptions( ...
    'SampleTime', 0.02, ...
    'ExperienceBufferLength', 1e6, ...
    'DiscountFactor', 0.99);

agent = rlSACAgent(actor, critic, agentOpts);

% Train in simulation
trainOpts = rlTrainingOptions( ...
    'MaxEpisodes', 500, ...
    'MaxStepsPerEpisode', 200, ...
    'ScoreAveragingWindowLength', 20, ...
    'StopTrainingCriteria', 'AverageReward', ...
    'StopTrainingValue', -1.0, ...
    'Verbose', true);

trainingStats = train(agent, env);

% Transfer test: change link masses in the Simulink model
set_param([mdl '/TwoLinkDynamics'], 'LinkMassScale', '1.3'); % e.g. mask parameter
simOptions = rlSimulationOptions('MaxSteps', 200);
simResults = sim(env, agent, simOptions);
      

This MATLAB/Simulink setup closely matches the Python SAC formulation but delegates the system dynamics to a Simulink block diagram, which is often preferred in control engineering workflows.

Mathematica is well-suited to verify analytic expressions for gradients. Consider a 1-DOF simplified manipulator with state \( s_t = (q_t, \dot{q}_t) \), torque action \( a_t \), and linear policy \( a_t = \theta_1 q_t + \theta_2 \dot{q}_t \). We model one-step dynamics as

\[ q_{t+1} = q_t + \Delta t \, \dot{q}_t, \quad \dot{q}_{t+1} = \dot{q}_t + \Delta t \, a_t. \]

The reward is

\[ r_t = - q_t^2 - \lambda a_t^2. \]

The code below symbolically computes the gradient of the one-step return with respect to \( \theta_1, \theta_2 \) and validates the DPG structure.

ClearAll["Global`*"];

(* Parameters *)
dt = Symbol["dt"];
lambda = Symbol["lambda"];

(* State and policy parameters *)
q = Symbol["q"];
dq = Symbol["dq"];
theta1 = Symbol["theta1"];
theta2 = Symbol["theta2"];

(* Deterministic policy a = theta . s *)
a[q_, dq_] := theta1*q + theta2*dq;

(* One-step dynamics *)
qNext[q_, dq_] := q + dt*dq;
dqNext[q_, dq_] := dq + dt*a[q, dq];

(* Reward *)
r[q_, dq_] := -q^2 - lambda*a[q, dq]^2;

(* One-step return J(theta) = r(s0, a0) + gamma * V(s1) with V(s1) ~ -q1^2 *)
gamma = Symbol["gamma"];
V[q_, dq_] := -q^2;

J[theta1_, theta2_] := Module[
  {q0, dq0, q1, dq1},
  q0 = q; dq0 = dq;
  q1 = qNext[q0, dq0];
  dq1 = dqNext[q0, dq0];
  r[q0, dq0] + gamma*V[q1, dq1]
];

gradTheta1 = D[J[theta1, theta2], theta1] // Simplify;
gradTheta2 = D[J[theta1, theta2], theta2] // Simplify;

Print["dJ/dtheta1 = ", gradTheta1];
Print["dJ/dtheta2 = ", gradTheta2];

(* Numerical gradient check *)
numGrad[q0_, dq0_, t1_, t2_, dtVal_, lambdaVal_, gammaVal_] := Module[
  {f, g1, g2},
  f[t1p_, t2p_] := 
    J[t1p, t2p] /. {
      q -> q0, dq -> dq0, dt -> dtVal,
      lambda -> lambdaVal, gamma -> gammaVal
    };
  g1 = (f[t1 + 1.*^-5, t2] - f[t1 - 1.*^-5, t2])/(2.*^-5);
  g2 = (f[t1, t2 + 1.*^-5] - f[t1, t2 - 1.*^-5])/(2.*^-5);
  {g1, g2}
];

(* Example numeric check *)
gradSym = {
  gradTheta1 /. {q -> 0.5, dq -> 0.1, dt -> 0.05, lambda -> 0.01, gamma -> 0.99,
    theta1 -> 0.2, theta2 -> -0.1},
  gradTheta2 /. {q -> 0.5, dq -> 0.1, dt -> 0.05, lambda -> 0.01, gamma -> 0.99,
    theta1 -> 0.2, theta2 -> -0.1}
};

gradNum = numGrad[0.5, 0.1, 0.2, -0.1, 0.05, 0.01, 0.99];
Print["Symbolic gradient: ", gradSym];
Print["Numeric gradient:   ", gradNum];
      

The resulting symbolic expressions confirm that \( \nabla_{\theta} J(\theta) \) is proportional to the state features times a term related to \( \nabla_a Q^{\mu}(s,a) \), as predicted by the deterministic policy gradient theorem.

Problem 1 (Soft Bellman Operator Contraction): Consider the soft Bellman operator used in SAC, acting on a bounded function \( Q : \mathcal{S}\times\mathcal{A} \to \mathbb{R} \):

\[ (\mathcal{T}^{\pi} Q)(s,a) = r(s,a) + \gamma \, \mathbb{E}_{s' \sim p(\cdot \mid s,a),\, a' \sim \pi(\cdot \mid s')} \big[ Q(s',a') - \alpha \log \pi(a' \mid s') \big]. \]

Show that \( \mathcal{T}^{\pi} \) is a \( \gamma \)-contraction in the sup norm and hence has a unique fixed point \( Q^{\pi} \).

Solution:

For two bounded functions \( Q_1, Q_2 \) we have

\[ \big| (\mathcal{T}^{\pi} Q_1)(s,a) - (\mathcal{T}^{\pi} Q_2)(s,a) \big| = \gamma \left| \mathbb{E} \big[ Q_1(s',a') - Q_2(s',a') \big] \right| \leq \gamma \, \| Q_1 - Q_2 \|_{\infty}. \]

Taking the supremum over \( (s,a) \) yields \( \| \mathcal{T}^{\pi} Q_1 - \mathcal{T}^{\pi} Q_2 \|_{\infty} \leq \gamma \, \| Q_1 - Q_2 \|_{\infty} \), so \( \mathcal{T}^{\pi} \) is a contraction and possesses a unique fixed point by the Banach fixed-point theorem.

Problem 2 (Deterministic Policy Gradient Derivation): For a deterministic policy \( \mu_{\theta} : \mathcal{S} \to \mathcal{A} \), show that

\[ \nabla_{\theta} J(\theta) = \mathbb{E}_{s \sim \rho^{\mu}} \left[ \nabla_{\theta} \mu_{\theta}(s) \, \nabla_a Q^{\mu}(s,a) \big|_{a = \mu_{\theta}(s)} \right], \]

where \( \rho^{\mu} \) is the discounted state visitation distribution.

Solution:

Starting from \( J(\theta) = \mathbb{E}_{s \sim \rho^{\mu}}[V^{\mu}(s)] \), we use the chain rule and the relation \( V^{\mu}(s) = Q^{\mu}(s, \mu_{\theta}(s)) \). Then

\[ \nabla_{\theta} J(\theta) = \int_{\mathcal{S}} \rho^{\mu}(s) \, \nabla_{\theta} V^{\mu}(s) \, \mathrm{d}s = \int_{\mathcal{S}} \rho^{\mu}(s) \, \nabla_{\theta} Q^{\mu}\big(s, \mu_{\theta}(s)\big) \,\mathrm{d}s. \]

Applying the chain rule inside the integral gives

\[ \nabla_{\theta} Q^{\mu}\big(s, \mu_{\theta}(s)\big) = \nabla_a Q^{\mu}(s,a)\big|_{a = \mu_{\theta}(s)} \, \nabla_{\theta} \mu_{\theta}(s), \]

which yields the desired expression after substituting back into the integral and interpreting it as an expectation over \( s \sim \rho^{\mu} \).

Problem 3 (Gradient under Domain Randomization): Suppose dynamics parameters \( d \) are sampled from a training distribution \( p_{\mathrm{train}}(d) \). Show that

\[ \nabla_{\theta} J_{\mathrm{sim}}(\theta) = \mathbb{E}_{d \sim p_{\mathrm{train}}} \big[ \nabla_{\theta} J(\theta; d) \big], \]

where \( J_{\mathrm{sim}}(\theta) \) is defined in Section 4.

Solution:

By definition, \( J_{\mathrm{sim}}(\theta) = \mathbb{E}_{d}[J(\theta; d)] \). Assuming we can interchange integration and differentiation (dominated convergence), we obtain

\[ \nabla_{\theta} J_{\mathrm{sim}}(\theta) = \nabla_{\theta} \int J(\theta; d) \, p_{\mathrm{train}}(d) \,\mathrm{d}d = \int \nabla_{\theta} J(\theta; d) \, p_{\mathrm{train}}(d) \,\mathrm{d}d = \mathbb{E}_{d} \big[ \nabla_{\theta} J(\theta; d) \big]. \]

Problem 4 (Quadratic Cost Manipulator, One-Step Optimal Torque): For the 1-DOF manipulator from Section 9, consider a single time step with state \( s_t = (q_t, \dot{q}_t) \) and torque \( a_t \). Ignore future rewards and find the torque \( a_t^{\star} \) minimizing

\[ r_t = - q_t^2 - \lambda a_t^2. \]

Solution:

Maximizing \( r_t \) is equivalent to minimizing \( q_t^2 + \lambda a_t^2 \) with respect to \( a_t \). Since \( q_t \) does not depend on \( a_t \), the minimum is obtained by minimizing \( \lambda a_t^2 \), which occurs at \( a_t^{\star} = 0 \). This is consistent with the intuition that with only instantaneous quadratic penalties and no future reward, the optimal action is to apply no torque.

Problem 5 (Sample Complexity for Transfer Evaluation): Let \( R_1, \dots, R_N \) be i.i.d. episode returns of a fixed policy under perturbed dynamics (transfer setting), with mean \( \mu \) and range bounded in \( [R_{\min}, R_{\max}] \). How many episodes \( N \) are sufficient to guarantee that the empirical mean \( \hat{\mu} \) satisfies \( |\hat{\mu} - \mu| \leq \varepsilon \) with probability at least \( 1 - \delta \)?

Solution:

By Hoeffding's inequality for bounded random variables,

\[ \mathbb{P} \big( |\hat{\mu} - \mu| \geq \varepsilon \big) \leq 2 \exp \left( - \frac{2 N \varepsilon^2}{(R_{\max} - R_{\min})^2} \right). \]

To ensure this probability is at most \( \delta \) we require

\[ 2 \exp \left( - \frac{2 N \varepsilon^2}{(R_{\max} - R_{\min})^2} \right) \leq \delta, \]

which, after rearranging, yields

\[ N \geq \frac{(R_{\max} - R_{\min})^2}{2 \varepsilon^2} \, \log \frac{2}{\delta}. \]

Thus this many transfer episodes suffice to estimate the mean return within \( \varepsilon \) with confidence \( 1 - \delta \).

In this lab you instantiated the abstract RL concepts from previous lessons on a concrete robotic manipulation task: defining a continuous MDP for a planar arm, specifying actor–critic objectives (including the soft Bellman operator), and implementing full training pipelines in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica. You also carried out transfer tests by perturbing dynamics parameters, quantified the resulting performance gap, and connected these empirical results to theoretical properties such as Bellman contractions, deterministic policy gradients, and the behavior of gradients under domain randomization.

The next chapter will deepen the sim-to-real perspective, introducing systematic domain randomization and calibration techniques designed to reduce the observed transfer gap \( \Delta J(\theta) \).