In many robotic tasks, the action space is continuous (e.g., joint torques, desired joint velocities, or end-effector twists). Value-based methods (like tabular Q-learning) do not scale naturally to such settings. Policy-based methods directly parametrize a stochastic control law, the policy \( \pi_{\boldsymbol{\theta}}(a \mid s) \), and optimize its parameters \( \boldsymbol{\theta} \) by gradient ascent on an objective function.

Consider an episodic MDP representing a robot executing a task from reset to termination. Let a trajectory be \( \tau = (s_0,a_0,s_1,a_1,\dots,s_T) \) with return \( R(\tau) = \sum_{t=0}^{T-1} \gamma^t r_t \). The performance objective is

\[ J(\boldsymbol{\theta}) \;=\; \mathbb{E}_{\tau \sim p_{\boldsymbol{\theta}}} \big[ R(\tau) \big] \;=\; \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau), \]

where \( p_{\boldsymbol{\theta}}(\tau) \) is the distribution over trajectories induced by the policy \( \pi_{\boldsymbol{\theta}} \) and the environment (robot + world) dynamics. We will derive unbiased estimators of \( \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \) and study their variance properties, which are crucial for sample-hungry robotic domains.

flowchart TD
  S0["Start: init theta"] --> COL["Collect rollouts from robot/simulator"]
  COL --> RET["Compute returns R(t) or advantages A(t)"]
  RET --> GRAD["Estimate grad J(theta) from samples"]
  GRAD --> UPD["Update theta by gradient ascent"]
  UPD --> COND{"Converged or budget?"}
  COND -->|no| COL
  COND -->|yes| POL["Deploy learned policy on robot"]
        

We first derive the likelihood-ratio (REINFORCE) form of the policy gradient. The key idea is to differentiate the trajectory distribution with respect to the policy parameters.

The trajectory distribution factorizes as

\[ p_{\boldsymbol{\theta}}(\tau) \;=\; \rho_0(s_0) \prod_{t=0}^{T-1} \pi_{\boldsymbol{\theta}}(a_t \mid s_t) P(s_{t+1} \mid s_t, a_t), \]

where \( \rho_0 \) is the initial state distribution and \( P \) is the environment transition kernel, which we assume does not depend on \( \boldsymbol{\theta} \). Then

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \;=\; \nabla_{\boldsymbol{\theta}} \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau) \;=\; \sum_{\tau} R(\tau)\, \nabla_{\boldsymbol{\theta}} p_{\boldsymbol{\theta}}(\tau). \]

Using the identity \( \nabla_{\boldsymbol{\theta}} p = p \nabla_{\boldsymbol{\theta}} \log p \), we obtain

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \;=\; \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau) \nabla_{\boldsymbol{\theta}}\log p_{\boldsymbol{\theta}}(\tau) \;=\; \mathbb{E}_{\tau \sim p_{\boldsymbol{\theta}}} \big[ R(\tau)\, \nabla_{\boldsymbol{\theta}} \log p_{\boldsymbol{\theta}}(\tau) \big]. \]

Since only the policy depends on \( \boldsymbol{\theta} \), we have

\[ \log p_{\boldsymbol{\theta}}(\tau) \;=\; \log \rho_0(s_0) + \sum_{t=0}^{T-1} \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) + \sum_{t=0}^{T-1} \log P(s_{t+1}\mid s_t,a_t), \]

and thus

\[ \nabla_{\boldsymbol{\theta}} \log p_{\boldsymbol{\theta}}(\tau) \;=\; \sum_{t=0}^{T-1} \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a_t \mid s_t). \]

Plugging this into the gradient expression yields the REINFORCE gradient:

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \;=\; \mathbb{E}_{\tau \sim p_{\boldsymbol{\theta}}} \left[ R(\tau) \sum_{t=0}^{T-1} \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \right]. \]

For robotics, this result means we can treat the robot and its dynamics as a black box simulator and still obtain unbiased policy gradient estimates by sampling trajectories and computing \( \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}} \).

In practice, we approximate the expectation with Monte Carlo samples. Let \( G_t = \sum_{k=t}^{T-1} \gamma^{k-t} r_k \) be the return from time \( t \). A common estimator is

