Consider a robot manipulator or mobile base with generalized configuration \( \mathbf{q} \in \mathbb{R}^n \) and generalized velocities \( \dot{\mathbf{q}} \in \mathbb{R}^n \). A standard state choice is \( \mathbf{x} = (\mathbf{q}, \dot{\mathbf{q}}) \in \mathbb{R}^{2n} \), and the control input (joint torques, wheel torques, or velocity commands) is \( \mathbf{u} \in \mathbb{R}^m \).

In continuous time, rigid-body dynamics may be written as

\[ \mathbf{M}(\mathbf{q}) \ddot{\mathbf{q}} + \mathbf{C}(\mathbf{q}, \dot{\mathbf{q}})\dot{\mathbf{q}} + \mathbf{g}(\mathbf{q}) = \mathbf{B}(\mathbf{q})\mathbf{u} + \mathbf{w}(t), \]

where \( \mathbf{w}(t) \) is a disturbance term. For reinforcement learning (RL) and MDP theory we almost always work with a discrete-time model, obtained for example by Euler integration with step size \( \Delta t > 0 \):

\[ \mathbf{x}_{t+1} = f(\mathbf{x}_t, \mathbf{u}_t) + \boldsymbol{\eta}_t, \]

where \( \boldsymbol{\eta}_t \) is a discrete-time noise term that captures modeling error, unmodeled dynamics, and sensor noise (if we embed estimation into the state).

A continuous control task specifies:

• A state space \( \mathcal{X} \subseteq \mathbb{R}^{d_x} \)
• An action space \( \mathcal{U} \subseteq \mathbb{R}^{d_u} \)
• Dynamics \( f : \mathcal{X} \times \mathcal{U} \to \mathcal{X} \)
• A scalar performance measure, typically a cumulative cost or reward

The MDP framework takes these ingredients and adds a precise probabilistic semantics for transitions and returns, enabling convergence proofs and algorithm design in later lessons.

flowchart TD
  A["Continuous control task (tracking, balancing, etc.)"]
    --> B["Choose state x = (q, dq, extra)"]
  B --> C["Choose action u (torques, velocities, etc.)"]
  C --> D["Specify dynamics x_next = f(x, u) + noise"]
  D --> E["Define reward r(x, u) and discount gamma"]
  E --> F["Obtain MDP (X, U, p, r, gamma) for RL"]
        

Formally, a (discounted) MDP suitable for robot control is a tuple \( \mathcal{M} = (\mathcal{X}, \mathcal{U}, p, r, \gamma, \rho_0) \), where:

• \( \mathcal{X} \) is a measurable state space (often a subset of \( \mathbb{R}^{d_x} \)).
• \( \mathcal{U} \) is a measurable action space (subset of \( \mathbb{R}^{d_u} \)).
• \( p(\cdot \mid \mathbf{x}, \mathbf{u}) \) is a stochastic kernel on \( \mathcal{X} \), giving the distribution of \( \mathbf{x}_{t+1} \) when control \( \mathbf{u}_t = \mathbf{u} \) is applied in state \( \mathbf{x}_t = \mathbf{x} \).
• \( r : \mathcal{X} \times \mathcal{U} \to \mathbb{R} \) is a (possibly unbounded) reward function.
• \( \gamma \in (0,1) \) is a discount factor.
• \( \rho_0 \) is an initial state distribution on \( \mathcal{X} \).

A (stationary) policy is a measurable mapping \( \pi : \mathcal{X} \to \mathcal{P}(\mathcal{U}) \), where \( \mathcal{P}(\mathcal{U}) \) denotes the set of probability measures on the action space. At time \( t \), the controller samples \( \mathbf{u}_t \sim \pi(\cdot \mid \mathbf{x}_t) \).

The discounted return of a policy \( \pi \) from time \( t \) is

\[ G_t = \sum_{k=0}^{\infty} \gamma^{k} r(\mathbf{x}_{t+k}, \mathbf{u}_{t+k}). \]

The associated state-value and action-value functions are

\[ V^{\pi}(\mathbf{x}) = \mathbb{E}\!\left[ G_0 \mid \mathbf{x}_0 = \mathbf{x}, \pi \right], \quad Q^{\pi}(\mathbf{x}, \mathbf{u}) = \mathbb{E}\!\left[ G_0 \mid \mathbf{x}_0 = \mathbf{x}, \mathbf{u}_0 = \mathbf{u}, \pi \right]. \]

