Integrated Energy Systems: Nuclear Case Study

In this notebook, we formulate a multi-period optimization problem that determines the optimal design of an integrated energy system producing electricity and hydrogen such that its net present value (NPV) is maximized for a given market signal. We consider an existing nuclear power plant of 1000 MW capcity as the genator. In general, due to safety restrictions, nuclear power plants are not amenable for ramping. Therefore, we assume that the power plant always operates at its base load (i.e., 1000 MW). We are interested in determining whether the following options are attractive (i.e., maximize NPV):

  • Produce hydrogen using a portion of the electricity, especially during periods of low electricity demand, and sell it.

  • Produce hydrogen, store it in a tank, and combust it to produce electrity during the periods of high electricity demand.

Towards this, we formulate the superstructure shown in the figure below.

image0

In the figure, arrows in red and blue denote that electricity and material, respectively, flows through them. ID denotes the name of the object containing the corresponding unit’s model in the code, and DD and OD stand for design decision and operating decision, respectively.

Here, we use a PEM (polymer electrolyte membrane) electrolyzer to produce hydrogen via water electrolysis. Current PEM technology requires ~54 kWh/kg of hydrogen (see https://www.h-tec.com/en/products/detail/h-tec-pem-electrolyser-me100-350/me100-350/ , last accessed on February 18, 2022). The produced hydrogen is stored in a tank, which can either be sold to the market via a pipeline, or combusted in a hydrogen turbine to produce electricity during the periods of high demand for electricity. Note that

  • The hydrogen tank is modeled as a simple inventory model i.e., only mass balance is enforced. The enthalpy balance is not included in the model. If needed, the SimplifiedHydrogenTank model in the flowsheet can be replaced with a more detailed tank model (the DISPATCHES repository has a more detailed tank model that enforces enthalpy balance as well).

  • The hydrogen turbine is modeled as a compressor, a stochiometric reactor, and a turbine connected in series.

In the flowsheet model, we use two different thermodynamic packages: one for the PEM electrolyzer and the hydrogen tank, and the other for the mixer and the hydrogen turbine. This is because, the material stream(s) to and from the former units only contain hydrogen. Whereas, the material stream(s) to and from the latter units contain nitrogen, oxygen, water and argon, in addition to hydrogen. The translator block facilitates the use of multiple thermodynamic packages by connecting properties across the packages.

Objective

For a given market signal, our objective is to determine the optimal design decisions,

  • Size of the PEM electrolyzer (pem_capacity: maximum rated capacity of the PEM electrolyzer, in MW)

  • Size of the hydrogen tank (tank_capacity: maximum amount of hydrogen that can be stored in the tank, in kg)

  • Size of the hydrogen turbine (turbine_capacity: maximum power the hydrogen turbine can produce, in MW)

and the optimal operating decisions,

  • Split fraction of electricity to the pem electrolyzer (m.fs.np_power_split.split_fraction["np_to_pem", 0])

  • Molar flowrate of hydrogen to the pipeline (m.fs.h2_tank.outlet_to_pipeline.flow_mol[0])

  • Molar flowrate of hydrogen to the turbine (m.fs.h2_tank.outlet_to_turbine.flow_mol[0])

which maximize the NPV.

First, we import the required packages and functions needed for the formulation of the optimization problem.

[1]:
# General python imports
import json
import logging

# Pyomo imports
from pyomo.environ import (
    ConcreteModel,
    RangeSet,
    Var,
    NonNegativeReals,
    Constraint,
    Expression,
    Objective,
    value,
    maximize,
    units as pyunits,
)

# Register currency as a unit
pyunits.load_definitions_from_strings(["USD = [currency]"])

# IDAES imports
from idaes.core.solvers import get_solver
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.apps.grid_integration import MultiPeriodModel

# Nuclear flowsheet function imports
from dispatches.case_studies.nuclear_case.nuclear_flowsheet import (
    build_ne_flowsheet,
    fix_dof_and_initialize,
)

Simulation of the Flowsheet

Before we formulate the multiperiod optimization problem, we first simulate the nuclear flowsheet and print some results. This helps the reader to familiarize themselves with the names, along with their default units, of a few important variables. The readers are encouraged to go through the nuclear_flowsheet.py script. It contains the function build_ne_flowsheet which assembles models for all the units in the flowsheet, and connects them via Arc objects. It also contains the function fix_dof_and_initialize which fixes the degrees of freedom and initializes the entire flowsheet.

The flowsheet has four degrees of freedom viz.,

  • Split fraction of electricity to the grid in the power splitter (variable m.fs.np_power_split.split_fraction["np_to_grid", 0] in the model). We fix this variable to 0.8.

  • Molar flowrate of hydrogen to the pipeline (variable m.fs.h2_tank.outlet_to_pipeline.flow_mol[0]). We fix this variable to 10 mol/s.

  • Molar flowrate of hydrogen to the turbine (variable m.fs.h2_tank.outlet_to_turbine.flow_mol[0]). We fix this variable to 10 mol/s.

  • Initial holdup of hydrogen in the tank (variable m.fs.h2_tank.tank_holdup_previous[0]). We fix this variable to 0 mol.

These variables are fixed in the fix_dof_and_initialize function. In addition to the above three variables, we also fix the variables shown in the table below. We do not refer these variables as degrees of freedom, because they remain fixed at these values in the multiperiod optimization model.

Name

Variable

Value

Pressure difference across h2_turbine’s compressor

m.fs.h2_turbine.compressor.deltaP

24.01 bar

Isentropic efficiency of h2_turbine’s compressor

m.fs.h2_turbine.compressor.efficiency_isentropic

0.86

Conversion of hydrogen in h2_turbine’s reactor

m.fs.h2_turbine.stoic_reactor.conversion

0.99

Pressure deifference across h2_turbine’s turbine

m.fs.h2_turbine.turbine.deltaP

-24.01 bar

Isentripic efficiency of h2_turbine’s turbine

m.fs.h2_turbine.turbine.efficiency_isentropic

0.89

Molar flow rate of air to h2_turbine

m.fs.mixer.air_feed.flow_mol[0]

10.76 * molar flowrate of hydrogen to turbine

Temperature of air

m.fs.mixer.air_feed.temperature[0]

300 K

Pressure of air

m.fs.mixer.air_feed_pressure[0]

1.01325 bar

Duration of the simulation for h2_tank

m.fs.h2_tank.dt

3600 s

[2]:
# Build the nuclear flowsheet with the capacity of the power plant as 1000 MW
m = build_ne_flowsheet(np_capacity=1000)

# Fix the degrees of freedom and initialize
fix_dof_and_initialize(
    m,
    split_frac_grid=0.8,
    tank_holdup_previous=0,
    flow_mol_to_pipeline=10,
    flow_mol_to_turbine=10,
)

# Ensure that the resulting model is a square problem i.e., its degrees of freedom must be 0
print("Degrees of freedom: ", degrees_of_freedom(m))
assert degrees_of_freedom(m) == 0

# Create a solver object with the default solver (IPOPT)
solver = get_solver()

# Simulate the entire flowsheet
solver.solve(m, tee=True)
2022-10-08 11:01:50 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
Degrees of freedom:  0
Ipopt 3.14.7: nlp_scaling_method=gradient-based
tol=1e-06


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.7, running with linear solver MUMPS 5.5.0.

Number of nonzeros in equality constraint Jacobian...:      593
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      204

Total number of variables............................:      215
                     variables with only lower bounds:       20
                variables with lower and upper bounds:      180
                     variables with only upper bounds:        0
Total number of equality constraints.................:      215
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 5.55e+04 1.00e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  0.0000000e+00 3.49e-08 2.36e-01  -1.0 9.80e-03    -  9.90e-01 1.00e+00h  1

Number of Iterations....: 1

                                   (scaled)                 (unscaled)
Objective...............:   0.0000000000000000e+00    0.0000000000000000e+00
Dual infeasibility......:   0.0000000000000000e+00    0.0000000000000000e+00
Constraint violation....:   2.3283064365386963e-10    3.4924596548080444e-08
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.3283064365386963e-10    3.4924596548080444e-08


Number of objective function evaluations             = 2
Number of objective gradient evaluations             = 2
Number of equality constraint evaluations            = 2
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 1
Total CPU secs in IPOPT (w/o function evaluations)   =      0.000
Total CPU secs in NLP function evaluations           =      0.000

EXIT: Optimal Solution Found.
[2]:
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 215, 'Number of variables': 215, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 0.18897509574890137}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

