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

HAL
open science

# First-Order Continuation for Heterogeneous Multi-Agent Systems on a Ring: From PDE Modeling to Feedback Control

Sébastien Fueyo, Carlos Canudas de Wit

► To cite this version:

Sébastien Fueyo, Carlos Canudas de Wit. First-Order Continuation for Heterogeneous Multi-Agent Systems on a Ring: From PDE Modeling to Feedback Control. 2026. hal-05533076

HAL Id: hal-05533076

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

Preprint submitted on 2 Mar 2026

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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# First-Order Continuation for Heterogeneous Multi-Agent Systems on a Ring: From PDE Modeling to Feedback Control

S. Fueyo and C. Canudas de Wit

***Abstract***—In large-scale multi-agent systems, analyzing individual dynamics quickly becomes impractical, which motivates the use of continuous models to describe collective behavior. We introduce a continuation method on a ring that systematically derives macroscopic partial differential equations (PDEs) from microscopic dynamics, while allowing for heterogeneous interaction structures between agents. Control laws designed at the PDE level can then be consistently transferred to the underlying microscopic system, rendering the agents effectively autonomous under feedback. The approach is illustrated on the first-order Optimal Velocity traffic model, where the vehicle distribution is successfully steered toward a prescribed stationary density profile despite the presence of non-uniform interactions.

***Index Terms***—Multi-agent systems, heterogeneous interactions, PDE continuation, feedback control, traffic flow.

## I. INTRODUCTION

In many practical situations, the number of interacting particles in a system is so large that a direct study of their individual dynamics is infeasible. To overcome this difficulty, one typically resorts to continuous models that summarize the collective behavior of these numerous entities into a macroscopic description. Interestingly, discrete and continuous descriptions often mirror each other in meaningful ways. For instance, in traffic-flow theory, the microscopic and macroscopic models derived in [1] display closely aligned dynamical behaviors. A similar correspondence appears in synchronization phenomena: the stability properties of phase-locked solutions in networks of Kuramoto oscillators remain consistent whether one analyzes the discrete model or its continuous analogue [11]. A related perspective is adopted in [2], where the macroscopic traffic density is reconstructed from sparse discrete vehicle data, highlighting once again how continuous descriptions can be consistently inferred from underlying microscopic dynamics.

Although discretizing partial differential equations (PDE) into systems of ordinary differential equations (ODE) is a classical and widely adopted procedure, the reverse task of deriving a systematic PDE representation from an ODE network has received considerably less attention. In particular settings, such as oscillator networks or traffic-flow models [1], [11], PDEs have been obtained through ad-hoc derivations. Other

works propose methodologies for special classes of systems, yet these approaches remain far from offering a general and systematic analogue to standard discretization schemes. A recent contribution by [10] introduced a continuation method capable of generating PDE representations for multi-agent systems with homogeneous interaction laws.

The continuation method in [10] assumes homogeneous interactions, meaning all agents follow identical coupling laws. In many realistic systems, however, interactions are heterogeneous: the strength, range, or functional form may depend on the agent or its position, as in traffic with non-identical drivers, oscillator networks with non-uniform couplings, or multi-agent systems with varying capabilities. Such heterogeneity breaks the symmetries exploited in homogeneous models, preventing direct application of continuation or mean-field techniques and invalidating simple translation-invariant PDE representations.

The first main objective of this paper is to address this limitation by developing a continuation procedure that remains valid in the presence of heterogeneous interactions among agents. A key contribution of our work is the derivation of a first-order continuation method suited for heterogeneous microscopic multi-agent systems evolving on a ring, with either fixed or moving agent positions. Such models are widely employed in the study of oscillator networks [11] and traffic dynamics [3], yet existing formulations typically assume homogeneous interactions and therefore fail to capture the variability that may arise between different agents.

The continuation procedure is particularly valuable because PDE formulations often provide a more compact and analytically tractable description of the underlying dynamics than their discrete counterparts. This continuous representation can be used to establish stability criteria and to design control laws directly at the PDE level. These controls can then be discretized and applied to the original microscopic model, as illustrated in [10] through classical PDE-based control techniques. A similar approach has also been employed for homogeneous multi-agent systems in [6], based on a mean-field approach, where the collective effect of interactions is approximated by an averaged density. Recent work by [7] addresses density control of heterogeneous multi-agent systems within a macroscopic mean-field framework. In that contribution, the precise form of the macroscopic dynamics is not assumed to be fully known; instead, heterogeneity and unmodeled effects are treated as bounded drift terms in a density evolution equation obtained through a mean-field description. The overall continuation-based control paradigm

