Electric Vehicle Charge Scheduling on Highway Networks from an Aggregate Cost Perspective

TL;DR

This paper presents an optimal scheduling method for electric vehicle charging on highways to minimize total system costs, considering both user and infrastructure perspectives using a hybrid systems approach.

Contribution

It introduces a hybrid systems based formulation for joint optimization of EV charging schedules from multiple stakeholder viewpoints.

Findings

Effective scheduling reduces overall system costs.

Balances interests of EV users and infrastructure providers.

Provides a framework for integrated EV charging management.

Abstract

In this paper, we attempt to optimally schedule the charging of long-range battery electric vehicles (BEVs) along highway networks, in order to minimize aggregate costs to the overall system consisting of utilities or electricity providers, station operators and other infrastructure, as well as EV users. Thus, we approach the problem from the perspective of both customers (EV car owners), as well as charging station operators and utilities using a hybrid systems based formulation.

Tables

Table 1. Table I: Definitions of general notation

Variable and Symbol	Notes
Car index $c\in \{1,\dots \dots ,p\}$	p = Total no. of cars
Time step $k\in \{1,\dots ..,H_p\}$	$H_p$ = Prediction horizon
Station index $n\in \{1,\dots .N\}$	N = Total no. of stations
${\nu }_0$	Starting node for a car
${\nu }_{H_p}$	Ending node for a car
$d_k$	Discrete trip distance counter
${ϵ}_k$	Discrete local (edge) distance counter
$\gamma \in \{0,1\}$	Constrained input: at a node or on an edge
$y\in \{0,1\}$	Constrained input: charging or not
$E_k$	Energy level of an EV at time k

Equations

$$C_{electricity}$$

$$\displaystyle\begin{split}C_{station}&=\sum_{k=1}^{H_{p}}\sum_{n=1}^{N}\bar{C}\;*\mid sgn\bigg{(}\frac{U(k,n)}{S_{n}}-0.5\bigg{)}\mid*\\ &\qquad\;(U(k,n)-S_{n})\end{split}$$

$$U(k,n)$$

$$\displaystyle sgn\bigg{(}\frac{U(k,n)}{S_{n}}-0.5\bigg{)}=\begin{cases}1,&U(k,n)>50\%\;of\;S_{n}\\ 0,&U(k,n)=S_{n}\\ -1,&U(k,n)<50\%\;of\;S_{n}\end{cases}$$

$$E_{k+1}$$

$$\sum _{i=1}^3{\beta }_k(i)=1$$

$$d_{k+1}^c$$

$$\sum _{i=1}^3{\beta }_k(i)=1$$

$${ϵ}_{k+1}^c$$

$$\sum _{i=1}^3{\beta }_k(i)=1$$

$$C_{customer}:=$$

$$+C^{degrade}(\Delta E)+C^{t_{waiting}}$$

$$\displaystyle C^{t_{driving}}(e^{ij}):=t_{e^{ij}}=t_{e^{ij}}^{0}*\bigg{(}1+0.15*\bigg{(}\frac{v_{e^{ij}}^{k}}{c_{e^{ij}}}\bigg{)}\bigg{)}$$

$$\overrightarrow{v_{e^{ij}}^k}=\sum _{c=1}^p{\xi }_k^c(e^{ij})$$