Verify that IPOPT converges to the optimal solution. Next, we print some results. Note the names of the variables and their default units (we use Pyomo’s units, imported as pyunits, to obtain the units of a variable). Units are needed to correctly define the cash flow expressions later.

[3]:
def get_units(obj):
    return str(pyunits.get_units(obj))

# Print results: power splitter
print("Nuclear plant power production : ", m.fs.np_power_split.electricity[0].value,
      get_units(m.fs.np_power_split.electricity[0]))
print("Electricity to grid            : ", m.fs.np_power_split.np_to_grid_port.electricity[0].value,
      get_units(m.fs.np_power_split.np_to_grid_port.electricity[0]))
print("Electricity to PEM             : ", m.fs.np_power_split.np_to_pem_port.electricity[0].value,
      get_units(m.fs.np_power_split.np_to_pem_port.electricity[0]))
print()

# Print results: PEM electrolyzer
print("Flowrate of H2 from pem        : ", m.fs.pem.outlet.flow_mol[0].value,
      get_units(m.fs.pem.outlet.flow_mol[0]))
print()

# Print results: Hydrogen tank
print("Flowrate of H2 to tank         : ", m.fs.h2_tank.inlet.flow_mol[0].value,
      get_units(m.fs.h2_tank.inlet.flow_mol[0]))
print("Flowrate of H2 to pipeline     : ", m.fs.h2_tank.outlet_to_pipeline.flow_mol[0].value,
      get_units(m.fs.h2_tank.outlet_to_pipeline.flow_mol[0]))
print("Flowrate of H2 to turbine      : ", m.fs.h2_tank.outlet_to_turbine.flow_mol[0].value,
      get_units(m.fs.h2_tank.outlet_to_turbine.flow_mol[0]))
print("Initial tank holdup            : ", m.fs.h2_tank.tank_holdup_previous[0].value,
      get_units(m.fs.h2_tank.tank_holdup_previous[0]))
print("Tank holdup at the end of 1 hr : ", m.fs.h2_tank.tank_holdup[0].value,
      get_units(m.fs.h2_tank.tank_holdup[0]))

# Print results: Hydrogen Turbine
print("H2 Turbine's compressor work   : ", m.fs.h2_turbine.compressor.work_mechanical[0].value,
      get_units(m.fs.h2_turbine.compressor.work_mechanical[0]))
print("H2 Turbine's turbine work      : ", m.fs.h2_turbine.turbine.work_mechanical[0].value,
      get_units(m.fs.h2_turbine.turbine.work_mechanical[0]))
print("Net power produced by turbine  : ", -value(m.fs.h2_turbine.work_mechanical[0]),
      get_units(m.fs.h2_turbine.work_mechanical[0]))
Nuclear plant power production :  1000000.0 kW
Electricity to grid            :  800000.0 kW
Electricity to PEM             :  199999.99999999997 kW

Flowrate of H2 from pem        :  505.4811999999999 mol/s

Flowrate of H2 to tank         :  505.4811999999999 mol/s
Flowrate of H2 to pipeline     :  10 mol/s
Flowrate of H2 to turbine      :  10 mol/s
Initial tank holdup            :  0 mol
Tank holdup at the end of 1 hr :  1747732.3199999996 mol
H2 Turbine's compressor work   :  1764604.5955977994 kg*m**2/s**3
H2 Turbine's turbine work      :  -2672169.126200441 kg*m**2/s**3
Net power produced by turbine  :  907564.5306026414 kg*m**2/s**3

Observe that the power variables in the power splitter model and the PEM electrolyzer model use kW. Whereas, the power variables (mechanical work) use in the hydrogen turbine model use W. This complete the simulation of the flowsheet. Before proceeding further, we delete the object m to avoid confusion

[4]:
# Delete the object containing the flowsheet
del m

Multiperiod Optimization Model: Deterministic

LMP Signal

Our objective is to determine the optimal size of the PEM, tank and turbine maximizing the NPV for a given ‘market signal’. Throughout the notebook, by market signal or price signal, we refer to the locational marginal price (LMP) (selling price of electricity, in $/MWh) as a function of time. The LMP depends on several factors such as weather, demand, generator mix of the grid, and so forth. Due to the uncertainty associated with some of these factors, it is not possible to predict the exact value of the LMP way into the future. Nevertheless, in this section, we assume that the LMP signal is accurate (i.e., there is no uncertainty in the price). Later, we show how the same framework/workflow can be used to easily formulate a stochastic program to handle the uncertainty in the price signal.

Here, we use the LMP data contained in the file lmp_signal.json. This dataset is generated by RAVEN/FORCE (https://github.com/idaholab/raven) using the 2019/2020 New York Independent System Operator (NYISO) price data. FORCE divides the 365 days of a year into a specified number of clusters (20, in our case), and generates the LMP signal for each cluster (Note that the LMP signal is the same for all the days of a cluster). The figure below plots the LMP signal for cluster 1 (left) and cluster 8 (right).

image0

As evident from the figure, the price can vary significantly in a day, and from cluster to cluster. Given the variation, we are interested in determining if producing hydrogen, especially during the periods when LMP is low, is attractive or not.

For the demonstration, we do not use the entire dataset in lmp_signal.json. Instead, we use the LMP signal only for the years 2022 and 2032. We assume that the plant lifetime is twenty years. We use the 2022 LMP signal for the first ten years (i.e., 2022 - 2031), and the 2032 LMP signal for the next ten years (i.e., 2032 - 2041). If desired, the LMP data for the intermediate years can be easily included in the model with a slight modification to the code.

[5]:
# Load the LMP dataset
with open("lmp_signal.json") as fp:
    lmp_dataset = json.load(fp)

# Gather the LMP data needed for the deterministic case
# Notation: lmp_dataset[scenario][year][cluster/day][time/hour]
lmp_deterministic = {year: {cluster: {hour: lmp_dataset["0"][str(year)][str(cluster)][str(hour)]
                                      for hour in range(1, 25)}
                            for cluster in range(1, 21)}
                     for year in [2022, 2032]}

# Size of each cluster/number of days in a year represented by the cluster
weights_days = {year: {cluster: lmp_dataset[str(0)][str(year)][str(cluster)]["num_days"]
                       for cluster in range(1, 21)}
                for year in [2022, 2032]}

Now, we formulate the multi-period price-taker problem to determine the optimal design and operating decisions maximizing the NPV. The optimization problem is of the form

\[\begin{split} \begin{aligned} \max_{D, u_{t, d}} \quad & \text{NPV}(D, u_{t, d})\\ & g(u_{t, d}) = 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C} \\ & h(u_{t, d}) \le 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C} \\ & f(u_{t-1,d}, u_{t,d}) = 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C} \\ & u_{t,d} \le D, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C} \end{aligned}\end{split}\]

Here, the sets \(\mathcal{T} = \{1, \dots, 24\}\) and \(\mathcal{C} = \{1, \dots, 20\}\) denote the set of hours in a day and the set of clusters/days in a year, respectively. \(D\) denotes the design decisions (or, first-stage variables) viz. size of the PEM electrolyzer (pem_capacity), tank (tank_capacity) and the turbine (turbine_capacity). \(u_{t, d}\) denotes the operating decisions (or, second-stage decisions) at time \(t\) of day/cluster \(d\). \(g(u_{t, d}) = 0\) and \(h(u_{t,d}) \le 0\) denote the flowsheet model, \(f(u_{t-1,d}, u_{t,d}) = 0\) connects the operating decisions at time \(t-1\) and at time \(t\) (e.g., tank holdup), and \(u_{t,d} \le D\) ensures that the operating decision values never exceed the design capacity values (e.g., tank holdup at any time instant must not exceed the tank capacity).

We use the MultiPeriodModel class in IDAES to formulate the multi-period optimization problem and specify the following arguments - n_time_points: Number of elements in \(\mathcal{T}\) - set_days: The set of clusters/representative days - set_years: The set of years - process_model_func: Function that builds an instance of the process flowsheet - initialization_func: Function that initializes an instance of the process flowsheet - unfix_dof_func: Function that unfixes a degrees of freedom for optimization - linking_variable_func: Function that yields pairs of variables that are connected across two consecutive time periods - flowsheet_options: Dictionary containing arguments required for process_model_func - initialization_options: Dictionary containing arguments required for initialization_func - use_stochastic_build: If True, uses the generalized method, build_stochastic_multi_period, for formulating the multi-period problem. - outlvl: Output level

