# Convergent NMPC-based Reinforcement Learning Using Deep Expected Sarsa and Nonlinear Temporal Difference Learning

Amine Salaje, Thomas Chevet, Nicolas Langlois

► To cite this version:

Amine Salaje, Thomas Chevet, Nicolas Langlois. Convergent NMPC-based Reinforcement Learning Using Deep Expected Sarsa and Nonlinear Temporal Difference Learning. European Control Conference, Jun 2025, Thessaloniki, Greece. pp.2644-2649, 10.23919/ECC65951.2025.11187042. hal-04983821

HAL Id: hal-04983821

https://hal.science/hal-04983821v1

Submitted on 22 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Convergent NMPC-based Reinforcement Learning Using Deep Expected Sarsa and Nonlinear Temporal Difference Learning

Amine Salaje, Thomas Chevet, and Nicolas Langlois

**Abstract**— In this paper, we present a learning-based nonlinear model predictive controller (NMPC) using an original reinforcement learning (RL) method to learn the optimal weights of the NMPC scheme, for which two methods are proposed. Firstly, the controller is used as the current action-value function of a deep Expected Sarsa where the subsequent action-value function, usually obtained with a secondary NMPC, is approximated with a neural network (NN). With respect to existing methods, we add to the NN's input the current value of the NMPC's learned parameters so that the network is able to approximate the action-value function and stabilize the learning performance. Additionally, with the use of the NN, the real-time computational burden is approximately halved without affecting the closed-loop performance. Secondly, we combine gradient temporal difference methods with a parametrized NMPC as a function approximator of the Expected Sarsa RL method to overcome the potential parameters' divergence and instability issues when nonlinearities are present in the function approximation. The simulation results show that the proposed approach converges to a locally optimal solution without instability problems.

## I. INTRODUCTION

Nonlinear model predictive control (NMPC) has proven effective in addressing complex control challenges [1], particularly in systems with nonlinear dynamics and constraints, such as mobile robotics [2]. It generates control signals by solving a constrained optimization problem – the optimal control problem (OCP) – at each time step. However, fine-tuning the parameters of the OCP remains a challenge, with ongoing research seeking reliable tuning methods [3]. Reinforcement learning (RL) offers a promising solution as it enables an agent to optimize these parameters by learning optimal value functions and policies within the scope of a Markov decision process (MDP) [4].

In problems with large MDPs or high-dimensional state spaces, modern RL methods often use deep neural networks (DNN) to approximate value functions or policies, allowing them to handle the complexity of nonlinear systems. However, DNN-based RL approaches often face difficulties in closed-loop stability analysis and can be challenging for formal verification and constraint management [5], thus raising concerns about ensuring safety in control systems. To address these issues, NMPC-based RL has been proposed and justified in [5] as a promising solution. By using the OCP as a function

approximator for the optimal policy in RL, this approach naturally satisfies state and input constraints while ensuring safety requirements.

In [5], the NMPC parameters are tuned to enhance closed-loop performance using Q-learning and policy gradient RL methods. Since then, various studies have deployed Q-learning with MPC as a function approximator to realize this task, with notable contributions found in [6]–[10]. However, one of the current challenges of these approaches is the significant computational demand for real-time implementation [11]. Indeed, in the context of temporal difference (TD) learning methods, it is required to solve at least two constrained nonlinear optimization problems at each sampling instant – one for estimating the current action-value function and another for approximating the subsequent action-value function.

To deal with this issue, a new data-driven approach was proposed in [12], where a parametrized NMPC has been used as an approximator for a deep double Expected Sarsa (ES) RL method. To estimate the current action-value function, a parametrized OCP is solved at each time step. However, the subsequent action-value function is approximated using a target neural network (NN) instead of solving a second OCP at the next time step. The NN is trained with inputs and outputs of the primary NMPC obtained at previous sampling times. This approach shows a clear advantage in reducing online computation time compared to [5]. However, as the second NMPC of [5] is replaced by a NN in [12], the simulation results show instability and potential divergence. As a first contribution, we propose in the present paper to tackle this problem by adding the parametrization vector of the NMPC at previous time steps to the training input of the NN, leading to a clear improvement in the convergence stability and control performance compared with [12].

While simple to use, conventional TD learning methods such as Q-learning lack formal guarantees for the closed-loop optimality of the learned policy [4]. Furthermore, they lack convergence guarantees when using nonlinear function approximation, which may cause potential divergence of the approximator parameters. To overcome these challenges, the authors of [13] proposed gradient temporal difference (GTD) for off-policy Q-learning methods as an effective solution, a method ensuring almost-sure convergence for any finite MDP and any smooth value function approximator to a locally optimal solution or equilibrium point. As a second contribution, we combine in this paper GTD methods with parametrized NMPC as a function approximator in the Expected Sarsa (ES) RL method to overcome the potential parameters divergence and instability issues. To the best of

<footnote>*This work was supported by ANR and Région Normandie through the HAISCoDe (ANR-20-THIA-0021) and PAMAP projects. For this reason and the purpose of Open Access, the authors have applied a CC BY public copyright licence to any Author Accepted Manuscript (AAM) version arising from this submission.</footnote>

<footnote>The authors are with Université de Rouen, ESIGELEC, IRSEEM, 76000 Rouen, France amine.salaje@groupe-esigelec.org, {thomas.chevet, nicolas.langlois}@esigelec.fr</footnote>
the authors' knowledge, this is the first work where an action-value method is used to tune the NMPC approximator using GTD methods. The simulation results show that the proposed approach converges to a locally optimal solution without instability problems compared to previous methods.