In continuous spaces, expectations above become integrals over \( \mathcal{X} \) and \( \mathcal{U} \). The existence of optimal policies and the well-posedness of the problem require technical assumptions (bounded rewards, suitable continuity of the kernel, etc.), which we will invoke but not fully prove in this introductory lesson.

For a robot whose discrete-time dynamics are \( \mathbf{x}_{t+1} = f(\mathbf{x}_t, \mathbf{u}_t) + \boldsymbol{\eta}_t \), where \( \boldsymbol{\eta}_t \) is independent noise with density \( \phi(\cdot) \), we can define the transition kernel

\[ p(\mathbf{x}' \mid \mathbf{x}, \mathbf{u}) = \phi\big(\mathbf{x}' - f(\mathbf{x}, \mathbf{u})\big). \]

In the deterministic case (no noise), the kernel is a Dirac measure

\[ p(\mathbf{x}' \mid \mathbf{x}, \mathbf{u}) = \delta\big(\mathbf{x}' - f(\mathbf{x}, \mathbf{u})\big), \]

meaning that \( \mathbf{x}_{t+1} \) is almost surely equal to \( f(\mathbf{x}_t, \mathbf{u}_t) \). Both models fit perfectly into the MDP definition.

To relate this to the underlying manipulator dynamics, define \( \mathbf{x} = (\mathbf{q}, \dot{\mathbf{q}}) \) and use a time step \( \Delta t \). A simple explicit Euler scheme gives

\[ \begin{aligned} \ddot{\mathbf{q}}_t &= \mathbf{M}(\mathbf{q}_t)^{-1} \Big(\mathbf{B}(\mathbf{q}_t)\mathbf{u}_t - \mathbf{C}(\mathbf{q}_t,\dot{\mathbf{q}}_t)\dot{\mathbf{q}}_t - \mathbf{g}(\mathbf{q}_t)\Big),\\ \dot{\mathbf{q}}_{t+1} &= \dot{\mathbf{q}}_t + \Delta t\,\ddot{\mathbf{q}}_t,\\ \mathbf{q}_{t+1} &= \mathbf{q}_t + \Delta t\,\dot{\mathbf{q}}_{t+1}. \end{aligned} \]

This defines \( f(\mathbf{x}_t, \mathbf{u}_t) \) and thus the kernel \( p(\cdot \mid \mathbf{x}_t, \mathbf{u}_t) \). The Markov property holds because the next state depends only on the current state and action, not on the full history:

\[ \mathbb{P}\big(\mathbf{x}_{t+1} \in A \mid \mathbf{x}_0, \mathbf{u}_0, \ldots, \mathbf{x}_t, \mathbf{u}_t\big) = \mathbb{P}\big(\mathbf{x}_{t+1} \in A \mid \mathbf{x}_t, \mathbf{u}_t\big) \]

for any measurable set \( A \subseteq \mathcal{X} \). This is a crucial structural property for applying MDP tools.

In many robotics texts, we formulate performance as cost minimization:

\[ J^{\pi}(\mathbf{x}_0) = \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^{t}\, c(\mathbf{x}_t, \mathbf{u}_t) \,\Big|\, \mathbf{x}_0, \pi \right], \]

where \( c \) is a nonnegative running cost. In RL, it is conventional to work with rewards \( r = -c \). The optimization problems

\[ \min_{\pi} J^{\pi}(\mathbf{x}_0) \quad \text{and} \quad \max_{\pi} V^{\pi}(\mathbf{x}_0) \]

are equivalent under the transformation \( r = -c \). For instance, tracking a desired configuration trajectory \( \mathbf{q}^{\star}_t \) with penalized torque usage might use

\[ c(\mathbf{x}_t, \mathbf{u}_t) = (\mathbf{q}_t - \mathbf{q}^{\star}_t)^{\top}\mathbf{Q} (\mathbf{q}_t - \mathbf{q}^{\star}_t) + \mathbf{u}_t^{\top}\mathbf{R}\mathbf{u}_t, \]

with positive semidefinite matrices \( \mathbf{Q} \) and positive definite \( \mathbf{R} \). Then

\[ r(\mathbf{x}_t, \mathbf{u}_t) = -c(\mathbf{x}_t,\mathbf{u}_t) \]

is a natural reward. Detailed reward engineering strategies for sparse or delayed signals will appear in a later lesson; here we only need the structural mapping from robotics objectives to MDP rewards.

For a fixed policy \( \pi \), the value function satisfies the Bellman equation

