# Mean-Field Control of Heterogeneous Piecewise Deterministic Markov Processes under Grid Constraints

Thomas Le Corre, Ana Bušić

► To cite this version:

Thomas Le Corre, Ana Bušić. Mean-Field Control of Heterogeneous Piecewise Deterministic Markov Processes under Grid Constraints. SmartGridComm 2025 - IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids, Sep 2025, Toronto, Canada. hal-05375490

HAL Id: hal-05375490

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

Submitted on 21 Nov 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
# Mean-Field Control of Heterogeneous Piecewise Deterministic Markov Processes under Grid Constraints

Thomas Le Corre

Inria, Paris, France

DI ENS, ENS, PSL University

Paris, France

thomas.le-corre@inria.fr

Ana Bušić

Inria, Paris, France

DI ENS, ENS, PSL University

Paris, France

ana.busic@inria.fr

***Abstract***—This paper investigates mean-field control problems for large populations of agents modeled by piecewise deterministic Markov processes. Motivated by the growing need for demand-side management in power systems, we develop a decentralized control framework that accounts for key operational grid constraints such as maximum aggregate power and ramping limits. These constraints are critical to ensure the feasibility and reliability of control strategies in real-world power grids. Furthermore, we address the limitation of classical mean-field approaches that assume agent homogeneity by extending the formulation to heterogeneous populations, where agent-specific characteristics (such as charging power or battery capacity) are encoded through fixed parameters. The resulting framework enables the scalable and feasible coordination of various flexible loads. Theoretical developments are illustrated through applications to electric vehicle charging, demonstrating the impact of the proposed constraints and heterogeneity modeling on the system's aggregate behavior.

***Index Terms***—mean field control, electric vehicles, demand side management, piecewise-deterministic Markov processes

## I. INTRODUCTION

Mean-field control problems have attracted considerable attention due to their wide range of applications in fields such as economics, energy systems, and networked infrastructures [1], [2]. When a large population of agents interacts through a global cost function, solving the corresponding control problem directly becomes computationally intractable due to the high dimensionality. The mean-field approximation addresses this issue by modeling the behavior of a representative agent subjected to the influence of the aggregate population state. This formulation naturally leads to decentralized control strategies aimed at aligning the behavior of individual agents with a shared global objective [3]. Piecewise deterministic Markov processes (PDMPs) offer a robust framework for modeling agents characterized by sudden transitions and deterministic dynamics between these transitions [4]. Such models have already been applied to Electric Vehicles (EVs), with a decentralized approach based on a Lagrangian formulation [5]. This methodology leverages a stochastic Uzawa algorithm, enabling efficient handling of the non-convex nature of the problem.

This work was carried out in the framework of the AI-NRGY project, funded by France 2030 (Grant No: ANR-22-PETA-0004).

This research is particularly relevant to demand-side management in power systems, which aims to maintain a balance between energy production and consumption. With the increasing integration of intermittent renewable energy sources, such as wind and solar power, managing electricity demand has become crucial to mitigate the variability in power generation and various strategies have been explored [6]–[8].

To ensure operational reliability, demand-side control strategies must account for physical grid constraints. Among the most critical are capacity limitations, such as bounds on the total consumption and ramping constraints, which limit the rate at which aggregate demand can change over time. These constraints are essential to prevent network congestion, reduce the risk of voltage instability, and guarantee secure operation, particularly in the presence of large populations of controllable loads. As emphasized in [9]–[11], failing to incorporate such network-aware constraints may lead to infeasible or suboptimal control strategies. In this work, we explicitly address these issues by integrating both maximum power and demand ramping constraints into the mean-field control framework.

A key challenge in demand-side control lies in the inherent heterogeneity of the controlled population. In the context of electric vehicle (EV) charging, individual units may differ widely in terms of battery capacity, initial state-of-charge, charging power capabilities, and user preferences. Classical mean-field control frameworks often assume a homogeneous population, where all agents are identically distributed and governed by the same dynamics. This simplifying assumption, while analytically convenient, limits the applicability of the model to realistic systems. Addressing this limitation requires extending the mean-field paradigm to incorporate agent heterogeneity in a tractable manner. This problem was addressed early on by approaches that added parameters to the state space to take account of this heterogeneity [12]. Recent works have begun to explore such extensions, either through structured population models [1], heterogeneous mean-field games [13], or parametric distributions over agent types [14]. In this paper, we build on this line of research by introducing a heterogeneous mean-field control formulation, where agent-specific characteristics are encoded through fixed parameters, and integrate it within a PDMP-based model.
**Contributions:** This paper aims to make PDMP-based mean-field control approaches more realistic, by extending the theoretical results of this approach. There are two major contributions:

