{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Large problems: Recipes and EAO strategy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this notebook we explore large problems (in terms of memory or run-time). We explain the main drivers for size & complexity and describe simplifications we found useful.\n", "\n", "EAO generates mathematical problems (LP or MIP) out of given assets or portfolios - then passing them to solvers (such as Gurobi, CPLEX or free solvers such as SCPI or HighS). In that sense it is very similar to other commercial or open source approaches.\n", "\n", "No matter what solver (or asset optimization framework) is used, once we model complex assets and/or large time spans, problems may hit memory or run-time limits. **EAO's strategy is not to hide away simplifications** but to keep the exact formulation of the problem. We leave it to the user, deciding where she can live with simplifications. \n", "\n", "In the following we give a **recipe, which levers to pull** to keep problems solvable." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Basics - setting things up" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import datetime as dt\n", "from copy import deepcopy\n", "import eaopack as eao\n", "import tracemalloc ### trace memory usage\n", "import time \n", "\n", "### load input data (price samples)\n", "data = pd.read_excel(\"data2024.xlsx\")\n", "data.set_index('datum', inplace=True)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "## keeping track of peak memory\n", "import psutil # import Python psutil module\n", "\n", "def memory_usage():\n", " process = psutil.Process()\n", " usage = process.memory_info().rss \n", " # Using memory_info() to check consumption\n", " return usage # Returning the memory in bytes\n", "\n", "def memory_used():\n", " return np.round(tracemalloc.get_traced_memory()[1]/(1024.**2),0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Large problems - focus on memory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Seemingly \"large\" problem, that actually isn't" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with setting up a typical \"large\" problem:\n", "\n", "* Let us optimize across one year - hourly resolution (for demo purpose - arguments also hold for longer periods)\n", "* 20 generation assets in the portfolio (we will start with very simple assets - that means no ramping restrictions or startups)\n", "* plus (later on) a storage\n", "* fixed demand to be covered" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "### setup \n", "number_of_generating_assets = 20\n", "Start = dt.datetime(2024,1,1)\n", "End = dt.datetime(2025,1,1) # one year (end, i.e. 01/01:2025 00:00) is not included)\n", "resolution = 'h' # hourly resolution" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we can easily steer the number of time steps modeled using the resolution (or frequency). For power, prices mostly live an a scale of 15min, 30min or hours and the resolution is needed. For gas or other commodities, days, weeks or months are enough. \n", "\n", "**Recommendation: Always to choose the lowest resolution/ frequency that captures your problem**. Trivial example: modelling a year on hourly level leads to 8.760 time steps, on daily level 365. A big difference! It may be worthwhile to compare hourly, 15min or daily optimization to decide what is enough. May save computation or waiting times." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us generate the portfolio. We use a simplistic way of generating generators with random capacity and generation costs - note that we start tracking memory usage" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Peak memory usage: 2 MB\n" ] } ], "source": [ "tracemalloc.start() ## start memory tracking\n", "# overall structure\n", "node_power = eao.assets.Node(name = 'power')\n", "\n", "# demand\n", "demand = eao.assets.SimpleContract(name = 'demand',\n", " nodes = node_power,\n", " min_cap = \"demand_curve\",\n", " max_cap = \"demand_curve\")\n", "\n", "all_assets = [demand] # collect for portfolio\n", "# supply - various assets\n", "## using random generation costs\n", "np.random.seed(42)\n", "for i in range(number_of_generating_assets):\n", " gen = eao.assets.SimpleContract(name = f'gen_{i}',\n", " nodes = node_power,\n", " min_cap = 0,\n", " max_cap = np.random.uniform(0.5,2),\n", " price = f'price_gen_{i}'\n", " )\n", " all_assets.append(deepcopy(gen))\n", " # add random prices to input data\n", " data[f'price_gen_{i}'] = np.random.uniform(5,20,len(data)) \n", "portf = eao.Portfolio(assets = all_assets)\n", "timegrid = eao.Timegrid(Start, End, freq=resolution)\n", "# stopping memory tracking\n", "print(\"Peak memory usage: \", f'{memory_used():,.0f}', \"MB\")\n", "tracemalloc.stop()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Setting up the portfolio does not require much memory. Basically little has happened but storing parameters in the objects.\n", "\n", "Now: We do the optimization. EAO will\n", "* Create the LP/MIP (potentially needing much memory)\n", "* solving the problem (potentially requiring much time)" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Peak memory usage: 149 MB\n" ] } ], "source": [ "tracemalloc.start() ## start memory tracking\n", "out = eao.optimize(portf, timegrid, data)\n", "# stopping memory tracking\n", "print(\"Peak memory usage: \", f'{memory_used():,.0f}', \"MB\")\n", "tracemalloc.stop()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that solving this setup does not constitute an issue in terms of memory used or time required. \n", "\n", "Why is that?\n", "* We used relatively simple assets with few restrictions. Therefore **little memory** is used despite a relatively large number of assets and many time steps\n", "* Relatively **quick to solve**. The problem is an LP (no binary variables e.g. for start ups). Such problems can typically solved efficiently\n", "\n", "Here, a few lines of code that allow us to check the nature of the optimization problem. Note that the above call \"out = eao.optimize(portf, timegrid, data)\" is a shortcut. We take it apart into the single steps to obtain the (1) optimization problem (2) do the optimization and (3) to extract the results in neat dataframes:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Is problem a MIP? False\n" ] } ], "source": [ "# (1) generate optim. problem\n", "opt = portf.setup_optim_problem(data, timegrid)\n", "### get info on problem type:\n", "print(f\"Is problem a MIP? {opt.is_MIP}\")\n", "#### to show the remaining steps for sake of completeness (not needed here)\n", "# (2) optimize\n", "# res = opt.optimize() ### commented out to save us some time\n", "# (3) extract result into neat dataframes - also portfolio and optim problem needed)\n", "#out = eao.io.extract_output(portf, opt, res)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... so: no MIP and pretty easy to solve. Let's remember for later how to do the check" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Storage: The most problematic asset when it comes to memory" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us add a storage to the portfolio. It's just one asset more, but we will see that a storage requires much memory when the optimzation problem is set up" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "\n", "# storage\n", "battery = eao.assets.Storage(name = 'battery',\n", " nodes= node_power,\n", " cap_out = 1, \n", " cap_in = 1, \n", " size = 2, \n", " eff_in = 0.9, \n", " no_simult_in_out = False, # creates MIP of true! To be analyzed later\n", " block_size = None) # optimization over whole time horizon. To be analyzed later\n" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Peak memory usage: 3,534 MB\n" ] } ], "source": [ "tracemalloc.start() ## start memory tracking\n", "battery.setup_optim_problem(prices = None, timegrid=eao.Timegrid(Start, End, freq=resolution))\n", "# stopping memory tracking\n", "print(\"Peak memory usage: \", f'{memory_used():,.0f}', \"MB\")\n", "tracemalloc.stop()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Much memory is needed - even to generate the LP part of the storage only!\n", "\n", "To provide some detail on the \"why\":\n", "* in each point in time (every hour!) the storage is at least empty and at most full\n", "* that means, the sum of all previous dispatches must be >= 0 and <= storage size\n", "\n", "So we have two restrictions for each time step, involving all previous time steps. A densely populated matrix\n", "\n", "**But wait: There's a very simple way around**: Is it realistic to assume the dispatch (e.g.) on Jan 1 affects (e.g.) my decision on Dec 31? In reality we can securely assume that a short-term storage is dispatched on a daily or weekly basis. This should even reflect trading reality better, than optimizing against the full time span. In the following we try it, using the parameter **block_size**." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Peak memory usage: 158 MB\n" ] } ], "source": [ "battery.block_size = 'W' # weekly optimization. In each block, the battery goes back to start_level\n", "\n", "tracemalloc.start() ## start memory tracking\n", "battery.setup_optim_problem(prices = None, timegrid=eao.Timegrid(Start, End, freq=resolution))\n", "# stopping memory tracking\n", "print(\"Peak memory usage: \", f'{memory_used():,.0f}', \"MB\")\n", "tracemalloc.stop()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Recommendation:** Always utilize the block_size and set it to a low value. We found that daily (\"d\") is often a good choice for batteries or heat storage. Most available optimization tools utilize such simplification \"under the hood\". However, we believe it should be up to the user to decide whether or when to use it.\n", "\n", "**Long-term storage:** An alternative for long-term storage such as seasonal storages for gas or large hydro storages: We can define different intrinsic frequencies for single assets. For a specific sample check out https://rivacon.github.io/EAO/samples/mixing_long_and_short_term/mixing_long_and_short_term.html" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Solvers\n", "\n", "For MIPs, solvers typically implement complex iterations using heuristics - and widely differ in efficiency. In EAO, for MIPs, we set the standard to the open source solver SCIP. For LPs we let CVXPY decide If you have a license for Gurobi or CPLEX, it may be worthwhile to switch. Simply set the parameter \"solver = ...\" when calling the optimization (see samples below).\n", "\n", "**We are using CVXPY as interface to solvers**. See https://www.cvxpy.org/tutorial/solvers/index.html for available solvers. You may need to install solvers to your environment." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Run-time - divide & conquer and LPs vs. MIPs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us now focus on run-time. While the size of a problem (roughly number of assets $\\times$ number of time steps, storage block size) determines the problem size, its complexity is not necessarily linked to it.\n", "\n", "A main lever for complexity is the question whether we're dealing with a **mixed integer problem (MIP) or a linear problem (LP)**:\n", "* An LP only contains continuous variables. An LP is convex, with guaranteed convergence and efficient solution algorithms.\n", "* A MIP contains \"switch\" variables such as start-up decisions. MIPs are not necessarily convex, no convergence guaranteed. Run-time may be long, even for smaller problems.\n", "\n", "Unfortunately there's no easy way out and it turns out to be difficult knowing in advance, which problem is hard and which is not. In the coming sections we explore how to deal with complex problems." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Long simulations - splitting in chunks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solving MIPs over long time frames is often not possible in one go (due to long run-time or memory usage). However, solving them in chunks of, say, weeks or days, is mostly doable. Particularly in power or gas setups, there is also no reason to believe that decisions early the year will influence decisions late in the year (e.g. starting a plant in January will hardly have an effect on whether the plant is running in December).\n", "\n", "Splitting a problem into chunks manually would be cumbersome. Therefore we introduce \"split_optimization\" in EAO. This will cut the problem into chunks of user-defined length (e.g. weeks) and set together results automatically.\n", "\n", "Attention in cases where we are using restrictions across a longer time frame - such as \"minimum take\". In this case, split-optimization will also divide down this restriction. E.g. a minimum take of 12 GWh across the year will give rise to 12 chunks of 1 GWh minimum take per month. In those case consider manual break-down or going for an LP (see below)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us define a sample portfolio: The main asset is the power plant (we're adding two of them with different parameters), where we use some of the more complex restrictions. Those are mainly start up costs and specific ramps for starting and shutting down as well as minimum down time. We add a battery with a maximum number of cycles per day to add complexity." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "Start = dt.datetime(2024,1,1)\n", "End = dt.datetime(2025,1,1) # one month - for purpose of this demo a bit shorter \n", "resolution = 'h' # hourly resolution" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "# define timegrid\n", "timegrid = eao.assets.Timegrid(Start, End, freq=resolution)\n", "\n", "### our sample here: optimized power plant\n", "# define portfolio\n", "node_power = eao.assets.Node('node_power')\n", "power_market = eao.assets.SimpleContract(name = 'market', nodes=node_power, price='price',\n", " min_cap = -110, # we restrict grid connection so battery & plant cannot\n", " max_cap = 110 ) # run full load together\n", "plant_1 = eao.assets.Plant(name = 'plant_1',\n", " nodes = node_power,\n", " extra_costs = 100,\n", " min_cap = 20.,\n", " max_cap = 100.,\n", " ramp = 20,\n", " start_costs = 1000, \n", " min_downtime = 8,\n", " time_already_running = 0,\n", " time_already_off = 1) \n", "\n", "plant_2 = eao.assets.Plant(name = 'plant_2',\n", " nodes = node_power,\n", " extra_costs = 90,\n", " min_cap = 50.,\n", " max_cap = 90.,\n", " ramp = 10,\n", " start_costs = 2000, \n", " min_downtime = 10,\n", " time_already_running = 0,\n", " time_already_off = 1) \n", "\n", "\n", "battery = eao.assets.Storage(name = 'battery',\n", " nodes = node_power,\n", " cap_out = 50, # MW\n", " cap_in = 50, # MW\n", " size = 100, # MWh\n", " eff_in = 0.9, # %\n", " start_level = 0.5, # MWh\n", " end_level = 0.5, # MWh\n", " max_cycles_freq = 'd',\n", " max_cycles_no = 1,\n", " no_simult_in_out = True, # creates MIP if true! If false, simultaneous\n", " # charge/discharge for negative prices possible\n", " block_size = 'd') # daily optimization\n", "portf = eao.portfolio.Portfolio([power_market, plant_1, plant_2, battery])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's try different variants: Optimization in one piece, split into weeks, relaxed to LP as well as different solvers." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Run-time: 10 min 14 sec\n", "Total dispatch plant (no splitting using SCIP as EAO default solver): 210,996 MWh\n" ] } ], "source": [ "start_time = time.time()\n", "out = eao.optimize(portf, timegrid, data)\n", "end_time = time.time()\n", "minutes, seconds = divmod(end_time - start_time, 60)\n", "print(f\"Run-time: {int(minutes)} min {seconds:.0f} sec\")\n", "print(f\"Total dispatch plant (no splitting using SCIP as EAO default solver): {out['dispatch'][['plant_1', 'plant_2']].sum().sum():,.0f} MWh\")" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Run-time: 11 min 56 sec\n", "Total dispatch plant (no splitting using HighS as solver): 210,896 MWh\n" ] } ], "source": [ "start_time = time.time()\n", "out = eao.optimize(portf, timegrid, data, solver = 'HIGHS')\n", "end_time = time.time()\n", "minutes, seconds = divmod(end_time - start_time, 60)\n", "print(f\"Run-time: {int(minutes)} min {seconds:.0f} sec\")\n", "print(f\"Total dispatch plant (no splitting using HighS as solver): {out['dispatch'][['plant_1', 'plant_2']].sum().sum():,.0f} MWh\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see little difference between SCIP and HighS (we found both to be very efficient for EAO purposes). Note that plant output is very similar, but not identical (you will hardly achieve perfect accuracy in MIPs). It's worthwhile trying different solvers for your problems. Differences may be big. \n", "\n", "Particularly if you have a license for CPLEX, Gurobi or another commercial solver, it's worth trying it. Performance may be much better.\n", "\n", "Now the run in weekly chunks:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Run-time: 2 min 24 sec\n", "Total dispatch plant (with splitting): 208,496 MWh\n" ] } ], "source": [ "start_time = time.time()\n", "out = eao.optimize(portf, timegrid, data, split_interval_size='W')\n", "end_time = time.time()\n", "minutes, seconds = divmod(end_time - start_time, 60)\n", "print(f\"Run-time: {int(minutes)} min {seconds:.0f} sec\")\n", "print(f\"Total dispatch plant (with splitting): {out['dispatch'][['plant_1', 'plant_2']].sum().sum():,.0f} MWh\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Check out the comparison of overall generation: the split optimization leads to similar results. A close look on hourly level may, of course, show greater differences in detail. For example: running through a weekend instead of shutting down for a Sunday night.\n", "\n", "Note: In this illustration we chose examples with run-times that do not hurt too much in an interactive notebook. But note that for complex assets which require a MIP, split optimization easily makes a great difference (also between no convergence at all and quick results).\n", "\n", "**Recommendation:** Use split optimization as far as possible for simulation across long time frames. Check your asset definitions and your analysis targets to guide you. Particularly where assets do not have long-range restrictions (such as minimum take or long-term storage), split optimization is a good choice." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Making the MIP an LP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A MIP is present if, for example, we have a binary decision such as whether to start oa plant that is tied to costs such as start-up costs. We utilize a boolean variables (0 = off, 1 = on) or (1 = is starting, 0 = is not starting).\n", "\n", "A drastic simplification is to allow all values [0,1], making the problem an LP. Sometimes cutting down drastically on run-time. The main use-case, though, is for testing purpose in cases where a MIP is not feasible or does not converge at all. For example, decision variables such as for start-up to a minimum load may lead to cases where a small demand cannot be fulfilled (ill formulated problem). In those cases, a relaxed problem may help to identify mistakes or workarounds.\n", "\n", "However, also for some use-cased, we believe this can be a valid approximation, particularly in long-term valuation.\n", "\n", "Let us test, using the parameter make_soft_problem:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Run-time: 9 min 3 sec\n", "Total dispatch plant (softened problem using default solver): 297,740 MWh\n" ] } ], "source": [ "start_time = time.time()\n", "out = eao.optimize(portf, timegrid, data, make_soft_problem=True)\n", "end_time = time.time()\n", "minutes, seconds = divmod(end_time - start_time, 60)\n", "print(f\"Run-time: {int(minutes)} min {seconds:.0f} sec\")\n", "print(f\"Total dispatch plant (softened problem using default solver): {out['dispatch'][['plant_1', 'plant_2']].sum().sum():,.0f} MWh\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our simple illustration, we see a limited impact run-time - yet in real-life problems differences may be much higher. The approximation appears somewhat ok, at least in terms of overall plant dispatch.\n", "\n", "**Recommendation:** Use make_soft_problem with care. When a problem does not converge or when it is not feasible (it easily happens that restrictions cannot be met due to errors in the setup), start trying the soft version. When using for a productive purposes, do thorough tests. Start-up processes, for example, will not be modelled correctly." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### " ] } ], "metadata": { "kernelspec": { "display_name": "eao", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.0" }, "orig_nbformat": 2 }, "nbformat": 4, "nbformat_minor": 2 }