\[ V^{\pi}(\mathbf{x}) = \int_{\mathcal{U}} \pi(\mathrm{d}\mathbf{u} \mid \mathbf{x}) \left[ r(\mathbf{x}, \mathbf{u}) + \gamma \int_{\mathcal{X}} p(\mathrm{d}\mathbf{x}' \mid \mathbf{x}, \mathbf{u}) V^{\pi}(\mathbf{x}') \right]. \]

Define the Bellman operator \( T^{\pi} \) acting on bounded measurable functions \( V : \mathcal{X} \to \mathbb{R} \) by

\[ (T^{\pi} V)(\mathbf{x}) = \int_{\mathcal{U}} \pi(\mathrm{d}\mathbf{u} \mid \mathbf{x}) \left[ r(\mathbf{x}, \mathbf{u}) + \gamma \int_{\mathcal{X}} p(\mathrm{d}\mathbf{x}' \mid \mathbf{x}, \mathbf{u}) V(\mathbf{x}') \right]. \]

Under bounded rewards, \( T^{\pi} \) is a contraction on the Banach space of bounded functions with sup norm. Specifically,

\[ \lVert T^{\pi} V_1 - T^{\pi} V_2 \rVert_{\infty} \leq \gamma \lVert V_1 - V_2 \rVert_{\infty}. \]

Sketch of proof. For any \( \mathbf{x} \),

\[ \begin{aligned} (T^{\pi}V_1)(\mathbf{x}) - (T^{\pi}V_2)(\mathbf{x}) &= \gamma \int_{\mathcal{U}} \pi(\mathrm{d}\mathbf{u}\mid \mathbf{x}) \int_{\mathcal{X}} p(\mathrm{d}\mathbf{x}' \mid \mathbf{x},\mathbf{u}) \big(V_1(\mathbf{x}')-V_2(\mathbf{x}')\big) \\ &\leq \gamma \lVert V_1 - V_2 \rVert_{\infty} \int_{\mathcal{U}} \pi(\mathrm{d}\mathbf{u}\mid \mathbf{x}) \int_{\mathcal{X}} p(\mathrm{d}\mathbf{x}' \mid \mathbf{x},\mathbf{u}) \\ &= \gamma \lVert V_1 - V_2 \rVert_{\infty}, \end{aligned} \]

and similarly for the negative difference, giving the sup norm inequality. Banach's fixed-point theorem then guarantees existence and uniqueness of \( V^{\pi} \) solving \( V^{\pi} = T^{\pi} V^{\pi} \).

The optimal Bellman operator is

\[ (T V)(\mathbf{x}) = \sup_{\mathbf{u} \in \mathcal{U}} \left[ r(\mathbf{x}, \mathbf{u}) + \gamma \int_{\mathcal{X}} p(\mathrm{d}\mathbf{x}'\mid \mathbf{x},\mathbf{u}) V(\mathbf{x}') \right], \]

and its unique fixed point is the optimal value function \( V^{\star} \). Computing or approximating this fixed point in continuous spaces motivates the RL algorithms to be introduced next (policy gradients, actor–critic, and off-policy methods).

The MDP description can be viewed as an abstraction of the physical control loop that your robot executes: sense state, choose control, apply dynamics, receive a scalar signal summarizing task progress, and iterate.

flowchart TD
  S["State x_t (q, dq, sensors)"] --> P["Controller or policy pi(u | x_t)"]
  P --> A["Apply action u_t to robot"]
  A --> D["Dynamics give x_{t+1} = f(x_t, u_t) + noise"]
  D --> R["Compute reward r(x_t, u_t)"]
  R --> L["Learning or design updates controller (later lessons)"]
  D --> S2["Next state x_{t+1}"]
  S2 --> P
        

In this lesson we formalize \( \mathcal{X} \), \( \mathcal{U} \), \( p \), and \( r \); later lessons will focus on how to design or learn the policy \( \pi \) using gradient-based and off-policy algorithms.

We now build a minimal Python abstraction for a continuous-state robot MDP. This is intentionally structural and does not yet include any learning. In practice, such an environment would be exposed via interfaces compatible with libraries like gymnasium, mujoco, or pybullet.

import numpy as np