For our problem, n_time_points=24,set_days=[1, 2, ..., 20], set_years=[2022, 2032], process_model_func=build_ne_flowsheet, initialization_func=fix_dof_and_initialize. We now construct two functions, unfix_dof and get_linking_variable_pairs, for unfix_dof and linking_variable_func arguments.

unfix_dof unfixes a few degrees of freedom for optimization. In particular, it - Unfixes the split fractions of the power splitter. The optimizer then determines the optimal split fraction of the electricity to the grid, and to the PEM electrolyzer. - Unfixes the initial tank holdup. The initial tank holdup is governed by the final tank holdup at the previous hour. - Unfixes the molar flow rate of hydrogen to the turbine and to the pipeline. - Unfixes the molar flow rate of air to the turbine. - Adds a constraint to fix the ratio of molar flow rate of air to molar flowrate of hydrogen (fuel-air ratio). - Specifies a small non-zero bounds on a few flow variables to avoid convergence issues.

[6]:
def unfix_dof(m, air_h2_ratio=10.76):
    """
    This function unfixes a few degrees of freedom for optimization
    """
    # Unfix the electricity split in the electrical splitter
    m.fs.np_power_split.split_fraction["np_to_grid", 0].unfix()

    # Unfix the holdup_previous and outflow variables
    m.fs.h2_tank.tank_holdup_previous.unfix()
    m.fs.h2_tank.outlet_to_turbine.flow_mol.unfix()
    m.fs.h2_tank.outlet_to_pipeline.flow_mol.unfix()

    # Unfix the flowrate of air to the mixer
    m.fs.mixer.air_feed.flow_mol.unfix()

    # Add a constraint to maintain the air to hydrogen flow ratio
    m.fs.mixer.air_h2_ratio = Constraint(
        expr=m.fs.mixer.air_feed.flow_mol[0] ==
             air_h2_ratio * m.fs.mixer.hydrogen_feed.flow_mol[0])

    # Set bounds on variables. A small non-zero value is set as the lower
    # bound on molar flowrates to avoid convergence issues
    m.fs.pem.outlet.flow_mol[0].setlb(0.001)

    m.fs.h2_tank.inlet.flow_mol[0].setlb(0.001)
    m.fs.h2_tank.outlet_to_turbine.flow_mol[0].setlb(0.001)
    m.fs.h2_tank.outlet_to_pipeline.flow_mol[0].setlb(0.001)

    m.fs.translator.inlet.flow_mol[0].setlb(0.001)
    m.fs.translator.outlet.flow_mol[0].setlb(0.001)

    m.fs.mixer.hydrogen_feed.flow_mol[0].setlb(0.001)
[7]:
def get_linking_variable_pairs(ct, ft):
    """Yield pairs of variables that need to be connected across time periods.

    The only variable that is connected across two time periods is the tank holdup.

    Args:
        ct: current time step
        ft: the next time step
    """
    return [(ct.fs.h2_tank.tank_holdup[0], ft.fs.h2_tank.tank_holdup_previous[0])]
[8]:
def add_capacity_variables(m):

    """
    This function declares the first-stage variables or design decisions, and
    adds constraints that ensure that the operational variables never exceed their
    design values.
    """
    if hasattr(m, "set_period"):
        set_period = m.set_period
    else:
        set_period = m.parent_block().set_period

    # Declare first-stage variables (Design decisions)
    m.pem_capacity = Var(
        within=NonNegativeReals,
        doc="Maximum capacity of the PEM electrolyzer (in kW)",
        units=pyunits.kW,
    )
    m.tank_capacity = Var(
        within=NonNegativeReals,
        doc="Maximum holdup of the tank (in mol)",
        units=pyunits.mol
    )
    m.h2_turbine_capacity = Var(
        within=NonNegativeReals,
        doc="Maximum power output from the turbine (in W)",
        units=pyunits.W
    )

    m.pem_capacity_constraint = Constraint(set_period)
    m.tank_capacity_constraint = Constraint(set_period)
    m.turbine_capacity_constraint = Constraint(set_period)

    for t in set_period:
        # Ensure that the electricity to the PEM elctrolyzer does not exceed the PEM capacity
        m.pem_capacity_constraint.add(
            t, m.period[t].fs.pem.electricity[0] <= m.pem_capacity
        )
        # Ensure that the final tank holdup does not exceed the tank capacity
        m.tank_capacity_constraint.add(
            t, m.period[t].fs.h2_tank.tank_holdup[0] <= m.tank_capacity
        )
        # Ensure that the power generated by the turbine does not exceed the turbine capacity
        m.turbine_capacity_constraint.add(
            t, - m.period[t].fs.h2_turbine.work_mechanical[0] <= m.h2_turbine_capacity
        )

Next, we use the MultiPeriodModel class to build the constraints \(g(u_{t,d}) = 0\) and \(h(u_{t,d}) \le 0\) \(\forall \; t \in \mathcal{T}\) and \(\forall \; d \in \mathcal{C}\). It creates an instance of the flowsheet (using the build_ne_flowsheet function) for each time instance, fixes its degrees of freedom and initializes it (using the fix_dof_and_initialize function), unfixes operational degrees of freedom (using the unfix_dof function), and adds constraints of the form \(f(u_{t-1,d}, u_{t,d}) = 0\) (using the get_linking_variable_pairs functions). Then, we add the capacity constraints of the form \(u_{t,d} \le D\) using add_capacity_variables function.

Tip: Building the multiperiod model can take some time. To monitor the progress, set outlvl=logging.INFO

[9]:
# Formulate the multiperiod optimization problem
m = MultiPeriodModel(
    n_time_points=24,
    set_days=[i for i in range(1, 21)],
    set_years=[2022, 2032],
    process_model_func=build_ne_flowsheet,
    initialization_func=fix_dof_and_initialize,
    unfix_dof_func=unfix_dof,
    linking_variable_func=get_linking_variable_pairs,
    flowsheet_options={"np_capacity": 1000},
    initialization_options={
        "split_frac_grid": 0.8,
        "tank_holdup_previous": 0,
        "flow_mol_to_pipeline": 10,
        "flow_mol_to_turbine": 10,
    },
    use_stochastic_build=True,
    outlvl=logging.WARNING,
)

# Add first-stage variables
add_capacity_variables(m)
[+   0.00] Beginning the formulation of the multiperiod optimization problem.
[+ 168.69] Completed the formulation of the multiperiod optimization problem.
2022-10-08 11:05:23 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
[+   7.48] Created an instance of the flowsheet and initialized it.
[+  17.89] Initialized the entire multiperiod optimization model.
[+   0.42] Unfixed the degrees of freedom from each period model.

Next, we use the numbers in the table below for constructing the cash flow expressions.

Unit

CAPEX

Fixed O&M

Variable O&M

PEM Electrolyzer

$1630 /kW

$47.9 /kW

$1.3 /MWh

Hydrogen tank

$29 /kWh

0

0

Hydrogen turbine

$947 /kW

$7 /kW

$4.25 /Mwh

We refer the reader to the Simulation of the Flowsheet section for the default units of the decision variables. We convert the cost numbers so that the units are consistent with that of the variable.

[10]:
CAPEX = {
    # "component": value * units
    "pem": 1630 * (pyunits.USD / pyunits.kW),
    "tank": 29 * (pyunits.USD / pyunits.kWh),
    "turbine": 947 * (pyunits.USD / pyunits.kW),
}

FOM = {
    # "component": value * units
    "pem": 47.9 * (pyunits.USD / pyunits.kW),
    "turbine": 7 * (pyunits.USD / pyunits.kW),
}

VOM = {
    # "component": value * units
    "pem": 1.3 * (pyunits.USD / pyunits.MWh),
    "turbine": 4.25 * (pyunits.USD / pyunits.MWh)
}
LHV_H2 = 33.3 * pyunits.kWh / pyunits.kg
MW_H2 = 2.016e-3 * pyunits.kg / pyunits.mol

