From f546f2f2b6a2a3eba859111dbd444f8ab6761854 Mon Sep 17 00:00:00 2001 From: Fabian Neumann Date: Mon, 27 May 2024 14:38:54 +0200 Subject: [PATCH] simplify linopy tutorial --- data-science-for-esm/14-workshop-linopy.ipynb | 4696 ++++++++++++++++- data-science-for-esm/_toc.yml | 4 +- 2 files changed, 4427 insertions(+), 273 deletions(-) diff --git a/data-science-for-esm/14-workshop-linopy.ipynb b/data-science-for-esm/14-workshop-linopy.ipynb index 47dcf300..5a675eb9 100644 --- a/data-science-for-esm/14-workshop-linopy.ipynb +++ b/data-science-for-esm/14-workshop-linopy.ipynb @@ -61,7 +61,7 @@ "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", + "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", @@ -84,13 +84,12 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "d94292b1-8073-49a8-a2e9-a2c769145b2d", "metadata": {}, "outputs": [], "source": [ "import linopy\n", - "import numpy as np\n", "\n", "m = linopy.Model()" ] @@ -103,19 +102,29 @@ "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 variables\n", + "### Adding decision variables\n", "\n", - "Variables in a linear optimization problem represent the decision variables. 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:" + "**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": null, + "execution_count": 2, "id": "64b35f1c-5398-4eb3-9135-aafce270b79b", "metadata": {}, "outputs": [], @@ -134,14 +143,74 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "19d4b5f5-b9ae-48f6-8ce6-cd4a418dba6b", "metadata": {}, - "outputs": [], + "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", @@ -149,15 +218,28 @@ "source": [ "### Adding Constraints\n", "\n", - "Constraints define the feasible region of the optimization problem. 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:" + "**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": null, + "execution_count": 6, "id": "1df4589f-4c99-46b5-a530-63f99b4515cc", "metadata": {}, - "outputs": [], + "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" ] @@ -172,10 +254,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "id": "688e0f65-9198-49b2-915a-2131084224c1", "metadata": {}, - "outputs": [], + "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" ] @@ -185,7 +280,7 @@ "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." + "… and linopy will automatically take over the separation of variables expression on the LHS, and constant values on the RHS." ] }, { @@ -193,12 +288,12 @@ "id": "f18d18a3-671f-4190-9257-4541e6c88c61", "metadata": {}, "source": [ - "The constraint is currently not assigned to the model. We assign it by calling the function `m.add_constraints`" + "The constraint is currently not assigned to the model. We assign it by calling the `add_constraints()` function:" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "id": "fa233f5c-38ec-4a5a-9008-502932df968d", "metadata": {}, "outputs": [], @@ -207,6 +302,53 @@ "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", @@ -214,12 +356,12 @@ "source": [ "### Adding the Objective \n", "\n", - "The objective function defines what you want to optimize. You can set the objective function of a Model in Linopy using the add_objective method. For our example that would be" + "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": null, + "execution_count": 11, "id": "18c6e6ce-37f9-4527-83a4-e5764b67b34c", "metadata": {}, "outputs": [], @@ -227,12 +369,73 @@ "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." + "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" ] }, { @@ -242,47 +445,1280 @@ "source": [ "### Solving the Model\n", "\n", - "Once you've defined your Model with variables, constraints, and an objective function, you can solve it using the `solve` method:" + "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": null, + "execution_count": 14, "id": "712a0033-d027-478a-b8f2-201ec4ed1cc1", "metadata": {}, - "outputs": [], + "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": [ - "The solution of the linear problem assigned to the variables under `solution` in form of a `xarray.Dataset`." + "### 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": null, + "execution_count": 15, "id": "1f28aa81-d16f-4b7e-9d8e-f661d4b8bcd8", "metadata": {}, - "outputs": [], + "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": null, + "execution_count": 16, "id": "470e2d25-9386-438d-932e-113464907726", "metadata": {}, - "outputs": [], + "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", @@ -333,7 +1769,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "id": "788956e7-1bc9-4298-85a1-9235486b8cad", "metadata": {}, "outputs": [], @@ -346,12 +1782,12 @@ "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." + "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": null, + "execution_count": 20, "id": "d477e41c-a89f-4d3b-a1af-385820638a75", "metadata": {}, "outputs": [], @@ -368,23 +1804,78 @@ "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 rhs have to match. \n", + "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. Like that, things remain clear and no unexpected broadcasting (which we show later) will happen. " + ":::{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": null, + "execution_count": 22, "id": "00a64dec-26a5-4f80-97c0-de0fe00188a2", "metadata": {}, - "outputs": [], + "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", @@ -401,10 +1892,36 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "id": "8d3a89e8-0b0e-480d-9cb8-f29931fb3559", "metadata": {}, - "outputs": [], + "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", @@ -421,21 +1938,94 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "id": "9cbf16b4-99ed-4e33-9ea9-e8eff3503b05", "metadata": {}, "outputs": [], "source": [ - "obj = (x + 2 * y).sum()\n", - "m.add_objective(obj)" + "obj = (x + 2 * y).sum()\n" ] }, { "cell_type": "code", - "execution_count": null, - "id": "37da40ed-71f5-4c3c-825b-26e0e2da46e2", + "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()" ] @@ -445,15 +2035,149 @@ "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. This can sometimes be helpful." + "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": null, + "execution_count": 28, "id": "33cdfad9-0ff3-4211-afaf-b0e27fa33d5a", "metadata": {}, - "outputs": [], + "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\");" ] @@ -518,14 +2242,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 30, "id": "30a08a5c-a63c-4c0e-9ec0-8fbbd2e26dbb", "metadata": {}, - "outputs": [], - "source": [ - "marginals_df = pd.DataFrame(\n", - " {\"Generator\": [\"Wind\", \"Coal\", \"Gas\", \"Oil\"], \"MarginalCost\": [0, 30, 60, 80]}\n", - ")" + "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" ] }, { @@ -538,14 +2278,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 31, "id": "5864cf82-59bb-4ec2-9767-4739fd675637", "metadata": {}, - "outputs": [], - "source": [ - "capacities_df = pd.DataFrame(\n", - " {\"Generator\": [\"Wind\", \"Coal\", \"Gas\", \"Oil\"], \"Capacity\": [3000, 35000, 8000, 2000]}\n", - ")" + "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" ] }, { @@ -558,7 +2314,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 32, "id": "9b018aa6-715e-4319-ae83-ce7a0fd3d890", "metadata": {}, "outputs": [], @@ -576,7 +2332,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 33, "id": "89d61ab1-078b-406d-b9db-1a6377843c79", "metadata": {}, "outputs": [], @@ -596,13 +2352,29 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 34, "id": "cd87f66b-e2f0-4676-966b-51c2f180efb1", "metadata": {}, - "outputs": [], + "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_df.Capacity, coords=[capacities_df.Generator], name=\"g\"\n", + " lower=0, upper=capacities, coords=[capacities.index], name=\"g\"\n", ")\n", "g" ] @@ -612,18 +2384,33 @@ "id": "7eceb706-c0b4-490e-ad65-e335793d4a68", "metadata": {}, "source": [ - "And and the objective:\n", + "And and the objective to minimize total operational costs:\n", "$$\\min_{g_s} \\sum_s o_s g_s$$" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 35, "id": "55cfe8c4-9816-4375-8979-a136219a96fa", "metadata": {}, - "outputs": [], - "source": [ - "m.add_objective(marginals_df.MarginalCost.values * g, sense=\"min\")\n", + "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" ] }, @@ -639,10 +2426,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 36, "id": "c37f66f5-fce5-4e54-a122-56b992ed2f95", "metadata": {}, - "outputs": [], + "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\")" ] @@ -652,15 +2452,42 @@ "id": "36ad37d3-5cbb-4097-9f9f-0440832359fa", "metadata": {}, "source": [ - "It always helps to write out the constraints before adding them to the model. Since they look good, let’s assign them." + "Then, we can solve the model:" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 37, "id": "1bd321b2-760f-48e4-b734-1c73daaf520d", "metadata": {}, - "outputs": [], + "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()" ] @@ -675,10 +2502,73 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 38, "id": "876e38fd-008a-4717-ab3b-0aadac8d41bc", "metadata": {}, - "outputs": [], + "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()" ] @@ -693,12 +2583,391 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 39, "id": "a5870398-3420-4443-aed4-ce407ea3ad2d", "metadata": {}, - "outputs": [], - "source": [ - "m.dual[\"energy_balance\"].item()" + "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\"]" ] }, { @@ -734,7 +3003,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 40, "metadata": {}, "outputs": [], "source": [ @@ -744,12 +3013,81 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 41, "id": "4aaaa1bd-112e-4502-92fe-6291c2998046", "metadata": {}, - "outputs": [], - "source": [ - "capacities_df = pd.DataFrame(\n", + "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", @@ -759,36 +3097,71 @@ " },\n", " index=countries,\n", ")\n", + "capacities.index.name = \"countries\"\n", + "capacities.columns.name = \"generators\"\n", "\n", - "capacities_df" + "capacities" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 42, "id": "8822d16f-1eab-4e7d-af3b-b791113629e0", "metadata": {}, - "outputs": [], + "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_series = pd.Series([30, 0, 60, 80, 0], index=generators)\n", - "marginal_series" + "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": null, + "execution_count": 43, "id": "01214bbf-786c-4dcf-b3ea-179b5285e51d", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "countries\n", + "South_Africa 42000\n", + "Mozambique 650\n", + "dtype: int64" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "loads_series = pd.Series([42000, 650], index=countries)\n", - "loads_series" + "load = pd.Series([42000, 650], index=countries)\n", + "load.index.name = \"countries\"\n", + "load" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 44, "id": "20b478e0-26e3-470b-82e6-ed0da0e0217f", "metadata": {}, "outputs": [], @@ -806,7 +3179,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 45, "id": "d2996e54-fd7e-4c99-b52d-02ec6c4ec2f1", "metadata": {}, "outputs": [], @@ -824,25 +3197,113 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "4c4bb906-4fc6-4264-8cac-47ef3fec6cd3", - "metadata": {}, - "outputs": [], - "source": [ - "countries_index = pd.Index(countries, name=\"countries\")\n", - "generators_index = pd.Index(generators, name=\"generators\")" + "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": null, + "execution_count": 47, "id": "d548c6ae-ae63-4a24-84f4-60cdc0537953", "metadata": {}, - "outputs": [], - "source": [ - "g = m.add_variables(\n", - " lower=0, upper=capacities_df, coords=[countries_index, generators_index], name=\"g\"\n", - ")\n", + "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" ] }, @@ -867,13 +3328,26 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 48, "id": "44d964fb-df92-410d-8465-cb134a0e8f9e", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "Variable\n", + "--------\n", + "flow_MZ_SA ∈ [-500, 500]" + ] + }, + "execution_count": 48, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "f_line = m.add_variables(lower=-transmission, upper=transmission, name=\"line_limit\")\n", - "f_line" + "f = m.add_variables(lower=-transmission, upper=transmission, name=\"flow_MZ_SA\")\n", + "f" ] }, { @@ -891,7 +3365,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 49, "id": "d4d7a22b-0b0d-43cc-bf15-658fa4e945ff", "metadata": {}, "outputs": [], @@ -899,43 +3373,70 @@ "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_line == loads_series.loc[country],\n", + " g.loc[country].sum() + sign * f == load[country],\n", " name=f\"{country}_KCL\",\n", " )" ] }, { - "cell_type": "markdown", - "id": "f24c3c2a-f6f9-4cfd-9732-ecfab10bcb48", + "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": [ - "The objective can be written as:\n", - "$$\\min_{g_{i,s}, f_\\ell} \\sum_s o_{i,s} g_{i,s}$$" + "m.constraints[\"Mozambique_KCL\"]" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", + "id": "f24c3c2a-f6f9-4cfd-9732-ecfab10bcb48", "metadata": {}, - "outputs": [], "source": [ - "o_is = [marginal_series] * len(countries)" + "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": null, + "execution_count": 51, "id": "75132441-0386-45bf-8f16-e2a6659922aa", "metadata": {}, - "outputs": [], - "source": [ - "obj = (g * o_is).sum()\n", + "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": null, + "execution_count": 52, "id": "d411da42-7782-40b8-829c-bff883fc2b4b", "metadata": {}, "outputs": [], @@ -953,10 +3454,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 53, "id": "6e41147e-6771-4c04-a1f0-b9ebbe8b0edc", "metadata": {}, - "outputs": [], + "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()" ] @@ -971,13 +3499,917 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 54, + "id": "a163fa60", "metadata": {}, - "outputs": [], + "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\": [35000] * 4,\n", - " \"Wind\": wind_ts,\n", - " \"Gas\": [8000] * 4,\n", - " \"Oil\": [2000] * 4,\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", - "capacities_ts" + "capacity_factors.index.name = \"time\"\n", + "capacity_factors.columns.name = \"generators\"\n", + "capacity_factors" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 62, "id": "74658603-cb8f-4cef-9070-ff9b2ed28c78", "metadata": {}, "outputs": [], "source": [ - "load_ts = [42000, 43000, 45000, 46000]" + "load = pd.Series(\n", + " [42000, 43000, 45000, 46000], index=time_index\n", + ")\n", + "load.index.name = \"time\"" ] }, { @@ -1104,7 +4615,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 63, "id": "46e65d7c-3113-43a3-a091-5f2f8bd3fdfe", "metadata": {}, "outputs": [], @@ -1126,13 +4637,40 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 64, "id": "54c6e92d-df1a-4d81-8d06-1c47ed1ca445", "metadata": {}, - "outputs": [], + "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_ts, coords=[time_index, generators_index], name=\"g\"\n", + " lower=0, upper=capacities * capacity_factors, name=\"g\"\n", ")\n", "g" ] @@ -1150,28 +4688,27 @@ }, { "cell_type": "code", - "execution_count": null, - "id": "2fe3cf97-7da2-4f18-96d4-b09e314413b6", - "metadata": {}, - "outputs": [], - "source": [ - "o_s = m.add_variables(\n", - " lower=marginal_series,\n", - " upper=marginal_series,\n", - " coords=[generators_index],\n", - " name=\"o_s\",\n", - ")\n", - "o_s" - ] - }, - { - "cell_type": "code", - "execution_count": null, + "execution_count": 65, "id": "62ba1212-347d-4b04-901d-f0cae3167e50", "metadata": {}, - "outputs": [], - "source": [ - "m.add_objective((o_s * g).sum(), sense=\"min\")\n", + "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" ] }, @@ -1188,14 +4725,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 66, "id": "cfc40d59-0780-422f-9c21-f3a59336a230", "metadata": {}, - "outputs": [], + "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_ts,\n", - " name=f\"energy_balance\",\n", + " g.sum(\"generators\") == load,\n", + " name=\"energy_balance\",\n", ")" ] }, @@ -1209,10 +4762,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 67, "id": "f3b63ce4-d7c6-4b09-8e57-151c0f0040b6", "metadata": {}, - "outputs": [], + "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()" ] @@ -1227,11 +4807,192 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 68, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "6082000.0" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "g.solution.round(2).to_dataframe()" + "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()" ] }, { @@ -1260,7 +5021,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 71, "id": "30b6d318-91e9-4074-b749-b2d1393e503e", "metadata": {}, "outputs": [], @@ -1273,10 +5034,34 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 72, "id": "596f621c-70e1-4bb4-8701-e706bbd756e3", "metadata": {}, - "outputs": [], + "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" ] @@ -1291,7 +5076,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 73, "id": "2d32118c-eeda-4ab4-8400-292aca5fa0ef", "metadata": {}, "outputs": [], @@ -1318,33 +5103,58 @@ "\n", "For the initial period, we set the state of charge to zero.\n", "\n", - "$$e_{0} = 0$$\n", - "\n", - "We also set the charging power and discharging power to zero.\n", - "\n", - "$$g_{discharge, 0} = 0$$\n", - "$$g_{charge, 0} = 0$$" + "$$e_{0} = 0$$" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 74, "id": "37f57529-a84c-4553-86ee-180211a53af8", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "Constraint `soc_initial`\n", + "------------------------\n", + "+1 battery_soc[0] = -0.0" + ] + }, + "execution_count": 74, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "for time_value in time_index:\n", - " if time_value == 0:\n", - " m.add_constraints(battery_soc.loc[0] == 0)\n", - " m.add_constraints(battery_charge.loc[0] == 0)\n", - " m.add_constraints(battery_discharge.loc[0] == 0)\n", - " else:\n", - " m.add_constraints(\n", - " battery_soc.loc[time_value]\n", - " == (1 - standing_loss) * battery_soc.loc[time_value - 1]\n", - " + efficiency * battery_charge.loc[time_value]\n", - " - 1 / efficiency * battery_discharge.loc[time_value]\n", - " )" + "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", + ")" ] }, { @@ -1359,40 +5169,41 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 76, "id": "e1f690ff-55d0-45ff-901e-58c8adb43c18", "metadata": {}, "outputs": [], "source": [ - "m.remove_constraints(f\"energy_balance\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "fcad7972-b06f-42d2-adab-80dd248bbd17", - "metadata": {}, - "outputs": [], - "source": [ - "g.loc[3, :].sum() + battery_discharge.loc[3] - battery_charge.loc[3] == load_ts[3]\n", - "# constraint at time \"3\"" + "m.remove_constraints(\"energy_balance\")" ] }, { "cell_type": "code", - "execution_count": null, - "id": "ffb5ad87-2484-4bd9-be58-5c544325bb3d", - "metadata": {}, - "outputs": [], + "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": [ - "for time_value in time_index:\n", - " m.add_constraints(\n", - " g.loc[time_value, :].sum()\n", - " + battery_discharge.loc[time_value]\n", - " - battery_charge.loc[time_value]\n", - " == load_ts[time_value],\n", - " name=f\"energy_balance_{time_value}\",\n", - " )" + "m.add_constraints(\n", + " g.sum(\"generators\") + battery_discharge - battery_charge == load,\n", + " name=\"energy_balance\",\n", + ")" ] }, { @@ -1405,10 +5216,44 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 78, "id": "5700dbf1-8ced-48f4-94b0-d991f5d1b4b0", "metadata": {}, - "outputs": [], + "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()" ] @@ -1423,40 +5268,337 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 79, "id": "b92c2f1e-113e-4dfc-9ec0-30c3ae680123", "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "6017200.655993455" + ] + }, + "execution_count": 79, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "g.solution.to_dataframe()" + "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": null, + "execution_count": 81, "id": "ed4ffe65-82ac-4565-a11a-fa5d99f977f5", "metadata": {}, - "outputs": [], + "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": null, + "execution_count": 82, "id": "2ff75605-f715-4f1c-9b39-0e2b46d0072b", "metadata": {}, - "outputs": [], + "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": null, + "execution_count": 83, "id": "41b539b7-1011-4fd9-b30b-b474c4e365cb", "metadata": {}, - "outputs": [], + "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()" ] @@ -1480,13 +5622,25 @@ "\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": "", + "display_name": "esm-2024", "language": "python", - "name": "" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1498,7 +5652,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.12" + "version": "3.11.8" } }, "nbformat": 4, diff --git a/data-science-for-esm/_toc.yml b/data-science-for-esm/_toc.yml index 749a45f8..07960525 100644 --- a/data-science-for-esm/_toc.yml +++ b/data-science-for-esm/_toc.yml @@ -11,10 +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: 14-workshop-linopy +- file: 08-workshop-pyomo