diff --git a/optimization201/Modeling_Session_2/completed_modeling2.ipynb b/optimization201/Modeling_Session_2/completed_modeling2.ipynb index 77961a1..fa2b773 100644 --- a/optimization201/Modeling_Session_2/completed_modeling2.ipynb +++ b/optimization201/Modeling_Session_2/completed_modeling2.ipynb @@ -1 +1 @@ -{"cells":[{"cell_type":"markdown","metadata":{"id":"ezGOS9cRadzN"},"source":["# Intermediate Modeling - Part 2: Back for More Sun\n","Just as we expanded on the first modeling example from Opti 101 for the first example of Opti 201, we'll do the same here. This example is motivated by [IEEE's Predict+Optimize Technical Challenge](https://ieee-dataport.org/competitions/ieee-cis-technical-challenge-predictoptimize-renewable-energy-scheduling) in 2021.\n","\n","In this example we are trying to minimize the amount of electricity purchased from the grid to supply a building with electricity over a 6-day period while meeting predicted demand. The building is outfitted with solar panels and two batteries to store and provide energy when appropriate. Decisions made at each 30-minute interval are:\n","1. Buy electricity from the grid\n","2. Charge/discharge the one battery\n","3. Charge/discharge the other battery\n","\n","The output of the optimization model is a charge/discharge schedule for the batteries as well as the time periods where electricity will be purchased. The demand is forecasted using [Prophet](https://facebook.github.io/prophet/). The notebook for this part of the problem can be found [here](https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/optimization101/Modeling_Session_2/energy_storage_ML.ipynb) If you want to have the first notebook handy, you can find it [here](https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/optimization101/Modeling_Session_2/completed_modeling2.ipynb).\n","\n","This model is a lot more complicated that the first from Opti 101, so we can copy the formulation and code in a concise way and highlight the parts relevant to the features discussed in the second modeling session video.\n","\n","In subsequent sections we'll see how to adapt this model to incorporate\n","- Multi-objective models,\n","- Multi-scenario models, and\n","- Solutions pools"]},{"cell_type":"markdown","metadata":{"id":"6uEfx_64u0jY"},"source":["## Opti 101 Model"]},{"cell_type":"markdown","metadata":{"id":"E5lo-xCHu4Jk"},"source":["### Formulation"]},{"cell_type":"markdown","metadata":{"id":"CXPsSNehu8Nu"},"source":["#### Sets and Input Parameters\n","\n","Batteries: $b \\in B = \\{\\texttt{Battery0, Battery1}\\}$\n","\n","Time periods: $t \\in T = \\{0,1,...,179\\}$ for the first Monday through Saturday of October, 2020\n","\n","$c_{b}$: capacity of battery $b \\in B \\quad\\quad \\texttt{capacity[b]}$\n","\n","$p_{b}$: loss of energy (as a percentage) during transfer into battery $b\\in B \\quad\\quad \\texttt{p}\\_\\texttt{loss[b]}$\n","\n","$q_{b}$: quantity of initial energy in battery $b \\in B \\quad\\quad \\texttt{initial[b]}$\n","\n","$solar_{t}$: solar power generation of the panel for time period $t \\in T \\quad\\quad \\texttt{solar}\\_\\texttt{values[b]}$\n","\n","$d_{t}$: total building and class energy demand for time period $t\\in T \\quad\\quad \\texttt{total}\\_\\texttt{deamnd[t]}$"]},{"cell_type":"markdown","metadata":{"id":"bOQnhuPkvIcg"},"source":["#### Decision Variables\n","Let $f^{in}_{b,t}$ be the amount each battery, $b$, `charges` at time period, $t$, $\\forall b\\in B, t\\in T$. $\\quad\\quad \\texttt{flow}\\_\\texttt{in[b,t]}$\n","\n","Let $f^{out}_{b,t}$ be the amount each battery `discharges`, similarly. $\\quad\\quad \\texttt{flow}\\_\\texttt{out[b,t]}$\n","\n","The max amount that each battery can charge or discharge in a single period to be 20 kW.\n","\n","$s_{b,t}$ is the current amount of energy in battery, $b$, at the end of time period, $t, \\forall b\\in B, t\\in T$. $\\quad\\texttt{state[b,t]}$\n","\n","$gen_{t}$ is the amount of available solar energy that is used in period $t$, $\\forall t \\in T$. $\\quad\\texttt{gen[t]}$\n","\n","$grid_{t}$ : This variable indicates the amount of energy purchased from the grid at time period, $t$, $\\forall t \\in T$. $\\quad\\texttt{grid[t]}$\n","\n","We are not considering the ability to \"sell back\" electricity, so $grid_t \\ge 0$."]},{"cell_type":"markdown","metadata":{"id":"jZ5Zbv6PF6xU"},"source":["#### Constraints\n","**Power balance:**\n","$$\n","\\begin{align}\n","\\sum_b(f^{out}_{b,t}-p_bf^{in}_{b,t}) + gen_t + grid_t = d_t \\quad \\forall t \\in T\n","\\end{align}\n","$$\n","\n","**Battery state:**\n","\n","First time period's charge/discharge decision.\n","\n","\\begin{equation}\n","s_{b,0} = q_b + p_bf^{in}_{b,0} - f^{out}_{b,0}\n","\\end{equation}\n","\n","For each time period after (i.e. $t\\ge1$), the state of a battery is found by:\n","\n","\\begin{equation}\n","s_{b,t} = s_{b,t-1} + p_bf^{in}_{b,t} - f^{out}_{b,t}, t \\ge 1\n","\\end{equation}\n","\n","**Solar availability:**\n","$$\n","\\begin{equation}\n","f^{in}_{\\texttt{Battery0},t} + f^{in}_{\\texttt{Battery1},t} + gen_t \\le solar_t, \\quad \\forall t \\in T\n","\\end{equation}\n","$$\n","\n","**Charge/discharge:**\n","We cannot charge and discharge the same battery in the same time period.\n","Define $z_t \\in \\{0,1\\}, t \\in T$.\n","$$\n","\\begin{align*}\n","f^{in}_{b,t} &\\leq 20*z_{b,t} &\\forall b \\in B, t \\in T \\\\\n","f^{out}_{b,t} &\\leq 20*(1-z_{b,t}) &\\forall b \\in B, t \\in T\n","\\end{align*}\n","$$\n","Hmm, this set of constraints looks a little familiar, huh?\n","\n","**Battery capacity:**\n","$$\n","\\begin{align*}\n","s_{b,t} \\le c_b, \\quad \\forall b \\in B, t\\in T\n","\\end{align*}\n","$$"]},{"cell_type":"markdown","metadata":{"id":"7j--9VCkuqCt"},"source":["#### Objective Function\n","The objective is to minimize the total amount of energy purchased from the grid over all time periods.\n","\n","\\begin{equation}\n","{\\rm min} \\space \\sum_{t} g_{t}\n","\\end{equation}"]},{"cell_type":"markdown","metadata":{"id":"dpyvjhp4oT4N"},"source":["### Gurobipy Code"]},{"cell_type":"markdown","metadata":{"id":"8EBFFKNPvft2"},"source":["#### Sets and Input parameters"]},{"cell_type":"code","execution_count":1,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":15980,"status":"ok","timestamp":1699599941524,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"iPw18bDQF9Lc","outputId":"4f3ddf9a-403a-4536-d766-0086b80c0423"},"outputs":[{"name":"stdout","output_type":"stream","text":["Collecting gurobipy\n"," Downloading gurobipy-10.0.3-cp310-cp310-manylinux2014_x86_64.whl (12.7 MB)\n","\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m12.7/12.7 MB\u001b[0m \u001b[31m20.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n","\u001b[?25hInstalling collected packages: gurobipy\n","Successfully installed gurobipy-10.0.3\n","Total Solar Generation: 4939.172 \n","Total Demand: 5250.6\n"]}],"source":["# install/import packages\n","%pip install gurobipy\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import gurobipy as gp\n","from gurobipy import GRB\n","\n","# read in solar forecast data\n","path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_2/'\n","solar_values_read = pd.read_csv(path+'pred_solar_values.csv')\n","solar_values = round(solar_values_read.yhat,3)\n","solar_values.reset_index(drop = True, inplace = True)\n","\n","# read in demaind data, where total demand is a fixed building demand and\n","# an estimated demand based on a proposed schedule for the building\n","schedule = pd.read_csv(path+'schedule_demand.csv')\n","avg_building = pd.read_csv(path+'building_demand.csv')\n","total_demand = schedule.sched_demand + avg_building.build_demand\n","print(f\"Total Solar Generation: {solar_values.sum()} \\nTotal Demand: {total_demand.sum()}\")\n","\n","# create data for batteries, including capity, efficency and initial charges\n","# also define time periods,\n","batteries = [\"Battery0\", \"Battery1\"]\n","capacity = {\"Battery0\": 60, \"Battery1\": 80} # in Kw\n","p_loss = {\"Battery0\": 0.95, \"Battery1\": 0.9} # proportion\n","initial = {\"Battery0\": 0, \"Battery1\": 0} # in kW\n","time_periods = range(len(solar_values_read))"]},{"cell_type":"markdown","metadata":{"id":"20u3pAiAulYD"},"source":["#### Decison Variables"]},{"cell_type":"code","execution_count":113,"metadata":{"executionInfo":{"elapsed":101,"status":"ok","timestamp":1699602343575,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"mpi8sWtKqEMw"},"outputs":[],"source":["m = gp.Model() #this defines the model that we'll add to as we finish the formulation\n","\n","flow_in = m.addVars(batteries, time_periods, name=\"flow_in\")\n","flow_out = m.addVars(batteries, time_periods, name=\"flow_out\")\n","\n","grid = m.addVars(time_periods, name=\"grid\")\n","state = m.addVars(batteries, time_periods, ub = [capacity[b] for b in batteries for t in time_periods], name=\"state\")\n","gen = m.addVars(time_periods, name=\"gen\")\n","\n","zwitch = m.addVars(batteries, time_periods, vtype=GRB.BINARY, name=\"zwitch\")"]},{"cell_type":"markdown","metadata":{"id":"pJIGt-hLqgWi"},"source":["#### Constraints\n","Note than the **capacity constraint** has been taken care of by putting an `ub` on the `state` decision variable."]},{"cell_type":"code","execution_count":114,"metadata":{"executionInfo":{"elapsed":138,"status":"ok","timestamp":1699602344717,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"7vnVisjTqiHL"},"outputs":[],"source":["# Power balance\n","m.addConstrs((gp.quicksum(flow_out[b,t] - p_loss[b]*flow_in[b,t] for b in batteries) + gen[t] + grid[t] == total_demand[t]\n"," for t in time_periods), name=\"power_balance\");\n","# Battery state\n","m.addConstrs((state[b,0] == initial[b] + p_loss[b]*flow_in[b,0] - flow_out[b,0] for b in batteries), name=\"initial_state\")\n","m.addConstrs((state[b,t] == state[b,t-1] + p_loss[b]*flow_in[b,t] - flow_out[b,t] for b in batteries for t in time_periods if t >= 1), name=\"subsequent_states\");\n","\n","#Solar availability\n","m.addConstrs((flow_in['Battery0',t] + flow_in['Battery1',t] + gen[t] <= solar_values[t] for t in time_periods), name = \"solar_avail\");\n","\n","#Charge/discharge\n","m.addConstrs((flow_in[b,t] <= 20*zwitch[b,t] for b in batteries for t in time_periods), name = \"to_charge\")\n","m.addConstrs((flow_out[b,t] <= 20*(1-zwitch[b,t]) for b in batteries for t in time_periods), name = \"or_not_to_charge\");\n"]},{"cell_type":"markdown","metadata":{"id":"AqsU88pnqpiQ"},"source":["#### Objective Function\n","The original objective is to minimize the amount of electricity purchased from the grid. There was also a section where we considered a different objective -- given expected electricity prices for each period, minimize the total cost. So we can compare the cost and amount of electricity purchased across scenarios, let's start to capture these values."]},{"cell_type":"code","execution_count":115,"metadata":{"executionInfo":{"elapsed":253,"status":"ok","timestamp":1699602346323,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"q-yYYpvEq0OJ"},"outputs":[],"source":["# read in estimated price of electricity for each time period\n","avg_price = pd.read_csv(path+'expected_price.csv')\n","price = avg_price.price\n","\n","# define a linear expression for total energy purchased from the grid\n","total_grid = grid.sum()\n","# define a linear expression for total cost\n","total_cost = gp.quicksum(avg_price.price[t]*grid[t] for t in time_periods)\n","\n","# set the objective to the original of min grid purchase\n","m.setObjective(total_grid, GRB.MINIMIZE)"]},{"cell_type":"markdown","metadata":{"id":"uYDpiO-2v7dD"},"source":["### Solve the Model"]},{"cell_type":"code","execution_count":116,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":686},"executionInfo":{"elapsed":190,"status":"ok","timestamp":1699602348004,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"NbNVPTBZwC8E","outputId":"6ed9a448-59d5-46ba-b589-385f29380ce6"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0xaa31bb6a\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [1e+00, 1e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Found heuristic solution: objective 5250.6000000\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Found heuristic solution: objective 4986.0960000\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 1.281678e+03, 594 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n","* 0 0 0 1281.6776848 1281.67768 0.00% - 0s\n","\n","Explored 1 nodes (594 simplex iterations) in 0.08 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 3: 1281.68 4986.1 5250.6 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.281677684779e+03, best bound 1.281677684779e+03, gap 0.0000%\n"]},{"data":{"text/html":["\n","
\n","
\n","\n","\n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n","
CostGridAmount
0803.181281.68
\n","
\n","
\n","\n","
\n"," \n","\n"," \n","\n"," \n","
\n","\n","
\n","
\n"],"text/plain":[" Cost GridAmount\n","0 803.18 1281.68"]},"execution_count":116,"metadata":{},"output_type":"execute_result"}],"source":["m.optimize()\n","\n","results = pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]],\n"," columns = ['Cost','GridAmount']\n"," )\n","\n","results"]},{"cell_type":"markdown","metadata":{"id":"TTBd1ssUpRxX"},"source":["The original objective is to minimize the amount of electricity purchased from the grid. There was also a section where we considered a different objective -- given expected electricity prices for each period, minimize the total cost. We'll read in the data and quickly update the objective and resolve."]},{"cell_type":"code","execution_count":117,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":228,"status":"ok","timestamp":1699602350218,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"vBeoWLq6pcpO","outputId":"ffdedc1b-0b30-4460-c837-639f5f8862fa"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x66bddde2\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","MIP start from previous solve produced solution with objective 617.569 (0.04s)\n","MIP start from previous solve produced solution with objective 617.569 (0.04s)\n","Loaded MIP start from previous solve with objective 617.569\n","\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.01s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 6.019761e+02, 754 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 601.97612 0 1 617.56942 601.97612 2.52% - 0s\n","H 0 0 601.9761241 601.97612 0.00% - 0s\n"," 0 0 601.97612 0 1 601.97612 601.97612 0.00% - 0s\n","\n","Explored 1 nodes (754 simplex iterations) in 0.11 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 2: 601.976 617.569 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240990e+02, gap 0.0000%\n"]}],"source":["m.setObjective(total_cost, GRB.MINIMIZE)\n","m.optimize()\n","\n","results = pd.concat([results,\n"," pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]], columns = ['Cost','GridAmount'])\n"," ],\n"," ignore_index=True\n"," )\n","\n","results\n","#create a copy for later\n","mp = m.copy()"]},{"cell_type":"markdown","metadata":{"id":"pU48SBLRxjun"},"source":["## Let's Find Some Hidden Gems"]},{"cell_type":"markdown","metadata":{"id":"_7L8k0t-0cgI"},"source":["### Handling Multiple Objectives\n","\n","The above results show two different scenarios, one that purely considers grid purchases, and the other purely prioritizes cost. Let's see how to use Gurobi's [multi-objective](https://www.gurobi.com/documentation/current/refman/working_with_multiple_obje.html) functionality to combine the two objectives by first implementing a `weighed` setup."]},{"cell_type":"markdown","metadata":{"id":"AOPH_t4eDAiV"},"source":["#### Weighted Objectives\n","In using a weighted objective, we need to obviously decide the `weights`. This can be a little tricky given units and scale."]},{"cell_type":"code","execution_count":82,"metadata":{"executionInfo":{"elapsed":119,"status":"ok","timestamp":1699601541125,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"H4ntqJhexqZY"},"outputs":[],"source":["m.setObjectiveN(total_cost, index=0, weight = 1, name=\"cost\")\n","m.setObjectiveN(total_grid, index=1, weight = 10, name= \"grid\")"]},{"cell_type":"markdown","metadata":{"id":"G6muNC5D0cUT"},"source":["Now we can run the optimization with our weighted objective."]},{"cell_type":"code","execution_count":83,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":1000},"executionInfo":{"elapsed":473,"status":"ok","timestamp":1699601542850,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"frdBl5DM0Zg4","outputId":"312fbaf3-fa93-4efe-c68a-8382b8cfd049"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x9cf62986\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: starting optimization with 2 objectives (1 combined) ...\n","---------------------------------------------------------------------------\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 1 (weighted) ...\n","---------------------------------------------------------------------------\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0xf11f89cc\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [1e+01, 1e+01]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Found heuristic solution: objective 55228.437000\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Found heuristic solution: objective 52452.392480\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 1.343869e+04, 703 iterations, 0.04 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 13438.6902 0 2 52452.3925 13438.6902 74.4% - 0s\n","H 0 0 13438.690179 13438.6902 0.00% - 0s\n"," 0 0 13438.6902 0 2 13438.6902 13438.6902 0.00% - 0s\n","\n","Explored 1 nodes (703 simplex iterations) in 0.14 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 3: 13438.7 52452.4 55228.4 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.343869017947e+04, best bound 1.343869017947e+04, gap 0.0000%\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: solved in 0.16 seconds (0.03 work units), solution count 3\n","\n"]},{"data":{"text/html":["\n","
\n","
\n","\n","\n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n","
CostGridAmount
0803.181281.68
1601.981317.69
2616.701282.20
\n","
\n","
\n","\n","
\n"," \n","\n"," \n","\n"," \n","
\n","\n","\n","
\n"," \n","\n","\n","\n"," \n","
\n","
\n","
\n"],"text/plain":[" Cost GridAmount\n","0 803.18 1281.68\n","1 601.98 1317.69\n","2 616.70 1282.20"]},"execution_count":83,"metadata":{},"output_type":"execute_result"}],"source":["m.ModelSense = GRB.MINIMIZE\n","m.optimize()\n","\n","results = pd.concat([results,\n"," pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]], columns = ['Cost','GridAmount'])\n"," ],\n"," ignore_index=True\n"," )\n","\n","results"]},{"cell_type":"markdown","metadata":{"id":"f1HiiJcZBLws"},"source":["Using this multi-objective approach, we have been able to significantly reduce the cost from the first solution* while adding on a very small amount of grid purchase.\n","\n","*Note that since we didn't specify anything in the first model that pertained to price there may be a solution that uses the same amout of purchased energy at a different cost.\n"]},{"cell_type":"markdown","metadata":{"id":"UcoCaW7zDH-H"},"source":["#### Hierarchical Objectives\n","This kind of multi-objective approach gives priority to a given objective in a given order, and for subsequent objectives, constraints are added to make sure the expressions of prior objectives to not degrade beyond given values.\n","\n","In the `hierarchical` approach the two main arguments are `abstol` and `reltol`. The former allows the first objective to degrade by a *raw* amount, while the latter is a *percentage*.\n","\n","For example, we can set first optimize for minimal total cost, and then allow that objective value to *degrade* by $50 or by a factor of 0.1 (i.e. increase by 10%) when we then minimize for the lowest amount of electricity purchased.\n","\n","An important concept in energy storage in maintaining battery health. One way this may be handled is to limit a battery's [depth of discharge](https://en.wikipedia.org/wiki/Depth_of_discharge), which can be expressed as the percentage of a battery's capacity that has been discharged. The way we decide to do this is by limiting the number of time periods our batteries are discharged below a certain depth.\n","\n","To do this we'll need binary variables to indicate when this occurs.\n","\n","Let $v_{b,t} = 1$ when battery $b$ is below a depth of discharge α (say 0.3) and 0 otherwise. We want to model\n","\n","\"If the battery gets below 30% it's capacity, then indicate depth of discharge violation.\"\n","\n","Let's replace some what comes after \"If\" and \"then\" with decision variables and (in)equalities. As a reminder $s_{b,t}$ is the current charge in the battery and $c_b$ is its capacity.\n","\n","$$\n","If \\space s_{b,t} < 0.3 \\times c_b, \\space then \\space v_{b,t} = 1\n","$$\n","\n","We can easily model this by taking the contrapositive of the conditional statement above and using indicator constraints.\n","\n","$$\n","If \\space v_{b,t} = 0, \\space then \\space s_{b,t} \\ge 0.3 \\times c_b,\n","$$\n","\n","````\n","To keep the problem within the 2,000 variable limit that comes with the `pip` install license,\n","this will only be implemented for `Battery0`.\n","````"]},{"cell_type":"code","execution_count":84,"metadata":{"executionInfo":{"elapsed":201,"status":"ok","timestamp":1699601547810,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"XgMclZuoAq3q"},"outputs":[],"source":["# Reminder that s_{b,t} is defined by state[b,t] in the code\n","v = m.addVars(time_periods, vtype=GRB.BINARY, name = 'v')\n","\n","# define a linear expression for total depth count\n","total_depth_count = v.sum()\n","\n","m.addConstrs((((v[t] == 0) >> (state['Battery0',t] >= 0.3*capacity['Battery0'])) for b in batteries for t in time_periods), name = \"discharge_depth\")\n","m.update()"]},{"cell_type":"markdown","metadata":{"id":"bXnwpJx1Ex85"},"source":["Now we'll dive into a three-piece objective."]},{"cell_type":"code","execution_count":85,"metadata":{"executionInfo":{"elapsed":97,"status":"ok","timestamp":1699601549020,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"jdVLoqVfHi4C"},"outputs":[],"source":["#the higher priority number is used first as an objective\n","m.setObjectiveN(total_cost, index=0, priority=2, reltol = 0.05, name=\"cost\")\n","m.setObjectiveN(total_depth_count, index=1, priority=1, reltol = 0.2, name= \"depth\")\n","m.setObjectiveN(total_grid, index=2, priority=0, name= \"grid\")"]},{"cell_type":"code","execution_count":86,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":1164,"status":"ok","timestamp":1699601551154,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"8EddQ72UVif6","outputId":"0419d2ca-b72d-4817-c5f8-281eacb3f1dc"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1980 columns and 4498 nonzeros\n","Model fingerprint: 0x4086b7a9\n","Model has 360 general constraints\n","Variable types: 1440 continuous, 540 integer (540 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n"," GenCon rhs range [2e+01, 2e+01]\n"," GenCon coe range [1e+00, 1e+00]\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: starting optimization with 3 objectives ... \n","---------------------------------------------------------------------------\n","\n","Multi-objectives: applying initial presolve ...\n","---------------------------------------------------------------------------\n","\n","Presolve added 305 rows and 302 columns\n","Presolve time: 0.02s\n","Presolved: 1745 rows and 2282 columns\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 1 (cost) ...\n","---------------------------------------------------------------------------\n","\n","Presolve removed 389 rows and 671 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1611 columns, 4160 nonzeros\n","Variable types: 1275 continuous, 336 integer (336 binary)\n","Found heuristic solution: objective 1021.0681360\n","\n","Root relaxation: objective 6.019761e+02, 697 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 601.97612 0 3 1021.06814 601.97612 41.0% - 0s\n","H 0 0 601.9761241 601.97612 0.00% - 0s\n"," 0 0 601.97612 0 3 601.97612 601.97612 0.00% - 0s\n","\n","Explored 1 nodes (697 simplex iterations) in 0.15 seconds (0.04 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 2: 601.976 1021.07 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240991e+02, gap 0.0000%\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 2 (depth) ...\n","---------------------------------------------------------------------------\n","\n","\n","Loaded user MIP start with objective 180\n","\n","Presolve removed 218 rows and 398 columns\n","Presolve time: 0.04s\n","Presolved: 1528 rows, 1884 columns, 4802 nonzeros\n","Variable types: 1373 continuous, 511 integer (511 binary)\n","Found heuristic solution: objective 66.0000000\n","\n","Root relaxation: objective 9.475247e+00, 1721 iterations, 0.05 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 9.47525 0 8 66.00000 9.47525 85.6% - 0s\n","H 0 0 14.0000000 9.47525 32.3% - 0s\n"," 0 0 10.03868 0 8 14.00000 10.03868 28.3% - 0s\n"," 0 0 10.26646 0 7 14.00000 10.26646 26.7% - 0s\n"," 0 0 10.30054 0 9 14.00000 10.30054 26.4% - 0s\n"," 0 0 10.31001 0 10 14.00000 10.31001 26.4% - 0s\n"," 0 0 10.31001 0 10 14.00000 10.31001 26.4% - 0s\n"," 0 0 10.41407 0 10 14.00000 10.41407 25.6% - 0s\n"," 0 0 10.45727 0 10 14.00000 10.45727 25.3% - 0s\n"," 0 0 10.45727 0 10 14.00000 10.45727 25.3% - 0s\n","H 0 0 13.0000000 10.45727 19.6% - 0s\n","H 0 0 12.0000000 10.45727 12.9% - 0s\n"," 0 2 10.71604 0 4 12.00000 10.71604 10.7% - 0s\n","\n","Cutting planes:\n"," Gomory: 2\n"," MIR: 12\n"," Flow cover: 9\n","\n","Explored 5 nodes (1802 simplex iterations) in 0.54 seconds (0.14 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 5: 12 13 14 ... 180\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.200000000000e+01, best bound 1.200000000000e+01, gap 0.0000%\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 3 (grid) ...\n","---------------------------------------------------------------------------\n","\n","\n","Loaded user MIP start with objective 1343.89\n","\n","Presolve removed 218 rows and 398 columns\n","Presolve time: 0.04s\n","Presolved: 1529 rows, 1884 columns, 4977 nonzeros\n","Variable types: 1373 continuous, 511 integer (511 binary)\n","\n","Root relaxation: objective 1.292876e+03, 1450 iterations, 0.03 seconds (0.01 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 1292.87646 0 13 1343.88759 1292.87646 3.80% - 0s\n"," 0 0 1302.07834 0 10 1343.88759 1302.07834 3.11% - 0s\n"," 0 0 1302.29363 0 12 1343.88759 1302.29363 3.10% - 0s\n"," 0 0 1302.91971 0 11 1343.88759 1302.91971 3.05% - 0s\n"," 0 0 1302.91971 0 12 1343.88759 1302.91971 3.05% - 0s\n"," 0 0 1303.20394 0 9 1343.88759 1303.20394 3.03% - 0s\n"," 0 0 1303.22949 0 9 1343.88759 1303.22949 3.03% - 0s\n"," 0 0 1303.23331 0 10 1343.88759 1303.23331 3.03% - 0s\n"," 0 0 1303.60134 0 12 1343.88759 1303.60134 3.00% - 0s\n"," 0 0 1303.60694 0 12 1343.88759 1303.60694 3.00% - 0s\n"," 0 0 1303.63927 0 12 1343.88759 1303.63927 2.99% - 0s\n"," 0 0 1303.63927 0 12 1343.88759 1303.63927 2.99% - 0s\n"," 0 2 1304.18118 0 12 1343.88759 1304.18118 2.95% - 0s\n","* 11 9 8 1313.8029332 1307.01963 0.52% 12.6 0s\n","* 43 7 8 1313.7843075 1309.46527 0.33% 13.1 1s\n","* 45 0 6 1309.5801997 1309.58020 0.00% 12.8 1s\n","\n","Cutting planes:\n"," Gomory: 4\n"," Implied bound: 2\n"," MIR: 12\n"," Flow cover: 13\n","\n","Explored 46 nodes (2108 simplex iterations) in 1.02 seconds (0.28 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 4: 1309.58 1313.78 1313.8 1343.89 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.309580199681e+03, best bound 1.309580199681e+03, gap 0.0000%\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: solved in 1.04 seconds (0.28 work units), solution count 9\n","\n"]}],"source":["m.optimize()"]},{"cell_type":"markdown","metadata":{"id":"TIMY-FXjn4LA"},"source":["We defined linear expressions for each part of the objective so we can use `getValue` to get their respective results. You can also query objectives as shown here."]},{"cell_type":"code","execution_count":87,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":89,"status":"ok","timestamp":1699601553769,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"jC4X0Bab_B5H","outputId":"ba1e01cd-ca50-4949-9274-0546ef669592"},"outputs":[{"name":"stdout","output_type":"stream","text":[" 632.14 14.0 1309.58"]}],"source":["for o in range(m.NumObj):\n"," # set which objective we will query\n"," m.params.ObjNumber = o\n"," # query the o-th objective value\n"," print(' ',round(m.ObjNVal,2), end='')"]},{"cell_type":"markdown","metadata":{"id":"nT1gJ5w4YJO4"},"source":["Play around with the priorities and tolerances to see how those change the overall solution."]},{"cell_type":"markdown","metadata":{"id":"fSw8RpaTY2CT"},"source":["### Handling Multiple Scenarios\n","In this section we will start with out **base model** with the objective being to minimize the total electricity purchase cost. To make is easier to edit constraints we can store."]},{"cell_type":"code","execution_count":88,"metadata":{"executionInfo":{"elapsed":283,"status":"ok","timestamp":1699601555555,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"XGsWkqaUnTbj"},"outputs":[],"source":["mm = gp.Model('mulit_scen') #this defines the model that we'll add to as we finish the formulation\n","flow_in = mm.addVars(batteries, time_periods, name=\"flow_in\")\n","flow_out = mm.addVars(batteries, time_periods, name=\"flow_out\")\n","grid = mm.addVars(time_periods, name=\"grid\")\n","state = mm.addVars(batteries, time_periods, ub = [capacity[b] for b in batteries for t in time_periods], name=\"state\")\n","gen = mm.addVars(time_periods, name=\"gen\")\n","zwitch = mm.addVars(batteries, time_periods, vtype=GRB.BINARY, name=\"zwitch\")\n","\n","# define a linear expression for total energy purchased from the grid\n","total_grid = grid.sum()\n","# define a linear expression for total cost\n","total_cost = gp.quicksum(avg_price.price[t]*grid[t] for t in time_periods)\n","\n","power_balance = mm.addConstrs((gp.quicksum(flow_out[b,t] - p_loss[b]*flow_in[b,t] for b in batteries) + gen[t] + grid[t] == total_demand[t]\n"," for t in time_periods), name=\"power_balance\")\n","initial_state = mm.addConstrs((state[b,0] == initial[b] + p_loss[b]*flow_in[b,0] - flow_out[b,0] for b in batteries), name=\"initial_state\")\n","subsequent_states = mm.addConstrs((state[b,t] == state[b,t-1] + p_loss[b]*flow_in[b,t] - flow_out[b,t] for b in batteries for t in time_periods if t >= 1), name=\"subsequent_states\")\n","solar_avail = mm.addConstrs((flow_in['Battery0',t] + flow_in['Battery1',t] + gen[t] <= solar_values[t] for t in time_periods), name = \"solar_avail\")\n","to_charge = mm.addConstrs((flow_in[b,t] <= 20*zwitch[b,t] for b in batteries for t in time_periods), name = \"to_charge\")\n","or_not_to_charge = mm.addConstrs((flow_out[b,t] <= 20*(1-zwitch[b,t]) for b in batteries for t in time_periods), name = \"or_not_to_charge\")\n","\n","mm.setObjective(total_cost, GRB.MINIMIZE)\n","mm.update()"]},{"cell_type":"markdown","metadata":{"id":"AXTW-_xBiopN"},"source":["#### Compare Four Scenarios\n","Here we'll outline four different scenarios that, outside of the base scenario, change the objective coefficients, variable bounds, and constraint RHS, respectively.\n","0. Base Scenario as above\n","> We'll set the number of scenarios and let the first be the base case.\n","\n"]},{"cell_type":"code","execution_count":89,"metadata":{"executionInfo":{"elapsed":231,"status":"ok","timestamp":1699601559348,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"BTzYqm6e6jX2"},"outputs":[],"source":["mm.NumScenarios=4\n","mm.Params.ScenarioNumber = 0\n","mm.ScenNName = 'Base model'"]},{"cell_type":"markdown","metadata":{"id":"dTUJg9_H6jf6"},"source":["1. (Mostly) higher costs\n",">Another column was added to the price dataframe where the original price was altered so high prices have become even higher and low prices have become even lower. The exact multiplier was random so the prices may seem to jump up and down a little more."]},{"cell_type":"code","execution_count":90,"metadata":{"executionInfo":{"elapsed":91,"status":"ok","timestamp":1699601560462,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"8YmIbGgr6Mwt"},"outputs":[],"source":["price2 = avg_price.price2\n","# change obj\n","mm.Params.ScenarioNumber = 1\n","mm.ScenNName = 'Increased price'\n","for t in time_periods:\n"," grid[t].ScenNObj = price2[t]"]},{"cell_type":"markdown","metadata":{"id":"hC44jopW6M9c"},"source":["2. Lower Battery Capacity\n","> Reduced battery capacities from 60 $→$ 48 and 80 $→$ 64.\n"]},{"cell_type":"code","execution_count":91,"metadata":{"executionInfo":{"elapsed":204,"status":"ok","timestamp":1699601561843,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"g-BQPRwY6NGT"},"outputs":[],"source":["capacity2 = {'Battery0': 48, 'Battery1': 64}\n","# change variable bound\n","mm.Params.ScenarioNumber = 2\n","mm.ScenNName = 'Low Battery'\n","for b in batteries:\n"," for t in time_periods:\n"," state[b,t].ScenNUB = capacity2[b]"]},{"cell_type":"markdown","metadata":{"id":"ndrbqFXO6NM0"},"source":["3. Higher Solar\n",">The solar prediction data frame has many other columns as standard output for a Prophet model (this was used to predict the solar availability). We'll take a convex combination of a lower estimae, point estimate, and upper estimate to create a new value for solar availability. The weights chosen will increase the forecast."]},{"cell_type":"code","execution_count":92,"metadata":{"executionInfo":{"elapsed":97,"status":"ok","timestamp":1699601563437,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"mO6tNo2j6NTD"},"outputs":[],"source":["solar_values2 = round(0.1*solar_values_read.yhat_lower +\n"," 0.6*solar_values_read.yhat +\n"," 0.3*solar_values_read.yhat_upper,3)\n","# just in case any forecasted value is negative, set it to 0\n","solar_values2[solar_values2 < 0] = 0\n","\n","# change constraint RHS\n","mm.Params.ScenarioNumber = 3\n","mm.ScenNName = 'High Solar'\n","for t in time_periods:\n"," solar_avail[t].ScenNRhs = solar_values2[t]"]},{"cell_type":"markdown","metadata":{"id":"VrpVF6XZ6NZl"},"source":["#### Querying Scenarios\n","We can now write an `.lp` file and optimize the scenarios. What do you see in the log that shows how Gurobi tries to solve these scenarios efficiently?"]},{"cell_type":"code","execution_count":93,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":438,"status":"ok","timestamp":1699601565611,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"VYfQt63wwU0R","outputId":"96499dcb-3f62-4898-f615-2bb9fc80cbe4"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x00b35a48\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 3e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","Solving a multi-scenario model with 4 scenarios...\n","\n","Found heuristic solution: objective 2722.4370000\n","Presolve removed 103 rows and 89 columns\n","Presolve time: 0.03s\n","Presolved: 1872 rows, 1889 columns, 5668 nonzeros\n","Found heuristic solution: objective 2713.0375100\n","Variable types: 1549 continuous, 340 integer (340 binary)\n","Found heuristic solution: objective 868.9568220\n","\n","Root relaxation: objective 4.448080e+02, 1364 iterations, 0.03 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 444.80800 0 5 868.95682 444.80800 48.8% - 0s\n","H 0 0 445.0857385 444.80800 0.06% - 0s\n","Scenario 3 has been solved (gap 0.0037%). 3/4 scenarios left.\n"," 0 0 445.06940 0 1 445.08574 445.06940 0.00% - 0s\n","H 0 0 445.0694028 445.06940 0.00% - 0s\n"," 0 0 578.28849 0 2 445.06940 578.28849 0.00% - 0s\n","\n","Optimal solution found at node 0 - now completing multiple scenarios...\n","\n"," Nodes | Current Node | Scenario Obj. Bounds | Work\n"," | | Worst |\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 2 579.16759 0 2 - 579.16759 - - 0s\n","Scenario 0 has been solved. 2/4 scenarios left.\n","Scenario 2 has been solved. 1/4 scenarios left.\n","Scenario 1 has been solved. 0/4 scenarios left.\n","\n","Cutting planes:\n"," Implied bound: 7\n"," MIR: 19\n"," Flow cover: 26\n"," Relax-and-lift: 1\n","\n","Explored 3 nodes (1835 simplex iterations) in 0.31 seconds (0.08 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 10: 445.069 445.069 445.086 ... 670.871\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 4.450694027691e+02, best bound 4.450694027691e+02, gap 0.0000%\n"]}],"source":["mm.write('ms.lp')\n","mm.optimize()\n","\n","#print obj value from each scenario\n","for s in range(m.NumScenarios):\n"," # Set the scenario number to query the information for this scenario\n"," mm.Params.ScenarioNumber = s\n"," print('\\nTotal cost for '+mm.ScenNName+' is '+'$'+str(round(mm.ScenNObjVal,2)))"]},{"cell_type":"markdown","metadata":{"id":"KM6o-kdx3YbO"},"source":["### Solution Pools\n","Dive in, the water's fine!\n","\n","A solution pool allows you to find multiple solutions that satisfy some criteria base off of search parameters set."]},{"cell_type":"code","execution_count":145,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":209,"status":"ok","timestamp":1699603570635,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"tYvUI5HeAVg9","outputId":"6c40ac97-549c-4bf8-8d9c-fe3c03043100"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x66bddde2\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Presolved: 1381 rows, 1741 columns, 4306 nonzeros\n","\n","Continuing optimization...\n","\n","\n","Explored 3964 nodes (8111 simplex iterations) in 0.02 seconds (0.00 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 500: 601.976 601.976 601.976 ... 602.035\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240991e+02, gap 0.0000%\n"]}],"source":["# Limit how many solutions to collect\n","mp.setParam(GRB.Param.PoolSolutions, 250)\n","# Limit the search space by setting a gap for the worst possible solution that will be accepted\n","mp.setParam(GRB.Param.PoolGap, 0.05)\n","# do a systematic search for the k-best solutions\n","mp.setParam(GRB.Param.PoolSearchMode, 2)\n","\n","mp.optimize()"]},{"cell_type":"markdown","metadata":{"id":"vOaMTUBrRRqW"},"source":["Now we can query all 250 solutions for variables of interest and look for anything interesting."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"DmSx15LsMs9Z"},"outputs":[],"source":["n = mp.SolCount\n","flow_250 = pd.DataFrame()\n","for i in range(n):\n"," mp.setParam(GRB.Param.SolutionNumber, i)\n"," for b in batteries:\n"," tmp = pd.DataFrame([[flow_in[b,t].Xn, flow_out[b,t], i] for t in time_periods], columns=[b+'_out',b+'_in','Scen'+str(i)])\n"," flow_250 = pd.concat([flow_250, tmp])"]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":245},"executionInfo":{"elapsed":149,"status":"error","timestamp":1699514549043,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"uyG2mtFOBHIl","outputId":"d254ff84-5808-4900-e498-57fe2d4c925a"},"outputs":[{"ename":"NameError","evalue":"ignored","output_type":"error","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)","\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0ms0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msol_level\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Battery0'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'orange'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0ms1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msol_level\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Battery1'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'blue'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Battery State (kWhr)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Time Period'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"]}],"source":["plt.figure(figsize=(12,5))\n","s0, = plt.plot(sol_level['Battery0'], c = 'orange')\n","s1, = plt.plot(sol_level['Battery1'], c = 'blue')\n","plt.ylabel('Battery State (kWhr)')\n","plt.xlabel('Time Period')\n","plt.legend([s0,s1],[\"Battery0\", \"Battery1\"])\n","plt.axhline(y=capacity['Battery0'], c='orange', linestyle='--', alpha = 0.5)\n","plt.axhline(y=capacity['Battery1'], c='blue', linestyle='--', alpha = 0.5)\n","print(f\"Periods at Battery0 at Full Capacity: {sum(sol_level['Battery0']==capacity['Battery0'])}\")\n","print(f\"Periods at Battery1 at Full Capacity: {sum(sol_level['Battery1']==capacity['Battery1'])}\");"]},{"cell_type":"markdown","metadata":{"id":"6RFo7kk39cil"},"source":[]},{"cell_type":"markdown","metadata":{"id":"iPu22TkK9cRE"},"source":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Yf6YU88f9buP"},"outputs":[],"source":[]}],"metadata":{"colab":{"authorship_tag":"ABX9TyN80moD6M6voTt1zw7joXtQ","provenance":[]},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":0} +{"cells":[{"cell_type":"markdown","metadata":{"id":"ezGOS9cRadzN"},"source":["# Intermediate Modeling - Part 2: Back for More Sun\n","Just as we expanded on the first modeling example from Opti 101 for the first example of Opti 201, we'll do the same here. This example is motivated by [IEEE's Predict+Optimize Technical Challenge](https://ieee-dataport.org/competitions/ieee-cis-technical-challenge-predictoptimize-renewable-energy-scheduling) in 2021.\n","\n","In this example we are trying to minimize the amount of electricity purchased from the grid to supply a building with electricity over a 6-day period while meeting predicted demand. The building is outfitted with solar panels and two batteries to store and provide energy when appropriate. Decisions made at each 30-minute interval are:\n","1. Buy electricity from the grid\n","2. Charge/discharge the one battery\n","3. Charge/discharge the other battery\n","\n","The output of the optimization model is a charge/discharge schedule for the batteries as well as the time periods where electricity will be purchased. The demand is forecasted using [Prophet](https://facebook.github.io/prophet/). The notebook for this part of the problem can be found [here](https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/optimization101/Modeling_Session_2/energy_storage_ML.ipynb) If you want to have the first notebook handy, you can find it [here](https://colab.research.google.com/github/Gurobi/modeling-examples/blob/master/optimization101/Modeling_Session_2/completed_modeling2.ipynb).\n","\n","This model is a lot more complicated that the first from Opti 101, so we can copy the formulation and code in a concise way and highlight the parts relevant to the features discussed in the second modeling session video.\n","\n","In subsequent sections we'll see how to adapt this model to incorporate\n","- Multi-objective models,\n","- Multi-scenario models, and\n","- Solutions pools"]},{"cell_type":"markdown","metadata":{"id":"6uEfx_64u0jY"},"source":["## Opti 101 Model"]},{"cell_type":"markdown","metadata":{"id":"E5lo-xCHu4Jk"},"source":["### Formulation"]},{"cell_type":"markdown","metadata":{"id":"CXPsSNehu8Nu"},"source":["#### Sets and Input Parameters\n","\n","Batteries: $b \\in B = \\{\\texttt{Battery0, Battery1}\\}$\n","\n","Time periods: $t \\in T = \\{0,1,...,179\\}$ for the first Monday through Saturday of October, 2020\n","\n","$c_{b}$: capacity of battery $b \\in B \\quad\\quad \\texttt{capacity[b]}$\n","\n","$p_{b}$: loss of energy (as a percentage) during transfer into battery $b\\in B \\quad\\quad \\texttt{p}\\_\\texttt{loss[b]}$\n","\n","$q_{b}$: quantity of initial energy in battery $b \\in B \\quad\\quad \\texttt{initial[b]}$\n","\n","$solar_{t}$: solar power generation of the panel for time period $t \\in T \\quad\\quad \\texttt{solar}\\_\\texttt{values[b]}$\n","\n","$d_{t}$: total building and class energy demand for time period $t\\in T \\quad\\quad \\texttt{total}\\_\\texttt{deamnd[t]}$"]},{"cell_type":"markdown","metadata":{"id":"bOQnhuPkvIcg"},"source":["#### Decision Variables\n","Let $f^{in}_{b,t}$ be the amount each battery, $b$, `charges` at time period, $t$, $\\forall b\\in B, t\\in T$. $\\quad\\quad \\texttt{flow}\\_\\texttt{in[b,t]}$\n","\n","Let $f^{out}_{b,t}$ be the amount each battery `discharges`, similarly. $\\quad\\quad \\texttt{flow}\\_\\texttt{out[b,t]}$\n","\n","The max amount that each battery can charge or discharge in a single period to be 20 kW.\n","\n","$s_{b,t}$ is the current amount of energy in battery, $b$, at the end of time period, $t, \\forall b\\in B, t\\in T$. $\\quad\\texttt{state[b,t]}$\n","\n","$gen_{t}$ is the amount of available solar energy that is used in period $t$, $\\forall t \\in T$. $\\quad\\texttt{gen[t]}$\n","\n","$grid_{t}$ : This variable indicates the amount of energy purchased from the grid at time period, $t$, $\\forall t \\in T$. $\\quad\\texttt{grid[t]}$\n","\n","We are not considering the ability to \"sell back\" electricity, so $grid_t \\ge 0$."]},{"cell_type":"markdown","metadata":{"id":"jZ5Zbv6PF6xU"},"source":["#### Constraints\n","**Power balance:**\n","$$\n","\\begin{align}\n","\\sum_b(f^{out}_{b,t}-p_bf^{in}_{b,t}) + gen_t + grid_t = d_t \\quad \\forall t \\in T\n","\\end{align}\n","$$\n","\n","**Battery state:**\n","\n","First time period's charge/discharge decision.\n","\n","\\begin{equation}\n","s_{b,0} = q_b + p_bf^{in}_{b,0} - f^{out}_{b,0}\n","\\end{equation}\n","\n","For each time period after (i.e. $t\\ge1$), the state of a battery is found by:\n","\n","\\begin{equation}\n","s_{b,t} = s_{b,t-1} + p_bf^{in}_{b,t} - f^{out}_{b,t}, t \\ge 1\n","\\end{equation}\n","\n","**Solar availability:**\n","$$\n","\\begin{equation}\n","f^{in}_{\\texttt{Battery0},t} + f^{in}_{\\texttt{Battery1},t} + gen_t \\le solar_t, \\quad \\forall t \\in T\n","\\end{equation}\n","$$\n","\n","**Charge/discharge:**\n","We cannot charge and discharge the same battery in the same time period.\n","Define $z_{b,t} \\in \\{0,1\\}, \\forall b \\in B, t \\in T$.\n","$$\n","\\begin{align*}\n","f^{in}_{b,t} &\\leq 20*z_{b,t} &\\forall b \\in B, t \\in T \\\\\n","f^{out}_{b,t} &\\leq 20*(1-z_{b,t}) &\\forall b \\in B, t \\in T\n","\\end{align*}\n","$$\n","Hmm, this set of constraints looks a little familiar, huh?\n","\n","**Battery capacity:**\n","$$\n","\\begin{align*}\n","s_{b,t} \\le c_b, \\quad \\forall b \\in B, t\\in T\n","\\end{align*}\n","$$"]},{"cell_type":"markdown","metadata":{"id":"7j--9VCkuqCt"},"source":["#### Objective Function\n","The objective is to minimize the total amount of energy purchased from the grid over all time periods.\n","\n","\\begin{equation}\n","{\\rm min} \\space \\sum_{t} g_{t}\n","\\end{equation}"]},{"cell_type":"markdown","metadata":{"id":"dpyvjhp4oT4N"},"source":["### Gurobipy Code"]},{"cell_type":"markdown","metadata":{"id":"8EBFFKNPvft2"},"source":["#### Sets and Input parameters"]},{"cell_type":"code","execution_count":1,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":15980,"status":"ok","timestamp":1699599941524,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"iPw18bDQF9Lc","outputId":"4f3ddf9a-403a-4536-d766-0086b80c0423"},"outputs":[{"name":"stdout","output_type":"stream","text":["Collecting gurobipy\n"," Downloading gurobipy-10.0.3-cp310-cp310-manylinux2014_x86_64.whl (12.7 MB)\n","\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m12.7/12.7 MB\u001b[0m \u001b[31m20.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n","\u001b[?25hInstalling collected packages: gurobipy\n","Successfully installed gurobipy-10.0.3\n","Total Solar Generation: 4939.172 \n","Total Demand: 5250.6\n"]}],"source":["# install/import packages\n","%pip install gurobipy\n","import pandas as pd\n","import matplotlib.pyplot as plt\n","import gurobipy as gp\n","from gurobipy import GRB\n","\n","# read in solar forecast data\n","path = 'https://raw.githubusercontent.com/Gurobi/modeling-examples/master/optimization101/Modeling_Session_2/'\n","solar_values_read = pd.read_csv(path+'pred_solar_values.csv')\n","solar_values = round(solar_values_read.yhat,3)\n","solar_values.reset_index(drop = True, inplace = True)\n","\n","# read in demaind data, where total demand is a fixed building demand and\n","# an estimated demand based on a proposed schedule for the building\n","schedule = pd.read_csv(path+'schedule_demand.csv')\n","avg_building = pd.read_csv(path+'building_demand.csv')\n","total_demand = schedule.sched_demand + avg_building.build_demand\n","print(f\"Total Solar Generation: {solar_values.sum()} \\nTotal Demand: {total_demand.sum()}\")\n","\n","# create data for batteries, including capity, efficency and initial charges\n","# also define time periods,\n","batteries = [\"Battery0\", \"Battery1\"]\n","capacity = {\"Battery0\": 60, \"Battery1\": 80} # in Kw\n","p_loss = {\"Battery0\": 0.95, \"Battery1\": 0.9} # proportion\n","initial = {\"Battery0\": 0, \"Battery1\": 0} # in kW\n","time_periods = range(len(solar_values_read))"]},{"cell_type":"markdown","metadata":{"id":"20u3pAiAulYD"},"source":["#### Decison Variables"]},{"cell_type":"code","execution_count":113,"metadata":{"executionInfo":{"elapsed":101,"status":"ok","timestamp":1699602343575,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"mpi8sWtKqEMw"},"outputs":[],"source":["m = gp.Model() #this defines the model that we'll add to as we finish the formulation\n","\n","flow_in = m.addVars(batteries, time_periods, name=\"flow_in\")\n","flow_out = m.addVars(batteries, time_periods, name=\"flow_out\")\n","\n","grid = m.addVars(time_periods, name=\"grid\")\n","state = m.addVars(batteries, time_periods, ub = [capacity[b] for b in batteries for t in time_periods], name=\"state\")\n","gen = m.addVars(time_periods, name=\"gen\")\n","\n","zwitch = m.addVars(batteries, time_periods, vtype=GRB.BINARY, name=\"zwitch\")"]},{"cell_type":"markdown","metadata":{"id":"pJIGt-hLqgWi"},"source":["#### Constraints\n","Note than the **capacity constraint** has been taken care of by putting an `ub` on the `state` decision variable."]},{"cell_type":"code","execution_count":114,"metadata":{"executionInfo":{"elapsed":138,"status":"ok","timestamp":1699602344717,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"7vnVisjTqiHL"},"outputs":[],"source":["# Power balance\n","m.addConstrs((gp.quicksum(flow_out[b,t] - p_loss[b]*flow_in[b,t] for b in batteries) + gen[t] + grid[t] == total_demand[t]\n"," for t in time_periods), name=\"power_balance\");\n","# Battery state\n","m.addConstrs((state[b,0] == initial[b] + p_loss[b]*flow_in[b,0] - flow_out[b,0] for b in batteries), name=\"initial_state\")\n","m.addConstrs((state[b,t] == state[b,t-1] + p_loss[b]*flow_in[b,t] - flow_out[b,t] for b in batteries for t in time_periods if t >= 1), name=\"subsequent_states\");\n","\n","#Solar availability\n","m.addConstrs((flow_in['Battery0',t] + flow_in['Battery1',t] + gen[t] <= solar_values[t] for t in time_periods), name = \"solar_avail\");\n","\n","#Charge/discharge\n","m.addConstrs((flow_in[b,t] <= 20*zwitch[b,t] for b in batteries for t in time_periods), name = \"to_charge\")\n","m.addConstrs((flow_out[b,t] <= 20*(1-zwitch[b,t]) for b in batteries for t in time_periods), name = \"or_not_to_charge\");\n"]},{"cell_type":"markdown","metadata":{"id":"AqsU88pnqpiQ"},"source":["#### Objective Function\n","The original objective is to minimize the amount of electricity purchased from the grid. There was also a section where we considered a different objective -- given expected electricity prices for each period, minimize the total cost. So we can compare the cost and amount of electricity purchased across scenarios, let's start to capture these values."]},{"cell_type":"code","execution_count":115,"metadata":{"executionInfo":{"elapsed":253,"status":"ok","timestamp":1699602346323,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"q-yYYpvEq0OJ"},"outputs":[],"source":["# read in estimated price of electricity for each time period\n","avg_price = pd.read_csv(path+'expected_price.csv')\n","price = avg_price.price\n","\n","# define a linear expression for total energy purchased from the grid\n","total_grid = grid.sum()\n","# define a linear expression for total cost\n","total_cost = gp.quicksum(avg_price.price[t]*grid[t] for t in time_periods)\n","\n","# set the objective to the original of min grid purchase\n","m.setObjective(total_grid, GRB.MINIMIZE)"]},{"cell_type":"markdown","metadata":{"id":"uYDpiO-2v7dD"},"source":["### Solve the Model"]},{"cell_type":"code","execution_count":116,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":686},"executionInfo":{"elapsed":190,"status":"ok","timestamp":1699602348004,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"NbNVPTBZwC8E","outputId":"6ed9a448-59d5-46ba-b589-385f29380ce6"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0xaa31bb6a\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [1e+00, 1e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Found heuristic solution: objective 5250.6000000\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Found heuristic solution: objective 4986.0960000\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 1.281678e+03, 594 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n","* 0 0 0 1281.6776848 1281.67768 0.00% - 0s\n","\n","Explored 1 nodes (594 simplex iterations) in 0.08 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 3: 1281.68 4986.1 5250.6 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.281677684779e+03, best bound 1.281677684779e+03, gap 0.0000%\n"]},{"data":{"text/html":["\n","
\n","
\n","\n","\n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n","
CostGridAmount
0803.181281.68
\n","
\n","
\n","\n","
\n"," \n","\n"," \n","\n"," \n","
\n","\n","
\n","
\n"],"text/plain":[" Cost GridAmount\n","0 803.18 1281.68"]},"execution_count":116,"metadata":{},"output_type":"execute_result"}],"source":["m.optimize()\n","\n","results = pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]],\n"," columns = ['Cost','GridAmount']\n"," )\n","\n","results"]},{"cell_type":"markdown","metadata":{"id":"TTBd1ssUpRxX"},"source":["The original objective is to minimize the amount of electricity purchased from the grid. There was also a section where we considered a different objective -- given expected electricity prices for each period, minimize the total cost. We'll read in the data and quickly update the objective and resolve."]},{"cell_type":"code","execution_count":117,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":228,"status":"ok","timestamp":1699602350218,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"vBeoWLq6pcpO","outputId":"ffdedc1b-0b30-4460-c837-639f5f8862fa"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x66bddde2\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","MIP start from previous solve produced solution with objective 617.569 (0.04s)\n","MIP start from previous solve produced solution with objective 617.569 (0.04s)\n","Loaded MIP start from previous solve with objective 617.569\n","\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.01s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 6.019761e+02, 754 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 601.97612 0 1 617.56942 601.97612 2.52% - 0s\n","H 0 0 601.9761241 601.97612 0.00% - 0s\n"," 0 0 601.97612 0 1 601.97612 601.97612 0.00% - 0s\n","\n","Explored 1 nodes (754 simplex iterations) in 0.11 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 2: 601.976 617.569 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240990e+02, gap 0.0000%\n"]}],"source":["m.setObjective(total_cost, GRB.MINIMIZE)\n","m.optimize()\n","\n","results = pd.concat([results,\n"," pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]], columns = ['Cost','GridAmount'])\n"," ],\n"," ignore_index=True\n"," )\n","\n","results\n","#create a copy for later\n","mp = m.copy()"]},{"cell_type":"markdown","metadata":{"id":"pU48SBLRxjun"},"source":["## Let's Find Some Hidden Gems"]},{"cell_type":"markdown","metadata":{"id":"_7L8k0t-0cgI"},"source":["### Handling Multiple Objectives\n","\n","The above results show two different scenarios, one that purely considers grid purchases, and the other purely prioritizes cost. Let's see how to use Gurobi's [multi-objective](https://www.gurobi.com/documentation/current/refman/working_with_multiple_obje.html) functionality to combine the two objectives by first implementing a `weighed` setup."]},{"cell_type":"markdown","metadata":{"id":"AOPH_t4eDAiV"},"source":["#### Weighted Objectives\n","In using a weighted objective, we need to obviously decide the `weights`. This can be a little tricky given units and scale."]},{"cell_type":"code","execution_count":82,"metadata":{"executionInfo":{"elapsed":119,"status":"ok","timestamp":1699601541125,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"H4ntqJhexqZY"},"outputs":[],"source":["m.setObjectiveN(total_cost, index=0, weight = 1, name=\"cost\")\n","m.setObjectiveN(total_grid, index=1, weight = 10, name= \"grid\")"]},{"cell_type":"markdown","metadata":{"id":"G6muNC5D0cUT"},"source":["Now we can run the optimization with our weighted objective."]},{"cell_type":"code","execution_count":83,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":1000},"executionInfo":{"elapsed":473,"status":"ok","timestamp":1699601542850,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"frdBl5DM0Zg4","outputId":"312fbaf3-fa93-4efe-c68a-8382b8cfd049"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x9cf62986\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: starting optimization with 2 objectives (1 combined) ...\n","---------------------------------------------------------------------------\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 1 (weighted) ...\n","---------------------------------------------------------------------------\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0xf11f89cc\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [1e+01, 1e+01]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Found heuristic solution: objective 55228.437000\n","Presolve removed 84 rows and 265 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1535 columns, 4109 nonzeros\n","Found heuristic solution: objective 52452.392480\n","Variable types: 1199 continuous, 336 integer (336 binary)\n","\n","Root relaxation: objective 1.343869e+04, 703 iterations, 0.04 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 13438.6902 0 2 52452.3925 13438.6902 74.4% - 0s\n","H 0 0 13438.690179 13438.6902 0.00% - 0s\n"," 0 0 13438.6902 0 2 13438.6902 13438.6902 0.00% - 0s\n","\n","Explored 1 nodes (703 simplex iterations) in 0.14 seconds (0.03 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 3: 13438.7 52452.4 55228.4 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.343869017947e+04, best bound 1.343869017947e+04, gap 0.0000%\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: solved in 0.16 seconds (0.03 work units), solution count 3\n","\n"]},{"data":{"text/html":["\n","
\n","
\n","\n","\n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n","
CostGridAmount
0803.181281.68
1601.981317.69
2616.701282.20
\n","
\n","
\n","\n","
\n"," \n","\n"," \n","\n"," \n","
\n","\n","\n","
\n"," \n","\n","\n","\n"," \n","
\n","
\n","
\n"],"text/plain":[" Cost GridAmount\n","0 803.18 1281.68\n","1 601.98 1317.69\n","2 616.70 1282.20"]},"execution_count":83,"metadata":{},"output_type":"execute_result"}],"source":["m.ModelSense = GRB.MINIMIZE\n","m.optimize()\n","\n","results = pd.concat([results,\n"," pd.DataFrame([[round(v,2) for v in [total_cost.getValue(),total_grid.getValue()]]], columns = ['Cost','GridAmount'])\n"," ],\n"," ignore_index=True\n"," )\n","\n","results"]},{"cell_type":"markdown","metadata":{"id":"f1HiiJcZBLws"},"source":["Using this multi-objective approach, we have been able to significantly reduce the cost from the first solution* while adding on a very small amount of grid purchase.\n","\n","*Note that since we didn't specify anything in the first model that pertained to price there may be a solution that uses the same amout of purchased energy at a different cost.\n"]},{"cell_type":"markdown","metadata":{"id":"UcoCaW7zDH-H"},"source":["#### Hierarchical Objectives\n","This kind of multi-objective approach gives priority to a given objective in a given order, and for subsequent objectives, constraints are added to make sure the expressions of prior objectives to not degrade beyond given values.\n","\n","In the `hierarchical` approach the two main arguments are `abstol` and `reltol`. The former allows the first objective to degrade by a *raw* amount, while the latter is a *percentage*.\n","\n","For example, we can set first optimize for minimal total cost, and then allow that objective value to *degrade* by $50 or by a factor of 0.1 (i.e. increase by 10%) when we then minimize for the lowest amount of electricity purchased.\n","\n","An important concept in energy storage in maintaining battery health. One way this may be handled is to limit a battery's [depth of discharge](https://en.wikipedia.org/wiki/Depth_of_discharge), which can be expressed as the percentage of a battery's capacity that has been discharged. The way we decide to do this is by limiting the number of time periods our batteries are discharged below a certain depth.\n","\n","To do this we'll need binary variables to indicate when this occurs.\n","\n","Let $v_{b,t} = 1$ when battery $b$ is below a depth of discharge α (say 0.3) and 0 otherwise. We want to model\n","\n","\"If the battery gets below 30% it's capacity, then indicate depth of discharge violation.\"\n","\n","Let's replace some what comes after \"If\" and \"then\" with decision variables and (in)equalities. As a reminder $s_{b,t}$ is the current charge in the battery and $c_b$ is its capacity.\n","\n","$$\n","If \\space s_{b,t} < 0.3 \\times c_b, \\space then \\space v_{b,t} = 1\n","$$\n","\n","We can easily model this by taking the contrapositive of the conditional statement above and using indicator constraints.\n","\n","$$\n","If \\space v_{b,t} = 0, \\space then \\space s_{b,t} \\ge 0.3 \\times c_b,\n","$$\n","\n","````\n","To keep the problem within the 2,000 variable limit that comes with the `pip` install license,\n","this will only be implemented for `Battery0`.\n","````"]},{"cell_type":"code","execution_count":84,"metadata":{"executionInfo":{"elapsed":201,"status":"ok","timestamp":1699601547810,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"XgMclZuoAq3q"},"outputs":[],"source":["# Reminder that s_{b,t} is defined by state[b,t] in the code\n","v = m.addVars(time_periods, vtype=GRB.BINARY, name = 'v')\n","\n","# define a linear expression for total depth count\n","total_depth_count = v.sum()\n","\n","m.addConstrs((((v[t] == 0) >> (state['Battery0',t] >= 0.3*capacity['Battery0'])) for b in batteries for t in time_periods), name = \"discharge_depth\")\n","m.update()"]},{"cell_type":"markdown","metadata":{"id":"bXnwpJx1Ex85"},"source":["Now we'll dive into a three-piece objective."]},{"cell_type":"code","execution_count":85,"metadata":{"executionInfo":{"elapsed":97,"status":"ok","timestamp":1699601549020,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"jdVLoqVfHi4C"},"outputs":[],"source":["#the higher priority number is used first as an objective\n","m.setObjectiveN(total_cost, index=0, priority=2, reltol = 0.05, name=\"cost\")\n","m.setObjectiveN(total_depth_count, index=1, priority=1, reltol = 0.2, name= \"depth\")\n","m.setObjectiveN(total_grid, index=2, priority=0, name= \"grid\")"]},{"cell_type":"code","execution_count":86,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":1164,"status":"ok","timestamp":1699601551154,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"8EddQ72UVif6","outputId":"0419d2ca-b72d-4817-c5f8-281eacb3f1dc"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1980 columns and 4498 nonzeros\n","Model fingerprint: 0x4086b7a9\n","Model has 360 general constraints\n","Variable types: 1440 continuous, 540 integer (540 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n"," GenCon rhs range [2e+01, 2e+01]\n"," GenCon coe range [1e+00, 1e+00]\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: starting optimization with 3 objectives ... \n","---------------------------------------------------------------------------\n","\n","Multi-objectives: applying initial presolve ...\n","---------------------------------------------------------------------------\n","\n","Presolve added 305 rows and 302 columns\n","Presolve time: 0.02s\n","Presolved: 1745 rows and 2282 columns\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 1 (cost) ...\n","---------------------------------------------------------------------------\n","\n","Presolve removed 389 rows and 671 columns\n","Presolve time: 0.02s\n","Presolved: 1356 rows, 1611 columns, 4160 nonzeros\n","Variable types: 1275 continuous, 336 integer (336 binary)\n","Found heuristic solution: objective 1021.0681360\n","\n","Root relaxation: objective 6.019761e+02, 697 iterations, 0.02 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 601.97612 0 3 1021.06814 601.97612 41.0% - 0s\n","H 0 0 601.9761241 601.97612 0.00% - 0s\n"," 0 0 601.97612 0 3 601.97612 601.97612 0.00% - 0s\n","\n","Explored 1 nodes (697 simplex iterations) in 0.15 seconds (0.04 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 2: 601.976 1021.07 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240991e+02, gap 0.0000%\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 2 (depth) ...\n","---------------------------------------------------------------------------\n","\n","\n","Loaded user MIP start with objective 180\n","\n","Presolve removed 218 rows and 398 columns\n","Presolve time: 0.04s\n","Presolved: 1528 rows, 1884 columns, 4802 nonzeros\n","Variable types: 1373 continuous, 511 integer (511 binary)\n","Found heuristic solution: objective 66.0000000\n","\n","Root relaxation: objective 9.475247e+00, 1721 iterations, 0.05 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 9.47525 0 8 66.00000 9.47525 85.6% - 0s\n","H 0 0 14.0000000 9.47525 32.3% - 0s\n"," 0 0 10.03868 0 8 14.00000 10.03868 28.3% - 0s\n"," 0 0 10.26646 0 7 14.00000 10.26646 26.7% - 0s\n"," 0 0 10.30054 0 9 14.00000 10.30054 26.4% - 0s\n"," 0 0 10.31001 0 10 14.00000 10.31001 26.4% - 0s\n"," 0 0 10.31001 0 10 14.00000 10.31001 26.4% - 0s\n"," 0 0 10.41407 0 10 14.00000 10.41407 25.6% - 0s\n"," 0 0 10.45727 0 10 14.00000 10.45727 25.3% - 0s\n"," 0 0 10.45727 0 10 14.00000 10.45727 25.3% - 0s\n","H 0 0 13.0000000 10.45727 19.6% - 0s\n","H 0 0 12.0000000 10.45727 12.9% - 0s\n"," 0 2 10.71604 0 4 12.00000 10.71604 10.7% - 0s\n","\n","Cutting planes:\n"," Gomory: 2\n"," MIR: 12\n"," Flow cover: 9\n","\n","Explored 5 nodes (1802 simplex iterations) in 0.54 seconds (0.14 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 5: 12 13 14 ... 180\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.200000000000e+01, best bound 1.200000000000e+01, gap 0.0000%\n","---------------------------------------------------------------------------\n","\n","Multi-objectives: optimize objective 3 (grid) ...\n","---------------------------------------------------------------------------\n","\n","\n","Loaded user MIP start with objective 1343.89\n","\n","Presolve removed 218 rows and 398 columns\n","Presolve time: 0.04s\n","Presolved: 1529 rows, 1884 columns, 4977 nonzeros\n","Variable types: 1373 continuous, 511 integer (511 binary)\n","\n","Root relaxation: objective 1.292876e+03, 1450 iterations, 0.03 seconds (0.01 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 1292.87646 0 13 1343.88759 1292.87646 3.80% - 0s\n"," 0 0 1302.07834 0 10 1343.88759 1302.07834 3.11% - 0s\n"," 0 0 1302.29363 0 12 1343.88759 1302.29363 3.10% - 0s\n"," 0 0 1302.91971 0 11 1343.88759 1302.91971 3.05% - 0s\n"," 0 0 1302.91971 0 12 1343.88759 1302.91971 3.05% - 0s\n"," 0 0 1303.20394 0 9 1343.88759 1303.20394 3.03% - 0s\n"," 0 0 1303.22949 0 9 1343.88759 1303.22949 3.03% - 0s\n"," 0 0 1303.23331 0 10 1343.88759 1303.23331 3.03% - 0s\n"," 0 0 1303.60134 0 12 1343.88759 1303.60134 3.00% - 0s\n"," 0 0 1303.60694 0 12 1343.88759 1303.60694 3.00% - 0s\n"," 0 0 1303.63927 0 12 1343.88759 1303.63927 2.99% - 0s\n"," 0 0 1303.63927 0 12 1343.88759 1303.63927 2.99% - 0s\n"," 0 2 1304.18118 0 12 1343.88759 1304.18118 2.95% - 0s\n","* 11 9 8 1313.8029332 1307.01963 0.52% 12.6 0s\n","* 43 7 8 1313.7843075 1309.46527 0.33% 13.1 1s\n","* 45 0 6 1309.5801997 1309.58020 0.00% 12.8 1s\n","\n","Cutting planes:\n"," Gomory: 4\n"," Implied bound: 2\n"," MIR: 12\n"," Flow cover: 13\n","\n","Explored 46 nodes (2108 simplex iterations) in 1.02 seconds (0.28 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 4: 1309.58 1313.78 1313.8 1343.89 \n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 1.309580199681e+03, best bound 1.309580199681e+03, gap 0.0000%\n","\n","---------------------------------------------------------------------------\n","Multi-objectives: solved in 1.04 seconds (0.28 work units), solution count 9\n","\n"]}],"source":["m.optimize()"]},{"cell_type":"markdown","metadata":{"id":"TIMY-FXjn4LA"},"source":["We defined linear expressions for each part of the objective so we can use `getValue` to get their respective results. You can also query objectives as shown here."]},{"cell_type":"code","execution_count":87,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":89,"status":"ok","timestamp":1699601553769,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"jC4X0Bab_B5H","outputId":"ba1e01cd-ca50-4949-9274-0546ef669592"},"outputs":[{"name":"stdout","output_type":"stream","text":[" 632.14 14.0 1309.58"]}],"source":["for o in range(m.NumObj):\n"," # set which objective we will query\n"," m.params.ObjNumber = o\n"," # query the o-th objective value\n"," print(' ',round(m.ObjNVal,2), end='')"]},{"cell_type":"markdown","metadata":{"id":"nT1gJ5w4YJO4"},"source":["Play around with the priorities and tolerances to see how those change the overall solution."]},{"cell_type":"markdown","metadata":{"id":"fSw8RpaTY2CT"},"source":["### Handling Multiple Scenarios\n","In this section we will start with out **base model** with the objective being to minimize the total electricity purchase cost. To make is easier to edit constraints we can store."]},{"cell_type":"code","execution_count":88,"metadata":{"executionInfo":{"elapsed":283,"status":"ok","timestamp":1699601555555,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"XGsWkqaUnTbj"},"outputs":[],"source":["mm = gp.Model('mulit_scen') #this defines the model that we'll add to as we finish the formulation\n","flow_in = mm.addVars(batteries, time_periods, name=\"flow_in\")\n","flow_out = mm.addVars(batteries, time_periods, name=\"flow_out\")\n","grid = mm.addVars(time_periods, name=\"grid\")\n","state = mm.addVars(batteries, time_periods, ub = [capacity[b] for b in batteries for t in time_periods], name=\"state\")\n","gen = mm.addVars(time_periods, name=\"gen\")\n","zwitch = mm.addVars(batteries, time_periods, vtype=GRB.BINARY, name=\"zwitch\")\n","\n","# define a linear expression for total energy purchased from the grid\n","total_grid = grid.sum()\n","# define a linear expression for total cost\n","total_cost = gp.quicksum(avg_price.price[t]*grid[t] for t in time_periods)\n","\n","power_balance = mm.addConstrs((gp.quicksum(flow_out[b,t] - p_loss[b]*flow_in[b,t] for b in batteries) + gen[t] + grid[t] == total_demand[t]\n"," for t in time_periods), name=\"power_balance\")\n","initial_state = mm.addConstrs((state[b,0] == initial[b] + p_loss[b]*flow_in[b,0] - flow_out[b,0] for b in batteries), name=\"initial_state\")\n","subsequent_states = mm.addConstrs((state[b,t] == state[b,t-1] + p_loss[b]*flow_in[b,t] - flow_out[b,t] for b in batteries for t in time_periods if t >= 1), name=\"subsequent_states\")\n","solar_avail = mm.addConstrs((flow_in['Battery0',t] + flow_in['Battery1',t] + gen[t] <= solar_values[t] for t in time_periods), name = \"solar_avail\")\n","to_charge = mm.addConstrs((flow_in[b,t] <= 20*zwitch[b,t] for b in batteries for t in time_periods), name = \"to_charge\")\n","or_not_to_charge = mm.addConstrs((flow_out[b,t] <= 20*(1-zwitch[b,t]) for b in batteries for t in time_periods), name = \"or_not_to_charge\")\n","\n","mm.setObjective(total_cost, GRB.MINIMIZE)\n","mm.update()"]},{"cell_type":"markdown","metadata":{"id":"AXTW-_xBiopN"},"source":["#### Compare Four Scenarios\n","Here we'll outline four different scenarios that, outside of the base scenario, change the objective coefficients, variable bounds, and constraint RHS, respectively.\n","0. Base Scenario as above\n","> We'll set the number of scenarios and let the first be the base case.\n","\n"]},{"cell_type":"code","execution_count":89,"metadata":{"executionInfo":{"elapsed":231,"status":"ok","timestamp":1699601559348,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"BTzYqm6e6jX2"},"outputs":[],"source":["mm.NumScenarios=4\n","mm.Params.ScenarioNumber = 0\n","mm.ScenNName = 'Base model'"]},{"cell_type":"markdown","metadata":{"id":"dTUJg9_H6jf6"},"source":["1. (Mostly) higher costs\n",">Another column was added to the price dataframe where the original price was altered so high prices have become even higher and low prices have become even lower. The exact multiplier was random so the prices may seem to jump up and down a little more."]},{"cell_type":"code","execution_count":90,"metadata":{"executionInfo":{"elapsed":91,"status":"ok","timestamp":1699601560462,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"8YmIbGgr6Mwt"},"outputs":[],"source":["price2 = avg_price.price2\n","# change obj\n","mm.Params.ScenarioNumber = 1\n","mm.ScenNName = 'Increased price'\n","for t in time_periods:\n"," grid[t].ScenNObj = price2[t]"]},{"cell_type":"markdown","metadata":{"id":"hC44jopW6M9c"},"source":["2. Lower Battery Capacity\n","> Reduced battery capacities from 60 $→$ 48 and 80 $→$ 64.\n"]},{"cell_type":"code","execution_count":91,"metadata":{"executionInfo":{"elapsed":204,"status":"ok","timestamp":1699601561843,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"g-BQPRwY6NGT"},"outputs":[],"source":["capacity2 = {'Battery0': 48, 'Battery1': 64}\n","# change variable bound\n","mm.Params.ScenarioNumber = 2\n","mm.ScenNName = 'Low Battery'\n","for b in batteries:\n"," for t in time_periods:\n"," state[b,t].ScenNUB = capacity2[b]"]},{"cell_type":"markdown","metadata":{"id":"ndrbqFXO6NM0"},"source":["3. Higher Solar\n",">The solar prediction data frame has many other columns as standard output for a Prophet model (this was used to predict the solar availability). We'll take a convex combination of a lower estimae, point estimate, and upper estimate to create a new value for solar availability. The weights chosen will increase the forecast."]},{"cell_type":"code","execution_count":92,"metadata":{"executionInfo":{"elapsed":97,"status":"ok","timestamp":1699601563437,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"mO6tNo2j6NTD"},"outputs":[],"source":["solar_values2 = round(0.1*solar_values_read.yhat_lower +\n"," 0.6*solar_values_read.yhat +\n"," 0.3*solar_values_read.yhat_upper,3)\n","# just in case any forecasted value is negative, set it to 0\n","solar_values2[solar_values2 < 0] = 0\n","\n","# change constraint RHS\n","mm.Params.ScenarioNumber = 3\n","mm.ScenNName = 'High Solar'\n","for t in time_periods:\n"," solar_avail[t].ScenNRhs = solar_values2[t]"]},{"cell_type":"markdown","metadata":{"id":"VrpVF6XZ6NZl"},"source":["#### Querying Scenarios\n","We can now write an `.lp` file and optimize the scenarios. What do you see in the log that shows how Gurobi tries to solve these scenarios efficiently?"]},{"cell_type":"code","execution_count":93,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":438,"status":"ok","timestamp":1699601565611,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"VYfQt63wwU0R","outputId":"96499dcb-3f62-4898-f615-2bb9fc80cbe4"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x00b35a48\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 3e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","\n","Solving a multi-scenario model with 4 scenarios...\n","\n","Found heuristic solution: objective 2722.4370000\n","Presolve removed 103 rows and 89 columns\n","Presolve time: 0.03s\n","Presolved: 1872 rows, 1889 columns, 5668 nonzeros\n","Found heuristic solution: objective 2713.0375100\n","Variable types: 1549 continuous, 340 integer (340 binary)\n","Found heuristic solution: objective 868.9568220\n","\n","Root relaxation: objective 4.448080e+02, 1364 iterations, 0.03 seconds (0.02 work units)\n","\n"," Nodes | Current Node | Objective Bounds | Work\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 0 444.80800 0 5 868.95682 444.80800 48.8% - 0s\n","H 0 0 445.0857385 444.80800 0.06% - 0s\n","Scenario 3 has been solved (gap 0.0037%). 3/4 scenarios left.\n"," 0 0 445.06940 0 1 445.08574 445.06940 0.00% - 0s\n","H 0 0 445.0694028 445.06940 0.00% - 0s\n"," 0 0 578.28849 0 2 445.06940 578.28849 0.00% - 0s\n","\n","Optimal solution found at node 0 - now completing multiple scenarios...\n","\n"," Nodes | Current Node | Scenario Obj. Bounds | Work\n"," | | Worst |\n"," Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n","\n"," 0 2 579.16759 0 2 - 579.16759 - - 0s\n","Scenario 0 has been solved. 2/4 scenarios left.\n","Scenario 2 has been solved. 1/4 scenarios left.\n","Scenario 1 has been solved. 0/4 scenarios left.\n","\n","Cutting planes:\n"," Implied bound: 7\n"," MIR: 19\n"," Flow cover: 26\n"," Relax-and-lift: 1\n","\n","Explored 3 nodes (1835 simplex iterations) in 0.31 seconds (0.08 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 10: 445.069 445.069 445.086 ... 670.871\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 4.450694027691e+02, best bound 4.450694027691e+02, gap 0.0000%\n"]}],"source":["mm.write('ms.lp')\n","mm.optimize()\n","\n","#print obj value from each scenario\n","for s in range(m.NumScenarios):\n"," # Set the scenario number to query the information for this scenario\n"," mm.Params.ScenarioNumber = s\n"," print('\\nTotal cost for '+mm.ScenNName+' is '+'$'+str(round(mm.ScenNObjVal,2)))"]},{"cell_type":"markdown","metadata":{"id":"KM6o-kdx3YbO"},"source":["### Solution Pools\n","Dive in, the water's fine!\n","\n","A solution pool allows you to find multiple solutions that satisfy some criteria base off of search parameters set."]},{"cell_type":"code","execution_count":145,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":209,"status":"ok","timestamp":1699603570635,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"tYvUI5HeAVg9","outputId":"6c40ac97-549c-4bf8-8d9c-fe3c03043100"},"outputs":[{"name":"stdout","output_type":"stream","text":["Gurobi Optimizer version 10.0.3 build v10.0.3rc0 (linux64)\n","\n","CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]\n","Thread count: 1 physical cores, 2 logical processors, using up to 2 threads\n","\n","Optimize a model with 1440 rows, 1800 columns and 4498 nonzeros\n","Model fingerprint: 0x66bddde2\n","Variable types: 1440 continuous, 360 integer (360 binary)\n","Coefficient statistics:\n"," Matrix range [9e-01, 2e+01]\n"," Objective range [6e-02, 2e+00]\n"," Bounds range [1e+00, 8e+01]\n"," RHS range [1e-01, 9e+01]\n","Presolved: 1381 rows, 1741 columns, 4306 nonzeros\n","\n","Continuing optimization...\n","\n","\n","Explored 3964 nodes (8111 simplex iterations) in 0.02 seconds (0.00 work units)\n","Thread count was 2 (of 2 available processors)\n","\n","Solution count 500: 601.976 601.976 601.976 ... 602.035\n","\n","Optimal solution found (tolerance 1.00e-04)\n","Best objective 6.019761240991e+02, best bound 6.019761240991e+02, gap 0.0000%\n"]}],"source":["# Limit how many solutions to collect\n","mp.setParam(GRB.Param.PoolSolutions, 250)\n","# Limit the search space by setting a gap for the worst possible solution that will be accepted\n","mp.setParam(GRB.Param.PoolGap, 0.05)\n","# do a systematic search for the k-best solutions\n","mp.setParam(GRB.Param.PoolSearchMode, 2)\n","\n","mp.optimize()"]},{"cell_type":"markdown","metadata":{"id":"vOaMTUBrRRqW"},"source":["Now we can query all 250 solutions for variables of interest and look for anything interesting."]},{"cell_type":"code","execution_count":null,"metadata":{"id":"DmSx15LsMs9Z"},"outputs":[],"source":["n = mp.SolCount\n","flow_250 = pd.DataFrame()\n","for i in range(n):\n"," mp.setParam(GRB.Param.SolutionNumber, i)\n"," for b in batteries:\n"," tmp = pd.DataFrame([[flow_in[b,t].Xn, flow_out[b,t], i] for t in time_periods], columns=[b+'_out',b+'_in','Scen'+str(i)])\n"," flow_250 = pd.concat([flow_250, tmp])"]},{"cell_type":"code","execution_count":null,"metadata":{"colab":{"base_uri":"https://localhost:8080/","height":245},"executionInfo":{"elapsed":149,"status":"error","timestamp":1699514549043,"user":{"displayName":"Jerry Yurchisin","userId":"04917630893127362098"},"user_tz":300},"id":"uyG2mtFOBHIl","outputId":"d254ff84-5808-4900-e498-57fe2d4c925a"},"outputs":[{"ename":"NameError","evalue":"ignored","output_type":"error","traceback":["\u001b[0;31m---------------------------------------------------------------------------\u001b[0m","\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)","\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfigure\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0ms0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msol_level\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Battery0'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'orange'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0ms1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msol_level\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Battery1'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'blue'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Battery State (kWhr)'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Time Period'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n","\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"]}],"source":["plt.figure(figsize=(12,5))\n","s0, = plt.plot(sol_level['Battery0'], c = 'orange')\n","s1, = plt.plot(sol_level['Battery1'], c = 'blue')\n","plt.ylabel('Battery State (kWhr)')\n","plt.xlabel('Time Period')\n","plt.legend([s0,s1],[\"Battery0\", \"Battery1\"])\n","plt.axhline(y=capacity['Battery0'], c='orange', linestyle='--', alpha = 0.5)\n","plt.axhline(y=capacity['Battery1'], c='blue', linestyle='--', alpha = 0.5)\n","print(f\"Periods at Battery0 at Full Capacity: {sum(sol_level['Battery0']==capacity['Battery0'])}\")\n","print(f\"Periods at Battery1 at Full Capacity: {sum(sol_level['Battery1']==capacity['Battery1'])}\");"]},{"cell_type":"markdown","metadata":{"id":"6RFo7kk39cil"},"source":[]},{"cell_type":"markdown","metadata":{"id":"iPu22TkK9cRE"},"source":[]},{"cell_type":"code","execution_count":null,"metadata":{"id":"Yf6YU88f9buP"},"outputs":[],"source":[]}],"metadata":{"colab":{"authorship_tag":"ABX9TyN80moD6M6voTt1zw7joXtQ","provenance":[]},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":0}