The remainder of this article is organized as follows. Section II exposes prerequisites on predictive control and reinforcement learning. Then, Section III presents our contributions. Finally, Section IV showcases the efficiency of our methods on numerical simulation examples.

## II. PREREQUISITES ON LEARNING AND CONTROL

We start by giving some context on the considered control strategy and background on the RL methods we use.

### A. Markov decision process for RL

Consider the tuple $(S, A, P, \ell, \rho, \gamma)$ which defines an MDP with continuous state space $S \subseteq \mathbb{R}^n$ and action space $A \subseteq \mathbb{R}^m$. The state transition dynamics have the underlying probability law $P$ such that $s_{k+1} \sim P(\cdot|a_k, s_k)$, with $s_{k+1}, s_k \in S$ and $a_k \in A$. In a control approach, these state transition dynamics can be expressed as $s_{k+1} = f(a_k, s_k)$ where $f$ is the discretized real system dynamics. The remaining elements of the tuple are the reward $\ell$, or stage cost in optimal control, which is associated to each MDP transition, the initial state distribution $\rho$, and the discount factor $\gamma \in (0, 1)$ determining the importance of future rewards.

The goal of RL is to find a policy, denoted by $\pi: S \to A$, minimizing a closed-loop performance objective [4]. Consider then a deterministic policy delivering the control input $a_k = \pi(s_k)$ resulting in state distribution $\rho^\pi$. Given an objective $J(\pi) = \mathbb{E}_{s \sim \rho^\pi} \left[ \sum_{k=0}^K \gamma^k \ell(s_k, a_k) | a_k = \pi(s_k) \right]$ on the expected discounted cumulative cost, with $K \in \mathbb{N} \cup \{+\infty\}$ a (possibly infinite) horizon, the RL problem aims to find the optimal policy $\pi^*$ satisfying $\pi^* = \arg\min J(\pi)$.

The optimal action-value function $Q^*(s_k, a_k)$ and value function $V^*(s_k)$ associated to the optimal policy $\pi^*$ can be determined using the Bellman equations

$$ 
\begin{aligned}
Q^*(s_k, a_k) &= \ell(s_k, a_k) + \gamma \mathbb{E} [V^*(s_{k+1}) | s_k, a_k], \\
V^*(s_k) &= Q^*(s_k, \pi^*(s_k)) = \min Q^*(s_k, a_k).
\end{aligned} 
\quad (1) 
$$

### B. NMPC as a function approximator

As in [5], we use a parametric optimization problem as a function approximator for the reinforcement learning. We therefore consider the parametrized NMPC

$$ 
\begin{aligned}
\min_{\mathbf{u}_i, \mathbf{x}_{i, \sigma_i}, \forall i \in \overline{0, N}} & \sum_{i=0}^{N-1} \gamma^i (L(\mathbf{x}_i, \mathbf{u}_i, \vartheta) + \boldsymbol{\omega}^\top \boldsymbol{\sigma}_i) \\
& + \gamma^N (V(\mathbf{x}_N, \vartheta_f) + \boldsymbol{\omega}_f^\top \boldsymbol{\sigma}_N) \quad (2a)
\end{aligned} 
$$

$$ \text{s.t.} \quad \mathbf{x}_0 = s_k, \quad (2b) $$

$$ \mathbf{x}_{i+1} = f(\mathbf{u}_i, \mathbf{x}_i), \quad \forall i \in \overline{1, N}, \quad (2c) $$

$$ \mathbf{x}_i \in S, \quad \forall i \in \overline{1, N}, \quad (2d) $$

$$ \mathbf{u}_i \in A, \quad \forall i \in \overline{0, N-1}, \quad (2e) $$

$$ g(\mathbf{x}_i) + \theta_c \le \boldsymbol{\sigma}_i, \quad \boldsymbol{\sigma}_i \ge 0, \forall i \in \overline{0, N}. \quad (2f) $$

In this formulation, $\forall i \in \overline{0, N}$, the variables $\mathbf{u}_i, \mathbf{x}_i$, and $\boldsymbol{\sigma}_i$ represent the primal decision variables, $L$ the NMPC stage cost and $V$ the terminal cost. The prediction horizon $N \in \mathbb{N}$ may be shorter than $K$ from the performance measure $J(\pi)$, and $\gamma$ is the discount factor referenced in Paragraph II-A. The function $f$ denotes the deterministic dynamic model driven by the NMPC scheme (2) as in Paragraph II-A, while the function $g$ imposes additional constraints on the state at each time step over the prediction horizon $N$. The variables $\boldsymbol{\sigma}_i$, for $i \in \overline{0, N}$, are slack variables, introduced to relax state constraints, with their violation penalized by the weights $\boldsymbol{\omega}$ and $\boldsymbol{\omega}_f$, preventing the potential infeasibility of the NMPC due to the additional constraints. Finally, $\vartheta$ and $\vartheta_f$ are vectors of NMPC weights, usually selected a priori based, e.g., on the dynamics $f$ or the control objective for the system. In this paper, these weights are treated as parameter vectors to be adjusted through an RL algorithm.

