Workflow for a stoachstic opimization¶

In this example of the ETHOS.FINE framework, a stochastic optimization is performed. It allows a energy system optimization to consider several years of operation input (e.g. different wind power due to different weather conditions) in the optimization of the energy system design. By optimizing a system considering several years of input data for operation, a more robust energy system is achieved.

1. Import packages¶

In [1]:

import fine as fn
import numpy as np
import pandas as pd

import fine as fn import numpy as np import pandas as pd

2. Initialize energy system model¶

Initialize an esM which has two years.

In [2]:

numberOfTimeSteps = 4
hoursPerTimeStep = 2190
numberOfYears = 2  # new test, before =1
investmentPeriodInterval = 1

# Create an energy system model instance and set stochastic model optimization
esM = fn.EnergySystemModel(
    stochasticModel=True,
    locations={"PerfectLand"},
    commodities={"electricity"},
    numberOfTimeSteps=numberOfTimeSteps,
    commodityUnitsDict={"electricity": r"kW$_{el}$"},
    hoursPerTimeStep=hoursPerTimeStep,
    costUnit="1 Euro",
    numberOfInvestmentPeriods=numberOfYears,
    investmentPeriodInterval=investmentPeriodInterval,
    lengthUnit="km",
    verboseLogLevel=2,
)

numberOfTimeSteps = 4 hoursPerTimeStep = 2190 numberOfYears = 2 # new test, before =1 investmentPeriodInterval = 1 # Create an energy system model instance and set stochastic model optimization esM = fn.EnergySystemModel( stochasticModel=True, locations={"PerfectLand"}, commodities={"electricity"}, numberOfTimeSteps=numberOfTimeSteps, commodityUnitsDict={"electricity": r"kW$_{el}$"}, hoursPerTimeStep=hoursPerTimeStep, costUnit="1 Euro", numberOfInvestmentPeriods=numberOfYears, investmentPeriodInterval=investmentPeriodInterval, lengthUnit="km", verboseLogLevel=2, )

3. Add a wind source¶

Either pass the parameters in a dict per investment period (modification per investment period possible) or constant (only one dataframe, value ...)

In [3]:

windOperationRateMax = {}
windOperationRateMax[0] = pd.DataFrame(
    [
        np.array(
            [
                0.4,
                0.4,
                0.6,
                0.6,
            ]
        )
    ],
    index=["PerfectLand"],
).T
windOperationRateMax[1] = pd.DataFrame(
    [
        np.array(
            [
                0.2,
                0.2,
                0.3,
                0.5,
            ]
        )
    ],
    index=["PerfectLand"],
).T
# different opexPerOperation per investmentperiod
windOpexPerOperation = {}
windOpexPerOperation[0] = 0.01
windOpexPerOperation[1] = 0.02

esM.add(
    fn.Source(
        esM=esM,
        name="Wind",
        commodity="electricity",
        hasCapacityVariable=True,
        operationRateMax=windOperationRateMax,
        capacityMax=4e6,
        investPerCapacity=2 * 2190,
        opexPerCapacity=0,
        interestRate=0,
        opexPerOperation=windOpexPerOperation,  # 0.01,
        economicLifetime=1,
    )
)

windOperationRateMax = {} windOperationRateMax[0] = pd.DataFrame( [ np.array( [ 0.4, 0.4, 0.6, 0.6, ] ) ], index=["PerfectLand"], ).T windOperationRateMax[1] = pd.DataFrame( [ np.array( [ 0.2, 0.2, 0.3, 0.5, ] ) ], index=["PerfectLand"], ).T # different opexPerOperation per investmentperiod windOpexPerOperation = {} windOpexPerOperation[0] = 0.01 windOpexPerOperation[1] = 0.02 esM.add( fn.Source( esM=esM, name="Wind", commodity="electricity", hasCapacityVariable=True, operationRateMax=windOperationRateMax, capacityMax=4e6, investPerCapacity=2 * 2190, opexPerCapacity=0, interestRate=0, opexPerOperation=windOpexPerOperation, # 0.01, economicLifetime=1, ) )

3. Add a electricity market¶

In [4]:

costs = pd.DataFrame(
    [
        np.array(
            [
                1,
                1,
                1,
                1,
            ]
        )
    ],
    index=["PerfectLand"],
).T


revenues = {}
revenues[0] = pd.DataFrame(
    [
        np.array(
            [
                0.0,
                0.0,
                0.0,
                0.0,
            ]
        )
    ],
    index=["PerfectLand"],
).T
revenues[1] = pd.DataFrame(
    [
        np.array(
            [
                0.0,
                0.0,
                0.0,
                0.0,
            ]
        )
    ],
    index=["PerfectLand"],
).T

maxpurchase = {}
maxpurchase[0] = pd.DataFrame(
    [
        np.array(
            [
                0.5e3,
                0.5e3,
                4e3,
                4e3,
            ]
        )
    ],
    index=["PerfectLand"],
).T
maxpurchase[1] = pd.DataFrame(
    [
        np.array(
            [
                0.5e3,
                0.5e3,
                4e3,
                4e3,
            ]
        )
    ],
    index=["PerfectLand"],
).T