Next, we define

  • A function append_costs_and_revenue to calculate the total variable operating cost and the total revenue

  • A function append_cashflows to build expressions for all the cash flows (CAPEX, fixed O&M, depreciation, etc.)

Observe that we assume (Line 39 in the function below) that each day begins with an empty tank (i.e., zero holdup). This is one of the limitations of working with representative days/clusters. The excess storage at the end of each day cannot be transferred to the following day. This limitation can be overcome by formulating an optimization problem for an entire year, instead of a few representative days. This, however, increases the problem size substantially.

[11]:
def append_op_costs_and_revenue(m, lmp, h2_price):

    lmp = lmp * pyunits.USD / pyunits.MWh
    h2_price = h2_price * pyunits.USD / pyunits.kg

    m.fs.vom = Expression(
        expr=pyunits.convert(VOM["pem"] * m.fs.pem.electricity[0] * pyunits.hr, to_units=pyunits.USD) -
        pyunits.convert(VOM["turbine"] * m.fs.h2_turbine.work_mechanical[0] * pyunits.hr, to_units=pyunits.USD)
    )

    m.fs.electricity_revenue = Expression(
        expr= lmp * pyunits.convert(m.fs.np_power_split.np_to_grid_port.electricity[0] * pyunits.hr, to_units=pyunits.MWh) -
        lmp * pyunits.convert(m.fs.h2_turbine.work_mechanical[0] * pyunits.hr, to_units=pyunits.MWh)
    )

    m.fs.hydrogen_revenue = Expression(
        expr=pyunits.convert(h2_price * MW_H2 * m.fs.h2_tank.outlet_to_pipeline.flow_mol[0] * pyunits.hr,
        to_units=pyunits.USD)
    )
[12]:
def append_cashflows(m):

    if hasattr(m, "set_time"):
        ps = m                         # Object containing information on sets and parameters
    else:
        ps = m.parent_block()

    set_time = ps.set_time             # Set of hours
    set_years = ps.set_years           # Set of years
    plant_life = ps.plant_life         # Plant lifetime
    tax_rate = ps.tax_rate             # Corporate tax rate
    discount_rate = ps.discount_rate   # Discount rate

    representative_days = False
    if hasattr(ps, "set_days"):
        representative_days = True
        set_days = ps.set_days         # Set of days/clusters
        weights_days = ps.weights_days # Weights associated with each cluster

    years_vec = [y - set_years.at(1) + 1 for y in set_years]
    years_vec.append(plant_life + 1)
    weights_years = {y: sum(1 / (1 + discount_rate) ** i
                            for i in range(years_vec[j], years_vec[j + 1]))
                     for j, y in enumerate(set_years)}

    # PEM CAPEX: $1630/kWh and pem_capacity is in kW,
    # Tank CAPEX: $29/kWh, the LHV of hydrogen is 33.3 kWh/kg,
    # the molecular mass of hydrogen is 2.016e-3 kg/mol and
    # tank_capacity is in moles
    # Turbine CAPEX: $947/kWh and turbine_capacity is in W
    m.capex = Expression(
        expr=(pyunits.convert(CAPEX["pem"] * m.pem_capacity, to_units=pyunits.USD) +
              pyunits.convert(CAPEX["tank"] * LHV_H2 * MW_H2 * m.tank_capacity, to_units=pyunits.USD) +
              pyunits.convert(CAPEX["turbine"] * m.h2_turbine_capacity, to_units=pyunits.USD)),
        doc="Total capital cost (in USD)"
    )

    # Fixed O&M of PEM: $47.9/kW
    # Fixed O&M of turbine: $7/kW
    @m.Expression(set_years,
                  doc="Fixed O&M cost per year (in USD)")
    def fixed_om_cost(blk, y):
        return (
            pyunits.convert(FOM["pem"] * m.pem_capacity, to_units=pyunits.USD) +
            pyunits.convert(FOM["turbine"] * m.h2_turbine_capacity, to_units=pyunits.USD)
        )

    # Variable O&M: PEM: $1.3/MWh and turbine: $4.25/MWh
    @m.Expression(set_years,
                  doc="Total variable O&M cost per year (in USD)")
    def variable_om_cost(blk, y):
        if representative_days:
            return (
                sum(weights_days[y][d] * blk.period[t, d, y].fs.vom
                    for t in set_time for d in set_days)
            )

        else:
            return sum(blk.period[t, y].fs.vom for t in set_time)

    @m.Expression(set_years,
                  doc="Revenue generated by selling electricity per year (in USD)")
    def electricity_revenue(blk, y):
        if representative_days:
            return (
                sum(weights_days[y][d] * blk.period[t, d, y].fs.electricity_revenue
                    for t in set_time for d in set_days)
            )

        else:
            return sum(blk.period[t, y].fs.electricity_revenue for t in set_time)


    @m.Expression(set_years,
                  doc="Revenue generated by selling hydrogen per year (in USD)")
    def h2_revenue(blk, y):
        if representative_days:
            return (
                sum(weights_days[y][d] * blk.period[t, d, y].fs.hydrogen_revenue
                    for t in set_time for d in set_days)
            )

        else:
            return sum(blk.period[t, y].fs.hydrogen_revenue for t in set_time)

    @m.Expression(set_years,
                  doc="Depreciation value per year (in USD)")
    def depreciation(blk, y):
        return (
            blk.capex / plant_life
        )

    @m.Expression(set_years,
                  doc="Net profit per year (in USD)")
    def net_profit(blk, y):
        return (
            blk.depreciation[y] + (1 - tax_rate) * (+ blk.h2_revenue[y]
                                                    + blk.electricity_revenue[y]
                                                    - blk.fixed_om_cost[y]
                                                    - blk.variable_om_cost[y]
                                                    - blk.depreciation[y])
        )

    m.npv = Expression(
        expr=sum(weights_years[y] * m.net_profit[y] for y in set_years) - m.capex,
        doc="Net present value (in USD)"
    )

Note: It is not necessary to define the above functions for the deterministic case. Nevertheless, we do so because they will be useful for the stochastic case, as we shall see in the next section.

We now build the remaining constraints and the objective function and solve the optimization problem.

[13]:
# Define parameters
m.plant_life = 20                            # Plant lifetime: 20 years
m.tax_rate = 0.2                             # Corporate tax rate: 20%
m.discount_rate = 0.08                       # Discount rate: 8%
m.h2_price = 3                               # Selling price of hydrogen: $3/kg
m.h2_demand = 1                              # Maximum amount of hydrogen that can be sold: 1 kg/s
m.weights_days = weights_days                # number of days represented by each cluster
[14]:
# Set initial holdup for each day (Assumed to be zero at the beginning of each day)
for y in m.set_years:
    for d in m.set_days:
        m.period[1, d, y].fs.h2_tank.tank_holdup_previous.fix(0)

# Hydrogen demand constraint. divide the RHS by the molecular mass to convert kg/s to mol/s
@m.Constraint(m.set_time, m.set_days, m.set_years)
def hydrogen_demand_constraint(blk, t, d, y):
    return blk.period[t, d, y].fs.h2_tank.outlet_to_pipeline.flow_mol[0] <= m.h2_demand / MW_H2

# Operating costs and revenue expressions at each hour
for (t, d, y) in m.set_period:
    append_op_costs_and_revenue(
        m.period[t, d, y],
        lmp=lmp_deterministic[y][d][t],
        h2_price=m.h2_price,
    )

# Append overall cashflows
append_cashflows(m)

# Define the objective function
m.obj = Objective(expr=m.npv, sense=maximize)

# Define the solver object. Using IPOPT
solver = get_solver()

# Solver the optimization problem
solver.solve(m, tee=True)
Ipopt 3.14.7: nlp_scaling_method=gradient-based
tol=1e-06


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.7, running with linear solver MUMPS 5.5.0.

Number of nonzeros in equality constraint Jacobian...:   580680
Number of nonzeros in inequality constraint Jacobian.:     7680
Number of nonzeros in Lagrangian Hessian.............:   195840

Total number of variables............................:   211163
                     variables with only lower bounds:    20123
                variables with lower and upper bounds:   176640
                     variables with only upper bounds:        0