Indeed, an RL algorithm can advantageously adjust the parameters of the NMPC cost function, model, and constraints in a way to enhance closed-loop performance. In theory, as stated in [5, Theorem 1], under specific assumptions and with a sufficiently expressive parametrization, the MPC scheme can approximate the optimal policy $\pi^*$.

When solving (2) at time $k \in \mathbb{N}$, we obtain the optimal control sequence $\mathbf{u}^*(s_k, \theta) = \{\mathbf{u}_0^*(s_k, \theta), \dots, \mathbf{u}_{N-1}^*(s_k, \theta)\}$, with $\theta = [\vartheta \ \vartheta_f]$ driving the system with dynamics (2c) towards its objective over the prediction horizon $N$.

### C. Expected Sarsa with NMPC as function approximator

ES is a TD control method used in model-free reinforcement learning [4]. To handle large MDPs, the ES algorithm employs function approximators parametrized by unknown parameters $\theta$ to approximate the optimal action-value function $Q^*$.

Under some assumptions given in [5], the NMPC scheme (2) approximates the optimal action-value function $Q^*$. Thus, for any $k \in \mathbb{N}$, we define

$$ Q_{\theta}(s_k, a_k) = (2), \quad (3) $$

an approximator of the optimal action-value function $Q^*$. Unlike Q-learning, when combining NMPC with the ES method for on-policy learning, the action-value function (3) does not need to be evaluated at the first action $\mathbf{u}_0 = \mathbf{a}_k$ [14] (the behavioral policy is the same as the target policy). Furthermore, we keep the same notation for $Q_\theta(s_k, a_k)$ (rather than $V(s_k)$) with respect to our comparison with [12], [14].

The policy $\pi_\theta(s_k)$ and the action $\mathbf{a}_k$ are obtained by solving (2), so that, for any $s_k \in S$,

$$ \pi_{\theta}(s_k) = \arg \min_{\mathbf{a}_k \in A} Q_{\theta}(s_k, \mathbf{a}_k) = \mathbf{u}_0^*(s_k, \theta), \quad (4) $$

$$ \mathbf{a}_k = \pi_{\theta}(s_k) + c_{\varepsilon}^k \varepsilon \quad (5) $$

satisfying the Bellman equation [4]. Here, $\varepsilon \sim \mathcal{N}(0, 1)$ represents a Gaussian term introducing weighted exploration of the continuous action space $A$. This exploration decreases exponentially during training due to the factor $c_\varepsilon \in (0, 1)$.
During the learning process, the action $\mathbf{a}_k$ defined by (5), is applied to the system at time step $k$, resulting in $\mathbf{s}_{k+1}$.

To simplify the notations, the time dependence is dropped for states $\mathbf{s} \in \mathcal{S}$ and actions $\mathbf{a} \in \mathcal{A}$. Then, $\mathbf{s}^+ \in \mathcal{S}$ denotes the state at the next time instant.

In order to find $Q^*$, we characterize the optimal parameters $\theta^*$ as those minimizing the least-squares problem

$$ \theta^* = \arg\min_{\theta} \mathbb{E}_{\pi_{\theta}} \left[ \frac{1}{2} (Q^*(\mathbf{s}, \mathbf{a}) - Q_{\theta}(\mathbf{s}, \mathbf{a}))^2 \right]. \quad (6) $$

Since the optimal $Q^*$ is usually unknown, we cannot solve (6) directly. To do so, we use here the continuous ES algorithm [14], [15].

The parametrized action-value function $Q_{\theta}$ is estimated using the NMPC scheme (2) parametrized by the vector $\theta$. At each sampling time, the parameters are updated based on the TD error

$$ \delta = \ell(\mathbf{s}, \mathbf{a}) + \gamma Q_{\theta}(\mathbf{s}^{+}, \mathbf{a}^{+}) - Q_{\theta}(\mathbf{s}, \mathbf{a}) \quad (7) $$

where $\mathbf{a}^+$ is obtained by applying $\mathbf{a}$, obtained as in (5) by solving the OCP (2) for state $\mathbf{s}$ (behavioral policy), to the controlled system and solving again the OCP (2) for the resulting state $\mathbf{s}^+$ (target policy). The parameter updates are determined by minimizing the mean-square TD error

$$ \min_{\theta} \mathbb{E}_{\pi_{\theta}} \left[ \frac{1}{2} \left( \ell(\mathbf{s}, \mathbf{a}) + \gamma Q_{\theta}(\mathbf{s}^{+}, \mathbf{a}^{+}) - Q_{\theta}(\mathbf{s}, \mathbf{a}) \right)^2 \right]. \quad (8) $$

We apply the semi gradient step to perform this minimization, thus giving the parameters' update rule

$$ \theta \leftarrow \theta + \alpha \delta \nabla_{\theta} Q_{\theta}(\mathbf{s}, \mathbf{a}) \quad (9) $$

where $\alpha > 0$ is the learning rate. By updating the unknown parameters, the action-value functions can converge towards its optimal value. The gradient $\nabla_{\theta}Q_{\theta}(\mathbf{s}, \mathbf{a})$ is obtained from sensitivity analysis. For more details the reader is referred to [5].

### III. PROPOSED METHODS

This section presents the main contributions of this paper, aiming to improve the overall performances of the learning algorithm. Our methods are then summarized in Algorithm 1.

#### A. Reducing computational complexity of NMPC-based ES