* Typical grid constraints are taken into account: ramping and maximum power constraints.

* The PDMP-based mean-field approach is extended to heterogeneous agents.

**Structure:** Section II introduces the main mathematical tools used in the article, such as the PDMP definition. The mean field control problem is introduced in Section III. This problem with ramping and maximum power constraints is illustrated through an application to electric vehicle (EV) charging, where these are shown to play a critical role in shaping charging strategies. In Section IV, we further adapt the framework to account for agent heterogeneity. A second example demonstrates this extension in a setting where EVs differ in their charging power.

## II. PRELIMINARIES

### A. Notations

To facilitate readability, we introduce the main notations used throughout the paper.

* Let $\mathcal{X}$ be the state space defined as $\mathcal{X} := \mathcal{I} \times [0, 1]$, where $\mathcal{I}$ is a finite set.

* The time horizon is the interval $[0, T]$ with $T > 0$.

* The space $C^1(\mathcal{X})$ (resp. $C^1([0, T] \times \mathcal{X})$) represents the set of continuously differentiable real-valued functions defined on $\mathcal{X}$ (resp. $[0, T] \times \mathcal{X}$).

* The space of Borel probability measures on $\mathcal{X}$ is denoted by $\mathcal{P}(\mathcal{X})$.

* The space $C_0([0, T] \times \mathcal{X}, \mathbb{R}_+^d)$ (resp. $C_0([0, T], \mathcal{P}(\mathcal{X}))$) denotes the set of continuous functions from $[0, T] \times \mathcal{X}$ to $\mathbb{R}_+^d$ (resp. $\mathcal{P}(\mathcal{X})$).

* $L^2(0, T)$ is the space of square-integrable functions on $[0, T]$ and $L^2_+(0, T)$ is the space of non-negative square-integrable functions on $[0, T]$.

* We define the following function space:

$$A := \{ \alpha \in C_0([0, T] \times \mathcal{X}, \mathbb{R}_+^d) \mid \forall i \in \mathcal{I}, \alpha_i(\cdot, i, \cdot) = 0 \},$$

### B. Modelization by Piecewise Deterministic Markov Processes

In this article, agents will be modeled by Piecewise Deterministic Markov Process (PDMP) [4] controlled by a common intensity $\alpha \in A$. The state of an agent at time $t$ is given by $X_t^\alpha = (I_t^\alpha, S_t^\alpha)$, where:

* $S_t^\alpha \in \mathbb{R}$ represents a continuous state variable (e.g., level battery of an EV), evolving according to the ordinary differential equation:

$$\frac{d}{dt} S_t^\alpha = b(X_t^\alpha),$$

where $b: \mathcal{X} \to \mathbb{R}$ describes the speed of $S_t$. For any $(\tau, t, j, s) \in [0, T] \times [\tau, T] \times \mathcal{X}$, we define the flow $\phi_{\tau,s}$

as the unique solution of this ordinary differential equation $\partial_t \phi_{\tau,s}(j, t) = b(j, \phi_{\tau,s}(j, t))$ with the initial condition $\phi_{\tau,s}(j, \tau) = s$.

* $I_t^\alpha \in \mathcal{I}$ represents a discrete mode (e.g., charging mode for an EV), switching between states in $I$ at jump times $\{T_k\}_{k \ge 0}$, which are determined by a state-dependent intensity function $\alpha \in A$:

$$P(T_{k+1} - T_k > s | X_{T_k} = (Y, s)) = e^{-\int_0^s \sum_{j \neq Y} \alpha_j(t, Y, s) dt}.$$

At each jump time $T_k$, the process transitions from state $(I_{T_k}^-, S_{T_k}^-)$ to a new state $(I_{T_k}, S_{T_k}^-)$:

$$I_{T_k} = i \quad \text{with probability} \quad \frac{\alpha_i(T_k, I_{T_k}^-, S_{T_k}^-)}{\sum_{j' \neq I_{T_k}^-} \alpha_{i'}(T_k, I_{T_k}^-, S_{T_k}^-)}.$$

The PDMP thus alternates between deterministic evolution phases and stochastic jumps, making it a suitable model for systems with sudden transitions and deterministic dynamics between jumps.

**Continuity equation (CE$\alpha$)** We denote $m^\alpha$ the law of a PDMP with control $\alpha$. For a given function $\alpha \in A$, a distribution $m$ satisfies the continuity equation in the weak sense if $m$ is in $C_0([0, T], \mathcal{P}(\mathcal{X}))$ and for any test function $\varphi \in C_c^\infty([0, T] \times \mathcal{X})$, we have:

$$\begin{aligned} & \int_0^1 \sum_i \varphi(T, i, s) m(T, i, ds) - \varphi(0, i, s) m_i^0(ds) \\ &= \int_0^T \int_0^1 \sum_i (\partial_t \varphi(t, i, s) + b(t, i, s) \partial_s \varphi(t, i, s)) m(t, i, ds) dt \\ &+ \int_0^T \int_0^1 \sum_i \sum_{\substack{j \in \mathcal{I} \\ j \neq i}} (\varphi(t, j, s) - \varphi(t, i, s)) \alpha_j(t, i, s) m(t, i, ds) dt. \end{aligned} \quad (1)$$

## III. MEAN FIELD PROBLEM OF PDMP UNDER GRID CONSTRAINTS

### A. Mean Field Control problem

We consider a population of $N$ independent and identically distributed agents. The goal is to minimize the expected cost function $J_N$ over all admissible controls $\alpha \in A$:

$$J_N(\alpha) := \mathbb{E} \left[ \frac{1}{N} \sum_{n=1}^{N} G(\alpha, X^{n, \alpha}) \right] \quad (2)$$

where $G(\alpha, X^{n, \alpha})$ is the individual cost of an agent, given by:

$$G(\alpha, x) := \int_{0}^{T} \left[ c(t, x_t) + \sum_{j \in \mathcal{I}} L(\alpha_j(t, x_t)) \right] dt + g(x_T). \quad (3)$$

We consider in this article, the **optimal mean field control problem** associated to $J_N$, defined as:

$$\begin{aligned} \min_{\alpha \in A} & J(\alpha) = \mathbb{E}[G(\alpha, X^\alpha)] \\ \text{s.t. } & X_t^\alpha \text{ is a PDMP of intensity } \alpha \\ & \text{and } \forall t, \mathbb{E}[p(X_t^\alpha)] \le P^{\text{max}}(t) \\ & \text{and } \forall t, \frac{\partial}{\partial t} \mathbb{E}[p(X_t^\alpha)] \le R^{\text{max}}(t) \end{aligned} \quad (4)$$
This problem can be intuitively interpreted as the limiting case where the number of agents $N$ tends to infinity. It is important to note here that the control parameter is not directly $I_t$ (the mode of heating for water heaters) but the intensity $\alpha$.

We make the following assumptions:

(A1) The continuous component is bounded: $\forall t, S_t \in [0, 1]$ and the vector field $b$ vanishes on the bounds.

(A2) The feasible set in (4) is non empty.

(A3) The function $p$ is in $C^1(\mathcal{X})$.

(A4) $c \in C^1([0, T] \times \mathcal{X})$, $g \in C^1(\mathcal{X})$ and $b$ is continuously differentiable with respect to the second variable.

(A5) The function $L: \mathbb{R} \to \mathbb{R} \cup \{+\infty\}$ is defined as:

$$L(x) := \begin{cases} l(x), & \text{if } x > 0, \\ 0, & \text{if } x = 0, \\ +\infty, & \text{otherwise,} \end{cases}$$

where $l \in C^1(\mathbb{R}_+, \mathbb{R}_+)$ is an increasing strongly convex function bounded from above by a quadratic function. We denote by $H$ the convex conjugate of $L$, i.e. $H$ is defined for any $x \in \mathbb{R}$ by $H(x) := \sup_{y \in \mathbb{R}} (xy - L(y))$.

**Remark:** Assumption (A2) is generally verified by the fact that control $\alpha = 0$ is feasible. In practice, and notably in the use cases considered in this article, this control will be feasible.

## B. Dualization

Let us introduce the Lagrangian $\mathcal{L}: A \times L_+^2(0, T) \times L_+^2(0, T) \to \mathbb{R}$:

$$\begin{aligned} \mathcal{L}(\alpha, \lambda, \eta) := & \int_0^T \left( \mathbb{E}[p(X_t^\alpha)] - P^{\text{max}}(t) \right) \lambda(t) dt \\ & + \int_0^T \left( \frac{\partial}{\partial t} \mathbb{E}[p(X_t^\alpha)] - R^{\text{max}}(t) \right) \eta(t) dt + \mathbb{E}[G(\alpha, X^\alpha)] \end{aligned} \quad (5)$$

and the associated dual function $W: L_+^2(0, T) \times L_+^2(0, T) \to \mathbb{R}$:

$$W(\lambda, \eta) := \inf_{\alpha \in A} \mathcal{L}(\alpha, \lambda, \eta). \quad (P)$$