Sébastien Fueyo and Carlos Canudas-de-Wit are with Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, GIPSA-Lab, F-38000 Grenoble, France (emails sebastien.fueyo@gipsa-lab.fr and carlos.canudas-de-wit@gipsa-lab.fr). This work was supported by grant ANR-22-CE48-0011 COCOON.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Continuation control paradigm</figcaption></figure>

is shown in Fig. 1.

The second contribution of this paper is the design of a stabilization feedback law at the macroscopic level for heterogeneous interactions within the framework of the first-order Optimal Velocity (OV) driver model, i.e. with different driver behavior. The control design relies on a Lyapunov functional approach formulated at the PDE level. Through the proposed continuation procedure, the resulting macroscopic feedback is consistently transferred to the underlying microscopic closed loop dynamics, yielding a distributed control strategy that steers the system toward a prescribed stationary density profile on a ring road. While the present work focuses specifically on the OV driver model and on time-invariant density targets, the continuation-based methodology may be extended to more general models and to time-varying model references, at the expense of additional PDE computations.

## II. PRELIMINARIES AND MODELLING FRAMEWORK

### A. Fixed agents

We consider a network of $N$ agents positioned on a ring of length $2\pi$. Each agent $i$ has a fixed spatial location $x_i \in [0, 2\pi)$ and an associated state $z_i(t) \in \mathbb{R}$ at time $t$. The positions $x_i$ are assumed to be uniformly spaced, with $\Delta x := x_{i+1} - x_i$ for all $i = 1, \dots, N$. The evolution of each agent's state $z_i(t)$ is described in Eulerian coordinates by

$$z_i(t) = F_i(z_1(t), \dots, z_N(t)), \quad i = 1, \dots, N, \quad (1)$$

where $z_i(t)$ represents the state of agent $i$ (e.g., phase, velocity, or another relevant quantity), and $F_i : \mathbb{R}^N \to \mathbb{R}$ captures the interactions of agent $i$ with all other agents in the network. This formulation allows for heterogeneous and possibly nonlinear interactions. As a representative example, we consider the Kuramoto oscillator model.

**Example 1 (Kuramoto Oscillators).** We consider a network of coupled oscillators, a central model in the study of collective synchronization. In this context, the state variable $z_i$ used in the general formulation (1) corresponds specifically to the phase of an oscillator. For notational clarity and following standard conventions, we denote the phase of oscillator $i$ by $\theta_i$, whose dynamics is given by

$$\dot{\theta}_i = \omega_i + \sum_{j=1}^{N} H_{ij}(\theta_j - \theta_i), \quad i = 1, \dots, N, \quad (2)$$

where $\omega_i$ is the natural frequency of oscillator $i$, and $H_{ij} : \mathbb{R} \to \mathbb{R}$ represents the heterogeneous coupling function between oscillators $i$ and $j$.

Due to the complexity of dealing with a large number of agents, it is common to adopt macroscopic models in the form of partial differential equations (PDEs). To the best of our knowledge, no general systematic method exists when dealing with heterogeneous microscopic models to derive such macroscopic PDEs, although several approaches have been proposed for particular instances of System (1). As an illustrative example, we consider a continuous formulation of the Kuramoto oscillator on a ring. Following [11], let us consider the discrete Kuramoto model with nearest-neighbor coupling for $i = 1, \dots, N$:

$$\dot{\theta}_i = \omega_i + F(\sin(\theta_{i+1} - \theta_i) - \sin(\theta_i - \theta_{i-1})), \quad (3)$$

with positions $x_i = i\Delta x$ and $\Delta x = 2\pi/N$. The continuum limit of (3) can be understood as

$$\partial_t \theta(t, x) = \omega(x) + \partial_x (G(\partial_x \theta(t, x))), \quad (4)$$

where $G(u) = F(\Delta x \sin(u))$.

We now turn to System (1) expressed in Lagrangian coordinates, meaning that the agents are no longer fixed in space.

### B. Moving agents

In the previous section, we considered agents with fixed positions $x_i$ on a ring, with their dynamics described by the state variables $z_i(t)$. In many applications, however, the positions of the agents themselves evolve in time. In such cases, the state of each agent corresponds directly to its position, and we denote it by $x_i(t)$ to emphasize its temporal evolution. Let $N \in \mathbb{N}$ and consider the collection of variables $\{x_i(t)\}_{i=1}^N$, representing the positions of $N$ vehicles moving along a circular road of length $2\pi$. Each position $x_i(t)$ is defined modulo $2\pi$, reflecting the periodic geometry of the ring. The evolution of each vehicle is governed by the general heterogeneous interaction model