$$C^{t_{charging}}(e^{ij}):=\{\begin{array}{l}\frac{(E_c-E_i)}{P_{max}}+\frac{1}{n}log(\frac{m-n{E}_c}{m-n{E}_f}),E_i<E_c\\ \frac{1}{n}log(\frac{m-n{E}_i}{m-n{E}_f}),else\end{array}$$

$$P(E):=\{\begin{array}{l}{P}_{max},E<E_c\\ m-n\cdot E,else\end{array}$$

$$\displaystyle C_{degrade}:=\bigg{(}a(\mathbb{E}(P_{ow}))^{2}+b*\mathbb{E}(P_{ow})+c\bigg{)}\mathbb{E}(t_{charging})$$

$$C_{waiting}:=f(\lambda ,\nu )$$

$${\delta }_k^e$$

$$x_{k+1}^l$$

$$i_k$$

$$x_{k+1}^{ct}$$

$${\delta }_{k+1}^e(1)=\{\begin{array}{ll}1,& {ϵ}_k\ge e_k^{ij}\\ 0,& {ϵ}_k\in [0,e^{ij})\end{array}$$

$${\delta }_{k+1}^e(2)=\{\begin{array}{ll}1,& {ϵ}_k\in (0,e^{ij})\\ 0,& {ϵ}_k\notin (0,e^{ij})\end{array}$$

$${\gamma }_{k+1}^l=\{\begin{array}{ll}1,& {\delta }_k^e(1)\\ 0/1,& \neg {\delta }_k^e(1)\\ 0,& {\delta }_k^e(2)\end{array}$$

$$y_{k+1}^l$$

$$\sum _{c=1}^p{y}_{k+1}^{ct}(c)\le ca{p}_{k+1}^n$$

$${\xi }_{k+1}^l(h)$$

$$\sum _{h=1}^{N^2}{\xi }_k^l(h)$$

$${\mu }_{k+1}^l$$

$$\sum _{i=1}^3{\beta }_k(i)$$

$$e_{k+1}$$

$$\underset{\{y,\gamma \in \{0,1{\}}^{H_p},\xi \in \{0,1{\}}^{H_p,N^2}{\}}^p}{\min }{C}_{total}$$

$$C_{total}=C_{stations}+C_{customers}+C_{electricity}$$

$$\text{eqref}eq:mlda$$

$$\underline{x}\le x_k\le \bar{x},\forall c,\forall k\in \{0,1,...H_p\}$$

$${ϵ}_0=0,\forall c,\forall k\in \{0,1,...H_p\}$$

$${ϵ}_{H_p}=0,\forall c,\forall k\in \{0,1,...H_p\}$$

Full text

I Introduction

I-A Motivation

The popularity of electric vehicles has been rising rapidly in recent years, driven mainly by the greater driving comfort, improved safety, cheaper maintenance and lower exhaust emissions they offer. As battery and materials costs continue to decline and we transition towards cleaner, more diversified and renewable electricity generation, the use of both plug-in hybrid and battery electric vehicles is likely to expand even more.

Increased EV utilization rates also mean that we will need more sophisticated and effective techniques to schedule the charging of both plug-in and battery electric vehicles in real-time, and allocate charging capacity or power in the most efficient manner. This is in order to minimize negative effects like congestion, long waiting times for customers and high infrastructure costs like uneven demands or loads for utilities and station operators. On a broader level, optimal EV charging will also play an important role in realizing smart grids of the future.

I-B Literature Review

The work proposed thus far for EV scheduling has been focused on four areas: optimization with respect to load on power distribution networks [5],[9],[12],[13],[14], queuing models for input-output dynamics of charging stations [1],[2],[6],[17], game-based approaches for trajectory optimization [2],[4],[8], and more traditional trajectory optimal control methods [3],[7],[10],[11],[15]. Note that power distribution problems either use game-based approaches or traditional optimal control methods but are considered distinct in their objectives: the former looks at power quality on networks whereas the latter has objective functions related to time or expected financial costs.

The work done in power distribution considers the effect of EV charging on the grid, mainly with regard to base load. In this sense EVs are treated as dispatchable loads that can fill valleys during otherwise low demand. Other work has considered topological effects of charging station power draws. The work described in [5] uses statistical estimation to predict future EV arrivals and the associated charging demands in order to perform load-shifting via model predictive control. [9] proposes a mean-field game-based approach for large EV populations. The objective is to fill the overnight valley in base load and does not consider EV routing. The work in [12] coordinates plug-in hybrid EVs (PHEV) in order to reduce system power losses, maximize load factor, and minimize load variance. A decentralized approach is proposed in [13] that focuses on load shifting via dispatch signals provided by a centralized authority. The issue of load shifting is considered in [14] with emphasis on circumventing load forecasting limitations. A locally optimal distributed approach is demonstrated. The work mentioned here focuses on optimization in the time domain–load shifting for valley filling throughout a day.

Many different queuing model approaches have been proposed for optimal EV scheduling. These works generally focus on modeling the arrival and departure/service rates of the EVs. While some approaches utilize probability theory as the foundation for the scheduling, others rely on optimization. The work in [1] presents a first-come first-served queuing model with arrival rates corresponding to an exponential distribution and Poisson departure rates in a Markov chain. This does not account for location and path planning. The dynamics of the EV battery are not considered. A similar approach is shown in [2] where the objective is to minimize EV waiting time at charging stations. This makes more considerations around path planning between nodes but limits the road network to a one-dimensional, unidirectional flow. [6] also minimizes EV waiting time using a similar queuing model. Theoretical bounds on waiting time are presented as well as a communication strategy (vehicle-to-charging station). A battery-swapping station approach in [17] utilizes queues and also defines probabilities around customer dissatisfaction.

Trajectory optimization problems have been formulated for electric vehicle charging with an emphasis on finding efficient algorithms and sensible objective functions. In [3] a hierarchical approach is proposed that considers power system network dynamics and battery dynamics. Given a fleet of EVs the highest level problem is station location. The provisioning problem sits below this in which the number of chargers per station is allocated according to a queuing model and the given demand. The final layer consists of a global cost minimization. The work in [7] focuses on a graph theoretic framework for EV scheduling. The work exploits previous classical work done in path planning to reduce the complexity of scheduling EV charging in large networks. [10] details heuristic search algorithms for solving the EV-routing problem on a road network with charging stations located at specific nodes. A mixed approach is described in [11] whereby probabilistic tools, power distribution concerns, and an objective function formulated for routing efficiency are used. More specifically, the problem maximizes the station profits while maintaining a high equality of service. Finally, the work in [15] utilizes finite horizon and receding horizon approaches to consider the two cases of 1) a priori knowledge of all users and 2) new inputs (vehicles) to the system during the horizon for vehicle routing on a charging network. The objective function takes into consideration the customer’s perspective in an attempt to satisfy demand while keeping the cost low. Here a variable charging rate is implemented, which adds complexity to the required infrastructure and solution space but enables more interesting solutions.

II Contributions

In this paper, we propose looking at the scheduling problem for EV charging from the perspective of minimizing aggregate costs to the system as a whole. This addresses an important null space in the literature since much of the previous work in this area has focused on optimizing the allocation of charging power with respect to only one of the players i.e. utilities, charging station operators, or customers.

A lot of the past research has also focused on determining optimal locations for charging stations in a given area. However, in this paper, we take the station locations (i.e. the nodes) as a fixed parameter for a particular city or urban area’s highway network. Instead, we aim to design a central controller or algorithm that determines the optimal trajectory that an EV can take given their starting position and desired destination, in a manner that maximizes benefits to all stakeholders involved.

Furthermore, a great deal of effort has gone into minimizing losses to utilities and the grid - and utilizing EVs to provide services like load flattening, voltage and frequency regulation. However, we place more emphasis on the opportunity costs paid by customers as well, in terms of both time spent and physical degradation of their car batteries. Following the prescribed route will also help alleviate range anxiety among users, and potentially contribute towards increased EV penetration.

Additionally, this work defines a hybrid systems theoretic model and uses finite state machines to describe the dynamics of electric vehicles on a highway network under two different frameworks. The model is first described under a framework that is more easily converted into code. The second description allows for more thorough analysis of the model’s properties.

III Model

The following model definitions and formulations formalize the Mixed Integer Quadratic Programming (MIQP) problem.

III-A Definitions and Notation

Table I defines the basic terms that are used throughout the work. The remainder of the notation is defined in-text when first introduced.

III-B Key Assumptions of Formulation

• •

All nodes or intersections in our highway network (represented as a graph) have a charging station.

• •

All EVs in our network have identical battery energy level capacities and the same discharging power dynamics while driving.

• •

All charging stations in our network provide identical identical charging power profiles for all EV customers.

III-C Highway Network Model

A highway network can be described as a connected graph $\mathcal{G}=(\mathcal{V},\mathcal{E},A)$ where $\mathcal{V}=\{1,..,N\}$ is a non-empty set of N nodes that represent road intersections and the potential location of a EV charging station, $\mathcal{E}$ is the set of edges representing the road segment lengths, and $A$ is the adjacency matrix where $a^{ii}=0$ and $a^{ij}>0$ for all connected nodes $i\neq j$. The adjacency matrix defines whether nodes are connected. For the set of edges, $\mathcal{E}$, an entry $e^{ij}>0$ indicates the transition from node $i$ to node $j$ with a positive weight. Generally, a graph can be directed or undirected. Considering that the graph represents a highway, traffic is assumed to be able to flow both ways without loss of generality. This can be generalized to cases with only one-way traffic flows, because if $\exists$ $e^{ij}>0$ then and physically $\nexists$ $e_{ji}>0$, then $e_{ji}:=0$ where the convention is taken to be that only positive-valued edges can be traversed. Equivalently, the zero-valued edges can be assigned to $NaN$ or any placeholder value. Finally, the neighborhood of a node is defined where $\mathcal{N}_{i}:=\{j|(i,j)\in\mathcal{E}\text{ s.t. }e^{ij}>0\}$.

The graph representation helps to formalize the work in notation but no novel or extensive graph theory is used here.

The network used in the simulations is presented in figure 1. It is included here to help motivate the application of graph theory to the EV network.

III-D Costs for the Network

[TABLE]

Let $P^{elec}_{k}(n)$ represent the electricity price (dollars/kWh) at a particular station n, as a function of location in the network and according to the time step k. The difference in energy when it is charging (i.e. when $y_{k}$ is true) is taken to be the amount of electricity, in kWh, that was just injected into the battery. This cost term then is in terms of dollars.

We originally considered including electricity costs using Time-of-Use (TOU) rates as part of stations’ costs, as described in equation (III-D) above. However, we later decided against this since previous studies [8] have shown that electricity costs are much smaller than costs associated with travel times. Thus, we chose time minimization for customers as the dominant decision criterion, where travelers choose the fastest feasible path while also satisfying the other constraints.

[TABLE]

In order to avoid either excessive under or over-utilization of charging stations and ensure a uniform burden on all nodes across the network, we decided to penalize stations whenever they use above or below 50% of their maximum rated capacity, as shown in (3). We imposed a flat linear cost $\bar{C}$ on the difference between the current station utilization and the preferred capacity of that station $S_{n}$. For simplicity, we defined station utilization $U(k,n)$ in terms of the number of cars charging at a particular station $n$ at any given time $k$. This can be easily computed by summing over the boolean variables $y_{k}^{c}(n)$ in our logic FSM described in Section IV-C, as shown in (3) above. In general the following definition holds:

[TABLE]

III-E Electric Vehicle Model

The EV model is defined as three-mode finite-state machine (FSM) (Fig. 2). Namely, there is a charging mode, driving mode, and a waiting mode (at the charging station). The three piecewise-continuous states are the energy in the battery, the trip distance as a scalar value, and the current edge distance traversed.

While less accurate than state of charge, the energy in the battery can be manipulated much more easily and with enough accuracy for the scheduling of trajectories. The dynamics of the battery are then a simplified version of Coulomb counting where the change in energy is a function of the power injected or consumed. Given $\frac{\partial E}{\partial t}:=P(t)$ as the most basic form of Coulomb counting, let it be discretized under the one-step Euler to $\Delta E:=E_{k+1}-E_{k}=P_{k}*T_{s}$, such that $E_{k+1}=E_{k}+P_{k}*T_{s}$ where $k$ denotes the discrete time and $T_{s}$ is the sampling time. This describes both the charge and discharge dynamics in this work. The charging power is a function of the energy in the battery as presented in [3] and the discharge power is a function of the incremental velocity, $\dot{\delta}_{k}$. Thus the switched dynamics can defined on $[\underline{E},\bar{E}]$ where Waiting = $\beta(1)$, Charging = $\beta(2)$, and Driving = $\beta(3)$:

[TABLE]

The second continuous state, the discrete distance ($d_{k}\in[0,d_{max}$), is defined as a trip distance counter. $d_{max}$ is the maximum feasible distance that a car would reasonably travel on the given network. This is later useful for forming the convex hull of the continuous-valued variables. The discrete local distance $\epsilon_{k}$ has the same evolution as the trip distance when in $[0,e^{ij})$. Again, the dynamics here in continuous time are $\frac{\partial d}{\partial t}:=\dot{\delta}(t)$. In discrete-time, under the one-step Euler it becomes $d_{k+1}=d_{k}+\dot{\delta}_{k}T_{s}$. The state $d_{k}$ is primarily useful for keeping track of whether the vehicle is moving or not. Given the incremental velocity, $\dot{\delta}_{k}$, the distance counter is defined:

[TABLE]

Similarly, $\epsilon_{k}^{c}\in[0,{e^{ij}}^{+}]$ but with a reset at each node:

[TABLE]

The guard for each of the three states is defined by boolean algebra. The extended logic is defined in Section IV. For the rest of this work, since the logic is the same for the three continuous states, the state vector will be denoted $x_{k}$.

III-F Costs to customers (EVs)

The costs to the customers will be the time spent charging, the time to drive along a specified path, the degradation of the battery, and the time to wait in a queue at a full charging station.

[TABLE]

The cost of travel time along a specific edge $e^{ij}$ was obtained from [8]. Here $t_{e^{ij}}^{0}$ and $c_{e^{ij}}$ are the free flow travel time (in the absence of traffic) and capacity of link $e^{ij}\in\mathcal{A}$ respectively, which are both predetermined. Then, the travel time is a strictly increasing function of the traffic flow $v_{e^{ij}}^{k}$:

[TABLE]

The traffic flow $v_{e^{ij}}^{k}$ along a particular edge can be calculated by summing over our vector-valued exogenous input variable $\vec{\xi^{c}_{k}}$ (explained in Section IV-C) over all cars at any given time $k$, which gives us the total number of cars currently on that highway.

[TABLE]

It is easiest to instead assign a quadratic cost to $C^{t_{driving}}$ that acts directly on the mode selector as a time-series vector instead of penalizing each edge. The time dependence on traffic/congestion can be separated out to be penalized independently. This will be a squared two-norm penalty on (13) i.e. $\|\vec{v_{e^{ij}}^{k}}\|_{2}^{2}$ for the prediction horizon.

Both the charging time costs as well as the battery degradation model were pulled from [3], where $E_{i}$ and $E_{f}$ are the battery energies before and after the charging step.

[TABLE]

The charging time can be directly penalized with a quadratic cost since there is direct access to the time spent charging, $\beta(1)$. The function in (14) is useful for understanding how the charging time increases with increased final state of charge, SOC.

The charging power can be derived from the charging power function in [3]:

[TABLE]

Where m and n are constant parameters for all the charging stations, since we assume they all use identical batteries.

As seen in Fig. 3 above, the rate of charging power decreases and charging duration increases as the SOC of the battery increases beyond the threshold value $E_{c}$. Physically, this implies that it may often be optimal to not charge cars to 100% of their capacity since this would greatly slow down the charging time. Instead, the controller would have them move on to the next available station/node and let other vehicles waiting in queue charge at that station.

Once the mean charging time $\mathbb{E}(t_{charging})$ and mean charging power $\mathbb{E}(P_{ow})$ are determined using the above relations, the expected battery degradation cost can then be calculated for each EV can then be calculated as follows:

[TABLE]

Where parameters a, b, and c are related to the Li-ion battery, which is the most common battery type used in EV applications.

Most of the previous literature encodes the time spent by EVs waiting at a charging outlet, as some function of the arrival rate $\lambda$ at that node and its blocking probability $\nu$, to characterize the random process.

[TABLE]

However, we decided to simplify matters by just penalizing the total waiting time with a quadratic cost. This was partly because neither of us had much prior experience working with queuing theory or Markov chains and also since the customer’s actual waiting time can be obtained directly from our two finite state machines described below in Section IV, if necessary. This does away with the need to estimate wait-times probabilistically.

IV Hybrid systems theoretic formulation

The EV scheduling problem is demonstrated to be piecewise continuous in the two states, $E_{k}andd_{k}$. The problem can be posed as a hybrid systems theoretic MIQP where the equality constraints are linear, the inequality constraints are affine, and the cost quadratic. Furthermore, the position of each vehicle on the highway graph can be represented as a finite-state machine in addition to the dynamics FSM. Thus, there exists two finite-state machines in this formulation.

IV-A Framework

In accordance with the work presented in [21], it is convenient to define a general framework to describe mixed discrete-time horizon automaton trajectories. The superscripts $ct$ denote continuous valued variables while $l$ denotes logical/binary variables. These superscripts are dropped outside of this section in order to simplify notation and allow for indexing cars, times, and nodes. The switched affine system, $x_{k+1}^{ct}=A^{i_{k}}x^{ct}_{k}+B^{i_{k}}\gamma^{ct}_{k}$ is defined in Section III-E. This generally encodes piecewise-affine dynamics. The event generator, $h$, creates a binary output vector, $\delta^{e}$, that is a function of state events that satisfy some affine condition. This is considered to be the endogenous input to the other binary functions as it is an input that is decided based on the evolution of the dynamics in relation to the domains. In the framework developed in EECS C291E/ME290S, it is equivalent to a guard. Next, define the FSM, $f^{l}$, that encodes state transitions and that outputs $x^{l}$. Lastly, the mode selector, $\mu$, is defined to decide which discrete state, $x_{k+1}^{l}$, the finite-state machine is switching to. The general input-output relations for this work’s modeling purposes are shown as follows and are not exhaustive of the work presented in [21]:

[TABLE]

To summarize these equations, the event generator is a function of the continuous valued state, $x^{ct}_{k}$, and input, $\gamma^{ct}_{k}$, where this exogenous input is zero for this work. This can also be time-variant for cases where there exist time-dependent guards. The logic FSM encodes the previously outputted value, the discrete input,$\gamma^{l}_{k}$, and the output from the event generator, $\delta^{e}_{k}$. The mode selector takes in the information from the inputs and decides what mode of the switched affine system to be in. The dynamics then propagate forward one time-step.

The work can be similarly presented in the framework from EECS C291E/ME C290S; however, for formalizing the work for coding purposes with a controller, (18) seemed to be a canonical formulation worth adopting. The analysis (Section V) of the resulting automaton is conducted with the tools developed in EECS C291E.

IV-B Event Generator and Guards

The event generator for this automaton is fairly simple. The domains of the continuous valued variables are box constraints, $x^{ct}\in[\underline{x},\bar{x}]$ except for $\epsilon_{k}\in[0,{e^{ij}}^{+}]$, where the $+$ denotes that depending on the step value, the length of $e^{ij}$ can be overstepped in one time-step. If the domain of $x^{ct}$ is violated, then the problem is infeasible. The first element of the event generator returns $x^{l}=true$ when the upper bound of $\epsilon$ is reached. The first element of the event generator returns $false$ for $\epsilon_{k+1}\in[0,e^{ij})$.

[TABLE]

The second element of the output of the event generator indicates whether the vehicle is moving on an edge, or has arrived at a node:

[TABLE]

IV-C Logic FSM

The logic FSM (boolean output) keeps track of what edge a car is currently on. In particular it encodes the event generator as well as the exogenous input variable, $\gamma_{k}^{l}$. This variable decides whether the vehicle is at a node or on an edge. It is bounded under certain conditions:

[TABLE]

This indicates that the exogenous input is $true$ when the event generator is also $true$ and can be $\{0,1\}$ when the event generator returns $false$ since this physically means that the car can stay or leave a node when it wants to. The input is further restricted when it is on an edge, which is the third case; it says that it has to be zero when it is traversing an edge.

The value of the second input variable, $y_{k}^{l}$, indicates whether a car is charging or not:

[TABLE]

If the car is at a node, indicated by $\gamma_{k+1}^{l}$, then the car may charge. If it is not at a node, then it cannot charge. An additional charging station capacity constraint, denoted $cap^{n}_{k}$ ($\in\mathbb{Z}$) is included as stations can only provide service to a certain number of EVs at a time. For continuity of charging, the condition $y_{k}=1$ is included to allow for cars to keep charging at the next time step if they already were in the previous time step, or $y_{k}=0$ prevents cars from charging if they hadn’t been previously and a station is at capacity. This is where the time indexing becomes muddled. The logic FSM is generally defined to output one-step ahead. However, since the first input affects the second input in a disjunctive condition, it should presumably still be labeled as $k+1$. To check this, the event generator at $k$ could be introduced directly into the condition for the second input. Thus, the time index for this is correct. The time index for the capacity condition is also non-obvious. The way it is written indicates that if the station is currently at capacity, then new cars can’t start charging until one time step has passed during which the station is not at capacity. While this is not optimal for large time steps, it creates an appropriate condition. (27) makes sure that numerous EVs don’t start charging in the next time step and exceed the station’s capacity constraints.

The third, vector-valued exogenous input variable is $\xi_{k}^{l}$. This encodes which edge a vehicle is currently on. Specifically, $\xi_{k}^{l}\in\{0,1\}^{N^{2}}$, or can be represented in a $NxN$ matrix corresponding to whether an edge is selected, since there are a maximum of $N^{2}$ edges in a graph with $N$ nodes. For example, in the vector form, if $\xi_{k}^{l}(5)$ is selected in a five-node graph, then it indicates the vehicle is traveling from node 1 to node 5. This allows for the next edge to be selected each time a car reaches a node. Under the FSM framework with $h$ acting as the index in $[1,N^{2}]$:

[TABLE]

The notation for the latter two cases is poor but essentially indicate that $\xi_{k+1}^{l}$ can only be true if that edge is connected to the current node. Here the constraints on the exogenous inputs under specific constraints were shown but the unconstrained conditions (as it is not deterministic) allow for a solver to choose values that minimize the objective function. The output of the FSM is then a binary vector $x_{k+1}^{l}:=[\gamma_{k+1}^{l},y_{k}^{l},(\xi_{k}^{l})^{T}]^{T}$.

IV-D Mode Selector

The discrete mode is then selected as a function of the results from the event generator, the exogenous input, and the logic FSM. Let Waiting = $\beta(1)$, Charging = $\beta(2)$, and Driving = $\beta(3)$ be binary variables in this illustration of the mode selection:

[TABLE]

The utility of $\xi_{k}^{l}$ becomes apparent now as the evolution of the current edge must be considered. If a node $j\in\mathcal{N}_{i}$, then $\xi_{k}^{l}(h)$ can be positive. The algebra for this relation is inconvenient and non-intuitive (Section XI) but the relation is as follows:

[TABLE]

The incremental velocity mentioned in Section III-E is defined as a non-constant increment as it decreases with congestion on an edge and is highest in free-flow conditions. In general it is of the form $\dot{\delta}_{k}:=\frac{e^{ij}}{t_{ij}}$ where $t_{ij}$ is the time to traverse an edge. For future work this can be exploited but for now it is left as constant.

Thus, the mode for the switched dynamics, as well as the next edge can be chosen. This is also illustrated in the left-hand FSM of Fig. 2.

V Analysis from a Hybrid Systems Perspective

The EV and highway network models described in the previous sections incorporates several elements of hybrid systems, according to the framework developed in the ME 290S/EE291E class. These can be thought of as both 1) continuous, embedded systems (in distance travelled and battery energy level) controlled by discrete logic and 2) multi-agent subsystems, with many subsystems (EVs) interacting with one another.

We can explicitly define the following entities for the hybrid automaton $\mathcal{H}(Q,X,Init,f,Dom,R)$ of the EV model described above:

• •

Discrete states of the EV $q\in Q$ where $Q=\{waiting\;\beta(1),charging\;\beta(2),driving\;\beta(3)\}$ is a finite collection of states possible for the two FSMs shown in Fig. 2.

• •

Continuous states $X=\{d_{k}^{c},E_{k}^{c},\epsilon_{k}^{c}\}\subseteq\mathbb{R}^{3}_{+++}$ for the distance traveled from the source $d_{k}^{c}$ and energy level of the battery $E_{k}^{c}$ or equivalently state of charge $SOC=\frac{E_{k}^{c}}{E^{c}_{max}}$, and $\epsilon_{k}^{c}$ is the counter for distance traveled along a specific edge.

• •

Discrete inputs $\Sigma=\{\gamma_{k}^{c},y_{k}^{c},\vec{\xi}_{k}^{c}$ }, where $\gamma_{k}^{c}\in\{0,1\}$ switches the EV between nodes and edges, $y_{k}^{c}\in\{0,1\}$ switches it between charging and waiting, and $\xi_{k}^{c}\in\{0,1\}^{N^{2}}$ chooses between edges at a node.

• •

Continuous inputs $V=\o$

• •

Initial state $Init\subset Q\;x\;X$ where the car always starts off by waiting at a source node $v_{0,i}$ at a certain distance $d_{k}^{c}$ from the origin, with an initial energy level $E_{k}^{c}\geq 0$ and $\epsilon_{k}^{c}=0$

• •

Continuous dynamics $f$: Both the battery energy level and distance traveled vary continuously and can be represented using continuous state, continuous time differential equations as a vector field. However, one step Euler discretizations were used for this digital control system instead, using the charging/discharging powers and the incremental velocity $\dot{\delta}_{k}$ respectively.

• •

The event generator $\delta^{e}_{k}$ is the guard that switches the EV between edges and nodes in the first FSM, and thus keeps the EV automaton inside the domain (e.g. this ensures that the EV never violates the prescribed domain by running out of charge while traveling along an edge).

• •

A transition relation or reset map $R:Q\;X\;\Sigma\mapsto 2^{Q}$ that describes the switching logic of the finite automaton. In our automaton, this is the mode selector $\mu$ that switches between the discrete states of charging, waiting, and driving states in the second FSM.

• •

Domain for the endogenous inputs,

$Dom=(\beta(1),\{d_{k}^{c}\geq 0,E_{k}^{c}\geq 0,\epsilon_{k}^{c}=0\})\\ \cup\;(\beta(2),\{d_{k}^{c}\geq 0,E_{k}^{c}\geq 0,\epsilon_{k}^{c}=0\})\\ \cup\;(\beta(3),\{d_{k}^{c}\geq 0,E_{k}^{c}>0,\epsilon_{k}^{c}\geq 0\})$

where $\cup$ denotes the union of sets

We attempted to analyze this system from a reachability standpoint but it is computationally challenging and time intensive to explicitly calculate the feasible trajectories/paths as well as all their possible permutations given a particular starting node and initial energy level. However, we can clearly see that all nodes and states in our EV network’s domain are reachable with finite executions, even without explicitly computing the reachable sets. This is because we can assume that any realistic highway road network will most likely be a connected graph i.e. there exists at least one path that you can take between every pair of vertices. Thus given a particular starting node, an EV can reach any other node in finite time and without running out of charge.

Furthermore, we can also argue that our automaton is domain preserving (i.e. the domain of possible discrete states (modes) and continuous states (state of charge and distance traveled $\in\mathbb{R_{+}}$) is an invariant set). Since our initial set of states is in the domain, and all the reachable states in the network from these initial states is in the domain, the hybrid automaton is domain preserving i.e

$Init\subseteq Dom\implies Reach\subseteq Q\;X\;Dom.$

Similar to what we learned in the class, we can then proceed to also define hybrid time sets $\tau$ and trajectories/executions $(\tau,q,x)$ over which our EV system (automaton) is allowed to evolve. For the EV system, we can envision hybrid trajectories with either finite, finite-open or Zeno hybrid time sets. Since $R(q,x)=\mu_{k}^{c}\neq\phi\;\forall$ states $(q,x)\in Dom$, we can conclude that $\mathcal{H}$ is non-blocking. This means that for all reachable states for which continuous evolution (in distance or energy) is impossible, a discrete transition (i.e. switching to a different mode) is possible. Of course, this is only true when the edges are short enough such that a car can traverse an edge without running out of energy and the horizon is long enough such that cars can wait to charge before driving further.

The automaton is also non-deterministic since 1) continuous evolution in $d_{k}^{c}$ and $E_{k}^{c}$ is possible even when discrete transitions between the modes $\{\beta(1),\beta(2),\beta(3)\}$ occur and 2) each discrete transition between any two of these modes could lead to multiple different destinations, implying that $|R(q,x)|\nleq 1$ in general. This makes intuitive sense since given a particular starting node (initial condition), the EV could reach several different destinations and thus there are a whole family of state trajectories (solutions) possible. There could also be multiple ways to travel between the same start and end nodes while meeting all the constraints and minimizing our objective (especially for larger networks than our 5-node example).

