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.
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
SimplifiedHydrogenTankmodel 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
# Pyomo imports
from pyomo.environ import (ConcreteModel,
RangeSet,
Var,
NonNegativeReals,
Constraint,
Expression,
Objective,
maximize,
units as pyunits)
# IDAES imports
from idaes.core.solvers import get_solver
from idaes.core.util.model_statistics import degrees_of_freedom
# Nuclear flowsheet function imports
from nuclear_flowsheet import (build_ne_flowsheet,
fix_dof_and_initialize)
# Import function for the construction of the multiperiod model
from multiperiod import (build_multiperiod_design,
plot_lmp_and_schedule)
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
|
|
24.01 bar |
Isentropic
efficiency of
|
|
0.86 |
Conversion of
hydrogen in
|
|
0.99 |
Pressure deifference
across
|
|
-24.01 bar |
Isentripic
efficiency of
|
|
0.89 |
Molar flow rate of
air to
|
|
10.76 * molar flowrate of hydrogen to turbine |
Temperature of air |
|
300 K |
Pressure of air |
|
1.01325 bar |
Duration of the
simulation for
|
|
3600 s |
[2]:
# Create a concrete model object
m = ConcreteModel()
# Build the nuclear flowsheet
build_ne_flowsheet(m)
# Fix the degrees of freedom and initialize
fix_dof_and_initialize(m)
# 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-07-01 17:07:21 [INFO] idaes.init.fs.pem.outlet_state: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.pem.outlet_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_in: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_tank: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator.properties_in: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator.properties_out: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.translator: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.air_feed_state: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.air_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.mixed_state: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.mixed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer.mixed_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.mixer: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine: Starting initialization...
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.compressor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Starting initialization
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:21 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume: Initialization Complete
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.stoic_reactor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Starting initialization
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Starting initialization
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Starting initialization
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine.turbine: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:07:22 [INFO] idaes.init.fs.h2_turbine: Initialization complete
Degrees of freedom: 0
Ipopt 3.13.2: 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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
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 8.83e+00 -1.0 4.05e+02 - 1.96e-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....: 4.6566128730773926e-10 3.4924596548080444e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 4.6566128730773926e-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.002
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.024416446685791016}], '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 : ", (- m.fs.h2_turbine.turbine.work_mechanical[0].value
- m.fs.h2_turbine.compressor.work_mechanical[0].value),
get_units(m.fs.h2_turbine.compressor.work_mechanical[0]))
Nuclear plant power production : 1000000.0 kW
Electricity to grid : 800000.0 kW
Electricity to PEM : 200000.0 kW
Flowrate of H2 from pem : 505.48119999999994 mol/s
Flowrate of H2 to tank : 505.48119999999994 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.3199999998 mol
H2 Turbine's compressor work : 1764604.5955977982 kg*m**2/s**3
H2 Turbine's turbine work : -2672169.126200439 kg*m**2/s**3
Net power produced by turbine : 907564.5306026407 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. Owing to the uncertain nature of some those 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 FORCE (link goes here) 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.
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
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).
First, we create a ConcreteModel object, and define sets and parameters
[6]:
# Create a ConcreteModel object
m = ConcreteModel()
# Define sets
m.set_time = RangeSet(24) # twenty fours in a day
m.set_days = RangeSet(20) # twenty clusters/days per year
m.set_years = [2022, 2032] # Set of years
# 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_deterministic # LMP signal
m.weights_days = weights_days # number of days represented by each cluster
As we saw earlier, functions build_ne_flowsheet and fix_dof_and_initialize yield an initialized square model of the nuclear flowsheet. For the multi-period optimization problem, we now construct a function that can be used to unfix a few degrees of freedom for optimization. The function below
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.
[7]:
def unfix_dof(m, **kwargs):
"""
This function unfixes a few degrees of freedom for optimization
"""
# Set defaults in case options are not passed to the function
options = kwargs.get("options", {})
air_h2_ratio = options.get("air_h2_ratio", 10.76)
# 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)
Next, we use the build_multiperiod_design function 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}\). This function 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), and unfixes operational degrees of freedom (using the unfix_dof
function) described above. Since we are dealing with multiple days/clusters and multiple years (2022 and 2032 LMP signals), we set the multiple_days and multiyear argument to True (see below. Note that the default value of these arguments is False).
Note: The function, by default, looks for m.set_time, m.set_days and m.set_years objects for the set of hours, days/clusters and years, respectively. However, if the sets are created using a different set of names, then they must be passed to the function explicitly. E.g. build_multiperiod_design(m, ..., set_time=m.set_time_object, set_days=m.set_days_object, set_years=m.set_years_object).
Tip: Building the multiperiod model can take some time. So, set verbose=True to moniter the progress
[8]:
# Build initialized multiperiod model: Constraints g(u(t,d)) = 0 and h(u(t,d)) <= 0
build_multiperiod_design(m,
flowsheet=build_ne_flowsheet,
initialization=fix_dof_and_initialize,
unfix_dof=unfix_dof,
multiple_days=True,
multiyear=True,
verbose=False)
[ 0.00] Processing input information.
[+ 0.00] Beginning the formulation of the multiperiod problem.
[+ 226.00] Completed the formulation of the multiperiod problem
2022-07-01 17:11:08 [INFO] idaes.init.fs.pem.outlet_state: Starting initialization
2022-07-01 17:11:08 [INFO] idaes.init.fs.pem.outlet_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:08 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:08 [INFO] idaes.init.fs.h2_tank.properties_in: Starting initialization
2022-07-01 17:11:08 [INFO] idaes.init.fs.h2_tank.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:08 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_tank: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator.properties_in: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator.properties_out: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.translator: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.air_feed_state: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.air_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.mixed_state: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.mixed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer.mixed_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.mixer: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_turbine: Starting initialization...
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Starting initialization
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:09 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.compressor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume: Initialization Complete
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.stoic_reactor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Starting initialization
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:10 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:11:11 [INFO] idaes.init.fs.h2_turbine.turbine: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:11:11 [INFO] idaes.init.fs.h2_turbine: Initialization complete
[+ 2.85] Created an instance of the flowsheet and initialized it.
[+ 24.06] Initialized the entire multiperiod optimization model.
[+ 0.47] Unfixed the degrees of freedom from each period model.
Next, we define
Function
build_connecting_constraintsto declare the design decisions and to build the constraints \(f(u_{t-1,d}, u_{t,d}) = 0\) and \(u_{t,d} \le D\).Function
append_costs_and_revenueto build expressions for all the cash flows (electricity revenue, hydrogen revenue, fixed O&M, variable 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.
Note: It is not necessary to define the above functions for the deterministic case. Nevertheless, we do so because it will be useful for the stochastic case, as we shall see in the next section.
[9]:
def build_connecting_constraints(m, set_time, set_days, set_years):
"""
This function declares the first-stage variables or design decisions,
adds constraints that ensure that the operational variables never exceed their
design values, and adds constraints connecting variables at t - 1 and t
"""
# Declare first-stage variables (Design decisions)
m.pem_capacity = Var(within=NonNegativeReals,
doc="Maximum capacity of the PEM electrolyzer (in kW)")
m.tank_capacity = Var(within=NonNegativeReals,
doc="Maximum holdup of the tank (in mol)")
m.h2_turbine_capacity = Var(within=NonNegativeReals,
doc="Maximum power output from the turbine (in W)")
# Ensure that the electricity to the PEM elctrolyzer does not exceed the PEM capacity
@m.Constraint(set_time, set_days, set_years)
def pem_capacity_constraint(blk, t, d, y):
return blk.period[t, d, y].fs.pem.electricity[0] <= m.pem_capacity
# Ensure that the final tank holdup does not exceed the tank capacity
@m.Constraint(set_time, set_days, set_years)
def tank_capacity_constraint(blk, t, d, y):
return blk.period[t, d, y].fs.h2_tank.tank_holdup[0] <= m.tank_capacity
# Ensure that the power generated by the turbine does not exceed the turbine capacity
@m.Constraint(set_time, set_days, set_years)
def turbine_capacity_constraint(blk, t, d, y):
return (
- blk.period[t, d, y].fs.h2_turbine.turbine.work_mechanical[0]
- blk.period[t, d, y].fs.h2_turbine.compressor.work_mechanical[0] <=
m.h2_turbine_capacity
)
# Connect the initial tank holdup at time t with the final tank holdup at time t - 1
@m.Constraint(set_time, set_days, set_years)
def tank_holdup_constraints(blk, t, d, y):
if t == 1:
# Each day begins with an empty tank
return (
blk.period[t, d, y].fs.h2_tank.tank_holdup_previous[0] == 0
)
else:
# Initial holdup at time t = final holdup at time t - 1
return (
blk.period[t, d, y].fs.h2_tank.tank_holdup_previous[0] ==
blk.period[t - 1, d, y].fs.h2_tank.tank_holdup[0]
)
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]:
def append_costs_and_revenue(m, ps, LMP):
"""
ps: Object containing information on sets and parameters
LMP: Dictionary containing the LMP data
"""
set_time = ps.set_time # Set of hours
set_days = ps.set_days # Set of days/clusters
set_years = ps.set_years # Set of years
weights_days = ps.weights_days # Weights associated with each cluster
h2_sp = ps.h2_price # Selling price of hydrogen
plant_life = ps.plant_life # Plant lifetime
tax_rate = ps.tax_rate # Corporate tax rate
discount_rate = ps.discount_rate # Discount rate
years_vec = [y - set_years[0] + 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=(1630 * m.pem_capacity +
(29 * 33.3 * 2.016e-3) * m.tank_capacity +
(947 / 1000) * m.h2_turbine_capacity),
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 (
47.9 * m.pem_capacity + 7e-3 * m.h2_turbine_capacity
)
# 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):
return (
(1.3 * 1e-3) * sum(weights_days[y][d] * blk.period[t, d, y].fs.pem.electricity[0]
for t in set_time for d in set_days) +
(4.25 * 1e-6) * sum(weights_days[y][d] * (
- blk.period[t, d, y].fs.h2_turbine.turbine.work_mechanical[0]
- blk.period[t, d, y].fs.h2_turbine.compressor.work_mechanical[0])
for t in set_time for d in set_days)
)
@m.Expression(set_years,
doc="Revenue generated by selling electricity per year (in USD)")
def electricity_revenue(blk, y):
return (
sum(weights_days[y][d] * LMP[y][d][t] *
(blk.period[t, d, y].fs.np_power_split.np_to_grid_port.electricity[0] * 1e-3 -
blk.period[t, d, y].fs.h2_turbine.turbine.work_mechanical[0] * 1e-6 -
blk.period[t, d, y].fs.h2_turbine.compressor.work_mechanical[0] * 1e-6)
for t in set_time for d in set_days)
)
@m.Expression(set_years,
doc="Revenue generated by selling hydrogen per year (in USD)")
def h2_revenue(blk, y):
return (
h2_sp * 2.016e-3 * 3600 *
sum(weights_days[y][d] *
blk.period[t, d, y].fs.h2_tank.outlet_to_pipeline.flow_mol[0]
for t in set_time for d in set_days)
)
@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)"
)
We now build the remaining constraints and the objective function by calling build_connecting_constraints and append_cost_and_revenue functions, and solve the optimization problem.
[11]:
# Build the connecting constraints
build_connecting_constraints(m,
set_time=m.set_time,
set_days=m.set_days,
set_years=m.set_years)
# 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 / 2.016e-3
# Append cash flow expressions
append_costs_and_revenue(m,
ps=m,
LMP=m.LMP)
# 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.13.2: 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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 580760
Number of nonzeros in inequality constraint Jacobian.: 7680
Number of nonzeros in Lagrangian Hessian.............: 195840
Total number of variables............................: 211203
variables with only lower bounds: 20163
variables with lower and upper bounds: 176640
variables with only upper bounds: 0
Total number of equality constraints.................: 208320
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 6.52e+03 -1.0 9.98e+03 - 2.59e-01 9.71e-01h 1
4 -1.5315313e+09 1.30e+02 8.88e+04 -1.0 1.98e+02 - 1.43e-01 4.23e-01h 1
5 -1.5315387e+09 9.95e+01 1.20e+05 -1.0 7.29e+02 - 2.34e-01 3.57e-01f 1
6 -1.5315650e+09 1.09e+00 6.39e+05 -1.0 1.87e+02 - 1.67e-01 1.00e+00f 1
7 -1.5315922e+09 3.97e-01 3.07e+06 -1.0 1.89e+02 - 6.84e-01 1.00e+00f 1
8 -1.5316800e+09 1.23e+00 1.03e+06 -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 8.72e+05 -1.0 1.83e+02 -4.0 1.18e-01 1.00e+00f 1
11 -1.5321112e+09 3.78e-01 4.32e+05 -1.0 5.99e+03 - 5.05e-01 1.00e+00f 1
12 -1.5332880e+09 4.00e-03 4.26e+05 -1.0 3.66e+03 - 1.31e-02 1.00e+00f 1
13 -1.5350114e+09 1.45e-06 3.95e+05 -1.0 1.44e+04 - 7.24e-02 1.00e+00f 1
14 -1.5377561e+09 2.18e-10 3.48e+05 -1.0 1.69e+04 - 1.19e-01 1.00e+00f 1
15 -1.5426917e+09 1.27e-10 3.04e+05 -1.0 3.00e+04 - 1.27e-01 1.00e+00f 1
16 -1.5524243e+09 3.43e-10 2.67e+05 -1.0 5.34e+04 - 1.22e-01 1.00e+00f 1
17 -1.5727455e+09 5.98e-10 2.35e+05 -1.0 9.82e+04 - 1.19e-01 1.00e+00f 1
18 -1.5945127e+09 6.02e-10 2.08e+05 -1.0 1.91e+05 - 1.13e-01 5.05e-01f 1
19 -1.5947343e+09 6.46e-10 1.77e+05 -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 7.89e-10 1.52e+05 -1.0 1.02e+05 - 1.43e-01 6.31e-02f 1
21 -1.5954440e+09 5.61e-10 1.34e+05 -1.0 1.11e+05 - 1.17e-01 7.52e-01f 1
22 -1.5955818e+09 5.49e-10 1.18e+05 -1.0 8.63e+04 - 1.15e-01 2.29e-01f 1
23 -1.5959467e+09 4.95e-10 9.17e+04 -1.0 8.59e+04 - 2.25e-01 6.63e-01f 1
24 -1.5962946e+09 5.55e-10 6.85e+04 -1.0 8.62e+04 - 2.53e-01 9.94e-01f 1
25 -1.5964117e+09 5.31e-10 4.99e+04 -1.0 8.40e+04 - 2.73e-01 1.00e+00f 1
26 -1.5964845e+09 6.88e-10 3.56e+04 -1.0 8.50e+04 - 2.86e-01 1.00e+00f 1
27 -1.5965329e+09 5.42e-10 2.50e+04 -1.0 8.48e+04 - 2.98e-01 1.00e+00f 1
28 -1.5965524e+09 5.34e-10 1.80e+04 -1.0 8.53e+04 - 2.80e-01 5.67e-01f 1
29 -1.5965611e+09 5.25e-10 1.46e+04 -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.27e-10 1.13e+04 -1.0 9.18e+04 - 2.21e-01 3.60e-01f 1
31 -1.5965878e+09 5.91e-10 8.75e+03 -1.0 9.17e+04 - 2.29e-01 3.99e-01f 1
32 -1.5965937e+09 5.47e-10 6.56e+03 -1.0 9.06e+04 - 2.50e-01 4.54e-01f 1
33 -1.5965974e+09 4.95e-10 4.78e+03 -1.0 7.97e+04 - 2.71e-01 5.12e-01f 1
34 -1.5966027e+09 4.91e-10 3.42e+03 -1.0 6.07e+04 - 2.84e-01 5.30e-01f 1
35 -1.5966013e+09 5.30e-10 2.41e+03 -1.0 8.43e+04 - 2.95e-01 6.23e-01f 1
36 -1.5965947e+09 5.78e-10 1.61e+03 -1.0 6.27e+04 - 3.32e-01 7.29e-01f 1
37 -1.5965872e+09 5.77e-10 1.01e+03 -1.0 6.84e+04 - 3.75e-01 9.34e-01f 1
38 -1.5965795e+09 5.52e-10 5.39e+02 -1.0 3.69e+04 - 4.64e-01 1.00e+00f 1
39 -1.5965702e+09 6.13e-10 1.83e+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 5.71e-10 5.90e+01 -1.0 2.92e+04 - 9.90e-01 1.00e+00f 1
41 -1.5965644e+09 5.93e-10 2.28e+00 -1.0 7.96e+03 - 1.00e+00 1.00e+00f 1
42 -1.5967508e+09 5.55e-10 8.78e+05 -2.5 1.54e+05 - 8.91e-01 9.82e-01f 1
43 -1.5967560e+09 5.64e-10 2.11e+00 -2.5 6.73e+03 - 1.00e+00 1.00e+00f 1
44 -1.5967560e+09 5.99e-10 4.61e-01 -2.5 3.47e+01 - 1.00e+00 1.00e+00f 1
45 -1.5967560e+09 5.30e-10 1.09e-02 -2.5 5.35e+00 - 1.00e+00 1.00e+00f 1
46 -1.5967588e+09 5.16e-10 1.19e+05 -3.8 4.45e+03 - 1.00e+00 5.56e-01f 1
47 -1.5967599e+09 5.34e-10 5.36e+04 -3.8 1.97e+03 - 1.00e+00 5.49e-01f 1
48 -1.5967604e+09 5.05e-10 1.91e+04 -3.8 8.72e+02 - 9.22e-01 6.43e-01f 1
49 -1.5967605e+09 5.54e-10 1.18e+03 -3.8 2.83e+02 - 1.00e+00 9.38e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 -1.5967605e+09 6.23e-10 1.09e-01 -3.8 5.05e+01 - 1.00e+00 1.00e+00f 1
51 -1.5967605e+09 6.44e-10 5.86e-03 -3.8 1.33e+01 - 1.00e+00 1.00e+00f 1
52 -1.5967605e+09 5.46e-10 2.15e-05 -3.8 8.11e-01 - 1.00e+00 1.00e+00h 1
53 -1.5967606e+09 5.37e-10 6.82e+03 -5.7 1.54e+02 - 8.46e-01 3.87e-01f 1
54 -1.5967606e+09 5.33e-10 4.79e+03 -5.7 4.71e+01 - 1.92e-01 2.71e-01f 1
55 -1.5967606e+09 4.94e-10 3.23e+03 -5.7 2.95e+01 - 3.27e-01 3.26e-01f 1
56 -1.5967606e+09 7.18e-10 3.53e+03 -5.7 1.90e+01 - 6.07e-01 1.02e-01f 1
57 -1.5967606e+09 5.37e-10 2.36e+03 -5.7 1.38e+01 - 7.60e-02 2.99e-01f 1
58 -1.5967606e+09 7.03e-10 1.40e+03 -5.7 9.27e+00 - 2.21e-01 3.76e-01f 1
59 -1.5967606e+09 5.23e-10 8.52e+02 -5.7 5.52e+00 - 4.97e-01 4.14e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
60 -1.5967606e+09 5.57e-10 5.14e+02 -5.7 3.07e+00 - 8.68e-01 4.77e-01f 1
61 -1.5967606e+09 5.48e-10 1.33e+02 -5.7 1.54e+00 - 1.00e+00 7.53e-01f 1
62 -1.5967606e+09 6.23e-10 1.71e-05 -5.7 3.76e-01 - 1.00e+00 1.00e+00f 1
63 -1.5967606e+09 6.89e-10 1.01e+02 -7.0 5.79e-01 - 9.95e-01 4.19e-01f 1
64 -1.5967606e+09 5.30e-10 5.01e-01 -7.0 3.32e-01 - 1.00e+00 9.95e-01f 1
65 -1.5967606e+09 5.49e-10 3.49e-09 -7.0 1.61e-03 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 65
(scaled) (unscaled)
Objective...............: -2.6787832953601621e+07 -1.5967606334020669e+09
Dual infeasibility......: 3.4906158039623750e-09 2.0806751750887140e-07
Constraint violation....: 1.1641532182693481e-10 5.4893600776040330e-10
Complementarity.........: 9.1160146931116271e-08 5.4338450671570302e-06
Overall NLP error.......: 9.1160146931116271e-08 5.4338450671570302e-06
Number of objective function evaluations = 66
Number of objective gradient evaluations = 66
Number of equality constraint evaluations = 66
Number of inequality constraint evaluations = 66
Number of equality constraint Jacobian evaluations = 66
Number of inequality constraint Jacobian evaluations = 66
Number of Lagrangian Hessian evaluations = 65
Total CPU secs in IPOPT (w/o function evaluations) = 49.355
Total CPU secs in NLP function evaluations = 21.167
EXIT: Optimal Solution Found.
[11]:
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 212160, 'Number of variables': 211203, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 72.74183440208435}], '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.
[12]:
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.013e-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.013e-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.013e-3 * 3600
for t in set_time]
# Plot the results
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)})
[13]:
# 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=m.LMP)
Optimal PEM capacity : 196.2616021908122 MW
Optimal tank capacity : 2.0756993646370207e-05 kg
OPtimal turbine capacity: 9.076695580817415e-05 MW
[14]:
# 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=m.LMP)
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.
[15]:
# 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).
[16]:
# 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
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\).
[17]:
# Create concrete model object
m = ConcreteModel()
# Define sets
m.set_time = RangeSet(24) # twenty fours in a day
m.set_days = RangeSet(20) # twenty clusters/days per year
m.set_years = [2022] # Set of years
m.set_scenarios = [0, 1] # Set of scenarios
# 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, we call the build_multiperiod_design function to formulate the multi-period optimization problem. However, this time, we also set stochastic=True.
[18]:
# Build initialized multiperiod model: Constraints g(u(t,d)) = 0 and h(u(t,d)) <= 0
build_multiperiod_design(m,
flowsheet=build_ne_flowsheet,
initialization=fix_dof_and_initialize,
unfix_dof=unfix_dof,
multiple_days=True,
multiyear=True,
stochastic=True,
verbose=False)
[ 0.00] Processing input information.
[+ 0.00] Beginning the formulation of the multiperiod problem.
[+ 253.45] Completed the formulation of the multiperiod problem
2022-07-01 17:18:34 [INFO] idaes.init.fs.pem.outlet_state: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.pem.outlet_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.pem.outlet_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_in: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_turbine: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank.properties_out_pipeline: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.h2_tank: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator.properties_in: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator.properties_out: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.translator: Initialization Complete optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.mixer.air_feed_state: Starting initialization
2022-07-01 17:18:34 [INFO] idaes.init.fs.mixer.air_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:34 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.mixer.hydrogen_feed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.mixer.mixed_state: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.mixer.mixed_state: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.mixer.mixed_state: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.mixer: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine: Starting initialization...
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.compressor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Starting initialization
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:35 [INFO] idaes.init.fs.h2_turbine.stoic_reactor.control_volume: Initialization Complete
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.stoic_reactor: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Starting initialization
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_in: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Starting initialization
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.control_volume.properties_out: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Starting initialization
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine.properties_isentropic: Property package initialization: optimal - Optimal Solution Found.
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine.turbine: Initialization Complete: optimal - Optimal Solution Found
2022-07-01 17:18:36 [INFO] idaes.init.fs.h2_turbine: Initialization complete
[+ 2.91] Created an instance of the flowsheet and initialized it.
[+ 24.67] Initialized the entire multiperiod optimization model.
[+ 0.43] Unfixed the degrees of freedom from each period model.
Next, add the connecting constraints and hydrogen demand constraints for each of the scenarios.
[19]:
for s in m.set_scenarios:
# Build the connecting constraints
build_connecting_constraints(m.scenario[s],
set_time=m.set_time,
set_days=m.set_days,
set_years=m.set_years)
# Append cash flow expressions
append_costs_and_revenue(m.scenario[s],
ps=m,
LMP=m.LMP[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 / 2.016e-3
[20]:
# Add non-anticipativity constraints
m.pem_capacity = Var(within=NonNegativeReals,
doc="Design PEM capacity (in kW)")
m.tank_capacity = Var(within=NonNegativeReals,
doc="Design tank capacity (in mol)")
m.h2_turbine_capacity = Var(within=NonNegativeReals,
doc="Design turbine capacity (in 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
[21]:
# 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.13.2: 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 http://projects.coin-or.org/Ipopt
This version of Ipopt was compiled from source code available at
https://github.com/IDAES/Ipopt as part of the Institute for the Design of
Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.
This version of Ipopt was compiled using HSL, a collection of Fortran codes
for large-scale scientific computation. All technical papers, sales and
publicity material resulting from use of the HSL codes within IPOPT must
contain the following acknowledgement:
HSL, a collection of Fortran codes for large-scale scientific
computation. See http://www.hsl.rl.ac.uk.
******************************************************************************
This is Ipopt version 3.13.2, running with linear solver ma27.
Number of nonzeros in equality constraint Jacobian...: 580772
Number of nonzeros in inequality constraint Jacobian.: 7680
Number of nonzeros in Lagrangian Hessian.............: 195840
Total number of variables............................: 211209
variables with only lower bounds: 20169
variables with lower and upper bounds: 176640
variables with only upper bounds: 0
Total number of equality constraints.................: 208326
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 4.40e+03 -1.0 1.01e+04 - 4.84e-01 9.71e-01f 1
4 -1.5737296e+09 2.21e+02 1.17e+04 -1.0 2.04e+02 - 2.69e-02 5.23e-02f 4
5 -1.5737373e+09 1.69e+02 5.96e+04 -1.0 1.39e+02 - 7.09e-02 2.32e-01h 2
6 -1.5737516e+09 9.33e+01 1.20e+05 -1.0 9.11e+01 - 8.23e-01 4.43e-01h 1
7 -1.5738919e+09 9.07e+00 3.52e+05 -1.0 5.28e+02 - 4.10e-01 9.39e-01f 1
8 -1.5741711e+09 2.33e+01 3.18e+06 -1.0 1.02e+03 - 3.86e-01 1.00e+00f 1
9 -1.5741718e+09 2.03e+01 8.01e+05 -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 1.76e+06 -1.0 2.88e+02 -4.5 1.95e-01 5.00e-01f 2
11 -1.5742405e+09 1.20e+00 2.98e+06 -1.0 4.64e+02 -5.0 1.67e-01 1.00e+00f 1
12 -1.5758575e+09 1.61e+00 2.91e+06 -1.0 1.25e+04 - 2.49e-02 1.00e+00f 1
13 -1.5779187e+09 3.22e-02 2.58e+06 -1.0 1.69e+04 - 1.13e-01 6.84e-01f 1
14 -1.5835775e+09 1.34e-02 2.32e+06 -1.0 2.97e+04 - 1.00e-01 9.54e-01f 1
15 -1.5978308e+09 4.90e-03 2.08e+06 -1.0 6.24e+04 - 1.06e-01 1.00e+00f 1
16 -1.6303589e+09 2.75e-04 1.87e+06 -1.0 1.53e+05 - 9.86e-02 8.53e-01f 1
17 -1.6306844e+09 2.73e-04 1.61e+06 -1.0 2.14e+05 - 1.40e-01 6.28e-03f 1
18 -1.6306914e+09 2.64e-04 9.51e+05 -1.0 4.69e+04 - 4.10e-01 3.19e-02f 1
19 -1.6307897e+09 3.96e-05 6.54e+05 -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 3.68e+05 -1.0 6.58e+04 - 4.37e-01 8.21e-01f 1
21 -1.6309292e+09 2.98e-06 1.64e+05 -1.0 6.00e+04 - 5.55e-01 1.00e+00f 1
22 -1.6309356e+09 1.34e-06 1.32e+05 -1.0 5.99e+04 - 1.96e-01 4.56e-01f 1
23 -1.6309409e+09 7.57e-07 9.41e+04 -1.0 5.16e+04 - 2.85e-01 4.97e-01f 1
24 -1.6309464e+09 2.32e-07 6.55e+04 -1.0 5.17e+04 - 3.04e-01 8.46e-01f 1
25 -1.6309419e+09 4.16e-08 4.27e+04 -1.0 4.01e+04 - 3.49e-01 9.11e-01f 1
26 -1.6309368e+09 1.38e-08 2.71e+04 -1.0 3.75e+04 - 3.65e-01 1.00e+00f 1
27 -1.6309299e+09 5.02e-09 1.61e+04 -1.0 3.24e+04 - 4.05e-01 1.00e+00f 1
28 -1.6309221e+09 1.48e-09 8.34e+03 -1.0 1.99e+04 - 4.82e-01 1.00e+00f 1
29 -1.6309164e+09 6.35e-10 3.43e+03 -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 6.09e-10 3.41e+01 -1.0 1.17e+04 - 9.90e-01 1.00e+00f 1
31 -1.6309083e+09 5.93e-10 8.19e-01 -1.0 4.00e+03 - 9.97e-01 1.00e+00f 1
32 -1.6310468e+09 6.13e-10 6.61e+05 -2.5 1.13e+05 - 9.25e-01 9.93e-01f 1
33 -1.6310486e+09 6.36e-10 4.06e+00 -2.5 2.83e+03 - 1.00e+00 1.00e+00f 1
34 -1.6310486e+09 5.88e-10 4.79e-01 -2.5 3.44e+01 - 1.00e+00 1.00e+00f 1
35 -1.6310486e+09 5.56e-10 1.13e-02 -2.5 5.29e+00 - 1.00e+00 1.00e+00f 1
36 -1.6310507e+09 5.19e-10 1.17e+05 -3.8 3.21e+03 - 1.00e+00 5.64e-01f 1
37 -1.6310513e+09 4.69e-10 5.31e+04 -3.8 1.39e+03 - 1.00e+00 5.46e-01f 1
38 -1.6310516e+09 5.42e-10 1.89e+04 -3.8 6.08e+02 - 1.00e+00 6.43e-01f 1
39 -1.6310517e+09 5.32e-10 2.78e+03 -3.8 1.72e+02 - 1.00e+00 8.53e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40 -1.6310517e+09 5.45e-10 5.47e-02 -3.8 2.77e+01 - 1.00e+00 1.00e+00f 1
41 -1.6310517e+09 4.97e-10 1.40e-03 -3.8 6.39e+00 - 1.00e+00 1.00e+00f 1
42 -1.6310517e+09 5.09e-10 3.08e+03 -5.7 8.60e+01 - 5.85e-01 3.77e-01f 1
43 -1.6310518e+09 5.31e-10 2.08e+03 -5.7 3.12e+01 - 3.37e-01 3.32e-01f 1
44 -1.6310518e+09 6.33e-10 3.28e+03 -5.7 1.98e+01 - 4.51e-01 1.05e-01f 1
45 -1.6310518e+09 5.33e-10 1.99e+03 -5.7 1.53e+01 - 9.80e-02 2.73e-01f 1
46 -1.6310518e+09 4.98e-10 7.84e+02 -5.7 9.71e+00 - 2.71e-01 4.37e-01f 1
47 -1.6310518e+09 5.06e-10 7.81e+02 -5.7 5.01e+00 - 5.83e-01 3.82e-01f 1
48 -1.6310518e+09 4.89e-10 6.89e+02 -5.7 2.88e+00 - 9.20e-01 4.72e-01f 1
49 -1.6310518e+09 5.72e-10 2.99e+02 -5.7 1.45e+00 - 1.00e+00 5.95e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50 -1.6310518e+09 5.27e-10 6.73e+01 -5.7 5.87e-01 - 1.00e+00 7.75e-01f 1
51 -1.6310518e+09 5.96e-10 2.96e-06 -5.7 1.32e-01 - 1.00e+00 1.00e+00f 1
52 -1.6310518e+09 5.05e-10 6.64e+01 -7.0 4.10e-01 - 1.00e+00 6.22e-01f 1
53 -1.6310518e+09 6.55e-10 3.09e-06 -7.0 1.53e-01 - 1.00e+00 1.00e+00f 1
54 -1.6310518e+09 5.51e-10 2.57e-12 -7.0 2.88e-04 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 54
(scaled) (unscaled)
Objective...............: -3.7401907403570741e+07 -1.6310517839609087e+09
Dual infeasibility......: 2.5730085157677933e-12 1.1220577829110959e-10
Constraint violation....: 1.1641532182693481e-10 5.5087978623191702e-10
Complementarity.........: 9.0916704111854490e-08 3.9647671129020865e-06
Overall NLP error.......: 9.0916704111854490e-08 3.9647671129020865e-06
Number of objective function evaluations = 63
Number of objective gradient evaluations = 55
Number of equality constraint evaluations = 63
Number of inequality constraint evaluations = 63
Number of equality constraint Jacobian evaluations = 55
Number of inequality constraint Jacobian evaluations = 55
Number of Lagrangian Hessian evaluations = 54
Total CPU secs in IPOPT (w/o function evaluations) = 45.865
Total CPU secs in NLP function evaluations = 18.631
EXIT: Optimal Solution Found.
[21]:
{'Problem': [{'Lower bound': -inf, 'Upper bound': inf, 'Number of objectives': 1, 'Number of constraints': 212166, 'Number of variables': 211209, 'Sense': 'unknown'}], 'Solver': [{'Status': 'ok', 'Message': 'Ipopt 3.13.2\\x3a Optimal Solution Found', 'Termination condition': 'optimal', 'Id': 0, 'Error rc': 0, 'Time': 67.20030307769775}], 'Solution': [OrderedDict([('number of solutions', 0), ('number of solutions displayed', 0)])]}
[22]:
# 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.417908572877676e-05 kg
OPtimal turbine capacity: 9.076374619487129e-05 MW
[23]:
generate_plots(m.scenario[0], d=8, y=2022, set_time=m.set_time, lmp=m.LMP[0])
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.