**1) Problem definition:** In [5] and [14], the parameter vector $\theta$ is tuned using Q-learning and an ES algorithm, respectively. In both cases, it is required to solve, at each learning step, a first OCP of the form (2) to obtain $Q_{\theta}(\mathbf{s}, \mathbf{a})$ as well as a second OCP to obtain the subsequent action-value function $Q_{\theta}(\mathbf{s}^{+}, \mathbf{a}^{+})$. These two action-value functions are then used to compute the TD error to be minimized. However, solving two constrained nonlinear optimization problems at each time step significantly increases the real-time computational complexity, resulting in longer training times. This issue worsens with large prediction horizons $N$, when the OCP (2) involves multiple constraints, etc. To address this problem, [12] replace the second OCP of the aforementioned

methods by a neural network approximator. However, this method may diverge due to an inaccurate approximation of the subsequent action-value function.

**2) NMPC-based off-policy regularized deep ES (RDES):** To better approximate the subsequent action-value function obtained with a second NMPC [14] using an NN, we propose to adjust the input data used to train the NN. The idea is to learn the action-value function $Q_{\theta}$ of (3) using a neural network $Q_{\varpi_{NN}}$ parametrized by the weights $\varpi_{NN}$ so that $Q_{\varpi_{NN}}(\mathbf{s}, \mathbf{a}, \theta) \approx Q_{\theta}(\mathbf{s}, \mathbf{a})$. To do so, the NN training input should be chosen carefully.

In order to train the NN, we store, at each learning step, the tuple $\{\mathbf{s}, \mathbf{a}, \theta, Q_{\theta}(\mathbf{s}, \mathbf{a})\}$ into a replay buffer $\mathcal{D}$. The current state $\mathbf{s}$, the action $\mathbf{a}$ and the parametrization vector $\theta$ are the NN's inputs while $Q_{\theta}(\mathbf{s}, \mathbf{a})$ is the regression target. Then, a mini-batch $\mathcal{B} \subset \mathcal{D}$ of $n \in \mathbb{N}$ elements is randomly sampled and we update $\varpi_{NN}$ using the Adam algorithm [16] with a learning rate $\zeta$ to minimize the cost function

$$ C_{NN}(\theta) = \frac{1}{n} \sum_{\{\mathbf{s}, \eta, \vartheta, q\} \in \mathcal{B}} (Q_{\varpi_{NN}}(\mathbf{s}, \eta, \vartheta) - q)^2 \quad (10) $$

where $Q_{\varpi_{NN}}(\mathbf{s}, \eta, \vartheta)$ for $\{\mathbf{s}, \eta, \vartheta, q\} \in \mathcal{B} \subset \mathcal{D}$ is the output of the NN with parameters $\varpi_{NN}$ and inputs $\mathbf{s}, \eta$, and $\vartheta$, whereas $q$ is the regression target.

After the NN's update, the subsequent action-value function in (7) is predicted by passing to the NN the next state $\mathbf{s}^+$ and the new target policy defined as $\mathbf{a}^+ = \pi_{\theta}(\mathbf{s}^+) = \mathbf{u}_1^*(\mathbf{s}, \theta)$ the optimal value of $\mathbf{u}_1$, the second input vector over the prediction horizon when solving the primary NMPC (4). By doing so, only one constrained optimization problem needs to be solved at each time step as well as the NN optimization process. Finally, we define the new TD error to be minimized, at each learning step, by

$$ \varrho = \ell(\mathbf{s}, \mathbf{a}) + \gamma Q_{\varpi_{NN}}(\mathbf{s}^{+}, \mathbf{a}^{+}, \theta) - Q_{\theta}(\mathbf{s}, \mathbf{a}). \quad (11) $$

We note that even when replacing the subsequent action-value function by the NN approximator, the constraint satisfaction is always maintained since the policy is derived from the OCP scheme (2). Compared with [14], our contribution then consists in adjusting the input data used to train the NN. We incorporate the current value of the parametrization vector $\theta$ of the action-value function defined in (3) to the NN's input. This is based on the assumption that the approximated subsequent action-value function should depend on $\theta$, as the current action-value function takes vector $\theta$ as input at each time step. Thus, by adding the parameter vector $\theta$ to the training data, the NN gains more knowledge about the behavior of the primary NMPC (2). Therefore, stabilizing the learning performance and improving the approximation of the action-value function, which is usually computed by solving a second NMPC, under the assumption

$$ Q_{\varpi_{NN}}(\mathbf{s}, \mathbf{a}, \theta) \rightarrow Q_{\theta}(\mathbf{s}, \mathbf{a}). \quad (12) $$

#### B. Nonlinear temporal difference learning

**1) Problem definition:** Traditional TD methods like Q-learning and ES have been successfully applied with function
approximation in various fields. However, it is well-known
that these algorithms can become unstable, particularly when
nonlinear function approximation is used, which can lead to
parameter divergence in the approximator. To address this
issue, we build on the results from [13], [17], and our idea is
to use a GTD approach to address instability and convergence
issues when using NMPC as function approximator of RL
methods. Indeed, the GTD methods are proven to converge
almost surely for any MDP, even with nonlinear function
approximators, to a locally optimal solution or equilibrium
point.

2) **Convergent NMPC-based gradient ES (GES):** In TD learning control methods, most works use an objective function representing how closely the approximate action-value function satisfies the Bellman equations. The true action-value function satisfies the Bellman equation in model-free RL.

