{
"cells": [
{
"cell_type": "markdown",
"id": "261dd435-3512-4a17-ae65-8c852d5ace69",
"metadata": {},
"source": [
"```{index} disjunctive programming\n",
"```\n",
"```{index} single: application; electric vehicles\n",
"```\n",
"```{index} single: solver; cbc\n",
"```\n",
"```{index} single: AMPL; Set\n",
"```\n",
"```{index} pandas dataframe\n",
"```\n",
"# Recharging strategy for an electric vehicle\n",
"\n",
"Whether it is to visit family, take a sightseeing tour or call on business associates, planning a road trip is a familiar and routine task. Here we consider a road trip on a pre-determined route for which need to plan rest and recharging stops. This example demonstrates use of AMPL disjunctions to model the decisions on where to stop. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "700a2a08-73ab-411e-a5e0-8dafaf2fbc30",
"metadata": {},
"outputs": [],
"source": [
"# install dependencies and select solver\n",
"%pip install -q amplpy matplotlib\n",
"\n",
"SOLVER = \"cbc\"\n",
"\n",
"from amplpy import AMPL, ampl_notebook\n",
"\n",
"ampl = ampl_notebook(\n",
" modules=[\"cbc\"], # modules to install\n",
" license_uuid=\"default\", # license to use\n",
") # instantiate AMPL object and register magics"
]
},
{
"cell_type": "markdown",
"id": "8d93181c-3ee4-42f5-899b-3aef314bf1ea",
"metadata": {},
"source": [
"## Problem Statement\n",
"\n",
"Given the current location $x$, battery charge $c$, and planning horizon $D$, the task is to plan a series of recharging and rest stops. Data is provided for the location and the charging rate available at each charging stations. The objective is to drive from location $x$ to location $x + D$ in as little time as possible subject to the following constraints:\n",
"\n",
"* To allow for unforeseen events, the state of charge should never drop below 20% of the maximum capacity.\n",
"* The the maximum charge is $c_{max} = 80$ kWh.\n",
"* For comfort, no more than 4 hours should pass between stops, and that a rest stop should last at least $t^{rest}$.\n",
"* Any stop includes a $t^{lost} = 10$ minutes of \"lost time\".\n",
"\n",
"For this first model we make several simplifying assumptions that can be relaxed as a later time.\n",
"\n",
"* Travel is at a constant speed $v = 100$ km per hour and a constant discharge rate $R = 0.24$ kWh/km\n",
"* The batteries recharge at a constant rate determined by the charging station.\n",
"* Only consider stops at the recharging stations."
]
},
{
"cell_type": "markdown",
"id": "04e6c32e-7deb-4ea7-8df6-a24a5185ad4a",
"metadata": {},
"source": [
"## Modeling\n",
"\n",
"The problem statement identifies four state variables.\n",
"\n",
"* $c$ the current battery charge\n",
"* $r$ the elapsed time since the last rest stop\n",
"* $t$ elapsed time since the start of the trip\n",
"* $x$ the current location\n",
"\n",
"The charging stations are located at positions $d_i$ for $i\\in I$ with capacity $C_i$. The arrival time at charging station $i$ is given by\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"c_i^{arr} & = c_{i-1}^{dep} - R (d_i - d_{i-1}) \\\\\n",
"r_i^{arr} & = r_{i-1}^{dep} + \\frac{d_i - d_{i-1}}{v} \\\\\n",
"t_i^{arr} & = t_{i-1}^{dep} + \\frac{d_i - d_{i-1}}{v} \\\\\n",
"\\end{align*}\n",
"$$\n",
"\n",
"where the script $t_{i-1}^{dep}$ refers to departure from the prior location. At each charging location there is a decision to make of whether to stop, rest, and recharge. If the decision is positive, then\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"c_i^{dep} & \\leq c^{max} \\\\\n",
"r_i^{dep} & = 0 \\\\\n",
"t_i^{dep} & \\geq t_{i}^{arr} + t_{lost} + \\frac{c_i^{dep} - c_i^{arr}}{C_i} \\\\\n",
"t_i^{dep} & \\geq t_{i}^{arr} + t_{rest}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"which account for the battery charge, the lost time and time required for battery charging, and allows for a minimum rest time. On the other hand, if a decision is make to skip the charging and rest opportunity,\n",
"\n",
"$$\n",
"\\begin{align*}\n",
"c_i^{dep} & = c_i^{arr} \\\\\n",
"r_i^{dep} & = r_i^{arr} \\\\\n",
"t_i^{dep} & = t_i^{arr}\n",
"\\end{align*}\n",
"$$\n",
"\n",
"The latter sets of constraints have an exclusive-or relationship. That is, either one or the other of the constraint sets hold, but not both. \n",
"\n",
"$$\n",
"\\begin{align*}\n",
"\\min \\quad & t_{n+1}^{arr} \\\\\n",
"\\text{s.t.} \\quad\n",
" & r_i^{arr} \\leq r^{max} & \\forall \\, i \\in I \\\\\n",
" & c_i^{arr} \\geq c^{min} & \\forall \\,i \\in I \\\\\n",
" & c_i^{arr} = c_{i-1}^{dep} - R (d_i - d_{i-1}) & \\forall \\,i \\in I \\\\\n",
" & r_i^{arr} = r_{i-1}^{dep} + \\frac{d_i - d_{i-1}}{v} & \\forall \\,i \\in I \\\\\n",
" & t_i^{arr} = t_{i-1}^{dep} + \\frac{d_i - d_{i-1}}{v} & \\forall \\,i \\in I \\\\\n",
"& \\begin{bmatrix}\n",
" c_i^{dep} & \\leq & c^{max} \\\\\n",
" r_i^{dep} & = & 0 \\\\\n",
" t_i^{dep} & \\geq & t_{i}^{arr} + t_{lost} + \\frac{c_i^{dep} - c_i^{arr}}{C_i} \\\\\n",
" t_i^{dep} & \\geq & t_{i}^{arr} + t_{rest}\n",
"\\end{bmatrix}\n",
"\\veebar\n",
"\\begin{bmatrix}\n",
" c_i^{dep} = c_i^{arr} \\\\\n",
" r_i^{dep} = r_i^{arr} \\\\\n",
" t_i^{dep} = t_i^{arr}\n",
"\\end{bmatrix} & \\forall \\, i \\in I.\n",
"\\end{align*}\n",
"$$\n"
]
},
{
"cell_type": "markdown",
"id": "82aee043-4d9a-4108-8649-f3bea3dcf210",
"metadata": {},
"source": [
"## Charging Station Information"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "fff98145-1b94-42cf-9b0d-ec3593495a8c",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" name \n",
" location \n",
" kw \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" S_00 \n",
" 191.6 \n",
" 150 \n",
" \n",
" \n",
" 1 \n",
" S_01 \n",
" 310.6 \n",
" 100 \n",
" \n",
" \n",
" 2 \n",
" S_02 \n",
" 516.0 \n",
" 50 \n",
" \n",
" \n",
" 3 \n",
" S_03 \n",
" 683.6 \n",
" 50 \n",
" \n",
" \n",
" 4 \n",
" S_04 \n",
" 769.9 \n",
" 50 \n",
" \n",
" \n",
" 5 \n",
" S_05 \n",
" 869.7 \n",
" 100 \n",
" \n",
" \n",
" 6 \n",
" S_06 \n",
" 1009.1 \n",
" 150 \n",
" \n",
" \n",
" 7 \n",
" S_07 \n",
" 1164.7 \n",
" 100 \n",
" \n",
" \n",
" 8 \n",
" S_08 \n",
" 1230.8 \n",
" 100 \n",
" \n",
" \n",
" 9 \n",
" S_09 \n",
" 1350.8 \n",
" 250 \n",
" \n",
" \n",
" 10 \n",
" S_10 \n",
" 1508.4 \n",
" 100 \n",
" \n",
" \n",
" 11 \n",
" S_11 \n",
" 1639.8 \n",
" 100 \n",
" \n",
" \n",
" 12 \n",
" S_12 \n",
" 1809.4 \n",
" 150 \n",
" \n",
" \n",
" 13 \n",
" S_13 \n",
" 1947.3 \n",
" 250 \n",
" \n",
" \n",
" 14 \n",
" S_14 \n",
" 2145.2 \n",
" 150 \n",
" \n",
" \n",
" 15 \n",
" S_15 \n",
" 2337.5 \n",
" 100 \n",
" \n",
" \n",
" 16 \n",
" S_16 \n",
" 2415.6 \n",
" 100 \n",
" \n",
" \n",
" 17 \n",
" S_17 \n",
" 2590.0 \n",
" 100 \n",
" \n",
" \n",
" 18 \n",
" S_18 \n",
" 2691.2 \n",
" 100 \n",
" \n",
" \n",
" 19 \n",
" S_19 \n",
" 2896.2 \n",
" 100 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# current location (km) and charge (kw)\n",
"x = 0\n",
"\n",
"# planning horizon\n",
"D = 2000\n",
"\n",
"# visualize\n",
"fig, ax = plt.subplots(1, 1, figsize=(15, 3))\n",
"\n",
"\n",
"def plot_stations(stations, x, D, ax=ax):\n",
" for station in stations.index:\n",
" xs = stations.loc[station, \"location\"]\n",
" ys = stations.loc[station, \"kw\"]\n",
" ax.plot([xs, xs], [0, ys], \"b\", lw=10, solid_capstyle=\"butt\")\n",
" ax.text(xs, 0 - 30, stations.loc[station, \"name\"], ha=\"center\")\n",
"\n",
" ax.plot([x, x + D], [0, 0], \"r\", lw=5, solid_capstyle=\"butt\", label=\"plan horizon\")\n",
" ax.plot([x, x + D], [0, 0], \"r.\", ms=20)\n",
"\n",
" ax.axhline(0)\n",
" ax.set_ylim(-50, 300)\n",
" ax.set_xlabel(\"Distance\")\n",
" ax.set_ylabel(\"kw\")\n",
" ax.set_title(\"charging stations\")\n",
" ax.legend()\n",
"\n",
"\n",
"plot_stations(stations, x, D)"
]
},
{
"cell_type": "markdown",
"id": "91ebd206-ca91-436c-86f9-71f89bb1f34d",
"metadata": {},
"source": [
"## Car Information"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "2d6dad34-d50b-4acf-bb44-2fbf3e044086",
"metadata": {},
"outputs": [],
"source": [
"# charge limits (kw)\n",
"c_max = 150\n",
"c_min = 0.2 * c_max\n",
"c = c_max\n",
"\n",
"# velocity km/hr and discharge rate kwh/km\n",
"v = 100.0\n",
"R = 0.24\n",
"\n",
"# lost time\n",
"t_lost = 10 / 60\n",
"t_rest = 10 / 60\n",
"\n",
"# rest time\n",
"r_max = 3"
]
},
{
"cell_type": "markdown",
"id": "c797f5b1-c0a2-498c-9f94-2beaa947d961",
"metadata": {},
"source": [
"## AMPL Model"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "5cabe50f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Overwriting ev_plan.mod\n"
]
}
],
"source": [
"%%writefile ev_plan.mod\n",
"\n",
"param n;\n",
"\n",
"# locations and road segments between location x and x + D\n",
"set STATIONS; # 1..n\n",
"set LOCATIONS; # 0, 1..n, D\n",
"set SEGMENTS; # 1..n + 1\n",
"\n",
"param C{STATIONS};\n",
"param D;\n",
"param c_min;\n",
"param c_max;\n",
"param v;\n",
"param R;\n",
"\n",
"param r_max;\n",
"param location{LOCATIONS};\n",
"param dist{SEGMENTS};\n",
"param t_lost;\n",
"\n",
"# distance traveled\n",
"var x{LOCATIONS} >= 0, <= 10000;\n",
"\n",
"# arrival and departure charge at each charging station\n",
"var c_arr{LOCATIONS} >= c_min, <= c_max;\n",
"var c_dep{LOCATIONS} >= c_min, <= c_max;\n",
"\n",
"# arrival and departure times from each charging station\n",
"var t_arr{LOCATIONS} >= 0, <= 100;\n",
"var t_dep{LOCATIONS} >= 0, <= 100;\n",
"\n",
"# arrival and departure rest from each charging station\n",
"var r_arr{LOCATIONS} >= 0, <= r_max;\n",
"var r_dep{LOCATIONS} >= 0, <= r_max;\n",
"\n",
"minimize min_time: t_arr[n + 1];\n",
" \n",
"s.t. drive_time {i in SEGMENTS}: t_arr[i] == t_dep[i-1] + dist[i]/v;\n",
"s.t. rest_time {i in SEGMENTS}: r_arr[i] == r_dep[i-1] + dist[i]/v;\n",
"s.t. drive_distance {i in SEGMENTS}: x[i] == x[i-1] + dist[i];\n",
"s.t. discharge {i in SEGMENTS}: c_arr[i] == c_dep[i-1] - R * dist[i];\n",
"\n",
"s.t. recharge {i in STATIONS}:\n",
" # list of constraints that apply if there is no stop at station i\n",
" ((c_dep[i] == c_arr[i] and t_dep[i] == t_arr[i] and r_dep[i] == r_arr[i])\n",
" or\n",
" # list of constraints that apply if there is a stop at station i\n",
" (t_dep[i] == t_lost + t_arr[i] + (c_dep[i] - c_arr[i])/C[i] and\n",
" c_dep[i] >= c_arr[i] and r_dep[i] == 0))\n",
" and not\n",
" ((c_dep[i] == c_arr[i] and t_dep[i] == t_arr[i] and r_dep[i] == r_arr[i])\n",
" and\n",
" (t_dep[i] == t_lost + t_arr[i] + (c_dep[i] - c_arr[i])/C[i] and\n",
" c_dep[i] >= c_arr[i] and r_dep[i] == 0));"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "d061b189-5ccd-4949-a1b9-f736533e8389",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"cbc 2.10.7: \b\b\b\b\b\b\b\b\b\b\b\bcbc 2.10.7: optimal solution; objective 24.142507\n",
"12091 simplex iterations\n",
"12091 barrier iterations\n",
"102 branching nodes\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" location \n",
" t_arr \n",
" t_dep \n",
" c_arr \n",
" c_dep \n",
" t_stop \n",
" \n",
" \n",
" index \n",
" \n",
" \n",
" \n",
" \n",
" 0 \n",
" 0.0 \n",
" 0.000000 \n",
" 0.000000 \n",
" 30.0000 \n",
" 150.0000 \n",
" 0.000000 \n",
" \n",
" \n",
" 1 \n",
" 191.6 \n",
" 1.916000 \n",
" 2.389227 \n",
" 104.0160 \n",
" 150.0000 \n",
" 0.473227 \n",
" \n",
" \n",
" 2 \n",
" 310.6 \n",
" 3.579227 \n",
" 4.031493 \n",
" 121.4400 \n",
" 150.0000 \n",
" 0.452267 \n",
" \n",
" \n",
" 3 \n",
" 516.0 \n",
" 6.085493 \n",
" 6.535840 \n",
" 100.7040 \n",
" 114.8880 \n",
" 0.450347 \n",
" \n",
" \n",
" 4 \n",
" 683.6 \n",
" 8.211840 \n",
" 8.211840 \n",
" 74.6640 \n",
" 74.6640 \n",
" 0.000000 \n",
" \n",
" \n",
" 5 \n",
" 769.9 \n",
" 9.074840 \n",
" 9.241507 \n",
" 53.9520 \n",
" 53.9520 \n",
" 0.166667 \n",
" \n",
" \n",
" 6 \n",
" 869.7 \n",
" 10.239507 \n",
" 10.740733 \n",
" 30.0000 \n",
" 63.4560 \n",
" 0.501227 \n",
" \n",
" \n",
" 7 \n",
" 1009.1 \n",
" 12.134733 \n",
" 12.848120 \n",
" 30.0000 \n",
" 112.0080 \n",
" 0.713387 \n",
" \n",
" \n",
" 8 \n",
" 1164.7 \n",
" 14.404120 \n",
" 14.570787 \n",
" 74.6640 \n",
" 74.6640 \n",
" 0.166667 \n",
" \n",
" \n",
" 9 \n",
" 1230.8 \n",
" 15.231787 \n",
" 15.231787 \n",
" 58.8000 \n",
" 58.8000 \n",
" 0.000000 \n",
" \n",
" \n",
" 10 \n",
" 1350.8 \n",
" 16.431787 \n",
" 17.078453 \n",
" 30.0000 \n",
" 150.0000 \n",
" 0.646667 \n",
" \n",
" \n",
" 11 \n",
" 1508.4 \n",
" 18.654453 \n",
" 18.654453 \n",
" 112.1760 \n",
" 112.1760 \n",
" -0.000000 \n",
" \n",
" \n",
" 12 \n",
" 1639.8 \n",
" 19.968453 \n",
" 20.135121 \n",
" 80.6400 \n",
" 80.6401 \n",
" 0.166668 \n",
" \n",
" \n",
" 13 \n",
" 1809.4 \n",
" 21.831121 \n",
" 22.236507 \n",
" 39.9361 \n",
" 75.7440 \n",
" 0.405386 \n",
" \n",
" \n",
" 14 \n",
" 1947.3 \n",
" 23.615507 \n",
" 23.615507 \n",
" 42.6480 \n",
" 42.6480 \n",
" 0.000000 \n",
" \n",
" \n",
" 15 \n",
" 2000.0 \n",
" 24.142507 \n",
" 0.000000 \n",
" 30.0000 \n",
" 30.0000 \n",
" -24.142507 \n",
" \n",
" \n",
"
\n",
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# visualize\n",
"\n",
"\n",
"def visualize(m):\n",
" D = m.param[\"D\"].value()\n",
"\n",
" results = get_results(m)\n",
"\n",
" results[\"t_stop\"] = results[\"t_dep\"] - results[\"t_arr\"]\n",
"\n",
" fig, ax = plt.subplots(2, 1, figsize=(15, 8), sharex=True)\n",
"\n",
" # plot stations\n",
" for station in stations.index:\n",
" xs = stations.loc[station, \"location\"]\n",
" ys = stations.loc[station, \"kw\"]\n",
" ax[0].plot([xs, xs], [0, ys], \"b\", lw=10, solid_capstyle=\"butt\")\n",
" ax[0].text(xs, 0 - 30, stations.loc[station, \"name\"], ha=\"center\")\n",
"\n",
" # plot planning horizon\n",
" ax[0].plot(\n",
" [x, x + D], [0, 0], \"r\", lw=5, solid_capstyle=\"butt\", label=\"plan horizon\"\n",
" )\n",
" ax[0].plot([x, x + D], [0, 0], \"r.\", ms=20)\n",
"\n",
" # annotations\n",
" ax[0].axhline(0)\n",
" ax[0].set_ylim(-50, 300)\n",
" ax[0].set_ylabel(\"kw\")\n",
" ax[0].set_title(\"charging stations\")\n",
" ax[0].legend()\n",
"\n",
" SEGMENTS = m.set[\"SEGMENTS\"].to_list()\n",
"\n",
" # plot battery charge\n",
" for i in SEGMENTS:\n",
" xv = [results.loc[i - 1, \"location\"], results.loc[i, \"location\"]]\n",
" cv = [results.loc[i - 1, \"c_dep\"], results.loc[i, \"c_arr\"]]\n",
" ax[1].plot(xv, cv, \"g\")\n",
"\n",
" STATIONS = m.set[\"STATIONS\"].to_list()\n",
"\n",
" # plot charge at stations\n",
" for i in STATIONS:\n",
" xv = [results.loc[i, \"location\"]] * 2\n",
" cv = [results.loc[i, \"c_arr\"], results.loc[i, \"c_dep\"]]\n",
" ax[1].plot(xv, cv, \"g\")\n",
"\n",
" # mark stop locations\n",
" for i in STATIONS:\n",
" if results.loc[i, \"t_stop\"] > 0:\n",
" ax[1].axvline(results.loc[i, \"location\"], color=\"r\", ls=\"--\")\n",
"\n",
" # show constraints on battery charge\n",
" ax[1].axhline(c_max, c=\"g\")\n",
" ax[1].axhline(c_min, c=\"g\")\n",
" ax[1].set_ylim(0, 1.1 * c_max)\n",
" ax[1].set_ylabel(\"Charge (kw)\")\n",
"\n",
"\n",
"visualize(ev_plan(stations, 0, 2000))"
]
},
{
"cell_type": "markdown",
"id": "479452ec-be4e-46e4-baa8-4f1e9de69968",
"metadata": {},
"source": [
"## Suggested Exercises\n",
"\n",
"1. Does increasing the battery capacity $c^{max}$ significantly reduce the time required to travel 2000 km? Explain what you observe.\n",
"\n",
"2. \"The best-laid schemes of mice and men go oft awry\" (Robert Burns, \"To a Mouse\"). Modify this model so that it can be used to update a plans in response to real-time measurements. How does the charging strategy change as a function of planning horizon $D$?"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}