Thus, the non-blocking property implies the existence of finite and infinite executions for all initial states i.e. local existence of solutions is guaranteed. Non-determinism implies that even though all executions $(\tau,q,x)$ can be extended to infinite executions, these are not necessarily unique.

By construction, the EV system defined above is also non-Zeno since our model only allows one discrete transition per time step. Thus, $\mathcal{H}$ can never accept Zeno executions since it’s impossible to complete infinite discrete transitions within finite time.

VI Optimization Problem

The modeling framework proposed in the previous section is for each car. In order to account for each car, the conditions hold for each vehicle on the network. Interactions between the cars are only considered insofar as they affect the incremental velocity, $\dot{\delta}_{k}$, and the availability of a charger at a station.

The idea behind this formulation is that it will be run as a receding horizon control (RHC) problem. The prediction horizon, Hp, will be computed, applied, and then calculated again at the next sampling time in the control horizon $H_{p}$.

[TABLE]

The constraints on the problem are the mixed logic dynamics defined in Section III-E and IV. The box constraints here are reiterated for clarity, and the initial and terminal constraints are added. The counter variable, $\epsilon_{k}$, is zero at the start and the end. The starting and ending nodes are dictated by $\xi_{k}$; the actual indexing on $\xi$ is non-trivial as it relates the coming-from and going-to a node in index arrays. It is included in the Appendix (Section XI) for brevity here but are fully-defined in how a node is selected.