$$
\Gamma^{\pi} Q(\mathbf{s}, \mathbf{a}) = \mathbb{E}_{\pi} \left[ \ell(\mathbf{s}, \mathbf{a}) + \gamma Q(\mathbf{s}^{+}, \mathbf{a}^{+}) \right] \quad (13)
$$

where $\Gamma^\pi$ is known as the Bellman operator. In GTD approaches, the objective function that measures how closely the approximation $Q_\theta(\mathbf{s}, \mathbf{a})$ satisfies the Bellman equation (13) is known as the mean-square projected Bellman error (MSPBE), where the goal is to minimize

$$
\mathcal{J}(\theta) = \|\Pi(\Gamma^{\pi\theta} Q_{\theta}) - Q_{\theta}\|_{\rho^{\pi}}^{2} \quad (14)
$$

using approximate gradient decent, Π is a projection operator
which projects action-value functions into the space $\mathcal{F} =$
$\{Q_{\theta} : \theta \in \mathbb{R}^d\}$ w.r.t. $\|\cdot\|_{\rho^{\pi}} : \Pi(\hat{Q}) = \arg \min_{f \in \mathcal{F}} \|\hat{Q} - f\|_{\rho^{\pi}}$.

Following [13], the projected objective can be written in
terms of expectations as

$$
\mathcal{J}(\theta) = \mathbb{E}[\delta\phi]^{\top}\mathbb{E}[\phi\phi^{\top}]^{-1}\mathbb{E}[\delta\phi], \quad (15)
$$

where $\delta$ is the TD error as in (7), and $\phi$ is the feature
vector used in approximating the action-value function. In
[13], the authors used an unrestricted linear form for the
action-value function approximation, specifically $Q_{\theta}(\mathbf{s}, \mathbf{a}) =
\phi^{\top}\theta$. However, in this paper, we use a nonlinear function
approximation defined by the OCP scheme (2). Therefore,
we represent $\phi$ as the gradient of the action-value function,
expressed as $\phi = \nabla_{\theta}Q_{\theta}(\mathbf{s}, \mathbf{a})$.

Finally, by applying the negative gradient on the MSPBE
objective function, and introducing a new set of weights $\mathbf{w}$
making use of the weight doubling trick [13], [17], [18], we
get the update rules

$$
\begin{align}
\theta &\leftarrow \theta + \alpha [\delta\phi - \gamma (\mathbf{w}^{\top}\phi)\phi^{+}], \nonumber \\
\mathbf{w} &\leftarrow \mathbf{w} + \beta [\delta - \phi^{\top}\mathbf{w}]\phi, \tag{16}
\end{align}
$$

where $\phi^+$ is the subsequent action-value features
$\nabla_{\theta}Q_{\theta}(\mathbf{s}^{+}, \mathbf{a}^{+})$. Under some assumptions [13] on the
positive step sizes $\alpha, \beta > 0$, it follows that $\theta$ converges
almost surely to a stable equilibrium point.

<table><caption>**Algorithm 1:** Learning algorithm.</caption>
  <tbody>
    <tr>
      <td>1</td>
      <td>Input: <b>&theta;</b>, <b>w</b>, <i>&omega;</i><sub>NN</sub>, dynamical system <i>f</i>, NMPC (2)</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Output: <i>&theta;</i><sup>*</sup></td>
    </tr>
    <tr>
      <td>3</td>
      <td>Initialization: <i>D</i> &larr; &empty;</td>
    </tr>
    <tr>
      <td>4</td>
      <td>for <i>k</i> &larr; 1 to <i>N</i><sub>episode</sub> do</td>
    </tr>
    <tr>
      <td>5</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;<i>s</i> &larr; <i>s</i><sub>0</sub></td>
    </tr>
    <tr>
      <td>6</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;for <i>t</i> &larr; 0 to <i>K</i> do</td>
    </tr>
    <tr>
      <td>7</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Solve (2) in <i>s</i> to obtain <b>a</b>, and <i>Q</i><sub>&theta;</sub>(<b>s</b>, <b>a</b>)</td>
    </tr>
    <tr>
      <td>8</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Get the action <b>a</b> from (5)</td>
    </tr>
    <tr>
      <td>9</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Observe <b>s</b><sup>+</sup> = <i>f</i>(<b>a</b>, <i>s</i>) and get the reward <i>l</i></td>
    </tr>
    <tr>
      <td>10</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if RDES used then</td>
    </tr>
    <tr>
      <td>11</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Store <i>D</i> &larr; <i>D</i> &cup; {<b>s</b>, <b>a</b>, <i>&theta;</i>, <i>Q</i><sub>&theta;</sub>(<b>s</b>, <b>a</b>)}</td>
    </tr>
    <tr>
      <td>12</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;if Card(<i>D</i>) &ge; <i>n</i> then</td>
    </tr>
    <tr>
      <td>13</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Randomly sample mini-batch of <i>n</i></td>
    </tr>
    <tr>
      <td>14</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;elements of <i>D</i> and get <i>&omega;</i><sub>NN</sub> (10)</td>
    </tr>
    <tr>
      <td>15</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the TD error <i>&delta;</i> with (11)</td>
    </tr>
    <tr>
      <td>16</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update parameters <i>&theta;</i> by (9)</td>
    </tr>
    <tr>
      <td>17</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;else if GES used then</td>
    </tr>
    <tr>
      <td>18</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Solve (2) in <b>s</b><sup>+</sup> to obtain <i>Q</i><sub>&theta;</sub>(<b>s</b><sup>+</sup>, <b>a</b><sup>+</sup>)</td>
    </tr>
    <tr>
      <td>19</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Compute the TD error <i>&delta;</i> with (7)</td>
    </tr>
    <tr>
      <td>20</td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Update parameters <i>&theta;</i> and <b>w</b> by (16)</td>
    </tr>
    <tr>
      <td></td>
      <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;<i>s</i> &larr; <b>s</b><sup>+</sup></td>
    </tr>
  </tbody>