**Proposition 1.** There exists a unique $(\bar{\lambda}, \bar{\eta}) \in L_+^2(0, T) \times L_+^2(0, T)$ maximizing $W$. And there is no duality gap, i.e.,

$$\max_{\lambda, \eta \in L_+^2(0, T)} W(\lambda, \eta) = \min_{\alpha \in A} \bar{J}(\alpha).$$

*Proof.* The proof is an adaptation of the Lemma 3.2 and Theorem 3.4 in [5]. □

We introduce, for any $\lambda, \eta \in L_+^2(0, T)$, the objective function $\bar{\mathcal{L}}(\cdot, \cdot, \lambda, \eta)$, defined on $A \times C_0([0, T], \mathcal{P}(\mathcal{X}))$, by

$$\begin{aligned} \bar{\mathcal{L}}(\alpha, m, \lambda, \eta) := & \int_0^1 \sum_i g(i, s) m(T, i, ds) + \int_0^T \int_0^1 \sum_i \left( c(t, i, s) \right. \\ & \left. \qquad + \sum_{j \in \mathcal{I}} L(\alpha_j(t, i, s)) + \eta(t) (p(j, s) - p(i, s)) \alpha_j(t, i, s) \right. \\ & \left. \qquad + \eta(t) b(i, s) p(i, s) + \lambda(t) p(i, s) \right) m(t, i, ds) dt \\ & - \int_0^T R^{\text{max}}(t) \eta(t) dt - \int_0^T P^{\text{max}}(t) \lambda(t) dt \end{aligned}$$

and the optimization problem:

$$\inf_{(\alpha, m) \in A \times C_0([0, T], \mathcal{P}(\mathcal{X}))} \bar{\mathcal{L}}(\alpha, m, \lambda, \eta) \quad (\bar{P})$$

where $m$ is a weak solution of $\text{CE}(\alpha)$ defined in 1.

**Proposition 2.** For any $\lambda, \eta \in L_+^2(0, T)$, Problem (P) is equivalent to Problem $(\bar{P})$.

*Proof.* The equivalence $(X^\alpha$ is a PDMP) $\Leftrightarrow$ $(m^\alpha$ is a weak solution of $(\text{CE}(\alpha))$) was established in [15]. The main difference in the proof lies in the expression of the term due to the maximum slope constraint, in this new setting:

$$\begin{aligned} & \int_0^T \left( \frac{\partial}{\partial t} \mathbb{E}[p(X_t^\alpha)] \right) \eta(t) dt \\ &= \int_0^T \left( \int_0^1 \sum_i \frac{\partial}{\partial t} (pm^\alpha)(t, i, s) \right) \eta(t) dt \\ &= \int_0^T \eta(t) \left( \int_0^1 \sum_i \frac{\partial p}{\partial s}(i, s) b(i, s) \right. \\ & \qquad \left. + \sum_{j \neq i} p(i, s) (\alpha_{j,i}(t, s) m_j^\alpha(t, s) - \alpha_{i,j}(t, s) m_i^\alpha(t, s)) \right) dt \\ &= \int_0^T \int_0^1 \sum_i \sum_{j \neq i} (p(j, s) - p(i, s)) \alpha_{i,j}(t, s) m_i^\alpha(t, s) \eta(t) dt \\ & \qquad + \int_0^T \int_0^1 \sum_i p(i, s) b(i, s) m_i^\alpha(t, s) \eta(t) dt \end{aligned}$$

where we used the continuity equation (1) and the assumption that the vector field vanishes on the boundaries. The equivalence between the two problems immediately follows. □

Now, with this proposition, we can obtain solutions $(\alpha, m^\alpha)$ and thus $\alpha$ to our minimization problem. For any $\lambda, \eta \in L_+^2(0, T)$, we define $\varphi[\lambda, \eta]$, the unique function in $C$ (defined in (1.1)), satisfying for any $(t, i, s) \in [0, T] \times \mathcal{X}$,

$$\begin{aligned} \varphi[\lambda, \eta](t, i, s) = & g(i, \phi_{t,s}(i, T)) + \int_t^T c(\tau, i, \phi_{t,s}(i, \tau)) d\tau \\ & + \int_t^T \eta(\tau) b(\tau, i, \phi_{t,s}(i, \tau)) \frac{\partial p}{\partial s}(\tau, i, \phi_{t,s}(i, \tau)) d\tau \\ & - \int_t^T \sum_{\substack{j \in \mathcal{I} \\ j \neq i}} H(\varphi[\eta](\tau, i, \phi_{t,s}(i, \tau)) - \varphi[\eta](\tau, j, \phi_{t,s}(i, \tau))) \\ & - (p(j, \phi_{t,s}(i, \tau)) - p(i, \phi_{t,s}(i, \tau)))\eta(\tau)) d\tau \\ & + \int_t^T \lambda(\tau) p(\tau, i, \phi_{t,s}(i, \tau)) d\tau. \end{aligned}$$

