examples/notebooks/optimization-with-linopy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "raw_mimetype": "text/restructuredtext",
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# Optimization with Linopy\n",
    "\n",
    "In PyPSA `v0.22`, an additional optimization module was introduced to the package. It is built on [Linopy](https://github.com/PyPSA/linopy) and aims at \n",
    "\n",
    "* **performance** as we know it from the native PyPSA optimization (`lopf` with `pyomo=False`, now deprecated)\n",
    "* **flexibility** as we know from the Pyomo implementation\n",
    "* **usability** as we know from pandas/xarray\n",
    "\n",
    "\n",
    "Linopy is a stand-alone package and works similar to Pyomo, but without the memory overhead and much faster. If you come from an older PyPSA version, you might be familiar with the pyomo-based ``n.lopf`` function. The Linopy-based optimization is the successor, is now the default optimization method in PyPSA, and is accessed via the [pypsa.Network.optimize](https://pypsa.readthedocs.io/en/latest/api/optimization.html) accessor. See the [migration guide](https://pypsa.readthedocs.io/en/stable/examples/optimization-with-linopy-migrate-extra-functionalities.html) for more information.\n",
    "\n",
    "For additional information on the Linopy package, have a look at the [documentation](https://linopy.readthedocs.io/en/latest/).\n",
    "\n",
    "## Let's get started\n",
    "\n",
    "Now, we demonstrate the behaviour of the optimization with linopy. The core functions for the optimization can be called via the [pypsa.Network.optimize](https://pypsa.readthedocs.io/en/latest/api/optimization.html) accessor. The accessor is used for creating, solving, modifying the optimization problem. Further, it supports to run different optimization formulations and provides helper functions. \n",
    "\n",
    "At first, we run the ordinary linearized optimal power flow (LOPF). We then extend the formulation by some additional constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pypsa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = pypsa.examples.ac_dc_meshed(from_master=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to make the network a bit more interesting, we modify its data: We set gas generators to non-extendable,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators.loc[n.generators.carrier == \"gas\", \"p_nom_extendable\"] = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... add ramp limits,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators.loc[n.generators.carrier == \"gas\", \"ramp_limit_down\"] = 0.2\n",
    "n.generators.loc[n.generators.carrier == \"gas\", \"ramp_limit_up\"] = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... add additional storage units (cyclic and non-cyclic) and fix one state_of_charge,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.add(\n",
    "    \"StorageUnit\",\n",
    "    \"su\",\n",
    "    bus=\"Manchester\",\n",
    "    marginal_cost=10,\n",
    "    inflow=50,\n",
    "    p_nom_extendable=True,\n",
    "    capital_cost=10,\n",
    "    p_nom=2000,\n",
    "    efficiency_dispatch=0.5,\n",
    "    cyclic_state_of_charge=True,\n",
    "    state_of_charge_initial=1000,\n",
    ")\n",
    "\n",
    "n.add(\n",
    "    \"StorageUnit\",\n",
    "    \"su2\",\n",
    "    bus=\"Manchester\",\n",
    "    marginal_cost=10,\n",
    "    p_nom_extendable=True,\n",
    "    capital_cost=50,\n",
    "    p_nom=2000,\n",
    "    efficiency_dispatch=0.5,\n",
    "    carrier=\"gas\",\n",
    "    cyclic_state_of_charge=False,\n",
    "    state_of_charge_initial=1000,\n",
    ")\n",
    "\n",
    "n.storage_units_t.state_of_charge_set.loc[n.snapshots[7], \"su\"] = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and add an additional store."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.add(\"Bus\", \"storebus\", carrier=\"hydro\", x=-5, y=55)\n",
    "n.add(\n",
    "    \"Link\",\n",
    "    [\"battery_power\", \"battery_discharge\"],\n",
    "    \"\",\n",
    "    bus0=[\"Manchester\", \"storebus\"],\n",
    "    bus1=[\"storebus\", \"Manchester\"],\n",
    "    p_nom=100,\n",
    "    efficiency=0.9,\n",
    "    p_nom_extendable=True,\n",
    "    p_nom_max=1000,\n",
    ")\n",
    "n.add(\n",
    "    \"Store\",\n",
    "    [\"store\"],\n",
    "    bus=\"storebus\",\n",
    "    e_nom=2000,\n",
    "    e_nom_extendable=True,\n",
    "    marginal_cost=10,\n",
    "    capital_cost=10,\n",
    "    e_nom_max=5000,\n",
    "    e_initial=100,\n",
    "    e_cyclic=True,\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Optimization\n",
    "\n",
    "The optimization based on linopy mimics the well-known [`n.lopf`](https://pypsa.readthedocs.io/en/v0.28.0/api_reference.html#pypsa.Network.lopf) optimization. We run it by calling the `optimize` accessor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a model instance attached to our network. It is a container of all variables, constraints and the objective function. You can modify this as much as you please, by directly adding or deleting variables or constraints etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modify model, optimize and feed back to network "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When you have a fresh network and you just want to create the model instance, run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.optimize.create_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Through the model instance we gain a lot of flexibility. Let's say for example we want to remove the Kirchhoff Voltage Law constraint, thus convert the model to a transport model. This can be done via"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model.constraints.remove(\"Kirchhoff-Voltage-Law\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we want to optimize the altered model and feed to solution back to the network. Here again, we use the `optimize` accessor."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.optimize.solve_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we followed the recommended way to run altered models:\n",
    "\n",
    "1. **Create the model instance** - `n.optimize.create_model()`\n",
    "2. **Modify the model to your needs**\n",
    "3. **Solve and feed back** - `n.optimize.solve_model()`\n",
    "\n",
    "\n",
    "For compatibility reasons the `optimize` function, also allows passing a `extra_funcionality` argument, as we know it from the `lopf` function. The above behaviour with use of the extra functionality is obtained through"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def remove_kvl(n, sns):\n",
    "    print(\"KVL removed!\")\n",
    "    n.model.constraints.remove(\"Kirchhoff-Voltage-Law\")\n",
    "\n",
    "\n",
    "n.optimize(extra_functionality=remove_kvl)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following, we exemplarily present a set of additional constraints. Note, the dual values of the additional constraints won't be stored in default data fields in the `PyPSA` network. But in any case, they are stored in the `linopy.Model`. \n",
    "\n",
    "Again, we **first build** the optimization model, **add our constraints** and finally **solve the network**. For the first step, we use again our accessor `optimize` to access the function `create_model`. This returns the `linopy` model that we can modify."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = n.optimize.create_model()  # the return value is the model, let's use it directly!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. **Minimum for state of charge**\n",
    "\n",
    "Assume we want to set a minimum state of charge of 50 MWh in our storage unit. This is done by: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sus = m.variables[\"StorageUnit-state_of_charge\"]\n",
    "m.add_constraints(sus >= 50, name=\"StorageUnit-minimum_soc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The return value of the `add_constraints` function is a array with the labels of the constraints. You can access the constraint now through: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.constraints[\"StorageUnit-minimum_soc\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and inspects its attributes like `lhs`, `sign` and `rhs`, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.constraints[\"StorageUnit-minimum_soc\"].rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. **Fix the ratio between ingoing and outgoing capacity of the Store**\n",
    "\n",
    "The battery in our system is modelled with two links and a store. We should make sure that its charging and discharging capacities, meaning their links, are somehow coupled. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "capacity = m.variables[\"Link-p_nom\"]\n",
    "eff = n.links.at[\"battery_power\", \"efficiency\"]\n",
    "lhs = capacity.loc[\"battery_power\"] - eff * capacity.loc[\"battery_discharge\"]\n",
    "m.add_constraints(lhs == 0, name=\"Link-battery_fix_ratio\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. **Every bus must in total produce the 20% of the total demand**\n",
    "\n",
    "For this, we use the linopy function `groupby_sum` which follows the pattern from `pandas`/`xarray`'s `groupby` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_demand = n.loads_t.p_set.sum().sum()\n",
    "buses = n.generators.bus.to_xarray()\n",
    "prod_per_bus = m.variables[\"Generator-p\"].groupby(buses).sum().sum(\"snapshot\")\n",
    "m.add_constraints(prod_per_bus >= total_demand / 5, name=\"Bus-minimum_production_share\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = prod_per_bus >= total_demand / 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and now let's solve the network again. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.optimize.solve_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing the constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if the system got our own constraints. We look at `n.constraints` which combines summarises constraints going into the linear problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model.constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last three entries show our constraints. Let's check whether out two custom constraint are fulfilled:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.links.loc[[\"battery_power\", \"battery_discharge\"], [\"p_nom_opt\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.storage_units_t.state_of_charge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.generators_t.p.groupby(n.generators.bus, axis=1).sum().sum() / n.loads_t.p.sum().sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks good! Now, let's see which dual values were parsed. Therefore we have a look into `n.model.dual` \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model.dual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model.dual[\"StorageUnit-minimum_soc\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n.model.dual[\"Link-battery_fix_ratio\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n.model.dual[\"Bus-minimum_production_share\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are the basic functionalities of the `optimize` accessor. There are many more functions like abstract optimziation formulations (security constraint optimization, iterative transmission expansion optimization, etc.) or helper functions (fixing optimized capacities, adding load shedding). Try them out if you want!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"\\n\".join([func for func in n.optimize.__dir__() if not func.startswith(\"_\")]))"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "",
   "language": "python",
   "name": ""
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}