</table>

IV. NUMERICAL APPLICATIONS

In this section we apply both our RL-based methods to
the case of setpoint tracking problem with static obstacle
avoidance for a diff-drive mobile robot. These simulations
aim to show the clear improvement in the methods we propose
compared with the methods from [14] and [12].

A. *Simulation model*

We consider a diff-drive robot with continuous-time dynamics [19]

$$
\dot{x} = v \cos \varphi, \quad (17a)
$$

$$
\dot{y} = v \sin \varphi, \qquad (17b)
$$

$$
\dot{\varphi} = \nu. \tag{17c}
$$

The robot's state is defined by its coordinates (x, y) in a 2D
plane and its orientation ϕ with inputs its linear and angular
speeds v and v. For training, we discretize the dynamics
(17) with the Runge-Kutta 4 method and a sampling period
T_s = 0.2 s, so that

$$
\mathbf{s}_{k+1} = f(\mathbf{u}_k, \mathbf{s}_k), \quad (18)
$$

where $\mathbf{s}_k^\top = [x_k \quad y_k \quad \varphi_k]$ and $\mathbf{u}_k^\top = [v_k \quad \nu_k]$.

The robot objective is to find the shortest path to reach an
equilibrium point (s<sup>ref</sup>, u<sup>ref</sup>) of the dynamics (18), starting
from a given initial position s<sub>0</sub>. Here s<sup>ref</sup> denotes the desired
target state vector. Along its trajectory, the robot must
maintain a desired safe distance from obstacles. To this end,
the robot is controlled using an NMPC as defined in (2). The
<table><caption>TABLE I
CONTROLLER AND ALGORITHM CONFIGURATION.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>α, β, ζ</td><td>10<sup>-7</sup>, 10<sup>-8</sup>, 10<sup>-2</sup></td></tr><tr><td>w<sub>init</sub></td><td>10<sup>-4</sup></td></tr><tr><td>Robot θ<sub>init</sub></td><td>[1 1 0.05 0.05 0.05 1 1 0.1 0.001]</td></tr><tr><td>Law of ρ, γ</td><td>N(0, 1), 0.97</td></tr><tr><td>N, K, N<sub>ep</sub></td><td>10, 129, 300</td></tr><tr><td>ω, ω<sub>f</sub></td><td>[100, 100, 100, 100]<sup>T</sup></td></tr><tr><td>Mini-batch size card(B)</td><td>128</td></tr><tr><td>Hidden layers, neurons per layer</td><td>2, 64</td></tr></tbody></table>

cost function (2a) is designed such that, for any $i \in [0, N-1]$,

$$L(\mathbf{s}_i, \mathbf{u}_i, \boldsymbol{\theta}) = \|\mathbf{s}_i - \mathbf{s}^{\text{ref}}\|_{\text{diag}(\boldsymbol{\vartheta}_x)}^2 + \|\mathbf{u}_i - \mathbf{u}^{\text{ref}}\|_{\text{diag}(\boldsymbol{\vartheta}_u)}^2,$$

$$\Gamma(\mathbf{s}_N, \boldsymbol{\theta}) = \|\mathbf{s}_N - \mathbf{s}^{\text{ref}}\|_{\text{diag}(\boldsymbol{\vartheta}_f)}^2$$

where $\text{diag}(\boldsymbol{\theta})$ is a diagonal matrix having for elements the components of $\boldsymbol{\theta}$, with $\boldsymbol{\vartheta}_x = [\theta_x \ \theta_y \ \theta_\varphi]$, $\boldsymbol{\vartheta}_u = [\theta_v \ \theta_\nu]$, and $\boldsymbol{\vartheta}_f = [\theta_{x_f} \ \theta_{y_f} \ \theta_{\varphi_f}]$, and $\|\mathbf{x}\|_Q^2 = \mathbf{x}^\top Q \mathbf{x}$.

The action space is $\mathcal{A} = [-0.6 \text{ m} \cdot \text{s}^{-1}, 0.6 \text{ m} \cdot \text{s}^{-1}] \times [-\pi/2, \pi/2]$. The state space is $\mathcal{S} = [-9 \text{ m}, 3 \text{ m}] \times [-9 \text{ m}, 3 \text{ m}] \times \mathbb{R}$. The robot has to stay at a given distance from 4 obstacles by satisfying the inequality constraint (2f), defined, for $n \in [1, 4]$, as

$$1 - 4 \frac{(x_i - x_{\text{obs},n})^2 + (y_i - y_{\text{obs},n})^2}{(d_{\text{rob}} + d_{\text{obs}})^2} + \theta_c \le 0, \quad (19)$$