**Proposition 3.** For any $\lambda, \eta \in L_+^2(0, T)$, the control $\alpha[\lambda, \eta]$ defined for any $(t, i, j, s) \in [0, T] \times I^2 \times [0, 1]$ by

$$\alpha_j[\lambda, \eta](t, i, s) := H' \left( \varphi[\eta](t, i, s) - \varphi[\eta](t, j, s) \right. \\ \left. - (p(j, s) - p(i, s))\eta(t) \right) \quad (6)$$

is a solution of the problem (P).
**Proof.** The proof follows the path of the proof of Proposition 4.1 in [5], by replacing their term due to congestion to our terms due to ramping constraint and maximum power. Then, by [16], Theorem 2.1, $(\alpha[\eta], m^{\alpha[\eta]})$ is a solution of problem $(\bar{P})$ and thus to problem $(P)$ (by Proposition 2). □

**Proposition 4.** W is Gâteaux differentiable. For any $\lambda, \eta \in L_+^2(0, T)$, the Gâteaux differential of W at $\eta$ is denoted by $DW(\lambda, \eta)(\cdot)$ and is given in the direction $\mu \in L_+^2(0, T)$ by:

$$DW(\lambda, \eta)(\mu) = \int_{0}^{T} \left( \frac{\partial}{\partial t} \mathbb{E} \left[ p(X_{\tau}^{\alpha[\lambda, \eta]}) \right] - R^{\text{max}}(\tau) \right) \mu(\tau) d\tau$$

and the Gâteaux differential of W at $\lambda$ is given in the direction $\mu \in L_+^2(0, T)$ by:

$$DW(\lambda, \eta)(\mu) = \int_{0}^{T} \left( \mathbb{E} \left[ p(X_{\tau}^{\alpha[\lambda, \eta]}) \right] - P^{\text{max}}(\tau) \right) \mu(\tau) d\tau$$

**Proof.** Gâteaux derivability is given by adapting Lemma 3.1, in [5]. □

### C. Algorithm

We have in Proposition 4, an expression for the differential that allows us to use a stochastic projected gradient descent algorithm by generating PDMP trajectories to obtain an approximation of the gradient. The step size of the gradient descent is chosen as $\rho_k = \frac{a}{k+1}$, with a chosen empirically to improve the speed of convergence. The time horizon is discretized in $N_t$ points and the subscript $t$ will denote the value of a variable at time $t$.

In this algorithm, we denote the positive part of a number $x$ by $(x)^+$. In the numerical results presented later in the article, Step 5 is the most computationally demanding. Its complexity grows linearly with the number $M$ of generated trajectories. There is therefore a trade-off between choosing a large value of $M$, which ensures accurate Monte Carlo estimation, and keeping the computational cost reasonable. All subsequent results were produced in less than 10 minutes on a standard desktop machine, confirming the practical efficiency of the proposed approach.

### D. Numerical results

A large number of electric vehicles (EVs) arrive at the beginning of the time interval $t = 0$ and leave at the end of the interval $t = 10$. They are modeled by PDMPs $X = (I, S)$,

<figure><img src="image_2.png" /><figcaption>Fig. 1: Normalized global consumptions of $10^4$ WHs for different set of constraints, $N_t = 200$ and $N_s = 40$</figcaption></figure>

### Algorithm 1

1: Initialize $\lambda^0, \eta^0$, set $\{\rho^k\}$ and $M \in \mathbb{N}^*$

2: $k \leftarrow 0$

3: **for** each $k = 0, 1, \dots$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp; $\alpha^k \leftarrow \alpha[\lambda^k, \eta^k]$ where $\alpha[\lambda^k, \eta^k]$ is defined in (6)

5: &nbsp;&nbsp;&nbsp;&nbsp; Generate $M$ PDMPs $X_1^{\alpha^k}, \dots, X_M^{\alpha^k}$

6: &nbsp;&nbsp;&nbsp;&nbsp; $\forall t, U_t^{k+1} \leftarrow \frac{1}{M} \sum_{m=1}^{M} p(X_{m,t}^{\alpha^k}) - P^{\text{max}}(t)$