Total number of equality constraints.................:   208280
Total number of inequality constraints...............:     3840
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:     3840

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -1.2218545e+09 1.75e+06 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.4944642e+09 1.91e+05 2.08e+04  -1.0 2.80e+06    -  8.10e-04 8.91e-01f  1
   2 -1.5310591e+09 7.77e+03 1.27e+04  -1.0 5.00e+05    -  1.82e-03 9.90e-01f  1
   3 -1.5315194e+09 2.28e+02 5.69e+02  -1.0 9.98e+03    -  2.59e-01 9.71e-01h  1
   4 -1.5315313e+09 1.30e+02 4.13e+02  -1.0 1.98e+02    -  1.43e-01 4.23e-01h  1
   5 -1.5315387e+09 9.95e+01 1.37e+03  -1.0 7.29e+02    -  2.34e-01 3.57e-01f  1
   6 -1.5315650e+09 1.09e+00 2.06e+03  -1.0 1.87e+02    -  1.67e-01 1.00e+00f  1
   7 -1.5315922e+09 3.97e-01 3.29e+03  -1.0 1.89e+02    -  6.84e-01 1.00e+00f  1
   8 -1.5316800e+09 1.23e+00 7.60e+03  -1.0 1.14e+03    -  6.66e-01 1.00e+00f  1
   9 -1.5317263e+09 1.90e+02 1.28e+06  -1.0 4.72e+04    -  3.52e-02 5.15e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -1.5317313e+09 1.36e+00 7.99e+03  -1.0 1.83e+02  -4.0 1.18e-01 1.00e+00f  1
  11 -1.5321112e+09 3.78e-01 1.25e+03  -1.0 5.99e+03    -  5.05e-01 1.00e+00f  1
  12 -1.5332880e+09 4.00e-03 6.05e+01  -1.0 3.66e+03    -  1.31e-02 1.00e+00f  1
  13 -1.5350114e+09 1.45e-06 5.12e+01  -1.0 1.44e+04    -  7.24e-02 1.00e+00f  1
  14 -1.5377561e+09 2.15e-10 2.40e+02  -1.0 1.69e+04    -  1.19e-01 1.00e+00f  1
  15 -1.5426917e+09 1.47e-10 1.06e+03  -1.0 3.00e+04    -  1.27e-01 1.00e+00f  1
  16 -1.5524243e+09 3.14e-10 1.01e+03  -1.0 5.34e+04    -  1.22e-01 1.00e+00f  1
  17 -1.5727455e+09 5.53e-10 8.05e+02  -1.0 9.82e+04    -  1.19e-01 1.00e+00f  1
  18 -1.5945127e+09 6.01e-10 7.26e+02  -1.0 1.91e+05    -  1.13e-01 5.05e-01f  1
  19 -1.5947343e+09 6.01e-10 6.15e+02  -1.0 1.76e+05    -  1.50e-01 7.45e-03f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -1.5947823e+09 6.64e-10 5.26e+02  -1.0 1.02e+05    -  1.43e-01 6.31e-02f  1
  21 -1.5954440e+09 5.01e-10 4.69e+02  -1.0 1.11e+05    -  1.17e-01 7.52e-01f  1
  22 -1.5955818e+09 5.97e-10 4.14e+02  -1.0 8.63e+04    -  1.15e-01 2.29e-01f  1
  23 -1.5959467e+09 4.94e-10 3.19e+02  -1.0 8.59e+04    -  2.25e-01 6.63e-01f  1
  24 -1.5962946e+09 7.00e-10 3.24e+02  -1.0 8.62e+04    -  2.53e-01 9.94e-01f  1
  25 -1.5964117e+09 5.98e-10 3.92e+02  -1.0 8.40e+04    -  2.73e-01 1.00e+00f  1
  26 -1.5964845e+09 6.68e-10 2.86e+02  -1.0 8.50e+04    -  2.86e-01 1.00e+00f  1
  27 -1.5965329e+09 5.81e-10 1.99e+02  -1.0 8.48e+04    -  2.98e-01 1.00e+00f  1
  28 -1.5965524e+09 5.92e-10 1.44e+02  -1.0 8.53e+04    -  2.80e-01 5.67e-01f  1
  29 -1.5965611e+09 5.03e-10 2.31e+02  -1.0 8.72e+04    -  1.90e-01 2.19e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30 -1.5965750e+09 5.61e-10 3.28e+02  -1.0 9.18e+04    -  2.21e-01 3.60e-01f  1
  31 -1.5965878e+09 4.92e-10 2.90e+02  -1.0 9.17e+04    -  2.29e-01 3.99e-01f  1
  32 -1.5965937e+09 5.81e-10 2.11e+02  -1.0 9.06e+04    -  2.50e-01 4.54e-01f  1
  33 -1.5965974e+09 4.24e-10 1.66e+02  -1.0 7.97e+04    -  2.71e-01 5.12e-01f  1
  34 -1.5966027e+09 4.51e-10 1.54e+02  -1.0 6.07e+04    -  2.84e-01 5.30e-01f  1
  35 -1.5966013e+09 5.31e-10 1.57e+02  -1.0 8.43e+04    -  2.95e-01 6.23e-01f  1
  36 -1.5965947e+09 5.06e-10 1.71e+02  -1.0 6.27e+04    -  3.32e-01 7.29e-01f  1
  37 -1.5965872e+09 5.72e-10 2.12e+02  -1.0 6.84e+04    -  3.75e-01 9.34e-01f  1
  38 -1.5965795e+09 5.43e-10 1.65e+02  -1.0 3.69e+04    -  4.64e-01 1.00e+00f  1
  39 -1.5965702e+09 6.88e-10 1.26e+02  -1.0 3.58e+04    -  6.60e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40 -1.5965659e+09 6.04e-10 5.90e+01  -1.0 2.92e+04    -  9.90e-01 1.00e+00f  1
  41 -1.5965644e+09 5.18e-10 2.28e+00  -1.0 7.96e+03    -  1.00e+00 1.00e+00f  1
  42 -1.5965644e+09 5.34e-10 6.54e-02  -1.0 1.30e+03    -  1.00e+00 1.00e+00f  1
  43 -1.5967508e+09 5.69e-10 4.41e+01  -2.5 1.54e+05    -  8.91e-01 9.82e-01f  1
  44 -1.5967560e+09 5.19e-10 2.11e+00  -2.5 6.73e+03    -  1.00e+00 1.00e+00f  1
  45 -1.5967560e+09 6.09e-10 4.61e-01  -2.5 3.47e+01    -  1.00e+00 1.00e+00f  1
  46 -1.5967560e+09 5.40e-10 1.10e-02  -2.5 5.35e+00    -  1.00e+00 1.00e+00f  1
  47 -1.5967588e+09 5.93e-10 4.13e+01  -3.8 4.45e+03    -  1.00e+00 5.56e-01f  1
  48 -1.5967599e+09 5.29e-10 1.09e+02  -3.8 1.97e+03    -  1.00e+00 5.49e-01f  1
  49 -1.5967604e+09 5.43e-10 1.31e+02  -3.8 8.72e+02    -  9.22e-01 6.43e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50 -1.5967605e+09 6.01e-10 2.67e+01  -3.8 2.83e+02    -  1.00e+00 9.38e-01f  1
  51 -1.5967605e+09 6.18e-10 1.09e-01  -3.8 5.05e+01    -  1.00e+00 1.00e+00f  1
  52 -1.5967605e+09 5.80e-10 5.86e-03  -3.8 1.33e+01    -  1.00e+00 1.00e+00f  1
  53 -1.5967605e+09 6.08e-10 2.15e-05  -3.8 8.11e-01    -  1.00e+00 1.00e+00h  1
  54 -1.5967606e+09 5.74e-10 1.36e+02  -5.7 1.54e+02    -  8.46e-01 3.87e-01f  1
  55 -1.5967606e+09 6.02e-10 9.27e+01  -5.7 4.71e+01    -  1.92e-01 2.71e-01f  1
  56 -1.5967606e+09 5.03e-10 6.26e+01  -5.7 2.95e+01    -  3.27e-01 3.26e-01f  1
  57 -1.5967606e+09 5.81e-10 7.87e+01  -5.7 1.90e+01    -  6.07e-01 1.02e-01f  1
  58 -1.5967606e+09 5.48e-10 4.90e+01  -5.7 1.38e+01    -  7.60e-02 2.99e-01f  1
  59 -1.5967606e+09 5.15e-10 2.67e+01  -5.7 9.27e+00    -  2.21e-01 3.76e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  60 -1.5967606e+09 5.25e-10 1.68e+01  -5.7 5.52e+00    -  4.97e-01 4.14e-01f  1
  61 -1.5967606e+09 4.94e-10 1.06e+01  -5.7 3.07e+00    -  8.68e-01 4.77e-01f  1
  62 -1.5967606e+09 5.41e-10 2.59e+00  -5.7 1.54e+00    -  1.00e+00 7.53e-01f  1
  63 -1.5967606e+09 6.61e-10 1.71e-05  -5.7 3.76e-01    -  1.00e+00 1.00e+00f  1
  64 -1.5967606e+09 6.88e-10 1.33e+00  -7.0 5.79e-01    -  9.95e-01 4.19e-01f  1
  65 -1.5967606e+09 6.08e-10 7.28e-03  -7.0 3.32e-01    -  1.00e+00 9.95e-01f  1
  66 -1.5967606e+09 6.74e-10 3.49e-09  -7.0 1.61e-03    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 66

                                   (scaled)                 (unscaled)