$$\dot{x}_i(t) = F_i(x_1(t), \dots, x_N(t)), \quad i = 1, \dots, N. \quad (5)$$

To illustrate the variety of models encompassed by (5), we consider a large-scale multi-agent interaction model, from which the Optimal Velocity traffic model is derived.

**Example 2 (Multi-agent and OV traffic model).** We consider systems of interacting agents where each agent's velocity depends on relative positions. In a general setting, agent $i$ evolves as

$$\dot{x}_i = \sum_{j=1}^{N} F_{i,j}(x_i - x_j), \quad i = 1, \dots, N, \quad (6)$$

where $F_{i,j}$ models pairwise interactions. A classical special case is the Optimal Velocity (OV) model in traffic flow, where each vehicle reacts only to its leader:

$$\dot{x}_i = V_i(x_{i+1} - x_i), \quad i = 1, \dots, N, \quad (7)$$

with the convention that $x_{N+1} = x_1$ and

$$V_i(s) = v_{\max,i} \left( 1 - \frac{s_{\text{stop}}}{s} \right), \quad (8)$$

where $v_{\max,i} > 0$ is the maximal speed and $s_{\text{stop}} > 0$ the standstill spacing.
As in the case of fixed (Eulerian) positions, there is no general or systematic method to derive a macroscopic PDE from the heterogeneous system (5) for Lagrangian models. However, for the specific structures appearing in Example 2 with homogeneous structure, two different techniques can be employed to derive continuum approximations.

1) Consider the large-scale multi-agent system (6). Under the homogeneity assumption $F_i \equiv F$, a mean-field approach [6] yields, in the continuum limit, the mass conservation law

$$ \partial_t \rho + \partial_x (\rho \mathcal{V}) = 0, \quad (9) $$

with $\rho(x, t)$ denoting the agent density and $\mathcal{V}(x, t)$ the mean velocity field perceived at position $x$ and time $t$. This field is defined by

$$ \mathcal{V}(x, t) = \int_{0}^{2\pi} f(x - y) \rho(y, t) dy, \quad (10) $$

where $f$ is the interaction kernel that quantifies the influence of agents at position $y$ on the velocity at $x$, under periodic boundary conditions.

2) Under the homogeneity assumption $V_i \equiv V$ for the car- following model (7), the continuum limit leads to the classical Lighthill–Whitham–Richards (LWR) model [9]

$$ \partial_t \rho + \partial_x \phi = 0, \quad \phi := \rho V \left(\frac{1}{\rho}\right), \quad (11) $$

with periodic boundary conditions. These derivations, how- ever, rely on homogeneity and do not directly extend to heterogeneous interactions.

In general, when the functions $F_i$ are heterogeneous or depend on more than just the spacing $x_{i+1} - x_i$, deriving a macroscopic PDE is nontrivial. In the following, we present a systematic continuation method to handle such first-order agent systems and construct their PDE limits, starting with the Eulerian formulation, which is simpler due to the agents’ fixed positions and corresponds to the first main objective of this paper.

### III. CONTINUATION FOR FIXED AGENTS

To interpret the microscopic model (1) in a continuum framework, we introduce a function $z : \mathbb{R}_+ \times [0, 2\pi] \to \mathbb{R}$ satisfying

$$ z(t, x_i) = z_i(t), \quad i = 1, \dots, N. $$

Because the agents lie on a ring, two paths connect any pair $(i, j)$; the effective distance between $x_i$ and $x_j$ is defined as the length of the shorter arc,

$$ d_N(i, j) := \min (|j - i|, N - |j - i|). $$

Using this circular distance, a first-order Taylor expansion of $z(t, x)$ around $x_i$ yields

$$ z(t, x_j) \approx z(t, x_i) + d_N(i, j) \Delta x \partial_x z(t, x_i), \quad j = 1, \dots, N. \quad (12) $$

Substituting this approximation into (1), we define the contin- uous mapping

$$ \dot{z}(t, x_i) = \mathcal{F}_i(z(t, x_i), \partial_x z(t, x_i)), \quad (13) $$

where

$$ \mathcal{F}_i(z, \partial_x z) := F_i \left( z + d_N(i, 1) \Delta x \partial_x z, \dots, \right. \\ \left. z + d_N(i, N) \Delta x \partial_x z \right). \quad (14) $$