VI-A Disciplined formulation

In order to convert this MIQP into a problem for a numerical solver, it must be reformulated into a problem that has ”good” numerical properties and uses the tools available in disciplined convex programming.

The canonical way of expressing piecewise functions is via the big-M relaxation [21], [20]. The problem,

[TABLE]

can be formulated to be,

[TABLE]

The four inequality constraints represent the two equality constraints where the first statement of (40) is tight when the logical, $\delta_{k}$ is $false$ in the first two constraints and the second statement of (40) is tight when $\delta_{k}$ is $true$. The point of interest is when the bound is not tight. In this case the values of $m_{i},M_{i}$ are important as they define the box constraints on the continuous variable $x_{k}\in\mathbb{R}$. If the bounds are arbitrarily small, (for $m_{i}$ it means $|m_{i}|$ is large). Then the domain for $x_{k}$ is very large. This leads to bad numerics for solvers. For bounds that are on the order of the actual solutions bounds, the solve time can be longer than necessary due to extra branching [20] in branch-and-bound solvers.

Thus, ”good” means that the convex hull of the piecewise-polytopes is found. This is generally non-trivial as the feasible set at a specific timestep is a subset of the box constraints that define the domain of the states and the inputs for all time. Fortunately, the modeling languages YALMIP provides a tool that allows for direct implications to be made. YALMIP uses the function implies(a,b) for logic programming to imply that if a is $true$, then b is too. If a is $false$, then b may or may not be $false$ as well.

