Prior to the description of the mathematical models, assumptions are here concisely reported. Various assumptions will be elaborated on after the following list.

• Perfect purity and total conversion are associated with EL, for which dynamics are neglected. Therefore, EL is treated as a node, splitting its feed stream into two outlets according to the stoichiometry of water splitting: H₂O → H₂ + 0.5O₂.

The outlets from an electrolyzer will contain traces of unconverted water, which should be removed. In case EL runs at high temperatures, O₂ has to be cooled prior to compression and storage. Consequently, flash condensation constitutes a reasonable option. In addition, H₂O can be removed in steam-selective polymeric membrane separators.

Biogas may contain traces of ammonia and hydrogen sulfide which should be removed before any catalytic operation. These process steps are required. Nevertheless, considering that purification operations would constitute a decoupled optimization problem from the syngas stabilization, i.e., to minimize the concentration of catalyst poisons at their outlets, and considering that the time-scale related to the effects of catalyst deactivation due to poisoning differs from the one related to short-term power fluctuation here investigated, purified biogas is considered in the control volume of the process system analyzed in this study. Furthermore, biogas could vary in its CH₄ concentration, as it may be provided by different biomass feedstock over time, or collected within the plant from different geographical areas in the region. For this reason and for the parametrization of TRI, a lower and upper bound is considered in order to design the plant for the most suitable biogas composition. In case this composition were not perfectly aligned with the actual feed gas, the ratio could be easily adjusted via membrane modules, already implemented to generate a stream of CO₂ intended for RWGS. As a lower bound for the methane concentration in the binary mixture of CO₂ and CH₄ fed to TRI, 55% is selected, whereas 75% is the upper bound, a rather high but possible value, as reported in EBA (2021) and IEA (2021).

The choice of a stoichiometric feed composition at RWGS results from the aim of minimizing the thermal effect that a surplus of H₂ or CO₂ would introduce, as larger gas volumes would be required at the reactor, thus demanding higher heat duty and reactor volume. On the other hand, feeding a molar excess of H₂ would prove reasonable to contain carbon deposition unless the catalyst shows little or no tendency in this direction, which is the case for the rhodium-supported catalyst proposed in this contribution.

The shell-side temperature of RWGS is consistent with the upper threshold of the validation range for the implemented kinetics (873–1073 K). An estimation of the electricity demand for gas compression will be provided with the discussion of results in the following section.

This section reports the main governing equations of the reactors as well as further reactor-specific assumptions. Extensive details regarding the modeling of kinetics and heat transfer can be found in the supplementary material – see Section 7.

Introducing the set of chemical components S ≔ { C O 2 , H 2 , C O , H 2 O , C H 4 , O 2 } , the reactors are described by dynamic, mono-dimensional, pseudo-homogeneous material and energy balances. The material balance for component α in molar formulation reads: ∂ x α ∂ t = − v ∂ x α ∂ l + 1 − ϵ ϵ ρ c a t c T σ α ( x α ) − x α ∑ α ∈ S σ α ( x α ) , (1) where the molar fraction x _α is differentiated in time and space. Equation 1 has to be fulfilled for all (t, l) ∈ (0, T] × (0, L), where L > 0 specifies the length of the reactor (l is the axial coordinate), T > 0 is the final time (t is the temporal coordinate).

In the energy balance, the temperature T is differentiated as: ∂ T ∂ t = − v ϵ c T C p ̃ g a s ∂ T ∂ l − 4 U d T T − T e x t − 1 − ϵ ρ c a t ∑ k ∈ K H ̃ k R k η k ϵ c T C p ̃ g a s + 1 − ϵ ρ c a t C p ̂ c a t , (2) defined for all (t, l) ∈ (0, T] × (0, L), where axial dispersion is neglected. Equation 2 introduces a summation term defined for the relevant kinetics k ∈ K≔ {reverse water-gas shift (RWGS), water-gas shift (WGS), steam-reforming of methane (SR), reverse methantion (RMETH) and catalytic oxidation of methane (OX)}. The term σ _α , concurring to (1), is the overall molar generation rate for the single component α, which reads: σ α = ∑ k ∈ K ν α , k R k η k . (3) State-of-the-art kinetics reported by Xu and Froment (1989) and De Groote and Froment (1997) are implemented for TRI. Combustion kinetics within TRI were adapted from Trimm and Lam (1980) by De Smet et al. (2001) for supported Ni catalysts. Effectiveness factors for this heterogeneous model are constant and were retrieved from De Groote and Froment (1997). Reversible kinetics of RWGS were taken from Richardson and Paripatyadar (1990), thus based on a Rh/γ −Al₂O₃ catalyst. The effectiveness factor of 0.3 observed by the authors is also incorporated in the calculations for this contribution. For adiabatic operations, i.e., inside the TRI reactor, the overall heat transfer coefficient U is 0, whilst it is U ≠ 0 for the RWGS mildly-endothermic, multitubular reactor.

The axial mole-averaged interstitial velocity v is defined by a total mass balance in quasi-steady state assumption. Mole and mass-averaged velocities coincide if dispersion is neglected. Thus, the molar-based axial velocity defined between the reactor feed (superscript ⁰) and the generic reactor section reads v = ∑ α ∈ S N α 0 W α c T A c r o s s ϵ ∑ α ∈ S x α W α . (4) The momentum balance, typically dominated by friction, reduces to the Ergun equation d p d l = − 1 − ϵ 2 ϵ 3 150 v ϵ 1 − ϵ d p 2 μ m i x + 1.75 1 − ϵ ϵ 3 ρ g a s v 2 ϵ 2 d p . (5) The time-dependency of momentum is neglected. Equations 1, 2, and 5 are completed by Dirichlet-type boundary conditions at each reactor inlet, respectively x α t , l = 0 = x α , 0 , T t , l = 0 = T 0 , p l = 0 = p 0 , (6) where x _α,0, T 0 and p ₀ are the molar fraction of component α, temperature and pressure at the reactor feed, pre-specified and constant over time.

After discretizing the reactors by equally-spaced finite volumes in the spatial direction and approximating the advection contributions with the upwind scheme, the resulting system of modeling equations is a semi-explicit differential-algebraic equation (DAE) system of index 1 on the time horizon 0 , T . Thus, it can be written in the form y ̇ ( t ) = f ( t , y ( t ) , z ( t ) ) ,   0 = g ( t , y ( t ) , z ( t ) ) , (7) where y ̇ is an abbreviation for d y d t and suitable initial values for (7) will be discussed in Section 5. Here:

• y ( t ) ∈ R d × [ 0 , T ] collects all differential variables, i.e., all variables that are differentiated with respect to time. These are the mole fractions coming out of the discretizations inside the RWGS and TRI reactor.

A preliminary optimization problem is set to identify suitable values of the design parameters that maximize the selectivity towards syngas, that is min z , Φ D O B J T R I → min z , Φ D 1 S e l e c t i v i t y → min z , Φ D N C H 4 + N C O 2 + N H 2 O N H 2 + N C O o u t , T R I . (8) Here, Equations 1, 2 are implemented in their steady-state formulation such that time derivatives vanish. The resulting modeling equations serve as equality constraints in a nonlinear programming (NLP) problem, that is min z S , r , Φ D r O B J r , s.t. 0 = f S , r ( z S , r , Φ D r ) , l b r < { z S , r , Φ D r } < u b r , and r ∈ { T R I } , (9) where the vector Φ _D collects the optimization variables: tube length L, tube diameter d _T , feed molar flowrates N α i n ( α ∈ S ), temperature T 0 and pressure p ₀, and the vector z _S,r which includes all reactor states. The function f _S,r incorporates the conservation laws (1), (2), and (5), where the subscript S indicates the steady-state formulation. Additional constraints bound the syngas composition as H 2 / C O ∈ 1.95 , 2.1 , while the productivity of carbon monoxide is fixed. Problem parameters, constraints and bounds are defined in Matlab v2018b, combined with CasADi v3.5.3, a symbolic framework for algorithmic differentiation and numeric optimization developed by Andersson et al. (2019). The nonlinear program is solved by IPOPT v3.12.3, running with the linear solver mumps, see, e.g., Wächter and Biegler (2006). The solver IPOPT identifies local optima. Therefore, in order to identify a global solution, optimization variables are randomized within bounds. Preliminary simulations in steady-state based on the randomized guesses are then ran in order to provide IPOPT with reasonable initial guesses. Optimal solutions are then compared and the design set ensuring the best objective is selected. The reactor is discretized with 150 equally-spaced points. Optimization variables, bounds and optimal values are reported in Table 2.

The selected bounds for the feed pressure to TRI are in line with the findings on the role of coke formation in catalytic partial oxidation of methane by De Groote and Froment (1997) and with the prevailing literature on partial oxidation and tri-reforming. Above 25 bar and for elevated inlet temperatures, coke deposition is neglected.

Moreover, the selected lower bound for the feed pressure at TRI, which coincides with its optimal value reported in Table 2, is consistent with the requirement of a typical Fischer-Tropsch synthesis reactor, i.e., no pressurization is required before this downstream syngas application.

At the best local optimum, the molar selectivity towards H₂ and CO is 2.28. The feed ratio between CH₄ and the sum of CO₂ and CH₄ is 0.75, which corresponds to the upper bound defined for the CH₄ molar content in this binary mixture (see Section 3.1 for modeling assumptions).

The set of relevant design parameters and operating conditions is reported in Table 3.

Accounting for its drop along the axial coordinate, pressure is slightly above atmospheric, in agreement with the range of validation for the kinetic model by Richardson and Paripatyadar (1990). At the prevailing temperature, pressure and feed composition, the selected feed and shell-side temperature are high enough to prevent from a thermodynamically relevant contribution of the methanation reaction, therefore neglected. The number of discretization points for RWGS simulations and optimal control is 150.

Standard Runge-Kutta methods (RKMs) can be adapted to DAEs according to Hairer and Wanner (1996), but the matrix of the tableau has to be invertible such that implicit methods have to be chosen. Thus, a simple implicit Euler scheme is sufficient for the adjoint Equation 17 as it is a linear DAE. Regarding (14), the forward simulation is computationally more demanding, such that parareal is briefly reviewed here, see, e.g., Lions et al. (2001); Baffico et al. (2002); Garmatter et al. (2021), in order to speed up the solution time of the overall OCP.

The main idea of parareal is to decompose the global time domain [0, T] into N _t smaller subdomains. Given initial values on each of these subdomains, the global time-dependent problem in each iteration of the method splits up into N _t many local problems on these subdomains, which can then be solved in parallel. The initial values can be generated using a coarse integrator, which should be cheap but can be inaccurate, and the subproblems are then solved in parallel using a fine integrator, which has to be accurate and is thus more expensive. Afterwards, the next iterate of parareal is generated via a correction step, where the fine solutions of the subproblems are used to correct the coarse sequential integrator.

Introducing the general time grid 0 =: t 0 < t 1 … < t N t ≔ T with variable step sizes h _n ≔t _n+1 − t _n , n = 0, …N _t − 1, G ( t n + 1 , t n , h n , w n k ) denotes the coarse integrator, that integrates on the subdomain [t _n , t _n+1] with step-size h _n and provides an inaccurate approximation to w(t _n+1), the solution of (14), using the initial values w n k . Here, k indicates the iteration number of parareal. The fine integrator has to be more accurate and can thus be more expensive. This can be achieved by using a higher order integration method or by operating on a refinement of the time grid or a combination of both. In this article, the fine integrator will always be a higher order method and it can additionally operate on a refinement of the time grid such that F ( t n + 1 , t n , R r e f , w n k ) denotes the fine integrator, that integrates on [t _n , t _n+1] using R r e f ∈ N time steps and provides a more accurate approximation to w(t _n+1) using the initial values w n k . Thus, for R _ref = 1 both integrators use the same grid and for R _ref ≥ 2 the fine solver F uses a refined grid. With this notation at hand, parareal is described in Algorithm 1.

