Optimising EV Charging Station Selection using PuLP in Python
I am trying to work out an optimization problem using PuLP in Python.
OBJECTIVE: Find the set with the fewest # of charging stations, such that the vehicle is able to complete the route without running out of charge.
Further details below.
* SET UP *
e_c is the energy consumed in each trip segment
e_r is the energy recharged at the charging station at the end of that trip segment
So the vehicle starts at 100% state of charge, and consumes e_c[0] energy to reach the first charging node, where it recharges e_r[0] kWh. It then consumes e_c[1] energy to reach the second charging node, where it recharges e_r[1] kWh.
BC = 100 # battery capacity
SOC_in = 1.0 # (%) initial state of charge
SOC_max = 1.0 #(%) maximum SOC
# energy consumed in each trip segment
e_c = [23.5, 4.6, 4.4, 33.8, 89.5, 0. , 24.5, 4.7, 21.5, 0. , 0. ,
16.2, 4.8, 20.1, 0. , 62.8, 0.6, 40.8, 41.1, 5.1, 16.7, 0. ,
39.6, 9.6, 9.8, 23.7, 0. , 0. , 23.4, 5.3, 8.6, 23.7, 18.6,
5.8, 9.1, 0.1, 0. , 0.1, 0. , 0.3, 19.8, 0. , 4.9, 23.4,
7.7, 9.6, 8. , 1. , 0. , 0. , 0. , 0. , 28.1, 5.5, 9.3,
0. , 0.1, 27.9, 41.3, 2.8, 20.8, 6.4, 5.7, 28.4, 7.7, 35.2,
2.8, 21. , 11.8, 0.1, 21.6, 7.4, 13.5, 8.6, 7.5, 0. , 12.8,
10.5, 0.4, 0. , 0. , 0.5, 37.2, 5.8, 19. , 10.7, 21.7, 17.1,
0.9, 0. , 17.1, 12. , 7.9, 2.4, 0.2, 20.6, 4.9, 21.1]
# energy recharged in each trip segment
e_r = [30.7, 64.8, 31.2, 181.2, 2.2, 66.5, 30.3, 52.8, 43.8,
54.1, 151.1, 33.9, 22.2, 68.3, 66.4, 70.2, 170.5, 260.2,
32.7, 72.6, 89.1, 138.1, 54.6, 70.2, 24.3, 42.2, 44.4,
38. , 38.8, 30.1, 28.7, 179.4, 48.6, 96.6, 0.2, 38.5,
44.5, 0.4, 0.4, 28.8, 13.5, 25.5, 44.9, 196.8, 269.8,
215.7, 81.1, 166.6, 0.4, 0.4, 0.1, 22.3, 111.9, 36.9,
43.3, 38. , 712.2, 9.7, 25.2, 242.3, 24.5, 49.4, 31.6,
231.6, 0.4, 27.5, 228.6, 31.7, 685.1, 812.4, 84.5, 36.7,
28.4, 214.7, 0.5, 46. , 61.3, 0.4, 0.4, 22.4, 0.4,
0.5, 47.5, 27. , 43.8, 182.3, 24.2, 76.3, 54. , 748.8,
72. , 46. , 31.3, 113.2, 0.4, 22.7, 70.6, 78.9]
candidate_nodes = range(len(e_c))
n = len(candidate_nodes)
* SANITY CHECK *
To check that it is possible for the vehicle to make it all the way through the route without running out of charge, I did a simple simulation of the vehicle going through the route as if all charging stations were available. This code threw no error.
max_energy = 100 #kWh
energy = max_energy
for idx,value in enumerate(e_c):
# subtract the energy needed to complete this trip segment
energy -= value
# check if the vehicle's battery dropped below 0
if energy < 0:
raise ValueError("vehicle ran out of charge!")
# add the energy from charging at the node, but no more than max battery capacity
energy = np.minimum(energy + e_r[idx], max_energy)
* Optimization Model *
In plain English, my optimization problem is:
- Minimise the number of charging stations set up, such that the vehicle is able to complete the route without running out of charge.
Other constraints:
The vehicle starts at 100% battery capacity.
The vehicle can charge to no more than 100% battery capacity.
Decision variables:
X[i] = {0,1} - Binary variable indicating whether a station is built at node i
E[i] - Continuous variable indicating the amount of energy vehicle has at node i
Objective:
min sum(X)
* PuLP Model Set-Up *
Here is the code for the PuLP model set-up
from PuLP import *
# Define optimization problem
model = LpProblem("EV charging optimization modellem", LpMinimize)
### Decision variables
X = LpVariable.dicts("UseLocation", candidate_nodes, cat=LpBinary)
E = LpVariable.dicts("energy", range(len(e)), lowBound=0, upBound = 100, cat = LpContinuous)
### Objective function
model += lpSum(X[i] for i in candidate_nodes)
### Constraints
# Initial energy level
model += E[S_start] == SOC_in * BC
for i in range(n-1):
# Energy at next stop (only add charge if charging station is built there)
model += E[i+1] == E[i] - e_c[i] + (e_r[i] * X[i])
# Energy does not exceed battery capacity
model += E[i+1] <= BC * SOC_max
# Has enough energy to reach the next stop without going below minimum battery level
model += E[i] >= e_c[i]
# Solve
status = model.solve()
The code runs successfully but provides the following output:
command line - /Users/Orie4416/opt/anaconda3/lib/python3.9/site-packages/pulp/solverdir/cbc/osx/64/cbc /var/folders/pb/drvp7crx7qn3x5331z925m1c0000gt/T/04ea923ff7ad4c48a355ebf128973c32-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/pb/drvp7crx7qn3x5331z925m1c0000gt/T/04ea923ff7ad4c48a355ebf128973c32-pulp.sol (default strategy 1)
At line 2 NAME MODEL
At line 3 ROWS
At line 297 COLUMNS
At line 1078 RHS
At line 1371 BOUNDS
At line 1568 ENDATA
Problem MODEL has 292 rows, 196 columns and 486 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 5.47589 - 0.00 seconds
Cgl0000I Cut generators found to be infeasible! (or unbounded)
Cgl0000I Cut generators found to be infeasible! (or unbounded)
Pre-processing says infeasible or unbounded
Option for printingOptions changed from normal to all
Total time (CPU seconds): 0.01 (Wallclock seconds): 0.04
Optimization failed - problem is infeasible
/Users/Orie4416/opt/anaconda3/lib/python3.9/site-packages/pulp/pulp.py:1352: UserWarning: Spaces are not permitted in the name. Converted to '_'
warnings.warn("Spaces are not permitted in the name. Converted to '_'")
Anyone got any ideas?
The output you have says "unbounded or infeasible" so you should be considering both of those problems when troubleshooting...
Is it unbounded? Not likely... You have it set as a minimization and you have lower bounds on both of your variables and your objective is a positive valued sum by inspection.
Is it infeasible? Well, let's take a look at your 2 charging constraints. The first says that the E[i] is whatever it shows up with plus the full charge available at the station, if it is built. The following constraint says the E[i] can't be more than the max (good idea!). So the problem is that you need a less than or equals constraint in your first constraint so it doesn't force a potential over-charge, which can make the pair of two constraints infeasible--as written the car can only take a charge at a station where taking the full charge will result in less than the max capacity. Make sense? So that limitation is causing infeasibility.
One byproduct of this to think about is that the model (as written) doesn't care about carrying any "extra charge" in the vehicle. So if it takes 61 units to get to the next charging station, it is perfectly happy leaving the previous charging station with 61 units, even if a "full tank" where available. E[i] here may assume any value in between 61 and what is available, as they all produce the same objective value.
If that troubles you, you might consider a minuscule bonus in the objective for keeping as full of a tank as possible. How small? Well, what you don't want to happen is have the bonus be large enough to outweigh an additional station, so the total "bonus" possible must be less than 1 as you are counting stations. So, I would maybe use a weighting factor of 1/(100 * 100) which would be the bonus for having full charge at 100 nodes which ensures that the total of all the charge states could never be > 1. (See the code.) I solved it both ways by replacing the objective function and re-solving.
Picture explains it better, I think. As you may note, the selection of which nodes to build is different (the dots), but the number to build is the same (equivalent objective value minus the bonus). I'd be more comfortable riding in the blue-line car. ;)
Then, after fixing a couple typos, your model runs just fine.
Code:
from pulp import *
from matplotlib import pyplot as plt
BC = 100 # battery capacity
SOC_in = 1.0 # (%) initial state of charge
SOC_max = 1.0 #(%) maximum SOC
# energy consumed in each trip segment
e_c = [23.5, 4.6, 4.4, 33.8, 89.5, 0. , 24.5, 4.7, 21.5, 0. , 0. ,
16.2, 4.8, 20.1, 0. , 62.8, 0.6, 40.8, 41.1, 5.1, 16.7, 0. ,
39.6, 9.6, 9.8, 23.7, 0. , 0. , 23.4, 5.3, 8.6, 23.7, 18.6,
5.8, 9.1, 0.1, 0. , 0.1, 0. , 0.3, 19.8, 0. , 4.9, 23.4,
7.7, 9.6, 8. , 1. , 0. , 0. , 0. , 0. , 28.1, 5.5, 9.3,
0. , 0.1, 27.9, 41.3, 2.8, 20.8, 6.4, 5.7, 28.4, 7.7, 35.2,
2.8, 21. , 11.8, 0.1, 21.6, 7.4, 13.5, 8.6, 7.5, 0. , 12.8,
10.5, 0.4, 0. , 0. , 0.5, 37.2, 5.8, 19. , 10.7, 21.7, 17.1,
0.9, 0. , 17.1, 12. , 7.9, 2.4, 0.2, 20.6, 4.9, 21.1]
# energy recharged in each trip segment
e_r = [30.7, 64.8, 31.2, 181.2, 2.2, 66.5, 30.3, 52.8, 43.8,
54.1, 151.1, 33.9, 22.2, 68.3, 66.4, 70.2, 170.5, 260.2,
32.7, 72.6, 89.1, 138.1, 54.6, 70.2, 24.3, 42.2, 44.4,
38. , 38.8, 30.1, 28.7, 179.4, 48.6, 96.6, 0.2, 38.5,
44.5, 0.4, 0.4, 28.8, 13.5, 25.5, 44.9, 196.8, 269.8,
215.7, 81.1, 166.6, 0.4, 0.4, 0.1, 22.3, 111.9, 36.9,
43.3, 38. , 712.2, 9.7, 25.2, 242.3, 24.5, 49.4, 31.6,
231.6, 0.4, 27.5, 228.6, 31.7, 685.1, 812.4, 84.5, 36.7,
28.4, 214.7, 0.5, 46. , 61.3, 0.4, 0.4, 22.4, 0.4,
0.5, 47.5, 27. , 43.8, 182.3, 24.2, 76.3, 54. , 748.8,
72. , 46. , 31.3, 113.2, 0.4, 22.7, 70.6, 78.9]
candidate_nodes = range(len(e_c))
n = len(candidate_nodes)
# Define optimization problem
model = LpProblem("EV charging optimization modellem", LpMinimize)
### Decision variables
X = LpVariable.dicts("UseLocation", candidate_nodes, cat=LpBinary)
E = LpVariable.dicts("energy", candidate_nodes, lowBound=0, upBound = 100, cat = LpContinuous)
### Objective function
model += lpSum(X[i] for i in candidate_nodes)
### Constraints
# Initial energy level
model += E[0] == SOC_in * BC # typo fixed
for i in range(n-1):
# Energy at next stop (only add charge if charging station is built there)
model += E[i+1] <= E[i] - e_c[i] + (e_r[i] * X[i]) # less than or equal => prevent over-charge
# Energy does not exceed battery capacity
model += E[i+1] <= BC * SOC_max
# Has enough energy to reach the next stop without going below minimum battery level
model += E[i] >= e_c[i]
# Solve
status = model.solve()
# print the "build plan":
build_plan = [i for i in candidate_nodes if X[i].varValue > 0]
print('build at nodes: ')
print(build_plan)
plt.plot(candidate_nodes, [E[i].varValue for i in candidate_nodes], color='r', label='no bonus')
plt.scatter(build_plan, [0 for t in build_plan], s=10, alpha=0.6, color='orange', label='no-bonus charge')
# replace the objective function with this
model += lpSum(X[i] for i in candidate_nodes) - 1/(100*100)*lpSum(E[i] for i in candidate_nodes)
status2 = model.solve()
# QA the build plan has the same NUMBER of build points
build_plan_2 = [i for i in candidate_nodes if X[i].varValue > 0]
assert len(build_plan) == len(build_plan_2)
print('build at nodes: ', build_plan_2)
plt.plot(candidate_nodes, [E[i].varValue for i in candidate_nodes], color='b', label='bonus')
plt.scatter(build_plan_2, [0 for t in build_plan_2], s=10, alpha=0.6, color='g', label='charge')
plt.legend(loc='upper right')
plt.show()
