15 Design before demand is known
Problem 15 | Hydrogen Production and Storage | Level 3: Open-ended challenge
15.1 Problem Statement
An electrolyzer operates over six four-hour periods, representing one day. Electricity use is constant at 0.05 MWh/kg of hydrogen. Power is limited to 4 MW. Initial hydrogen inventory is 50 kg, and final inventory must also be 50 kg. Storage has no leakage and no charging-rate restriction. Production and withdrawal rates are each uniform within a period. Inventory therefore varies linearly between period boundaries.
| Period | Electricity price (USD/MWh) | Hydrogen demand (kg/period) |
|---|---|---|
| 1 | 30 | 120 |
| 2 | 25 | 150 |
| 3 | 80 | 180 |
| 4 | 100 | 220 |
| 5 | 45 | 180 |
| 6 | 35 | 140 |
Capacity is enforced at every period boundary, including the 50 kg initial stock. Under the uniform-rate assumption, this also bounds inventory throughout each period. If deliveries are discrete or rates vary within a period, an additional inventory profile or shorter time steps is necessary for physical sizing. A daily capacity charge of 0.15 USD/(kg d) is used where investment is optimized. It is an assumed daily equivalent, not a full capital-cost model.
Your task. Keep one shared tank capacity. Daily demand is either 0.9 or 1.15 times the table, each with probability 0.5. The demand state is fully observed before the first dispatch decision. In either state, the electrolyzer is initially off, consumes at least 1 MW when on, and incurs 10 USD per startup. Minimize capacity cost plus expected electricity and startup cost.
15.1.1 Before you code
Why may dispatch depend on the state while tank capacity may not? Would the same model be valid if the demand state were revealed only halfway through the day?
15.1.2 Deliverables
- Define decision variables, units, objective, and constraints. State which data and assumptions are fixed.
- Predict one feature of the solution and explain why it should occur.
- Build and solve a Pyomo model. Check balances and bounds independently.
- Interpret the result and investigate the follow-up question below. For an open-ended challenge, state and defend your chosen policy separately from the reference solution.
15.2 Model Formulation
Let \(P_{st}\) be power in MW, \(I_{st}\) end inventory in kg, and \(K\) capacity in kg. With \(\Delta t=4\) h and \(e=0.05\) MWh/kg, \[I_{st}=I_{s,t-1}+\frac{\Delta t}{e}P_{st}-D_{st},\qquad 0\le P_{st}\le4,\quad0\le I_{st}\le K.\] Use \(I_{s0}=50\) before the first period and \(I_{s6}=50\) at the end. Total production must equal total demand. The electricity cost is \(\sum_t c_t\Delta t P_{st}\) USD/d. Omitting \(\Delta t\) would understate cost by a factor of four.
Use state-specific binaries \(y_{st}\) for operation and \(b_{st}\) for startup: \[y_{st}\le P_{st}\le4y_{st},\qquad b_{st}\ge y_{st}-y_{s,t-1},\quad y_{s0}=0.\] The minimum-load coefficient is 1 MW. Minimize \[0.15K+\sum_s\pi_s\sum_t\left(c_t\Delta tP_{st}+10b_{st}\right).\] The positive startup cost makes unnecessary startups unattractive. Capacity is shared across the two states; dispatch is recourse after full daily demand revelation. There is no minimum up/down time, ramp limit, or commitment carried to the next day. The terminal stock is restored separately in both states.
15.3 Pyomo Implementation
The code below is the model-building portion of models/hydrogen.py. The level argument selects this problem’s assumptions. Run the complete module from the source bundle to reproduce the solution, including its independent checks:
~/.venvs/optim/bin/python -m models.hydrogen 3"""Daily hydrogen dispatch, tank design, and scenario-dependent operation."""
import pyomo.environ as pyo
from models.common import value
PRICE = [30, 25, 80, 100, 45, 35] # USD/MWh
DEMAND = [120, 150, 180, 220, 180, 140] # kg per four-hour period
DT = 4 # h
SPECIFIC_ENERGY = 0.05 # MWh/kg
TANK_CHARGE = 0.15 # USD/(kg capacity d)
def build(level):
states = (
{"nominal": (1, 1)} if level < 3 else {"low": (0.9, 0.5), "high": (1.15, 0.5)}
)
m = pyo.ConcreteModel(name="Hydrogen production and storage")
m.S = pyo.Set(initialize=list(states))
m.T = pyo.RangeSet(0, 5)
m.k = pyo.Var(bounds=(50, 500))
if level == 1:
m.k.fix(150)
m.power = pyo.Var(m.S, m.T, bounds=(0, 4))
m.stock = pyo.Var(m.S, m.T, domain=pyo.NonNegativeReals)
m.balance = pyo.Constraint(
m.S,
m.T,
rule=lambda m, s, t: (
m.stock[s, t]
== (50 if t == 0 else m.stock[s, t - 1])
+ DT * m.power[s, t] / SPECIFIC_ENERGY
- DEMAND[t] * states[s][0]
),
)
m.capacity = pyo.Constraint(m.S, m.T, rule=lambda m, s, t: m.stock[s, t] <= m.k)
m.terminal = pyo.Constraint(m.S, rule=lambda m, s: m.stock[s, 5] == 50)
if level == 3:
m.on = pyo.Var(m.S, m.T, domain=pyo.Binary)
m.start = pyo.Var(m.S, m.T, domain=pyo.Binary)
m.upper = pyo.Constraint(
m.S, m.T, rule=lambda m, s, t: m.power[s, t] <= 4 * m.on[s, t]
)
m.lower = pyo.Constraint(
m.S, m.T, rule=lambda m, s, t: m.power[s, t] >= m.on[s, t]
)
m.startup = pyo.Constraint(
m.S,
m.T,
rule=lambda m, s, t: (
m.start[s, t] >= m.on[s, t] - (0 if t == 0 else m.on[s, t - 1])
),
)
fixed = 0 if level == 1 else TANK_CHARGE * m.k
m.obj = pyo.Objective(
expr=fixed
+ sum(
states[s][1]
* sum(
DT * PRICE[t] * m.power[s, t]
+ (10 * m.start[s, t] if level == 3 else 0)
for t in m.T
)
for s in m.S
)
)
m._states = states
return mThe common solve function calls appsi_highs, checks optimal termination before loading a solution, and checks constraint residuals, bounds, and integer domains. After building the model, use:
from models.common import solve
m = solve(build(3))15.4 Optimal Solution
The reference objective is 1,785.0500 USD/d. This value is optimal for the explicitly stated model and data.
| State | Period | Power (MW) | Demand (kg) | End stock (kg) |
|---|---|---|---|---|
| low | 1 | 3.9375 | 108.0000 | 257.0000 |
| low | 2 | 4.0000 | 135.0000 | 442.0000 |
| low | 3 | 0.0000 | 162.0000 | 280.0000 |
| low | 4 | -0.0000 | 198.0000 | 82.0000 |
| low | 5 | 1.0000 | 162.0000 | 0.0000 |
| low | 6 | 2.2000 | 126.0000 | 50.0000 |
| high | 1 | 4.0000 | 138.0000 | 232.0000 |
| high | 2 | 4.0000 | 172.5000 | 379.5000 |
| high | 3 | 1.0062 | 207.0000 | 253.0000 |
| high | 4 | 0.0000 | 253.0000 | 0.0000 |
| high | 5 | 2.5875 | 207.0000 | 0.0000 |
| high | 6 | 2.6375 | 161.0000 | 50.0000 |
| Quantity | Value |
|---|---|
| Tank capacity (kg) | 442.0000 |
| Expected electricity (MWh/d) | 50.7375 |
15.4.1 Check your solution
Check power is either zero or between 1 and 4 MW. Check mass balance and terminal stock in each state. Count startups directly from the on/off sequence.
The automated residual, bound, and integrality audit passed with a maximum violation of 2.84e-14 in model units. Domain-specific checks passed. Model units differ across equations, so this numerical audit does not replace dimensional analysis.
15.5 Brief Discussion
The shared tank is 442.0 kg and expected daily cost is 1785.05 USD. Expected demand is 1.025 times the original profile, so the change from the deterministic case combines a higher mean demand, uncertainty, and commitment costs. It is not a pure value-of-uncertainty comparison.
15.5.1 Extension to investigate
Formulate the nonanticipativity constraints required if the state is learned after period 2. State which decisions must coincide before that time, and explain how uncertainty in earlier withdrawals affects the scenario definition.