To embed the discrete system (1) into a continuum descrip- tion, we introduce a set of smooth interpolation functions $\{\lambda_i(x)\}_{i=1}^N$ satisfying $\lambda_i(x_j) = \delta_{ij}$, for $i, j = 1, \dots, N$, where $\delta_{ij}$ is the Kronecker delta, i.e., $\delta_{ij} = 1$ if $i = j$ and $\delta_{ij} = 0$ otherwise. Each function $\lambda_i(x)$ is smooth and vanishes at all locations $x_j$ with $j \neq i$, so that it “activates” only the contribution of the $i$-th agent. Using these interpolation functions, we define a smooth, continuous function

$$ \mathcal{F}(x, z, \partial_x z) := \sum_{i=1}^{N} \lambda_i(x) \mathcal{F}_i(z, \partial_x z), \quad (15) $$

with $z$ and $\partial_x z$ evaluated at $(t, x)$. The function $\mathcal{F}$ interpo- lates the discrete agent dynamics, satisfying $\mathcal{F}(x_i, z, \partial_x z) = \mathcal{F}_i(z, \partial_x z)$ at each agent position $x_i$. This construction embeds the discrete dynamics into a continuous framework suitable for a PDE description. This allows us to interpret (1) as the one- dimensional PDE

$$ \partial_t z(t, x) = \mathcal{F}(x, z(t, x), \partial_x z(t, x)), \quad (16) $$

with the partial derivatives $\partial_t z := \frac{\partial z}{\partial t}$ and $\partial_x z := \frac{\partial z}{\partial x}$, and with periodic boundary conditions $z(t, 0) = z(t, 2\pi)$ due to the ring geometry. We summarize in the next lemma the first order continuation of a microscopic heterogeneous model with fixed agent.

**Lemma III.1.** Within a continuum framework, Equation (1) with fixed positions (i.e., in Eulerian coordinates, with agents uniformly spaced along the ring) can be regarded as an heterogeneous first-order PDE