While implies() is convenient in clarity of expressions, it does have existing bugs. For instance, through talking with the developer of YALMIP, it was discovered that the operator not() was not overloaded while all the other MATLAB logical expressions were. Errors such as this while using a somewhat black-box tool (based on big-M but with many numerical tools adjusting it) made debugging the program more difficult.

The problem was coded in both YALMIP using implies() and purely with the big-M formulation. While including more binary variables in both cases is cumbersome and intuitively more complicated for a solver, it actually creates cleaner big-M relaxations in both cases. The reason for this is that it creates disjunctive behavior between conditional cases involving multiple continuous-valued variables. For example, the trivial implication ”if $x>5$, then $y<3$” where y is the variable that should be constrained can result in the value of x being modified so that the condition is true. A disjunctive equivalent would then include $\delta\in\{0,1\}$ so that ”if $x>5$, then $\delta$. If $\delta$, then $y<3$”.

VII Process and Preliminary Results

VII-A Process

Various reformulations of the problem have been undertaken over the past two months. Originally, the problem used a model that was more consistent with previous work done [8],[3]. In order to keep it in the same space as more hybrid system theoretic approaches, a first draft of the current model was proposed.

The first draft incorporated the basic dynamics of keeping track of trip distance and energy in the battery. In order to reduce the order of computational complexity, it was thought that it would be best to pre-compute the $n$ shortest total travel time trajectories, for example $n=3$, under fully available charging stations and zero traffic congestion (i.e. free-flow). While this would seem to reduce the solution space when using a branch-and-bound solver as the number of ”branching” variables is reduced, the scaling of this was very poor because every combination of each car’s $n$ feasible trajectories would need to be computed with every other car’s $n$ feasible trajectories.