Objective...............:  -2.6787832953601610e+07   -1.5967606334020660e+09
Dual infeasibility......:   3.4906158125004671e-09    2.0806751801780723e-07
Constraint violation....:   1.1641532182693481e-10    6.7401506598230288e-10
Complementarity.........:   9.1160146931071113e-08    5.4338450671543375e-06
Overall NLP error.......:   9.1160146931071113e-08    5.4338450671543375e-06


Number of objective function evaluations             = 67
Number of objective gradient evaluations             = 67
Number of equality constraint evaluations            = 67
Number of inequality constraint evaluations          = 67
Number of equality constraint Jacobian evaluations   = 67
Number of inequality constraint Jacobian evaluations = 67
Number of Lagrangian Hessian evaluations             = 66
Total CPU secs in IPOPT (w/o function evaluations)   =     29.894
Total CPU secs in NLP function evaluations           =     17.501

EXIT: Optimal Solution Found.
[14]:
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 212120, 'Number of variables': 211163, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 50.57620286941528}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}

Verify that the solver converges to the optimal solution. We now print the results and answer the questions posed at the beginning of this notebook. To recall, we are interested in determining whether the following options are attractive (i.e., maximize NPV):

  • Produce hydrogen using a portion of the electricity, especially during periods of low electricity demand, and sell it.

  • Produce hydrogen, store it in a tank, and combust it to produce electrity during the periods of high electricity demand.

[15]:
def generate_plots(m, d, y, set_time, lmp):
    LMP = [lmp[y][d][t] for t in set_time]

    # Power from nuclear power plant to the grid (convert it to MW)
    np_to_grid = [m.period[t, d, y].fs.np_power_split.split_fraction["np_to_grid", 0].value
                  for t in set_time]
    # Hydrogen production rate in the PEM electrolyzer (convert it to kg/hr)
    h2_production = [m.period[t, d, y].fs.pem.outlet.flow_mol[0].value * 2.016e-3 * 3600
                     for t in set_time]
    # Hydrogen flowrate to pipeline (convert it to kg/hr)
    h2_to_pipeline = [m.period[t, d, y].fs.h2_tank.outlet_to_pipeline.flow_mol[0].value * 2.016e-3 * 3600
                      for t in set_time]
    # Hydrogen flowrate to turbine (convert it to kg/hr)
    h2_to_turbine = [m.period[t, d, y].fs.h2_tank.outlet_to_turbine.flow_mol[0].value * 2.016e-3 * 3600
                     for t in set_time]

    # Plot the results
    if hasattr(m, "plot_lmp_and_schedule"):
        plot_lmp_and_schedule = m.plot_lmp_and_schedule

    else:
        plot_lmp_and_schedule = m.parent_block().plot_lmp_and_schedule

    plot_lmp_and_schedule(
        lmp=LMP,
        schedule={"power_to_grid": np_to_grid,
                  "h2_production": h2_production,
                  "h2_to_pipeline": h2_to_pipeline,
                  "h2_to_turbine": h2_to_turbine},
        y_label={"power_to_grid": "Split fraction to grid [-]",
                 "h2_production": "Hydrogen production (kg/hr)",
                 "h2_to_pipeline": "Hydrogen to pipeline (kg/hr)",
                 "h2_to_turbine": "Hydrogen to turbine (kg/hr)"},
        y_range={"power_to_grid": (0.5, 1.02),
                 "h2_production": (0, 4000),
                 "h2_to_pipeline": (0, 4000),
                 "h2_to_turbine": (-0.5, 10)},
    )
[16]:
# Print Results
print("Optimal PEM capacity    : ", m.pem_capacity.value * 1e-3, "MW")
print("Optimal tank capacity   : ", m.tank_capacity.value * 2.016e-3, "kg")
print("OPtimal turbine capacity: ", m.h2_turbine_capacity.value * 1e-6, "MW")

# Using `generate_plots` we plot optimal operating conditions for a specific day
generate_plots(m, d=1, y=2022, set_time=m.set_time, lmp=lmp_deterministic)
Optimal PEM capacity    :  196.2616021908122 MW
Optimal tank capacity   :  2.0756993639026943e-05 kg
OPtimal turbine capacity:  9.076695580817342e-05 MW
../_images/examples_multiperiod_design_pricetaker_30_1.png
[17]:
# Using `generate_plots` we plot optimal operating conditions for a specific day
generate_plots(m, d=8, y=2022, set_time=m.set_time, lmp=lmp_deterministic)
../_images/examples_multiperiod_design_pricetaker_31_0.png
  • Clearly, when the selling price of hydrogen is $3/kg, it is attractive (maximizes NPV) to produce hydrogen using a portion of electricity when its price is sufficiently low.

  • However, storing hydrogen and combusting it later to produce electricity is not attractive. This is inferred from the optimal size of the hydrogen tank and the size of the turbine, respectively (both the values are zero. We see a small nonzero value because we imposed small nonzero lower bounds on flows to avoid convergence issues).

The optimal size of the PEM electrolyzer for the production of hydrogen is 196.26 MW for a hydrogen demand of 1 kg/s. The hydrogen tank and the hydrogen turbine must not be built to maximize the NPV. The hydrogen produced by the PEM electrolyzer must be sold immediately to the hydrogen market via the pipeline.

[18]:
# These optimization problems tend to be very large, so we
# delete the model after analyzing the results to save memory.
del m

Multiperiod Optimization Model: Stochastic

LMP Signal

As mentioned in the previous section, due to uncertainty in various factors such as weather, demand, etc., it is not possible to determine the future locational marginal price accurately. There are many approaches to take into account the uncertainty in the LMP signal during the decision making process. One such approach involves the generation of potential LMP scenarios along with their associated probabilities and use them to formulate a stochastic program.

Here, for demonstration, we consider two different realizations of the LMP signal (i.e., two scenarios: scenario 0 and scenario 1).

[19]:
# Gather the LMP data needed for the stochastic case
# Notation: lmp_dataset[scenario][year][cluster/day][time/hour]
lmp_stochastic = {scenario: {year: {cluster: {hour: lmp_dataset[str(scenario)][str(year)][str(cluster)][str(hour)]
                                              for hour in range(1, 25)}
                                    for cluster in range(1, 21)}
                             for year in [2022]}
                  for scenario in [0, 1]}

Next, we formulate an optimization problem of the form

\[\begin{split} \begin{aligned} \max_{D, D_s, u_{s, t, d}} \quad & \sum_s w_s \text{NPV}_s(D_s, u_{s, t, d})\\ & g(u_{s, t, d}) = 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C}; \forall\; s \in \mathcal{S} \\ & h(u_{s, t, d}) \le 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C}; \forall\; s \in \mathcal{S} \\ & f(u_{s, t-1, d}, u_{s, t, d}) = 0, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C}; \forall\; s \in \mathcal{S} \\ & u_{s, t, d} \le D_s, & & \forall \; t \in \mathcal{T}; \; \forall \; d \in \mathcal{C}; \forall\; s \in \mathcal{S}\\ & D_s = D, & & \forall\; s \in \mathcal{S} \end{aligned}\end{split}\]

