Power-to-Syngas: A Parareal Optimal Control Approach

A chemical plant layout for the production of syngas from renewable power, H2O and biogas, is presented to ensure a steady productivity of syngas with a constant H2-to-CO ratio under time-dependent electricity provision. An electrolyzer supplies H2 to the reverse water-gas shift reactor. The system compensates for a drop in electricity supply by gradually operating a tri-reforming reactor, fed with pure O2 directly from the electrolyzer or from an intermediate generic buffering device. After the introduction of modeling assumptions and governing equations, suitable reactor parameters are identified. Finally, two optimal control problems are investigated, where computationally expensive model evaluations are lifted via parareal and necessary objective derivatives are calculated via the continuous adjoint method. For the first time, modeling, simulation, and optimal control are applied to a combination of the reverse water-gas shift and tri-reforming reactor, exploring a promising pathway in the conversion of renewable power into chemicals.


INTRODUCTION
In the context of energy transition towards carbon-free drivers for the chemical and fuel industry, electricity from renewable sources plays a primary role. Availability constraints are associated with the use of renewable electrical power. Engeland et al. (2017) introduced the notion of climate related energy (CRE), which encompasses wind-, solar-power and, partly, hydro-power. Typically, its intermittent provision on the daily and hourly time-scale is attributable to the temporal trends of the related natural source and other factors, such as wind velocity, solar radiation, and atmospheric precipitation. The intermittent nature of the supply as well as the consumer demand contribute to price volatility. As reported by Schill (2014), positive residual power loads in the grid occur whenever demand exceeds the generation of CRE.
Currently, by means of relatively flexible operations and proper forecasts, state-of-the-art power plants can adjust their outputs to meet market demands and manage fluctuations in the provision of renewable power. Reservoir-type hydro-power, i.e., conversion of surplus electricity into potential energy by pumping water into elevated basins, can also contribute to this strategy. Nevertheless, the number of such reservoirs is limited due to morphological constraints. Moreover, as reported in Sinn (2017), seasonal fluctuations are expected to require most of the storage capacity if the share of renewable electricity was to become predominant.
In this framework, Power-to-X processes are relevant contributors to the energy transition: they are conceived to take advantage of unsteady power inputs, thus transformed into valuable chemicals and energy vectors. In particular, as highlighted in a review article by Wulf et al. (2020), the role of H 2 is of cardinal importance in the conversion of electricity. As a matter of fact, projects and industrial applications related to water electrolysis (EL) of increasing capacities are expected in the near future. This contribution shall focus on Power-to-Syngas towards Fischer-Tropsch, for which a molar ratio H 2 /CO of 2 must be ensured. H 2 and CO 2 can be converted to syngas via the reverse water-gas shift reaction in a non-adiabatic fixed-bed reactor (RWGS), possibly supported by a number of candidate catalysts, as thoroughly reviewed by Daza and Kuhn (2016). Its mild endothermicity allows for a moderate thermal input to sustain the generation of CO from CO 2 .
As the inherent complexity of a tri-phase Fischer-Tropsch synthesis reactor benefits from a steady syngas supply, shortterm, e.g., hourly, intermittency could be levelized via H 2 buffering, although resulting in higher investment and operating costs in terms of EL volumes and electricity supply. As an example, if the power load required to attain the nominal flowrate of syngas is 50 MW, 10 wind turbines of 5 MW each would fulfill the duty. Nevertheless, a seasonal shortage of renewable power would require the installation of X additional turbines and a dedicated buffering system for the storage of H 2 . Furthermore, Kaiser et al. (2013) report on the high energy intensity required for H 2 liquefaction, amounting up to 30% of its lower heating value, and on its low density, both in gaseous and liquefied state. Moreover, buffering may not guarantee sufficient productivity levels over long term, e.g., seasonal, shortage of renewable power supply. Besides H 2 , EL generates O 2 which could be used to sustain the endothermicity of RWGS. Assuming that the thermal duty of RWGS equals its standard enthalpy of reaction (35.3 kJmol −1 CO at 950 K), resulting in a conversion of 42.7% of an equimolar feed mixture in isothermal operations and for the lower heating value (LHV) of CH 4 (800 kJ mol −1 ), 0.0882 mol O 2 per mole of CO produced are required for heat generation. EL generates 1.5 mol O 2 mol −1 CO if a syngas ratio H 2 /CO of 2 is to be attained. Consequently, 94.1% of the O 2 produced within the plant does not contribute to the generation of syngas. In conclusion, the large majority of O 2 is not utilized within the process: even if the stabilization of syngas were attained via H 2 storage strategies, it would be reasonable to consider O 2 as a valuable by-product to be stored within suitable buffering devices in the plant.
Another route to syngas is tri-reforming of methane (TRI), traditionally meant to convert combustion off-gases and natural gas within a catalytic, adiabatic reactor fed with oxygen or air, as in Song and Pan (2004). The oxidation of CH 4 contributes to lower the levels of carbon deposition and to allow for a tunable outlet composition (Vita et al. (2018)). Biogas from anaerobic digestion and CO 2 can partially or completely substitute the traditional feedstock. The use of O 2 results in high adiabatic temperature peaks although, if compared with air, it requires smaller process volumes and separation costs.
Considering the aforementioned possible downsides of H 2 buffering and the possibility to use a TRI reactor, this contribution proposes the combined implementation of a RWGS and a TRI reactor as an alternative to a H 2 buffering strategy. Here, O 2 is provided either directly from the EL or indirectly, as suitable storage devices are assumed after the splitting of water. If the stored O 2 were not sufficient to feed TRI on a seasonal basis, it could be used to enrich air, then fed to the reactor. In order to ensure a steady syngas supply at the specified H 2 /CO molar ratio of 2 under variable electricity provision, the continuous switching between RWGS and TRI operations should be enforced via appropriate control strategies.
After the introduction of modeling assumptions and reactor models, design parameters and operating conditions for the reactors are identified. Afterwards, the overall optimal control strategy is outlined where the final goal is to meet syngas specifications under a time-varying electricity supply. From the mathematical side, parareal shall be employed to speed up expensive model evaluations during the solution of the optimization problems. parareal, see, e.g., Lions et al. (2001); Baffico et al. (2002), aims at decomposing the global time domain into several smaller domains and, given initial values on these subdomains, the global problem splits up into local subproblems that, in each iteration, can be solved in parallel. The initial values are generated using a cheap but possibly inaccurate coarse time integrator and the subproblems are then solved in parallel using an expensive but accurate fine time integrator. Initial values for the next iteration are generated in a correction step utilizing the information from the fine integrator. The derivatives of the objective function that are necessary for the optimization are calculated via the continuous adjoint method, see, e.g., Cao et al. (2002Cao et al. ( , 2003, since the discrete adjoint method is not available due to the usage of parareal as a forward integrator.
In a numerical investigation, two optimal control problems are considered: at first, an academic problem shall verify the approach of combining parareal with the continuous adjoint method. The second optimal control problem then targets the scenario previously described: a desired syngas ratio and flowrate shall be attained over a time horizon while a power fluctuation results in a varying temporal provision of H 2 from EL to the plant. Both problems will be solved once with parareal for the forward integration and once using the fine integrator previously used inside parareal. A time comparison as well as a comparison of the obtained optimal solutions is then of interest. Figure 1 represents the plant layout. EL converts H 2 O into H 2 and O 2 . H 2 can either contribute with CO 2 to feed the RWGS reactor, or be directly supplied to the product stream to adjust the H 2 /CO molar ratio in the plant outlet. Conversely, O 2 is either supplied to TRI to sustain the endothermicity of reforming reactions, or suitably buffered for later use. Purification strategies prior to Fischer Tropsch are not investigated in this article.

PLANT LAYOUT
In this framework, the most rational source of CO 2 is biogas, part of which could be split into CH 4 and CO 2 via a suitable stateof-the-art separation strategy as reviewed by Maggi et al. (2020), e.g., membranes or PSA. The resulting CH 4 stream would then be used to sustain the heat demand of RWGS. Surplus CH 4 could be injected into the natural gas grid.

MATHEMATICAL MODEL
Key modeling assumptions and main governing equations for the reactors are reported in this section. For the sake of clarity, important symbols are introduced within the text, whereas a broader list of symbols is reported in Table 1.

Modeling Assumptions
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 2 O → H 2 + 0.5O 2 . • H 2 can be integrated directly from EL into the product stream. • The difference between O 2 from EL and O 2 to TRI approximates the accumulation of O 2 in the buffering system, which is not directly modeled. • Biogas is a binary mixture of CH 4 and CO 2 for example provided by anaerobic digestion. • Feed composition, temperature and pressure at the each reactor feed are fixed, whereas their flowrates can vary. • RWGS is fed with an equimolar mixture of CO 2 and H 2 . The desired feed composition of TRI shall be identified via optimization in Section 4. • The shell-side temperature at RWGS, denoted as T ext , is assumed to be constant along the reactor axis and equal to 1073 K, which is assumed to coincide with the outer skin temperature. • The electricity demand at EL largely surpasses other contributions, which are therefore neglected in the optimal control application.
The outlets from an electrolyzer will contain traces of unconverted water, which should be removed. In case EL runs at high temperatures, O 2 has to be cooled prior to compression and storage. Consequently, flash condensation constitutes a reasonable option. In addition, H 2 O can be removed in steamselective 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 4 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 2 intended for RWGS. As a lower bound for the methane concentration in the binary mixture of CO 2 and CH 4 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 2 or CO 2 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 2 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 FIGURE 1 | Overall Plant layout: the RWGS and TRI reactors produce raw syngas, where H 2 and O 2 are supplied via EL. H 2 may be fed to RWGS and/or bypass the reactor to the syngas product stream whereas excess O 2 may be stored in a buffering device. CO 2 could be obtained from biogas using a suitable separation strategy.
Frontiers in Energy Research | www.frontiersin.org September 2021 | Volume 9 | Article 720489 for gas compression will be provided with the discussion of results in the following section.

Modeling of Reactors
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 materialsee Section 7. Introducing the set of chemical components Sd{CO 2 , H 2 , CO, H 2 O, CH 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: 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: 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 ∈ Kd {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: 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/c −Al 2 O 3 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 0 ) and the generic reactor section reads The momentum balance, typically dominated by friction, reduces to the Ergun equation The time-dependency of momentum is neglected. Equations 1, 2, and 5 are completed by Dirichlet-type boundary conditions at each reactor inlet, respectively where x α,0 , T 0 and p 0 are the molar fraction of component α, temperature and pressure at the reactor feed, pre-specified and constant over time.

Semidiscretized Plant Model
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)), where _ y is an abbreviation for dy dt 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.
T] collects all algebraic variables. These are the pressure values inside the reactors as well as the outlet molar flowrates from the reactors and the H 2 -integration into the plant outlet stream.