class ContinuousRobotMDP:
    """
    Minimal continuous-state, continuous-action MDP wrapper for a robot.

    State x = [q, dq] in R^{2n}
    Action u in R^{m} (e.g. joint torques)
    """

    def __init__(self, n_dof, m_act, dt=0.02):
        self.n = n_dof
        self.m = m_act
        self.dt = dt

        # Cost weights for quadratic running cost
        self.Q = np.eye(2 * n_dof)   # state cost
        self.R = 0.01 * np.eye(m_act)  # action cost

        # Simple joint limits and action limits
        self.q_limit = np.ones(n_dof) * np.pi
        self.u_limit = np.ones(m_act) * 10.0

    def dynamics(self, x, u):
        """
        Simple double-integrator style dynamics per joint:
          q_{t+1} = q_t + dt * dq_t
          dq_{t+1} = dq_t + dt * ddq_t
        where ddq_t approximates M^{-1} B u.
        """
        x = np.asarray(x, dtype=float)
        u = np.asarray(u, dtype=float)

        n = self.n
        q = x[:n]
        dq = x[n:]

        # Clip control to limits
        u_clipped = np.clip(u, -self.u_limit, self.u_limit)

        # Crude acceleration model (unit inertia)
        ddq = u_clipped  # in real robots, use full rigid-body dynamics

        q_next = q + self.dt * dq
        dq_next = dq + self.dt * ddq

        # Enforce simple joint limits
        q_next = np.clip(q_next, -self.q_limit, self.q_limit)

        x_next = np.concatenate([q_next, dq_next])
        return x_next

    def reward(self, x, u, x_des=None):
        """
        Quadratic tracking reward in state and control.
        For now, desired state x_des is the zero vector if not given.
        """
        x = np.asarray(x, dtype=float)
        u = np.asarray(u, dtype=float)

        if x_des is None:
            x_des = np.zeros_like(x)

        dx = x - x_des
        cost = dx.T @ self.Q @ dx + u.T @ self.R @ u
        return -float(cost)

    def step(self, x, u, x_des=None):
        """
        Apply one MDP transition:
          x_next ~ p(. | x, u)
        Here the dynamics are deterministic; stochasticity could be added
        by injecting Gaussian noise into x_next.
        """
        x_next = self.dynamics(x, u)
        r = self.reward(x, u, x_des)

        # Example termination condition: huge norm means failure
        done = bool(np.linalg.norm(x_next) > 1e3)

        return x_next, r, done

# Example usage:
if __name__ == "__main__":
    env = ContinuousRobotMDP(n_dof=2, m_act=2, dt=0.01)
    x = np.array([0.1, -0.1, 0.0, 0.0])  # [q1, q2, dq1, dq2]
    u = np.array([1.0, -0.5])
    x_next, r, done = env.step(x, u)
    print("x_next:", x_next)
    print("reward:", r, "done:", done)
      

In later lessons, this environment interface will be connected to RL algorithms (policy gradients, actor–critic, off-policy methods). For now, it should be understood as a concrete example of the abstract MDP \( (\mathcal{X}, \mathcal{U}, p, r, \gamma) \).

In C++, continuous control in robotics commonly uses Eigen for linear algebra and dedicated rigid-body libraries such as Pinocchio or RBDL. Below is a minimal header-only style abstraction of an MDP transition for a robot. The actual dynamics could call into such libraries.

#include <Eigen/Dense>
#include <cmath>

struct RobotState {
    Eigen::VectorXd x; // [q; dq]
};

struct RobotAction {
    Eigen::VectorXd u; // torques or commands
};

class RobotMDP {
public:
    RobotMDP(int n_dof, int m_act, double dt)
        : n_(n_dof), m_(m_act), dt_(dt)
    {
        Q_ = Eigen::MatrixXd::Identity(2 * n_, 2 * n_);
        R_ = 0.01 * Eigen::MatrixXd::Identity(m_, m_);
        q_limit_ = Eigen::VectorXd::Constant(n_, M_PI);
        u_limit_ = Eigen::VectorXd::Constant(m_, 10.0);
    }

    RobotState dynamics(const RobotState& s, const RobotAction& a) const {
        Eigen::VectorXd q = s.x.head(n_);
        Eigen::VectorXd dq = s.x.tail(n_);

        Eigen::VectorXd u = a.u;
        u = u.cwiseMax(-u_limit_).cwiseMin(u_limit_);

        Eigen::VectorXd ddq = u; // placeholder for M(q)^{-1} B u

        Eigen::VectorXd q_next = q + dt_ * dq;
        Eigen::VectorXd dq_next = dq + dt_ * ddq;
        q_next = q_next.cwiseMax(-q_limit_).cwiseMin(q_limit_);

        RobotState s_next;
        s_next.x.resize(2 * n_);
        s_next.x << q_next, dq_next;
        return s_next;
    }

