convert.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from otoole import ReadCsv\n",
    "import os\n",
    "\n",
    "\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "from numpy import inf\n",
    "\n",
    "from otoole.utils import (\n",
    "    _read_file,\n",
    "    validate_config\n",
    ")\n",
    "\n",
    "import logging\n",
    "\n",
    "from dask.distributed import Client\n",
    "\n",
    "from variables import add_variables\n",
    "from constraints import add_demand_constraints\n",
    "from validate import validate_user_data\n",
    "\n",
    "logger = logging.getLogger(__name__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_name = 'model_one'\n",
    "config_path = os.path.join('test', f'{model_name}.yaml')\n",
    "folder_path = os.path.join('test', model_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(f'{model_name}.nc'):\n",
    "    with open(config_path, \"r\") as config_file:\n",
    "        config = _read_file(config_file, '.yaml')\n",
    "    validate_config(config)\n",
    "\n",
    "    read_strategy = ReadCsv(user_config=config)\n",
    "\n",
    "    model, defaults = read_strategy.read(folder_path)\n",
    "\n",
    "    # Reshape the CSV files to create the dataset of sets and parameters\n",
    "    # If the sparse package is available (Python <3.11)\n",
    "    # data_vars = {x: xr.DataArray().from_series(y.VALUE, sparse=True) for x, y in model.items() if config[x]['type'] == 'param'}\n",
    "    data_vars = {x: y.VALUE.to_xarray() for x, y in model.items() if config[x]['type'] == 'param'}\n",
    "    coords = {x: y.values.T[0] for x, y in model.items() if config[x]['type'] == 'set'}\n",
    "    ds = xr.Dataset(data_vars=data_vars, coords=coords)\n",
    "    ds = ds.assign_coords({'_REGION': model['REGION'].values.T[0]})\n",
    "\n",
    "    for param, default in defaults.items():\n",
    "        if config[param]['type'] == 'param':\n",
    "            ds[param].attrs['default'] = default\n",
    "            if default != 0:\n",
    "                ds[param] = ds[param].fillna(default)\n",
    "    # ds.to_netcdf(f'{model_name}.nc', engine='netcdf4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ds = xr.open_dataset(f'{model_name}.nc')\n",
    "# ds = ds.drop_sel(YEAR=range(2025, 2071))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# client = Client()\n",
    "# client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from linopy import Model, solvers, available_solvers\n",
    "\n",
    "chunks = {\n",
    "    'YEAR': ds.YEAR.size,\n",
    "    'TIMESLICE': ds.TIMESLICE.size,\n",
    "    'REGION': 1,\n",
    "    'TECHNOLOGY': 100,\n",
    "    'FUEL': 100,\n",
    "    'MODE_OF_OPERATION': ds.MODE_OF_OPERATION.size,\n",
    "    'SEASON': ds.SEASON.size\n",
    "}\n",
    "\n",
    "m = Model(force_dim_names=True) #, chunk='auto')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = add_variables(ds, m)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discounting\n",
    "\n",
    "```ampl\n",
    "param DiscountRate{r in REGION};\n",
    "param DiscountRateIdv{r in REGION, t in TECHNOLOGY}, default DiscountRate[r];\n",
    "\n",
    "param DiscountFactor{r in REGION, y in YEAR} :=\n",
    "\t(1 + DiscountRate[r]) ^ (y - min{yy in YEAR} min(yy) + 0.0);\n",
    "param DiscountFactorMid{r in REGION, y in YEAR} :=\n",
    "\t(1 + DiscountRate[r]) ^ (y - min{yy in YEAR} min(yy) + 0.5);\n",
    "\n",
    "param OperationalLife{r in REGION, t in TECHNOLOGY};\n",
    "\n",
    "param CapitalRecoveryFactor{r in REGION, t in TECHNOLOGY} :=\n",
    "\t(1 - (1 + DiscountRateIdv[r,t])^(-1))/(1 - (1 + DiscountRateIdv[r,t])^(-(OperationalLife[r,t])));\n",
    "param PvAnnuity{r in REGION, t in TECHNOLOGY} :=\n",
    "\t(1 - (1 + DiscountRate[r])^(-(OperationalLife[r,t]))) * (1 + DiscountRate[r]) / DiscountRate[r];\n",
    "\n",
    "param DiscountRateStorage{r in REGION, s in STORAGE};\n",
    "param DiscountFactorStorage{r in REGION, s in STORAGE, y in YEAR} :=\n",
    "\t(1 + DiscountRateStorage[r, s]) ^ (y - min{yy in YEAR} min(yy) + 0.0);\n",
    "param DiscountFactorMidStorage{r in REGION, s in STORAGE, y in YEAR} :=\n",
    "\t(1 + DiscountRateStorage[r, s]) ^ (y - min{yy in YEAR} min(yy) + 0.5);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "discount_factor = ((1 + ds['DiscountRate']) ** (ds.coords['YEAR'] - min(ds.coords['YEAR'])))\n",
    "discount_factor_mid = ((1 + ds['DiscountRate']) ** (ds.coords['YEAR'] - min(ds.coords['YEAR']) + 0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "discount_factor_idv = ((1 + ds['DiscountRateIdv']) ** (ds.coords['YEAR'] - min(ds.coords['YEAR'])))\n",
    "discount_factor_mid_idv = ((1 + ds['DiscountRateIdv']) ** (ds.coords['YEAR'] - min(ds.coords['YEAR']) + 0.5))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pv_annuity = (1 - (1 + ds['DiscountRateIdv'])**(-(ds['OperationalLife']))) * (1 + ds['DiscountRateIdv']) / ds['DiscountRateIdv']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "capital_recovery_factor = (1 - (1 + ds['DiscountRateIdv'])**(-1))/(1 - (1 + ds['DiscountRateIdv'])**(-(ds['OperationalLife'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constraints"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Storage\n",
    "\n",
    "```ampl\n",
    "s.t. S1_RateOfStorageCharge{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:\n",
    "\tsum{t in TECHNOLOGY, m in MODE_OF_OPERATION, l in TIMESLICE: TechnologyToStorage[r,t,s,m] > 0}\n",
    "\tRateOfActivity[r,l,t,m,y] * TechnologyToStorage[r,t,s,m] * Conversionls[l,ls] * Conversionld[l,ld] * Conversionlh[l,lh]\n",
    "\t=\n",
    "\tRateOfStorageCharge[r,s,ls,ld,lh,y];\n",
    "\n",
    "s.t. S2_RateOfStorageDischarge{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:\n",
    "\tsum{t in TECHNOLOGY, m in MODE_OF_OPERATION, l in TIMESLICE: TechnologyFromStorage[r,t,s,m] > 0}\n",
    "\tRateOfActivity[r,l,t,m,y] * TechnologyFromStorage[r,t,s,m] * Conversionls[l,ls] * Conversionld[l,ld] * Conversionlh[l,lh]\n",
    "\t=\n",
    "\tRateOfStorageDischarge[r,s,ls,ld,lh,y];\n",
    "\n",
    "s.t. S3_NetChargeWithinYear{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:\n",
    "\tsum{l in TIMESLICE:Conversionls[l,ls]>0 && Conversionld[l,ld] > 0 && Conversionlh[l,lh] > 0}\n",
    "\t(RateOfStorageCharge[r,s,ls,ld,lh,y] - RateOfStorageDischarge[r,s,ls,ld,lh,y]) * YearSplit[l,y] *\n",
    "\tConversionls[l,ls] * Conversionld[l,ld] * Conversionlh[l,lh]\n",
    "\t=\n",
    "\tNetChargeWithinYear[r,s,ls,ld,lh,y];\n",
    "\n",
    "s.t. S4_NetChargeWithinDay{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:\n",
    "\t(RateOfStorageCharge[r,s,ls,ld,lh,y] - RateOfStorageDischarge[r,s,ls,ld,lh,y]) * DaySplit[lh,y]\n",
    "\t=\n",
    "\tNetChargeWithinDay[r,s,ls,ld,lh,y];\n",
    "\n",
    "s.t. S5_and_S6_StorageLevelYearStart{r in REGION, s in STORAGE, y in YEAR}:\n",
    "\tif y = min{yy in YEAR} min(yy)\n",
    "\tthen StorageLevelStart[r,s]\n",
    "\telse StorageLevelYearStart[r,s,y-1] + sum{ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET}\n",
    "\tNetChargeWithinYear[r,s,ls,ld,lh,y-1]\n",
    "\t=\n",
    "\tStorageLevelYearStart[r,s,y];\n",
    "\n",
    "s.t. S7_and_S8_StorageLevelYearFinish{r in REGION, s in STORAGE, y in YEAR}:\n",
    "\tif y < max{yy in YEAR} max(yy)\n",
    "\tthen StorageLevelYearStart[r,s,y+1]\n",
    "\telse StorageLevelYearStart[r,s,y] + sum{ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET}\n",
    "\tNetChargeWithinYear[r,s,ls,ld,lh,y]\n",
    "\t=\n",
    "\tStorageLevelYearFinish[r,s,y];\n",
    "\n",
    "s.t. S9_and_S10_StorageLevelSeasonStart{r in REGION, s in STORAGE, ls in SEASON, y in YEAR}:\n",
    "\tif ls = min{lsls in SEASON} min(lsls)\n",
    "\tthen StorageLevelYearStart[r,s,y]\n",
    "\telse StorageLevelSeasonStart[r,s,ls-1,y] + sum{ld in DAYTYPE, lh in DAILYTIMEBRACKET}\n",
    "\tNetChargeWithinYear[r,s,ls-1,ld,lh,y]\n",
    "\t=\n",
    "\tStorageLevelSeasonStart[r,s,ls,y];\n",
    "\n",
    "s.t. S11_and_S12_StorageLevelDayTypeStart{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, y in YEAR}:\n",
    "\tif ld = min{ldld in DAYTYPE} min(ldld)\n",
    "\tthen StorageLevelSeasonStart[r,s,ls,y]\n",
    "\telse StorageLevelDayTypeStart[r,s,ls,ld-1,y] + sum{lh in DAILYTIMEBRACKET}\n",
    "\tNetChargeWithinDay[r,s,ls,ld-1,lh,y] * DaysInDayType[ls,ld-1,y]\n",
    "\t=\n",
    "\tStorageLevelDayTypeStart[r,s,ls,ld,y];\n",
    "\n",
    "s.t. S13_and_S14_and_S15_StorageLevelDayTypeFinish{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, y in YEAR}:\n",
    "\tif ls = max{lsls in SEASON} max(lsls) && ld = max{ldld in DAYTYPE} max(ldld)\n",
    "\tthen StorageLevelYearFinish[r,s,y]\n",
    "\telse if ld = max{ldld in DAYTYPE} max(ldld)\n",
    "\tthen StorageLevelSeasonStart[r,s,ls+1,y]\n",
    "\telse StorageLevelDayTypeFinish[r,s,ls,ld+1,y] - sum{lh in DAILYTIMEBRACKET}\n",
    "\tNetChargeWithinDay[r,s,ls,ld+1,lh,y] * DaysInDayType[ls,ld+1,y]\n",
    "\t=\n",
    "\tStorageLevelDayTypeFinish[r,s,ls,ld,y];\n",
    "\n",
    "#\n",
    "##########\t\tStorage Constraints\t\t\t\t#############\n",
    "#\n",
    "s.t. SC1_LowerLimit_BeginningOfDailyTimeBracketOfFirstInstanceOfDayTypeInFirstWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: 0 <= (StorageLevelDayTypeStart[r,s,ls,ld,y]+sum{lhlh in DAILYTIMEBRACKET:lh-lhlh>0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageLowerLimit[r,s,y];\n",
    "s.t. SC1_UpperLimit_BeginningOfDailyTimeBracketOfFirstInstanceOfDayTypeInFirstWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: (StorageLevelDayTypeStart[r,s,ls,ld,y]+sum{lhlh in DAILYTIMEBRACKET:lh-lhlh>0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageUpperLimit[r,s,y] <= 0;\n",
    "s.t. SC2_LowerLimit_EndOfDailyTimeBracketOfLastInstanceOfDayTypeInFirstWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: 0 <= if ld > min{ldld in DAYTYPE} min(ldld) then (StorageLevelDayTypeStart[r,s,ls,ld,y]-sum{lhlh in DAILYTIMEBRACKET:lh-lhlh<0} NetChargeWithinDay[r,s,ls,ld-1,lhlh,y])-StorageLowerLimit[r,s,y];\n",
    "s.t. SC2_UpperLimit_EndOfDailyTimeBracketOfLastInstanceOfDayTypeInFirstWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: if ld > min{ldld in DAYTYPE} min(ldld) then (StorageLevelDayTypeStart[r,s,ls,ld,y]-sum{lhlh in DAILYTIMEBRACKET:lh-lhlh<0} NetChargeWithinDay[r,s,ls,ld-1,lhlh,y])-StorageUpperLimit[r,s,y] <= 0;\n",
    "s.t. SC3_LowerLimit_EndOfDailyTimeBracketOfLastInstanceOfDayTypeInLastWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:  0 <= (StorageLevelDayTypeFinish[r,s,ls,ld,y] - sum{lhlh in DAILYTIMEBRACKET:lh-lhlh<0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageLowerLimit[r,s,y];\n",
    "s.t. SC3_UpperLimit_EndOfDailyTimeBracketOfLastInstanceOfDayTypeInLastWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}:  (StorageLevelDayTypeFinish[r,s,ls,ld,y] - sum{lhlh in DAILYTIMEBRACKET:lh-lhlh<0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageUpperLimit[r,s,y] <= 0;\n",
    "s.t. SC4_LowerLimit_BeginningOfDailyTimeBracketOfFirstInstanceOfDayTypeInLastWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: \t0 <= if ld > min{ldld in DAYTYPE} min(ldld) then (StorageLevelDayTypeFinish[r,s,ls,ld-1,y]+sum{lhlh in DAILYTIMEBRACKET:lh-lhlh>0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageLowerLimit[r,s,y];\n",
    "s.t. SC4_UpperLimit_BeginningOfDailyTimeBracketOfFirstInstanceOfDayTypeInLastWeekConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: if ld > min{ldld in DAYTYPE} min(ldld) then (StorageLevelDayTypeFinish[r,s,ls,ld-1,y]+sum{lhlh in DAILYTIMEBRACKET:lh-lhlh>0} NetChargeWithinDay[r,s,ls,ld,lhlh,y])-StorageUpperLimit[r,s,y] <= 0;\n",
    "s.t. SC5_MaxChargeConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: RateOfStorageCharge[r,s,ls,ld,lh,y] <= StorageMaxChargeRate[r,s];\n",
    "s.t. SC6_MaxDischargeConstraint{r in REGION, s in STORAGE, ls in SEASON, ld in DAYTYPE, lh in DAILYTIMEBRACKET, y in YEAR}: RateOfStorageDischarge[r,s,ls,ld,lh,y] <= StorageMaxDischargeRate[r,s];\n",
    "#\n",
    "#########\t\tStorage Investments\t\t\t\t#############\n",
    "#\n",
    "s.t. SI1_StorageUpperLimit{r in REGION, s in STORAGE, y in YEAR}: AccumulatedNewStorageCapacity[r,s,y]+ResidualStorageCapacity[r,s,y] = StorageUpperLimit[r,s,y];\n",
    "s.t. SI2_StorageLowerLimit{r in REGION, s in STORAGE, y in YEAR}: MinStorageCharge[r,s,y]*StorageUpperLimit[r,s,y] = StorageLowerLimit[r,s,y];\n",
    "s.t. SI3_TotalNewStorage{r in REGION, s in STORAGE, y in YEAR}: sum{yy in YEAR: y-yy < OperationalLifeStorage[r,s] && y-yy>=0} NewStorageCapacity[r,s,yy]=AccumulatedNewStorageCapacity[r,s,y];\n",
    "s.t. SI4_UndiscountedCapitalInvestmentStorage{r in REGION, s in STORAGE, y in YEAR}: CapitalCostStorage[r,s,y] * NewStorageCapacity[r,s,y] = CapitalInvestmentStorage[r,s,y];\n",
    "s.t. SI5_DiscountingCapitalInvestmentStorage{r in REGION, s in STORAGE, y in YEAR}: CapitalInvestmentStorage[r,s,y]/(DiscountFactorStorage[r,s,y]) = DiscountedCapitalInvestmentStorage[r,s,y];\n",
    "s.t. SI6_SalvageValueStorageAtEndOfPeriod1{r in REGION, s in STORAGE, y in YEAR: (y+OperationalLifeStorage[r,s]-1) <= (max{yy in YEAR} max(yy))}: 0 = SalvageValueStorage[r,s,y];\n",
    "s.t. SI7_SalvageValueStorageAtEndOfPeriod2{r in REGION, s in STORAGE, y in YEAR: (DepreciationMethod[r]=1 && (y+OperationalLifeStorage[r,s]-1) > (max{yy in YEAR} max(yy)) && DiscountRateStorage[r,s]=0) || (DepreciationMethod[r]=2 && (y+OperationalLifeStorage[r,s]-1) > (max{yy in YEAR} max(yy)))}: CapitalInvestmentStorage[r,s,y]*(1-(max{yy in YEAR} max(yy) - y+1)/OperationalLifeStorage[r,s]) = SalvageValueStorage[r,s,y];\n",
    "s.t. SI8_SalvageValueStorageAtEndOfPeriod3{r in REGION, s in STORAGE, y in YEAR: DepreciationMethod[r]=1 && (y+OperationalLifeStorage[r,s]-1) > (max{yy in YEAR} max(yy)) && DiscountRateStorage[r,s]>0}: CapitalInvestmentStorage[r,s,y]*(1-(((1+DiscountRateStorage[r,s])^(max{yy in YEAR} max(yy) - y+1)-1)/((1+DiscountRateStorage[r,s])^OperationalLifeStorage[r,s]-1))) = SalvageValueStorage[r,s,y];\n",
    "s.t. SI9_SalvageValueStorageDiscountedToStartYear{r in REGION, s in STORAGE, y in YEAR}: SalvageValueStorage[r,s,y]/((1+DiscountRateStorage[r,s])^(max{yy in YEAR} max(yy)-min{yy in YEAR} min(yy)+1)) = DiscountedSalvageValueStorage[r,s,y];\n",
    "s.t. SI10_TotalDiscountedCostByStorage{r in REGION, s in STORAGE, y in YEAR}: DiscountedCapitalInvestmentStorage[r,s,y]-DiscountedSalvageValueStorage[r,s,y] = TotalDiscountedStorageCost[r,s,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Demand"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EQ_SpecifiedDemand{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\tSpecifiedAnnualDemand[r,f,y] <> 0}:\n",
    "\tSpecifiedAnnualDemand[r,f,y] * SpecifiedDemandProfile[r,f,l,y] / YearSplit[l,y]\n",
    "\t=\n",
    "\tRateOfDemand[r,l,f,y];\n",
    "\t```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mask = ~ds['SpecifiedAnnualDemand'].isnull()\n",
    "# con = m['RateOfDemand'] == (ds['SpecifiedAnnualDemand'] * ds['SpecifiedDemandProfile'] / ds['YearSplit'])\n",
    "# m.add_constraints(con, name='EQ_SpecifiedDemand', mask=mask)\n",
    "m = add_demand_constraints(ds, m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Capacity Adequacy A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CAa1_TotalNewCapacity{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tAccumulatedNewCapacity[r,t,y]\n",
    "\t=\n",
    "\tsum{yy in YEAR: y-yy < OperationalLife[r,t] && y - yy >= 0} NewCapacity[r,t,yy];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cap = m['NewCapacity'].rename(YEAR=\"BUILDYEAR\")\n",
    "mask = (ds.YEAR - new_cap.data.BUILDYEAR >= 0) & (ds.YEAR - new_cap.data.BUILDYEAR < ds.OperationalLife)\n",
    "con = m['AccumulatedNewCapacity'] - new_cap.where(mask).sum(\"BUILDYEAR\") == 0\n",
    "m.add_constraints(con, name='CAa1_TotalNewCapacity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CAa2_TotalAnnualCapacity{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tAccumulatedNewCapacity[r,t,y] + ResidualCapacity[r,t,y]\n",
    "\t=\n",
    "\tTotalCapacityAnnual[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalCapacityAnnual'] - m['AccumulatedNewCapacity'] == ds['ResidualCapacity'].fillna(0)\n",
    "m.add_constraints(con, name=\"CAa2_TotalAnnualCapacity\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CAa3_TotalActivityOfEachTechnology{r in REGION, t in TECHNOLOGY, l in TIMESLICE, y in YEAR}:\n",
    "\tsum{m in MODE_OF_OPERATION} RateOfActivity[r,l,t,m,y]\n",
    "\t=\n",
    "\tRateOfTotalActivity[r,t,l,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfActivity'].sum(dims='MODE_OF_OPERATION') - m['RateOfTotalActivity'] == 0\n",
    "m.add_constraints(con, name='CAa3_TotalActivityOfEachTechnology')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CAa4_Constraint_Capacity{r in REGION, l in TIMESLICE, t in TECHNOLOGY, y in YEAR}:\n",
    "\tRateOfTotalActivity[r,t,l,y]\n",
    "\t<=\n",
    "\tTotalCapacityAnnual[r,t,y] * CapacityFactor[r,t,l,y] * CapacityToActivityUnit[r,t];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfTotalActivity'] - (m['TotalCapacityAnnual'] * ds['CapacityFactor'] * ds['CapacityToActivityUnit']) <=0\n",
    "mask = ~ds['CapacityFactor'].isnull()\n",
    "m.add_constraints(con, name='CAa4_Constraint_Capacity', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CAa5_TotalNewCapacity{r in REGION, t in TECHNOLOGY, y in YEAR: CapacityOfOneTechnologyUnit[r,t,y]<>0}:\n",
    "\tCapacityOfOneTechnologyUnit[r,t,y] * NumberOfNewTechnologyUnits[r,t,y]\n",
    "\t=\n",
    "\tNewCapacity[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = ds['CapacityOfOneTechnologyUnit'] * m['NumberOfNewTechnologyUnits'] - m['NewCapacity'] == 0\n",
    "mask = ds['CapacityOfOneTechnologyUnit'].notnull()\n",
    "m.add_constraints(con, name='CAa5_TotalNewCapacity', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Capacity Adequacy B\n",
    "```ampl\n",
    "s.t. CAb1_PlannedMaintenance{r in REGION, t in TECHNOLOGY, y in YEAR: AvailabilityFactor[r,t,y] < 1}:\n",
    "\tsum{l in TIMESLICE} RateOfTotalActivity[r,t,l,y] * YearSplit[l,y]\n",
    "\t<=\n",
    "\tsum{l in TIMESLICE} (TotalCapacityAnnual[r,t,y] * CapacityFactor[r,t,l,y] * YearSplit[l,y])\n",
    "\t* AvailabilityFactor[r,t,y] * CapacityToActivityUnit[r,t];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ds['AvailabilityFactor'] < 1\n",
    "con = (m['RateOfTotalActivity'] * ds['YearSplit']).sum(dims='TIMESLICE') \\\n",
    "      - ((m['TotalCapacityAnnual'] * ds['CapacityFactor'] * ds['YearSplit']).sum(dims='TIMESLICE') * ds['AvailabilityFactor'] * ds['CapacityToActivityUnit']) <= 0\n",
    "m.add_constraints(con, name='CAb1_PlannedMaintenance', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy Balance A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa1_RateOfFuelProduction1{\n",
    "\tr in REGION, l in TIMESLICE, f in FUEL, t in TECHNOLOGY, m in MODE_OF_OPERATION, y in YEAR:\n",
    "\tOutputActivityRatio[r,t,f,m,y] <> 0}:\n",
    "\tRateOfActivity[r,l,t,m,y] * OutputActivityRatio[r,t,f,m,y]\n",
    "\t=\n",
    "\tRateOfProductionByTechnologyByMode[r,l,t,m,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ds['OutputActivityRatio'].notnull()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfActivity'] * ds['OutputActivityRatio'] - m['RateOfProductionByTechnologyByMode'] == 0\n",
    "m.add_constraints(con, name='EBa1_RateOfFuelProduction1', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa2_RateOfFuelProduction2{r in REGION, l in TIMESLICE, f in FUEL, t in TECHNOLOGY, y in YEAR}:\n",
    "\tsum{m in MODE_OF_OPERATION: OutputActivityRatio[r,t,f,m,y] <> 0} RateOfProductionByTechnologyByMode[r,l,t,m,f,y]\n",
    "\t=\n",
    "\tRateOfProductionByTechnology[r,l,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ds['OutputActivityRatio'].sum('MODE_OF_OPERATION') != 0\n",
    "con = m['RateOfProductionByTechnologyByMode'].where(mask).sum(dims='MODE_OF_OPERATION') - m['RateOfProductionByTechnology'] == 0\n",
    "m.add_constraints(con, name='EBa2_RateOfFuelProduction2', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa3_RateOfFuelProduction3{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t\t    (sum{t in TECHNOLOGY, m in MODE_OF_OPERATION} OutputActivityRatio[r,t,f,m,y]) <> 0}:\n",
    "\tsum{t in TECHNOLOGY} RateOfProductionByTechnology[r,l,t,f,y]\n",
    "\t=\n",
    "\tRateOfProduction[r,l,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfProductionByTechnology'].sum(dims='TECHNOLOGY') - m['RateOfProduction'] == 0\n",
    "mask = ds['OutputActivityRatio'].sum(dim=['TECHNOLOGY', 'MODE_OF_OPERATION']) != 0\n",
    "m.add_constraints(con, name='EBa3_RateOfFuelProduction3', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa4_RateOfFuelUse1{r in REGION, l in TIMESLICE, f in FUEL, t in TECHNOLOGY, m in MODE_OF_OPERATION, y in YEAR:\n",
    "\t\t\t\t\t\t InputActivityRatio[r,t,f,m,y] <> 0}:\n",
    "\tRateOfActivity[r,l,t,m,y] * InputActivityRatio[r,t,f,m,y]\n",
    "\t=\n",
    "\tRateOfUseByTechnologyByMode[r,l,t,m,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfActivity'] * ds['InputActivityRatio'] - m['RateOfUseByTechnologyByMode'] == 0\n",
    "mask = ds['InputActivityRatio'].notnull()\n",
    "m.add_constraints(con, name='EBa4_RateOfFuelUse1', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa5_RateOfFuelUse2{r in REGION, l in TIMESLICE, f in FUEL, t in TECHNOLOGY, y in YEAR:\n",
    "\t\t\t\t\t\t sum{m in MODE_OF_OPERATION} InputActivityRatio[r,t,f,m,y] <> 0}:\n",
    "\tsum{m in MODE_OF_OPERATION: InputActivityRatio[r,t,f,m,y] <> 0}\n",
    "\tRateOfUseByTechnologyByMode[r,l,t,m,f,y]\n",
    "\t=\n",
    "\tRateOfUseByTechnology[r,l,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfUseByTechnologyByMode'].sum(dims='MODE_OF_OPERATION') - m['RateOfUseByTechnology'] == 0\n",
    "mask = ds['InputActivityRatio'].sum(dim='MODE_OF_OPERATION') != 0\n",
    "m.add_constraints(con, name='EBa5_RateOfFuelUse2', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa6_RateOfFuelUse3{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t sum{t in TECHNOLOGY, m in MODE_OF_OPERATION} InputActivityRatio[r,t,f,m,y] <> 0}:\n",
    "\tsum{t in TECHNOLOGY} RateOfUseByTechnology[r,l,t,f,y]\n",
    "\t=\n",
    "\tRateOfUse[r,l,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfUseByTechnology'].sum(dims='TECHNOLOGY') - m['RateOfUse'] == 0\n",
    "mask = ds['InputActivityRatio'].sum(dim=['TECHNOLOGY', 'MODE_OF_OPERATION']) != 0\n",
    "m.add_constraints(con, name='EBa6_RateOfFuelUse3', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa7_EnergyBalanceEachTS1{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t\t   (sum{t in TECHNOLOGY, m in MODE_OF_OPERATION} OutputActivityRatio[r,t,f,m,y]) <> 0}:\n",
    "\tRateOfProduction[r,l,f,y] * YearSplit[l,y]\n",
    "\t=\n",
    "\tProduction[r,l,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfProduction'] * ds['YearSplit']) - m['Production'] == 0\n",
    "mask = ds['OutputActivityRatio'].sum(dim=['TECHNOLOGY', 'MODE_OF_OPERATION']) != 0\n",
    "m.add_constraints(con, name='EBa7_EnergyBalanceEachTS1', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa8_EnergyBalanceEachTS2{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t\t   (sum{t in TECHNOLOGY, m in MODE_OF_OPERATION} InputActivityRatio[r,t,f,m,y]) <> 0}:\n",
    "\tRateOfUse[r,l,f,y] * YearSplit[l,y]\n",
    "\t=\n",
    "\tUse[r,l,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfUse'] * ds['YearSplit']) - m['Use'] == 0\n",
    "mask = ds['InputActivityRatio'].sum(dim=['TECHNOLOGY', 'MODE_OF_OPERATION']) != 0\n",
    "m.add_constraints(con, name='EBa8_EnergyBalanceEachTS2', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa9_EnergyBalanceEachTS3{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t\t   SpecifiedAnnualDemand[r,f,y] <> 0}:\n",
    "\tRateOfDemand[r,l,f,y] * YearSplit[l,y]\n",
    "\t=\n",
    "\tDemand[r,l,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfDemand'] * ds['YearSplit']) - m['Demand'] == 0\n",
    "mask = ds['SpecifiedAnnualDemand'].notnull()\n",
    "m.add_constraints(con, name='EBa9_EnergyBalanceEachTS3', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa10_EnergyBalanceEachTS4{r in REGION, rr in REGION, l in TIMESLICE, f in FUEL, y in YEAR:\n",
    "\t\t\t\t\t\t\t\tTradeRoute[r,rr,f,y] <> 0}:\n",
    "\tTrade[r,rr,l,f,y]\n",
    "\t=\n",
    "\t-Trade[rr,r,l,f,y];\n",
    "```\n",
    "\n",
    "```python\n",
    "def EnergyBalanceEachTS4_rule(model, r, rr, l, f, y):\n",
    "    return model.Trade[r, rr, l, f, y] + model.Trade[rr, r, l, f, y] == 0\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr = ds['TradeRoute']\n",
    "from_mask = tr.where(tr.REGION != tr._REGION).notnull()\n",
    "to_mask = tr.where(tr._REGION != tr.REGION).notnull()\n",
    "con = m['Trade'].where(from_mask) + m['Trade'].where(to_mask) == 0\n",
    "m.add_constraints(con, name='EBa10_EnergyBalanceEachTS4')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBa11_EnergyBalanceEachTS5{r in REGION, l in TIMESLICE, f in FUEL, y in YEAR}:\n",
    "\tProduction[r,l,f,y]\n",
    "\t>=\n",
    "\tDemand[r,l,f,y] + Use[r,l,f,y] + sum{rr in REGION} Trade[r,rr,l,f,y] * TradeRoute[r,rr,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['Production'] - (m['Demand'] + m['Use'] + (m['Trade'] * ds['TradeRoute'].fillna(0)).sum(dims='_REGION')) >= 0\n",
    "m.add_constraints(con, name='EBa11_EnergyBalanceEachTS5')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy Balance B"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBb1_EnergyBalanceEachYear1{r in REGION, f in FUEL, y in YEAR}:\n",
    "\tsum{l in TIMESLICE} Production[r,l,f,y]\n",
    "\t=\n",
    "\tProductionAnnual[r,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['Production'].sum(dims='TIMESLICE') - m['ProductionAnnual'] == 0\n",
    "m.add_constraints(con, name='EBb1_EnergyBalanceEachYear1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBb2_EnergyBalanceEachYear2{r in REGION, f in FUEL, y in YEAR}:\n",
    "\tsum{l in TIMESLICE} Use[r,l,f,y]\n",
    "\t=\n",
    "\tUseAnnual[r,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['Use'].sum('TIMESLICE') - m['UseAnnual'] == 0\n",
    "m.add_constraints(con, name='EBb2_EnergyBalanceEachYear2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBb3_EnergyBalanceEachYear3{r in REGION, rr in REGION, f in FUEL, y in YEAR}:\n",
    "\tsum{l in TIMESLICE} Trade[r,rr,l,f,y]\n",
    "\t=\n",
    "\tTradeAnnual[r,rr,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['Trade'].sum('TIMESLICE') - m['TradeAnnual'] == 0\n",
    "mask = ds.coords['REGION'] != ds.coords['_REGION']\n",
    "m.add_constraints(con, name='EBb3_EnergyBalanceEachYear3', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. EBb4_EnergyBalanceEachYear4{r in REGION, f in FUEL, y in YEAR}:\n",
    "\tProductionAnnual[r,f,y]\n",
    "\t>=\n",
    "\tUseAnnual[r,f,y] + sum{rr in REGION} TradeAnnual[r,rr,f,y] * TradeRoute[r,rr,f,y] + AccumulatedAnnualDemand[r,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['ProductionAnnual'] - m['UseAnnual'] - (m['TradeAnnual'].sum('_REGION') * ds['TradeRoute'].sum('_REGION')) >= ds['AccumulatedAnnualDemand'].fillna(0)\n",
    "m.add_constraints(con, name='EBb4_EnergyBalanceEachYear4')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accounting Technology Production/Use"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. Acc1_FuelProductionByTechnology{r in REGION, l in TIMESLICE, t in TECHNOLOGY, f in FUEL, y in YEAR}:\n",
    "\tRateOfProductionByTechnology[r,l,t,f,y] * YearSplit[l,y]\n",
    "\t=\n",
    "\tProductionByTechnology[r,l,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfProductionByTechnology'] * ds['YearSplit']) - m['ProductionByTechnology'] == 0\n",
    "mask = ds['OutputActivityRatio'].sum('MODE_OF_OPERATION') != 0\n",
    "m.add_constraints(con, name='Acc1_FuelProductionByTechnology', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. Acc2_FuelUseByTechnology{r in REGION, l in TIMESLICE, t in TECHNOLOGY, f in FUEL, y in YEAR}:\n",
    "\tRateOfUseByTechnology[r,l,t,f,y] * YearSplit[l,y]\n",
    "\t=\n",
    "\tUseByTechnology[r,l,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfUseByTechnology'] * ds['YearSplit']) - m['UseByTechnology'] == 0\n",
    "mask = ds['InputActivityRatio'].sum('MODE_OF_OPERATION') != 0\n",
    "m.add_constraints(con, name='Acc2_FuelUseByTechnology', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. Acc3_AverageAnnualRateOfActivity{r in REGION, t in TECHNOLOGY, m in MODE_OF_OPERATION, y in YEAR}:\n",
    "\tsum{l in TIMESLICE} RateOfActivity[r,l,t,m,y]*YearSplit[l,y]\n",
    "\t=\n",
    "\tTotalAnnualTechnologyActivityByMode[r,t,m,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfActivity'] * ds['YearSplit']).sum('TIMESLICE') - m['TotalAnnualTechnologyActivityByMode'] == 0\n",
    "# mask = ds['OutputActivityRatio'].sum('FUEL') != 0\n",
    "m.add_constraints(con, name='Acc3_AverageAnnualRateOfActivity')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. Acc4_ModelPeriodCostByRegion{r in REGION}:\n",
    "\tsum{y in YEAR}TotalDiscountedCost[r,y] = ModelPeriodCostByRegion[r];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalDiscountedCost'].sum('YEAR') - m['ModelPeriodCostByRegion'] == 0\n",
    "m.add_constraints(con, name='Acc4_ModelPeriodCostByRegion')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Capital Costs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CC1_UndiscountedCapitalInvestment{r in REGION, t in TECHNOLOGY, y in YEAR}: \n",
    "        CapitalCost[r,t,y] * NewCapacity[r,t,y] * CapitalRecoveryFactor[r,t] * PvAnnuity[r,t] \n",
    "        = \n",
    "        CapitalInvestment[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = ds['CapitalCost'].fillna(0) * m['NewCapacity'] * capital_recovery_factor * pv_annuity - m['CapitalInvestment'] == 0\n",
    "m.add_constraints(con, name='CC1_UndiscountedCapitalInvestment')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. CC2_DiscountingCapitalInvestment{r in REGION, t in TECHNOLOGY, y in YEAR}: \n",
    "    CapitalInvestment[r,t,y]  / DiscountFactor[r,y] = DiscountedCapitalInvestment[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['CapitalInvestment'] / discount_factor) - m['DiscountedCapitalInvestment'] == 0\n",
    "m.add_constraints(con, name='CC2_DiscountingCapitalInvestment')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Salvage Value\n",
    "\n",
    "### GNU MathProg implementation\n",
    "```ampl\n",
    "s.t. SV1_SalvageValueAtEndOfPeriod1{\n",
    "    r in REGION, t in TECHNOLOGY, y in YEAR: \n",
    "        DepreciationMethod[r]=1 && \n",
    "        (y + OperationalLife[r,t]-1) > (max{yy in YEAR} max(yy)) && \n",
    "        DiscountRate[r]>0}: \n",
    "    SalvageValue[r,t,y] \n",
    "    = \n",
    "    CapitalCost[r,t,y] * NewCapacity[r,t,y] * CapitalRecoveryFactor[r,t] * PvAnnuity[r,t] * \n",
    "    (1-(((1+DiscountRate[r])^(max{yy in YEAR} max(yy) - y+1)-1)/((1+DiscountRate[r])^OperationalLife[r,t]-1)));\n",
    "```\n",
    "### Pyomo implementation\n",
    "```python\n",
    "def SalvageValueAtEndOfPeriod1_rule(model, r, t, y):\n",
    "    if (\n",
    "        model.DepreciationMethod[r] == 1\n",
    "        and ((y + model.OperationalLife[r, t] - 1) > max(model.YEAR))\n",
    "        and model.DiscountRate[r] > 0\n",
    "    ):\n",
    "        return model.SalvageValue[r, t, y] == model.CapitalCost[\n",
    "            r, t, y\n",
    "        ] * model.NewCapacity[r, t, y] * (\n",
    "            1\n",
    "            - (\n",
    "                ((1 + model.DiscountRate[r]) ** (max(model.YEAR) - y + 1) - 1)\n",
    "                / ((1 + model.DiscountRate[r]) ** model.OperationalLife[r, t] - 1)\n",
    "            )\n",
    "        )\n",
    "    elif (\n",
    "        model.DepreciationMethod[r] == 1\n",
    "        and ((y + model.OperationalLife[r, t] - 1) > max(model.YEAR))\n",
    "        and model.DiscountRate[r] == 0\n",
    "    ) or (\n",
    "        model.DepreciationMethod[r] == 2\n",
    "        and (y + model.OperationalLife[r, t] - 1) > (max(model.YEAR))\n",
    "    ):\n",
    "        return model.SalvageValue[r, t, y] == model.CapitalCost[\n",
    "            r, t, y\n",
    "        ] * model.NewCapacity[r, t, y] * (\n",
    "            1 - (max(model.YEAR) - y + 1) / model.OperationalLife[r, t]\n",
    "        )\n",
    "    else:\n",
    "        return model.SalvageValue[r, t, y] == 0\n",
    "\n",
    "\n",
    "model.SalvageValueAtEndOfPeriod1 = Constraint(\n",
    "    model.REGION, model.TECHNOLOGY, model.YEAR, rule=SalvageValueAtEndOfPeriod1_rule\n",
    ")\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def numerator(y: int):\n",
    "    return ((1 + ds['DiscountRateIdv']) ** (max(ds.coords['YEAR']) - y + 1) - 1)\n",
    "\n",
    "def denominator():\n",
    "    return ((1 + ds['DiscountRateIdv']) ** ds['OperationalLife'] - 1)\n",
    "\n",
    "def salvage_cost(ds):\n",
    "    return ds['CapitalCost'].fillna(0) * (1 - (numerator(ds.coords['YEAR']) / denominator()))\n",
    "\n",
    "con = m['SalvageValue'] - (m['NewCapacity'] * salvage_cost(ds)) == 0\n",
    "mask = (ds['DepreciationMethod'] == 1) & ((ds.coords['YEAR'] + ds['OperationalLife'] - 1) > max(ds.coords['YEAR'])) & (ds['DiscountRateIdv'] > 0)\n",
    "m.add_constraints(con, name='SV1_SalvageValueAtEndOfPeriod1', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. SV2_SalvageValueAtEndOfPeriod2{r in REGION, t in TECHNOLOGY, y in YEAR: \n",
    "        (DepreciationMethod[r]=1 && \n",
    "        (y + OperationalLife[r,t]-1) > (max{yy in YEAR} max(yy)) && \n",
    "        DiscountRate[r]=0) \n",
    "        || (DepreciationMethod[r]=2 && \n",
    "        (y + OperationalLife[r,t]-1) > (max{yy in YEAR} max(yy)))}: \n",
    "    SalvageValue[r,t,y] = CapitalCost[r,t,y] * NewCapacity[r,t,y] * CapitalRecoveryFactor[r,t] * PvAnnuity[r,t] *(1-(max{yy in YEAR} max(yy) - y+1)/OperationalLife[r,t]);\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def numerator(y: int):\n",
    "    return 1 - (max(ds.coords['YEAR']) - y + 1) - 1\n",
    "\n",
    "def denominator():\n",
    "    return ds['OperationalLife']\n",
    "\n",
    "def salvage_cost(ds):\n",
    "    return ds['CapitalCost'].fillna(0) * (1 - (numerator(ds.coords['YEAR']) / denominator()))\n",
    "\n",
    "con = m['SalvageValue'] - (m['NewCapacity'] * salvage_cost(ds)) == 0\n",
    "mask = ((ds['DepreciationMethod'] == 1) & ((ds.coords['YEAR'] + ds['OperationalLife'] - 1) > max(ds.coords['YEAR'])) & (ds['DiscountRateIdv'] == 0)) | ((ds['DepreciationMethod'] == 2) & ((ds.coords['YEAR'] + ds['OperationalLife'] - 1) > max(ds.coords['YEAR'])))\n",
    "m.add_constraints(con, name='SV2_SalvageValueAtEndOfPeriod2', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. SV3_SalvageValueAtEndOfPeriod3{r in REGION, t in TECHNOLOGY, y in YEAR: (y + OperationalLife[r,t]-1) <= (max{yy in YEAR} max(yy))}: \n",
    "    SalvageValue[r,t,y] = 0;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['SalvageValue'] == 0\n",
    "mask = ((ds.coords['YEAR'] + ds['OperationalLife'] - 1) <= max(ds.coords['YEAR']))\n",
    "m.add_constraints(con, name='SV3_SalvageValueAtEndOfPeriod3', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. SV4_SalvageValueDiscountedToStartYear{r in REGION, t in TECHNOLOGY, y in YEAR}: \n",
    "    DiscountedSalvageValue[r,t,y] = SalvageValue[r,t,y]/((1+DiscountRate[r])^(1+max{yy in YEAR} max(yy)-min{yy in YEAR} min(yy)));\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def discounting(ds):\n",
    "    return (1 + ds['DiscountRateIdv']) ** (1 + max(ds.coords['YEAR']) - min(ds.coords['YEAR']))\n",
    "con = m['DiscountedSalvageValue'] - m['SalvageValue'] / discounting(ds) == 0\n",
    "m.add_constraints(con, name='SV4_SalvageValueDiscountedToStartYear')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operating Costs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. OC1_OperatingCostsVariable{r in REGION, t in TECHNOLOGY, y in YEAR: sum{m in MODE_OF_OPERATION} VariableCost[r,t,m,y] <> 0}:\n",
    "\tsum{m in MODE_OF_OPERATION}\n",
    "\tTotalAnnualTechnologyActivityByMode[r,t,m,y] * VariableCost[r,t,m,y]\n",
    "\t=\n",
    "\tAnnualVariableOperatingCost[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['TotalAnnualTechnologyActivityByMode'] * ds['VariableCost'].fillna(0)).sum(dims='MODE_OF_OPERATION') - m['AnnualVariableOperatingCost'] == 0\n",
    "mask = (ds['VariableCost'].sum(dim='MODE_OF_OPERATION') != 0) & (~ds['VariableCost'].sum(dim='MODE_OF_OPERATION').isnull())\n",
    "m.add_constraints(con, name='OC1_OperatingCostsVariable', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. OC2_OperatingCostsFixedAnnual{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tTotalCapacityAnnual[r,t,y]*FixedCost[r,t,y]\n",
    "\t=\n",
    "\tAnnualFixedOperatingCost[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['TotalCapacityAnnual'] * ds['FixedCost'].fillna(0)) - m['AnnualFixedOperatingCost'] == 0\n",
    "# mask = ~ds['FixedCost'].isnull()\n",
    "m.add_constraints(con, name='OC2_OperatingCostsFixedAnnual')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. OC3_OperatingCostsTotalAnnual{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tAnnualFixedOperatingCost[r,t,y] + AnnualVariableOperatingCost[r,t,y]\n",
    "\t=\n",
    "\tOperatingCost[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['AnnualFixedOperatingCost'] + m['AnnualVariableOperatingCost'] - m['OperatingCost'] == 0\n",
    "# mask = (ds['VariableCost'].sum(dim='MODE_OF_OPERATION') != 0) & (~ds['FixedCost'].isnull())\n",
    "m.add_constraints(con, name='OC3_OperatingCostsTotalAnnual')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. OC4_DiscountedOperatingCostsTotalAnnual{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tOperatingCost[r,t,y] / DiscountFactorMid[r, y]\n",
    "\t=\n",
    "\tDiscountedOperatingCost[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['OperatingCost'] / discount_factor_mid - m['DiscountedOperatingCost'] == 0\n",
    "m.add_constraints(con, name='OC4_DiscountedOperatingCostsTotalAnnual')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total Discounted Costs\n",
    "\n",
    "```ampl\n",
    "s.t. TDC1_TotalDiscountedCostByTechnology{r in REGION, t in TECHNOLOGY, y in YEAR}: \n",
    "    DiscountedOperatingCost[r,t,y] + DiscountedCapitalInvestment[r,t,y] \n",
    "    + DiscountedTechnologyEmissionsPenalty[r,t,y] - DiscountedSalvageValue[r,t,y] = TotalDiscountedCostByTechnology[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['DiscountedOperatingCost'] + m['DiscountedCapitalInvestment'] + m['DiscountedTechnologyEmissionsPenalty'] \\\n",
    "        - m['DiscountedSalvageValue'] - m['TotalDiscountedCostByTechnology'] == 0\n",
    "m.add_constraints(con, name='TDC1_TotalDiscountedCostByTechnology')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. TDC2_TotalDiscountedCost{r in REGION, y in YEAR}: \n",
    "    sum{t in TECHNOLOGY} TotalDiscountedCostByTechnology[r,t,y] + \n",
    "    sum{s in STORAGE} TotalDiscountedStorageCost[r,s,y] \n",
    "    = TotalDiscountedCost[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    con = m['TotalDiscountedCostByTechnology'].sum('TECHNOLOGY') + m['TotalDiscountedStorageCost'].sum('STORAGE') - m['TotalDiscountedCost'] == 0\n",
    "except KeyError as ex:\n",
    "    con = m['TotalDiscountedCostByTechnology'].sum('TECHNOLOGY') - m['TotalDiscountedCost'] == 0\n",
    "finally:\n",
    "    m.add_constraints(con, name='TDC2_TotalDiscountedCost')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total Capacity Constraints\n",
    "\n",
    "```ampl\n",
    "s.t. TCC1_TotalAnnualMaxCapacityConstraint{r in REGION, t in TECHNOLOGY, y in YEAR: TotalAnnualMaxCapacity[r,t,y] <> -1}: \n",
    "    TotalCapacityAnnual[r,t,y] <= TotalAnnualMaxCapacity[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalCapacityAnnual'] <= ds['TotalAnnualMaxCapacity']\n",
    "mask = ds['TotalAnnualMaxCapacity'] >= 0\n",
    "m.add_constraints(con, name='TCC1_TotalAnnualMaxCapacityConstraint', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. TCC2_TotalAnnualMinCapacityConstraint{r in REGION, t in TECHNOLOGY, y in YEAR: TotalAnnualMinCapacity[r,t,y]>0}: \n",
    "    TotalCapacityAnnual[r,t,y] >= TotalAnnualMinCapacity[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalCapacityAnnual'] >= ds['TotalAnnualMinCapacity']\n",
    "mask = ds['TotalAnnualMinCapacity'] > 0\n",
    "m.add_constraints(con, name='TCC2_TotalAnnualMinCapacityConstraint', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Annual Activity Constraints\n",
    "\n",
    "```ampl\n",
    "s.t. AAC1_TotalAnnualTechnologyActivity{r in REGION, t in TECHNOLOGY, y in YEAR}: \n",
    "    sum{l in TIMESLICE} RateOfTotalActivity[r,t,l,y]*YearSplit[l,y] = TotalTechnologyAnnualActivity[r,t,y];\n",
    "```\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (m['RateOfTotalActivity'] * ds['YearSplit']).sum('TIMESLICE') - m['TotalTechnologyAnnualActivity'] == 0\n",
    "m.add_constraints(con, name='AAC1_TotalAnnualTechnologyActivity')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. AAC2_TotalAnnualTechnologyActivityUpperLimit{r in REGION, t in TECHNOLOGY, y in YEAR: TotalTechnologyAnnualActivityUpperLimit[r,t,y] <> -1}:\n",
    "    TotalTechnologyAnnualActivity[r,t,y] <= TotalTechnologyAnnualActivityUpperLimit[r,t,y] ;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalTechnologyAnnualActivity'] <= ds['TotalTechnologyAnnualActivityUpperLimit']\n",
    "mask = ds['TotalTechnologyAnnualActivityUpperLimit'] >= 0\n",
    "m.add_constraints(con, name='AAC2_TotalAnnualTechnologyActivityUpperLimit', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. AAC3_TotalAnnualTechnologyActivityLowerLimit{r in REGION, t in TECHNOLOGY, y in YEAR: TotalTechnologyAnnualActivityLowerLimit[r,t,y]>0}: \n",
    "    TotalTechnologyAnnualActivity[r,t,y] >= TotalTechnologyAnnualActivityLowerLimit[r,t,y] ;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalTechnologyAnnualActivity'] >= ds['TotalTechnologyAnnualActivityLowerLimit']\n",
    "mask = ds['TotalTechnologyAnnualActivityLowerLimit'] > 0\n",
    "m.add_constraints(con, name='AAC3_TotalAnnualTechnologyActivityLowerLimit', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total Activity Constraints\n",
    "\n",
    "```ampl\n",
    "s.t. TAC1_TotalModelHorizonTechnologyActivity{r in REGION, t in TECHNOLOGY}: \n",
    "    sum{y in YEAR} TotalTechnologyAnnualActivity[r,t,y] = TotalTechnologyModelPeriodActivity[r,t];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalTechnologyAnnualActivity'].sum('YEAR') - m['TotalTechnologyModelPeriodActivity'] == 0\n",
    "m.add_constraints(con, name='TAC1_TotalModelHorizonTechnologyActivity')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. TAC2_TotalModelHorizonTechnologyActivityUpperLimit{r in REGION, t in TECHNOLOGY: TotalTechnologyModelPeriodActivityUpperLimit[r,t]<>-1}: \n",
    "    TotalTechnologyModelPeriodActivity[r,t] <= TotalTechnologyModelPeriodActivityUpperLimit[r,t] ;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalTechnologyModelPeriodActivity'] <= ds['TotalTechnologyModelPeriodActivityUpperLimit']\n",
    "mask = ds['TotalTechnologyModelPeriodActivityUpperLimit'] >= 0\n",
    "m.add_constraints(con, name='TAC2_TotalModelHorizonTechnologyActivityUpperLimit', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. TAC3_TotalModelHorizenTechnologyActivityLowerLimit{r in REGION, t in TECHNOLOGY: TotalTechnologyModelPeriodActivityLowerLimit[r,t]>0}: \n",
    "    TotalTechnologyModelPeriodActivity[r,t] >= TotalTechnologyModelPeriodActivityLowerLimit[r,t] ;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['TotalTechnologyModelPeriodActivity'] >= ds['TotalTechnologyModelPeriodActivityLowerLimit']\n",
    "mask = ds['TotalTechnologyModelPeriodActivityLowerLimit'] > 0\n",
    "m.add_constraints(con, name='TAC3_TotalModelHorizenTechnologyActivityLowerLimit', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reserve Margin Constraint\n",
    "\n",
    "```ampl\n",
    "s.t. RM1_ReserveMargin_TechnologiesIncluded_In_Activity_Units{r in REGION, l in TIMESLICE, y in YEAR: ReserveMargin[r,y] > 0}:\n",
    "\tsum {t in TECHNOLOGY}\n",
    "\tTotalCapacityAnnual[r,t,y] * ReserveMarginTagTechnology[r,t,y] * CapacityToActivityUnit[r,t]\n",
    "\t=\n",
    "\tTotalCapacityInReserveMargin[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = ds['ReserveMarginTagTechnology'] * ds['CapacityToActivityUnit'] * m['TotalCapacityAnnual'] - m['TotalCapacityInReserveMargin'] == 0\n",
    "mask = (ds['ReserveMargin'] > 0) & (ds['ReserveMarginTagTechnology'] == 1)\n",
    "m.add_constraints(con, name='RM1_ReserveMargin_TechnologiesIncluded_In_Activity_Units', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```amlp\n",
    "s.t. RM2_ReserveMargin_FuelsIncluded{r in REGION, l in TIMESLICE, y in YEAR: ReserveMargin[r,y] > 0}:\n",
    "\tsum {f in FUEL}\n",
    "\tRateOfProduction[r,l,f,y] * ReserveMarginTagFuel[r,f,y]\n",
    "\t=\n",
    "\tDemandNeedingReserveMargin[r,l,y];\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = ds['ReserveMarginTagFuel'] * m['RateOfProduction'] - m['DemandNeedingReserveMargin'] == 0\n",
    "mask = (ds['ReserveMargin'] > 0) & (ds['ReserveMarginTagFuel'] == 1)\n",
    "m.add_constraints(con, name='RM2_ReserveMargin_FuelsIncluded', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```amlp\n",
    "s.t. RM3_ReserveMargin_Constraint{r in REGION, l in TIMESLICE, y in YEAR: ReserveMargin[r,y] > 0}:\n",
    "\tDemandNeedingReserveMargin[r,l,y] * ReserveMargin[r,y]\n",
    "\t<=\n",
    "\tTotalCapacityInReserveMargin[r,y];\n",
    "\t```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = ds['ReserveMargin'] * m['DemandNeedingReserveMargin'] - m['TotalCapacityInReserveMargin'] == 0\n",
    "mask = ds['ReserveMargin'] > 0\n",
    "m.add_constraints(con, name='RM3_ReserveMargin_Constraint', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RE Production Target NTS: Should change demand for production\n",
    "```ampl\n",
    "s.t. RE1_FuelProductionByTechnologyAnnual{r in REGION, t in TECHNOLOGY, f in FUEL, y in YEAR}: \n",
    "    sum{l in TIMESLICE} ProductionByTechnology[r,l,t,f,y] = ProductionByTechnologyAnnual[r,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['ProductionByTechnology'].sum('TIMESLICE') - m['ProductionByTechnologyAnnual'] == 0\n",
    "mask = ds['OutputActivityRatio'].sum('MODE_OF_OPERATION') != 0\n",
    "m.add_constraints(con, name='RE1_FuelProductionByTechnologyAnnual', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. RE2_TechIncluded{r in REGION, y in YEAR}: sum{t in TECHNOLOGY, f in FUEL} ProductionByTechnologyAnnual[r,t,f,y]*RETagTechnology[r,t,y] = TotalREProductionAnnual[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['ProductionByTechnologyAnnual'] * ds['RETagTechnology'].fillna(0) - m['TotalREProductionAnnual'] == 0\n",
    "mask = ds['RETagTechnology'] == 1\n",
    "m.add_constraints(con, name='RE2_TechIncluded', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. RE3_FuelIncluded{r in REGION, y in YEAR}: sum{l in TIMESLICE, f in FUEL} \n",
    "    RateOfProduction[r,l,f,y]*YearSplit[l,y]*RETagFuel[r,f,y] = RETotalProductionOfTargetFuelAnnual[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RateOfProduction'] * ds['YearSplit'] * ds['RETagFuel'].fillna(0) - m['RETotalProductionOfTargetFuelAnnual'] == 0\n",
    "mask = ds['RETagFuel'] == 1\n",
    "m.add_constraints(con, name='RE3_FuelIncluded', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. RE4_EnergyConstraint{r in REGION, y in YEAR}:\n",
    "    REMinProductionTarget[r,y] * RETotalProductionOfTargetFuelAnnual[r,y] <= TotalREProductionAnnual[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['RETotalProductionOfTargetFuelAnnual'] * ds['REMinProductionTarget'].fillna(0) - m['TotalREProductionAnnual'] <= 0\n",
    "mask = ds['REMinProductionTarget'] > 0\n",
    "m.add_constraints(con, name='RE4_EnergyConstraint', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. RE5_FuelUseByTechnologyAnnual{r in REGION, t in TECHNOLOGY, f in FUEL, y in YEAR}: \n",
    "    sum{l in TIMESLICE} RateOfUseByTechnology[r,l,t,f,y] * YearSplit[l,y] = UseByTechnologyAnnual[r,t,f,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = (ds['YearSplit'] * m['RateOfUseByTechnology']).sum('TIMESLICE') - m['UseByTechnologyAnnual'] == 0\n",
    "mask = ds['InputActivityRatio'].sum('MODE_OF_OPERATION') > 0\n",
    "m.add_constraints(con, name='RE5_FuelUseByTechnologyAnnual', mask=mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Emissions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E1_AnnualEmissionProductionByMode{r in REGION, t in TECHNOLOGY, e in EMISSION, m in MODE_OF_OPERATION, y in YEAR:\n",
    "\t\t\t\t\t\t\t\t\t   EmissionActivityRatio[r,t,e,m,y] <> 0}:\n",
    "\tEmissionActivityRatio[r,t,e,m,y] * TotalAnnualTechnologyActivityByMode[r,t,m,y]\n",
    "\t=\n",
    "\tAnnualTechnologyEmissionByMode[r,t,e,m,y];\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = ds['EmissionActivityRatio'].notnull()\n",
    "con = ds['EmissionActivityRatio'] * m['TotalAnnualTechnologyActivityByMode'] - m['AnnualTechnologyEmissionByMode'] == 0\n",
    "m.add_constraints(con, name='E1_AnnualEmissionProductionByMode', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E2_AnnualEmissionProduction{r in REGION, t in TECHNOLOGY, e in EMISSION, y in YEAR}:\n",
    "\tsum{m in MODE_OF_OPERATION}\n",
    "\tAnnualTechnologyEmissionByMode[r,t,e,m,y]\n",
    "\t=\n",
    "\tAnnualTechnologyEmission[r,t,e,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ds['EMISSION'].size > 0:\n",
    "    mask = ds['EmissionActivityRatio'].sum('MODE_OF_OPERATION') != 0\n",
    "    con = m['AnnualTechnologyEmissionByMode'].sum(dims='MODE_OF_OPERATION') - m['AnnualTechnologyEmission'] == 0\n",
    "    m.add_constraints(con, name='E2_AnnualEmissionProduction', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E3_EmissionsPenaltyByTechAndEmission{r in REGION, t in TECHNOLOGY, e in EMISSION, y in YEAR: EmissionsPenalty[r,e,y] <> 0}:\n",
    "\tAnnualTechnologyEmission[r,t,e,y] * EmissionsPenalty[r,e,y]\n",
    "\t=\n",
    "\tAnnualTechnologyEmissionPenaltyByEmission[r,t,e,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = (ds['EmissionsPenalty'].notnull()) & (ds['EmissionActivityRatio'].sum('MODE_OF_OPERATION') != 0)\n",
    "con = (m['AnnualTechnologyEmission'] * ds['EmissionsPenalty']) - m['AnnualTechnologyEmissionPenaltyByEmission'] == 0\n",
    "m.add_constraints(con, name='E3_EmissionsPenaltyByTechAndEmission', mask=mask)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E4_EmissionsPenaltyByTechnology{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tsum{e in EMISSION} AnnualTechnologyEmissionPenaltyByEmission[r,t,e,y]\n",
    "\t=\n",
    "\tAnnualTechnologyEmissionsPenalty[r,t,y];\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['AnnualTechnologyEmissionPenaltyByEmission'].sum('EMISSION') - m['AnnualTechnologyEmissionsPenalty'] == 0\n",
    "m.add_constraints(con, name='E4_EmissionsPenaltyByTechnology')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E5_DiscountedEmissionsPenaltyByTechnology{r in REGION, t in TECHNOLOGY, y in YEAR}:\n",
    "\tAnnualTechnologyEmissionsPenalty[r,t,y] / DiscountFactorMid[r,y]\n",
    "\t=\n",
    "\tDiscountedTechnologyEmissionsPenalty[r,t,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['AnnualTechnologyEmissionsPenalty'] / discount_factor_mid - m['DiscountedTechnologyEmissionsPenalty'] == 0\n",
    "m.add_constraints(con, name='E5_DiscountedEmissionsPenaltyByTechnology')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E6_EmissionsAccounting1{r in REGION, e in EMISSION, y in YEAR}:\n",
    "\tsum{t in TECHNOLOGY}\n",
    "\tAnnualTechnologyEmission[r,t,e,y]\n",
    "\t=\n",
    "\tAnnualEmissions[r,e,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ds['EMISSION'].size > 0:\n",
    "    con = m['AnnualTechnologyEmission'].sum(dims=['TECHNOLOGY']) - m['AnnualEmissions'] == 0\n",
    "    m.add_constraints(con, name='E6_EmissionsAccounting1')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E7_EmissionsAccounting2{r in REGION, e in EMISSION}:\n",
    "\tsum{y in YEAR} AnnualEmissions[r,e,y]\n",
    "\t=\n",
    "\tModelPeriodEmissions[r,e] - ModelPeriodExogenousEmission[r,e];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if ds['EMISSION'].size > 0:\n",
    "    con = m['ModelPeriodEmissions'] - m['AnnualEmissions'].sum('YEAR') == ds['ModelPeriodExogenousEmission'].fillna(0)\n",
    "    m.add_constraints(con, name='E7_EmissionsAccounting2')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E8_AnnualEmissionsLimit{r in REGION, e in EMISSION, y in YEAR: AnnualEmissionLimit[r, e, y] <> -1}:\n",
    "\tAnnualEmissions[r,e,y] + AnnualExogenousEmission[r,e,y]\n",
    "\t<=\n",
    "\tAnnualEmissionLimit[r,e,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if 'E8_AnnualEmissionsLimit' in m.constraints:\n",
    "    m.remove_constraints('E8_AnnualEmissionsLimit')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "con = m['AnnualEmissions'] <= ds['AnnualEmissionLimit'] - ds['AnnualExogenousEmission'].fillna(0)\n",
    "mask = ds['AnnualEmissionLimit'] != -1\n",
    "m.add_constraints(con, name='E8_AnnualEmissionsLimit', mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```ampl\n",
    "s.t. E9_ModelPeriodEmissionsLimit{r in REGION, e in EMISSION: ModelPeriodEmissionLimit[r, e] <> -1}:\n",
    "\tModelPeriodEmissions[r,e]\n",
    "\t<=\n",
    "\tModelPeriodEmissionLimit[r,e];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# con = m['ModelPeriodEmissions'] <= ds['ModelPeriodEmissionLimit']\n",
    "mask = ds['ModelPeriodEmissionLimit'] != -1\n",
    "m.add_constraints(lhs=m['ModelPeriodEmissions'], \n",
    "                  sign='<=', \n",
    "                  rhs=ds['ModelPeriodEmissionLimit'], \n",
    "                  name='E9_ModelPeriodEmissionsLimit', \n",
    "                  mask=mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Objective Function\n",
    "```ampl\n",
    "minimize cost: sum{r in REGION, y in YEAR} TotalDiscountedCost[r,y];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = m['TotalDiscountedCost'].sum(dims=['REGION', 'YEAR'])\n",
    "m.add_objective(expr=objective, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.to_file(f'{model_name}.lp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(available_solvers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solve(io_api='direct', log_fn='gurobi.log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for variable, data in m.solution.items():\n",
    "    data = data.to_dataframe(name=variable)\n",
    "    data[data[variable] != 0.0].to_csv(f\"results/linopy/{variable}.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# m.variables.get_name_by_label(645700)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# m.constraints.get_name_by_label(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "variable = re.compile('x[0-9]+')\n",
    "constraint = re.compile('c[0-9]+')\n",
    "\n",
    "with open(f'annotated_{model_name}.lp', 'w') as annotated:\n",
    "    with open(f'{model_name}.lp', 'r') as raw_lp_file:\n",
    "        for idx, line in enumerate(raw_lp_file):\n",
    "            vars = variable.search(line)\n",
    "            cons = constraint.match(line)\n",
    "            if cons:\n",
    "                label = int(cons.group(0)[1:])\n",
    "                name = m.constraints.get_name_by_label(label)\n",
    "                line = f\"{name}:\\n\"\n",
    "            if vars:\n",
    "                var = vars.group(0)\n",
    "                label = int(vars.group(0)[1:])\n",
    "                name = m.variables.get_name_by_label(label)\n",
    "                line = line.replace(var, name)\n",
    "            annotated.write(line)\n",
    "            if idx > 1789554:\n",
    "                break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare LP files\n",
    "import re\n",
    "constraints = [x for x in m.constraints]\n",
    "variables = [x for x in m.variables]\n",
    "with open(f'{model_name}_gnu.lp', 'r') as gnu:\n",
    "    with open(f'annotated_{model_name}.lp', 'r') as linopy:\n",
    "        file = gnu.read()\n",
    "        file2 = linopy.read()\n",
    "        for x in constraints:\n",
    "            num_gnu = file.count(x)\n",
    "            num_linopy = file2.count(x)\n",
    "            if num_gnu != num_linopy:\n",
    "                print(f\"{x}: {num_gnu} {num_linopy}\")          \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = m.solution['AnnualEmissions'].to_dataframe()\n",
    "df[df['AnnualEmissions'] != 0].plot(kind='bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solution['UseByTechnologyAnnual'].sum(\"TECHNOLOGY\").to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solution['ProductionByTechnologyAnnual'].sum(\"TECHNOLOGY\").to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solution['Demand'].sum(\"TIMESLICE\").to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solution['TotalDiscountedCost'].to_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = m.solution['ProductionByTechnologyAnnual'].to_dataframe()\n",
    "df[df != 0].dropna().reset_index().pivot_table(values='ProductionByTechnologyAnnual', index=[\"TECHNOLOGY\",\"FUEL\"], columns='YEAR', aggfunc=sum)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m.solution['TotalCapacityAnnual'].sum('REGION').to_dataframe().plot(kind='bar', stacked=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kwargs = dict(broadcast_like=\"vars\", filter_missings=True)\n",
    "vars = m.constraints.iter_ravel(\"vars\", **kwargs)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[x for x in vars]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "5434a038c19fa95fe8846d13fe13863eee5c1bf63c39831b2cc7a52a89e337c7"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}