The next step was to allow for the model to switch edges itself. While this would increase the complexity of each problem being solved, it removed all of the pre-computation. Since moving to this scheme, the basic model has been modified multiple times. The main modifications have been the guard and reset relations. The guards for the problem were originally defined based on the relation of the counter variable $\epsilon$ to the current edge length in both transitions to and from a node. This often led to blocking cases where once in the driving state, the vehicle could not stop again to charge.

Throughout this process, it was difficult to understand how relaxing a constraint would affect the feasible set. Generally, visualizing N-step invariant sets is straightforward with time-invariant dynamics that are continuous-valued functions. With mixed-integer programs it becomes more difficult. YALMIP has convenient plotting functions: the constraints can be formed into n-dimensional polyhedra and projected down onto arbitrary lower dimensions. Unfortunately, the computation time for using this for mixed-integer constraints is prohibitive due to the large number of binary variables.

One of the key areas that is misunderstood in this formulation is the evolution of the logical operators as defined in Section IV-C. Somewhat counter-intuitively, there are numerous cases where ”if $z_{k+1}$, then $w_{k+1}$”, with $z_{k+1},w_{k+1}\in\{0,1\}$, is different (i.e. feasible/infeasible) than ”if $z_{k}$, then $w_{k}$” where the entire disjunctive logic’s time indices is changed appropriately.