7: &nbsp;&nbsp;&nbsp;&nbsp; $\lambda_{k+1} \leftarrow (\lambda^k + \rho^k U^{k+1})^+$

8: &nbsp;&nbsp;&nbsp;&nbsp; $\forall t, V_t^{k+1} \leftarrow \frac{1}{M} \sum_{m=1}^{M} \sum_j (p(j, S_{m,t}^{\alpha^k}) - p(X_{m,t}^{\alpha^k}) \alpha_j^k(t, X_{m,t}^{\alpha^k}) + \frac{\partial p}{\partial s}(X_{m,t}^{\alpha^k}) b(X_{m,t}^{\alpha^k}) - R^{\text{max}}(t))$

9: &nbsp;&nbsp;&nbsp;&nbsp; $\eta_{k+1} \leftarrow (\eta^k + \rho^k V^{k+1})^+$

10: **end for**

where $I$ is the mode of charge (1 for ON and 0 for OFF) and $S$ is the state of charge of the battery. We penalize the fact that the vehicles are not sufficiently charged by choosing the final cost $g(i, z) := 30 \times (1 - e^{5(z-0.75)})^+$. The penalization of high intensity is quadratic: $L(x) = x^2$. The charging speed is 0 in mode $i=0$ and

$$b(1, z) = \begin{cases} a, & \text{if } z \le 0.75, \\ a(1-z), & \text{if } z > 0.75, \end{cases}$$

with $a = 0.4h^{-1}$ in mode 1. The cost $c$ models the electricity cost and is high between 2h and 4h and between 6h and 8h, as shown in Fig 1. This creates sharp peaks in the nominal consumption when prices are low, rapidly increasing at 0h and 6h. We want to limit this type of sudden increase by imposing a strong constraint on the slope of consumption. We therefore define $p(X_t^\alpha)$ as the consumption of an EV $X^{\alpha[\lambda, \eta]}$ at time $t$ and the expression $\frac{\partial}{\partial t} E \left[ p(X_t^{\alpha[\lambda, \eta]}) \right]$ is thus the slope of the global consumption at time $t$. We choose $R^{\text{max}} = 0.15/h$, which limits the increase in overall consumption to 0.15 per hour. In the same way, we impose a strong constraint on the decrease in consumption by limiting it to the same value of $0.15/h$. The maximum power is $P^{\text{max}} = 0.2$.

The green curve on Fig. 1(a) shows the result of the algorithm with the Maximum Power constraint. The constraint is met and a part of the global consumption is shifted to the second low-price period (5h-7h30). The green curve in Fig. 1(b) shows the result of the algorithm with the slope constraint,
where global consumption is no longer subject to sharp peaks. In Fig. 1(c), both constraints are implemented. Between Fig. 1(b) and 1 (c), introducing the maximum power constraint shifts a portion of the consumption to later time periods. Since the algorithm seeks to minimize electricity costs, we can observe that in the 3 figures, consumption occurs mainly during low-cost periods. However, the constraints force it to consume outside, particularly on Figures 1(b,c).

## IV. HETEROGENEITY

In this section, we consider that the population is composed of E classes of agents. In the case of EVs, these classes can be seen as different types of vehicles (different battery sizes, different charging power, etc.). The proportion of a class e is denoted $v_e$ and the intensity associated with this class, $\alpha_e$.

### A. Heterogeneous Mean Field Problem

The heterogeneous mean-field control problem is defined as:

$$ \begin{array}{ll} \min_{\{\alpha_e\}_{1 \le e \le E} \in A^E} & \displaystyle\sum_{e=1}^{E} v_e J(\alpha_e) \\ \text{s.t. } X_t^{\alpha_e} \text{ is a PDMP of intensity } \alpha_e & \\ \text{and } \forall t, \displaystyle\sum_{e=1}^{E} v_e \mathbb{E}[p_e(X_t^{\alpha_e})] \le P^{\text{max}}(t) & (7) \\ \text{and } \forall t, \displaystyle\sum_{e=1}^{E} v_e \frac{\partial}{\partial t} \mathbb{E}[p_e(X_t^{\alpha_e})] \le R^{\text{max}}(t) & \end{array} $$

The Lagrangian of this problem $\mathcal{L}^h: A^E \times L_+^2(0,T) \times L_+^2(0,T) \to \mathbb{R}$ is easy to obtain:

$$ \mathcal{L}^h(\{\alpha_e\}_{1 \le e \le E}, \lambda, \eta) := \sum_{e=1}^{E} v_e \mathcal{L}(\alpha_e, \lambda, \eta) $$