\[ \hat{g} \;=\; \frac{1}{N} \sum_{i=1}^N \sum_{t=0}^{T_i-1} \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t^{(i)} \mid s_t^{(i)}) \, G_t^{(i)}, \]

where index \( i \) ranges over sampled episodes. This estimator is unbiased but often has large variance, especially in long-horizon robotic tasks.

3.1. Baselines

A standard variance reduction technique is to subtract a baseline \( b(s_t) \) that does not depend on \( a_t \):

\[ \hat{g}_b \;=\; \frac{1}{N} \sum_{i=1}^N \sum_{t=0}^{T_i-1} \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t^{(i)} \mid s_t^{(i)}) \, \big(G_t^{(i)} - b(s_t^{(i)})\big). \]

Claim. For any baseline \( b \) independent of the action, the estimator \( \hat{g}_b \) is still unbiased: \( \mathbb{E}[\hat{g}_b] = \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \).

Proof sketch. Consider a single time step:

\[ \mathbb{E}_{a \sim \pi_{\boldsymbol{\theta}}(\cdot \mid s)} \big[ \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a \mid s) \, b(s) \big] \;=\; b(s) \sum_a \pi_{\boldsymbol{\theta}}(a \mid s) \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a \mid s). \]

Since \( \sum_a \pi_{\boldsymbol{\theta}}(a \mid s) = 1 \), we have \( \sum_a \nabla_{\boldsymbol{\theta}} \pi_{\boldsymbol{\theta}}(a \mid s) = 0 \). Using \( \nabla_{\boldsymbol{\theta}}\log \pi = \nabla_{\boldsymbol{\theta}}\pi / \pi \), it follows that the expectation is zero. Summing over time and trajectories preserves unbiasedness.

The optimal baseline in a scalar setting is the value function \( V^{\pi}(s_t) = \mathbb{E}[G_t \mid s_t] \), which motivates actor–critic methods.

3.2. Gaussian Policies for Continuous Actions

For torque control, a common policy is Gaussian:

\[ \pi_{\boldsymbol{\theta}}(a \mid s) \;=\; \mathcal{N}\big(a \mid \mu_{\boldsymbol{\theta}}(s), \sigma^2\big), \]

where the mean \( \mu_{\boldsymbol{\theta}}(s) \) is the output of a neural network (or linear function of features) and \( \sigma^2 \) is either fixed or learned. For a scalar action, the log-density is

\[ \log \pi_{\boldsymbol{\theta}}(a \mid s) \;=\; -\frac{1}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2} \big(a - \mu_{\boldsymbol{\theta}}(s)\big)^2. \]

Differentiating gives

\[ \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a \mid s) \;=\; \frac{a - \mu_{\boldsymbol{\theta}}(s)}{\sigma^2} \,\nabla_{\boldsymbol{\theta}} \mu_{\boldsymbol{\theta}}(s). \]

This expression is easy to compute via automatic differentiation in Python (PyTorch, JAX) and forms the backbone of many continuous-control algorithms for robot arms and legged robots.

Actor–critic algorithms introduce a parametric critic that approximates value functions and provides a low-variance baseline for the actor (policy). Let \( V_{\mathbf{w}}(s) \) be a differentiable function approximator with parameters \( \mathbf{w} \). The actor uses an advantage estimate \( \hat{A}_t \) instead of the raw return \( G_t \).

4.1. TD Error as Advantage

A simple one-step advantage estimate uses the temporal-difference (TD) error:

\[ \delta_t \;=\; r_t + \gamma V_{\mathbf{w}}(s_{t+1}) - V_{\mathbf{w}}(s_t). \]

For small step sizes, the update

\[ \mathbf{w} \leftarrow \mathbf{w} + \alpha_c \,\delta_t \,\nabla_{\mathbf{w}} V_{\mathbf{w}}(s_t) \]

performs stochastic gradient descent on the squared TD error. The actor then updates

\[ \boldsymbol{\theta} \leftarrow \boldsymbol{\theta} + \alpha_a \, \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \,\delta_t, \]