$$ \left\{ \begin{array}{l} \partial_t z = \mathcal{F}(x, z, \partial_x z), \\ z(t, 0) = z(t, 2\pi). \end{array} \right. \quad (17) $$

Let us show on a simple heterogeneous Kuramoto oscillator the result of the continuation method.

**Example 3.** Consider, for $i = 1, \dots, N$, the system:

$$ \dot{\theta}_i = \omega_i + a_i (\sin(\theta_{i+1} - \theta_i) - \sin(\theta_i - \theta_{i-1})), \quad (18) $$

where $a_i \in \mathbb{R}$. Let us define smooth functions such that $a(x_i) = a_i$ and $\omega(x_i) = \omega_i$. Put:

$$ \mathcal{F}(x, \theta, \partial_x \theta) := \omega(x) + 2a(x) \sin \partial_x \theta. \quad (19) $$

Then the continuum limit (18) is given by Equation (17), with $\mathcal{F}$ given by (19).

We now turn to the continuation for moving agents.

### IV. CONTINUATION FOR MOVING AGENTS

We aim at interpreting the microscopic model (5) in a con- tinuum framework. To derive the macroscopic model, we build upon the continuation method introduced in [10]. While [10] focuses on homogeneous systems, we extend this approach here to handle heterogeneous microscopic models, providing a
systematic procedure for such cases. The discrete vehicle index $i \in \{1, \dots, N\}$ is replaced by a continuously varying index $M \in [0, N]$, interpreted as a virtual continuum of vehicles or Moskowitz function. This construction extends the discrete positions $\{x_i\}_{i=1}^N$ to a continuous function $x(t, M)$ satisfying

$$x(t, M) = x_i(t) \quad \text{for} \quad M = i \in \{1, \dots, N\},$$

from which the macroscopic quantities such as density and velocity can be defined.

### A. Conservation law

Following [8], we introduce the Moskowitz function $M(t, x)$, defined for $(t, x) \in \mathbb{R}_+ \times [0, 2\pi]$, as the cumulative number of vehicles located upstream of position $x$ at time $t$, measured from the reference vehicle $x(t, 0)$. It satisfies

$$M(t, x_i(t)) = i, \quad M(t, 0) = M(t, 2\pi),$$

so that $M$ coincides with the continuous vehicle index introduced earlier. The density $\rho$ and the flow rate $\phi$ are obtained by differentiation:

$$\rho := \partial_x M, \quad \phi := -\partial_t M. \qquad (20)$$

The first relation expresses how the cumulative count varies along the road, while the second measures its rate of change in time. Combining both and using the commutation of mixed derivatives yields the standard conservation law

$$\partial_t \rho + \partial_x \phi = 0. \qquad (21)$$

Defining the velocity field as

$$v(t, x) := \frac{\phi(t, x)}{\rho(t, x)}, \qquad (22)$$

we recover the classical transport form

$$\partial_t \rho + \partial_x (v \rho) = 0, \qquad (23)$$

which expresses the conservation of vehicles along the ring road. It remains to exploit the equation of (5) to show that $v$ is a nonlinear function depending on $\rho$.

### B. Speed law

First notice that the functions $x(t, \cdot)$ and $M(t, \cdot)$ are inverse to each other where, for $z \in [0, 2\pi]$ and $z_0 \in [0, N]$, it holds that $x(t, M(t, z)) = z$ and $M(t, x(t, z_0)) = z_0$. Differentiating these relations gives

$$\partial_M x \cdot \partial_x M = 1, \quad \partial_t M + \partial_t x \cdot \partial_x M = 0. \qquad (24)$$

Hence, from Equations (20)-(22)-(24), it follows that

$$\partial_M x = \frac{1}{\rho}, \quad v = \partial_t x. \qquad (25)$$

These identities highlight the correspondence between the microscopic trajectories $x_i(t)$ and the macroscopic fields $(\rho, v)$, linking vehicle density, flow, and speed through the Moskowitz function. First, we note that the velocity field at the discrete positions satisfies $v(t, x_i(t)) \approx \dot{x}_i(t)$, leading to the approximation

$$v \approx \dot{x}_i. \qquad (26)$$

It remains to approximate the $x_i(t)$. By definition:

$$x_j(t) - x_i(t) = (j - i) \frac{x(t, M(t, x_{i+1}(t))) - x(t, M(t, x_i(t)))}{M(t, x_j(t)) - M(t, x_i(t))}. \qquad (27)$$

Thus, applying a first-order Taylor expansion, the finite difference $x_j - x_i$ can then be approximated by the partial derivative of $x$ with respect to $M$, yielding

$$x_j(t) - x_i(t) \approx (j - i) \partial_M x(t, i). \qquad (28)$$

Consequently, from (25), we can approximate the positions themselves as

$$x_j(t) \approx x_i(t) + (j - i) \frac{1}{\rho(t, x_i(t))}, \qquad (29)$$

where $\rho$ denotes the local agent density. Using the approximation given in (29), together with Equations (5)–(26), we obtain that the discrete velocity function can be expressed as

$$v(t, x_i) \approx F_i \left( x_i + (j-1) \frac{1}{\rho}, \dots, x_i + (j-N) \frac{1}{\rho} \right). \qquad (30)$$

Finally, by performing a spatial continuation of the discrete mappings $F_i$, we construct a smooth continuum velocity field at the agent positions:

$$v(t, x) = F_i(t, x, \rho). \qquad (31)$$

Here, interpolation refers to the process of extending the discrete agent dynamics $F_i$ to a smooth function of both space and time. Since the agents move, their positions $x_i(t)$ change with time, and each agent's contribution must be tracked dynamically. We achieve this by defining smooth spatio-temporal interpolation functions $\{\lambda_i(t, x)\}_{i=1}^N$ satisfying $\lambda_i(t, x_j(t)) = \delta_{ij}$, for $i, j = 1, \dots, N$, where $\delta_{ij}$ is the Kronecker delta. Each $\lambda_i(t, x)$ vanishes at all positions $x_j(t)$ with $j \neq i$, ensuring that the $i$-th agent contributes only around its current location, while the smoothness spreads this contribution continuously over space. Using these functions, the continuum velocity field is given by

$$\mathcal{F}(t, x, \rho) := \sum_{i=1}^{N} \lambda_i(t, x) \mathcal{F}_i(x, \rho). \qquad (32)$$

The smooth spatio-temporal interpolation functions $\{\lambda_i(t, x)\}_{i=1}^N$ exist under the mild assumption that the agent positions $x_i(t)$ remain distinct for all $t$ and vary smoothly in time. This ensures that each agent's contribution can be smoothly extended to a continuous function of space and time, embedding the discrete moving-agent dynamics into a smooth spatio-temporal framework suitable for PDE analysis. This construction provides a smooth, time-dependent interpolation of the discrete agent velocities, yielding a well-defined continuum representation. Thus, the discrete system (5), with the use of (31) and (23) can be interpreted as a 1-D PDE in the continuum limit:

$$\partial_t \rho + \partial_x (\rho \mathcal{F}(t, x, \rho)) = 0. \qquad (33)$$

Finally, due to the ring geometry of the system, the agent density satisfies the periodic boundary condition $\rho(t, 0) = \rho(t, 2\pi)$
for $t \ge 0$. The following lemma presents the first-order continuum approximation of a heterogeneous microscopic model with moving agents.

**Lemma IV.1.** *Within a continuum framework, Equation (5) with moving agents (i.e., in Lagrangian coordinates) can be regarded as an heterogeneous first-order PDE*

$$ \left\{ \begin{array}{l} \partial_t \rho + \partial_x (\rho \mathcal{F}(t, x, \rho)) = 0, \\ \rho(t, 0) = \rho(t, 2\pi). \end{array} \right. \qquad (34) $$

First notice that Lemma IV.1 allow us to recover the macroscopic system (11) from the microscopic (7) with $V_i \equiv V$, thus our work generalizes in the heterogeneous case the work in [9].

**Example 4.** Consider the car-following model of Example 2, where each driver adjusts their velocity according to Equations (7)–(8). Let $\{x_i(t)\}_{i=1}^N$ denote the vehicle positions and $\{v_{\max,i}\}_{i=1}^N$ their individual maximal speeds. We define a smooth, spatially continuous maximal velocity function $v_{\max}^c(t, x)$ by Gaussian smoothing of the discrete velocities:

$$ v_{\max}^c(t, x) := \sum_{i=1}^{N} v_{\max,i} \frac{\exp\left(-\frac{(x-x_i(t))^2}{2\sigma^2}\right)}{\sum_{j=1}^{N} \exp\left(-\frac{(x-x_j(t))^2}{2\sigma^2}\right)}, \quad (35) $$

where $\sigma > 0$ is an enough small smoothing parameter. The corresponding continuum flux is then

$$ \mathcal{F}(t, x, \rho) := v_{\max}^c(t, x) (1 - s_{\text{stop}}) \rho, \quad (36) $$

and the continuum limit of the car-following dynamics is given by the LWR-type PDE (18) with $\mathcal{F}(t, x, \rho)$ defined in (36).

The second objective is to use the macroscopic model to study control properties of the microscopic system.

## V. CONTINUATION-BASED FEEDBACK CONTROL AT THE MACROSCOPIC LEVEL FOR THE OPTIMAL VELOCITY MODEL

The continuation method has proven effective for the control of large-scale networks in both homogeneous [4], [6], [10] and heterogeneous [7] settings. Here, we illustrate its potential within the first-order Optimal Velocity (OV) traffic model, enabling the design of a macroscopic feedback law that stabilizes a prescribed stationary density on a ring road. While we focus on the OV model with time-invariant target densities, the method might be extended to time-dependent target densities. Such an extension, however, would require additional PDE computations.

We consider the OV model (7)–(8), which describes each vehicle’s natural headway-based driving behavior, to which we add a control input $u_i(t)$ representing additional driver or driver-assist actions (e.g., overtaking or speed regulation), while preserving the intrinsic dynamics. The resulting controlled dynamics is

$$ \dot{x}_i(t) = v_{\max,i} \left( 1 - \frac{s_{\text{stop}}}{x_{i+1}(t) - x_i(t)} \right) + u_i(t), \quad (37) $$

where $v_{\max,i}$ and $s_{\text{stop}}$, $i$ reflect driver heterogeneity, ensuring each vehicle retains its individual behavior under the applied control.

Applying the continuation procedure to the microscopic control inputs $u_i(t)$, we introduce a macroscopic control function $u(t, x)$ such that $u(t, x_i) = u_i(t)$, for $i = 1, \dots, N$. Following the same steps as in Section IV, the microscopic control problem (37) can then be reformulated in the continuum limit, where the macroscopic density $\rho(t, x)$ evolves according to

$$ \partial_t \rho(t, x) + \partial_x(\rho(t, x)(\mathcal{F}(t, x, \rho(t, x)) + u(t, x))) = 0, \quad (38) $$

with periodic boundary conditions and $\mathcal{F}(t, x, \rho)$ defined in (36). We aim at steering the system toward a prescribed target reference density, denoted by $\rho^r(x)$. For illustration purposes, we define $\rho^r(x)$ as a unimodal, time-independent distribution:

$$ \rho^r(x) = \frac{2N}{Z} e^{\kappa \cos(x - \mu_1)}, \quad (39) $$

with $\kappa > 0$ controlling the peak, $\mu_1 = \pi$, and normalization $Z = 4\pi I_0(\kappa)$. Defining the deviation error $\rho^e := \rho - \rho^r$, the error dynamics become

$$ \left\{ \begin{array}{l} \partial_t \rho^e + \partial_x(\rho(\mathcal{F}(t, x, \rho) + u)) = 0, \\ \rho^e(t, 0) = \rho^e(t, 2\pi). \end{array} \right. \qquad (40) $$

Motivated by the structure of the deviation equation (40), we introduce the following feedback control law, parameterized by $K > 0$:

$$ u(t, x) = \frac{-\rho^r(x) \mathcal{F}(t, x, \rho(t, x)) + K \int_0^x (\rho - \rho^r) ds}{\rho(x)}. \quad (41) $$

Under this feedback, the error equation (40) reduces to

$$ \partial_t \rho^e - \partial_x (\rho^e \mathcal{F}(t, x, \rho(t, x))) + K \rho^e = 0, \quad (42) $$

which is a scalar one-dimensional linear hyperbolic PDE. Let $\|f\|_\infty := \sup_{x \in [0, 2\pi]} |f(x)|$ and $\|f\|_{L^2} := (\int_0^{2\pi} |f(x)|^2 dx)^{1/2}$ denote the uniform and $L^2$ norms on $[0, 2\pi]$, respectively. From now on, we assume the existence of sufficiently regular solutions that remain uniformly bounded from below and above; we refer to [5] for a discussion on regularity issues in hyperbolic conservation laws. We show that the error density $\rho^e$ converges exponentially to zero in the $L^2$ sense.

**Theorem V.1.** Under the feedback law $u$ defined in (41), the density tracking error $\rho^e := \rho - \rho^r$ converges exponentially in $L^2$-norm.

*Proof.* Let $V(t)$ be the Lyapunov functional defined by $V(t) := \frac{1}{2} \|\rho^e(t, \cdot)\|_{L^2}^2$, for $t \ge 0$. Differentiating $V$ yields

$$ \dot{V}(t) = \int_{0}^{2\pi} \rho^e(t, x) \partial_t \rho^e(t, x) dx. \quad (43) $$

Substituting the error dynamics (40) into (43), we obtain

$$ \dot{V}(t) = - \int_{0}^{2\pi} \rho^e \partial_x(\rho^e \mathcal{F}(x, \rho)) dx - K \|\rho^e\|_{L^2}^2. \quad (44) $$

Using integration by parts and the periodic boundary condition, the transport term can be rewritten as

$$ \dot{V}(t) = -\frac{1}{2} \int_{0}^{2\pi} (\rho^e)^2 \partial_x(\mathcal{F}(t, x, \rho)) dx - K \|\rho^e\|_{L^2}^2. \quad (45) $$
Taking the spatial derivative of $\mathcal{F}$ defined in (36) and using the decomposition $\rho = \rho^e + \rho^r$ and $\partial_x \rho = \partial_x \rho^e + \partial_x \rho^r$, it yields that

$$ \dot{V}(t) = -\frac{1}{2}(I_1 + I_2 + I_3 + I_4) - K \|\rho^e\|_{L^2}^2, $$

where

$$ \begin{aligned} I_1 &= \int_0^{2\pi} (\rho^e)^2 \partial_x v_{\max}^c(t, x) dx, \\ I_2 &= -s_{\text{stop}} \int_0^{2\pi} (\rho^e)^2 \partial_x v_{\max}^c(t, x) (\rho^r + \rho^e) dx, \\ I_3 &= -s_{\text{stop}} \int_0^{2\pi} (\rho^e)^2 v_{\max}^c(t, x) \partial_x \rho^r dx, \\ I_4 &= -s_{\text{stop}} \int_0^{2\pi} (\rho^e)^2 v_{\max}^c(t, x) \partial_x \rho^e dx. \end{aligned} \quad (46) $$

By construction, we have:

$$ \begin{aligned} |I_1| &\le \|\partial_x v_{\max}^c\|_{L^\infty} \|\rho^e\|_{L^2}^2, \\ |I_2| &\le s_{\text{stop}} \|\partial_x v_{\max}^c\|_{L^\infty} (\|\rho^r\|_{L^\infty} + \|\rho^e\|_{L^\infty}) \|\rho^e\|_{L^2}^2, \\ |I_3| &\le s_{\text{stop}} \|v_{\max}^c\|_{L^\infty} \|\partial_x \rho^r\|_{L^\infty} \|\rho^e\|_{L^2}^2. \end{aligned} $$

Hence, these three terms are quadratic in $\|\rho^e\|_{L^2}$. Finally, the remaining term reads

$$ I_4 = -s_{\text{stop}} \int_0^{2\pi} (\rho^e)^2 v_{\max}^c(t, x) \partial_x \rho^e dx. $$

Integrating by parts and using periodicity yields

$$ \int_0^{2\pi} (\rho^e)^2 v_{\max}^c \partial_x \rho^e dx = -\frac{1}{3} \int_0^{2\pi} \partial_x v_{\max}^c(t, x) (\rho^e)^3 dx, $$

so that

$$ |I_4| \le \frac{s_{\text{stop}}}{3} \|\partial_x v_{\max}^c\|_{L^\infty} \|\rho^e\|_{L^\infty} \|\rho^e\|_{L^2}. \quad (47) $$

To sum up, collecting the above estimates and using the assumption that the solutions are uniformly bounded from above, there exist constants $C_1 > 0$ such that

$$ \dot{V}(t) \le -(K - C_1) \|\rho^e\|_{L^2}^2, \quad (48) $$

so that for $K$ enough large, we have the exponential stability. □

We now turn to the numerical illustration of Theorem V.1.

a) **Microscopic feedback law:** At the microscopic level, $u_i(t)$ is obtained by interpolating the macroscopic control $u(t, x)$ at the agent positions, i.e., $u_i(t) = u(t, x_i(t))$. Implementation aspects are discussed below.

b) **Micro-macro coupling and implementation:** The method relies on a dynamic coupling between the microscopic OV model and its macroscopic continuum approximation. At each time step, the macroscopic density $\rho(t, \cdot)$ evolves according to (38), the feedback control $u(t, x)$ is computed on the grid and evaluated at vehicle positions to define $u_i(t)$, and the microscopic dynamics (37) are updated. The new agent positions are then used to reconstruct the macroscopic velocity field $\mathcal{F}(t, x, \rho)$ for the next time step. This closed-loop micro-macro procedure ensures that the macroscopic control steers the density toward $\rho^r$, while the microscopic dynamics provide the instantaneous configuration required for a consistent continuation.

c) **Illustration:** Fig. 2 shows the empirical density $\rho(x)$ and vehicle positions at the initial time and at the final time $T = 9$ under the feedback law. The macroscopic control drives the density toward the target profile $\rho^r$, and the vehicles rearrange on the ring to match this target.

