Constraint-coupled: charging of Plug-in Electric Vehicles (PEVs)¶
For the constraint-coupled set-up, we consider the problem of determining an optimal overnight charging schedule for a fleet of Plug-in Electric Vehicles (PEVs) [FalMa17]. The model described in this page is inspired by the model in the paper [VuEs16] (we consider the “charge-only” case without the integer constraints on the input variables). The complete code of this example is given at the end of this page.
Problem formulation¶
Suppose there is a fleet of \(N\) PEVs (agents) that must be charged by drawing power from the same electricity distribution network. Assuming the vehicles are connected to the grid at a certain time (e.g., at midnight), the goal is to determine an optimal overnight schedule to charge the vehicles, since the electricity price varies during the charging period.
Formally, we divide the entire charging period into a total of \(T = 24\) time slots, each one of duration \(\Delta T = 20\) minutes. For each PEV \(i \in \{1, \ldots, N\}\), the charging power at time step \(k\) is equal to \(P_i u_i(k)\), where \(u_i(k) \in [0, 1]\) is the input to the system and \(P_i\) is the maximum charging power that can be fed to the \(i\)-th vehicle.
System model¶
The state of charge of the \(i\)-th battery is denoted by \(e_i(k)\), its initial state of charge is \(E_i^\text{init}\), which by the end of the charging period has to attain at least \(E_i^\text{ref}\). The charging conversion efficiency is denoted by \(\zeta_i^\text{u} \triangleq 1 - \zeta_i\), where \(\zeta_i > 0\) encodes the conversion losses. The battery’s capacity limits are denoted by \(E_i^\text{min}, E_i^\text{max} \ge 0\). The system’s dynamics are therefore given by
To model congestion avoidance of the power grid, we further consider the following (linear) coupling constraints among all the variables
where \(P^\text{max}\) is the maximum power that the be drawn from the grid.
Optimization problem¶
We assume that, at each time slot \(k\), electricity has unit cost equal to \(C^\text{u}(k)\). Since the goal is to minimize the overall consumed energy price, the global optimization problem can be posed as
The problem is recognized to be a constraint-coupled problem, with local variables \(x_i\) equal to the stack of \(e_i(k), u_i(k)\) for \(k \in \{0, \ldots, T-1\}\), plus \(e_i(T)\). The local objective function is equal to
the local constraint set is equal to
and the local coupling constraint function \(g_i : \mathbb{R}^{2T+1} \rightarrow \mathbb{R}^T\) has components
The goal is to make each agent compute its portion \(x_i^\star = (e_i^\star, u_i^\star)\) of an optimal solution \((x_1^\star, \ldots, x_N^\star)\) of the optimization problem, so that all of them can know their own assignment of the optimal charging schedule, given by \((u_i^\star(0), \ldots, u_i^\star(T-1))\).
Data generation model¶
The data are generated according to table in [VuEs16] (see Appendix).
Complete code¶
###############################################################
# CONSTRAINT-COUPLED Example
# Charging of Plug-in Electric Vehicles (PEVs)
#
# The problem consists of finding an optimal overnight schedule to
# charge electric vehicles. See [Vu16] for the problem model
# and generation parameters.
#
# Note: in this example we consider the "charging-only" case
###############################################################
# Compared Algorithms:
#
# - Distributed Dual Subgradient
# - Distributed Primal Decomposition
###############################################################
import dill as pickle
import numpy as np
from numpy.random import uniform as rnd
from mpi4py import MPI
from disropt.agents import Agent
from disropt.algorithms import DualSubgradientMethod, PrimalDecomposition
from disropt.functions import Variable
from disropt.utils.graph_constructor import binomial_random_graph, metropolis_hastings, binomial_random_graph_sequence
from disropt.problems import ConstraintCoupledProblem
# get MPI info
NN = MPI.COMM_WORLD.Get_size()
agent_id = MPI.COMM_WORLD.Get_rank()
# Generate a common graph (everyone uses the same seed)
Adj = binomial_random_graph(NN, p=0.2, seed=1)
W = metropolis_hastings(Adj)
# generate edge probabilities
edge_prob = np.random.uniform(0.3, 0.9, (NN, NN))
edge_prob[Adj == 0] = 0
i_lower = np.tril_indices(NN)
edge_prob[i_lower] = edge_prob.T[i_lower] # symmetrize
#####################
# Generate problem parameters
#####################
# Problem parameters are generated according to the table in [Vu16]
#### Common parameters
TT = 24 # number of time windows
DeltaT = 20 # minutes
# PP_max = 3 * NN # kWh
PP_max = 0.5 * NN # kWh
CC_u = rnd(19,35, (TT, 1)) # EUR/MWh - TT entries
#### Individual parameters
np.random.seed(10*agent_id)
PP = rnd(3,5) # kW
EE_min = 1 # kWh
EE_max = rnd(8,16) # kWh
EE_init = rnd(0.2,0.5) * EE_max # kWh
EE_ref = rnd(0.55,0.8) * EE_max # kWh
zeta_u = 1 - rnd(0.015, 0.075) # pure number
#####################
# Generate problem object
#####################
# normalize unit measures
DeltaT = DeltaT/60 # minutes -> hours
CC_u = CC_u/1e3 # Euro/MWh -> Euro/KWh
# optimization variables
z = Variable(2*TT + 1) # stack of e (state of charge) and u (input charging power)
e = np.vstack((np.eye(TT+1), np.zeros((TT, TT+1)))) @ z # T+1 components (from 0 to T)
u = np.vstack((np.zeros((TT+1, TT)), np.eye(TT))) @ z # T components (from 0 to T-1)
# objective function
obj_func = PP * (CC_u @ u)
# coupling function
coupling_func = PP*u - (PP_max/NN)
# local constraints
e_0 = np.zeros((TT+1, 1))
e_T = np.zeros((TT+1, 1))
e_0[0] = 1
e_T[TT] = 1
constr = [e_0 @ e == EE_init, e_T @ e >= EE_ref] # feedback and reference constraints
for kk in range(0, TT):
e_cur = np.zeros((TT+1, 1))
u_cur = np.zeros((TT, 1))
e_new = np.zeros((TT+1, 1))
e_cur[kk] = 1
u_cur[kk] = 1
e_new[kk+1] = 1
constr.append(e_new @ e == e_cur @ e + PP*DeltaT*zeta_u*u_cur @ u) # dynamics
constr.extend([u_cur @ u <= 1, u_cur @ u >= 0]) # input constraints
constr.extend([e_new @ e <= EE_max, e_new @ e >= EE_min]) # state constraints
#####################
# Distributed algorithms
#####################
# local agent and problem
agent = Agent(
in_neighbors=np.nonzero(Adj[agent_id, :])[0].tolist(),
out_neighbors=np.nonzero(Adj[:, agent_id])[0].tolist(),
in_weights=W[agent_id, :].tolist())
pb = ConstraintCoupledProblem(obj_func, constr, coupling_func)
agent.set_problem(pb)
# instantiate the algorithms
# y0 = np.zeros((TT, 1))
# mu0 = np.zeros((TT, 1))
y0 = 10*np.random.rand(TT, 1)
mu0 = 10*np.random.rand(TT, 1)
theothers = [i for i in range(NN) if i != agent_id]
y_others = agent.communicator.neighbors_exchange(y0, theothers, theothers, False)
y_others[agent_id] = y0
y_mean = sum([x for _, x in y_others.items()])/NN
y0 -= y_mean
dds = DualSubgradientMethod(agent=agent,
initial_condition=mu0,
enable_log=True)
dpd = PrimalDecomposition (agent=agent,
initial_condition=y0,
enable_log=True)
############################################
num_iterations = 1000
# generate sequence of adjacency matrices
Adj_seq = binomial_random_graph_sequence(Adj, num_iterations, edge_prob, NN, 1)
def step_gen(k): # define a stepsize generator
return 1/((k+1)**0.6)
# def step_gen(k): # define a stepsize generator
# return 0.01
def update_graph(k):
Adj_k = Adj_seq[k]
W_k = metropolis_hastings(Adj_k)
agent.set_neighbors(in_neighbors=np.nonzero(Adj_k[agent_id, :])[0].tolist(),
out_neighbors=np.nonzero(Adj_k[:, agent_id])[0].tolist())
agent.set_weights(in_weights=W_k[agent_id, :].tolist())
# run the algorithms
if agent_id == 0:
print("Distributed dual subgradient")
_, dds_seq = dds.run(iterations=num_iterations, stepsize=step_gen, verbose=True)
if agent_id == 0:
print("Distributed primal decomposition")
dpd_seq, _, _ = dpd.run(iterations=num_iterations, stepsize=step_gen, M=30.0, verbose=True)
# save information
if agent_id == 0:
with open('info.pkl', 'wb') as output:
pickle.dump({'N': NN, 'iterations': num_iterations, 'n_coupling': TT}, output, pickle.HIGHEST_PROTOCOL)
with open('agent_{}_objective_func.pkl'.format(agent_id), 'wb') as output:
pickle.dump(obj_func, output, pickle.HIGHEST_PROTOCOL)
with open('agent_{}_coupling_func.pkl'.format(agent_id), 'wb') as output:
pickle.dump(coupling_func, output, pickle.HIGHEST_PROTOCOL)
with open('agent_{}_local_constr.pkl'.format(agent_id), 'wb') as output:
pickle.dump(constr, output, pickle.HIGHEST_PROTOCOL)
np.save("agent_{}_seq_dds.npy".format(agent_id), dds_seq)
np.save("agent_{}_seq_dpd.npy".format(agent_id), dpd_seq)
import numpy as np
import matplotlib.pyplot as plt
from disropt.problems import Problem
from disropt.functions import ExtendedFunction
from disropt.constraints import ExtendedConstraint
import dill as pickle
# initialize
with open('info.pkl', 'rb') as inp:
info = pickle.load(inp)
NN = info['N']
iters = info['iterations']
SS = info['n_coupling']
# load agent data
seq_dds = {}
seq_dpd = {}
local_obj_func = {}
local_coup_func = {}
local_constr = {}
for i in range(NN):
seq_dds[i] = np.load("agent_{}_seq_dds.npy".format(i))
seq_dpd[i] = np.load("agent_{}_seq_dpd.npy".format(i))
with open('agent_{}_objective_func.pkl'.format(i), 'rb') as inp:
local_obj_func[i] = pickle.load(inp)
with open('agent_{}_coupling_func.pkl'.format(i), 'rb') as inp:
local_coup_func[i] = pickle.load(inp)
with open('agent_{}_local_constr.pkl'.format(i), 'rb') as inp:
local_constr[i] = pickle.load(inp)
# solve centralized problem
dim = sum([f.input_shape[0] for _, f in local_obj_func.items()]) # size of overall variable
centr_obj_func = 0
centr_constr = []
centr_coupling_func = 0
pos = 0
for i in range(NN):
n_i = local_obj_func[i].input_shape[0]
centr_obj_func += ExtendedFunction(local_obj_func[i], n_var = dim-n_i, pos=pos)
centr_constr.extend(ExtendedConstraint(local_constr[i], n_var = dim-n_i, pos=pos))
centr_coupling_func += ExtendedFunction(local_coup_func[i], n_var = dim-n_i, pos=pos)
pos += n_i
centr_constr.append(centr_coupling_func <= 0)
global_pb = Problem(centr_obj_func, centr_constr)
x_centr = global_pb.solve()
cost_centr = centr_obj_func.eval(x_centr).flatten()
x_centr = x_centr.flatten()
# compute costs and coupling values
cost_seq_dds = np.zeros(iters) - cost_centr
cost_seq_dpd = np.zeros(iters) - cost_centr
coup_seq_dds = np.zeros((iters, SS))
coup_seq_dpd = np.zeros((iters, SS))
for t in range(iters):
for i in range(NN):
cost_seq_dds[t] += local_obj_func[i].eval(seq_dds[i][t, :])
cost_seq_dpd[t] += local_obj_func[i].eval(seq_dpd[i][t, :])
coup_seq_dds[t] += np.squeeze(local_coup_func[i].eval(seq_dds[i][t, :]))
coup_seq_dpd[t] += np.squeeze(local_coup_func[i].eval(seq_dpd[i][t, :]))
# plot cost
plt.figure()
plt.title('Evolution of cost error')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$| (\sum_{i=1}^N f_i(x_i^k) - f^\star)/f^\star |$")
plt.semilogy(np.arange(iters), np.abs(cost_seq_dds/cost_centr), label='Distributed Dual Subgradient')
plt.semilogy(np.arange(iters), np.abs(cost_seq_dpd/cost_centr), label='Distributed Primal Decomposition')
plt.legend()
# plot maximum coupling constraint value
plt.figure()
plt.title('Evolution of maximum coupling constraint value')
plt.xlabel(r"iteration $k$")
plt.ylabel(r"$\max_{s} \: \sum_{i=1}^N g_{is}(x_i^k)$")
plt.plot(np.arange(iters), np.amax(coup_seq_dds, axis=1), label='Distributed Dual Subgradient')
plt.plot(np.arange(iters), np.amax(coup_seq_dpd, axis=1), label='Distributed Primal Decomposition')
plt.legend()
plt.show()
The two files can be executed by issuing the following commands in the example folder:
> mpirun -np 50 --oversubscribe python launcher.py
> python results.py
References
- VuEs16(1,2)
Vujanic, R., Esfahani, P. M., Goulart, P. J., Mariéthoz, S., & Morari, M. (2016). A decomposition method for large scale MILPs, with performance guarantees and a power system application. Automatica, 67, 144-156.
- FalMa17
Falsone, A., Margellos, K., Garatti, S., & Prandini, M. (2017). Dual decomposition for multi-agent distributed optimization with coupling constraints. Automatica, 84, 149-158.