which can be seen as a policy gradient step where \( \delta_t \approx Q^{\pi}(s_t,a_t) - V^{\pi}(s_t) \) is an estimate of the advantage \( A^{\pi}(s_t,a_t) \). When the critic is accurate, this yields the correct gradient direction but with much reduced variance compared to Monte Carlo returns.

4.2. Architecture for Robotic Control

In robotic applications, the actor and critic are typically neural networks taking as input the robot state (joint angles and velocities, base pose, end-effector pose, etc.) and possibly task parameters (desired goal pose). The actor outputs the mean of a Gaussian over torques or desired velocities; the critic outputs a scalar value.

flowchart TD
  ENV["Robot + world"] --> ST["State s_t"]
  ST --> ACT["Actor pi_theta"]
  ACT --> ACTN["Action a_t (torques/velocities)"]
  ACTN --> ENV
  ENV --> REW["Reward r_t, next state s_{t+1}"]
  ST --> CRT["Critic V_w"]
  REW --> CRT
  CRT --> DELTA["TD error delta_t"]
  DELTA --> UPD_ACT["Update theta using grad log pi * delta_t"]
  DELTA --> UPD_CRT["Update w using delta_t * grad V_w"]
        

This loop can run either in simulation (e.g., PyBullet, MuJoCo) or on a real robot with appropriate safety constraints and reset mechanisms.

5.1. Action Scaling and Torque Limits

Robots have strict actuator limits. Typically, the actor outputs an unconstrained vector \( u \in \mathbb{R}^m \), which is squashed with \( \tanh \) and scaled:

\[ a_t \;=\; \mathbf{u}_{\max} \odot \tanh(u_t), \]

where \( \mathbf{u}_{\max} \) encodes torque or velocity limits per joint and \( \odot \) denotes elementwise product. The policy gradient must account for the Jacobian of this squashing transformation, which is automatically handled by modern autodiff frameworks.

5.2. Reward Shaping

For tasks such as reaching, a typical reward is

\[ r_t \;=\; - \lambda_{\text{pos}} \big\| x_{\text{ee}}(q_t) - x_{\text{goal}} \big\|_2^2 - \lambda_{\text{torque}} \| \tau_t \|_2^2, \]

where \( x_{\text{ee}}(q_t) \) is the end-effector position (obtained via forward kinematics from prior courses), \( x_{\text{goal}} \) is the target, and \( \tau_t \) the commanded torques. The quadratic regularization stabilizes learning.

5.3. Sample Efficiency vs Bias

Monte Carlo policy gradients are unbiased but sample-inefficient. Actor–critic using TD introduces bias (due to bootstrapping) but significantly improves sample efficiency, which is crucial when rollouts are expensive (e.g., with high-fidelity dynamic simulators and complex manipulators).

We illustrate an on-policy actor–critic method for a torque-controlled pendulum, which is a standard surrogate for a single-link robot arm. The environment can be instantiated via gymnasium or PyBullet (e.g., Pendulum-v1).

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim
import gymnasium as gym

env = gym.make("Pendulum-v1")  # continuous torque control
obs_dim = env.observation_space.shape[0]
act_dim = env.action_space.shape[0]

class Actor(nn.Module):
    def __init__(self, obs_dim, act_dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(obs_dim, 64),
            nn.Tanh(),
            nn.Linear(64, 64),
            nn.Tanh(),
            nn.Linear(64, act_dim)
        )
        # log standard deviation (state-independent for simplicity)
        self.log_std = nn.Parameter(torch.zeros(act_dim))

    def forward(self, obs):
        mu = self.net(obs)
        std = torch.exp(self.log_std)
        return mu, std

    def sample(self, obs):
        mu, std = self.forward(obs)
        dist = torch.distributions.Normal(mu, std)
        a = dist.rsample()          # reparameterized sample
        logp = dist.log_prob(a).sum(axis=-1)
        # squash to [-2, 2] (Pendulum action bounds)
        squashed = 2.0 * torch.tanh(a)
        return squashed, logp