where $\mathcal{L}$ is defined in (5), and the associated dual function is $W: L_+^2(0,T) \times L_+^2(0,T) \to \mathbb{R}$:

$$ W^h(\lambda, \eta) := \inf_{\{\alpha_e\}_{1 \le e \le E} \in A^E} \mathcal{L}^h(\{\alpha_e\}_{1 \le e \le E}, \lambda, \eta). \quad (P^h) $$

**Proposition 5.** (i) There exists a unique $(\tilde{\lambda}, \tilde{\eta}) \in L_+^2(0,T) \times L_+^2(0,T)$ maximizing $W^h$ and there is no duality gap.

(ii) For any $\lambda, \eta \in L_+^2(0,T)$, the control $\alpha[\lambda, \eta]$ defined for any $(t, i, j, s) \in [0, T] \times I^2 \times [0, 1]$ by

$$ \alpha_{j,e}[\lambda, \eta](t, i, s) := H' \left( \varphi_e[\eta](t, i, s) - \varphi_e[\eta](t, j, s) \right. \\ \left. - (p_e(j, s) - p_e(i, s))\eta(t) \right) \quad (8) $$

is a minimizer of problem (P$^h$).

(iii) $W^h$ is Gâteaux differentiable. For any $\lambda, \eta \in L_+^2(0,T)$, the Gâteaux differential of W at $\eta$ is denoted by $DW(\lambda, \eta)(\cdot)$ and is given in the direction $\mu \in L_+^2(0,T)$ by:

$$ DW^h(\lambda, \eta)(\mu) = \sum_{e=1}^{E} v_e \int_{0}^{T} \left( \frac{\partial}{\partial t} E \left[ p_e(X_{\tau}^{\alpha_e[\lambda, \eta]}) \right] - R^{\text{max}}(\tau) \right) \mu(\tau) d\tau $$

and the Gâteaux differential of W at $\lambda$ is given in the direction $\mu \in L_+^2(0,T)$ by:

$$ DW^h(\lambda, \eta)(\mu) = \sum_{e=1}^{E} v_e \int_{0}^{T} \left( E \left[ p_e(X_{\tau}^{\alpha_e[\lambda, \eta]}) \right] - P^{\text{max}}(\tau) \right) \mu(\tau) d\tau $$

*Proof.* The Lagrangian is a linear combination of E convex Lagrangians that we have already solved in the previous section. (i), (ii), and (iii) are obtained by rewriting Propositions 1,3 and 4. □

### B. Numerical results

It is easy to adapt Algorithm 1 to take into account the modification of the gradient: Line 5 is modified to generate $M_e$ PDMPs for each class e. and Line 6 and 8 are modified to sum over all classes. Since the computationally expensive step is the gradient evaluation, that must be computed for each class in the previous algorithm, the overall complexity of the algorithm scales proportionally with E.

We consider the same use case as in section III-D but with four class of EVs defined by their charging power: 3,7 kW, 7,4 kW, 11 kW and 22 kW. Electricity cost is set to 0. The number of EVs is $10^4$, equally divided between the classes.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Total consumption and repartition of the consumption depending on the power class when maximum power and ramping constraints are implemented ($N_t = 200$ and $N_s = 40$)</figcaption></figure>

In Fig. 2, the total consumption respects both the maximum power and ramping constraints. It is interesting to note how the consumption of the different classes is distributed. The algorithm first shifts the consumption of the largest consumers (the 22kW EVs) to the end of the day, to allow the smaller ones to consume before. This can be explained by the fact that L penalizes the abrupt switch of EVs from one state to another. The algorithm therefore tries to connect the EVs progressively, but the constraints force it to postpone some of these connections until later in the day. The cost induced by this shift, in the dual, is proportional to power, whereas the price associated with L in the dual is not, which explains why large consumers are shifted more easily. The slight constraint violations in Fig.3 comes from the mean-field approximation. While such discrepancies vanish in the theoretical limit of an infinite population, they may appear in practice due to the
<figure><img src="image_4.png" /><figcaption>Fig. 3: Intensity $\alpha$ to switch from mode 0 to mode 1 for different EV types</figcaption></figure>

finite number of agents. Here, these deviations remain minimal given the large number of agents ($10^4$).

This shift can be seen in the control intensity of the different classes, in Fig 3. For the first 3 classes, control is similar, while for the 22kW class, there is less intensity to switch the EVs to mode ON at the start of the day. For all classes, there is a high-intensity area at the end of the day. This is what forces EVs to charge at the end of the day: an EV that is not sufficiently charged has a high probability of switching to ON mode at that point. On the other hand, there is little or no intensity above 0.75%, as there is a strong penalty $g$ for being below this value.