    double reward(const RobotState& s, const RobotAction& a,
                  const RobotState* s_des = nullptr) const {
        Eigen::VectorXd x_des;
        if (s_des) {
            x_des = s_des->x;
        } else {
            x_des = Eigen::VectorXd::Zero(2 * n_);
        }

        Eigen::VectorXd dx = s.x - x_des;
        double cost = dx.transpose() * Q_ * dx
                    + a.u.transpose() * R_ * a.u;
        return -cost;
    }

    // One MDP step: (s, a) -> (s_next, r)
    void step(const RobotState& s, const RobotAction& a,
              RobotState& s_next, double& r, bool& done) const {
        s_next = dynamics(s, a);
        r = reward(s, a, nullptr);

        done = (s_next.x.norm() > 1e3);
    }

private:
    int n_, m_;
    double dt_;
    Eigen::MatrixXd Q_, R_;
    Eigen::VectorXd q_limit_, u_limit_;
};
      

This C++ abstraction mirrors the mathematical structure of the MDP. Learning algorithms in C++ (e.g., custom policy gradient implementations or bindings to existing RL libraries) would interact with step in exactly the way the Bellman equations formalize.

Java is less commonly used for high-end robotics, but the same MDP abstractions apply. Numerical linear algebra can be handled by libraries such as EJML or Apache Commons Math. Below is a simple Java MDP interface and one concrete implementation with a quadratic tracking reward.

public interface ContinuousMDP {
    int stateDim();
    int actionDim();

    double[] step(double[] x, double[] u, StepResult result);
}

class StepResult {
    public double reward;
    public boolean done;
}

public class SimpleRobotMDP implements ContinuousMDP {
    private final int n;   // dof
    private final int m;   // actuation dimension
    private final double dt;
    private final double[] qLimit;
    private final double[] uLimit;

    public SimpleRobotMDP(int nDof, int mAct, double dt) {
        this.n = nDof;
        this.m = mAct;
        this.dt = dt;
        this.qLimit = new double[n];
        this.uLimit = new double[m];
        for (int i = 0; i < n; ++i) qLimit[i] = Math.PI;
        for (int j = 0; j < m; ++j) uLimit[j] = 10.0;
    }

    @Override
    public int stateDim() { return 2 * n; }

    @Override
    public int actionDim() { return m; }

    private double[] dynamics(double[] x, double[] u) {
        double[] q = new double[n];
        double[] dq = new double[n];
        for (int i = 0; i < n; ++i) {
            q[i] = x[i];
            dq[i] = x[i + n];
        }

        double[] uClip = new double[m];
        for (int j = 0; j < m; ++j) {
            uClip[j] = Math.max(-uLimit[j], Math.min(uLimit[j], u[j]));
        }

        double[] ddq = new double[n];
        for (int i = 0; i < n; ++i) {
            ddq[i] = uClip[Math.min(i, m - 1)];
        }

        double[] qNext = new double[n];
        double[] dqNext = new double[n];
        for (int i = 0; i < n; ++i) {
            qNext[i] = q[i] + dt * dq[i];
            dqNext[i] = dq[i] + dt * ddq[i];
            qNext[i] = Math.max(-qLimit[i], Math.min(qLimit[i], qNext[i]));
        }

        double[] xNext = new double[2 * n];
        for (int i = 0; i < n; ++i) {
            xNext[i] = qNext[i];
            xNext[i + n] = dqNext[i];
        }
        return xNext;
    }

    private double reward(double[] x, double[] u) {
        // Quadratic penalty toward zero state and small actions
        double cost = 0.0;
        for (int i = 0; i < x.length; ++i) {
            cost += x[i] * x[i];
        }
        for (int j = 0; j < u.length; ++j) {
            cost += 0.01 * u[j] * u[j];
        }
        return -cost;
    }

    @Override
    public double[] step(double[] x, double[] u, StepResult result) {
        double[] xNext = dynamics(x, u);
        result.reward = reward(x, u);
        double norm = 0.0;
        for (int i = 0; i < xNext.length; ++i) {
            norm += xNext[i] * xNext[i];
        }
        result.done = (norm > 1e6);
        return xNext;
    }
}
      

This mirrors the MDP semantics: step implements the kernel \( p \) and reward \( r \). A separate policy class would implement the mapping \( \pi(\mathbf{u}\mid \mathbf{x}) \), and an RL algorithm (not yet introduced) would update the policy based on observed transitions.