VII-B Preliminary Results

VII-B1 Improper constraints

The current version of the code exhibits somewhat reasonable behavior when constrained under the switched dynamics and logic defined in Section III but there exists faulty edge conditions due to unknown errors in the discrete time logic.

In particular, Fig. 4 shows the trajectory of a vehicle under all of the constraints save for the station capacity. The objective here is to simply maximize the distance traveled via a quadratic penalty on the difference between the terminal distance and some arbitrarily large distance. This is the test cost function to exploit the modes and the dynamics. In the example, it can be seen from Fig. 6 that the first car does not charge despite reaching a switching condition (a node) where it can charge under the relaxed constraints. Instead, it waits at a node. This reduces the cumulative distance traveled and thus exposes an issue in the logic defined above. The second car, however, is able to charge, and does so for one time-step. Ostensibly, neither should be waiting in the beginning but do so. This is thought to be the manifestation of issues with initialization and time indexing errors at this point in the project.

The local distance counter, $\epsilon$, is shown in Fig. 5 to be increasing over $H_{p}=10$ and resets when a node is reached. The counter logic is thus illustrated to function well. The fourth quadrant indicates that the current edge changes. This is indicated in terms of the vector-value of the indicator function $\xi$. The edges take values defined by the set of edges, $\mathcal{E}$, then and can change at a node.