class Critic(nn.Module):
    def __init__(self, obs_dim):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(obs_dim, 64),
            nn.Tanh(),
            nn.Linear(64, 64),
            nn.Tanh(),
            nn.Linear(64, 1)
        )

    def forward(self, obs):
        return self.net(obs).squeeze(-1)

actor = Actor(obs_dim, act_dim)
critic = Critic(obs_dim)
actor_opt = optim.Adam(actor.parameters(), lr=3e-4)
critic_opt = optim.Adam(critic.parameters(), lr=3e-3)
gamma = 0.99

def run_episode():
    obs, _ = env.reset()
    traj = []
    done = False
    while not done:
        obs_t = torch.as_tensor(obs, dtype=torch.float32)
        with torch.no_grad():
            a, logp = actor.sample(obs_t)
        a_np = a.numpy()
        next_obs, reward, terminated, truncated, _ = env.step(a_np)
        done = terminated or truncated
        traj.append((obs, a_np, reward, logp.item(), next_obs))
        obs = next_obs
    return traj

for episode in range(2000):
    traj = run_episode()

    # collect tensors
    obs = torch.as_tensor([t[0] for t in traj], dtype=torch.float32)
    logp = torch.as_tensor([t[3] for t in traj], dtype=torch.float32)
    rewards = [t[2] for t in traj]
    next_obs = torch.as_tensor([t[4] for t in traj], dtype=torch.float32)

    # compute returns and critic values
    returns = []
    G = 0.0
    for r in reversed(rewards):
        G = r + gamma * G
        returns.append(G)
    returns.reverse()
    returns = torch.as_tensor(returns, dtype=torch.float32)

    V = critic(obs)
    with torch.no_grad():
        V_next = critic(next_obs)
    # TD error as advantage
    r_t = torch.as_tensor(rewards, dtype=torch.float32)
    delta = r_t + gamma * V_next - V
    advantages = delta.detach()

    # critic loss (TD error squared)
    critic_loss = (delta ** 2).mean()
    critic_opt.zero_grad()
    critic_loss.backward()
    critic_opt.step()

    # actor loss: gradient ascent on J => minimize -E[log pi * advantage]
    actor_loss = -(logp * advantages).mean()
    actor_opt.zero_grad()
    actor_loss.backward()
    actor_opt.step()

    if episode % 100 == 0:
        print(f"Episode {episode}, return = {returns[0].item():.2f}")
      

This code uses Gaussian policies with squashing and a shared observation encoder. In a full robotic manipulator, the observation vector would be joint angles, velocities, and possibly end-effector pose and contact information, taken from the robot middleware (ROS, PyBullet, etc.).

For real-time control on embedded controllers, a lightweight linear-Gaussian actor–critic can be implemented using C++ with Eigen for linear algebra and a physics engine such as Bullet or ODE providing the dynamics.

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

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

struct Actor {
    MatrixXd W;   // weights: act_dim x obs_dim
    VectorXd b;   // bias: act_dim
    double log_std;

    Actor(int obs_dim, int act_dim)
        : W(MatrixXd::Zero(act_dim, obs_dim)),
          b(VectorXd::Zero(act_dim)),
          log_std(std::log(0.1)) {}

    VectorXd mean(const VectorXd& s) const {
        return W * s + b;
    }

    VectorXd sample(const VectorXd& s, double& logp,
                    std::mt19937& gen) const {
        VectorXd mu = mean(s);
        double std = std::exp(log_std);
        std::normal_distribution<double> dist(0.0, std);
        VectorXd eps(mu.size());
        for (int i = 0; i < mu.size(); ++i) {
            eps(i) = dist(gen);
        }
        VectorXd a = mu + eps;
        // log probability
        double log_det = mu.size() * std::log(std);
        double quad = (a - mu).squaredNorm() / (2.0 * std * std);
        logp = -0.5 * mu.size() * std::log(2.0 * M_PI)
               - log_det - quad;
        return a;
    }

    // gradient of log pi wrt W for linear Gaussian (no squashing)
    MatrixXd grad_logp_W(const VectorXd& s,
                         const VectorXd& a) const {
        double std = std::exp(log_std);
        VectorXd mu = mean(s);
        VectorXd diff = (a - mu) / (std * std);
        // outer product: diff * s'
        return diff * s.transpose();
    }

    VectorXd grad_logp_b(const VectorXd& s,
                         const VectorXd& a) const {
        double std = std::exp(log_std);
        VectorXd mu = mean(s);
        return (a - mu) / (std * std);
    }
};

struct Critic {
    VectorXd w;  // linear value: V(s) = w' s

    explicit Critic(int obs_dim) : w(VectorXd::Zero(obs_dim)) {}

    double value(const VectorXd& s) const {
        return w.dot(s);
    }
};

int main() {
    const int obs_dim = 4;  // example: pendulum state
    const int act_dim = 1;
    Actor actor(obs_dim, act_dim);
    Critic critic(obs_dim);
    double alpha_actor = 1e-3;
    double alpha_critic = 5e-3;
    double gamma = 0.99;
    std::mt19937 gen(42);

    // Pseudo-code loop:
    // for each episode:
    //   rollout states, actions, rewards, next_states using physics engine
    std::vector<VectorXd> states, next_states, actions;
    std::vector<double> rewards, logps;

    // ... fill buffers from simulator ...

    // Online update
    for (size_t t = 0; t < states.size(); ++t) {
        const VectorXd& s = states[t];
        const VectorXd& a = actions[t];
        const VectorXd& s_next = next_states[t];
        double r = rewards[t];

        double V_s = critic.value(s);
        double V_next = critic.value(s_next);
        double delta = r + gamma * V_next - V_s;

        // critic update
        critic.w += alpha_critic * delta * s;

        // actor update
        MatrixXd grad_W = actor.grad_logp_W(s, a);
        VectorXd grad_b = actor.grad_logp_b(s, a);
        actor.W += alpha_actor * delta * grad_W;
        actor.b += alpha_actor * delta * grad_b;
    }

    return 0;
}
      

This skeleton abstracts away the robot dynamics, which are provided by a simulation environment or actual hardware interface. The gradients use the closed-form expression for a linear-Gaussian policy.

8.1. Java Skeleton (Using Simple Linear Function Approximators)

In Java, an actor–critic loop can be built on top of a robotics simulator such as jMonkeyEngine or custom middleware. Below is a minimalistic sketch using plain arrays for a linear-Gaussian policy.

public class LinearActor {
    private final int obsDim;
    private final int actDim;
    private final double[][] W;  // actDim x obsDim
    private final double[] b;
    private double logStd = Math.log(0.1);

    public LinearActor(int obsDim, int actDim) {
        this.obsDim = obsDim;
        this.actDim = actDim;
        this.W = new double[actDim][obsDim];
        this.b = new double[actDim];
    }

    public double[] mean(double[] s) {
        double[] mu = new double[actDim];
        for (int i = 0; i < actDim; i++) {
            double v = b[i];
            for (int j = 0; j < obsDim; j++) {
                v += W[i][j] * s[j];
            }
            mu[i] = v;
        }
        return mu;
    }

    public double[] sample(double[] s, java.util.Random rng,
                           double[] outLogp) {
        double[] mu = mean(s);
        double std = Math.exp(logStd);
        double[] a = new double[actDim];
        double logp = 0.0;
        for (int i = 0; i < actDim; i++) {
            double eps = rng.nextGaussian() * std;
            a[i] = mu[i] + eps;
            double diff = a[i] - mu[i];
            logp += -0.5 * Math.log(2.0 * Math.PI * std * std)
                    - 0.5 * diff * diff / (std * std);
        }
        outLogp[0] = logp;
        return a;
    }

    public void update(double[] s, double[] a, double advantage,
                       double alpha) {
        double std = Math.exp(logStd);
        double invVar = 1.0 / (std * std);
        double[] mu = mean(s);
        for (int i = 0; i < actDim; i++) {
            double diff = (a[i] - mu[i]) * invVar;
            for (int j = 0; j < obsDim; j++) {
                W[i][j] += alpha * advantage * diff * s[j];
            }
            b[i] += alpha * advantage * diff;
        }
    }
}
      

8.2. MATLAB/Simulink — Using Reinforcement Learning Toolbox