where $(x_{\text{obs},n}, y_{\text{obs},n})$ is the center of the n-th obstacle, while $d_{\text{rob}} = 0.5 \text{ m}$ and $d_{\text{obs}} = 2 \text{ m}$ represent the diameters of a safety circle surrounding the robot and an obstacle, respectively. Parameter $\theta_c$ serves as a tuning variable, adjusting the intensity of the collision avoidance constraints. Overall, the parameters vector is $\boldsymbol{\theta} = [\theta_x \ \theta_y \ \theta_\varphi \ \theta_v \ \theta_\nu \ \theta_{x_f} \ \theta_{y_f} \ \theta_{\varphi_f} \ \theta_c]$.

To train the RL agent, the robot's initial state is set as $\mathbf{s}_0^\top = [-2.5 \text{ m} \ 1.5 \text{ m} \ 0 \text{ rad}]$, and the target state is defined as $\mathbf{s}^{\text{ref}} = [8.5 \text{ m} \ 2 \text{ m} \ \pi/2]^\top$, with the reference input $\mathbf{u}^{\text{ref}} = \mathbf{0}$. The vectors $\mathbf{s}_0$ and $\mathbf{s}^{\text{ref}}$ remain fixed throughout all learning episodes. The RL stage cost $\ell$ used to evaluate the RL performance defined in Paragraph II is

$$\ell = \begin{cases} \|\mathbf{s}_k - \mathbf{s}^{\text{ref}}\|_Q^2 + \|\mathbf{u}_k\|_2 & \text{if } \|\mathbf{s}_k - \mathbf{s}^{\text{ref}}\|_2 < d_t \\ \|\mathbf{s}_k - \mathbf{s}^{\text{ref}}\|_Q^2 + \Upsilon(\mathbf{s}_k) & \text{if } \|\mathbf{s}_k - \mathbf{s}^{\text{ref}}\|_2 \ge d_t \end{cases}, \quad (20)$$

where $Q = \text{diag}(1, 1, 0.1)$ is the objective weighting matrix, $d_t = 1.5 \text{ m}$ is a constant distance threshold, and $\Upsilon(\mathbf{s}_k) = \sum_{n=1}^4 c \max\{0, \Xi(\mathbf{s}_k) + d_o\}$ is the obstacle penalty collision where $\Xi$ is as defined as in (19), and $d_o = 0.35 \text{ m}$ is a desired safe distance between the robot and the obstacle, and $c = 50$ is a penalty weight. If a collision occurs, $\Upsilon(\mathbf{s}_k)$ adds a positive penalty to the RL stage cost $\ell$.

## B. Simulation results

To showcase the advantages of our algorithm, we train $\boldsymbol{\theta}$ for 300 episodes in four ways: one using the method from [14], one using the method from [12], and the last two using

<figure><img src="image_2.png" /><figcaption>Fig. 1. RL stage cost $\ell$ sum over learning episodes.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. RL stage cost $\ell$ sum over learning episodes.</figcaption></figure>

the proposed methods from paragraphs III-A and III-B. The numerical computation is performed using the Ipopt solver provided by the CasADi software framework [20] on a PC equipped of 16 GB of RAM.

For fair comparison, we compare the deep ES method from [12] with our RDES method from paragraph III-A. Indeed, for both these methods, only one NMPC must be solved to compute the current action-value function, whereas the subsequent action-value function is approximated by a neural network. Fig. 1 shows the sum of RL stage cost across training episodes for both approaches. In the blue graph, we can observe that the RL performance decreases faster, with the learning process remaining stable throughout training episodes, highlighting the effect of adding the parametrization vector into the NN's training input for better approximation of the subsequent action-value function. In contrast, the black dashed graph from [12] shows divergent RL performance, likely due to inaccuracies in approximating the subsequent action-value function.

To demonstrate the effectiveness of our second contribution from Paragraph III-B, we compare our proposed GES with the ES method from [14]. In both approaches, two NMPCs are solved to compute the TD error in (7). As we already stated, conventional temporal-difference learning methods can become unstable or even diverge when nonlinear function approximation is involved. In Fig. 2, we see that the blue graph reflects an unstable learning process as training episodes progress, attributed to the nonlinearities in the NMPC function approximator. In contrast, the red graph displays stable RL performance, with the GES method converging toward a stable equilibrium point. This demonstrates the effectiveness of GTD learning, which is a true gradient descent method,
<figure><img src="image_4.png" /><figcaption>Fig. 3. Robot control signals under the learned policies.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Robot states' trajectories under the learned policies.</figcaption></figure>

<table><caption>TABLE II
TRAINING TIME COMPARISON.</caption><thead><tr><td>Method</td><td>Training time (h)</td></tr></thead><tbody><tr><td>NMPC-based ES from [14]</td><td>≈ 18</td></tr><tr><td>NMPC-based GES from III-B</td><td>≈ 18</td></tr><tr><td>NMPC-based RDES from III-A</td><td>≈ 9</td></tr></tbody></table>

ensuring convergence.

Table II presents the real-time training duration of the three methods for 300 learning episodes. The deep ES method from [12] is excluded, as it diverges before reaching the 300th episode. As shown, our RDES method completes the training episodes in half the time required for both the method from [14] and our GES method.

## V. CONCLUSION

This paper proposes two reinforcement learning-based methods for tuning a nonlinear model predictive controller (NMPC). In the first method, a parametrized NMPC is used as the current action-value function in a regularized deep Expected Sarsa (RDES) algorithm, while a neural network (NN) with a novel feature representation serves as the subsequent action-value function. By incorporating the NMPC parameter vector into the NN's input training data, the learning process becomes more stable, and training time is reduced without compromising closed-loop performance. The second method introduces a combination of NMPC-based gradient Expected Sarsa (GES) with convergence guarantees. This approach stabilizes learning and ensures convergence to an equilibrium point, even in the presence of nonlinearities. Finally, we find that both methods achieve the same optimal policy. However, RDES reaches this outcome in half the training time, whereas GES provides convergence guarantees to a locally stable equilibrium.