Here, the sets \(\mathcal{S} = \{0, 1\}\), \(\mathcal{T} = \{1, \dots, 24\}\) and \(\mathcal{C} = \{1, \dots, 20\}\) denote the the set of scenarios, the set of hours in a day and the set of clusters/days in a year, respectively. \(D\) denotes the design decisions (or, first-stage variables) viz. size of the PEM electrolyzer (pem_capacity), tank (tank_capacity) and the turbine (turbine_capacity), and \(D_s\) denotes the design decisions in scenario \(s\). \(u_{s, t, d}\) denotes the operating decisions (or, second-stage decisions) in scenario \(s\), at time \(t\) of day/cluster \(d\). \(g(u_{s, t, d}) = 0\) and \(h(u_{s, t,d}) \le 0\) denote the flowsheet model, \(f(u_{s, t-1,d}, u_{s, t,d}) = 0\) connects the operating decisions at time \(t-1\) and at time \(t\) (e.g., tank holdup), and \(u_{s, t, d} \le D_s\) ensures that the operating decision values never exceed the design capacity values (e.g., tank holdup at any time instant must not exceed the tank capacity) in scenario \(s\). Finally, \(D_s = D\) ensures that the design decisions are maintained the same in all scenarios (non-anticipativity constraints).

In the objective function, \(w_s\) denotes the probability associated with scenario \(s\), and \(\text{NPV}_s\) stands for the net present value calculated in scenario \(s\).

As before, we use the MultiPeriodModel class to formulate the stochastic multi-period optimization problem. In addition to all the arguments mentioned before, we also specify set_scenarios, the set of scenarios, as [0, 1].

[20]:
# Formulate the stochastic multiperiod optimization problem
m = MultiPeriodModel(
    n_time_points=24,
    set_days=[i for i in range(1, 21)],
    set_years=[2022],
    set_scenarios=[0, 1],
    process_model_func=build_ne_flowsheet,
    initialization_func=fix_dof_and_initialize,
    unfix_dof_func=unfix_dof,
    linking_variable_func=get_linking_variable_pairs,
    flowsheet_options={"np_capacity": 1000},
    initialization_options={
        "split_frac_grid": 0.8,
        "tank_holdup_previous": 0,
        "flow_mol_to_pipeline": 10,
        "flow_mol_to_turbine": 10,
    },
    use_stochastic_build=True,
    outlvl=logging.WARNING,
)
[+   0.00] Beginning the formulation of the multiperiod optimization problem.
[+ 190.23] Completed the formulation of the multiperiod optimization problem.
2022-10-08 11:13:14 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
[+   7.36] Created an instance of the flowsheet and initialized it.
[+  13.78] Initialized the entire multiperiod optimization model.
[+   0.44] Unfixed the degrees of freedom from each period model.
[21]:
# Define parameters
m.plant_life = 20                            # Plant lifetime: 20 years
m.tax_rate = 0.2                             # Corporate tax rate: 20%
m.discount_rate = 0.08                       # Discount rate: 8%
m.h2_price = 3                               # Selling price of hydrogen: $3/kg
m.h2_demand = 1                              # Maximum amount of hydrogen that can be sold: 1 kg/s
m.LMP = lmp_stochastic                       # LMP signal
m.weights_days = weights_days                # number of days represented by each cluster
m.weights_scenarios = {0: 0.5, 1: 0.5}       # Equal probability for both the scenarios

Next, add the connecting constraints and hydrogen demand constraints for each of the scenarios.

[22]:
for s in m.set_scenarios:
    # Add first-stage variables to each scenario
    add_capacity_variables(m.scenario[s])

    # Set initial holdup for each day (Assumed to be zero at the beginning of each day)
    for y in m.set_years:
        for d in m.set_days:
            m.scenario[s].period[1, d, y].fs.h2_tank.tank_holdup_previous.fix(0)

    # Operating costs and revenue expressions at each hour
    for (t, d, y) in m.set_period:
        append_op_costs_and_revenue(
            m.scenario[s].period[t, d, y],
            lmp=lmp_stochastic[s][y][d][t],
            h2_price=m.h2_price,
        )

    # Append overall cashflows
    append_cashflows(m.scenario[s])

    # Hydrogen demand constraint.
    # Divide the RHS by the molecular mass to convert kg/s to mol/s
    obj = m.scenario[s]
    @obj.Constraint(m.set_time, m.set_days, m.set_years)
    def hydrogen_demand_constraint(blk, t, d, y):
        return blk.period[t, d, y].fs.h2_tank.outlet_to_pipeline.flow_mol[0] <= m.h2_demand / MW_H2
[23]:
# Add non-anticipativity constraints
m.pem_capacity = Var(within=NonNegativeReals,
                     doc="Design PEM capacity (in kW)",
                     units=pyunits.kW)
m.tank_capacity = Var(within=NonNegativeReals,
                      doc="Design tank capacity (in mol)",
                      units=pyunits.mol)
m.h2_turbine_capacity = Var(within=NonNegativeReals,
                            doc="Design turbine capacity (in W)",
                            units=pyunits.W)

@m.Constraint(m.set_scenarios)
def non_anticipativity_pem(blk, s):
    return blk.pem_capacity == blk.scenario[s].pem_capacity

@m.Constraint(m.set_scenarios)
def non_anticipativity_tank(blk, s):
    return blk.tank_capacity == blk.scenario[s].tank_capacity

@m.Constraint(m.set_scenarios)
def non_anticipativity_turbine(blk, s):
    return blk.h2_turbine_capacity == blk.scenario[s].h2_turbine_capacity
[24]:
# Define the objective function
m.obj = Objective(expr=sum(m.weights_scenarios[s] * m.scenario[s].npv
                           for s in m.set_scenarios),
                  sense=maximize)

# Define the solver object. Using the default solver: IPOPT
solver = get_solver()

# Solve the optimization problem
solver.solve(m, tee=True)
Ipopt 3.14.7: nlp_scaling_method=gradient-based
tol=1e-06


******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

This is Ipopt version 3.14.7, running with linear solver MUMPS 5.5.0.

Number of nonzeros in equality constraint Jacobian...:   580692
Number of nonzeros in inequality constraint Jacobian.:     7680
Number of nonzeros in Lagrangian Hessian.............:   195840

Total number of variables............................:   211169
                     variables with only lower bounds:    20129
                variables with lower and upper bounds:   176640
                     variables with only upper bounds:        0
Total number of equality constraints.................:   208286
Total number of inequality constraints...............:     3840
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:     3840

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0 -1.2559159e+09 1.75e+06 1.00e+02  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1 -1.5355493e+09 1.94e+05 2.06e+04  -1.0 2.80e+06    -  2.43e-03 8.89e-01f  1
   2 -1.5732609e+09 7.89e+03 1.24e+04  -1.0 5.04e+05    -  6.39e-04 9.90e-01f  1
   3 -1.5737278e+09 2.33e+02 5.62e+02  -1.0 1.01e+04    -  4.84e-01 9.71e-01f  1
   4 -1.5737296e+09 2.21e+02 5.34e+02  -1.0 2.04e+02    -  2.69e-02 5.23e-02f  4
   5 -1.5737373e+09 1.69e+02 4.32e+02  -1.0 1.39e+02    -  7.09e-02 2.32e-01h  2
   6 -1.5737516e+09 9.33e+01 2.73e+02  -1.0 9.11e+01    -  8.23e-01 4.43e-01h  1
   7 -1.5738919e+09 9.07e+00 5.54e+02  -1.0 5.28e+02    -  4.10e-01 9.39e-01f  1
   8 -1.5741711e+09 2.33e+01 8.82e+03  -1.0 1.02e+03    -  3.86e-01 1.00e+00f  1
   9 -1.5741718e+09 2.03e+01 7.71e+03  -1.0 1.02e+02  -4.0 3.95e-01 1.25e-01f  4
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10 -1.5741812e+09 8.83e+00 3.18e+03  -1.0 2.88e+02  -4.5 1.95e-01 5.00e-01f  2
  11 -1.5742405e+09 1.20e+00 7.18e+02  -1.0 4.64e+02  -5.0 1.67e-01 1.00e+00f  1
  12 -1.5758575e+09 1.61e+00 7.49e+02  -1.0 1.25e+04    -  2.49e-02 1.00e+00f  1
  13 -1.5779187e+09 3.22e-02 1.63e+02  -1.0 1.69e+04    -  1.13e-01 6.84e-01f  1
  14 -1.5835775e+09 1.34e-02 1.78e+03  -1.0 2.97e+04    -  1.00e-01 9.54e-01f  1
  15 -1.5978308e+09 4.90e-03 6.72e+03  -1.0 6.24e+04    -  1.06e-01 1.00e+00f  1
  16 -1.6303589e+09 2.75e-04 1.26e+03  -1.0 1.53e+05    -  9.86e-02 8.53e-01f  1
  17 -1.6306844e+09 2.73e-04 1.19e+03  -1.0 2.14e+05    -  1.40e-01 6.28e-03f  1
  18 -1.6306914e+09 2.65e-04 1.02e+03  -1.0 4.69e+04    -  4.10e-01 3.19e-02f  1
  19 -1.6307897e+09 3.96e-05 4.32e+02  -1.0 6.76e+04    -  3.12e-01 6.44e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -1.6308760e+09 6.12e-06 2.08e+02  -1.0 6.58e+04    -  4.37e-01 8.21e-01f  1
  21 -1.6309292e+09 2.98e-06 1.35e+02  -1.0 6.00e+04    -  5.55e-01 1.00e+00f  1
  22 -1.6309356e+09 1.34e-06 1.77e+02  -1.0 5.99e+04    -  1.96e-01 4.56e-01f  1
  23 -1.6309409e+09 7.57e-07 1.22e+02  -1.0 5.16e+04    -  2.85e-01 4.97e-01f  1
  24 -1.6309464e+09 2.32e-07 2.00e+02  -1.0 5.17e+04    -  3.04e-01 8.46e-01f  1
  25 -1.6309419e+09 4.16e-08 1.76e+02  -1.0 4.01e+04    -  3.49e-01 9.11e-01f  1
  26 -1.6309368e+09 1.38e-08 1.92e+02  -1.0 3.75e+04    -  3.65e-01 1.00e+00f  1
  27 -1.6309299e+09 5.02e-09 1.47e+02  -1.0 3.24e+04    -  4.05e-01 1.00e+00f  1
  28 -1.6309221e+09 1.50e-09 1.26e+02  -1.0 1.99e+04    -  4.82e-01 1.00e+00f  1
  29 -1.6309164e+09 6.38e-10 9.58e+01  -1.0 2.01e+04    -  5.89e-01 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  30 -1.6309105e+09 5.24e-10 1.97e+01  -1.0 1.17e+04    -  9.90e-01 1.00e+00f  1
  31 -1.6309083e+09 6.80e-10 8.19e-01  -1.0 4.00e+03    -  9.97e-01 1.00e+00f  1
  32 -1.6310231e+09 5.76e-10 1.11e+01  -1.7 9.28e+04    -  1.00e+00 1.00e+00f  1
  33 -1.6310238e+09 6.89e-10 2.90e-02  -1.7 1.40e+03    -  1.00e+00 1.00e+00f  1
  34 -1.6310487e+09 6.66e-10 4.20e+00  -2.5 2.05e+04    -  1.00e+00 9.98e-01f  1
  35 -1.6310485e+09 6.35e-10 1.15e+01  -2.5 1.67e+02    -  1.00e+00 1.00e+00f  1
  36 -1.6310486e+09 5.79e-10 2.05e+00  -2.5 7.12e+01    -  1.00e+00 1.00e+00f  1
  37 -1.6310486e+09 5.54e-10 1.59e-01  -2.5 1.98e+01    -  1.00e+00 1.00e+00f  1
  38 -1.6310486e+09 5.65e-10 1.35e-03  -2.5 1.83e+00    -  1.00e+00 1.00e+00h  1
  39 -1.6310507e+09 5.38e-10 5.20e+01  -3.8 3.21e+03    -  1.00e+00 5.64e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  40 -1.6310513e+09 4.66e-10 1.38e+02  -3.8 1.39e+03    -  1.00e+00 5.46e-01f  1
  41 -1.6310516e+09 5.36e-10 1.77e+02  -3.8 6.08e+02    -  1.00e+00 6.43e-01f  1
  42 -1.6310517e+09 6.60e-10 6.39e+01  -3.8 1.72e+02    -  1.00e+00 8.53e-01f  1
  43 -1.6310517e+09 5.15e-10 5.47e-02  -3.8 2.77e+01    -  1.00e+00 1.00e+00f  1
  44 -1.6310517e+09 5.44e-10 1.40e-03  -3.8 6.39e+00    -  1.00e+00 1.00e+00f  1
  45 -1.6310517e+09 5.62e-10 6.58e+01  -5.7 8.60e+01    -  5.85e-01 3.77e-01f  1
  46 -1.6310518e+09 5.11e-10 4.45e+01  -5.7 3.12e+01    -  3.37e-01 3.32e-01f  1
  47 -1.6310518e+09 5.34e-10 7.57e+01  -5.7 1.98e+01    -  4.51e-01 1.05e-01f  1
  48 -1.6310518e+09 5.69e-10 4.02e+01  -5.7 1.53e+01    -  9.80e-02 2.73e-01f  1
  49 -1.6310518e+09 5.18e-10 1.88e+01  -5.7 9.71e+00    -  2.71e-01 4.37e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  50 -1.6310518e+09 5.36e-10 1.45e+01  -5.7 5.01e+00    -  5.83e-01 3.82e-01f  1
  51 -1.6310518e+09 5.54e-10 1.12e+01  -5.7 2.88e+00    -  9.20e-01 4.72e-01f  1
  52 -1.6310518e+09 5.29e-10 4.71e+00  -5.7 1.45e+00    -  1.00e+00 5.95e-01f  1
  53 -1.6310518e+09 5.57e-10 1.07e+00  -5.7 5.87e-01    -  1.00e+00 7.75e-01f  1
  54 -1.6310518e+09 6.13e-10 2.96e-06  -5.7 1.32e-01    -  1.00e+00 1.00e+00f  1
  55 -1.6310518e+09 5.26e-10 8.59e-01  -7.0 4.10e-01    -  1.00e+00 6.22e-01f  1
  56 -1.6310518e+09 5.33e-10 3.09e-06  -7.0 1.53e-01    -  1.00e+00 1.00e+00f  1
  57 -1.6310518e+09 5.52e-10 2.57e-12  -7.0 2.88e-04    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 57

                                   (scaled)                 (unscaled)
Objective...............:  -3.7401907403570741e+07   -1.6310517839609084e+09
Dual infeasibility......:   2.5730052496087487e-12    1.1220563585787805e-10
Constraint violation....:   1.1641532182693481e-10    5.5245674701609460e-10
Complementarity.........:   9.0916704112864471e-08    3.9647671129461297e-06
Overall NLP error.......:   9.0916704112864471e-08    3.9647671129461297e-06


Number of objective function evaluations             = 66
Number of objective gradient evaluations             = 58
Number of equality constraint evaluations            = 66
Number of inequality constraint evaluations          = 66
Number of equality constraint Jacobian evaluations   = 58
Number of inequality constraint Jacobian evaluations = 58
Number of Lagrangian Hessian evaluations             = 57
Total CPU secs in IPOPT (w/o function evaluations)   =     27.447
Total CPU secs in NLP function evaluations           =     15.755

EXIT: Optimal Solution Found.
[24]:
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 212126, 'Number of variables': 211169, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 46.29058003425598}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
[25]:
# Print Results
print("Optimal PEM capacity    : ", m.pem_capacity.value * 1e-3, "MW")
print("Optimal tank capacity   : ", m.tank_capacity.value * 2.016e-3, "kg")
print("OPtimal turbine capacity: ", m.h2_turbine_capacity.value * 1e-6, "MW")
Optimal PEM capacity    :  196.26160219185448 MW
Optimal tank capacity   :  1.4179085730078668e-05 kg
OPtimal turbine capacity:  9.076374619487333e-05 MW
[26]:
generate_plots(m.scenario[0], d=8, y=2022, set_time=m.set_time, lmp=m.LMP[0])
../_images/examples_multiperiod_design_pricetaker_45_0.png

Even in the stochastic case, we observe the same result i.e., producing hydrogen and selling it during the periods of low electricity demand is profitable. Whereas, building a turbine to combust hydrogen and produce electricity is not attractive.