The terminal condition is numerically poor in this example where the end-node constraint cannot be written explicitly due to the construction focused on edges and adjacency instead of nodes. This results in the condition: At Hp, the car must be within $\Delta$ of the final node, j, (either pointing to j (i.e. ij) or from j (i.e. ji)). The formulation is non-intuitive making it difficult to debug. In this example the cars should point from the end destination node 4. This means that it can take the value (4,9,14,19,24) insofar as that edge is connected by a node that is incident with the previous edge. The resultant edge is in fact labeled 20, which corresponds to node 5 pointing at node 4. This outcome occurs because the terminal constraint on $\epsilon$ was relaxed to maintain feasibility. The terminal constraint here needs to be reformulated to be more intuitive.

VII-B2 Computational complexity- scaling

Another issue is that the problem scales poorly. The parameters that change with the problem parameters are $H_{p}$, $C$, and $N$. Fig. 7 indicates that the time to solve the problem increases dramatically with increased prediction horizon, Hp, over $H_{p}=10$. Increasing the number of EVs has diminished effects at some threshold and is thus not considered to be as much of a computational barrier. The number of nodes in the network does not seem to have a drastic effect on computation but larger networks have not been tested.

The computational issue lies in the complexity underlying the prediction horizon. In particular, there is a trade-off between the time-step size and the horizon length. With larger time-steps, the precision of the model begins to falter and cars can overshoot nodes (i.e. $\epsilon$ can attain a value much greater than its current edge length in one time-step). This results in inaccurate characterization of the trip distance of a vehicle. Smaller time-steps improve this. However, in order to reach a terminal node from a defined initial node, the time horizon must be sufficiently large. The requirements for these two parameters are inversely proportional. As a result of this complexity and desired small time-step, the current algorithm (i.e. the exact solution provided by Gurobi–a branch-and-bound solver) is not a good option.

VIII Next steps

Moving forward, we aim to polish and fix MATLAB code bugs in our most recent approach for programming the logic FSM and counter variables. In particular, the definitions for the initial and terminal conditions need to be refined or reformulated. This will allow the optimal trajectory to be computed correctly given starting and ending nodes on the small example network.

However, we recognize that our current solution does not scale well with larger networks as we increase the time horizon and number of participating agents (stations and EVs) greatly. Considering the relatively high-level nature of this scheduling formulation, it may make sense to explore heuristic approaches that scale very well and compromise the optimality of the solution. Further time should be spent visualizing the feasible sets and understand how they vary with time. Understanding the nuances of the variant sets will enable the logic to be smoothed out and the time indices to be coded with more rigorous reasoning.

The last step in mind is to simulate this using data from real highway networks that have high EV penetration, multiple charging stations that have finite capacity, and many, perhaps directed, edges. We will try to use data sets that have already been used by [8] and [3] in the past.

IX Conclusion

This work demonstrated the formulation of an EV scheduling model that is hybrid systems theoretic. The model is represented using two frameworks, one that is more canonical from a control standpoint, and one that is more rigorous when considering desirable properties of the model. The simulations indicate that the model has logical, or at least indexing, errors. The simulations illustrate the working principle of a switched dynamical system with extended decision-making logic.

Cost functions accounting for the aggregate cost to all players were introduced, many drawn from recent works in the field and cited as such. The utility of creating such a cost function is yet to be realized operationally.

The immediate next steps for this work are laid out so that the model can be implemented in simulations using data from real highway networks.

X Acknowledgements

We would like to thank Prof. Claire Tomlin and Forrest Laine for their informative in-class lectures and discussions related to modeling hybrid systems, as well as for offering advice and help during office hours. We also thank Prof. Francesco Borrelli since we drew a lot of inspiration for our latest hybrid systems theoretic approach from his book on Model Predictive Control. Finally, we would like to thank Johan Lofberg for offering some assistance in working through existing issues with implies() in YALMIP.

XI Appendix

XI-A Indexing with respect to edges and nodes

The constraints are built iteratively. In particular, in the YALMIP modeling language with the ”constraints” list recursively defined: