Dynamic Optimization of Proton Exchange Membrane Water Electrolzyers Considering Usage-Based Degradation (2024)

2.3.1 Outer Problem

The outer problem iterates over the total number of cells in the electrolyzer set up, $N_{c}$, and number of days of \chH2 storage $\Gamma_{\ch{H2},stor}$. Each cell is assumed to have an area of 450 $\frac{{\mathrm{cm}}^{2}}{}$ represented by $A_{cell}$. The cost of the stack and storage can then be calculated as

where $C^{mBoP}$ mechanical balance of plant capital cost, $C^{eBoP}$, electrical balance of plant capital cost, $p^{stack}$ is the capital cost of the stack in $/$\frac{{\mathrm{cm}}^{2}}{}$, $\dot{m}_{\ch{H2}}$ is the daily \chH2 production rate in kilograms as calculated by the inner problem, and $p^{storage}$ is the cost of storage in $/$\frac{\mathrm{kg}}{}$ \chH2.

The mechanical and electrical balance of plant are scaled based on the rate of \chH2 production and peak power consumption for the facility, as shown in Equations 4 and 5. The mechanical balance of plant represents the compressors, pumps, DI water system, and \chH2 separation system while the electrical balance of plant includes all electrical components necessary to make the interconnection to the grid including a rectifier. Here, $\alpha^{mBoP}$ is a cost factor in $/$\frac{\mathrm{kg}}{}$ \chH2, $\alpha^{eBoP}$ is a cost factor in $/$\frac{\mathrm{kW}}{}$e, and $P^{max}$ is peak power consumption of the system in $\frac{\mathrm{kW}}{}$e as calculated by the inner problem. $P_{max}$ is then a parameter passed to the outer problem to calculate the CAPEX for a given outer problem iteration.

2.3.2 Inner Problem

The inner problem takes a fixed $N_{c}$ and $\Gamma_{\ch{H2},stor}$ and minimizes operating cost while meeting the daily production target of $50,000$\frac{\mathrm{kg}}{\mathrm{d}}$$. The inner problem can be written as a minimization of the operating cost as defined in Equation 6.

where x represents the primary decision variables, $C^{elec}$ is the annual cost of electricity used by the electrolyzer, $C^{BoP,elec}$ is the annual cost of electricity consumed by the balance of plant, $C^{water}$ is the annual cost of the deionized water, and $C^{\ch{N2}}$ is the annual cost of the nitrogen gas used at the anode for purging.

The model is subject to the mass and energy balances outlined in the SI as well as a demand constraint of 50,000 $\frac{\mathrm{kg}}{\mathrm{d}}$ of \chH2. An operating temperature range of $60^{\circ}$C - $80^{\circ}$C was imposed on the stack. Additionally, a minimum and maximum current density were applied in this study. The minimum current density was chosen based on minimizing degradation at low voltages and the high current density was chosen based on the limit of commercialized PEM today [Weis2019-vj]. In order to ensure 4% \chH2 in \chO2 at the anode is never reached a safety factor of 2 was applied making another constraint of a maximum of 2% \chH2 in \chO2 at the anode.

Evaluating each variable in the system for 15 minute time periods over the course of a year quickly becomes computationally intractable. To maintain computational tractability, we approximate annual operating costs through modeling operation over representative days at 15 min resolution assuming the hourly price of electricity holds for the entire hour. Time domain reduction using k-means clustering was performed for each electricity price time series data in order to find 7 representative days. Further details of the clustering method can be found in SI 5.2.

The clustering created a mapping $f:d\rightarrow r$. Each real day, $d$, is mapped to a representative day, $r$. To distinguish between variables for real or representative days, any variable with a bar above it is a variable that is indexed by the set of representative days.

Beginning with the cost of electricity for the electrolyzer,

where $\tau$ represents the end of a representative day, $p^{elec}$ represents the price of electricity, $i$ is the current density, $A=N_{c}A_{cell}$ is the total electrolyzer area, $V^{undeg}$ is the voltage without considering degradation (based on polarization curve - see Equations S1 -S12), and $V^{deg,cuml}$ is the cumulative increase in voltage due to degradation.

The first term in the integral can be mapped to its respective representative day via $f(d)\rightarrow r$. However, the second term that accounts for the cumulative degradation that changes over the 365 day range and is inter-temporally coupled. So $C^{elec}_{d}$ becomes:

In accordance with the degradation correlation based on experiential literature values seen in Figure 4, we use the following model to represent the change in degradation over a representative day $r$:

Cumulative degradation for real day $d$ can then be written as

where $\delta\bar{V}_{deg,f(d)}(t)$ comes from solving the differential equation in Equation 10 up to the time point of interest. Initial degradation at the beginning of the year for a fresh stack is assumed to be 0 $\frac{\mathrm{V}}{}$.

Using the mapping $f(d)\rightarrow r$, $C^{elec}_{d}$ can be written