<figure><img src="image_4.png" /><figcaption>Fig. 2. Microscopic feedback control of the OV model with $N = 40$ vehicles on a ring of length $L = 2\pi$ and feedback gain $K = 15$. Top: initial vehicle density (red line) and vehicle positions on the ring (blue dots), with target density (black dashed). Bottom: empirical density and vehicle positions after feedback converge toward the target density. The maximal speeds $v_{\max,i}$ of the vehicles are chosen randomly in $[1, 5]$, and the minimal headway is $s_{\text{stop}} = 0.1$.</figcaption></figure>

REFERENCES

[1] A. Aw, A. Klar, M. Rascle, and T. Materne. Derivation of continuum traffic flow models from microscopic follow-the-leader models. *SIAM Journal on applied mathematics*, 63(1):259–278, 2002.

[2] N. Baloul, A. Hayat, T. Liard, and P. Lissy. Traffic flow reconstruction from limited collected data. In *2025 IEEE 64th Conference on Decision and Control (CDC)*, pages 7981–7986, 2025.

[3] S. Cui, B. Seibold, R. Stern, and D. B. Work. Stabilizing traffic flow via a single autonomous vehicle: Possibilities and limitations. In *2017 IEEE Intelligent Vehicles Symposium (IV)*, pages 1336–1341, 2017.