MATLAB's Reinforcement Learning Toolbox provides built-in actor–critic agents that can be connected to Simulink models of robot dynamics (e.g., a manipulator modeled with Simscape Multibody). A minimal script for a continuous-control actor–critic agent is:

% Define observation and action spaces
obsInfo = rlNumericSpec([obsDim 1], ...
    'LowerLimit', -inf, 'UpperLimit', inf);
actInfo = rlNumericSpec([actDim 1], ...
    'LowerLimit', -uMax, 'UpperLimit', uMax);

% Create Simulink environment
mdl = "RobotPendulumModel";
open_system(mdl);
env = rlSimulinkEnv(mdl, "RL_Agent_Block", obsInfo, actInfo);

% Critic network (state-value)
statePath = [
    featureInputLayer(obsDim, "Normalization","none","Name","state")
    fullyConnectedLayer(64, "Name","fc1")
    reluLayer("Name","relu1")
    fullyConnectedLayer(64, "Name","fc2")
    reluLayer("Name","relu2")
    fullyConnectedLayer(1, "Name","value")
];
criticNet = dlnetwork(layerGraph(statePath));
critic = rlValueFunction(criticNet, obsInfo, ...
    "ObservationInputNames","state");

% Actor network (Gaussian mean)
actorPath = [
    featureInputLayer(obsDim, "Normalization","none","Name","state")
    fullyConnectedLayer(64, "Name","fc1")
    reluLayer("Name","relu1")
    fullyConnectedLayer(64, "Name","fc2")
    reluLayer("Name","relu2")
    fullyConnectedLayer(actDim, "Name","mu")
    tanhLayer("Name","tanh")
];
actorNet = dlnetwork(layerGraph(actorPath));
actor = rlContinuousGaussianActor(actorNet, obsInfo, actInfo, ...
    "ObservationInputNames","state");

agentOpts = rlACAgentOptions( ...
    "DiscountFactor", 0.99, ...
    "ActorOptimizerOptions", rlOptimizerOptions("LearnRate",1e-4), ...
    "CriticOptimizerOptions", rlOptimizerOptions("LearnRate",1e-3));

agent = rlACAgent(actor, critic, agentOpts);

% Train
trainOpts = rlTrainingOptions( ...
    "MaxEpisodes", 2000, ...
    "MaxStepsPerEpisode", 400, ...
    "StopTrainingCriteria","AverageReward", ...
    "StopTrainingValue",-50);

trainingStats = train(agent, env, trainOpts);
      

Here RobotPendulumModel is a Simulink model with a torque input and pendulum state output; the RL Agent block interfaces with the actor and critic described above.

8.3. Wolfram Mathematica — REINFORCE for a Toy MDP

Mathematica can be used to prototype policy gradient algorithms for low-dimensional systems using symbolic or numerical differentiation.

(* Simple 1D state, 1D action Gaussian policy *)
Clear[mu, theta, sigma, piLogGrad];

mu[s_, theta_] := theta[[1]]*s + theta[[2]];
sigma = 0.5;

(* gradient of log pi wrt theta *)
piLogGrad[s_, a_, theta_List] := Module[
  {m = mu[s, theta], diff},
  diff = (a - m)/(sigma^2);
  {diff*s, diff}
];

(* rollout of finite-horizon episode *)
rollout[theta_List, envStep_, s0_, horizon_] := Module[
  {s = s0, data = {}, a, r, sNext},
  Do[
    a = RandomVariate[NormalDistribution[mu[s, theta], sigma]];
    {sNext, r} = envStep[s, a];
    AppendTo[data, <|
      "s" -> s, "a" -> a, "r" -> r, "sNext" -> sNext
    |>];
    s = sNext;
    ,
    {t, 1, horizon}
  ];
  data
];

(* return from t *)
discountedReturns[data_List, gamma_] := Module[
  {R = 0, returns = ConstantArray[0, Length[data]]},
  Do[
    R = data[[k, "r"]] + gamma*R;
    returns[[k]] = R;
    ,
    {k, Length[data], 1, -1}
  ];
  returns
];

