Dynamic PSA-CEM for EV Charging Pricing

• The paper proposes a novel hybrid algorithm integrating dynamic probabilistic sensitivity analysis with CEM to efficiently address high-dimensional, bilevel EV charging pricing challenges.
• It employs a rolling-horizon framework and the Method of Successive Averages to manage temporal dependencies and user equilibrium within a robust stochastic optimization setting.
• Empirical results demonstrate improved queue management and pricing performance in EV charging stations, highlighting the method’s scalability and adaptability.

The Dynamic Probabilistic Sensitivity Analysis-guided Cross-Entropy Method (PSA-CEM) is an adaptive, high-dimensional stochastic optimization approach developed for bilevel, nonconvex, and behaviorally heterogeneous dynamic pricing problems, as exemplified by electric vehicle charging systems. Integrating the Cross-Entropy Method (CEM) with a dynamic probabilistic sensitivity analysis, PSA-CEM targets efficient global optimization in problems characterized by vast decision spaces and intricate lower-level stochastic equilibria. In combination with rolling-horizon decomposition and the Method of Successive Averages (MSA), PSA-CEM enables tractable, scalable solution of time-coupled pricing strategies under queueing and multinomial logit (MNL) user choice models (Zhang et al., 20 Jan 2026).

1. Optimization Problem Structure

The target application involves dynamic pricing of EV charging, with the decision vector $$\boldsymbol\theta=(p_{11},p_{12},\dots,p_{N_sN_h})\in\Omega\subset\mathbb R^d$$ representing all station-hour price variables. The upper-level objective function to be maximized is

$$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$

subject to box constraints $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$. Lower-level user flows and queueing effects are resolved by an MNL model augmented with queuing-theoretic approximations. This construction yields a bilevel, high-dimensional stochastic program that is intractable to solve directly without specialized metaheuristics.

2. Dynamic Probabilistic Sensitivity Analysis Integration

PSA is interleaved into CEM iterations. At each iteration $\ell$, candidate price vectors $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$ are sampled from a parameterized density $f^{(\ell)}$, and their performance evaluated. To quantify the influence of each decision variable $\theta_k$ on the distribution of $F$, a “frozen” sample set is constructed for $k$ by fixing $\theta_k$ at its elite-mean $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$0 while retaining all other sample components, and recomputing the induced pdf $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$1 for $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$2. The probabilistic sensitivity index is given by the Kullback–Leibler divergence: $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$3 where $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$4 is the empirical density under normal sampling. The set of indices with $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$5 (with threshold $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$6) defines the active set $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$7. Only $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$8 undergoes adaptive distribution parameter updates, freezing all others, enabling dimensionality reduction and focused optimization.

3. Cross-Entropy Method Formulation

Each CEM iteration proceeds by sampling $$F(\boldsymbol\theta)=\sum_{t=1}^{N_h}\sum_{i=1}^{N_s}\sum_j Q_{ij}^\Delta(t)(p_i(t)-\nu_i(t))p_{ij}(t) + (1-\omega)[\sum_{j}\sum_{t,i}\kappa Q_{ij}^\Delta(t) p_{ij}(t) - \nu\sum_{t,i,j}T_{ij}^q(t)p_{ij}(t) -\eta\sum_{t,i}\lambda_i(t)\pi_{i,c_i}(t)]$$9 candidate $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$0 from a multivariate normal with independent marginals,

$\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$1

Performance scores $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$2 are computed for each sample, and the top $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$3 fraction (e.g., $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$4) by $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$5 value comprise the elite set $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$6. Update equations for means and variances of the elite samples are

$\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$7

For $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$8, the parameter update with smoothing $\nu_i(t)\leq p_i(t)\leq p_{\max,i}(t)$9 is

$\ell$0

while variables not in the active set remain frozen. This regime promotes rapid convergence on sensitive decision variables and avoids unnecessary update noise in insensitive ones.

4. Rolling-Horizon Integration

PSA-CEM is embedded in a rolling-horizon procedure, essential for addressing temporal dependencies of queue states and time-varying user demand. The overall optimization horizon $\ell$1 is divided into windows of length $\ell$2 (typically one hour). For each window:

• Queue-states and initial prices are set.
• PSA-CEM is run to convergence (stopping after a relative elite-mean change $\ell$3 for 2 consecutive iterations).
• The optimized period's prices are fixed, queue transitions are computed via MSA, and the window is advanced, with updated queue states and sampling parameters carried forward.

This receding-horizon decomposition permits large-scale, multi-period optimization, accommodating system dynamics and state carryover between windows without exponential computational cost growth.

5. Algorithmic Steps and Update Logic

The PSA-CEM algorithm can be summarized as follows:

1. Initialize $\ell$4, $\ell$5 for each $\ell$6; set iteration $\ell$7.
2. For each rolling window:
   • Set queue-states using the previous solution.
   • Iterate until convergence:
     • Draw $\ell$8 samples from current parameterization.
     • Compute lower-level equilibrium (MNL+MSA) and $\ell$9 for each.
     • Form elite set of the best $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$0 samples.
     • Every $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$1 iterations, recalculate sensitivity indices $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$2 and active set $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$3.
     • Update $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$4, $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$5 for $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$6; freeze others.
     • Check for convergence in variances or elite-mean relative change.
   • Fix the period's optimal prices and advance queue states.

6. Implementation Parameters and Practical Considerations

Typical parameterization for effective operation:

Parameter	Typical Value	Role
Population size $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$7	$\{\boldsymbol\theta^{(m)}\}_{m=1}^N$8	Number of samples per iteration
Elite fraction $\{\boldsymbol\theta^{(m)}\}_{m=1}^N$9	$f^{(\ell)}$0	Fraction for elite selection
Smoothing $f^{(\ell)}$1	$f^{(\ell)}$2	Balance update/new estimate
PSA update $f^{(\ell)}$3	$f^{(\ell)}$4 iterations	Sensitivity recalculation frequency
Sensitivity $f^{(\ell)}$5	$f^{(\ell)}$6	PSA threshold for activeness
Convergence $f^{(\ell)}$7	$f^{(\ell)}$8	Elite mean relative change tolerance
Rolling window $f^{(\ell)}$9	$\theta_k$0 h	Horizon per optimization window

Auxiliary implementation recommendations include the use of Gaussian moment-matched densities for pdf fitting, variance bounds to prevent collapse, parallelization of batch evaluations, and warm-starting of $\theta_k$1 and queue-states between windows. The lower-level equilibrium in each sample is efficiently resolved with MSA using a diminishing stepsize $\theta_k$2 (Zhang et al., 20 Jan 2026).

7. Computational Impact and Application Scope

The PSA-CEM method offers dimensionality screening and targeted search, enabling solution of high-dimensional pricing problems with high fidelity to behavioral and system uncertainties. Empirical studies in the context of EV charging networks show marked improvements over fixed and time-of-use pricing, both in user utility and queuing performance, as evidenced in real-world trials with 22 stations. The adaptive freezing of insensitive variables, coupled with rolling-horizon decomposition, ensures both computational efficiency and robustness. A plausible implication is applicability to a broader class of bilevel stochastic programs with similar structural properties, especially where sensitivity varies across dimensions and temporal coupling is significant (Zhang et al., 20 Jan 2026).