In MATLAB, we can represent the MDP transition as a function and use Simulink to capture the continuous-time dynamics. The Reinforcement Learning Toolbox provides rlSimulinkEnv, which wraps Simulink models into MDP-like environments. Below is a purely MATLAB function that implements a discrete-time MDP step.

function [x_next, r, done] = robot_mdp_step(x, u, dt, Q, R)
% robot_mdp_step  One transition of a robot MDP.
%   x     : state vector [q; dq]
%   u     : control vector
%   dt    : time step
%   Q, R  : cost weighting matrices

n = numel(x) / 2;
q = x(1:n);
dq = x(n+1:end);

% Joint and action limits
q_limit = pi * ones(n, 1);
u_limit = 10 * ones(size(u));

u = max(-u_limit, min(u_limit, u));

% Simple double-integrator dynamics
ddq = u;  % placeholder for full rigid-body dynamics

dq_next = dq + dt * ddq;
q_next = q + dt * dq_next;

q_next = max(-q_limit, min(q_limit, q_next));

x_next = [q_next; dq_next];

% Quadratic cost and reward
dx = x_next;  % track zero state
cost = dx' * Q * dx + u' * R * u;
r = -cost;

done = norm(x_next) > 1e3;
end
      

In Simulink, a typical structure is:

• A block computing the manipulator dynamics (e.g. via Robotics System Toolbox).
• Integrator blocks to obtain \( \mathbf{q} \) and \( \dot{\mathbf{q}} \).
• A block for the reward calculation \( r(\mathbf{x}_t, \mathbf{u}_t) \).
• Input ports for the action \( \mathbf{u}_t \) and output ports for state and reward.

This diagram, when wrapped by rlSimulinkEnv, furnishes an MDP interface where each simulation step corresponds to a transition drawn from the kernel \( p \) derived from the underlying dynamics.

Wolfram Mathematica is well suited to symbolic manipulations of dynamics and numerical simulation. Below, we define a simple one-degree-of-freedom robot as an MDP with Euler integration and a quadratic reward.

(* State: x = {q, dq}; action: u *)
Clear[robotStep]

robotStep[x_List, u_, dt_, qLimit_, uLimit_, qWeight_, dqWeight_, rWeight_] :=
 Module[{q, dq, uClip, ddq, qNext, dqNext, xNext, cost, reward, done},
  {q, dq} = x;
  uClip = Clip[u, {-uLimit, uLimit}];

  (* Simple double-integrator model: ddq = uClip *)
  ddq = uClip;

  dqNext = dq + dt*ddq;
  qNext = q + dt*dqNext;
  qNext = Clip[qNext, {-qLimit, qLimit}];

  xNext = {qNext, dqNext};

  cost = qWeight*qNext^2 + dqWeight*dqNext^2 + rWeight*uClip^2;
  reward = -cost;

  done = Norm[xNext] > 10.^3;

  << Return as association >>
  << Association[
    "NextState" -> xNext,
    "Reward" -> reward,
    "Done" -> done
  ]
 ]

(* Example usage *)
dt = 0.01;
x0 = {0.1, 0.0};
u0 = 1.0;

params = {dt, Pi, 10.0, 1.0, 0.1, 0.01}; (* dt, qLimit, uLimit, qW, dqW, rW *)
result = robotStep[x0, u0, Sequence @@ params];
result["NextState"]
result["Reward"]
result["Done"]
      

Mathematica can also be used to symbolically compute the Bellman optimality equation for low-dimensional problems or to check contraction properties analytically for specific choices of dynamics and reward, making it a powerful companion to the theoretical material of this lesson.

Problem 1 (Markov Property from Manipulator Dynamics). Let \( \mathbf{x}_t = (\mathbf{q}_t, \dot{\mathbf{q}}_t) \) and suppose the discrete-time dynamics are given by the Euler update derived from the rigid-body equations. Show that \( \{\mathbf{x}_t\}_{t \geq 0} \) is a Markov process when the control at time \( t \) is \( \mathbf{u}_t \) and the disturbance \( \boldsymbol{\eta}_t \) is independent and identically distributed.

Solution. By construction,

\[ \mathbf{x}_{t+1} = f(\mathbf{x}_t, \mathbf{u}_t) + \boldsymbol{\eta}_t. \]

For any measurable set \( A \subseteq \mathcal{X} \),