The annual cost of electricity is then written:

where $w_{r}$ represents the weight of representative day $r$ as determined by the clustering.

$C^{BoP,elec}$, $C^{water}$, and $C^{\ch{N2}}$ can be calculated in a similar manner using the associated representative days and mapping $f$.

where $W^{compress}$ is the work of the \chH2 storage compressor, $p^{water}$ is the price of deionized water, and $p^{\ch{N2}}$ is the price of nitrogen. Compressor work is evaluated using a single stage isentropic compression with an isentropic factor of 0.7 as outlined in Kahn et. al and Chung et. al [Chung2024-vr, Kahn2021-jv].

2.3.3 H₂ Storage

Long duration storage as well as intra-day storage are formulated as follows:

where $\delta_{r}$ is the positive or negative change in storage levels that occurs over representative day $r$ and allows for \chH2 storage to carry over between representative periods. Equation 18 represents a wrapping constraint between the last day in the year and the first day in the year. Equation 19 allows \chH2 storage to carry over between real days. Equation 20 relates the beginning of a real day of storage to the associated beginning of the representative day. Equations 21 - 22 ensure the mass balance across the \chH2 storage is closed over real and representative days. Finally, Equation 23 ensures all storage state of charge variables are less than the maximum amount of \chH2 storage. A visual representation of the storage formulation is included in the SI in Figure S1.

2.3.4 Implementation Approach

The outer and inner problems can be combined in an iterative manner using the algorithm described in Figure 2.

Upper and lower bounds for the electrolyzer area and number of days of \chH2 storage are chosen based on physical limits of the electrolyzer in order to satisfy the 50,000 $\frac{\mathrm{kg}}{\mathrm{d}}$ \chH2 demand. For this study, an area range of $40,000-300,000$ 450 $\frac{{\mathrm{cm}}^{2}}{}$ cells and a storage range of $0.1-14$ days were chosen where one day of storage is equivalent to 50,000 $\frac{\mathrm{kg}}{}$ of \chH2. From the GSS algorithm two trial areas and two trial storages are chosen as is done by Sandhya et al assuming unimodality in each of the search directions [Sandhya-Rani2019-pw]. This creates a total of four trial points in this 2D search space. The operational optimization is run for those 4 trial points which results in four different variable operating costs. For the PEM system, in the limit of very large electrolyzer areas, high capital costs would result in a high present value. In the limit of small areas, the system would need to operate at high current densities which would result in more stack replacement and a high present value. Storage also exhibits a similar unimodal trend at the extremes.

With the size of each system and the variable operating cost estimated, the present value for each of the four trial points can be calculated. If the trial areas and number of storage days are within the user defined tolerance (0.1% in this study), then the algorithm terminates and the final solution corresponds to the solution of trial point problems with the lowest present value. The variable operating routine from the associated minimum present value is then the optimal operation routine and the area and number of days of storage from the minimum present value are the optimal sizing of the system.

If still outside the tolerance, regions in the 2D space where the minimum cannot lie are eliminated. For example, for $A$ and $\Gamma_{\ch{H2},stor}$ where $A_{lb}\leq A\leq A_{ub}$ and $\Gamma_{\ch{H2},stor,lb}\leq\Gamma_{\ch{H2},stor}\leq\Gamma_{\ch{H2},stor,ub}$ and given four present values that are functions of the two trial areas and two trial storages $PV_{A}(A_{1},\Gamma_{\ch{H2},stor,1}),PV_{B}(A_{2},\Gamma_{\ch{H2},stor,1}),PV$, we can eliminate a portion of the search region by finding the minimum of the present values. If $PV_{A}$ is the smallest present value, we know the minimum will not lie in $A_{2}\leq A\leq A_{ub}$ and $\Gamma_{\ch{H2},stor,2}\leq\Gamma_{\ch{H2},stor}\leq\Gamma_{\ch{H2},stor,ub}$. Thus $A_{2}$ and $\Gamma_{\ch{H2},stor,2}$ become the new upper bounds. Similar eliminations can be made if the smallest present value is one of the other three values. The algorithm repeats until converged on a locally optimal solution. Due to the non-linearity of the problem, a global optimum is not guaranteed.

The NLP to solve the operational problem was implemented in Python using Pyomo, specifically Pyomo DAE with IPOPT as the non-linear solver [Nicholson2018-mw, Wachter2006-la]. Linear solver MA86 was chosen to use in IPOPT for its ability to handle large problems well [UnknownUnknown-hg]. The bilevel optimization approach with the degradation correlation results in a problem with $\approx 200,000$ continuous variables, $\approx 200,000$ equality constraints, and $\approx 2,500$ inequality constraints. Full details of each run can be found in Table 2 and in the SI in Table S4. The problem was run using resources on MIT’s SuperCloud HPC [Reuther2018-xn].