(* REINFORCE update *)
reinforceUpdate[theta_List, envStep_, s0_, horizon_, alpha_, gamma_] :=
 Module[{data, G, grad = {0., 0.}},
  data = rollout[theta, envStep, s0, horizon];
  G = discountedReturns[data, gamma];
  Do[
    grad += G[[t]] * piLogGrad[
      data[[t, "s"]], data[[t, "a"]], theta
    ];
    ,
    {t, 1, Length[data]}
  ];
  theta + alpha * grad
];

(* Example environment: quadratic cost *)
envStep[s_, a_] := Module[
  {sNext, r},
  sNext = s + 0.1*a;
  r = - (sNext^2 + 0.01*a^2);
  {sNext, r}
];

theta0 = {0., 0.};
theta = Nest[
  reinforceUpdate[#, envStep, 0.5, 50, 0.001, 0.99] &,
  theta0, 500];

theta
      

This code treats the dynamics as a black box envStep and approximates the policy gradient via Monte Carlo. Extending this to multi-dimensional robot states is straightforward conceptually, though computationally heavier.

Problem 1 (Trajectory Policy Gradient Derivation). Starting from the definition \( J(\boldsymbol{\theta}) = \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau) \), derive the REINFORCE policy gradient expression \( \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) = \mathbb{E}_{\tau} \big[ R(\tau) \sum_t \nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \big] \) .

Solution.

We compute \( \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \):

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) = \nabla_{\boldsymbol{\theta}} \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau) = \sum_{\tau} R(\tau) \nabla_{\boldsymbol{\theta}} p_{\boldsymbol{\theta}}(\tau). \]

Using \( \nabla p = p \nabla \log p \), we obtain

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) = \sum_{\tau} p_{\boldsymbol{\theta}}(\tau) R(\tau) \nabla_{\boldsymbol{\theta}} \log p_{\boldsymbol{\theta}}(\tau) = \mathbb{E}_{\tau} \big[ R(\tau) \nabla_{\boldsymbol{\theta}} \log p_{\boldsymbol{\theta}}(\tau) \big]. \]

Since \( \log p_{\boldsymbol{\theta}}(\tau) = \sum_t \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) + \text{const} \) , differentiating gives \( \nabla_{\boldsymbol{\theta}}\log p_{\boldsymbol{\theta}}(\tau) = \sum_t \nabla_{\boldsymbol{\theta}}\log\pi_{\boldsymbol{\theta}}(a_t \mid s_t) \) . Substituting yields the desired result.

Problem 2 (Baseline Invariance). Let \( b(s_t) \) be any function of the state only. Show that

\[ \mathbb{E}\left[ \sum_t \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \, b(s_t) \right] = 0. \]

Solution.

Condition on \( s_t \). Then

\[ \mathbb{E} \big[ \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a_t \mid s_t) \, b(s_t) \mid s_t \big] = b(s_t) \sum_a \pi_{\boldsymbol{\theta}}(a \mid s_t) \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a \mid s_t). \]

Using \( \nabla_{\boldsymbol{\theta}}\log \pi(a \mid s_t) = \nabla_{\boldsymbol{\theta}} \pi(a \mid s_t) / \pi(a \mid s_t) \), the sum becomes \( \sum_a \nabla_{\boldsymbol{\theta}}\pi(a \mid s_t) \), which is \( \nabla_{\boldsymbol{\theta}} \sum_a \pi(a \mid s_t) = \nabla_{\boldsymbol{\theta}} 1 = 0 \). Therefore the conditional expectation is zero for each \( s_t \), and so the unconditional expectation is zero, proving baseline invariance.

Problem 3 (Gaussian Policy Gradient). For the scalar Gaussian policy \( \pi_{\boldsymbol{\theta}}(a \mid s) = \mathcal{N}(a \mid \mu_{\boldsymbol{\theta}}(s), \sigma^2) \) , derive \( \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a \mid s) = \frac{a - \mu_{\boldsymbol{\theta}}(s)}{\sigma^2} \nabla_{\boldsymbol{\theta}} \mu_{\boldsymbol{\theta}}(s) \) .

Solution.

We have

