The explosive growth of battery electric vehicle (BEV) adoption presents a dual-faceted challenge and opportunity for modern power systems. While the uncoordinated charging of a large-scale battery electric vehicle fleet can strain grid infrastructure, their collective, controllable charging capacity represents a vast and flexible resource. Effectively harnessing this flexibility necessitates the aggregation of individual battery electric vehicle (BEV) charging profiles into a single, manageable entity for participation in grid services, economic dispatch, and electricity markets. The fundamental mathematical object describing this collective capability is the Aggregate Feasible Region (AFR), which precisely defines the set of all possible aggregated power trajectories the fleet can follow while satisfying each vehicle’s individual constraints.
However, computing the exact AFR, defined as the Minkowski sum of individual feasible regions, is an NP-hard problem. Consequently, researchers have focused on developing approximate methods with tractable computational complexity. Existing approaches broadly fall into two categories: outer approximations and inner approximations. Outer approximations are computationally simple but inevitably include infeasible power profiles within the approximated region, leading to over-optimistic and potentially undeliverable dispatch signals. In contrast, inner approximations guarantee feasibility but tend to be overly conservative, underestimating the true flexibility of the battery electric vehicle (BEV) fleet and leaving valuable grid-supporting potential untapped. The core challenge lies in the highly random and dispersed connection periods of battery electric vehicles, which cause significant heterogeneity in the shape and dimensionality of individual feasible regions, severely degrading the performance of standard approximation techniques.
This article addresses this critical gap. We propose a novel inner-outer combined approximation method specifically designed for the AFR of battery electric vehicles. Our method first employs a dimension-trimming technique to construct a tailored inner approximation for each BEV, significantly improving accuracy. It then leverages an efficient outer approximation at the cluster level to refine the final aggregated region. The result is an approximate AFR that strikes a superior balance between feasibility and conservativeness while maintaining a simple, interpretable Virtual Battery (VB) model form favored by grid operators. Furthermore, to quantitatively evaluate and compare different AFR approximation techniques, we introduce two key performance indicators: Feasibility and Conservativeness. To compute these metrics, we develop an enhanced Hit-and-Run sampling algorithm capable of efficiently handling the high-dimensional polytopes with equality constraints that characterize BEV feasible regions.
Problem Formulation: The Feasible Region of a Battery Electric Vehicle
We consider a fleet of \(N\) battery electric vehicles participating in a coordinated charging program over a scheduling horizon \(\mathcal{K} = [t_0, t_{\text{end}}]\), discretized into intervals of length \(\Delta t\). A charging session for BEV \(i\) is characterized by a tuple \((t_{\text{in}}^i, t_{\text{out}}^i, E_{\text{exp}}^i, P_N^i)\), representing its plug-in time, plug-out time, target energy demand, and rated charging power, respectively. The average charging efficiency is denoted by \(\eta\). The scheduling horizon is partitioned into three segments relative to each BEV’s connection: period \(\mathcal{A}^i=[t_0, t_{\text{in}}^i)\) (before arrival), period \(\mathcal{B}^i=[t_{\text{in}}^i, t_{\text{out}}^i]\) (connected), and period \(\mathcal{C}^i=(t_{\text{out}}^i, t_{\text{end}}]\) (after departure).
The feasible charging region for BEV \(i\), denoted \(\Phi_i\), encompasses all hourly power vectors \(\mathbf{P}^i \in \mathbb{R}^{|\mathcal{K}|}\) that satisfy its operational constraints. It can be formulated as a convex polytope using the half-space representation:
$$\Phi_i = \{ \mathbf{P}^i \in \mathbb{R}^{|\mathcal{K}|} \mid \mathbf{M} \mathbf{P}^i \leq \mathbf{H}^i \}$$
where \(\mathbf{M}\) is a constant matrix defining the constraint normals, and \(\mathbf{H}^i\) is a vector of limits specific to BEV \(i\).
The matrix \(\mathbf{M}\) and vector \(\mathbf{H}^i\) are constructed from power and energy constraints. Let \(\mathbf{P}^i_{\mathcal{A}}, \mathbf{P}^i_{\mathcal{B}}, \mathbf{P}^i_{\mathcal{C}}}\) be the power vectors for periods \(\mathcal{A}, \mathcal{B}, \mathcal{C}\). The power constraints are:
$$
\begin{aligned}
\mathbf{P}^i_{\mathcal{A}} &= \mathbf{0}, \quad \mathbf{P}^i_{\mathcal{C}} = \mathbf{0}, \\
\mathbf{0} &\leq \mathbf{P}^i_{\mathcal{B}} \leq \overline{\mathbf{P}}^i_{\mathcal{B}}, \quad \text{where } \overline{\mathbf{P}}^i_{\mathcal{B}} = P_N^i \cdot \mathbf{1}_{|\mathcal{B}^i|}.
\end{aligned}
$$
The energy dynamics follow \(E_k^i = E_{k-1}^i + \eta \Delta t P_k^i\), with boundary conditions \(E^i(t_{\text{in}})=0\) and \(E^i(t_{\text{out}})=E_{\text{exp}}^i\). This leads to energy envelope constraints that define the fastest and slowest possible charging trajectories. The combined constraints can be written in the standard Virtual Battery (VB) form, a first-order dynamic model akin to a storage device:
$$
\begin{aligned}
\underline{\mathbf{P}}^i &\leq \mathbf{P}^i \leq \overline{\mathbf{P}}^i, \\
\underline{\mathbf{E}}^i &\leq \boldsymbol{\Gamma} \mathbf{P}^i \leq \overline{\mathbf{E}}^i,
\end{aligned}
$$
where \(\boldsymbol{\Gamma}\) is a lower-triangular matrix with non-zero elements \(\eta \Delta t\), and \(\underline{\mathbf{P}}^i, \overline{\mathbf{P}}^i, \underline{\mathbf{E}}^i, \overline{\mathbf{E}}^i\) are the power and energy envelopes. This VB representation is compact and widely accepted for modeling flexible resources.
The precise AFR for the fleet is the Minkowski sum of all individual regions:
$$\Phi_{\text{agg}} = \bigoplus_{i=1}^{N} \Phi_i = \left\{ \mathbf{P} \in \mathbb{R}^{|\mathcal{K}|} \mid \mathbf{P} = \sum_{i=1}^{N} \mathbf{P}^i, \ \mathbf{P}^i \in \Phi_i \right\}.$$
Direct computation of \(\Phi_{\text{agg}}\) is intractable for large \(N\), motivating the need for efficient approximation methods.
Existing Approximation Methods and Their Limitations
Several methods have been proposed to approximate the AFR. Their key characteristics and drawbacks are summarized below, particularly in the context of battery electric vehicle aggregation where connection times are stochastic.
| Method | Type | Basic Principle | Advantages | Limitations for BEVs | Output Model |
|---|---|---|---|---|---|
| Direct Summation [7-9] | Outer | Sums individual constraint envelopes: \(\overline{\mathbf{P}}_{\text{agg}} = \sum_i \overline{\mathbf{P}}^i\), etc. | Extremely fast. Yields a VB model. | Contains significant infeasible region. Accuracy drops with diverse \(\mathcal{B}^i\). | VB Model |
| Homothet-based Inner [16] | Inner | Fits each \(\Phi_i\) inside a scaled/translated baseline polytope \(\Phi_0\). | Guarantees feasible aggregation. Moderate computation. | Highly conservative if \(\Phi_0\) and \(\Phi_i\) have different effective dimensions (due to random \(\mathcal{B}^i\)). | VB Model |
| Structure-Preserving Inner [21] | Inner | Uses a matrix scaling factor to better match individual BEV constraints. | Improves accuracy over simple homothets for BEVs. | Computational complexity is prohibitive for large fleets (\(O(N|\mathcal{K}|^2)\) variables). | VB Model |
| Zonotope/Ellipsoid Methods [13-15] | Inner/Outer | Approximates regions with specific geometric shapes (zonotopes, ellipsoids). | Efficient Minkowski sum operations. | May not align well with BEV polytope structure. Loss of VB model form. | Non-VB Model |
The central issue for battery electric vehicle aggregation is dimensionality mismatch. A global baseline polytope \(\Phi_0\), typically an average of all \(\Phi_i\), has non-zero flexibility across most of \(\mathcal{K}\). An individual BEV’s feasible region \(\Phi_i\), however, is effectively confined to its much shorter connection period \(\mathcal{B}^i\). Attempting to fit a high-dimensional \(\Phi_0\) inside a lower-dimensional \(\Phi_i\) forces the scaling factor to near zero in the homothet method, resulting in a trivial, overly conservative inner approximation. While the structure-preserving method mitigates this, its computational cost is unacceptable for the large-scale battery electric vehicle fleets envisioned for grid services.
The Proposed Inner-Outer Combined Approximation Method
Our method synergistically combines inner and outer approximation steps to overcome the above limitations while maintaining computational efficiency and the desirable VB model form.
Step 1: Establish a Global Baseline Polytope.
We first compute a baseline polytope \(\Phi_0\) representing the “average” BEV in the fleet over the full horizon \(\mathcal{K}\):
$$\Phi_0 = \{ \mathbf{P} \in \mathbb{R}^{|\mathcal{K}|} \mid \mathbf{M} \mathbf{P} \leq \mathbf{H}^0 \}, \quad \text{where} \quad \mathbf{H}^0 = \frac{1}{N} \sum_{i=1}^{N} \mathbf{H}^i.$$
Here, \(\mathbf{H}^0\) contains the average power and energy envelopes of the fleet.
Step 2: Dimension-Trimmed Maximum Inner Approximation for Each BEV.
This is the core innovation to handle random connection periods. For each battery electric vehicle \(i\), we tailor the baseline to its specific active dimensions.
- Dimension Trimming: We construct a diagonal “trimming matrix” \(\mathbf{C}^i\) that zeros out entries in \(\mathbf{H}^0\) corresponding to periods when BEV \(i\) is not connected or has no energy flexibility. This yields a trimmed baseline polytope \(\Phi_0^i\):
$$\Phi_0^i = \{ \mathbf{P} \in \mathbb{R}^{|\mathcal{K}|} \mid \mathbf{M} \mathbf{P} \leq \mathbf{C}^i \mathbf{H}^0 \}.$$
Crucially, \(\Phi_0^i\) now has the same effective dimensionality (non-zero constraints) as \(\Phi_i\). - Maximum Inner Approximation: We find the largest homothet of the trimmed baseline \(\Phi_0^i\) that fits inside BEV \(i\)’s true feasible region \(\Phi_i\). This involves solving a small linear program for each BEV:
$$
\begin{aligned}
&\underset{\alpha_i, \boldsymbol{\beta}_i}{\text{maximize}} && \alpha_i \\
&\text{subject to} && \alpha_i \geq 0, \\
&&& \alpha_i \mathbf{C}^i \mathbf{H}^0 + \mathbf{M} \boldsymbol{\beta}_i \leq \mathbf{H}^i.
\end{aligned}
$$
Let \((\alpha_i^*, \boldsymbol{\beta}_i^*)\) be the optimal solution. The maximum inner approximation for BEV \(i\) is:
$$\Phi_{i,\text{inner}} = \alpha_i^* \Phi_0^i + \boldsymbol{\beta}_i^* = \{ \mathbf{P} \mid \mathbf{M} \mathbf{P} \leq \alpha_i^* \mathbf{C}^i \mathbf{H}^0 + \mathbf{M} \boldsymbol{\beta}_i^* \}.$$
Because \(\Phi_0^i\) aligns with \(\Phi_i\)’s dimensions, \(\alpha_i^*\) is typically positive, yielding a non-conservative inner approximation.
Step 3: Efficient Outer Approximation for Cluster Aggregation.
With the set of inner-approximated regions \(\{\Phi_{i,\text{inner}}\}\), we compute an outer approximation of their Minkowski sum. Since all \(\Phi_{i,\text{inner}}\) share the same constraint matrix \(\mathbf{M}\), their sum’s outer approximation is simply the sum of their right-hand sides:
$$\Phi_{\text{agg, in-out}}^* = \left\{ \mathbf{P} \mid \mathbf{M} \mathbf{P} \leq \sum_{i=1}^{N} \left( \alpha_i^* \mathbf{C}^i \mathbf{H}^0 + \mathbf{M} \boldsymbol{\beta}_i^* \right) \right\}.$$
This operation is computationally trivial—just a summation of vectors.
Step 4: Final Refinement with Exact Energy Envelopes.
We observe that the energy envelope constraints from the simple outer approximation (Direct Summation) are always feasible for the true AFR. Therefore, we replace the energy-related part of the constraint vector in \(\Phi_{\text{agg, in-out}}^*\) with the exact sum of individual energy envelopes from Step 1. The final approximated AFR, \(\Phi_{\text{agg, in-out}}\), retains the VB model form:
$$
\boxed{\begin{aligned}
\underline{\mathbf{P}}_{\text{agg}}^* &\leq \mathbf{P} \leq \overline{\mathbf{P}}_{\text{agg}}^* \\
\underline{\mathbf{E}}_{\text{agg}} &\leq \boldsymbol{\Gamma} \mathbf{P} \leq \overline{\mathbf{E}}_{\text{agg}}
\end{aligned}}
$$
where \(\underline{\mathbf{E}}_{\text{agg}} = \sum_i \underline{\mathbf{E}}^i\) and \(\overline{\mathbf{E}}_{\text{agg}} = \sum_i \overline{\mathbf{E}}^i\).
This method inherently performs a “soft-clustering”: BEVs with similar connection periods \(\mathcal{B}^i\) are approximated using similarly trimmed baselines, and their inner approximations aggregate precisely via the homothet property before the final outer summation. The computational complexity is linear in the number of battery electric vehicles, \(O(N)\), making it highly scalable.
Quantitative Evaluation Metrics and Sampling Methodology
To objectively compare AFR approximation methods, we propose two metrics based on Monte Carlo sampling. Let \(\Phi_{\text{agg}}\) denote the true (but unknown) AFR and \(\widehat{\Phi}_{\text{agg}}\) an approximation of it.
1. Feasibility (\(F\)): The probability that a randomly sampled point from the approximated region is actually feasible for the true aggregate.
$$F(\widehat{\Phi}_{\text{agg}}) = \frac{1}{|\mathcal{S}_{\widehat{\Phi}}|} \sum_{\mathbf{P} \in \mathcal{S}_{\widehat{\Phi}}} \mathbb{I}(\mathbf{P} \in \Phi_{\text{agg}})$$
where \(\mathcal{S}_{\widehat{\Phi}}\) is a uniform sample set from \(\widehat{\Phi}_{\text{agg}}\), and \(\mathbb{I}(\cdot)\) is the indicator function. For an outer approximation, \(F \leq 1\); a higher \(F\) is better.
2. Conservativeness (\(C\)): The probability that a randomly sampled point from the true region is contained within the approximation.
$$C(\widehat{\Phi}_{\text{agg}}) = \frac{1}{|\mathcal{S}_{\Phi}|} \sum_{\mathbf{P} \in \mathcal{S}_{\Phi}} \mathbb{I}(\mathbf{P} \in \widehat{\Phi}_{\text{agg}})$$
where \(\mathcal{S}_{\Phi}\) is a uniform sample set from \(\Phi_{\text{agg}}\). For an inner approximation, \(C \leq 1\); a higher \(C\) is better (less conservative).
| Approximation Type | Feasibility (F) | Conservativeness (C) | Ideal Target |
|---|---|---|---|
| Outer Approximation | \(\leq 1\) (Higher is better) | \(= 1\) | Maximize \(F\) |
| Inner Approximation | \(= 1\) | \(\leq 1\) (Higher is better) | Maximize \(C\) |
| Proposed Combined Method | \(\leq 1\) (should be > Outer’s F) | \(\leq 1\) (should be > Inner’s C) | Balance high \(F\) and high \(C\) |
Computing these metrics requires solving two challenges: uniformly sampling high-dimensional polytopes (\(\Phi_{\text{agg}}\) and \(\widehat{\Phi}_{\text{agg}}\)) and evaluating the indicator function \(\mathbb{I}(\mathbf{P} \in \Phi_{\text{agg}})\).
Enhanced Hit-and-Run Sampling: The standard Hit-and-Run algorithm performs poorly for BEV polytopes due to equality constraints (e.g., total energy demand) and “sharp corners.” Our enhancement modifies the direction-finding step. At a current point \(\mathbf{P}_k\), let \(\mathcal{U}\) be the set of “active” inequality constraints (within threshold \(\epsilon\)). We generate a random direction \(\mathbf{d}\) by solving:
$$
\begin{aligned}
&\underset{\mathbf{d}}{\text{minimize}} && ||\mathbf{d} + \boldsymbol{\rho}||_2 \\
&\text{subject to} && \mathbf{M}_j’ \mathbf{d} \leq 0, \quad \forall j \in \mathcal{U} \\
&&& \mathbf{Q}_l^T \mathbf{d} = 0, \quad \forall l \in \{\text{equality constraints}\}
\end{aligned}
$$
where \(\boldsymbol{\rho}\) is a random vector, \(\mathbf{M}_j’\) are the normals of active inequalities, and \(\mathbf{Q}_l\) are the normals of equalities. This ensures \(\mathbf{d}\) points inward and remains parallel to equality hyperplanes, enabling efficient exploration.
Feasibility Check: To evaluate \(\mathbb{I}(\mathbf{P} \in \Phi_{\text{agg}})\) for a candidate aggregate trajectory \(\mathbf{P}\), we solve a Constraint Satisfaction Problem (CSP) that checks if \(\mathbf{P}\) can be decomposed into feasible individual BEV schedules:
$$
\begin{aligned}
&\text{find} && \{\mathbf{P}^1, \mathbf{P}^2, …, \mathbf{P}^N\} \\
&\text{subject to} && \mathbf{M} \mathbf{P}^i \leq \mathbf{H}^i, \quad \forall i=1,…,N \\
&&& \sum_{i=1}^{N} \mathbf{P}^i = \mathbf{P}.
\end{aligned}
$$
If a feasible solution exists, the indicator is 1; otherwise, it is 0.
Numerical Case Studies and Analysis
We validate the proposed method using three real-world datasets representing distinct battery electric vehicle charging profiles: workplace charging (Caltech ACN [24]), electric bus depot charging, and residential charging. The scheduling horizon is 24 hours (\(\Delta t = 1\) hr). We compare our Inner-Outer Combined method against: 1) Direct Summation (Outer), 2) Homothet-based Inner Approximation, and 3) Structure-Preserving Inner Approximation.
Visual Comparison of Approximated AFRs
The power and energy envelopes of the approximated AFRs for a 20-BEV sample from the workplace dataset clearly illustrate the trade-offs. The Homothet-based method yields the narrowest power envelope, indicating high conservativeness, with flexibility limited primarily to midday hours. The Structure-Preserving method improves upon this. Our proposed method produces power envelopes that are almost as wide as the optimistic Direct Summation outer approximation, while we know from its construction that it contains far fewer infeasible points. All methods share the same (feasible) energy envelopes from the outer approximation in the final step.
Quantitative Performance over Multiple Days
We evaluate the metrics over five consecutive days for each dataset (N=20 per day). The table below shows average results for the workplace dataset.
| Method | Avg. Feasibility (F) (%) | Avg. Conservativeness (C) (%) | Avg. Computation Time (s) |
|---|---|---|---|
| Direct Summation (Outer) | 74.2 | 100.0 | <0.01 |
| Homothet-based Inner | 100.0 | ~0 (highly conservative) | ~0.05 |
| Structure-Preserving Inner | 100.0 | 45.7 | ~9.0 |
| Proposed Inner-Outer Combined | 99.8 | 78.3 | ~0.06 |
Our method achieves near-perfect feasibility (significantly higher than the Outer method) while maintaining a high conservativeness score (substantially higher than the Structure-Preserving Inner method). Crucially, it does this at a computational cost comparable to the simple Homothet method and orders of magnitude faster than the Structure-Preserving method.
Stability Under Heterogeneous Connection Times
To test robustness, we simulated a simple two-BEV scenario, varying the plug-in time of the second battery electric vehicle. The performance metrics as a function of the overlap in their connection periods reveal key insights. The feasibility of both the Outer and our method dips slightly when connection periods partially overlap but remains high. Most strikingly, the conservativeness of the Homothet method falls to zero as soon as connection periods differ, and the Structure-Preserving method’s conservativeness also declines steadily. Our proposed method demonstrates the most stable and gradual decline in conservativeness, showcasing its superior adaptability to the inherent randomness in battery electric vehicle behavior.
Scalability Analysis
We assess computation time as the fleet size \(N\) increases. The results are clear: the computational burden of the Structure-Preserving method grows rapidly, becoming impractical (over 500 seconds for N=40). In contrast, the computation time for both the Homothet method and our proposed Inner-Outer Combined method scales linearly and remains negligible (around 2 seconds for N=1000), confirming the high scalability essential for real-world battery electric vehicle aggregation.
Sampling Effectiveness
Applying our enhanced Hit-and-Run sampler to a high-dimensional AFR (24-dimensional for a 24-hour schedule) generates a diverse set of power trajectories that uniformly fill the envelope constraints. A 2-dimensional projection of the samples confirms their uniform distribution within the feasible space, validating the reliability of the Monte Carlo estimates used for the \(F\) and \(C\) metrics.
Conclusions and Future Work
This article addresses the critical problem of approximating the aggregate feasible region for large-scale battery electric vehicle fleets. The random and dispersed connection times of battery electric vehicles render traditional approximation methods either overly optimistic or excessively conservative. We propose an inner-outer combined approximation method that introduces a dimension-trimming technique to build accurate inner approximations tailored to each BEV’s active charging period, followed by an efficient outer aggregation. The resulting AFR maintains a simple Virtual Battery model form, facilitating its integration into standard grid dispatch and market frameworks.
To enable quantitative comparison, we defined two key metrics—Feasibility and Conservativeness—and developed an enhanced Hit-and-Run sampling algorithm capable of handling the high-dimensional, equality-constrained polytopes representative of BEV aggregation. Numerical simulations using real-world charging data from workplace, bus depot, and residential settings demonstrate that our method consistently outperforms existing approaches. It achieves feasibility close to 100% (significantly better than outer approximations) and conservativeness above 75% (markedly better than inner approximations), all while maintaining linear computational complexity and scalability to thousands of battery electric vehicles.
Future work will focus on extending this deterministic framework to account for the stochastic nature of BEV arrival, departure, and energy demand. Integrating these uncertainties into a probabilistic Virtual Battery model would provide aggregators and system operators with a more robust characterization of the available battery electric vehicle flexibility, further enhancing the reliability of grid services provided by this burgeoning distributed resource.