REACTOR DESIGN
Prior to optimal control, reactor designs for TRI and RWGS are identified separately, ensuring that the desired productivity of 270 ton CO /day and a syngas ratio of H 2 /CO 2, suitable for Fischer-Tropsch synthesis, are attained in each case.

Design of TRI
A preliminary optimization problem is set to identify suitable values of the design parameters that maximize the selectivity towards syngas, that is 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 zS,r ,ΦDr OBJ r , where the vector Φ D collects the optimization variables: tube length L, tube diameter d T , feed molar flowrates N in α (α ∈ S), temperature T 0 and pressure p 0 , 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 /CO ∈ [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 Frontiers in Energy Research | www.frontiersin.org September 2021 | Volume 9 | Article 720489 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 2 and CO is 2.28. The feed ratio between CH 4 and the sum of CO 2 and CH 4 is 0.75, which corresponds to the upper bound defined for the CH 4 molar content in this binary mixture (see Section 3.1 for modeling assumptions).

Design of RWGS
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.

TRI Reactor
Reactor profiles reflect the typical adiabatic operation: non-zero gradients are reported at the beginning, flattening out to zero once the energy input from the catalyzed combustion of CH 4 is exhausted by the endothermic reactions. The selected NLP formulation and objective determine a solution for the reactor length at its upper bound, although the gradients flatten out before its first half. Consequently, the reactor length is reduced by 90% of the upper bound, whereas the value of the remaining optimization variables is set to their optimum. Steady-state simulation profiles are reported in Figure 2: reforming contributions exploit the heat from catalyzed combustion within the very first reactor section. After the adiabatic peak, temperature decreases and stabilizes around 1200 K. The short reactor length determines a negligible pressure drop, in agreement with Chein et al. (2017), and the velocity reflects the temperature profile. Position and magnitude of the hightemperature peak at the very inlet of the reactor is consistent with scenarios presented by Chein et al. (2017), Arab Aboosadi et al. (2011, and Rezaei and Catalan (2020). For the reason that the maximum adiabatic temperature undermines the thermal stability of materials, recommendable studies concerning sidefeeding strategies of O 2 by means of membrane and sequential injection points have been proposed by Alipour-Dehkordi and Khademi (2020) and Rezaei and Catalan (2020), although disregarded from the current contribution for sake of mathematical simplicity in the optimal control application. Furthermore, the small ratio between the resulting reactor length and diameter could lead to important radial gradients and an uneven gas distribution, possibly resolved by partitioning the feed stream into a bundle, as proposed by Alipour-Dehkordi and Khademi (2020).

RWGS Reactor
Reactor profiles are shown in Figure 3. The temperature profile drops to a minimum within the first quarter of reactor length, where the mass action is at its maximum. Gradually, it shifts towards the shell-side temperature. Temperature and velocity are directly proportional, whereas the pressure drops almost linearly along the reactor length. Simulation results indicate a conversion of 44.7%.
Assuming that EL operates at 1073 K, the ideal electrical power demand, which equals the Gibbs free energy of reaction, is 188 kJmol −1 H 2 . Therefore, the nominal flowrate of H 2 to RWGS (247.5 mol H 2 s −1 ) is ensured if EL is supplied with 46.53 MW. Posdziech et al. (2019) identified an EL efficiency of 82%, based on the lower heating-value of H 2 (LHV H 2 ), corresponding to 240 kJmol −1 H2 , which outlines a power demand of 72.4 MW. For the given operation, a  syngas ratio H 2 /CO of 2 is attained if an integration of 85 mol H 2 s −1 is accounted for, which corresponds to a total power demand at EL of 62.5 MW, based on the Gibbs free energy, and 97.3 MW, based on the LHV H2 efficiency. CH 4 is a candidate fuel to sustain the heat demand of RWGS, possibly provided from biogas. Given that its lower heating value (LHV CH4 ) is 800 kJmol −1 , an estimate of the required flowrate is provided by the following calculation, accounting for the reactor discretization (N z points in axial coordinate and axial discretization segments of length Δz): resulting in 5.8 mol s −1 of CH 4 , 5% of the molar flowrate required for nominal operations of TRI towards the maximum selectivity to syngas (see Table 2). If biogas has a molar concentration of 60% in CH 4 , 9.7 mol s −1 must be separated. Consequently, nominal RWGS operations require considerably less CH 4 than TRI, the latter virtually demanding no electricity other than compression duties.

Oxygen Utilization
Given that the stoichiometric combustion of 5.8 mol CH 4 s −1 requires 11.6 mol O 2 s −1 , nominal RWGS operations generate an excess of 155 mol O 2 s −1 which can be buffered to allow for the operation of TRI in lack of renewable electricity. Assuming that RWGS (TRI is off) and TRI (RWGS is off) are respectively operated for 50% of a given time-horizon, the ratio of O 2 generated (not intended for combustion) to O 2 fed to TRI (

OPTIMAL CONTROL FOR THE PLANT MODEL
In this section, two different optimal control problems (OCPs) related to the plant model (7) are considered. The control vector u (u RWGS , u TRI ) u ∈ R 2 contains the inlet flowrates to the RWGS and TRI reactors, where the feed composition is fixed as assumed in Section 3.1. The inlet molar flowrate to EL, denoted as EL in , is not treated as a control variable but is predefined in the model. Thus, one obtains the H 2 -integration to the product stream as where the composition of u RWGS , the total inlet flowrate to RWGS, is an equimolar mixture of CO 2 and H 2 as described in Section 3.1. Finally, although the buffer device is not modeled as mentioned in Section 3.1, assuming that O 2 is exclusively intended for the feed to TRI, the accumulation of O 2 inside the buffer can be calculated as  is the oxygen required to burn the required amount of methane calculated in (10).
In the first OCP, the control u and the EL inlet flowrate EL in are constant in time in order to clearly describe the approach. In the second OCP these quantities are then time-dependent functions. As a result, they have to discretized in time such that the amount of control variables scales with the amount of time steps and the problem as a whole becomes more challenging from the optimization point of view.
Introducing the combined state w(t)d[y(t), z(t)] u ∈ R d+q and the matrix Md I d×d 0 d×q 0 q×d 0 q×q , where I d×d ∈ R d×d is the identity matrix, the system can be rewritten as where the nonlinear functionF : R × R d+q × R 2 → R d+q encodes f and g from (13). Based on this formulation, the abstract solution operator F : u 1w(t) maps a control u to the state solution w(t) of (14). Using this notation, a general OCP in reduced form inlcuding box constraints on the control is proposed as In this reduced form of the OCP, the state w is not an optimization variable and thus only implicitly known to the optimizer through the solution operator F . As a result, whenever the objective is to be evaluated, the solution operator has to be evaluated for the current value of the control, meaning that (14) has to be solved. Before the numerical solution of (14) is tackled and two specific OCPs are introduced, it is clear that any optimizer applied to (OCP) will require dΓ du , the gradient of the objective function with respect to the control. Thus, the continuous adjoint method from Cao et al. (2003) for its calculation is briefly reviewed.
Introducing a yet unspecified variable λ ∈ R d+q , the augmented objective function is formed as where here and in the following the dependency of λ, w, _ w, F, and associated quantities on time is often suppressed. Since F(w, _ w, t, u) 0 due to the reduced approach, the sensitivity of Γ with respect to u is where the subscripts indicate partial derivatives and the integrals are to be taken component-wise. Integration by parts for the term T 0 λ u F _ w _ w u dt together with the fact that F _ w ≡ M and further letting the adjoint state λ(t) ∈ R d+q satisfy the adjoint equation the final sensitivity reads Thus, the computation of the gradient inside an optimizer requires • the solution of the forward problem (14), • the solution of the adjoint problem (17), • the evaluation of the gradient (18).

Regarding the adjoint Equation 17
, which is a linear DAE, it is sufficient to choose initial values such that λ(T) u M 0, which is easily obtained according to (Cao et al., 2003, Eq. (9)). Thus, the last term in (18) becomes λ(0) u Mw u (0) and one has to calculate w u (0), the sensitivity of the initial state with respect to the control. The initialization of the system will be covered separately for each OCP after the numerical solution of (14) and (17) is discussed in the following section.

Numerical Time Integration
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 dT with variable step sizes h n dt n+1 − t n , n 0, . . .N t − 1, G(t n+1 , t n , h n , w k n ) 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 k n . 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 ref , w k n ) denotes the fine integrator, that integrates on [t n , t n+1 ] using R ref ∈ N time steps and provides a more accurate approximation to w(t n+1 ) using the initial values w k n . 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 k n+1 G(t n+,1 , t n , h n , w k n ) 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.

First OCP
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 This resembles the scenario: given a desired state x d (t) and an initial control u init ∈ R 2 , problem (OCP1) wants to find an optimal control u opt ∈ R 2 for which the corresponding state F (u opt ) 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 vectorw 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 ofw 0 corresponding to the reactor outlet flowrates and the H 2 -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.

Numerical Results
Regarding the physical discretization and various model parameters, the same values as in Section 4 are chosen in this section and in the later Section 5.3.1. Furthermore, Matlab v2021a was invoked to obtain these results.
For the solution of (OCP1), Matlab's fmincon routine with its SQP solver is called, where the objective gradient is calculated via (18), the adjoint state is calculated via (17), and this adjoint equation is solved via the implicit-Euler method. The necessary derivatives ofF are calculated via ADiMat, see Bischof et al. (2002). The integrals in (OCP1) and (18) are approximated via the trapezoidal rule. Regarding fmincon, default options and thus default termination criteria are used.
In this academic example, it is T 1 s and the time horizon [0, 1] is discretized with N t 100 equidistant points resulting in h 0 / h Nt −1 1 Nt . For the evaluation of F (u) during the optimization, either the fine integrator F or parareal introduced in Section 5.1 are chosen where parareal is invoked with K max 10 and ε tol 10 −6 . Furthermore, via R ref 2, the fine integrator always operates on a refined grid. The computations were carried out on a machine with 60 Intel(R) Xeon(R) CPU E7-4880 v2 @ 2.50 GHz cores, where 50 cores where used for the parallel step.
The experiment consists of solving (OCP1) once using the fine integrator F for the evaluations of F (u) and once using parareal and the achieved speed-up is then of interest. To keep the comparison fair, Matlab is in any case only allowed to use one computational thread inside the RKMs. As discussed in Garmatter et al. (2021), the obtainable speedup using parareal is reduced when using an expensive coarse integrator. To partly lift this downside, the scenario of a cheaper coarse solver is mimicked by selecting 10 −4 as the tolerance for the Newton-solver inside G where the tolerance inside F is 10 −8 .
The remaining experiment setup is (numbers associated to the control u always have the unit of measurement mol.s −1 ) • x d (t) F (u d ) with u d (30, 1) u and EL in 30 mol.s −1 is obtained via the coarse solver G to prevent inverse crime; • u init (15, 15) u simulating a generically active plant; • u l (0.1, 0.1) u , u u (60, ∞) u , and β 10 −6 . Here, the upper bound for the RWGS corresponds to two times the inlet flowrate to EL to prevent negative H 2 -integration values due to (11).
Regarding the solution quality, the solution of (OCP1) using F inside the SQP solver of fmincon was and the solution using parareal was u par ≈ (30.0003, 0.9946) u , with u par − u d ||u d || ≈ 1.8025 · 10 −4 , Regarding the solution time, the optimization using the fine integrator required 8186 s and the optimization using parareal required 3181 s resulting in a speed-up of 2.57. Here, parareal required on average only one iteration, where it was already discussed in Garmatter et al. (2021) that this quick convergence stems from the used implicit RKM as a coarse solver. The downside of this approach is that parareal on average did spend 70.70 s per iteration in the coarse solver and only 8.73 s in the fine solver. Thus, the obtained speed-up could be improved by using N t 100 cores (effectively halving the time spent in the fine solver) but it is clear that the coarse solver does indeed limit the maximum obtainable speed-up. Thus, the use of an implicit explicit (IMEX) numerical scheme, see, e.g., Boscarino (2007); Steiner et al. (2014) for IMEX methods for DAEs, could be investigated to leverage this downside.

Second OCP
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 Nt dT, the control function is defined to be piecewise constant in time, that is Frontiers in Energy Research | www.frontiersin.org September 2021 | Volume 9 | Article 720489 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 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 SR d ∈ R and flowrate SF d ∈ R in the product stream of the plant. Letting SR(t): R → R: t1f 1 (F (u(t))) and SF(t): R → R: t1f 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 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 2 -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 ofw 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 constantin-time EL inflow profile. The control values u 0 then again affect the algebraic variables corresponding to the reactor outlet flowrates and the H 2 -integration and the resultingw 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).

Numerical Results
A time horizon of 1 h (resulting in T 3600 s) is chosen to reflect the realistic scenario proposed in (OCP2). As a result, the number of equidistant time-steps in the time-grid is increased to N t 200.
Regarding the experiment setup, a steady-state simulation is performed where the control values and the EL inflow are both constant in time for the values of u (465, 20) u and EL 313 mol.s −1 (numbers associated to the control u again have the unit of measurement mol.s −1 ). This results in a desired syngas ratio SR d 1.9998 and desired syngas flowrate SF d 111.44 mol s −1 . As mentioned in the previous section, this steady-state simulation is involved in computing the initial values for the DAE system during the optimization. Finally, the EL inflow is now dynamic and the hard-coded flow profile can be seen in the top left of Figure 4.
The solution of (OCP2) is carried out as described in Section 5.2.1 and the comparison between the solution using parareal and the solution via the fine integrator F is again of interest. The only difference in the solution setup is that fmincon is now also allowed to terminate as soon as the objective function value is below 10 −3 . With sight on the scale of the desired values SF d and SR d inside the objective function J 2 , this is justified as more precision is not required from an application point of view.
It is the aim of the optimization to adjust the plant activity to still produce the desired syngas ratio and flowrate under this changing inflow to the EL. Finally, the upper and lower bounds u u (t) and u l (t) at each time step coincide with the values from Section 5.2.1 with the only difference being that the first value of u u (t) (the upper bound to the total RWGS inflow) now is two times the value of the EL-profile from the top left of Figure 4 to again prevent negative H 2integration values. Based on the control values used in the steady state simulation, the initial values for the optimization are constant in time as u init (t) (465, 20) u .
The solutions to (OCP2) are investigated in Figures 4, 5. Here, only the results obtained using parareal are depicted. This is justified, as the results obtained using the fine integrator F are qualitatively the same as u fine − u par ||u fine || ≈ 2.1183 · 10 −4 , and F (u fine ) − F (u par ) ||F (u fine )|| ≈ 3.8393 · 10 −6 .
Regarding the solution time, the optimization using the fine integrator required 6.56 h and the optimization using parareal required 3.85 h resulting in a speed-up of 1.70. Due to the time-dependent dynamics, parareal now requires on average two iterations although the iteration error in line 12 of Algorithm 1 after the first iteration is on average 2.4956 · 10 −6 and thus very close to ε tol 10 −6 . The average times spent in the coarse and fine solver per iteration were 233.46 s and 31.16 s. These results are in line with the results from the previous experiment: the coarse/fine time per parareal iteration stems from the increased amount of time steps as well as the second average parareal iteration and the overall increased solution time comes from fmincon requiring more iterations due to (OCP2) being fully time-dependent including the EL profile depicted in the top left of Figure 4. Finally, do note that fmincon was terminated here as soon as the objective function value was lower than 10 −3 (further note that at termination the first-order optimality was around 1 · 10 −3 in both cases). As already mentioned (OCP2) targets a desired syngas flowrate and ratio, such that, due to the scale of the problem, more accuracy for this objective function is not required from an application point of view.
Regarding the solution quality, it can be seen in the bottom row of Figure 4 that the aim of the optimization was achieved: the desired syngas ratio and flowrate are matched very well under the prescribed change to the inlet flowrate to EL. In the top right of Figure 4 it can be seen how this was achieved: the inflow to the RWGS unit (blue) is adjusted in accordance to the inflow to EL (top left of Figure 4) and TRI (red) is then activated to still produce the desired amount of syngas. The H 2 -integration (yellow) is then a result of (11). Finally, the accumulation rate of O 2 in the buffering device (purple) at each time point can be calculated via (12), where O 2,TRI,in 0.168 and this value is a result from Table 2.
A more detailed investigation of the optimal state is provided in Figure 5. Here, the flowrates (top row) and flow compositions (bottom row) of the species at the RWGS reactor outlet (left column), the TRI reactor outlet (middle column), and at the combined plant outlet (right column) are depicted. Regarding the reactors, it is observed that the outlet compositions are essentially constant in time, consistent with the fact that equilibrium conditions are essentially met for the given operations. The changes in the outlet flowrates of the reactors are in agreement with the temporal profiles of the controls (the inlet flows to the reactors) depicted in the top right of Figure 4. Finally, it can be observed that an almost complete conversion of CH 4 was achieved and O 2 has entirely been depleted. The flowrate of H 2 and CO is kept constant in time, which was one aim of the optimization. Furthermore, if the content of steam in the raw syngas is constant over time, CO 2 drops at the minimum power input, suggesting that following purification units must handle variable loads of CO 2 over time before feeding the Fischer-Tropsch reactor.

Estimation of the Power Demand for Gas Compression
In Section 3.1 it was assumed that the power demand at EL is the only relevant contribution within the plant. In order to verify this statement, it is assumed that: • raw substances are fed to the control volume at ambient pressure; • Fischer-Tropsch (FT) synthesis occurs in a reactor operated at 25 bar; • the compression work is based on an isothermal transformation at 100°C.
The work demand related to the reactors to attain the pressure level required at the FT reactor, starting from ambient pressure (amb), is estimated as  The total power demand of the plant is reported in Figure 4A. As it can be seen, the power demand is largely dominated by the electrolyzer, which justifies the assumption made in Section 3.1 related to power consumption in the context of optimal control calculations.

CONCLUSION
A layout for a chemical plant that aims at a syngas production at constant flowrate and a H 2 -to-CO-ratio of 2 suitable for Fischer-Tropsch synthesis from renewable power, H 2 O and biogas was presented. The plant included a reverse water-gas shift (RWGS) reactor, which was to be used unless power shortage prevented the operation of the electrolyzer and thus the H 2 supply for the RWGS unit. A tri-reforming reactor was thus included into the plant to still provide the desired syngas flowrate and ratio at the product stream. An optimal control problem set to ensure a desired syngas ratio and flowrate at the plant outlet over a prescribed time horizon, given a variable power provision to the electrolyzer, was presented. From the mathematical side, parareal was utilized to speed up the expensive forward time integrations during the optimization and the continuous adjoint method was employed to compute the derivatives of the objective function. In a numerical investigation this novel approach was compared to the scenario where a traditional Runge-Kutta-Method was used for the time integration. It could be seen that both methods yielded qualitatively the same optimal solution, while the novel approach was significantly faster due to parareal. Furthermore, the objective of the optimal control problem was achieved: the desired syngas flowrate and ratio were matched over time under a prescribed piecewise linear inflow profile to the electrolyzer. From an engineering perspective, the numerical results allowed to identify appropriate feeding strategies at the reactors responding to a fluctuation in power supply and to monitor the temporal composition in the outlet stream and the accumulation of oxygen in the buffering device.
Future research paths should set the proposed layout in the framework of a thorough cost-based technological comparison with alternative buffering strategies for renewable power to syngas (buffering of H 2 ). Besides, appropriate O 2 -dosing strategies to decrease the temperature levels at tri-reforming should be accounted for, as well as a suitable strategy for the buffering of oxygen. The inclusion of heat exchangers and compressors would surely set feasibility bounds to the flexibility of the plant, an important aspect to be addressed. From a process-control perspective, it would be desirable to relax the condition of fixed temperature, pressure and composition at the reactors inlets throughout the control horizon. Furthermore, a mathematical investigation regarding the coarse solver inside parareal should be accounted for. Here, a faster solver could provide even more speed-up for the overall approach. As a result, a more in-depth analysis of the engineering model or even a real-time application could then become possible.

SUPPLEMENTARY MATERIAL
In this section, modeling details and parametrization is reported.

Reaction Rates for Tri-reforming of Methane
The parameters are selected from Xu and Froment (1989) and reported in , where the kinetics parameters k i and K j related to reaction and adsorption result from the following Arrhenius-like relations Values of coefficients A(k i , K j ) and E(k i , K j ) are listed in Table 4.

Reaction Rates for RWGS
RWGS kinetics are adapted from Richardson and Paripatyadar (1990 Effectiveness factor is set to 0.3, as reported in Richardson and Paripatyadar (1990).

Overall Heat Transfer Coefficient RWGS
The overall heat transfer coefficient reads 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).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
AM contributed Sections 2, 3, 4, and 7: he provided the plant layout, the governing equations as well as the corresponding kinetics and the optimal parameters for the different reactors. DG contributed Section 5: he implemented the optimal control approach as well as parareal, developed the combination of the continuous adjoint method and parareal and provided the main numerical results. AM is supervised by KS, DG is mentored by MS. SS aided in the optimal control approach.
FUNDING DG, SS, and MS acknowledge the financial support by the Federal Ministry of Education and Research of Germany with in the project P2Chem (support code 05M18OCB).