In NMPC-based RL methods, adapting to dynamic environments, such as shifting obstacle positions during learning, depends on whether the changes are deterministic or stochastic. Deterministic changes allow pre-scheduled updates to the nonlinear programming (NLP) problem, adjusting constraints and costs to account for updated obstacle positions, while the RL method adapts its value function and policy accordingly. For stochastic changes, performance metrics or real-time detection of obstacle movements trigger dynamic updates to the NLP scheme, modifying the NLP structure to match the current environment. Meanwhile, for the RL component, one can consider adapting its learning rate to quickly update the value function and prioritizing recent data in the replay buffer to focus on the current environment's state. This ensures the policy remains effective and the system handles environmental shifts robustly, whether changes are predictable or random.

## REFERENCES

[1] L. Grüne and J. Pannek, *Nonlinear Model Predictive Control*, 2nd ed. Cham, Switzerland: Springer, 2017.

[2] T. P. Nascimento, C. E. T. Dórea, and L. M. G. Gonçalves, “Nonholonomic mobile robots’ trajectory tracking model predictive control: a survey,” *Robotica*, vol. 36, no. 5, pp. 676–696, 2018.

[3] Z. Zhang, O. Babayomi, T. Dragicevic, R. Heydari, C. Garcia, J. Rodriguez et al., “Advances and opportunities in the model predictive control of microgrids: Part i—primary layer,” *Int. J. Electr. Power Energy Syst.*, vol. 134, p. 107411, 2022.

[4] R. S. Sutton and A. G. Barto, *Reinforcement Learning: An introduction*. Cambridge, MA: MIT Press, 2018.

[5] S. Gros and M. Zanon, “Data-driven economic NMPC using reinforcement learning,” *IEEE Trans. Autom. Control*, vol. 65, no. 2, pp. 636–648, 2019.

[6] M. Zanon, S. Gros, and A. Bemporad, “Practical reinforcement learning of stabilizing economic MPC,” in *Proc. IEEE ECC*. IEEE, 2019, pp. 2258–2263.

[7] S. Gros, M. Zanon, and A. Bemporad, “Safe reinforcement learning via projection on a safe set: How to achieve optimality?” *IFAC-PapersOnLine*, vol. 53, no. 2, pp. 8076–8081, 2020.

[8] A. B. Martinsen, A. M. Lekkas, and S. Gros, “Combining system identification with reinforcement learning-based MPC,” *IFAC-PapersOnLine*, vol. 53, no. 2, pp. 8130–8135, 2020.

[9] A. B. Kordabad, H. N. Esfahani, A. M. Lekkas, and S. Gros, “Reinforcement learning based on scenario-tree MPC for ASVs,” in *Proc. ACC*. IEEE, 2021, pp. 1985–1990.

[10] H. N. Esfahani, A. B. Kordabad, and S. Gros, “Approximate robust NMPC using reinforcement learning,” in *Proc. IEEE ECC*. IEEE, 2021, pp. 132–137.

[11] A. B. Kordabad, D. Reinhardt, A. S. Anand, and S. Gros, “Reinforcement learning for MPC: Fundamentals and current challenges,” *IFAC-PapersOnLine*, vol. 56, no. 2, pp. 5773–5780, 2023.

[12] H. MoradiMaryamnegari, M. Frego, and A. Peer, “Data-driven model predictive control using deep double expected Sarsa,” in *Proc. CODIT*. IEEE, 2023, pp. 345–350.

[13] H. R. Maei, C. Szepesvári, S. Bhatnagar, and R. S. Sutton, “Toward off-policy learning control with function approximation,” in *Proc. 27th ICML*, 2010, pp. 719–726.

[14] H. Moradimaryamnegari, M. Frego, and A. Peer, “Model predictive control-based reinforcement learning using expected Sarsa,” *IEEE Access*, vol. 10, pp. 8177–8191, 2022.

[15] H. Van Seijen, H. Van Hasselt, S. Whiteson, and M. Wiering, “A theoretical and empirical analysis of expected Sarsa,” in *Proc. IEEE ADPRL*. IEEE, 2009, pp. 177–184.

[16] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” *arXiv preprint*, p. 1412.6980, 2014.

[17] H. Maei, C. Szepesvari, S. Bhatnagar, D. Precup, D. Silver, and R. S. Sutton, “Convergent temporal-difference learning with arbitrary smooth function approximation,” *Adv. Neural Inf. Process.*, vol. 22, 2009.
[18] R. S. Sutton, H. R. Maei, D. Precup, S. Bhatnagar, D. Silver, C. Szepesvári, and E. Wiewiora, "Fast gradient-descent methods for temporal-difference learning with linear function approximation," in *Proc. 26th ICML*, 2009, pp. 993–1000.

[19] M. Sani, B. Robu, and A. Hably, "Dynamic obstacles avoidance using nonlinear model predictive control," in *Proc. IECON*. IEEE, 2021, pp. 1924–1929.

[20] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, "CasADi – A software framework for nonlinear optimization and optimal control," *Math. Program. Comput.*, vol. 11, no. 1, pp. 1–36, 2019.