[4] S. Fueyo and C. Canudas de Wit. A continuation-based control design for stabilizing second-order macroscopic traffic flow on circular roads. In *2025 IEEE 64th Conference on Decision and Control (CDC)*, pages 230–235, 2025.

[5] M. Garavello, K. Han, B. Piccoli, et al. *Models for vehicular traffic on networks*, volume 9. American Institute of Mathematical Sciences (AIMS), Springfield, MO, 2016.

[6] G. C. Maffettone, A. Boldini, M. Di Bernardo, and M. Porfiri. Continuation control of large-scale multiagent systems in a ring. *IEEE Control Systems Letters*, 7:841–846, 2023.

[7] G. C. Maffettone, D. Salzano, and M. di Bernardo. Robust macroscopic density control of heterogeneous multi-agent systems, 2026.

[8] G. Newell. A simplified theory of kinematic waves in highway traffic, part I: General theory. *Transportation Research Part B: Methodological*, 27(4):281–287, 1993.

[9] D. Nikitin. *Scalable Large-Scale Control of Network Aggregates*. PhD dissertation, Université Grenoble Alpes, Sep 2021.

[10] D. Nikitin, C. Canudas-de Wit, and P. Frasca. A continuation method for large-scale modeling and control: From ODEs to PDE, a round trip. *IEEE Transactions on Automatic Control*, 67(10):5118–5133, 2022.

[11] D. Nikitin, C. Canudas-de Wit, P. Frasca, and U. Ebels. Synchronization of spin-torque oscillators via continuation method. *IEEE Transactions on Automatic Control*, 68(11):6621–6635, 2023.