## V. FUTURE WORKS

A natural direction is to develop a Model Predictive Control (MPC) approach to tackle real-world data, such as the dataset provided in [17]. This approach will build upon the heterogeneous framework, with vehicle arrivals progressively revealed to the algorithm, enabling it to dynamically adjust its charging policy. In such settings, the number of classes is typically larger than in the application considered in this article, and it is necessary to compute a distinct control intensity $\alpha_e$ for each class $e$. This naturally raises the question of how to parallelize the computation associated with the resolution of the numerical scheme in Algorithm 1, to maintain tractability and enable real-time implementation.

A second perspective focuses on the integration of multiple charging modes, which would depend on the class of each agent, and particularly includes Vehicle-to-Grid (V2G) mode, where energy can be returned to the grid. The PDMP-based modeling approach is particularly well suited to capture these class-dependent dynamics, as demonstrated in [5]. Within the heterogeneous framework, it would allow different classes of agents to operate with different modes. This extension has the potential to provide greater flexibility to the power system, especially in highly constrained operational scenarios.

A third research direction would be to enhance the algorithm in two ways. First, by analyzing more carefully the impact of the parameter $M$ on both the computational cost and the estimation quality, with the goal of optimizing the trade-off between accuracy and efficiency. Second, by performing a more systematic study of the problem formulation: the

constraints used in this work were selected without tuning, and a more thorough investigation could reveal which types of constraints or agent dynamics is particularly challenging.

## REFERENCES

[1] Carmona, R., Delarue, F. (2018). *Probabilistic Theory of Mean Field Games with Applications I-II*. Springer.

[2] Achdou, Y., Buera, F.J., Lasry, J.M., Lions, P.L., Moll, B. (2014). *Partial differential equation models in macroeconomics*. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372(2028).

[3] Huang, Malhamé, Caines, *Social optima in mean field LQG control: centralized and decentralized strategies*. IEEE Transactions on Automatic Control, 2012.

[4] M. H. A. Davis, *Markov Models and Optimization*. Chapman and Hall, 1993.

[5] A. Séguret, T. Le Corre, and N. Oudjane, *A decentralized algorithm for a mean-field control problem of piecewise deterministic Markov processes*. ESAIM: Probability and Statistics, 28:22–45, 2024.

[6] P. Garabè and L. Russo, *Ensemble control for large-scale energy management*. IEEE Transactions on Control of Network Systems, 2022.

[7] M. Hochberg, P. R. Kumar, and P. Varaiya, *Demand dispatch: Using real-time measurements to manage electric demand*. IEEE Transactions on Smart Grid, 2006.

[8] A. Bušić and S. Meyn, *Distributed Control of Thermostatically Controlled Loads: Kullback-Leibler Optimal Control in Continuous Time*. Proceedings of the 2019 IEEE 58th Conference on Decision and Control (CDC), 2019.

[9] Callaway, D. S., Hiskens, I. A. (2011). *Achieving controllability of electric loads*. Proceedings of the IEEE, 99(1), 184–199.

[10] Mathieu, J. L., Koch, S., Callaway, D. S. (2013). *State estimation and control of electric loads to manage real-time energy imbalance*. IEEE Transactions on Power Systems, 28(1), 430–440.

[11] Bashash, S., Fathy, H. K. (2011). *Modeling and control of aggregate power of thermostatically controlled loads with demand-side bidding*. Proceedings of the American Control Conference (ACC), 5637–5642.

[12] Huang, Malhamé, Caines *Large population stochastic dynamic games: closed-loop McKean-Vlasov systems and the Nash certainty equivalence principle*. Communications in Information & Systems, 2006.

[13] Achdou, Y., Bardi, M., Cirant, M. (2020). Mean field games models of segregation. *Mathematical Models and Methods in Applied Sciences*, 30(06), 1119–1162.

[14] Bayraktar, E., Budhiraja, A., Cohen, A. (2013). *A stochastic differential game for a multiclass single-server queue*. Mathematics of Operations Research, 38(4), 655–678.

[15] C. Cocozza-Thivent, R. Eymard, S. Mercier, and M. Roussignol, *Characterization of the marginal distributions of Markov processes used in dynamic reliability*, Journal of Applied Mathematics and Stochastic Analysis, 2006.

[16] A. Séguret, *An optimal control problem for the continuity equation arising in smart charging*, J. Math. Anal. Appl., 2023, 127891.

[17] Elaad OpenDataset. *Elaad OpenDataset*, 2019. https://platform.elaad.io/analyses/ElaadNL_opendata.php.