Table S5 shows how the bi-level optimization approach scales with the number of representative days. As we increase the number of representative days we expect convergence to a value for storage, number of cells, and LCOH, although this will not be monotonic due to the different effect of each added day to the optimization. Although increasing number of representative days allows for capturing greater day-to-day variability in the underlying electricity price time series, it increases the overall problem size and the computational burden as indicated by the time to convergence and avg. run time per GSS iteration. We choose 7 representative days to strike a balance between computational tractability and accuracy.

2.4.1 Techno-economic Data

For the 2022 case, capital and operating cost assumptions were made consistent with the 2019 version of NREL’s H2A model adjusting for inflation from 2019 to 2022 [Brian_D_James_Daniel_A_DeSantis_Genevieve_Saur2016-wz]. The 2030 case capital costs come from Fraunhofer ISE’s cost forecast for PEM and alkaline electrolysis [Holst2021-nd]. Summary of the cost assumptions can be found in Table 1.

Using the parameters in Table 1, capital cost ($C^{CAPEX}=C^{stack}+C^{storage}$), unplanned replacement costs ($C^{UPR}$), planned replacement costs ($C^{PR}$), fixed operating cost ($C^{fOPEX}$), and variable operating cost ($C^{vOPEX}$) can be estimated. Assuming a plant life of 40 years and a discount rate of 8%, the present value of the project ($PV$) can be estimated as shown in Equation 24.

The direct capital needed to calculate the above costs is the sum of the cost of the stack and cost of the balance of plant. For the scenario run without useage-dependent stack degradation, the replacement rate for the stack was fixed at 7 years [Brian_D_James_Daniel_A_DeSantis_Genevieve_Saur2016-wz]. For the scenarios where degradation is accounted for as shown in Equation 10, a degradation rate-informed replacement rate is calculated as follows. Based on the fixed H2A degradation rate, approximately $\Delta^{max}_{deg}\approx 1$\frac{\mathrm{V}}{}$$ of degradation is allowed to occur before the stack needs to be replaced [Brian_D_James_Daniel_A_DeSantis_Genevieve_Saur2016-wz]. The operational optimization of the problem gives the degradation of a fresh stack after one year of operation ($\Delta^{1}_{deg}$). The replacement rate can then be defined as:

For implementation purposes, if the degradation of the stack exceeds 1 V in a given year, the stack is allowed to run until the end of the year before being replaced.

Fixed operating cost is the sum of the cost of labor, overhead, tax and insurance, and material. The labor rate was assumed to be $70 per hour per worker where 10 workers are needed for a system of this size, and the plant is assumed to operate 350 days per year.

Variable operating cost for the first year is simply obtained from the objective function of the operational optimization. Post-convergence of the operational optimization, the variable operating cost is updated to include degradation for the years leading up to replacement. This assumes that degradation in each future year is the same as the first year after stack replacement.

The amount of \chH2 produced in one year ($M_{\ch{H2}}$) can be discounted in a similar way as the present value of the project as seen in Equation 26. The levelized cost of \chH2 is then calculated as the ratio of the two present values.

2.4.2 Electrolyzer Performance and Model Assumptions

Validation of the electrochemical model without degradation at different temperatures and operating pressures is shown in Figure 3. The model fits well across varying current densities, pressures, and temperatures.

Many degradation studies have been done on PEM cells that attempt to characterize the degradation rate of the cell as a function of the applied current density [Suermann2019-kz, Rakousky2017-te, Papakonstantinou2020-hj, Alia2019-lw, Li2021-nj, Frensch2019-ff]. We use data from these studies to construct a usage-based degradation rate. The operating current density and degradation rate have been averaged over the lifetime of the cell and only square and hold wave patterns were taken from the degradation studies [Suermann2019-kz, Rakousky2017-te, Papakonstantinou2020-hj, Alia2019-lw, Li2021-nj, Frensch2019-ff]. Additionally some studies exhibited a negative degradation rate likely to a decrease in ohmic losses from membrane thinning. For this work, those data points were not considered.

Figure 4 illustrates the points taken from literature, and Equation 28 shows the correlation found. The degradation rate increases linearly with the square of the current density for $i>1$\frac{\mathrm{A}}{{\mathrm{cm}}^{2}}$$. The experiments at low current density do not follow this linear trend but instead have a relatively constant degradation rate. The degradation rate was thus modeled as a piece-wise function shown in Equation 28. The constant in the equation is a fitted parameter from the data.

Pushing electrolyzers to higher and higher current densities past the 4 $\frac{\mathrm{A}}{{\mathrm{cm}}^{2}}$ maximum presented in this work makes the economics more favorable when degradation is ignored [Chung2024-vr]. Bench scale experimentalists have already begun to push up to 10 $\frac{\mathrm{A}}{{\mathrm{cm}}^{2}}$ [Martin2022-ek]. However, operating at higher and higher current densities could potentially result in increased degradation rates if one assumes the extrapolation of the trend in Figure 4.