ROS 2 and MATLAB Co-Simulation for Autonomous Mobile Robot Path Tracking Using NMPC
QUICK ANSWER
Learn how to set up ROS 2 and MATLAB co-simulation for mobile robot path tracking using Nonlinear Model Predictive Control (NMPC). Full code included.
What you'll learn
- 1. Kinematic Model of a Differential Drive Robot
- 2. NMPC Problem Formulation
- State and Actuator Constraints:
- 3. Architecture: Bridging ROS 2 and MATLAB
- 4. Implementation Steps in MATLAB
Tracking complex trajectories with non-holonomic mobile robots requires handling physical constraints such as actuator saturation, wheel slip, and sharp cornering. Standard controllers like Pure Pursuit or classic PID work well on gentle paths, but they degrade when robots operate near motor limits or encounter sudden turns.
Nonlinear Model Predictive Control (NMPC) solves a constrained optimization problem at each sampling step. It predicts future robot states over a receding time horizon and calculates optimal control commands while strictly enforcing physical limits.
This tutorial demonstrates how to build an end-to-end co-simulation pipeline between ROS 2 and MATLAB. You will configure the kinematic model of a differential drive robot, formulate the NMPC cost function, set up ROS 2 publishers and subscribers, and run closed-loop trajectory tracking.
1. Kinematic Model of a Differential Drive Robot
A standard unicycle-type differential-drive mobile robot operates in a 2D plane with state vector:
x(t) = [p_x, p_y, θ]T
where p_x and p_y represent the 2D position coordinates of the robot center, and θ is the orientation (yaw angle).
The control input vector is:
u(t) = [v, ω]T
where v is linear velocity and ω is angular velocity.
The continuous kinematic equations of motion are:
dp_x / dt = v · cos(θ)dp_y / dt = v · sin(θ)dθ / dt = ω
Using Euler forward discretization with sampling period Ts, the discrete-time state update equations are:
p_x(k+1) = p_x(k) + v(k) · cos(θ(k)) · Tsp_y(k+1) = p_y(k) + v(k) · sin(θ(k)) · Tsθ(k+1) = θ(k) + ω(k) · Ts
2. NMPC Problem Formulation
At each time step k, the NMPC controller solves the following optimal control problem across prediction horizon N:
min Σj=0N-1 [ ||xj - xref,j||Q2 + ||uj - uref,j||R2 ] + ||xN - xref,N||Q_N2
State and Actuator Constraints:
- Linear velocity limits:
0.0 ≤ v ≤ vmax(forward motion only) - Angular velocity limits:
-ωmax ≤ ω ≤ ωmax - Acceleration limits:
|Δv| ≤ amax · Ts - Steering rate limits:
|Δω| ≤ αmax · Ts
The state weighting matrix Q = diag(q_x, q_y, q_θ) prioritizes Cartesian path accuracy, while the input penalty matrix R = diag(r_v, r_ω) damps aggressive control actions to prevent motor saturation and wheel slip.
3. Architecture: Bridging ROS 2 and MATLAB
| ROS 2 / Gazebo Node | MATLAB Workspace (Controller) |
|---|---|
Publishes /odom topic (Position, Yaw, Velocity) |
Subscriber callback updates current estimated robot state |
Subscribes to /cmd_vel topic |
NMPC solver computes optimal [v, ω] and publishes command |
Communication between MATLAB and ROS 2 runs over native DDS (Data Distribution Service). MATLAB acts as a node within the ROS 2 graph, reading robot odometry, executing the sequential quadratic programming (SQP) solver, and publishing velocity commands directly back to the simulator or physical robot.
4. Implementation Steps in MATLAB
Step 1: Initialize the ROS 2 Node and Topics
Create the MATLAB ROS 2 node, subscribe to the odometry topic, and configure the velocity publisher:
ros2node_matlab = ros2node("/matlab_nmpc_node");
cmd_pub = ros2publisher(ros2node_matlab, "/cmd_vel", "geometry_msgs/Twist");
cmd_msg = ros2message(cmd_pub);
odom_sub = ros2subscriber(ros2node_matlab, "/odom", "nav_msgs/Odometry");
Step 2: Define Prediction Horizon and Weights
Set a sampling time of 0.1 s and a prediction horizon of N = 15 steps. This balances tracking fidelity with solver computation time:
Ts = 0.1; % Sampling period (10 Hz)
N = 15; % Prediction horizon
Q = diag([10.0, 10.0, 1.5]); % State weights [x, y, theta]
R = diag([0.5, 0.2]); % Control weights [v, w]
lb = repmat([0.0; -1.2], N, 1);
ub = repmat([1.0; 1.2], N, 1);
Step 3: Execute the Real-Time Control Loop
At each iteration, parse the robot pose, solve the NMPC cost function using fmincon, apply the first control input, and publish to /cmd_vel.
5. Common Pitfalls and Engineering Solutions
1. Heading Angle Discontinuities
Heading errors can jump from +π to -π. Always normalize angular error using atan2(sin(err), cos(err)) before passing states to the solver. Without angle wrapping protection, the robot will spin erratically near the branch cut.
2. Solver Computation Time vs. Jitter
Keep the prediction horizon N between 10 and 20 steps. If the solver takes longer than the sampling period (Ts = 0.1 s), the robot will lag behind the reference trajectory. For real-time hardware execution, generate standalone C++ code using MATLAB Coder.
3. Obstacle Avoidance Feasibility
Hard distance constraints can cause the optimization solver to return infeasible status codes if the robot gets trapped. Use soft barrier functions with slack variables so the solver returns a safe braking command rather than terminating.
6. Project Support and Custom Engineering Services
Need assistance with advanced robotics controllers, custom vehicle kinematics, or ROS 2 hardware integration?
- For Students and Academic Researchers: We provide complete simulation models, thesis mentoring, and reproducible benchmark code.
- For Engineering Teams and Companies: We design custom NMPC path planners, Embedded Coder deployment pipelines, and Hardware-in-the-Loop (HIL) test rigs.
Contact the MatlabSolutions engineering team to discuss your project requirements or request custom MATLAB/Simulink source code.
Model Setup
%% NMPC Path Tracking with ROS 2 Co-Simulation
% Author: MatlabSolutions (https://matlabsolutions.com)
% Description: Closed-loop NMPC trajectory tracking for a differential drive
% robot communicating via ROS 2.
clear; clc; close all;
%% 1. Simulation and Vehicle Parameters
Ts = 0.1; % Sample time (seconds)
N = 15; % Prediction horizon
v_max = 1.0; % Max linear velocity (m/s)
v_min = 0.0; % Min linear velocity (forward only)
w_max = 1.2; % Max angular velocity (rad/s)
w_min = -1.2;
%% 2. Generate Reference Trajectory (Figure-8)
t_total = 40;
t_vec = 0:Ts:t_total;
num_steps = length(t_vec);
x_ref = 3 * sin(0.2 * t_vec);
y_ref = 3 * sin(0.4 * t_vec);
theta_ref = atan2(gradient(y_ref), gradient(x_ref));
ref_traj = [x_ref; y_ref; theta_ref];
%% 3. Configure ROS 2 Node and Topics
disp('Initializing ROS 2 Node...');
try
ros2node_matlab = ros2node("/matlab_nmpc_node");
catch
disp('Node already initialized or using default domain.');
end
% Publisher for velocity commands
cmd_pub = ros2publisher(ros2node_matlab, "/cmd_vel", "geometry_msgs/Twist");
cmd_msg = ros2message(cmd_pub);
% Subscriber for robot odometry
odom_sub = ros2subscriber(ros2node_matlab, "/odom", "nav_msgs/Odometry");
%% 4. NMPC Weight Matrices
Q = diag([10.0, 10.0, 1.5]); % State tracking weights (x, y, theta)
R = diag([0.5, 0.2]); % Input penalty (v, w)
%% 5. Simulation Initial Conditions
current_state = [x_ref(1) + 0.3; y_ref(1) - 0.2; theta_ref(1) + 0.1];
u_prev = [0.0; 0.0];
state_history = zeros(3, num_steps);
input_history = zeros(2, num_steps);
disp('Starting Co-Simulation Loop...');
%% 6. Main Control Loop
for k = 1:(num_steps - N)
% Read current state from ROS 2 odometry
[odom_msg, status] = receive(odom_sub, 0.05);
if status
pos = odom_msg.pose.pose.position;
orient = odom_msg.pose.pose.orientation;
% Convert quaternion to yaw angle
siny_cosp = 2 * (orient.w * orient.z + orient.x * orient.y);
cosy_cosp = 1 - 2 * (orient.y * orient.y + orient.z * orient.z);
yaw = atan2(siny_cosp, cosy_cosp);
current_state = [pos.x; pos.y; yaw];
end
state_history(:, k) = current_state;
% Extract current reference window
ref_window = ref_traj(:, k:(k + N - 1));
% Solve NMPC Optimization Step using fmincon
u_init = repmat(u_prev, N, 1);
lb = repmat([v_min; w_min], N, 1);
ub = repmat([v_max; w_max], N, 1);
options = optimoptions('fmincon', ...
'Display', 'off', ...
'Algorithm', 'sqp', ...
'MaxIterations', 25);
cost_func = @(u_seq) nmpc_cost(u_seq, current_state, ref_window, Q, R, Ts, N);
u_optimal = fmincon(cost_func, u_init, [], [], [], [], lb, ub, [], options);
% Extract first control action
u_apply = u_optimal(1:2);
input_history(:, k) = u_apply;
u_prev = u_apply;
% Publish to ROS 2 /cmd_vel
cmd_msg.linear.x = double(u_apply(1));
cmd_msg.angular.z = double(u_apply(2));
send(cmd_pub, cmd_msg);
% Discrete state propagation for internal validation
current_state(1) = current_state(1) + u_apply(1) * cos(current_state(3)) * Ts;
current_state(2) = current_state(2) + u_apply(1) * sin(current_state(3)) * Ts;
current_state(3) = current_state(3) + u_apply(2) * Ts;
% Keep heading bounded within [-pi, pi]
current_state(3) = atan2(sin(current_state(3)), cos(current_state(3)));
end
% Stop robot at the end of the run
cmd_msg.linear.x = 0.0;
cmd_msg.angular.z = 0.0;
send(cmd_pub, cmd_msg);
disp('Path tracking complete.');
%% 7. Plotting Trajectory Tracking Performance
figure('Name', 'NMPC Trajectory Tracking Results');
plot(x_ref, y_ref, 'k--', 'LineWidth', 1.5); hold on;
plot(state_history(1, 1:num_steps-N), state_history(2, 1:num_steps-N), 'b-', 'LineWidth', 2);
grid on; axis equal;
xlabel('X Position (m)');
ylabel('Y Position (m)');
legend('Reference Path', 'Robot Actual Path', 'Location', 'best');
title('Autonomous Mobile Robot NMPC Path Tracking');
%% Local Function: NMPC Cost Function
function J = nmpc_cost(u_seq, x0, ref_window, Q, R, Ts, N)
J = 0;
x = x0;
for i = 1:N
u = u_seq((2*i - 1):(2*i));
x_ref = ref_window(:, i);
% Predict next state using differential-drive kinematics
x(1) = x(1) + u(1) * cos(x(3)) * Ts;
x(2) = x(2) + u(1) * sin(x(3)) * Ts;
x(3) = x(3) + u(2) * Ts;
% Angular error correction
err = x - x_ref;
err(3) = atan2(sin(err(3)), cos(err(3)));
% Accumulate quadratic stage cost
J = J + err' * Q * err + u' * R * u;
end
end