\[ \log \pi_{\boldsymbol{\theta}}(a \mid s) = -\frac{1}{2}\log(2\pi\sigma^2) -\frac{1}{2\sigma^2} \big(a - \mu_{\boldsymbol{\theta}}(s)\big)^2. \]

The first term does not depend on \( \boldsymbol{\theta} \). Differentiating the second term,

\[ \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a \mid s) = -\frac{1}{2\sigma^2} \nabla_{\boldsymbol{\theta}} \big(a - \mu_{\boldsymbol{\theta}}(s)\big)^2 = -\frac{1}{2\sigma^2} \cdot 2 \big(a - \mu_{\boldsymbol{\theta}}(s)\big) \big(-\nabla_{\boldsymbol{\theta}} \mu_{\boldsymbol{\theta}}(s)\big), \]

which simplifies to \( \frac{a - \mu_{\boldsymbol{\theta}}(s)}{\sigma^2} \nabla_{\boldsymbol{\theta}} \mu_{\boldsymbol{\theta}}(s) \) as required.

Problem 4 (Actor–Critic as Policy Gradient). Assume the critic gives exact values \( V_{\mathbf{w}}(s) = V^{\pi}(s) \). Show that using the TD error \( \delta_t = r_t + \gamma V^{\pi}(s_{t+1}) - V^{\pi}(s_t) \) in the actor update yields an unbiased estimate of the policy gradient that uses the advantage \( A^{\pi}(s_t,a_t) = Q^{\pi}(s_t,a_t) - V^{\pi}(s_t) \).

Solution.

Under exact values, the expected TD error conditioned on \( s_t,a_t \) satisfies

\[ \mathbb{E}[\delta_t \mid s_t, a_t] = \mathbb{E}\big[ r_t + \gamma V^{\pi}(s_{t+1}) - V^{\pi}(s_t) \mid s_t, a_t \big] = Q^{\pi}(s_t,a_t) - V^{\pi}(s_t) = A^{\pi}(s_t,a_t). \]

Therefore, replacing \( G_t - V^{\pi}(s_t) \) with \( \delta_t \) in the policy gradient estimator preserves the expected value, yielding an unbiased estimator of \( \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) \).

Problem 5 (One-Step MDP and Gradient). Consider a one-step MDP modeling a single control decision for a manipulator: the state \( s \) is observed once, an action \( a \) is sampled from \( \pi_{\boldsymbol{\theta}}(a \mid s) \), and reward \( r(s,a) \) is received. Show that the policy gradient reduces to the gradient of the expected reward \( J(\boldsymbol{\theta}) = \mathbb{E}_{a \sim \pi_{\boldsymbol{\theta}}(\cdot \mid s)} [r(s,a)] \) .

Solution.

The expected reward is

\[ J(\boldsymbol{\theta}) = \sum_a \pi_{\boldsymbol{\theta}}(a \mid s) r(s,a). \]

Differentiating gives

\[ \nabla_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) = \sum_a r(s,a) \nabla_{\boldsymbol{\theta}} \pi_{\boldsymbol{\theta}}(a \mid s) = \sum_a r(s,a) \pi_{\boldsymbol{\theta}}(a \mid s) \nabla_{\boldsymbol{\theta}} \log \pi_{\boldsymbol{\theta}}(a \mid s), \]

which is exactly the REINFORCE form \( \mathbb{E}_{a \sim \pi_{\boldsymbol{\theta}}} [r(s,a)\nabla_{\boldsymbol{\theta}}\log \pi_{\boldsymbol{\theta}}(a \mid s)] \) . Thus, in the one-step case, policy gradient ascent coincides with gradient ascent on the expected reward.

In this lesson we formulated policy-based reinforcement learning for robotics, derived the trajectory policy gradient theorem, and established the REINFORCE estimator. We showed how baselines preserve unbiasedness while dramatically reducing variance, and built actor–critic methods where a learned critic approximates value functions and provides TD-based advantage estimates. Finally, we sketched practical implementations across Python, C++, Java, MATLAB/Simulink, and Mathematica, emphasizing continuous-action Gaussian policies appropriate for torque- and velocity-controlled robots. These tools form the basis for more advanced off-policy and entropy-regularized methods in the next lesson.