Algorithm 1 is terminated either after a fixed amount of K max ∈ N iterations or as soon as the relative change in the iterate w k + 1 − w k / w k is below some tolerance ɛ _tol > 0. Note that in line 10, the values of the coarse integrator G n + 1 k = G ( t n + , 1 , t n , h n , w n k ) have been calculated in the previous iteration already and thus can be reused. Regarding the integrators, G should be fast but at the same time an implicit RKM such that the implicit Euler method is chosen and for F a higher order method is required such that the Lobatto IIIC method is chosen.

A first optimal control problem based on a tracking-type objective function with an added regularization term is considered. That is, for a given desired state x d ( t ) : R → R d + q , the optimal control problem reads min u ∈ R 2 J 1 ( u ; β ) ≔ 1 2 ∫ 0 T F ( u ) − x d ( t ) 2 2 + β u 2 2 d t , and u l ≤ u ≤ u u . (OCP1) This resembles the scenario: given a desired state x _d (t) and an initial control u i n i t ∈ R 2 , problem (OCP1) wants to find an optimal control u o p t ∈ R 2 for which the corresponding state F ( u o p t ) best fits the desired state x _d (t). As the desired state will be calculated based on a chosen desired control u d ∈ R 2 , this problem becomes a parameter-identification problem and serves the purpose of validating the mathematical approach rather than resembling an actual application.

Regarding the initialization of the DAE system (14), the initial values w ( 0 , u ) ∈ R d + q are constructed as follows: an intermediate vector w ̃ 0 ∈ R d + q is formed that consists of fixed values for the mole fractions, the temperature and the pressure at each discretization point inside the reactors. Furthermore, the control u affects the algebraic variables of w ̃ 0 corresponding to the reactor outlet flowrates and the H₂-integration. Next, w ̃ 0 serves as the initial value for a Newton method that solves F ̃ ( 0 , w , u ) = 0 . The solution of this nonlinear problem then is w(0, u), the actual initial values for the DAE (14) such that these initial values are always consistent. Furthermore, w _u (0) has only nonzero values in the algebraic part, such that Mw _u (0) = 0 and the last term in (18) vanishes.

As mentioned in the beginning of Section 5 the control u and the EL inlet flowrate EL_in are now time-dependent. Thus, EL_in is a function of time predefined in the model and a time-dependent control has to be introduced. Remembering the time-grid 0 = : t 0 < t 1 < … < t N t ≔ T , the control function is defined to be piecewise constant in time, that is u : R → R 2 : t ↦ ∑ i = 1 N u i χ i ( t ) , with χ i ( t ) ≔ 1 , t ∈ ( t i − 1 , t i ] , 0 , else , (19) where the coefficients u i ∈ R 2 contain the control values on the time interval (t _i−1, t _i ]. Together with the control at the initial time point u ( t 0 ) = u ( 0 ) = u 0 ∈ R 2 and abusing Matlab notation, these coefficients can be collected as u ≔ [ u 0 ; u 1 ; … ; u N t ] ∈ R p ,  with  p ≔ 2 ( N t + 1 ) , (20) such that the overall control dimension p scales with the number of time steps. The OCP presented in the following shall now capture the scenario described in the Introduction: due to external influences (intermittency of electricity), the inlet flowrate to EL is changing according to a piecewise linear profile and the aim of the OCP then is to stabilize the plant to still yield a desired syngas ratio S R d ∈ R and flowrate S F d ∈ R in the product stream of the plant. Letting S R ( t ) : R → R : t ↦ f 1 ( F ( u ( t ) ) ) and S F ( t ) : R → R : t ↦ f 2 ( F ( u ( t ) ) ) be functions that compute the syngas ratio and flowrate at the outlet at time t, based on the value of the control u(t) and thus the state F ( u ( t ) ) , the optimal control problem reads min u ( t ) ∈ R 2 J 2 ( u ( t ) ) ≔ 1 2 ∫ 0 T c S R ( S R ( t ) − S R d ) 2 + c S F ( S F ( t ) − S F d ) 2 d t , and u l ≤ u ≤ u u . (OCP2) Here, c S R ≔ max { S R d , S F d } S R d and c S F ≔ max { S R d , S F d } S F d are scaling constants that ensure that both contributions to the objective function are equally weighted. Furthermore, the dynamic EL inlet flowrate profile is again hard-coded in the model such that it does not appear as a constraint in (OCP2). Together with the dynamic inflow to the RWGS unit (the first control component), the H₂-integration to the syngas outlet stream at each time point can then be determined according to (11).

The initialization of the system (14) for (OCP2) is similar to the initialization described in Section 5.2. The only difference being that the entries of w ̃ 0 corresponding to the mole fractions, the temperature and the pressure inside the reactors come from a steady-state-simulation that reflects the desired syngas ration and flowrate but uses a constant-in-time EL inflow profile. The control values u ₀ then again affect the algebraic variables corresponding to the reactor outlet flowrates and the H₂-integration and the resulting w ̃ 0 serves as the initial guess for the Newton method resulting in the consistent initial values.

From a theoretical perspective, it is clear that the formulation in (OCP2) differs from (OCP) due to the time-dependency of the control. Nonetheless, the continuous adjoint method can be applied to calculate the objective gradient, where the necessary details can be found in (Gerdts, 2012, Section 5.3.3).

The parameters are selected from Xu and Froment (1989) and reported in Table 4. The governing kinetic expressions are

• D E N = 1 + K C O p C O , [ b a r ] + K H 2 p H 2 , [ b a r ] + K C H 4 p C H 4 , [ b a r ] + K H 2 O x H 2 O x H 2 ,

where the kinetics parameters k _i and K _j related to reaction and adsorption result from the following Arrhenius-like relations k i = A ( k i ) exp − E ( k i ) R g a s T a n d K j = A ( K j ) exp − E ( K j ) R g a s T . (23) Values of coefficients A(k _i , K _j ) and E(k _i , K _j ) are listed in Table 4.

RWGS kinetics are adapted from Richardson and Paripatyadar (1990). D E N = 1 + K C O 2 p C O 2 , [ a t m ] + K H 2 p H 2 , [ a t m ] (24) and R R W G S = 1 0 3 k R W G S K C O 2 K H 2 p C O 2 , [ a t m ] p H 2 , [ a t m ] − p C O , [ a t m ] p H 2 O , [ a t m ] K e q R W G S D E N 2 , (25) Where kinetics parameters k _i and K _j related to reaction and adsorption, result from Arrhenius-like relations k R W G S = 350 ⁡ exp − 81030 R g a s T , K C O 2 = 0.5771 ⁡ exp 9262 R g a s T , K H 2 = 1.494 ⁡ exp 6025 R g a s T . (26) Effectiveness factor is set to 0.3, as reported in Richardson and Paripatyadar (1990).

The overall heat transfer coefficient reads U = 1 α e f f − 1 . (27) The shell-side heat transfer coefficient is not accounted for. Instead, a constant skin temperature is assumed along the pipe. The definition of α _eff for non-adiabatic packed-bed reactors is provided by Martin and Nilles (1993). 1 α e f f = 1 α w + d T 8 Λ r . (28) Here, α _w and Λ _r are retrieved from Bauer and Schlünder (1976) and discussed by Tsotsas (2010) as well as Martin and Nilles (1993). α w = 1.3 + 5 d T / d p λ b e d λ m i x + 0.19 ρ g a s v ϵ d P μ m i x 0.75 μ m i x C p ̃ g a s λ m i x 0.33 λ m i x d p (29) and Λ r = λ b e d + v ϵ c T C p ̃ g a s d p 8 2 − 1 − 2 d T / d p 2 . (30) Embedded in the definition of α _w , the heat conductivity across the packed-bed λ _bed is defined in Tsotsas (2010) by the following steps:

• λ _bed = k _bed λ _mix ,