\[ \mathbb{P}(\mathbf{x}_{t+1} \in A \mid \mathbf{x}_0, \mathbf{u}_0, \ldots, \mathbf{x}_t, \mathbf{u}_t) = \mathbb{P}\big(f(\mathbf{x}_t,\mathbf{u}_t)+\boldsymbol{\eta}_t \in A \mid \mathbf{x}_t, \mathbf{u}_t \big), \]

because \( \boldsymbol{\eta}_t \) is independent of the past and the only dependence on history in the right-hand side is through \( \mathbf{x}_t \) and \( \mathbf{u}_t \). Therefore,

\[ \mathbb{P}(\mathbf{x}_{t+1} \in A \mid \mathbf{x}_0,\mathbf{u}_0,\ldots,\mathbf{x}_t,\mathbf{u}_t) = \mathbb{P}(\mathbf{x}_{t+1} \in A \mid \mathbf{x}_t,\mathbf{u}_t), \]

which is exactly the Markov property.

Problem 2 (Robot LQR as an MDP). Consider the linear stochastic system

\[ \mathbf{x}_{t+1} = \mathbf{A}\mathbf{x}_t + \mathbf{B}\mathbf{u}_t + \boldsymbol{\xi}_t, \quad \boldsymbol{\xi}_t \sim \mathcal{N}(\mathbf{0},\mathbf{\Sigma}), \]

with quadratic cost

\[ c(\mathbf{x}_t, \mathbf{u}_t) = \mathbf{x}_t^{\top}\mathbf{Q}\mathbf{x}_t + \mathbf{u}_t^{\top}\mathbf{R}\mathbf{u}_t, \]

where \( \mathbf{Q} \) is positive semidefinite and \( \mathbf{R} \) is positive definite. Identify the corresponding MDP components \( (\mathcal{X},\mathcal{U},p,r,\gamma) \) when \( \gamma \in (0,1) \), and argue that \( J^{\pi}(\mathbf{x}_0) \) is finite for all stationary policies \( \pi \).

Solution. We take

• \( \mathcal{X} = \mathbb{R}^{d_x} \), \( \mathcal{U} = \mathbb{R}^{d_u} \).
• \( p(\cdot \mid \mathbf{x},\mathbf{u}) \) is Gaussian with mean \( \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{u} \) and covariance \( \mathbf{\Sigma} \).
• \( r(\mathbf{x},\mathbf{u}) = -c(\mathbf{x},\mathbf{u}) \).
• \( \gamma \in (0,1) \) as given.

Because \( c(\mathbf{x},\mathbf{u}) \geq 0 \) and is quadratic, it grows at most polynomially in \( \lVert \mathbf{x} \rVert \) and \( \lVert \mathbf{u} \rVert \). For any stationary policy with at most linear growth in \( \mathbf{x} \) (e.g. linear feedback), standard results on discounted linear systems imply that the state moments \( \mathbb{E}[\lVert \mathbf{x}_t \rVert^2] \) grow at most geometrically, while \( \gamma^{t} \to 0 \). Thus the discounted sum

\[ J^{\pi}(\mathbf{x}_0) = \mathbb{E}\!\left[\sum_{t=0}^{\infty} \gamma^t c(\mathbf{x}_t,\mathbf{u}_t)\right] \]

converges, and hence \( V^{\pi} \) and \( Q^{\pi} \) are finite everywhere.

Problem 3 (Contraction of the Optimal Bellman Operator). Let \( T \) be the optimal Bellman operator defined in Section 5. Assume rewards are bounded: there exists \( R_{\max} < \infty \) such that \( \lvert r(\mathbf{x},\mathbf{u}) \rvert \leq R_{\max} \) for all \( (\mathbf{x},\mathbf{u}) \). Show that for all bounded functions \( V_1, V_2 \),

\[ \lVert T V_1 - T V_2 \rVert_{\infty} \leq \gamma \lVert V_1 - V_2 \rVert_{\infty}. \]

Solution. For any \( \mathbf{x} \),

\[ \begin{aligned} (T V_1)(\mathbf{x}) - (T V_2)(\mathbf{x}) &= \sup_{\mathbf{u}} \left[ r(\mathbf{x},\mathbf{u}) + \gamma \int p(\mathrm{d}\mathbf{x}'\mid \mathbf{x},\mathbf{u}) V_1(\mathbf{x}') \right] \\ &\quad - \sup_{\mathbf{u}} \left[ r(\mathbf{x},\mathbf{u}) + \gamma \int p(\mathrm{d}\mathbf{x}'\mid \mathbf{x},\mathbf{u}) V_2(\mathbf{x}') \right]. \end{aligned} \]

