diff --git a/data-science-for-esm/14-workshop-linopy.ipynb b/data-science-for-esm/14-workshop-linopy.ipynb new file mode 100644 index 00000000..5a675eb9 --- /dev/null +++ b/data-science-for-esm/14-workshop-linopy.ipynb @@ -0,0 +1,5660 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "47098e73-2dbb-4b9d-8abf-c2de66e1f264", + "metadata": {}, + "source": [ + "# Introduction to `linopy`" + ] + }, + { + "cell_type": "markdown", + "id": "cdd93543-0b5c-4e03-adc5-5c9d44c7774a", + "metadata": {}, + "source": [ + ":::{note}\n", + "This material is in part adapted from the following resources:\n", + "- [Linopy Getting Started](https://linopy.readthedocs.io/en/latest/index.html)\n", + "- [PyPSA simple electricity market examples](https://pypsa.readthedocs.io/en/latest/examples/simple-electricity-market-examples.html)\n", + ":::\n", + "\n", + "\n", + "\n", + "[Linopy](https://linopy.readthedocs.io/en/latest/index.html) is an open-source framework for formulating, solving, and analyzing optimization problems with Python.\n", + "\n", + "With Linopy, you can create optimization models within Python that consist of decision variables, constraints, and optimization objectives. You can then solve these instances using a variety of commercial and open-source solvers (specialised software).\n", + "\n", + "[Linopy](https://linopy.readthedocs.io/en/latest/index.html) supports a wide range of problem types, including:\n", + "\n", + "- **Linear programming**\n", + "- Integer programming\n", + "- Mixed-integer programming\n", + "- Quadratic programming\n", + "\n", + "\n", + ":::{note}\n", + "Documentation for this package is available at https://linopy.readthedocs.io.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "id": "e6221514-355e-40eb-9194-b9026add6e40", + "metadata": {}, + "source": [ + ":::{note}\n", + "If you have not yet set up Python on your computer, you can execute this tutorial in your browser via [Google Colab](https://colab.research.google.com/). Click on the rocket in the top right corner and launch \"Colab\". If that doesn't work download the `.ipynb` file and import it in [Google Colab](https://colab.research.google.com/).\n", + "\n", + "Then install the following packages by executing the following command in a Jupyter cell at the top of the notebook.\n", + "\n", + "```sh\n", + "!pip install pandas linopy highspy\n", + "```\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "id": "72fb5873-a86d-45e4-9c85-eca1b96d1da6", + "metadata": {}, + "source": [ + "## Solve a Basic Model\n", + "\n", + "In this example, we explain the basic functions of the linopy `Model` class. First, we are setting up a very simple linear optimization model, given by\n", + "\n", + "Minimize:\n", + " $$x + 2y$$\n", + "subject to:\n", + " $$ x \\ge 0 $$\n", + " $$y \\ge 0 $$\n", + " $$3x + 7y \\ge 10 $$\n", + " $$5x + 2y \\ge 3 $$" + ] + }, + { + "cell_type": "markdown", + "id": "cc582211-9363-4107-932c-48c9c2279535", + "metadata": {}, + "source": [ + "### Initializing a `Model`\n", + "\n", + "The Model class in Linopy is a fundamental part of the library. It serves as a container for all the relevant data associated with a linear optimization problem. This includes variables, constraints, and the objective function." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d94292b1-8073-49a8-a2e9-a2c769145b2d", + "metadata": {}, + "outputs": [], + "source": [ + "import linopy\n", + "\n", + "m = linopy.Model()" + ] + }, + { + "cell_type": "markdown", + "id": "bf2b2258-f0f3-4118-ac39-6d902e1d61e3", + "metadata": {}, + "source": [ + "This creates a new Model object, which you can then use to define your optimization problem." + ] + }, + { + "cell_type": "markdown", + "id": "f0f00003", + "metadata": {}, + "source": [ + ":::{note}\n", + "It is good practice to choose a short variable name (like `m`) to reduce the verbosity of your code.\n", + ":::" + ] + }, + { + "cell_type": "markdown", + "id": "4844ce4c-4ba6-4dec-957a-69dce93c07f0", + "metadata": {}, + "source": [ + "### Adding decision variables\n", + "\n", + "**Variables** are the unknowns of an optimisation problems and are intended to be given values by solving an optimisation problem. A variable can always be assigned with a lower and an upper bound. In our case, both `x` and `y` have a lower bound of zero (default is unbouded). In linopy, you can add variables to a `Model` using the `add_variables()` method:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "64b35f1c-5398-4eb3-9135-aafce270b79b", + "metadata": {}, + "outputs": [], + "source": [ + "x = m.add_variables(lower=0, name=\"x\")\n", + "y = m.add_variables(lower=0, name=\"y\");" + ] + }, + { + "cell_type": "markdown", + "id": "57bd1e94-7ead-4c86-a590-f2c2998b5f92", + "metadata": {}, + "source": [ + "`x` and `y` are linopy variables of the class `linopy.Variable`. Each of them contain all relevant information that define it. The `name` parameter is optional but can be useful for referencing the variables later." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "19d4b5f5-b9ae-48f6-8ce6-cd4a418dba6b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable\n", + "--------\n", + "x ∈ [0, inf]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "f6c58398", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "linopy.model.Variables\n", + "----------------------\n", + " * x\n", + " * y" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.variables" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "c7f01250", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable\n", + "--------\n", + "x ∈ [0, inf]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.variables[\"x\"]" + ] + }, + { + "cell_type": "markdown", + "id": "d07cc8d8-1ea6-47e2-9b84-5326948a602b", + "metadata": {}, + "source": [ + "### Adding Constraints\n", + "\n", + "**Constraints** are equality or inequality expressions that define the *feasible* space of the decision variables. They consist of the left hand side (LHS) and the right hand side (RHS). The first constraint that we want to write down is $3x + 7y = 10$ which we write out exactly in the mathematical way:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "1df4589f-4c99-46b5-a530-63f99b4515cc", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint (unassigned)\n", + "-----------------------\n", + "+3 x + 7 y ≥ 10.0" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "3 * x + 7 * y >= 10" + ] + }, + { + "cell_type": "markdown", + "id": "2b7ddd9a-8601-4a06-9d40-b8760651bd2a", + "metadata": {}, + "source": [ + "Note, we can also mix the constant and the variable expression, like this" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "688e0f65-9198-49b2-915a-2131084224c1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint (unassigned)\n", + "-----------------------\n", + "+3 x + 7 y ≥ 10.0" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "3 * x + 7 * y - 10 >= 0" + ] + }, + { + "cell_type": "markdown", + "id": "a45ca5fa-0696-4688-a1e2-16c8eb06644b", + "metadata": {}, + "source": [ + "… and linopy will automatically take over the separation of variables expression on the LHS, and constant values on the RHS." + ] + }, + { + "cell_type": "markdown", + "id": "f18d18a3-671f-4190-9257-4541e6c88c61", + "metadata": {}, + "source": [ + "The constraint is currently not assigned to the model. We assign it by calling the `add_constraints()` function:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "fa233f5c-38ec-4a5a-9008-502932df968d", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_constraints(3 * x + 7 * y >= 10)\n", + "m.add_constraints(5 * x + 2 * y >= 3);" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "25c7a077", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "linopy.model.Constraints\n", + "------------------------\n", + " * con0\n", + " * con1" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.constraints" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "024a26c4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `con0`\n", + "-----------------\n", + "+3 x + 7 y ≥ 10.0" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.constraints[\"con0\"]" + ] + }, + { + "cell_type": "markdown", + "id": "f7fc48a8-ab7e-4baf-93d2-ff907e543db1", + "metadata": {}, + "source": [ + "### Adding the Objective \n", + "\n", + "The objective function defines what you want to optimize. It is a function of variables that a solver attempts to maximize or minimize. You can set the objective function of a `linopy.Model` using the `add_objective()` method. For our example that would be" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "18c6e6ce-37f9-4527-83a4-e5764b67b34c", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_objective(x + 2 * y, sense=\"min\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "4040afc5", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Objective:\n", + "----------\n", + "LinearExpression: +1 x + 2 y\n", + "Sense: min\n", + "Value: None" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.objective" + ] + }, + { + "cell_type": "markdown", + "id": "014ba3b0-afa4-46bf-a84a-0553eef830a5", + "metadata": {}, + "source": [ + "Note, we can either minimize or maximize in linopy. Per default, linopy applies `sense='min'` making it not necessary to explicitly define the optimization sense. In summary:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "a1e8788b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Linopy LP model\n", + "===============\n", + "\n", + "Variables:\n", + "----------\n", + " * x\n", + " * y\n", + "\n", + "Constraints:\n", + "------------\n", + " * con0\n", + " * con1\n", + "\n", + "Status:\n", + "-------\n", + "initialized" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "markdown", + "id": "b8ec9c21-de8c-40a1-8d30-2e3ee3f6d7ed", + "metadata": {}, + "source": [ + "### Solving the Model\n", + "\n", + "Once you've defined your `linopy.Model` with variables, constraints, and an objective function, you can solve it using the `solve` method:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "712a0033-d027-478a-b8f2-201ec4ed1cc1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "2 rows, 2 cols, 4 nonzeros\n", + "2 rows, 2 cols, 4 nonzeros\n", + "Presolve : Reductions: rows 2(-0); columns 2(-0); elements 4(-0) - Not reduced\n", + "Problem not reduced by presolve: solving the LP\n", + "Using EKK dual simplex solver - serial\n", + " Iteration Objective Infeasibilities num(sum)\n", + " 0 0.0000000000e+00 Pr: 2(13) 0s\n", + " 2 2.8620689655e+00 Pr: 0(0) 0s\n", + "Model status : Optimal\n", + "Simplex iterations: 2\n", + "Objective value : 2.8620689655e+00\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "c314737f", + "metadata": {}, + "source": [ + "Solvers are needed to compute solutions to the optimization models. There exists a large variety of solvers. In many cases, they specialise in certain problem types or solving algorithms, e.g. linear or nonlinear problems.\n", + "\n", + "- **open-source examples**: [CBC](https://www.coin-or.org/Cbc/), [GLPK](https://www.gnu.org/software/glpk/), [Ipopt](https://coin-or.github.io/Ipopt/), [HiGHS](https://highs.dev)\n", + "- **commercial examples**: [Gurobi](https://www.gurobi.com/), [CPLEX](https://www.ibm.com/de-de/analytics/cplex-optimizer), [FICO Xpress](https://www.fico.com/en/products/fico-xpress-optimization)\n", + "\n", + "The open-source solvers are sufficient to handle meaningful linopy models with hundreds to several thousand variables and constraints. However, as applications get large or more complex, there may be a need to turn to a commercial solvers (which often provide free academic licenses).\n", + "\n", + "For this course, we use HiGHS, which is already in the course environment `esm-2024`." + ] + }, + { + "cell_type": "markdown", + "id": "ca69649b-0fa6-4ca6-a557-8eb913895b19", + "metadata": {}, + "source": [ + "### Retrieving optimisation results\n", + "\n", + "The solution of the linear problem is assigned to the variables under `solution` in form of a `xarray.Dataset`." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "1f28aa81-d16f-4b7e-9d8e-f661d4b8bcd8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'solution' ()> Size: 8B\n",
+       "array(0.03448276)
" + ], + "text/plain": [ + " Size: 8B\n", + "array(0.03448276)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x.solution" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "470e2d25-9386-438d-932e-113464907726", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'solution' ()> Size: 8B\n",
+       "array(1.4137931)
" + ], + "text/plain": [ + " Size: 8B\n", + "array(1.4137931)" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "y.solution" + ] + }, + { + "cell_type": "markdown", + "id": "5de5e02c", + "metadata": {}, + "source": [ + "We can also read out the objective value:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "b30ceb74", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "2.8620689655172415" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.objective.value" + ] + }, + { + "cell_type": "markdown", + "id": "5f52e8fb", + "metadata": {}, + "source": [ + "And the dual values (or shadow prices) of the model's constraints: " + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "f6604954", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'con0' ()> Size: 8B\n",
+       "array(0.27586207)
" + ], + "text/plain": [ + " Size: 8B\n", + "array(0.27586207)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.dual[\"con0\"]" + ] + }, + { + "cell_type": "markdown", + "id": "34b0c100-b354-4e8c-8043-a45bfa78b999", + "metadata": {}, + "source": [ + "Well done! You solved your first linopy model!" + ] + }, + { + "cell_type": "markdown", + "id": "8d87cf68-8655-4f04-b87c-ab1e670f2884", + "metadata": {}, + "source": [ + "## Use Coordinates\n", + "\n", + "Now, the real power of the package comes into play! \n", + "\n", + "Linopy is structured around the concept that variables, and therefore expressions and constraints, have coordinates. That is, a `Variable` object actually contains multiple variables across dimensions, just as we know it from a `numpy` array or a `pandas.DataFrame`.\n", + "\n", + "Suppose the two variables `x` and `y` are now functions of time `t` and we would modify the problem according to: " + ] + }, + { + "cell_type": "markdown", + "id": "d40eb0b9-35ef-4867-b57c-dfca3b0da66c", + "metadata": {}, + "source": [ + "Minimize:\n", + "$$\\sum_t x_t + 2 y_t$$\n", + "\n", + "subject to:\n", + "\n", + "$$x_t \\ge 0 \\qquad \\forall t $$\n", + "$$y_t \\ge 0 \\qquad \\forall t $$\n", + "$$3x_t + 7y_t \\ge 10 t \\qquad \\forall t$$\n", + "$$5x_t + 2y_t \\ge 3 t \\qquad \\forall t$$\n", + "\n", + "whereas `t` spans all the range from 0 to 10." + ] + }, + { + "cell_type": "markdown", + "id": "53c08a98-05aa-4c0b-a046-1d16215ce382", + "metadata": {}, + "source": [ + "In order to formulate the new problem with linopy, we start again by initializing a model." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "788956e7-1bc9-4298-85a1-9235486b8cad", + "metadata": {}, + "outputs": [], + "source": [ + "m = linopy.Model()" + ] + }, + { + "cell_type": "markdown", + "id": "9cf08174-0d0a-473c-8ee2-7799958aef0b", + "metadata": {}, + "source": [ + "Again, we define `x` and `y` using the `add_variables()` function, but now we are adding a `coords` argument. This automatically creates optimization variables for all coordinates, in this case time-steps `t`." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "d477e41c-a89f-4d3b-a1af-385820638a75", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "time = pd.Index(range(10), name=\"time\")\n", + "\n", + "x = m.add_variables(\n", + " lower=0,\n", + " coords=[time],\n", + " name=\"x\",\n", + ")\n", + "y = m.add_variables(lower=0, coords=[time], name=\"y\")" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "b31a72ca", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable (time: 10)\n", + "-------------------\n", + "[0]: x[0] ∈ [0, inf]\n", + "[1]: x[1] ∈ [0, inf]\n", + "[2]: x[2] ∈ [0, inf]\n", + "[3]: x[3] ∈ [0, inf]\n", + "[4]: x[4] ∈ [0, inf]\n", + "[5]: x[5] ∈ [0, inf]\n", + "[6]: x[6] ∈ [0, inf]\n", + "[7]: x[7] ∈ [0, inf]\n", + "[8]: x[8] ∈ [0, inf]\n", + "[9]: x[9] ∈ [0, inf]" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x" + ] + }, + { + "cell_type": "markdown", + "id": "32dfa7f0-3f35-466d-8669-6134ca18a26d", + "metadata": {}, + "source": [ + "Following the previous example, we write the constraints out using the syntax from above, while multiplying the RHS with `t`. Note that the coordinates from the LHS and the RSH have to match. \n", + "\n", + ":::{note}\n", + "In the beginning, it is recommended to use explicit dimension names. In this way, things remain clear and no unexpected broadcasting (which we show later) will happen.\n", + ":::" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "00a64dec-26a5-4f80-97c0-de0fe00188a2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint (unassigned) (time: 10):\n", + "-----------------------------------\n", + "[0]: +3 x[0] + 7 y[0] ≥ -0.0\n", + "[1]: +3 x[1] + 7 y[1] ≥ 10.0\n", + "[2]: +3 x[2] + 7 y[2] ≥ 20.0\n", + "[3]: +3 x[3] + 7 y[3] ≥ 30.0\n", + "[4]: +3 x[4] + 7 y[4] ≥ 40.0\n", + "[5]: +3 x[5] + 7 y[5] ≥ 50.0\n", + "[6]: +3 x[6] + 7 y[6] ≥ 60.0\n", + "[7]: +3 x[7] + 7 y[7] ≥ 70.0\n", + "[8]: +3 x[8] + 7 y[8] ≥ 80.0\n", + "[9]: +3 x[9] + 7 y[9] ≥ 90.0" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "factor = pd.Series(time, index=time)\n", + "\n", + "3 * x + 7 * y >= 10 * factor" + ] + }, + { + "cell_type": "markdown", + "id": "946c5674-d5f5-4b78-97d3-337e6de10f73", + "metadata": {}, + "source": [ + "It always helps to write out the constraints before adding them to the model. Since they look good, let's assign them." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "8d3a89e8-0b0e-480d-9cb8-f29931fb3559", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Linopy LP model\n", + "===============\n", + "\n", + "Variables:\n", + "----------\n", + " * x (time)\n", + " * y (time)\n", + "\n", + "Constraints:\n", + "------------\n", + " * con1 (time)\n", + " * con2 (time)\n", + "\n", + "Status:\n", + "-------\n", + "initialized" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "con1 = m.add_constraints(3 * x + 7 * y >= 10 * factor, name=\"con1\")\n", + "con2 = m.add_constraints(5 * x + 2 * y >= 3 * factor, name=\"con2\")\n", + "m" + ] + }, + { + "cell_type": "markdown", + "id": "ff3e5770-383f-4011-aed7-7de526e69c7c", + "metadata": {}, + "source": [ + "Now, when it comes to the objective, we use the `sum` function of `linopy.LinearExpression`. This stacks all terms all terms of the `time` dimension and writes them into one big expression. " + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "9cbf16b4-99ed-4e33-9ea9-e8eff3503b05", + "metadata": {}, + "outputs": [], + "source": [ + "obj = (x + 2 * y).sum()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "075a8a0b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LinearExpression\n", + "----------------\n", + "+1 x[0] + 2 y[0] + 1 x[1] ... +2 y[8] + 1 x[9] + 2 y[9]" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "obj" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "8f4e2168", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_objective(obj, overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "37da40ed-71f5-4c3c-825b-26e0e2da46e2", + "metadata": {}, + "source": [ + "Then, we can solve:" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "436f52a8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "18 rows, 18 cols, 36 nonzeros\n", + "18 rows, 18 cols, 36 nonzeros\n", + "Presolve : Reductions: rows 18(-2); columns 18(-2); elements 36(-4)\n", + "Solving the presolved LP\n", + "Using EKK dual simplex solver - serial\n", + " Iteration Objective Infeasibilities num(sum)\n", + " 0 0.0000000000e+00 Pr: 18(585) 0s\n", + " 18 1.2879310345e+02 Pr: 0(0) 0s\n", + "Solving the original LP from the solution after postsolve\n", + "Model status : Optimal\n", + "Simplex iterations: 18\n", + "Objective value : 1.2879310345e+02\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "a5f4769a-98d0-446f-81cf-b41d7c4063e2", + "metadata": {}, + "source": [ + "In order to inspect the solution. You can go via the variables, i.e. `y.solution` or via the `solution` aggregator of the model, which combines the solution of all variables." + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "id": "33cdfad9-0ff3-4211-afaf-b0e27fa33d5a", + "metadata": {}, + "outputs": [ + { + "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", + " \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", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
xy
time
00.0000000.000000
10.0344831.413793
20.0689662.827586
30.1034484.241379
40.1379315.655172
50.1724147.068966
60.2068978.482759
70.2413799.896552
80.27586211.310345
90.31034512.724138
\n", + "
" + ], + "text/plain": [ + " x y\n", + "time \n", + "0 0.000000 0.000000\n", + "1 0.034483 1.413793\n", + "2 0.068966 2.827586\n", + "3 0.103448 4.241379\n", + "4 0.137931 5.655172\n", + "5 0.172414 7.068966\n", + "6 0.206897 8.482759\n", + "7 0.241379 9.896552\n", + "8 0.275862 11.310345\n", + "9 0.310345 12.724138" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solution.to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "id": "0c262072", + "metadata": {}, + "source": [ + "Sometimes it can be helpful to plot the solution:" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "id": "5ba03b54", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGwCAYAAABcnuQpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABLeElEQVR4nO3de5xM9ePH8dfM7J1d69Kuy65bknuIShR+sSEi0YW+lVC+KL663/ENqb5SKVJaffvmW7lFymVTSCW3VN8UlbtoCbtYdmdnzu+Pw2btYmfN7Dmz834+Hh7tZ3b2zHv3s5x353PmHIdhGAYiIiIiQcppdQARERGR86EyIyIiIkFNZUZERESCmsqMiIiIBDWVGREREQlqKjMiIiIS1FRmREREJKiFWR0g0LxeL7///juxsbE4HA6r44iIiEgRGIbB4cOHqVq1Kk7n2Y+9lPoy8/vvv5OcnGx1DBERESmGnTt3kpSUdNbnlPoyExsbC5g/jLi4OL9u2+12s2TJElJSUggPD/frtsV3mg970XzYi+bDXjQf55aZmUlycnLefvxsSn2ZObm0FBcXF5AyExMTQ1xcnH4ZbUDzYS+aD3vRfNiL5qPoinKKiE4AFhERkaCmMiMiIiJBTWVGREREglqpP2emqDweD26326evcbvdhIWFcfz4cTweT4CSnZ/w8HBcLpfVMURERAIm5MuMYRjs3buXQ4cOFetrK1euzM6dO219DZv4+HgqV65s64wiIiLFFfJl5mSRSUhIICYmxqcdvtfr5ciRI5QtW/acF/SxgmEYZGVlkZ6eDkCVKlUsTiQiIuJ/IV1mPB5PXpGpWLGiz1/v9XrJyckhKirKlmUGIDo6GoD09HQSEhK05CQiIqWOPffAJeTkOTIxMTEWJwmsk9+fr+cEiYiIBIOQLjMnlfZzSUr79yciIqFNZUZERESCmsqMiIiIBDWVGREREQlqKjMiIiJSPIYBmxaB12tpDJUZERER8V32YZg9AP57M3z1kqVRVGZOYxgGWTm5Rf5zLMfj0/PP9McwjCJn3LdvH5UrV2bs2LF5j33zzTdERESwZMmSQPxYRERE/rLnO3j9avjfLHC4zD8WCumL5hXmmNtDg6cWl/jrbhx9LTERRZuOCy64gLfeeosePXqQkpJCvXr1uO222xg8eDApKSkBTioiIiHLMGDNm7D4cfBkQ1wS9HoLql9uaSyVmSDVpUsXBg4cSN++fWnZsiVRUVE8++yzVscSEZHS6ngGzL8XNs4zx3U7Q4/XIKaCtblQmSkgOtzFxtHXFum5Xq+Xw5mHiY2LPe/bGUSH+36I7oUXXqBRo0Z88MEHrF27lqioqPPKICIiUqjd62HmnXBoOzjDoMMoaDUEbHJRVpWZ0zgcjiIv93i9XnIjXMREhFlyb6YtW7bw+++/4/V62b59O02aNCnxDCIiUooZBnwzBZY8CV43lKsOvVMhqYXVyfJRmQlSOTk59O3bl5tvvpl69erRv39/fvjhBxITE62OJiIipcGxgzBvKPy8wBzX6wrdJ0F0eWtzFUJlJkg9/vjjZGRk8PLLL1O2bFkWLlxI//79WbBggdXRREQk2O1cA7Pugowd4IqAlDFw2UDbLCudTm/NDkLLli1j4sSJvPPOO8TFxeF0OnnnnXdYuXIlkydPtjqeiIgEK68XvnwZUjuZRaZ8Lei/BC6/27ZFBnRkJii1a9cOt9ud77Hq1atz6NAhawKJiEjwyzoAcwfBLycuT9LwBuj2EkSVszZXEajMiIiIhLrtX8Ps/pC5G1yR0GkctLjL1kdjTqUyIyIiEqq8XvjyRfhsDBgeqFgHek+Hyo2tTuYTlRkREZFQdGQfzL0HfltqjhvfBF0nQGSstbmKQWVGREQk1GxbCbP6w5G9EBYNXZ6HZrcFzbLS6VRmREREQoXXAytegOXPguGFSheby0qJDaxOdl5UZkRERELB4T9gzgDYusIcN+1rHpGJKGNtLj+w9DozK1asoFu3blStWhWHw8GHH36Y9zm3283DDz9M48aNKVOmDFWrVuX222/n999/ty6wiIhIMPrtc5jSxiwy4THQY4p5k8hSUGTA4jJz9OhRLrnkEiZNmlTgc1lZWaxfv54nn3yS9evXM2fOHDZv3sz1119vQVIREZEg5MmFz56Bd26Ao+mQ0ADuXg5Nb7U6mV9ZuszUuXNnOnfuXOjnypUrR1paWr7HXnnlFS677DJ27NhB9erVSyKiiIhIcMr8HWYPgO1fmuPmd0Dn8RAebW2uAAiqc2YyMjJwOBzEx8ef8TnZ2dlkZ2fnjTMzMwFz2er0q+a63W4Mw8Dr9eL1en3OYxhG3n+L8/Ulxev1YhgGbrcbl8tldZyAOTm/p8+zWEPzYS+aD3sJ9Hw4fluKa/5gHFl/YkSUwdNlAkbDG0++eEBe0998+dk4jJN7ZIs5HA7mzp1Ljx49Cv388ePHadOmDfXq1eM///nPGbczcuRIRo0aVeDxGTNmEBMTk++xsLAwKleuTHJyMhEREeeV385ycnLYuXMne/fuJTc31+o4IiISIA4jl3p75lD3D/Omw4eiq7O25lCORlW2OJnvsrKy6NOnDxkZGcTFxZ31uUFRZtxuN71792bHjh0sW7bsrN9UYUdmkpOT2b9/f4GvO378ODt37qRmzZpERUX5nNkwDA4fPkxsbCyOEnxv/r///W/uv/9+du3aRWRkZN7jvXr1okyZMrz99tv5nn/8+HG2bdtGcnJysb7PYOF2u0lLS6Njx46Eh4dbHSfkaT7sRfNhLwGZj8zduOYOxLlrNQCeS+/C22E0hAXnv/uZmZlUqlSpSGXG9stMbrebm266ia1bt/LZZ5+d8xuKjIzMt4M/KTw8vMAvjMfjweFw4HQ6cTpPnAttGODOKlI2r9cL7iwcbtdfX19c4TFFvljRzTffzPDhw1mwYAG9e/cGYP/+/Xz88ccsWrSoQBan04nD4Sj0Z1Aahcr3GSw0H/ai+bAXv83HpkXw4SA4dhAi4+D6V3A17EEwn1jgy8/F1mXmZJH55Zdf+Pzzz6lYsWIJvGgWjK1apKc6gXh/ve5jvxf5LXLR0dH06dOH1NTUvDLz7rvvkpSURLt27fyVSERE7C43B5aOgq9PvCu4ajPolQoValmbq4RZWmaOHDnCr7/+mjfeunUrGzZsoEKFClStWpVevXqxfv16FixYgMfjYe/evQBUqFChVJ/jUhQDBw6kZcuW7N69m2rVqpGamsqdd95ZostdIiJioYPbYdZdsHutOb7879BxFIQVXJ0o7SwtM2vXrqV9+/Z54xEjRgBwxx13MHLkSObPnw9A06ZN833d559/HrgjEOEx5lGSIvB6vWQePkxcbKx/lpl80KxZMy655BL+/e9/c+211/LDDz/w0UcfnV8GEREJDj99BPOGwPEMiCoH3V+D+l2tTmUZS8tMu3btONv5x5acm+xwFP2KiF4vhHvM559vmSmGAQMG8OKLL7J79246dOhAcnJyiWcQEZESlJsNS56E1a+b46SW0OstiA/ta69ZegVgOT99+/Zl9+7dvPHGG9x1111WxxERkUA6sAWmpfxVZK68F/otDPkiAyozQS0uLo4bb7yRsmXLnvH6PCIiUgr8OBdebwt7NkB0BejzAaQ8Ay69Mw1s/m4mObc9e/bQt2/fQt+OLiIiQc59HBY/BmunmePkK8xlpXLVrM1lMyozQerAgQMsWbKEzz77rNAbdYqISJDb/yvMvBP++MEctxkB7R8Hl3bdp9NPJEg1b96cgwcPMn78eC6++GKr44iIiD99PxMWDIecIxBTCXq+DnU6WJ3KtlRmgtS2bdusjiAiIv6WkwWLHob1/zbHNa+Cnm9AXBVrc9mcyoyIiIgd7NtkLiulbwQc0PYhaPswOIP5pgQlQ2UGi65nU4JK+/cnIhL0NsyAj+83b6lTJgFufBNqt7U6VdAI6TJz8iZWWVlZREdHW5wmcLKyzBtn6uZyIiI2k3MUFjwK380wx7XbmctKZRMsjRVsQrrMuFwu4uPjSU9PByAmJsanext5vV5ycnI4fvz4+d/OIAAMwyArK4v09HTi4+NxuXSoUkTELmKP7STsrQ7w5y/gcEK7x+CqEVpWKoaQLjMAlStXBsgrNL4wDINjx44RHR1t6xs8xsfH532fIiJiMcPA8e07tN00Eofhhtgq5rJSzTZWJwtaIV9mHA4HVapUISEhAbfb7dPXut1uVqxYwdVXX23bJZzw8HAdkRERsYvsw7DgH4T9MBMAb+1rcN44FcpUsjhYcAv5MnOSy+XyeafvcrnIzc0lKirKtmVGRERsYs/35ruVDvyG4XCxscqN1L1lEs4IXcH9fKnMiIiIBJJhmLcjWPQYeLIhLgnPDVP59fv91HXY73zLYKSfooiISKAczzCPxnx8v1lk6naGQV9gJF1mdbJSRUdmREREAmH3epjVDw5uA2cYdBgFrYaAwwE+nqMpZ6cyIyIi4k+GAd+8DkueAK8bylWH3qmQ1MLqZKWWyoyIiIi/HDsI84bCzwvMcb2u0H0SRJe3NlcppzIjIiLiD7vWwsx+kLEDXBGQ8gxcdre5rCQBpTIjIiJyPgwDvp4En44Eby6Urwm9p0PVZhYHCx0qMyIiIsWVdQA+/DtsXmSOG94A3V6CqHLW5goxKjMiIiLFsWMVzLoLMneDKxI6jYMWd2lZyQIqMyIiIr7weuHLifDZM2B4oMKF5rJSlSZWJwtZKjMiIiJFdXQ/zL0Hfv3UHDfuDV1fhMhYa3OFOJUZERGRoti2EmYPgMN7ICwaujwHzf6mZSUbUJkRERE5G68HvvgXLBsHhhcqXWwuKyU2sDqZnKAyIyIiciaH/4A5A2HrcnPctC90eR4iylibS/JRmRERESnMlmUweyAcTYfwGLhuAjS91epUUgiVGRERkVN5PbDsWVjxPGBAQgNzWemCi61OJmegMiMiInJS5h7zJN/tK81x8zug83gIj7Y2l5yVyoyIiAiYb7eeczdk/QkRZaHrRGjS2+pUUgQqMyIiEto8ufD5M7DyRXOc2NhcVqpUx9JYUnQqMyIiEroydsGs/rBzlTluOQBSxkB4lLW5xCcqMyIiEpo2LYIPB8GxgxAZB9e/bN4oUoKOyoyIiISW3BxYOgq+nmSOqzSF3qlQobalsaT4VGZERCR0HNxu3ul691pzfPkg6DgawiKtzSXnRWVGRERCw08LYN5gOJ4BUeWg+2tQv6vVqcQPVGZERKR0y82GtKfgmynmuFoL6PUWlK9hbS7xG5UZEREpvQ5sgZn9YM8Gc9xqKFzzNIRFWBpL/EtlRkRESqcf58L8+yA7E6LLQ48pcHEnq1NJAKjMiIhI6eI+Dosfg7XTzHHyFdBrGpRLsjaXBIzKjIiIlB77f4WZd8IfP5jjNiOg/WPgCrc0lgSW08oXX7FiBd26daNq1ao4HA4+/PDDfJ83DIORI0dStWpVoqOjadeuHT/++KM1YUVExN6+nwlT25pFJqYi3DYbOjytIhMCLC0zR48e5ZJLLmHSpEmFfv65555jwoQJTJo0iTVr1lC5cmU6duzI4cOHSzipiIjYVk4WzL8X5gyAnCNQow0M+hLqdLA6mZQQS5eZOnfuTOfOnQv9nGEYTJw4kccff5yePXsC8Pbbb5OYmMiMGTO45557SjKqiIjY0b5N5rJS+kbAAVc/CG0fBpfOoggltp3trVu3snfvXlJSUvIei4yMpG3btnz11VdnLDPZ2dlkZ2fnjTMzMwFwu9243W6/Zjy5PX9vV4pH82Evmg97KY3z4fj+PVyLHsLhzsIok4Cn+2SMWm3Ba4DX3t9naZwPf/PlZ2PbMrN3714AEhMT8z2emJjI9u3bz/h148aNY9SoUQUeX7JkCTExMf4NeUJaWlpAtivFo/mwF82HvZSG+XB5smmy622qH1gJwL6yDVhXcxDZPx2Fnz6xOJ1vSsN8BEpWVlaRn2vbMnOSw+HINzYMo8Bjp3r00UcZMWJE3jgzM5Pk5GRSUlKIi4vzaza3201aWhodO3YkPFwnmFlN82Evmg97KTXzkf4TYXP74ziwGcPhxHvVQ8S3/gfXOF1WJ/NJqZmPADq5slIUti0zlStXBswjNFWqVMl7PD09vcDRmlNFRkYSGVnwhmHh4eEB+4UJ5LbFd5oPe9F82EvQzodhwPp/w8KHIPc4lK2Mo9c0XDXbEFw1Jr+gnY8S4MvPxdJ3M51NrVq1qFy5cr5DcDk5OSxfvpwrr7zSwmQiIlKisg/DnIHw0X1mkbnwGhi0Emq2sTqZ2ISlR2aOHDnCr7/+mjfeunUrGzZsoEKFClSvXp3hw4czduxYLrroIi666CLGjh1LTEwMffr0sTC1iIiUmD3fw6x+8Oev4HDB/z0BrYeD07b/Ly4WsLTMrF27lvbt2+eNT57rcscddzB9+nQeeughjh07xuDBgzl48CCXX345S5YsITY21qrIIiJSEgzDvB3BosfAkw1x1eDGaVCjldXJxIYsLTPt2rXDMIwzft7hcDBy5EhGjhxZcqFERMRaxzPgo2HmjSIBLroWbpgCMRWszSW2ZdsTgEVEJAT9/q15EbyD28AZBh1GwhVDtKwkZ6UyIyIi1jMMWD0VljwBnhwoVx16vQXJLa1OJkFAZUZERKx17CDMGwo/LzDH9bpC90kQXd7aXBI0VGZERMQ6u9aa71Y6tAOc4ZDyDFx+D5zl4qgip1OZERGRkmcY8PWr8OnT4M2F8jWhVypUa251MglCKjMiIlKysg7Ah3+HzYvMcYPucP0rEFXO2lwStFRmRESk5OxYBbP6Q+YucEVCp7HQor+WleS8qMyIiEjgeb3w5UT47BkwPFDhQug9Hao0sTqZlAIqMyIiElhH98Pce+DXT81x497Q9UWI1NXcxT9UZkREJHC2fQmz+8PhPRAWBZ2fg+a3a1lJ/EplRkRE/M/rgS8mwLKxYHihUl1zWSmxodXJpBRSmREREf86kg6zB8DW5eb4kj5w3QsQUcbaXFJqqcyIiIj/bFkGswfC0XQIj4Hr/gVN+1idSko5lRkRETl/Xg8sHw/LnwMMuKC+uayUUM/qZBICVGZEROT8ZO6BOQNh2xfmuPnt0Gk8RMRYm0tChsqMiIgU36+fwpx7IGs/RJSFrhOhSW+rU0mIUZkRERHfeXLh8zGwcoI5TmxsLitVqmNpLAlNKjMiIuKbjN3mtWN2fG2OW/SHa8dCeJS1uSRkqcyIiEjRbV4McwfBsQMQEQvXvwyNelqdSkKcyoyIiJybxw1LR8FXr5jjKk2hdypUqG1pLBFQmRERkXM5tANm9oPda83xZfdAyj8hLNLaXCInqMyIiMiZ/bQA5g2G4xkQVQ66vwr1u1mdSiQflRkRESkoNwfSnoJvJpvjapdCr1QoX8PaXCKFUJkREZH8DmyFWf3g92/NcauhcM3TEBZhbS6RM1CZERGRv/z4Icy/F7IzIbo89JgMF3e2OpXIWanMiIgIuI/DksdhzZvmOPly6PUWlEuyNpdIEajMiIiEuj9/g5l3wN4fzHGbf0D7x8EVbm0ukSJSmRERCWU/zIKPhkHOEYipCDdMhYs6WJ1KxCcqMyIioch9DBY+DOvfNsc1WsONb0JcVWtziRSDyoyISKjZtxlm3gnpPwIOuPoBaPsIuLRLkOCk31wRkRDi+P59WPQguLOgTAL0nAoXtrc6lsh5UZkREQkFOUdptv0Nwr79whzXuhp6vgmxidbmEvEDlRkRkdIu/SfCPriD6gc2YTicONo9ClfdD06X1clE/EJlRkSktDIM+PYd+OQhHLnHOB4WT9gt0wmro2UlKV1UZkRESqPsw7BgBPzwAQDe2u35POZGOtRoY3EwEf9zWh1ARET8bO8PMLWdWWQcLrjmaTy3vE9OeJzVyUQCQkdmRERKC8OAdamw8BHwZENcNbhxGtRoBW631elEAkZlRkSkNDieCR/dBz/ONccXXWveJLJMRWtziZQAlRkRkWD3+wbzIngHt4IzDDqMhCuGgFNnEkhoUJkREQlWhgGr3zDvdu3JgXLVzTtdJ7e0OplIiVKZEREJRscOwfyh8NNH5rheV+g+CaLLWxpLxAoqMyIiwWbXOph1JxzaAc5wSHkGLr8HHA6rk4lYwtYLqrm5uTzxxBPUqlWL6OhoateuzejRo/F6vVZHExEpeYYBX78Kb11rFpnyNaH/ErhikIqMhLRilZl33nmH1q1bU7VqVbZv3w7AxIkTmTdvnl/DjR8/nilTpjBp0iR++uknnnvuOZ5//nleeeUVv76OiIjtZR2A/94Kix8DrxsadId7VkC15lYnE7Gcz2Vm8uTJjBgxgi5dunDo0CE8Hg8A8fHxTJw40a/hvv76a7p37851111HzZo16dWrFykpKaxdu9avryMiYms7voEpV8HmheCKhOv+Bb3fhqhyVicTsQWfz5l55ZVXeOONN+jRowfPPvts3uMtWrTggQce8Gu4Nm3aMGXKFDZv3kzdunX57rvvWLly5VlLU3Z2NtnZ2XnjzMxMANxuN24/XzTq5Pb8vV0pHs2HvWg+/MDw4lw1CefnY3AYHowKtcm9YRpUbgy5uT5tSvNhL5qPc/PlZ+Nzmdm6dSvNmjUr8HhkZCRHjx71dXNn9fDDD5ORkUG9evVwuVx4PB7GjBnDrbfeesavGTduHKNGjSrw+JIlS4iJifFrvpPS0tICsl0pHs2HvWg+iifCnUnzHVNJzPwegF3lr+C7pH7krt8J7Cz2djUf9qL5OLOsrKwiP9fnMlOrVi02bNhAjRo18j2+cOFCGjRo4Ovmzur999/nP//5DzNmzKBhw4Zs2LCB4cOHU7VqVe64445Cv+bRRx9lxIgReePMzEySk5NJSUkhLs6/9yVxu92kpaXRsWNHwsPD/bpt8Z3mw140H8Xn2PE1rrkP4TiyFyMsCk/KOBKb3kbKeZzkq/mwF83HuZ1cWSkKn8vMgw8+yJAhQzh+/DiGYbB69Wr++9//Mm7cON58801fN3fO13rkkUe45ZZbAGjcuDHbt29n3LhxZywzkZGRREZGFng8PDw8YL8wgdy2+E7zYS+aDx94vbDyX/D5WDC8UKkujt7TCUts6LeX0HzYi+bjzHz5ufhcZvr160dubi4PPfQQWVlZ9OnTh2rVqvHSSy/llQ5/ycrKwnna5bhdLpfemi0ipc+RdJhzN2z53Bxfcit0eQEiy1qbSyQIFOuieQMHDmTgwIHs378fr9dLQkKCv3MB0K1bN8aMGUP16tVp2LAh3377LRMmTOCuu+4KyOuJiFhiy3KYMxCO/AHhMWaJadbX6lQiQeO8rgBcqVIlf+Uo1CuvvMKTTz7J4MGDSU9Pp2rVqtxzzz089dRTAX1dEZES4fXA8udg+XjAgAvqQ+/pkFDP6mQiQaVYJwA7znIS2pYtW84r0KliY2OZOHGi369fIyJiucN7YfYA2PaFOW72N+j8HEQE5l2XIqWZz2Vm+PDh+cZut5tvv/2WRYsW8eCDD/orl4hI6fXrUvP8mKz9EF4Guk2EJjdZnUokaPlcZoYNG1bo46+++qquzCsicjaeXFg2Fr6YABiQ2NhcVqpUx+pkIkHNbzea7Ny5M7Nnz/bX5kRESpeM3fB2N/jiX4ABLe6CAZ+qyIj4wXmdAHyqWbNmUaFCBX9tTkSk9Ni8BObeA8cOQEQsXP8yNOppdSqRUsPnMtOsWbN8JwAbhsHevXvZt28fr732ml/DiYgENY8blo6Gr142x1UuMZeVKtS2NJZIaeNzmenRo0e+sdPp5IILLqBdu3bUq6e3E4qIAHBoJ8y6C3atNseX3QMp/4SwglcoF5Hz43OZefrppwORQ0Sk9Pj5E/jw73D8EESWg+6ToMH1VqcSKbWKVGZ8udmTv2/mKCISNHJz4NOnYdWJJfdql0Kvt6B8TUtjiZR2RSoz8fHxZ71QHpjnzjgcDjwej1+CiYgElYPbYGY/+H29OW41FK55GsIiLI0lEgqKVGY+//zzQOcQEQleG+fDvKGQnQHR5aHHZLi4s9WpREJGkcpM27ZtA51DRCT4uI9D2pOweqo5Tr4cbpwG8cnW5hIJMcW+zkxWVhY7duwgJycn3+NNmjQ571AiIrb3528w807Y+705bj0c/u8JcIVbmUokJPlcZvbt20e/fv1YuHBhoZ/XOTMiUur9bzbMHwY5hyGmItwwFS7qYHUqkZDl8+0Mhg8fzsGDB1m1ahXR0dEsWrSIt99+m4suuoj58+cHIqOIiD24j8FHw8zrx+QchhqtYdBKFRkRi/l8ZOazzz5j3rx5tGzZEqfTSY0aNejYsSNxcXGMGzeO6667LhA5RUSstf8Xc1npj/8BDrj6AWj7CLj8dlcYESkmn/8WHj16lISEBAAqVKjAvn37qFu3Lo0bN2b9+vV+DygiYrnv3ocF/wD3USiTAD2nwoXtrU4lIif4vMx08cUXs2nTJgCaNm3K66+/zu7du5kyZQpVqlTxe0AREcvkZMG8ITD3brPI1LraXFZSkRGxFZ+PzAwfPpw9e/YA5q0Nrr32Wt59910iIiKYPn26v/OJiFgj/SdzWWnfz+BwmktKVz8ATpfVyUTkNEUuMz169GDAgAHceuutOJ3mAZ1mzZqxbds2fv75Z6pXr06lSpUCFlREpEQYBmx4Fz5+AHKPQdnKcOObUOsqq5OJyBkUeZnp2LFj9OjRg6SkJB577DF++eUXAGJiYmjevLmKjIgEv+wjMHeQubSUewwu/D9zWUlFRsTWilxmFi9ezLZt2/j73//OBx98QL169bj66qv597//zbFjxwKZUUQk8Pb+D6a2g+/fA4cLrnkK+s6GshdYnUxEzsGnE4CTkpJ48skn+fXXX/n000+pUaMGgwcPpnLlytxzzz188803gcopIhIYhgFrU+HNa+DPXyC2Ktz5MVx1Pzh9fo+EiFig2BdIaN++Pe3bt+fw4cPMmDGDxx57jGnTppGbm+vPfCIigXM8ExYMN6/oC3DRteZNIstUtDSWiPjmvK72tGXLFqZPn8706dPJyMigQwddBVNEgsSe78x3Kx3YAs4wuOZpaDVUR2NEgpDPZebYsWPMnDmT1NRUVqxYQfXq1RkwYAD9+vUjOVl3ihURmzMMWPMmLH4MPDlQLhl6pUJyS6uTiUgxFbnMfPXVV6SmpvLBBx+Qk5NDjx49WLx4sY7GiEjwOHYI5t8LP524j9zF10H3SRBTwdJYInJ+ilxm2rRpwyWXXMKYMWPo27cv5cuXD2QuERH/2r0OZvaDQ9vBGQ4p/4TLB4HDYXUyETlPRS4za9eupXnz5oHMIiLif4YBqyZD2lPgdUN8DeidCtUutTqZiPhJkcuMioyIBJ2sA+YF8DZ9Yo7rXw/XvwLR8ZbGEhH/0r3rRaR02rkaZt0FGTvBFQHXjoWWA7SsJFIKqcyISOni9cLXr8DS0eDNhQq1ofd0qHKJ1clEJEBUZkSk9Dj6J3w4CH5ZYo4b3QhdJ0JUnKWxRCSwVGZEpHTY/hXM6g+Hf4ewKOg8HprfoWUlkRBQpDLTrFkzHEX8B2H9+vXnFUhExCdeL6ycAJ+PBcMDFS8yl5UqN7I6mYiUkCKVmR49egQ4hohIMRzZB3MGwpbPzXGTW+C6f0FkWWtziUiJKlKZefrppwOdQ0TEN1tXwOwBcOQPCIuG616Apn21rCQSgnTOjIgEF68HVjwPy8eD4YUL6kHvtyGhntXJRMQiPpcZj8fDiy++yAcffMCOHTvIycnJ9/kDBw74LZyISD6H95rLSltXmONmt0Hn5yEixtpcImIpn+91P2rUKCZMmMBNN91ERkYGI0aMoGfPnjidTkaOHBmAiCIiwG+fwZQ2ZpEJLwM3TIXur6rIiIjvZebdd9/ljTfe4IEHHiAsLIxbb72VN998k6eeeopVq1YFIqOIhDJPLiz9J7zTE47ug8RGcPcyuORmq5OJiE34XGb27t1L48aNAShbtiwZGRkAdO3alY8//ti/6UQktGXshre7wRcvAAZc2g8GfAoX1LU6mYjYiM9lJikpiT179gBQp04dliwxr7S5Zs0aIiMj/ZtORELXL2nmstKOryAiFm6cBt0mQni01clExGZ8LjM33HADS5cuBWDYsGE8+eSTXHTRRdx+++3cddddfg+4e/dubrvtNipWrEhMTAxNmzZl3bp1fn8dEbEJjxvSnoJ3e8GxA1C5CdyzHBr3sjqZiNiUz+9mevbZZ/M+7tWrF0lJSXz11VfUqVOH66+/3q/hDh48SOvWrWnfvj0LFy4kISGB3377jfj4eL++jojYxKGd5p2ud602x5fdDR3/CeFR1uYSEVs77+vMXHHFFVxxxRX+yFLA+PHjSU5OJjU1Ne+xmjVrnvVrsrOzyc7OzhtnZmYC4Ha7cbvdfs13cnv+3q4Uj+bDXnydD8fmhbg+uhfH8UMYkXF4ur6EUa/byY0FKmbI0N8Pe9F8nJsvPxuHYRiGry+we/duvvzyS9LT0/F6vfk+d9999/m6uTNq0KAB1157Lbt27WL58uVUq1aNwYMHM3DgwDN+zciRIxk1alSBx2fMmEFMjN7CKWI3Dm8uDX5/nzr7FgNwMKY2a2sOJisyweJkImKlrKws+vTpQ0ZGBnFxcWd9rs9lJjU1lUGDBhEREUHFihXz3YDS4XCwZcuW4qUuRFSUeWh5xIgR9O7dm9WrVzN8+HBef/11br/99kK/prAjM8nJyezfv/+cPwxfud1u0tLS6NixI+Hh4X7dtvhO82EvRZqPQ9txzRmAc8+3AHguuwfv/z0NrogSTBoa9PfDXjQf55aZmUmlSpWKVGZ8XmZ66qmneOqpp3j00UdxOn0+f9gnXq+XFi1aMHbsWMC8e/ePP/7I5MmTz1hmIiMjC31XVXh4eMB+YQK5bfGd5sNezjgfG+fDvKGQnQFR8dBjMq56XXCVeMLQor8f9qL5ODNffi4+t5GsrCxuueWWgBcZgCpVqtCgQYN8j9WvX58dO3YE/LVFJEBys+GTB+GDv5lFJqklDPoC6nWxOpmIBCmfG0n//v2ZOXNmILIU0Lp1azZt2pTvsc2bN1OjRo0SeX0R8bM/f4NpHWH1VHPcehj0Wwjx1a3NJSJBzedlpnHjxtG1a1cWLVpE48aNCxwGmjBhgt/C/eMf/+DKK69k7Nix3HTTTaxevZqpU6cydepUv72GiJSQ/82G+cMg5zBEV4AbXoe6KVanEpFSwOcyM3bsWBYvXszFF18MUOAEYH9q2bIlc+fO5dFHH2X06NHUqlWLiRMn0rdvX7++jogEkPsYLHoQ1p24xEL1VubVfMtVszaXiJQaPpeZCRMm8NZbb3HnnXcGIE5BXbt2pWvXriXyWiLiX2WP7yFseidI/xFwwFUjoN1j4DrvS1yJiOTx+V+UyMhIWrduHYgsIlKKOP43k7abnsLhzYaYStBzKtS5xupYIlIK+XwC8LBhw3jllVcCkUVESoOcLJg3hLB5fyfMm423Rmv4+5cqMiISMD4fmVm9ejWfffYZCxYsoGHDhgVOAJ4zZ47fwolIkEn/GWbeAft+xsDBpsrdubDP6zgjdW8lEQkcn8tMfHw8PXv2DEQWEQlWhgEb3oWPH4DcY1A2EU/3yWzaeIQLnboMnogEls9l5tSbPoqIkH0EPr4fvn/PHNduBz3fwIgsDxs/sTSaiIQGvaVARIpv7/9gVj/YvxkcTmj/GLS5H5xO3elaREpMkcpM8+bNWbp0KeXLl6dZs2ZnvZ7M+vXr/RZORGzKMGDddFj0COQeh9gq5rVjauqdjiJS8opUZrp3755388bu3bv7/eJ4IhJEjmfCguHmFX0B6nSEG6ZAmUqWxhKR0FWkMvP000/nfTxy5MhAZRERu9vzHcy8Ew5sAYcLrnkKrrzPXFYSEbGIz/8C1a5dmz///LPA44cOHaJ27dp+CSUiNmMYsPoNeLODWWTikswbRLYZriIjIpbz+QTgbdu24fF4CjyenZ3Nrl27/BJKRGzk2CH46D7YOM8c1+0MPV6DmAqWxhIROanIZWb+/Pl5Hy9evJhy5crljT0eD0uXLqVWrVr+TSci1tq9Dmb2g0PbwRkOHUfBFYNB582JiI0Uucz06NEDMO+Mfccdd+T7XHh4ODVr1uRf//qXX8OJiEUMA1ZNhrSnwOuG+OrQazokXWp1MhGRAopcZrxeLwC1atVizZo1VKqkdy6IlEpZB2DeUNj0sTmu3w2unwTR8ZbGEhE5E5/Pmdm6dWsgcoiIHexcY14EL2MnuCIgZQxcNlDLSiJia8V6G8LSpUvp2rUrF154IXXq1KFr1658+umn/s4mIiXF64UvX4LUTmaRKV8L+qfB5XeryIiI7flcZiZNmkSnTp2IjY1l2LBh3HfffcTFxdGlSxcmTZoUiIwiEkhH/4T/3nzi/JhcaNgT7lkBVZtanUxEpEh8XmYaN24cL774IkOHDs177L777qN169aMGTMm3+MiYnPbv4JZ/eHw7+CKhM7PwqX9dDRGRIKKz0dmMjMz6dSpU4HHU1JSyMzM9EsoEQkwrxdWvADTu5pFpmIdGLgUWtylIiMiQcfnMnP99dczd+7cAo/PmzePbt26+SWUiATQkX3w7o3w2T/B8ECTm+Hu5VC5sdXJRESKxedlpvr16zNmzBiWLVtGq1atAFi1ahVffvkl999/Py+//HLec++77z7/JRWR87f1C5g9AI7shbBo6PI8NLtNR2NEJKj5XGamTZtG+fLl2bhxIxs3bsx7PD4+nmnTpuWNHQ6HyoyIXXg9sOJ5WD4eDC9cUA96T4eE+lYnExE5b7rOjEhpd/gPmDMAtq4wx01vgy7PQUQZa3OJiPiJz2XmpP379+NwOKhYsaI/84iIP/32OcwZCEf3QXgZ6DoBLrnF6lQiIn7l0wnAhw4dYsiQIVSqVInExEQSEhKoVKkSQ4cO5dChQwGKKCI+8+TC0n/COzeYRSahIdy9TEVGREqlIh+ZOXDgAK1atWL37t307duX+vXrYxgGP/30E9OnT2fp0qV89dVXlC9fPpB5ReRcMn83rx2z4ytzfOmd0OlZCI+2NJaISKAUucyMHj2aiIgIfvvtNxITEwt8LiUlhdGjR/Piiy/6PaSIFNEvaTD3Hsj6EyLKQreXoHEvq1OJiARUkZeZPvzwQ1544YUCRQagcuXKPPfcc4Vef0ZESoDHbd6O4N1eZpGp3Ni8JYGKjIiEgCIfmdmzZw8NGzY84+cbNWrE3r17/RJKRHxwaCfM7g87vzHHLQdCyjMQHmVtLhGRElLkMlOpUiW2bdtGUlJSoZ/funWr3tkkUtI2LYS5g+D4IYiMg+tfgYY9rE4lIlKiirzM1KlTJx5//HFycnIKfC47O5snn3yy0Hs2iUgA5ObA4sfhv7eYRaZqM3NZSUVGREJQkY/MjBo1ihYtWnDRRRcxZMgQ6tWrB8DGjRt57bXXyM7O5p133glYUBE54eA2mHUX7F5njq8YDB1GQViEpbFERKxS5DKTlJTE119/zeDBg3n00UcxDAMwb1vQsWNHJk2aRHJycsCCigiwcT7MGwrZGRBVDnpMhnrXWZ1KRMRSPl0BuFatWixcuJCDBw/yyy+/AFCnTh0qVKgQkHAickJuNix5AlZPNcdJLaHXWxBf3dpcIiI2UKzbGZQvX57LLrvM31lEpDB//gaz+sGe78zxlffBNU+BK9zaXCIiNlHsezOJSAn43xyYfx/kHIboCnDDFKh7rdWpRERsRWVGxI7cx2DRo7Au1RxXbwU3ToNy1azNJSJiQyozInaz/xeYeSf88T/AAVeNgHaPgUt/XUVECqN/HUXs5PsP4KPh4D4KMZWg51Soc43VqUREbE1lRsQOcrJg4YPw7X/Mcc2roOcbEFfF2lwiIkFAZUbEauk/m8tK+34CHND2YWj7EDhdVicTEQkKKjMiVvr2Xfj4fsg9BmUTzaMxtdtanUpEJKgU+d5MdjBu3DgcDgfDhw+3OorI+ck+Yt4gct5gs8jUbgeDVqrIiIgUQ9AcmVmzZg1Tp06lSZMmVkcROT9//GguK+3fDA4ntH8M2ozQspKISDEFRZk5cuQIffv25Y033uCZZ54563Ozs7PJzs7OG2dmZgLgdrtxu91+zXVye/7erhSP7efDMHBseAfXksdw5B7HKFsZzw1TMapfCR6v+acUsf18hBjNh71oPs7Nl5+Nwzh5x0gbu+OOO6hQoQIvvvgi7dq1o2nTpkycOLHQ544cOZJRo0YVeHzGjBnExMQEOKlI4cI8x7hkRypJh1YB8EdsE9bXuJuc8DiLk4mI2FNWVhZ9+vQhIyODuLiz/1tp+yMz7733HuvXr2fNmjVFev6jjz7KiBEj8saZmZkkJyeTkpJyzh+Gr9xuN2lpaXTs2JHwcN0nx2q2nY+93xM2pz+OQ1sxHC687R+nwhVD6eAIqlPWfGbb+QhRmg970Xyc28mVlaKwdZnZuXMnw4YNY8mSJURFRRXpayIjI4mMjCzweHh4eMB+YQK5bfGdbebDMGDNm7D4MfDkQFwSjl5v4ap+OaF0doxt5kMAzYfdaD7OzJefi63LzLp160hPT+fSSy/Ne8zj8bBixQomTZpEdnY2Llco7RYkaBzPgPn3wsZ55rhuZ+jxGsRUsDaXiEgpZOsyc8011/DDDz/ke6xfv37Uq1ePhx9+WEVG7Gn3evPdSoe2gzMMOo6GKwaDw2F1MhGRUsnWZSY2NpZGjRrle6xMmTJUrFixwOMiljMM+GYKLHkSvG6Irw69pkPSpef8UhERKT5blxmRoHHsIMwbCj8vMMf1ukL3VyE63tJYIiKhIOjKzLJly6yOIJLfzjUw6y7I2AGuCEgZA5cN1LKSiEgJCboyI2IbXi98PQmWjgJvLpSvBb1ToWozq5OJiIQUlRmR4sg6YN5b6ZfF5rjhDdDtZYjSRfBEREqayoyIr7Z/DbP7Q+ZucEVC52fh0n5aVhIRsYjKjEhReb3w5Yvw2RgwPFCxDvSeDpUbW51MRCSkqcyIFMWRfTD3bvjtM3Pc+CboOgEiY63NJSIiKjMi57T1C5g9AI7shbBo6PI8NLtNy0oiIjahMiNyJl4PrHgBlj8LhhcqXWwuKyU2sDqZiIicQmVGpDCH/4A5A2DrCnPc9Dbo8hxElLE2l4iIFKAyI3K63z6HOXfD0XQIj4GuL8Ilt1idSkREzkBlRuQkT665pLTiBcCAhIbmstIFda1OJiIiZ6EyIwKQ+bt5ku/2L81x8zug83gIj7Y2l4iInJPKjMgvn5pvu876EyLKQreXoHEvq1OJiEgRqcxI6PK44fMxsPJFc1y5MfR+GypeaG0uERHxicqMhKaMXeadrnd+Y45bDjDvdh0eZW0uERHxmcqMhJ5Ni+DDQXDsIETGwfWvQMMeVqcSEZFiUpmR0JGbA0tHwdeTzHHVZtArFSrUsjaXiIicF5UZCQ0Ht5vLSrvXmuMrBkOHkRAWaWksERE5fyozUvr9tADmDYbjGRBVDnpMhnrXWZ1KRET8RGVGSq/cbEh7Cr6ZYo6TWkKvtyC+urW5RETEr1RmpHQ6sAVm9oM9G8zxlffCNU+DK9zSWCIi4n8qM1LqOH6aBx//A7IzIboC3DAF6l5rdSwREQkQlRkpPXKP02TndMK+/cwcV28FN06DctWszSUiIgGlMiOlw/5fCZt5J7X2/2CO24yA9o+DS7/iIiKlnf6ll+D3/UxYMBxHzhGyw2Jx9ZpGWD0tK4mIhAqVGQleOVmw6GFY/28AvDVasyz2Jv7vwv+zOJiIiJQkp9UBRIpl3yZ485oTRcYBbR/G02cOx8PLW51MRERKmI7MSPDZMAM+vh/cWVAmAW58E2q3Bbfb6mQiImIBlRkJHjlH4eMH4LsZ5rh2O+j5BpRNsDSWiIhYS2VGgsMfP8LMO2H/ZnA4od1jcNUIcLqsTiYiIhZTmRF7MwzzvJiFD0HucYitYi4r1WxjdTIREbEJlRmxr+zDsOAf8MNMc1ynA9zwOpSpZG0uERGxFZUZsac935vLSgd+A4cLrnkSrhwGTr0BT0RE8lOZEXsxDFg7DRY9Bp5siKtm3um6+hVWJxMREZtSmRH7OJ4B8++DjR+a47qdoMdkiKlgaSwREbE3lRmxh93rYVY/OLgNnGHQYRS0GgIOh9XJRETE5lRmxFqGAd+8DkueAK8bylWH3qmQ1MLqZCIiEiRUZsQ6xw7CvKHw8wJzXK8rdJ8E0bolgYiIFJ3KjFhj11qY2Q8ydoArAlKegcvu1rKSiIj4TGVGSpZhwNeT4NOR4M2F8jWh93So2sziYCIiEqxUZqTkZB2AD/8OmxeZ44Y3QLeXIKqctblERCSoqcxIydixCmbdBZm7wRUJncZBi7u0rCQiIufN1pdTHTduHC1btiQ2NpaEhAR69OjBpk2brI4lvvB64YsJkNrFLDIVLoQBn0LL/ioyIiLiF7YuM8uXL2fIkCGsWrWKtLQ0cnNzSUlJ4ejRo1ZHk6I4uh9m9Ialo8DwQOPecM9yqNLE6mQiIlKK2HqZadGiRfnGqampJCQksG7dOq6++mqLUkmRbFsJswfA4T0QFgVdnodmf9PRGBER8Ttbl5nTZWRkAFChwpkvb5+dnU12dnbeODMzEwC3243b7fZrnpPb8/d2g5rXg/PLF3F+8RwOw4tRqS65N0yDhPqQmxvQl9Z82Ivmw140H/ai+Tg3X342DsMwjABm8RvDMOjevTsHDx7kiy++OOPzRo4cyahRowo8PmPGDGJiYgIZMeRFug9x6bYpXHBkIwA7KrTh+6Q78LgiLU4mIiLBJisriz59+pCRkUFcXNxZnxs0ZWbIkCF8/PHHrFy5kqSkpDM+r7AjM8nJyezfv/+cPwxfud1u0tLS6NixI+Hh4X7ddrBxbF2Ba94gHEfTMcJj8HR6DqPJLSWaQfNhL5oPe9F82Ivm49wyMzOpVKlSkcpMUCwz3XvvvcyfP58VK1actcgAREZGEhlZ8EhAeHh4wH5hArlt2/N6YNmzsOJ5wICEBjh6TyfsgostixTS82FDmg970XzYi+bjzHz5udi6zBiGwb333svcuXNZtmwZtWrVsjqSnCpzj3mS7/aV5rj57dBpPERoOU9EREqOrcvMkCFDmDFjBvPmzSM2Npa9e/cCUK5cOaKjoy1OF+J+/RTm3A1Zf0JEWeg6EZr0tjqViIiEIFuXmcmTJwPQrl27fI+npqZy5513lnwgAU8ufP4MrHzRHCc2Nu+tVKmOpbFERCR02brMBMm5yaEjYxfM6g87V5njFv3h2rEQHmVtLhERCWm2LjNiI5sWwYeD4NhBiIwzbxDZqKfVqURERFRm5Bxyc8zbEXw9yRxXaQq9U6FCbUtjiYiInKQyI2d2cLt5p+vda83x5YOg42gI00XwRETEPlRmpHA/LYB5g+F4BkSVg+6vQf2uVqcSEREpQGVG8svNhrSn4Jsp5rhaC+j1FpSvYW0uERGRM1CZkb8c2AIz+8GeDea41VC45mkIi7A0loiIyNmozIjpx7kw/z7IzoTo8tBjClzcyepUIiIi56QyE+rcx2HxY7B2mjlOvgJ6TYNyZ78HloiIiF2ozISy/b/CzDvhjx/McZt/QPvHwaWbnomISPBQmQlV38+EBcMh5wjEVIQbpsJFHaxOJSIi4jOVmVCTkwWLHob1/zbHNVrDjW9CXFVrc4mIiBSTykwo2bfJXFZK3wg44OoHoe3D4NKvgYiIBC/txULFhhnw8f3gzoIyCdBzKlzY3upUIiIi501lprTLOQofPwDfzTDHtdpCzzcgNtHaXCIiIn6iMlOa/bHRXFbavwkcTmj3KFx1PzhdVicTERHxG5WZ0sgw4Nt34JMHIfc4lK1sXjumZhurk4mIiPidykxpk30YFoyAHz4wxxdeAze8DmUvsDaXiIhIgKjMlCZ7fzCXlf78FRwu+L8noPVwcDqtTiYiIhIwKjOlgWHA2rdg0aPgyYa4anDjNKjRyupkIiIiAacyE+yOZ8BHw8wbRQLU7QQ9JkNMBWtziYiIlBCVmWD2+7cwsx8c3ArOMOgwEloNBYfD6mQiIiIlRmUmGBkGrJ4KS54ATw6Uqw693oLkllYnExERKXEqM8Hm2EGYNxR+XmCO63WF7pMgury1uURERCyiMhNMdq2DWXfCoR3gDIeUZ+Dye7SsJCIiIU1lJhgYBnz9Knz6NHhzoXxN6JUK1ZpbnUxERMRyKjN2l3UAPhwMmxea4wY94PqXIaqcpbFERETsQmXGznZ8A7Pugsxd4IqETmOhRX8tK4mIiJxCZcaOvF746iVY+k8wPFDhQug9Hao0sTqZiIiI7ajM2M3R/TB3EPyaZo4b94auL0JkrLW5REREbEplxk62fQmz+8PhPRAWBZ2fg+a3a1lJRETkLFRm7MDrgS8mwLKxYHihUl1zWSmxodXJREREbE9lxmpH0mHOQNiyzBxf0geuewEiylgaS0REJFiozFhpy3KYPQCOpkN4DFz3L2jax+pUIiIiQUVlxgpeDywfD8ufAwxIaGBeBC+hntXJREREgo7KTEnL3GMuK237whw3vx06jYeIGGtziYiIBCmVmZL061KYczdk7YeIstB1IjTpbXUqERGRoKYyUxI8ufD5GFg5wRwnNjbfrVSpjqWxREQktHm9Bm6vF4/XwO0xyPWc+NhrfpzrNcj1GLhPPJ7r9eL2GCee7yXXY5DrNah9QRnqJlp3PTSVmUDL2G1eO2bH1+a4RX+4diyER1mbS0REisRfO/xc718fZ+e42bDXwf5VOzBwnNjGKV934rln3ob5X7fXwHPi9fJynfz8iVx52zgtb67Xi9fwz89oSPsLefBa6877VJkJpM1LYO49cOwARMZBt5egUU+rU4mI+J33lB1wvp3m6TvpE497irDDz/V4zZ31ia8xn+8tsLPOv4Mu3g7/1Lx5pcRjfk+Gn3b4Bblg68+B2nixOR0Q5nIS7nTgcjoIdzkJczkIc5787ykfu5yEOR1UKRdtaWaVmUDwuGHpaPjqZXNcpSn0ToUKtS2NJSLW8py6oz1lZ+o+ZQdb2A7/1B3zqTtg92lHCE49cpC3Uz/x3Bx3Llu2OVn54Y94DP7aWRd6lCH//70XVlBKbodvT66TO3rnXzv0c+3ww04UA5fTgcsB+/f9QVLVKkSEuQpu48R2w11/FQrXadsIz/d6p7zuKUXk9G3kz/XXx+HOv7bvdAbfVedVZvzt0A7zTte71pjjywdBx9EQFmltLpEgUZQdfqE72BM7VY+n4GH+c+7wCxwJKOQwv/f01zvzDt88EvDXEYeTywbW7/Cd8MfuEnu1ouzw83bShR4JMJ/rcp1hG2fYWRe2jXw7bpcDl9PaHb7b7eaTTz6hS5dLCA8P99NPPHSpzPjTzx/Dh4Ph+CGIKgfdX4X63axOJaWQp5D/Mz7X4fpz7fDz7awL7JhPHK4/eZi/CDv8nFwP+/908eaOVeR6+WsbJ17/9PMNTpYW63f4JevkzvH0HXOY87Qd7Gk74MJ31mfahhMnBlt++4UG9S4mMjzstGJw9h1+WKFHAk4tIid29KXg//AlOAVFmXnttdd4/vnn2bNnDw0bNmTixIlcddVVVsf6iycHPn0Svplsjqtdal4Er3wNa3OFMMMw8Br8dbi+iDv8giffeYt0JOBMO/wz/d974UXk1BMLS8sO3wFHMs97K4UdXj/XDj//EYC/dtynH6L/a3vOAjvrU3f4eTvpQraR99xTikGB5YHTdvhhLgcuR8nt8N1uN5/kbKZL29o6EiClju3LzPvvv8/w4cN57bXXaN26Na+//jqdO3dm48aNVK9e3ep4xGSn43q7C+zZYD7Qaihc8zSERViaqyjOtMM/5/r8GU+cO/cOv/CT7/KXi9OXB862wz+9iOS4XTywOg23J2j29n5zpvX0vI8LPUR/6lp8wR3wGbdxhiMBeTtppwOH4eX7Dd9yWcsWREWEn3OHH3baNk59DYfuHC8iZ2H7MjNhwgT69+/PgAEDAJg4cSKLFy9m8uTJjBs3zrJch4+7ObZhNm1/fhKn9xieyHh2tZ3AgaT/w7PryFl3+Cd31oUdCcgtxg4/9wzbOLVgnDxacHoxKH0cwJm/r+Ic0vf1//B92eGHn14wClmfL1oRsd8O3+12Y+wwaH/xBToSICIBZesyk5OTw7p163jkkUfyPZ6SksJXX31V6NdkZ2eTnZ2dN87MNA9xu91u3G6337L99J8HuWxXKgBrvHW5L+Ne9syPAgrPFUzy7/AL7oBP3fEXOPx+po8L+9pTdsqn76QLlIF8xeD055ofG14Pq77+knZXX0V0ZESBolE61/C94IVcr9U5Cjr5982ff++k+DQf9qL5ODdffja2LjP79+/H4/GQmJiY7/HExET27t1b6NeMGzeOUaNGFXh8yZIlxMT47/5HfxyJp4Xh4A1vV1719sIbHkY5DFxO8z36Lsdf/z31Y6fDKOSxgs/L/7Hx18dOcGL+99zbMIqw7fwfO06MA8p74k8RGYD7xJ9zqRQF/1v9RfFySUCkpaVZHUFOofmwF83HmWVlZRX5ubYuMyedfvjcMIwzHlJ/9NFHGTFiRN44MzOT5ORkUlJSiIuL82OqLmTv7Unl9VtZ1bGjDqPbgNvtJi0tjY6aD1vQfNiL5sNeNB/ndnJlpShsXWYqVaqEy+UqcBQmPT29wNGakyIjI4mMLHhNl/DwcP//wlSuD2wNzLal2DQf9qL5sBfNh71oPs7Ml5+LM4A5zltERASXXnppgcNwaWlpXHnllRalEhERETux9ZEZgBEjRvC3v/2NFi1a0KpVK6ZOncqOHTsYNGiQ1dFERETEBmxfZm6++Wb+/PNPRo8ezZ49e2jUqBGffPIJNWrognQiIiISBGUGYPDgwQwePNjqGCIiImJDtj5nRkRERORcVGZEREQkqKnMiIiISFBTmREREZGgpjIjIiIiQU1lRkRERIKayoyIiIgENZUZERERCWoqMyIiIhLUguIKwOfDMAzAt1uJF5Xb7SYrK4vMzEzd9dQGNB/2ovmwF82HvWg+zu3kfvvkfvxsSn2ZOXz4MADJyckWJxERERFfHT58mHLlyp31OQ6jKJUniHm9Xn7//XdiY2NxOBx+3XZmZibJycns3LmTuLg4v25bfKf5sBfNh71oPuxF83FuhmFw+PBhqlatitN59rNiSv2RGafTSVJSUkBfIy4uTr+MNqL5sBfNh71oPuxF83F25zoic5JOABYREZGgpjIjIiIiQU1l5jxERkby9NNPExkZaXUUQfNhN5oPe9F82Ivmw79K/QnAIiIiUrrpyIyIiIgENZUZERERCWoqMyIiIhLUVGZEREQkqKnMFNNrr71GrVq1iIqK4tJLL+WLL76wOlJIGjduHC1btiQ2NpaEhAR69OjBpk2brI4lJ4wbNw6Hw8Hw4cOtjhLSdu/ezW233UbFihWJiYmhadOmrFu3zupYISk3N5cnnniCWrVqER0dTe3atRk9ejRer9fqaEFNZaYY3n//fYYPH87jjz/Ot99+y1VXXUXnzp3ZsWOH1dFCzvLlyxkyZAirVq0iLS2N3NxcUlJSOHr0qNXRQt6aNWuYOnUqTZo0sTpKSDt48CCtW7cmPDychQsXsnHjRv71r38RHx9vdbSQNH78eKZMmcKkSZP46aefeO6553j++ed55ZVXrI4W1PTW7GK4/PLLad68OZMnT857rH79+vTo0YNx48ZZmEz27dtHQkICy5cv5+qrr7Y6Tsg6cuQIzZs357XXXuOZZ56hadOmTJw40epYIemRRx7hyy+/1NFjm+jatSuJiYlMmzYt77Ebb7yRmJgY3nnnHQuTBTcdmfFRTk4O69atIyUlJd/jKSkpfPXVVxalkpMyMjIAqFChgsVJQtuQIUO47rrr6NChg9VRQt78+fNp0aIFvXv3JiEhgWbNmvHGG29YHStktWnThqVLl7J582YAvvvuO1auXEmXLl0sThbcSv2NJv1t//79eDweEhMT8z2emJjI3r17LUolYN5hdcSIEbRp04ZGjRpZHSdkvffee6xfv541a9ZYHUWALVu2MHnyZEaMGMFjjz3G6tWrue+++4iMjOT222+3Ol7Iefjhh8nIyKBevXq4XC48Hg9jxozh1ltvtTpaUFOZKSaHw5FvbBhGgcekZA0dOpTvv/+elStXWh0lZO3cuZNhw4axZMkSoqKirI4jgNfrpUWLFowdOxaAZs2a8eOPPzJ58mSVGQu8//77/Oc//2HGjBk0bNiQDRs2MHz4cKpWrcodd9xhdbygpTLjo0qVKuFyuQochUlPTy9wtEZKzr333sv8+fNZsWIFSUlJVscJWevWrSM9PZ1LL7007zGPx8OKFSuYNGkS2dnZuFwuCxOGnipVqtCgQYN8j9WvX5/Zs2dblCi0PfjggzzyyCPccsstADRu3Jjt27czbtw4lZnzoHNmfBQREcGll15KWlpavsfT0tK48sorLUoVugzDYOjQocyZM4fPPvuMWrVqWR0ppF1zzTX88MMPbNiwIe9PixYt6Nu3Lxs2bFCRsUDr1q0LXK5g8+bN1KhRw6JEoS0rKwunM/+u1+Vy6a3Z50lHZophxIgR/O1vf6NFixa0atWKqVOnsmPHDgYNGmR1tJAzZMgQZsyYwbx584iNjc07YlauXDmio6MtThd6YmNjC5yvVKZMGSpWrKjzmCzyj3/8gyuvvJKxY8dy0003sXr1aqZOncrUqVOtjhaSunXrxpgxY6hevToNGzbk22+/ZcKECdx1111WRwtuhhTLq6++atSoUcOIiIgwmjdvbixfvtzqSCEJKPRPamqq1dHkhLZt2xrDhg2zOkZI++ijj4xGjRoZkZGRRr169YypU6daHSlkZWZmGsOGDTOqV69uREVFGbVr1zYef/xxIzs72+poQU3XmREREZGgpnNmREREJKipzIiIiEhQU5kRERGRoKYyIyIiIkFNZUZERESCmsqMiIiIBDWVGREREQlqKjMiIiIS1FRmRMSWli1bhsPh4NChQ1ZHERGb0xWARcQW2rVrR9OmTZk4cSIAOTk5HDhwgMTERBwOh7XhRMTWdKNJEbGliIgIKleubHUMEQkCWmYSEcvdeeedLF++nJdeegmHw4HD4WD69On5lpmmT59OfHw8CxYs4OKLLyYmJoZevXpx9OhR3n77bWrWrEn58uW599578Xg8edvOycnhoYceolq1apQpU4bLL7+cZcuWWfONikhA6MiMiFjupZdeYvPmzTRq1IjRo0cD8OOPPxZ4XlZWFi+//DLvvfcehw8fpmfPnvTs2ZP4+Hg++eQTtmzZwo033kibNm24+eabAejXrx/btm3jvffeo2rVqsydO5dOnTrxww8/cNFFF5Xo9ykigaEyIyKWK1euHBEREcTExOQtLf38888Fnud2u5k8eTIXXnghAL169eKdd97hjz/+oGzZsjRo0ID27dvz+eefc/PNN/Pbb7/x3//+l127dlG1alUAHnjgARYtWkRqaipjx44tuW9SRAJGZUZEgkZMTExekQFITEykZs2alC1bNt9j6enpAKxfvx7DMKhbt26+7WRnZ1OxYsWSCS0iAacyIyJBIzw8PN/Y4XAU+pjX6wXA6/XicrlYt24dLpcr3/NOLUAiEtxUZkTEFiIiIvKduOsPzZo1w+PxkJ6ezlVXXeXXbYuIfejdTCJiCzVr1uSbb75h27Zt7N+/P+/oyvmoW7cuffv25fbbb2fOnDls3bqVNWvWMH78eD755BM/pBYRO1CZERFbeOCBB3C5XDRo0IALLriAHTt2+GW7qamp3H777dx///1cfPHFXH/99XzzzTckJyf7ZfsiYj1dAVhERESCmo7MiIiISFBTmREREZGgpjIjIiIiQU1lRkRERIKayoyIiIgENZUZERERCWoqMyIiIhLUVGZEREQkqKnMiIiISFBTmREREZGgpjIjIiIiQe3/AZcH96p2XDQ6AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "m.solution.to_dataframe().plot(grid=True, ylabel=\"Optimal Value\");" + ] + }, + { + "cell_type": "markdown", + "id": "603aa77f-4a13-4e02-bf0c-196e3cdfcbf0", + "metadata": {}, + "source": [ + "Alright! Now you learned how to set up linopy variables and expressions with coordinates. For more advanced `linopy` operations you can check out the [User Guide](https://linopy.readthedocs.io/en/latest/user-guide.html)." + ] + }, + { + "cell_type": "markdown", + "id": "f520ffce-6783-4c89-81f8-cc0c7822439c", + "metadata": {}, + "source": [ + "## Electricity Market Examples" + ] + }, + { + "cell_type": "markdown", + "id": "b1bc668c-c37c-4939-8890-d4e945ba2739", + "metadata": {}, + "source": [ + "### Single bidding zone, single period" + ] + }, + { + "cell_type": "markdown", + "id": "a2459e7a-d740-45dd-a834-bc098a5825cb", + "metadata": {}, + "source": [ + "We want to minimise operational cost of an example electricity system representing South Africa subject to generator limits and meeting the load:\n", + "\n", + "\\begin{equation}\n", + " \\min_{g_s} \\sum_s o_s g_s\n", + " \\end{equation}\n", + " such that\n", + " \\begin{align}\n", + " g_s &\\leq G_s \\\\\n", + " g_s &\\geq 0 \\\\\n", + " \\sum_s g_s &= d\n", + " \\end{align}" + ] + }, + { + "cell_type": "markdown", + "id": "b5c0cb03-aa1b-47f0-92d7-b98fea75bc4d", + "metadata": {}, + "source": [ + "We are given the following information on the South African electricity system:" + ] + }, + { + "cell_type": "markdown", + "id": "816b7ed4-fac6-45ee-8991-1057416f0983", + "metadata": {}, + "source": [ + "Marginal costs in EUR/MWh" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "id": "30a08a5c-a63c-4c0e-9ec0-8fbbd2e26dbb", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Wind 0\n", + "Coal 30\n", + "Gas 60\n", + "Oil 80\n", + "dtype: int64" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "marginal_costs = pd.Series(\n", + " [0, 30, 60, 80], index=[\"Wind\", \"Coal\", \"Gas\", \"Oil\"]\n", + ")\n", + "marginal_costs" + ] + }, + { + "cell_type": "markdown", + "id": "42b0c6e8-6699-4d9c-a958-727283e4ca0a", + "metadata": {}, + "source": [ + "Power plant capacities in MW" + ] + }, + { + "cell_type": "code", + "execution_count": 31, + "id": "5864cf82-59bb-4ec2-9767-4739fd675637", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Wind 3000\n", + "Coal 35000\n", + "Gas 8000\n", + "Oil 2000\n", + "dtype: int64" + ] + }, + "execution_count": 31, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "capacities = pd.Series(\n", + " [3000, 35000, 8000, 2000], index=[\"Wind\", \"Coal\", \"Gas\", \"Oil\"]\n", + ")\n", + "capacities" + ] + }, + { + "cell_type": "markdown", + "id": "fc235ee5-d77f-4c68-bfb8-85048393590e", + "metadata": {}, + "source": [ + "Inelastic demand in MW" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "9b018aa6-715e-4319-ae83-ce7a0fd3d890", + "metadata": {}, + "outputs": [], + "source": [ + "load = 42000" + ] + }, + { + "cell_type": "markdown", + "id": "b7177f81-88e0-42a0-afe1-64360856a2de", + "metadata": {}, + "source": [ + "We now start building the model" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "89d61ab1-078b-406d-b9db-1a6377843c79", + "metadata": {}, + "outputs": [], + "source": [ + "m = linopy.Model()" + ] + }, + { + "cell_type": "markdown", + "id": "aa66cf86-9836-421f-b507-343adc2d5c38", + "metadata": {}, + "source": [ + "Let's define the dispatch variables `g` with the `lower` and `upper` bound:\n", + "$$g_s \\leq G_s $$\n", + "$$g_s \\geq 0 $$" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "cd87f66b-e2f0-4676-966b-51c2f180efb1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable (dim_0: 4)\n", + "-------------------\n", + "[Wind]: g[Wind] ∈ [0, 3000]\n", + "[Coal]: g[Coal] ∈ [0, 3.5e+04]\n", + "[Gas]: g[Gas] ∈ [0, 8000]\n", + "[Oil]: g[Oil] ∈ [0, 2000]" + ] + }, + "execution_count": 34, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g = m.add_variables(\n", + " lower=0, upper=capacities, coords=[capacities.index], name=\"g\"\n", + ")\n", + "g" + ] + }, + { + "cell_type": "markdown", + "id": "7eceb706-c0b4-490e-ad65-e335793d4a68", + "metadata": {}, + "source": [ + "And and the objective to minimize total operational costs:\n", + "$$\\min_{g_s} \\sum_s o_s g_s$$" + ] + }, + { + "cell_type": "code", + "execution_count": 35, + "id": "55cfe8c4-9816-4375-8979-a136219a96fa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Objective:\n", + "----------\n", + "LinearExpression: +0 g[Wind] + 30 g[Coal] + 60 g[Gas] + 80 g[Oil]\n", + "Sense: min\n", + "Value: None" + ] + }, + "execution_count": 35, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_objective(marginal_costs.values * g, sense=\"min\")\n", + "m.objective" + ] + }, + { + "cell_type": "markdown", + "id": "61dfa4ce-d028-463e-a812-66e79ffb798b", + "metadata": {}, + "source": [ + "Which is subject to: \n", + "\n", + "$$\\sum_s g_s = d$$" + ] + }, + { + "cell_type": "code", + "execution_count": 36, + "id": "c37f66f5-fce5-4e54-a122-56b992ed2f95", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `energy_balance`\n", + "---------------------------\n", + "+1 g[Wind] + 1 g[Coal] + 1 g[Gas] + 1 g[Oil] = 42000.0" + ] + }, + "execution_count": 36, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_constraints(g.sum() == load, name=\"energy_balance\")" + ] + }, + { + "cell_type": "markdown", + "id": "36ad37d3-5cbb-4097-9f9f-0440832359fa", + "metadata": {}, + "source": [ + "Then, we can solve the model:" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "1bd321b2-760f-48e4-b734-1c73daaf520d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "1 rows, 3 cols, 3 nonzeros\n", + "0 rows, 0 cols, 0 nonzeros\n", + "Presolve : Reductions: rows 0(-1); columns 0(-4); elements 0(-4) - Reduced to empty\n", + "Solving the original LP from the solution after postsolve\n", + "Model status : Optimal\n", + "Objective value : 1.2900000000e+06\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 37, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "f00e76dd-6831-469e-8a5c-f797a6b7b39c", + "metadata": {}, + "source": [ + "This is the optimimal generator dispatch (MW)" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "876e38fd-008a-4717-ab3b-0aadac8d41bc", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + "
g
dim_0
Wind3000.0
Coal35000.0
Gas4000.0
Oil0.0
\n", + "
" + ], + "text/plain": [ + " g\n", + "dim_0 \n", + "Wind 3000.0\n", + "Coal 35000.0\n", + "Gas 4000.0\n", + "Oil 0.0" + ] + }, + "execution_count": 38, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solution.to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "id": "f77acf05-cf1e-48be-9318-4f4b33efa113", + "metadata": {}, + "source": [ + "And the market clearing price we can read from the shadow price of the energy balance constraint (i.e. the added cost of increasing electricity demand by one unit):" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "a5870398-3420-4443-aed4-ce407ea3ad2d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'energy_balance' ()> Size: 8B\n",
+       "array(60.)
" + ], + "text/plain": [ + " Size: 8B\n", + "array(60.)" + ] + }, + "execution_count": 39, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.dual[\"energy_balance\"]" + ] + }, + { + "cell_type": "markdown", + "id": "5cc105a5-7df0-4b76-90b7-85d8bafa2121", + "metadata": {}, + "source": [ + "### Two bidding zones with transmission" + ] + }, + { + "cell_type": "markdown", + "id": "b1c5e8dc-21fd-4961-a2c5-4b6dae862efe", + "metadata": {}, + "source": [ + "Let's add a spatial dimension, such that the optimisation problem is expanded to\n", + "\\begin{equation}\n", + " \\min_{g_{i,s}, f_\\ell} \\sum_s o_{i,s} g_{i,s}\n", + "\\end{equation}\n", + "such that\n", + "\\begin{align}\n", + " g_{i,s} &\\leq G_{i,s} \\\\\n", + " g_{i,s} &\\geq 0 \\\\\n", + " \\sum_s g_{i,s} - \\sum_\\ell K_{i\\ell} f_\\ell &= d_i & \\text{KCL} \\\\\n", + " |f_\\ell| &\\leq F_\\ell & \\text{line limits} \\\\\n", + " \\sum_\\ell C_{\\ell c} x_\\ell f_\\ell &= 0 & \\text{KVL} \n", + "\\end{align}\n", + "\n", + "In this example, we connect the previous South African electricity system with a hydro generation unit in Mozambique through a single transmission line. Note that because a single transmission line will not result in any cycles, we can neglect KVL in this case.\n", + "\n", + "We are given the following data (all in MW):" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "metadata": {}, + "outputs": [], + "source": [ + "generators = [\"Coal\", \"Wind\", \"Gas\", \"Oil\", \"Hydro\"]\n", + "countries = [\"South_Africa\", \"Mozambique\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "4aaaa1bd-112e-4502-92fe-6291c2998046", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
generatorsCoalWindGasOilHydro
countries
South_Africa350003000800020000
Mozambique00001200
\n", + "
" + ], + "text/plain": [ + "generators Coal Wind Gas Oil Hydro\n", + "countries \n", + "South_Africa 35000 3000 8000 2000 0\n", + "Mozambique 0 0 0 0 1200" + ] + }, + "execution_count": 41, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "capacities = pd.DataFrame(\n", + " {\n", + " \"Coal\": [35000, 0],\n", + " \"Wind\": [3000, 0],\n", + " \"Gas\": [8000, 0],\n", + " \"Oil\": [2000, 0],\n", + " \"Hydro\": [0, 1200],\n", + " },\n", + " index=countries,\n", + ")\n", + "capacities.index.name = \"countries\"\n", + "capacities.columns.name = \"generators\"\n", + "\n", + "capacities" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "id": "8822d16f-1eab-4e7d-af3b-b791113629e0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "generators\n", + "Coal 30\n", + "Wind 0\n", + "Gas 60\n", + "Oil 80\n", + "Hydro 0\n", + "dtype: int64" + ] + }, + "execution_count": 42, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# variable costs in EUR/MWh\n", + "marginal_costs = pd.Series([30, 0, 60, 80, 0], index=generators)\n", + "marginal_costs.index.name = \"generators\"\n", + "marginal_costs" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "01214bbf-786c-4dcf-b3ea-179b5285e51d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "countries\n", + "South_Africa 42000\n", + "Mozambique 650\n", + "dtype: int64" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "load = pd.Series([42000, 650], index=countries)\n", + "load.index.name = \"countries\"\n", + "load" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "20b478e0-26e3-470b-82e6-ed0da0e0217f", + "metadata": {}, + "outputs": [], + "source": [ + "transmission = 500" + ] + }, + { + "cell_type": "markdown", + "id": "e46e0a36-486f-48a9-be71-647208aaa1f5", + "metadata": {}, + "source": [ + "Let's start with a new model instance" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "id": "d2996e54-fd7e-4c99-b52d-02ec6c4ec2f1", + "metadata": {}, + "outputs": [], + "source": [ + "m = linopy.Model()" + ] + }, + { + "cell_type": "markdown", + "id": "436d75a0-8257-4c94-ae18-4462b899fbf5", + "metadata": {}, + "source": [ + "Now we create dispatch variables, as before, with the `upper` and `lower` bound for each countries and generators." + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "id": "83936e48", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
generatorsCoalWindGasOilHydro
countries
South_Africa350003000800020000
Mozambique00001200
\n", + "
" + ], + "text/plain": [ + "generators Coal Wind Gas Oil Hydro\n", + "countries \n", + "South_Africa 35000 3000 8000 2000 0\n", + "Mozambique 0 0 0 0 1200" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "capacities" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "id": "d548c6ae-ae63-4a24-84f4-60cdc0537953", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable (countries: 2, generators: 5)\n", + "--------------------------------------\n", + "[South_Africa, Coal]: g[South_Africa, Coal] ∈ [0, 3.5e+04]\n", + "[South_Africa, Wind]: g[South_Africa, Wind] ∈ [0, 3000]\n", + "[South_Africa, Gas]: g[South_Africa, Gas] ∈ [0, 8000]\n", + "[South_Africa, Oil]: g[South_Africa, Oil] ∈ [0, 2000]\n", + "[South_Africa, Hydro]: g[South_Africa, Hydro] ∈ [0, 0]\n", + "[Mozambique, Coal]: g[Mozambique, Coal] ∈ [0, 0]\n", + "[Mozambique, Wind]: g[Mozambique, Wind] ∈ [0, 0]\n", + "[Mozambique, Gas]: g[Mozambique, Gas] ∈ [0, 0]\n", + "[Mozambique, Oil]: g[Mozambique, Oil] ∈ [0, 0]\n", + "[Mozambique, Hydro]: g[Mozambique, Hydro] ∈ [0, 1200]" + ] + }, + "execution_count": 47, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g = m.add_variables(lower=0, upper=capacities, name=\"g\")\n", + "g" + ] + }, + { + "cell_type": "markdown", + "id": "8c60116f-5d1b-4cf6-b110-211f8d7bd28b", + "metadata": {}, + "source": [ + "We now define the line limit for the transmission line, assuming that power flowing from Mozambique\tto South Africa is positive." + ] + }, + { + "cell_type": "markdown", + "id": "e5d1fc38-c2c9-4fcb-9306-4a077b1cacff", + "metadata": {}, + "source": [ + "The line limit equation can be defined as \n", + "\\begin{align}\n", + "|f_\\ell| &\\leq F_\\ell & \\text{line limits}\n", + "\\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "id": "44d964fb-df92-410d-8465-cb134a0e8f9e", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable\n", + "--------\n", + "flow_MZ_SA ∈ [-500, 500]" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f = m.add_variables(lower=-transmission, upper=transmission, name=\"flow_MZ_SA\")\n", + "f" + ] + }, + { + "cell_type": "markdown", + "id": "3a23ed0a-2d31-470b-a263-e1d4c1bde616", + "metadata": {}, + "source": [ + "The energy balance constraint is replaced by KCL, where we take into account local generation as well as incoming or outgoing flows. The KCL equation can be defined as:\n", + "\\begin{align}\n", + " \\sum_s g_{i,s} - \\sum_\\ell K_{i\\ell} f_\\ell &= d_i & \\text{KCL} \\\\\n", + "\\end{align}\n", + "\n", + "We also need the incidence matrix $K_{i\\ell}$ of this network (here it's very simple!) and assume some direction for the flow variable. Here, we picked the orientation from South Africa to Mozambique. This means that if the values for the flow variable $f_\\ell$ are positive South Africa exports to Mozambique and vice versa if the variable takes negative values." + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "d4d7a22b-0b0d-43cc-bf15-658fa4e945ff", + "metadata": {}, + "outputs": [], + "source": [ + "for country in countries:\n", + " sign = -1 if country == \"Mozambique\" else 1 # minimal incidence matrix\n", + " m.add_constraints(\n", + " g.loc[country].sum() + sign * f == load[country],\n", + " name=f\"{country}_KCL\",\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "b3b94f78", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `Mozambique_KCL`\n", + "---------------------------\n", + "+1 g[Mozambique, Coal] + 1 g[Mozambique, Wind] + 1 g[Mozambique, Gas] + 1 g[Mozambique, Oil] + 1 g[Mozambique, Hydro] - 1 flow_MZ_SA = 650.0" + ] + }, + "execution_count": 50, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.constraints[\"Mozambique_KCL\"]" + ] + }, + { + "cell_type": "markdown", + "id": "f24c3c2a-f6f9-4cfd-9732-ecfab10bcb48", + "metadata": {}, + "source": [ + "The objective can be written as:\n", + "$$\\min_{g_{i,s}, f_\\ell} \\sum_s o_{i,s} g_{i,s}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "id": "75132441-0386-45bf-8f16-e2a6659922aa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "LinearExpression\n", + "----------------\n", + "+30 g[South_Africa, Coal] + 30 g[Mozambique, Coal] + 0 g[South_Africa, Wind] ... +80 g[Mozambique, Oil] + 0 g[South_Africa, Hydro] + 0 g[Mozambique, Hydro]" + ] + }, + "execution_count": 51, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "obj = (g * marginal_costs).sum()\n", + "obj" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "id": "d411da42-7782-40b8-829c-bff883fc2b4b", + "metadata": {}, + "outputs": [], + "source": [ + "m.add_objective(obj, sense=\"min\")" + ] + }, + { + "cell_type": "markdown", + "id": "277b43f9-1dd7-4f3e-a41e-7ad154e943ae", + "metadata": {}, + "source": [ + "We now solve the model." + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "id": "6e41147e-6771-4c04-a1f0-b9ebbe8b0edc", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "1 rows, 4 cols, 4 nonzeros\n", + "0 rows, 0 cols, 0 nonzeros\n", + "Presolve : Reductions: rows 0(-2); columns 0(-11); elements 0(-12) - Reduced to empty\n", + "Solving the original LP from the solution after postsolve\n", + "Model status : Optimal\n", + "Objective value : 1.2600000000e+06\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "7b0462ef-e412-497b-a591-18409deb9c93", + "metadata": {}, + "source": [ + "Now, we print the optimization results" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "id": "a163fa60", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1260000.0" + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.objective.value" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": {}, + "outputs": [ + { + "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", + " \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", + " \n", + "
solution
countriesgenerators
South_AfricaCoal35000.0
Wind3000.0
Gas3500.0
Oil0.0
Hydro0.0
MozambiqueCoal0.0
Wind0.0
Gas0.0
Oil0.0
Hydro1150.0
\n", + "
" + ], + "text/plain": [ + " solution\n", + "countries generators \n", + "South_Africa Coal 35000.0\n", + " Wind 3000.0\n", + " Gas 3500.0\n", + " Oil 0.0\n", + " Hydro 0.0\n", + "Mozambique Coal 0.0\n", + " Wind 0.0\n", + " Gas 0.0\n", + " Oil 0.0\n", + " Hydro 1150.0" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.solution.to_dataframe()" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "78afd5db", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'dual' ()> Size: 8B\n",
+       "array(60.)\n",
+       "Coordinates:\n",
+       "    countries  <U12 48B 'South_Africa'
" + ], + "text/plain": [ + " Size: 8B\n", + "array(60.)\n", + "Coordinates:\n", + " countries \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'dual' ()> Size: 8B\n",
+       "array(-0.)\n",
+       "Coordinates:\n",
+       "    countries  <U10 40B 'Mozambique'
" + ], + "text/plain": [ + " Size: 8B\n", + "array(-0.)\n", + "Coordinates:\n", + " countries \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", + " \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", + "
generatorsCoalWindGasOilHydro
time
010.3111
110.6111
210.4111
310.5111
\n", + "" + ], + "text/plain": [ + "generators Coal Wind Gas Oil Hydro\n", + "time \n", + "0 1 0.3 1 1 1\n", + "1 1 0.6 1 1 1\n", + "2 1 0.4 1 1 1\n", + "3 1 0.5 1 1 1" + ] + }, + "execution_count": 61, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "capacity_factors = pd.DataFrame(\n", + " {\n", + " \"Coal\": 4*[1],\n", + " \"Wind\": [0.3, 0.6, 0.4, 0.5],\n", + " \"Gas\": 4*[1],\n", + " \"Oil\": 4*[1],\n", + " \"Hydro\": 4*[1],\n", + " },\n", + " index=time_index,\n", + " columns=generators,\n", + ")\n", + "capacity_factors.index.name = \"time\"\n", + "capacity_factors.columns.name = \"generators\"\n", + "capacity_factors" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "74658603-cb8f-4cef-9070-ff9b2ed28c78", + "metadata": {}, + "outputs": [], + "source": [ + "load = pd.Series(\n", + " [42000, 43000, 45000, 46000], index=time_index\n", + ")\n", + "load.index.name = \"time\"" + ] + }, + { + "cell_type": "markdown", + "id": "c5e51151-1ff9-4012-898a-5ac93924b396", + "metadata": {}, + "source": [ + "We now start building the model:" + ] + }, + { + "cell_type": "code", + "execution_count": 63, + "id": "46e65d7c-3113-43a3-a091-5f2f8bd3fdfe", + "metadata": {}, + "outputs": [], + "source": [ + "m = linopy.Model()" + ] + }, + { + "cell_type": "markdown", + "id": "39a799b0-81a1-47f3-8dd4-8905cdc1d016", + "metadata": {}, + "source": [ + "Let's define the dispatch variables `g` with the `lower` and `upper` bound:\n", + " \\begin{align}\n", + " g_{s,t} &\\leq \\hat{g}_{s,t} G_{i,s} \\\\\n", + " g_{s,t} &\\geq 0 \\\\\n", + " \\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": 64, + "id": "54c6e92d-df1a-4d81-8d06-1c47ed1ca445", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Variable (time: 4, generators: 5)\n", + "---------------------------------\n", + "[0, Coal]: g[0, Coal] ∈ [0, 3.5e+04]\n", + "[0, Wind]: g[0, Wind] ∈ [0, 900]\n", + "[0, Gas]: g[0, Gas] ∈ [0, 8000]\n", + "[0, Oil]: g[0, Oil] ∈ [0, 2000]\n", + "[0, Hydro]: g[0, Hydro] ∈ [0, 0]\n", + "[1, Coal]: g[1, Coal] ∈ [0, 3.5e+04]\n", + "[1, Wind]: g[1, Wind] ∈ [0, 1800]\n", + "\t\t...\n", + "[2, Oil]: g[2, Oil] ∈ [0, 2000]\n", + "[2, Hydro]: g[2, Hydro] ∈ [0, 0]\n", + "[3, Coal]: g[3, Coal] ∈ [0, 3.5e+04]\n", + "[3, Wind]: g[3, Wind] ∈ [0, 1500]\n", + "[3, Gas]: g[3, Gas] ∈ [0, 8000]\n", + "[3, Oil]: g[3, Oil] ∈ [0, 2000]\n", + "[3, Hydro]: g[3, Hydro] ∈ [0, 0]" + ] + }, + "execution_count": 64, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g = m.add_variables(\n", + " lower=0, upper=capacities * capacity_factors, name=\"g\"\n", + ")\n", + "g" + ] + }, + { + "cell_type": "markdown", + "id": "6c7611d1-e77d-4b08-b8b3-50e3ec4ad91d", + "metadata": {}, + "source": [ + "Then, we add the objective:\n", + "\\begin{equation}\n", + " \\min_{g_{s,t}} \\sum_s o_{s} g_{s,t}\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": 65, + "id": "62ba1212-347d-4b04-901d-f0cae3167e50", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Objective:\n", + "----------\n", + "LinearExpression: +30 g[0, Coal] + 30 g[1, Coal] + 30 g[2, Coal] ... +0 g[1, Hydro] + 0 g[2, Hydro] + 0 g[3, Hydro]\n", + "Sense: min\n", + "Value: None" + ] + }, + "execution_count": 65, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_objective((g * marginal_costs).sum(), sense=\"min\")\n", + "m.objective" + ] + }, + { + "cell_type": "markdown", + "id": "f402a97e-9ab7-401e-9e2c-b569e7da35c5", + "metadata": {}, + "source": [ + "Which is subject to:\n", + "\\begin{align}\n", + " \\sum_s g_{s,t} &= d_t\n", + "\\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "cfc40d59-0780-422f-9c21-f3a59336a230", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `energy_balance` (time: 4):\n", + "--------------------------------------\n", + "[0]: +1 g[0, Coal] + 1 g[0, Wind] + 1 g[0, Gas] + 1 g[0, Oil] + 1 g[0, Hydro] = 42000.0\n", + "[1]: +1 g[1, Coal] + 1 g[1, Wind] + 1 g[1, Gas] + 1 g[1, Oil] + 1 g[1, Hydro] = 43000.0\n", + "[2]: +1 g[2, Coal] + 1 g[2, Wind] + 1 g[2, Gas] + 1 g[2, Oil] + 1 g[2, Hydro] = 45000.0\n", + "[3]: +1 g[3, Coal] + 1 g[3, Wind] + 1 g[3, Gas] + 1 g[3, Oil] + 1 g[3, Hydro] = 46000.0" + ] + }, + "execution_count": 66, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_constraints(\n", + " g.sum(\"generators\") == load,\n", + " name=\"energy_balance\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "f5997e56-a5b6-4ff7-bcf4-a4f4d5bcb10b", + "metadata": {}, + "source": [ + "We now solve the model:" + ] + }, + { + "cell_type": "code", + "execution_count": 67, + "id": "f3b63ce4-d7c6-4b09-8e57-151c0f0040b6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "4 rows, 12 cols, 12 nonzeros\n", + "0 rows, 0 cols, 0 nonzeros\n", + "Presolve : Reductions: rows 0(-4); columns 0(-20); elements 0(-20) - Reduced to empty\n", + "Solving the original LP from the solution after postsolve\n", + "Model status : Optimal\n", + "Objective value : 6.0820000000e+06\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 67, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "16db569f-41c0-4522-9867-9542922116d6", + "metadata": {}, + "source": [ + "We display the results. For ease of reading, we round the results to 2 decimals:" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6082000.0" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.objective.value" + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "0bd273c4", + "metadata": {}, + "outputs": [ + { + "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", + " \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", + "
generatorsCoalWindGasOilHydro
time
035000.0900.06100.00.00.0
135000.01800.06200.00.00.0
235000.01200.08000.0800.00.0
335000.01500.08000.01500.00.0
\n", + "
" + ], + "text/plain": [ + "generators Coal Wind Gas Oil Hydro\n", + "time \n", + "0 35000.0 900.0 6100.0 0.0 0.0\n", + "1 35000.0 1800.0 6200.0 0.0 0.0\n", + "2 35000.0 1200.0 8000.0 800.0 0.0\n", + "3 35000.0 1500.0 8000.0 1500.0 0.0" + ] + }, + "execution_count": 69, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.solution.round(2).to_dataframe().squeeze().unstack()" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "9d78b612", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + "
energy_balance
time
060.0
160.0
280.0
380.0
\n", + "
" + ], + "text/plain": [ + " energy_balance\n", + "time \n", + "0 60.0\n", + "1 60.0\n", + "2 80.0\n", + "3 80.0" + ] + }, + "execution_count": 70, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.dual.to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "id": "b5b9981c-73e0-41a3-a2ba-4d75a6ef2c25", + "metadata": {}, + "source": [ + "### Single bidding zone with several periods and storage" + ] + }, + { + "cell_type": "markdown", + "id": "af49627f-9d2c-4542-bd8d-4a501b4059c5", + "metadata": {}, + "source": [ + "Now, we want to expand the optimisation model with a storage unit to do price arbitrage to reduce oil consumption." + ] + }, + { + "cell_type": "markdown", + "id": "7882d6ab-3108-49c5-b904-6d55602c853d", + "metadata": {}, + "source": [ + "We have been given the following characteristics of the storage:" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "30b6d318-91e9-4074-b749-b2d1393e503e", + "metadata": {}, + "outputs": [], + "source": [ + "storage_energy = 6000 # MWh\n", + "storage_power = 1000 # MW\n", + "efficiency = 0.9 # discharge = charge\n", + "standing_loss = 0.00001 # per hour" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "596f621c-70e1-4bb4-8701-e706bbd756e3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Linopy LP model\n", + "===============\n", + "\n", + "Variables:\n", + "----------\n", + " * g (time, generators)\n", + "\n", + "Constraints:\n", + "------------\n", + " * energy_balance (time)\n", + "\n", + "Status:\n", + "-------\n", + "ok" + ] + }, + "execution_count": 72, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m" + ] + }, + { + "cell_type": "markdown", + "id": "087352e8-0e95-4120-a7b0-508ee148e034", + "metadata": {}, + "source": [ + "To model a storage unit, we need three additional variables for the discharging and charging of the storage unit and for its state of charge (energy filling level). We can directly define the bounds of these variables in the variable definition:" + ] + }, + { + "cell_type": "code", + "execution_count": 73, + "id": "2d32118c-eeda-4ab4-8400-292aca5fa0ef", + "metadata": {}, + "outputs": [], + "source": [ + "battery_discharge = m.add_variables(\n", + " lower=0, upper=storage_power, coords=[time_index], name=\"battery_discharge\"\n", + ")\n", + "battery_charge = m.add_variables(\n", + " lower=0, upper=storage_power, coords=[time_index], name=\"battery_charge\"\n", + ")\n", + "battery_soc = m.add_variables(\n", + " lower=0, upper=storage_energy, coords=[time_index], name=\"battery_soc\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "e874f80b-7482-4313-8418-87b213fa0441", + "metadata": {}, + "source": [ + "Then, we implement the storage consistency equations,\n", + "\n", + "$$e_{t} = (1-\\text{standing loss}) \\cdot e_{t-1} + \\eta \\cdot g_{charge, t} - \\frac{1}{\\eta} \\cdot g_{discharge, t}$$\n", + "\n", + "For the initial period, we set the state of charge to zero.\n", + "\n", + "$$e_{0} = 0$$" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "id": "37f57529-a84c-4553-86ee-180211a53af8", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `soc_initial`\n", + "------------------------\n", + "+1 battery_soc[0] = -0.0" + ] + }, + "execution_count": 74, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_constraints(battery_soc.loc[0] == 0, name=\"soc_initial\")" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "id": "2bad2277", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `soc_consistency` (time: 3):\n", + "---------------------------------------\n", + "[1]: +1 battery_soc[1] - 1 battery_soc[0] - 0.9 battery_charge[1] + 1.111 battery_discharge[1] = -0.0\n", + "[2]: +1 battery_soc[2] - 1 battery_soc[1] - 0.9 battery_charge[2] + 1.111 battery_discharge[2] = -0.0\n", + "[3]: +1 battery_soc[3] - 1 battery_soc[2] - 0.9 battery_charge[3] + 1.111 battery_discharge[3] = -0.0" + ] + }, + "execution_count": 75, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_constraints(\n", + " battery_soc.loc[1:] == (1 - standing_loss) * battery_soc.shift(time=1).loc[1:] + efficiency * battery_charge.loc[1:] - 1 / efficiency * battery_discharge.loc[1:],\n", + " name=\"soc_consistency\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "726ff9ca-900f-47e9-9216-b79ad3202a99", + "metadata": {}, + "source": [ + "And we also need to modify the energy balance to include the contributions of storage discharging and charging.\n", + "\n", + "For that, we should first remove the existing energy balance constraint, which we seek to overwrite." + ] + }, + { + "cell_type": "code", + "execution_count": 76, + "id": "e1f690ff-55d0-45ff-901e-58c8adb43c18", + "metadata": {}, + "outputs": [], + "source": [ + "m.remove_constraints(\"energy_balance\")" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "id": "7873dec7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `energy_balance` (time: 4):\n", + "--------------------------------------\n", + "[0]: +1 g[0, Coal] + 1 g[0, Wind] + 1 g[0, Gas] ... +1 g[0, Hydro] + 1 battery_discharge[0] - 1 battery_charge[0] = 42000.0\n", + "[1]: +1 g[1, Coal] + 1 g[1, Wind] + 1 g[1, Gas] ... +1 g[1, Hydro] + 1 battery_discharge[1] - 1 battery_charge[1] = 43000.0\n", + "[2]: +1 g[2, Coal] + 1 g[2, Wind] + 1 g[2, Gas] ... +1 g[2, Hydro] + 1 battery_discharge[2] - 1 battery_charge[2] = 45000.0\n", + "[3]: +1 g[3, Coal] + 1 g[3, Wind] + 1 g[3, Gas] ... +1 g[3, Hydro] + 1 battery_discharge[3] - 1 battery_charge[3] = 46000.0" + ] + }, + "execution_count": 77, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.add_constraints(\n", + " g.sum(\"generators\") + battery_discharge - battery_charge == load,\n", + " name=\"energy_balance\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "e9769874-f052-4b88-88a7-5e3d4334e4b5", + "metadata": {}, + "source": [ + "We now solve the model:" + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "id": "5700dbf1-8ced-48f4-94b0-d991f5d1b4b0", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d-dirty]\n", + "Copyright (c) 2023 HiGHS under MIT licence terms\n", + "Presolving model\n", + "7 rows, 21 cols, 29 nonzeros\n", + "4 rows, 10 cols, 15 nonzeros\n", + "Presolve : Reductions: rows 4(-4); columns 10(-22); elements 15(-26)\n", + "Solving the presolved LP\n", + "Using EKK dual simplex solver - serial\n", + " Iteration Objective Infeasibilities num(sum)\n", + " 0 5.3580000000e+06 Pr: 2(10300) 0s\n", + " 7 6.0172006560e+06 Pr: 0(0) 0s\n", + " 7 6.0172006560e+06 Pr: 0(0) 0s\n", + "Solving the original LP from the solution after postsolve\n", + "Model status : Optimal\n", + "Simplex iterations: 7\n", + "Objective value : 6.0172006560e+06\n", + "HiGHS run time : 0.00\n" + ] + }, + { + "data": { + "text/plain": [ + "('ok', 'optimal')" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.solve()" + ] + }, + { + "cell_type": "markdown", + "id": "80f712d8-b070-4eb3-9d3f-2e9ec560f9df", + "metadata": {}, + "source": [ + "We display the results:" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "b92c2f1e-113e-4dfc-9ec0-30c3ae680123", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "6017200.655993455" + ] + }, + "execution_count": 79, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "m.objective.value" + ] + }, + { + "cell_type": "code", + "execution_count": 80, + "id": "413c8dab", + "metadata": {}, + "outputs": [ + { + "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", + " \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", + "
generatorsCoalWindGasOilHydro
time
035000.0900.05100.00.00000.0
135000.01800.07200.00.00000.0
235000.01200.08000.00.00000.0
335000.01500.08000.01490.00820.0
\n", + "
" + ], + "text/plain": [ + "generators Coal Wind Gas Oil Hydro\n", + "time \n", + "0 35000.0 900.0 5100.0 0.0000 0.0\n", + "1 35000.0 1800.0 7200.0 0.0000 0.0\n", + "2 35000.0 1200.0 8000.0 0.0000 0.0\n", + "3 35000.0 1500.0 8000.0 1490.0082 0.0" + ] + }, + "execution_count": 80, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.solution.to_dataframe().squeeze().unstack()" + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "id": "ed4ffe65-82ac-4565-a11a-fa5d99f977f5", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + "
solution
time
01000.0000
10.0000
2800.0000
39.9918
\n", + "
" + ], + "text/plain": [ + " solution\n", + "time \n", + "0 1000.0000\n", + "1 0.0000\n", + "2 800.0000\n", + "3 9.9918" + ] + }, + "execution_count": 81, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "battery_discharge.solution.to_dataframe()" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "id": "2ff75605-f715-4f1c-9b39-0e2b46d0072b", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + "
solution
time
00.0
11000.0
20.0
30.0
\n", + "
" + ], + "text/plain": [ + " solution\n", + "time \n", + "0 0.0\n", + "1 1000.0\n", + "2 0.0\n", + "3 0.0" + ] + }, + "execution_count": 82, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "battery_charge.solution.to_dataframe()" + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "id": "41b539b7-1011-4fd9-b30b-b474c4e365cb", + "metadata": {}, + "outputs": [ + { + "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", + " \n", + " \n", + "
solution
time
0-0.000000
1900.000000
211.102111
30.000000
\n", + "
" + ], + "text/plain": [ + " solution\n", + "time \n", + "0 -0.000000\n", + "1 900.000000\n", + "2 11.102111\n", + "3 0.000000" + ] + }, + "execution_count": 83, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "battery_soc.solution.to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "id": "3c3386c4-946f-4f8e-8582-18730f60be83", + "metadata": {}, + "source": [ + "### Exercise" + ] + }, + { + "cell_type": "markdown", + "id": "c4d9eb93-073e-4725-957d-ebc2aa603b0c", + "metadata": {}, + "source": [ + "- Using the conversion efficiencies and specific emissions from the lecture slides, add a constraint that limits the total emissions in the four periods to 50% of the unconstrained optimal solution. How does the optimal objective value and the generator dispatch change?\n", + "\n", + "- Reimplement the storage consistency constraint such that the initial state of charge is not zero but corresponds to the state of charge in the final period of the optimisation horizon.\n", + "\n", + "- What parameters of the storage unit would have to be changed to reduce the objective? What's the sensitivity?" + ] + }, + { + "cell_type": "markdown", + "id": "231d90d1", + "metadata": {}, + "source": [] + }, + { + "cell_type": "markdown", + "id": "d97637fb", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "esm-2024", + "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.11.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/data-science-for-esm/_toc.yml b/data-science-for-esm/_toc.yml index 6e458888..07960525 100644 --- a/data-science-for-esm/_toc.yml +++ b/data-science-for-esm/_toc.yml @@ -11,9 +11,10 @@ chapters: - file: 06-workshop-atlite - file: 05-workshop-pysheds - file: 07-workshop-networkx -- file: 08-workshop-pyomo +- file: 14-workshop-linopy - file: 09-workshop-pypsa - file: 10-workshop-pypsa-cem - file: 11-workshop-groupwork - file: 12-workshop-pypsa-sector-coupling - file: 13-workshop-interactive-visualisation +- file: 08-workshop-pyomo