costs = pd.DataFrame( [ np.array( [ 1, 1, 1, 1, ] ) ], index=["PerfectLand"], ).T revenues = {} revenues[0] = pd.DataFrame( [ np.array( [ 0.0, 0.0, 0.0, 0.0, ] ) ], index=["PerfectLand"], ).T revenues[1] = pd.DataFrame( [ np.array( [ 0.0, 0.0, 0.0, 0.0, ] ) ], index=["PerfectLand"], ).T maxpurchase = {} maxpurchase[0] = pd.DataFrame( [ np.array( [ 0.5e3, 0.5e3, 4e3, 4e3, ] ) ], index=["PerfectLand"], ).T maxpurchase[1] = pd.DataFrame( [ np.array( [ 0.5e3, 0.5e3, 4e3, 4e3, ] ) ], index=["PerfectLand"], ).T

In [5]:

esM.add(
    fn.Source(
        esM=esM,
        name="Electricity market",
        commodity="electricity",
        hasCapacityVariable=False,
        operationRateMax=maxpurchase,
        commodityCostTimeSeries=costs,
        # commodityCost= 1,
        commodityRevenueTimeSeries=revenues,
    )
)  # eur/kWh

esM.add( fn.Source( esM=esM, name="Electricity market", commodity="electricity", hasCapacityVariable=False, operationRateMax=maxpurchase, commodityCostTimeSeries=costs, # commodityCost= 1, commodityRevenueTimeSeries=revenues, ) ) # eur/kWh

5. Add a industry sink¶

In [6]:

revenuesDemand = {}
revenuesDemand[0] = pd.DataFrame(
    [
        np.array(
            [
                0.1,
                0.1,
                0.1,
                0.1,
            ]
        )
    ],
    index=["PerfectLand"],
).T
revenuesDemand[1] = pd.DataFrame(
    [
        np.array(
            [
                0.2,
                0.2,
                0.2,
                0.2,
            ]
        )
    ],
    index=["PerfectLand"],
).T

demand = {}
demand[0] = pd.DataFrame(
    [
        np.array(
            [
                2e3,
                1e3,
                1e3,
                1e3,
            ]
        )
    ],
    index=["PerfectLand"],
).T  # first investmentperiod
demand[1] = pd.DataFrame(
    [
        np.array(
            [
                2e3,
                1e3,
                1e3,
                1e3,
            ]
        )
    ],
    index=["PerfectLand"],
).T  # second investmentperiod

esM.add(
    fn.Sink(
        esM=esM,
        name="EDemand",
        commodity="electricity",
        hasCapacityVariable=False,
        operationRateFix=demand,
        commodityRevenueTimeSeries=revenuesDemand,  # new compared to original model
    )
)

revenuesDemand = {} revenuesDemand[0] = pd.DataFrame( [ np.array( [ 0.1, 0.1, 0.1, 0.1, ] ) ], index=["PerfectLand"], ).T revenuesDemand[1] = pd.DataFrame( [ np.array( [ 0.2, 0.2, 0.2, 0.2, ] ) ], index=["PerfectLand"], ).T demand = {} demand[0] = pd.DataFrame( [ np.array( [ 2e3, 1e3, 1e3, 1e3, ] ) ], index=["PerfectLand"], ).T # first investmentperiod demand[1] = pd.DataFrame( [ np.array( [ 2e3, 1e3, 1e3, 1e3, ] ) ], index=["PerfectLand"], ).T # second investmentperiod esM.add( fn.Sink( esM=esM, name="EDemand", commodity="electricity", hasCapacityVariable=False, operationRateFix=demand, commodityRevenueTimeSeries=revenuesDemand, # new compared to original model ) )

6. Optimize¶

In [7]:

esM.optimize(
    timeSeriesAggregation=False,
    solver=fn.utils.ImplementedSolvers.STANDARD_SOLVER.value,
)

esM.optimize( timeSeriesAggregation=False, solver=fn.utils.ImplementedSolvers.STANDARD_SOLVER.value, )

GLPSOL: GLPK LP/MIP Solver, v4.65
Parameter(s) specified in the command line:
 --write C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmpjot3ih00.glpk.raw --wglp
 C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmp0rg09re0.glpk.glp --cpxlp C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmp21dx5j_a.pyomo.lp
Reading problem data from 'C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmp21dx5j_a.pyomo.lp'...
24 rows, 33 columns, 56 non-zeros
194 lines were read
Writing problem data to 'C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmp0rg09re0.glpk.glp'...
169 lines were written
GLPK Simplex Optimizer, v4.65
24 rows, 33 columns, 56 non-zeros
Preprocessing...
9 rows, 10 columns, 18 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  1.314e+03  ratio =  1.314e+03
GM: min|aij| =  9.834e-01  max|aij| =  1.017e+00  ratio =  1.034e+00
EQ: min|aij| =  9.691e-01  max|aij| =  1.000e+00  ratio =  1.032e+00
Constructing initial basis...
Size of triangular part is 9
      0: obj =   4.560000000e+03 inf =   7.075e+01 (4)
      4: obj =   3.456000000e+04 inf =   0.000e+00 (0)
*    11: obj =   2.914000000e+04 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (50671 bytes)
Writing basic solution to 'C:\Users\RAEF3F~1.MAI\AppData\Local\Temp\tmpjot3ih00.glpk.raw'...
66 lines were written

7. Get results of robust energy system, optimized for two years of operations¶

In [8]:

esM.getOptimizationSummary("SourceSinkModel").loc["Wind", "capacity", "[kW$_{el}$]"][
    "PerfectLand"
]

esM.getOptimizationSummary("SourceSinkModel").loc["Wind", "capacity", "[kW$_{el}$]"][ "PerfectLand" ]

Out[8]:

3.42465753424658