Using the inequality \( \sup_{\mathbf{u}} a_{\mathbf{u}} - \sup_{\mathbf{u}} b_{\mathbf{u}} \leq \sup_{\mathbf{u}} (a_{\mathbf{u}} - b_{\mathbf{u}}) \), we obtain

\[ (T V_1)(\mathbf{x}) - (T V_2)(\mathbf{x}) \leq \gamma \sup_{\mathbf{u}} \int p(\mathrm{d}\mathbf{x}' \mid \mathbf{x},\mathbf{u}) \big(V_1(\mathbf{x}')-V_2(\mathbf{x}')\big). \]

The integral is bounded by \( \lVert V_1 - V_2 \rVert_{\infty} \), since the transition kernel integrates to 1. Therefore

\[ (T V_1)(\mathbf{x}) - (T V_2)(\mathbf{x}) \leq \gamma \lVert V_1 - V_2 \rVert_{\infty}. \]

Swapping \( V_1 \) and \( V_2 \) gives the reverse inequality, hence taking the supremum over \( \mathbf{x} \) yields the desired contraction property.

Problem 4 (Equivalence of Cost and Reward Formulations). Let \( c(\mathbf{x},\mathbf{u}) \geq 0 \) be a running cost and define \( r(\mathbf{x},\mathbf{u}) = -c(\mathbf{x},\mathbf{u}) \). Show that for any policy \( \pi \) and initial state \( \mathbf{x}_0 \),

\[ J^{\pi}(\mathbf{x}_0) = \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^t c(\mathbf{x}_t,\mathbf{u}_t) \right] = - V^{\pi}(\mathbf{x}_0), \]

where \( V^{\pi} \) is the value function defined with rewards \( r \).

Solution. By definition of \( V^{\pi} \) in the reward convention,

\[ V^{\pi}(\mathbf{x}_0) = \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^t r(\mathbf{x}_t,\mathbf{u}_t) \right] = \mathbb{E}\!\left[ \sum_{t=0}^{\infty} \gamma^t (-c(\mathbf{x}_t,\mathbf{u}_t)) \right] = - J^{\pi}(\mathbf{x}_0). \]

Thus maximizing \( V^{\pi} \) over policies is exactly equivalent to minimizing \( J^{\pi} \), and we can switch freely between cost and reward viewpoints.

Problem 5 (State Design for Markovianity). Suppose a manipulator is controlled using only joint positions in the state: \( \tilde{\mathbf{x}}_t = \mathbf{q}_t \). Provide an example (in words or equations) showing why \( \tilde{\mathbf{x}}_t \) may violate the Markov property, and explain how augmenting the state restores it.

Solution. If the dynamics depend on joint velocities, e.g.

\[ \ddot{\mathbf{q}}_t = g(\mathbf{q}_t, \dot{\mathbf{q}}_t, \mathbf{u}_t), \]

then the next configuration \( \mathbf{q}_{t+1} \) depends on both \( \mathbf{q}_t \) and \( \dot{\mathbf{q}}_t \), even if the controller only observes \( \mathbf{q}_t \). The conditional distribution of \( \mathbf{q}_{t+1} \) given \( \tilde{\mathbf{x}}_t = \mathbf{q}_t \) and \( \mathbf{u}_t \) depends on the unobserved velocities, which in turn depend on the entire past. Thus \( \tilde{\mathbf{x}}_t \) is not Markov. Augmenting the state to \( \mathbf{x}_t = (\mathbf{q}_t,\dot{\mathbf{q}}_t) \) yields dynamics of the form \( \mathbf{x}_{t+1} = f(\mathbf{x}_t,\mathbf{u}_t) + \boldsymbol{\eta}_t \), which are Markov as shown in Problem 1.

In this lesson we embedded continuous robot control problems into the MDP framework. Starting from familiar manipulator and mobile robot dynamics, we defined continuous-state, continuous-action MDPs, specified transition kernels from the dynamics, and mapped standard quadratic tracking objectives into reward functions. We derived the Bellman equations in continuous spaces and showed their contraction properties, establishing the theoretical foundation on which RL algorithms operate. Finally, we implemented language-level abstractions in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, all reflecting the same mathematical structure.

Subsequent lessons in this chapter will build on this formalism to introduce policy gradient methods, actor–critic architectures, and off-policy algorithms that can handle the high-dimensional, continuous control problems typical in modern robotics.